text:
struct PolynomialMapCubedSphereTriangulation{T} <: Triangulation{2,3}
  model::T
end

Gridap.Geometry.get_cell_coordinates(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_coordinates(get_grid(trian.model))

Gridap.Geometry.get_reffes(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_reffes(get_grid(trian.model))

Gridap.Geometry.get_cell_type(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_type(get_grid(trian.model))

Gridap.Geometry.get_node_coordinates(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_node_coordinates(get_grid(trian.model))

Gridap.Geometry.get_cell_node_ids(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_node_ids(get_grid(trian.model))

Gridap.Geometry.get_cell_map(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_map(trian.model)

Gridap.Geometry.get_grid(trian::PolynomialMapCubedSphereTriangulation) = Gridap.Geometry.get_grid(trian.model)

struct AnalyticalMapCubedSphereTriangulation{T} <: Triangulation{2,3}
  model::T
end

const CSDTT = Union{PolynomialMapCubedSphereTriangulation,AnalyticalMapCubedSphereTriangulation}

# Triangulation API
function Gridap.Geometry.get_facet_normal(trian::CSDTT)
  function _unit_outward_normal(v::Gridap.Fields.MultiValue{Tuple{2,3}},sign_flip::Bool)
    n1 = v[1,2]*v[2,3] - v[1,3]*v[2,2]
    n2 = v[1,3]*v[2,1] - v[1,1]*v[2,3]
    n3 = v[1,1]*v[2,2] - v[1,2]*v[2,1]
    n = VectorValue(n1,n2,n3)
    (-1)^sign_flip*n/norm(n)
  end

  # Get the Jacobian of the cubed sphere mesh
  map   = get_cell_map(trian)
  Jt    = lazy_map(∇,map)

  # Get the index of the panel for each element
  fl = get_face_labeling(trian.model)
  panel_id  = fl.d_to_dface_to_entity[3]
  sign_flip = [panel_id[i] == 25 || panel_id[i] == 21 || panel_id[i] == 24 for i=1:length(panel_id)]
  fsign_flip = lazy_map(Gridap.Fields.ConstantField,sign_flip)

  lazy_map(Operation(_unit_outward_normal),Jt,fsign_flip)
end

function Gridap.CellData.get_normal_vector(trian::CSDTT)
  cell_normal = Gridap.Geometry.get_facet_normal(trian)
  Gridap.CellData.GenericCellField(cell_normal,trian,ReferenceDomain())
end

# Delegating to the underlying face Triangulation

Gridap.Geometry.get_cell_coordinates(trian::AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_coordinates(get_grid(trian.model.cubed_sphere_linear_model))

Gridap.Geometry.get_reffes(trian::AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_reffes(get_grid(trian.model.cubed_sphere_linear_model))

Gridap.Geometry.get_cell_type(trian::AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_type(get_grid(trian.model.cubed_sphere_linear_model))

Gridap.Geometry.get_node_coordinates(trian::AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_node_coordinates(get_grid(trian.model.cubed_sphere_linear_model))

Gridap.Geometry.get_cell_node_ids(trian::AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_node_ids(get_grid(trian.model.cubed_sphere_linear_model))

Gridap.Geometry.get_cell_map(trian::AnalyticalMapCubedSphereTriangulation) = trian.model.cell_map

Gridap.Geometry.get_grid(trian::AnalyticalMapCubedSphereTriangulation) = get_grid(trian.model.cubed_sphere_linear_model)

Gridap.Geometry.get_background_model(trian::AnalyticalMapCubedSphereTriangulation) = trian.model

function Gridap.Geometry.get_glue(a::AnalyticalMapCubedSphereTriangulation,D::Val{2})
  nc=num_cells(a.model)
  tface_to_mface=Gridap.Fields.IdentityVector(nc)
  tface_to_mface_map=Fill(Gridap.Fields.GenericField(identity),nc)
  mface_to_tface=tface_to_mface
  Gridap.Geometry.FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
@testset "subgraphs" begin

d = parse_dot("""
graph {
  subgraph {
    a -- b
    a -- c
  }
  d -- e
}
""")[1]
@test nodes(d) == Set(["a", "b", "c", "d", "e"])
@test edges(d) == Set([("d","e"),("a","c"),("a","b")])

d = parse_dot("""
graph {
  subgraph {
    a -- b
    a -- c
  }
  d; e;
}
""")[1]
@test nodes(d) == Set(["a", "b", "c", "d", "e"])
@test edges(d) == Set([("a","c"),("a","b")])

d = parse_dot("""
graph {
  a -- subgraph {
    b -- c
    c -- d
  } -- e
}
""")[1]
@test nodes(d) == Set(["a", "b", "c", "d", "e"])
@test edges(d) == Set([("a","b"),("a","c"),("a","d"),("b","c"),("b","e"),("c","d"),("c","e"),("d","e")])

d = parse_dot("""
graph {
  subgraph { a -- b } -- subgraph { c -- d }
}
""")[1]
@test nodes(d) == Set(["a", "b", "c", "d"])
@test edges(d) == Set([("a","b"),("a","c"),("a","d"),("b","c"),("b","d"),("c","d")])

println("subgraphs ok")

end

ext: jl


text:
using ACE, Printf, Test, LinearAlgebra
using ACE: PolyTransform, rbasis,
using JuLIP: evaluate, evaluate_d
using ACE.JacobiPolys: Jacobi
using PyPlot

trans = PolyTransform(1, 1)
B2 = rbasis(4, trans, 2, 0.5, 2.0)

rr = range(0, 2.5, length=100)
BB = zeros(100, length(B2))
for (i, r) in enumerate(rr)
   BB[i, :] = evaluate(B2, r)
end

##
clf()
for i = 1:length(B2)
   plot(rr, BB[:,i], label = "P_$i")
end
legend()
xlabel("r")
ylabel("P_i(r) = J_i(x) * w(x)")

PyPlot.display_figs()


text:
#===================================================================================================
Constraint handling differential evolution
===================================================================================================#
"""
    bounce_back!(t::AV{T}, b::AV{T}, bounds::AM{<:Real}) where T<:AF

Relocates randomly the parameter of the trial vector `t` in between the bound it exceeded and the
corresponding parameter from the base vector `b`.
"""
function bounce_back!(t::AV{T}, b::AV{T}, bounds::AM{<:Real}) where T<:AF
    @assert length(t) == length(b) == size(bounds, 1)
    @inbounds for i in eachindex(t)
        lower, upper = bounds[i, 1], bounds[i, 2]
        if t[i] < lower  # smaller than the lower bound
            t[i] = lower + rand(T) * (b[i] - lower)
        elseif t[i] > upper # exceed the upper bound
            t[i] = upper - rand(T) * (upper - b[i])
        end
    end
end


text:
# Nibble is a half byte (4 bits).
# This file defines the Nibble type and related stuff

"""
Nibble is basically an unsigned integer stored in a half byte, giving the data range 0..15

It is defined as an explicit type to allow for clear method signatures.
"""
#primitive type Nibble <: Number 8 end
primitive type Nibble <: Unsigned 8 end

function Nibble(v::Integer) ::Nibble
    @boundscheck checkulimit(v,0x0F)
    reinterpret(Nibble,v%UInt8)
end

"""
    nibble(x::T)

For types T which include a nibble field, define nibble(..)
to extract that field.

Recommended for bitfield-like types,  which have one exactly one Nibble field.

"""
function nibble end

Base.convert(::Type{UInt8}, x::Nibble) = reinterpret(UInt8,x)

Base.convert(::Type{Nibble}, x::UInt8) = Nibble(x)

Base.UInt8(x::Nibble) = reinterpret(UInt8,x)

Base.UInt64(x::Nibble) = reinterpret(UInt8,x)%UInt64

Base.Int64(x::Nibble) = reinterpret(UInt8,x)%Int64

Base.promote_rule(::Type{Nibble}, ::Type{U}) where U<:Unsigned  = UInt64

Base.promote_type(::Type{Nibble}, ::Type{Nibble}) = UInt8

Base.promote_rule(::Type{Nibble}, ::Type{U}) where U<:Signed = Int64

function Base.Char(v::Nibble)
    code = UInt8(v)
    if code>0x09
        code += 0x07
    end
    return '0'+code
end

function Base.show(io::IO, v::Nibble)
    print(io, Char(v))
end

#Base.promote_rule(::Type{T}, ::Type{Nibble}) where T<:Enum  = T

#Base.promote_rule(::Type{Nibble}, ::Type{T}) where T<:Enum  = Nibble

# Base.convert(::Type{T}, x::Nibble) where T<:Enum = T(UInt8(x))

Base.convert(::Type{Nibble}, x::T) where T<:Enum = Nibble(x)

Nibble(v::T) where T<:Enum = Nibble(UInt8(v))

Base.promote_rule(::Type{Nibble}, ::Type{T}) where T<:Enum  = Nibble


text:
@testset "Print/show" begin
    # Test whether we can reconstruct the objects from their show values

    "Simulates a user printing an object in the REPL, copying and pasting the result and pressing enter"
    function show_and_eval(object)
        io = IOBuffer()
        show(io, object)
        io |> take! |> String |> Meta.parse |> eval
    end

    acoustic = Acoustic(1.89, 1.5, 2)
    circle = Circle(2.0)
    particle = Particle(acoustic, circle)
    simulation = FrequencySimulation([particle], plane_source(acoustic))

    @test acoustic == show_and_eval(acoustic)
    @test circle == show_and_eval(circle)
    @test particle == show_and_eval(particle)

    # Currently can't do this with FrequencySimulations because of Source type
    @test show(devnull, simulation) == nothing

end


text:
@testset "Open-loop disturbance" begin

    nw = 4
    nx = 6
    H = 4
    N_sample = 3
    H_sim = N_sample*H
    w = [rand(nw) for t=1:H]
    d = open_loop_disturbances(w, N_sample)

    idx = [1,1,1,2,2,2,3,3,3,4,4,4]
    cnt = [1,2,3,1,2,3,1,2,3,1,2,3]
    for t = 1:H_sim
        x = rand(nx)
        wt = disturbances(d, x, t)
        @test d.idx == idx[t]
        @test d.cnt == cnt[t]
        @test norm(wt - d.w[d.idx] / d.N_sample) < 1e-10
    end

end

@testset "Open-loop policy" begin

    nu = 4
    nx = 6
    H = 4
    N_sample = 3
    H_sim = N_sample*H
    u = [rand(nu) for t=1:H]
    p = open_loop_policy(u, N_sample=N_sample)

    idx = [1,1,1,2,2,2,3,3,3,4,4,4]
    cnt = [1,2,3,1,2,3,1,2,3,1,2,3]
    for t = 1:H_sim
        x = nothing
        traj = nothing
        ut = policy(p, x, traj, t)
        @test p.idx == idx[t]
        @test p.cnt == cnt[t]
        @test norm(ut - p.u[p.idx] / p.N_sample) < 1e-10
    end

end


text: compute(n::Int)::Int = n * (n + 1) * (3n ^ 2 - n - 2) ÷ 12

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
export MaestroDataset, reader, @dataset_str

using Pkg.TOML
using DataDeps
using Random

using Transformers.Datasets
import Transformers.Datasets: trainfile, testfile, devfile, Mode

const dataset_configs = open(TOML.parse, joinpath(@__DIR__, "datasets.toml"))

# For clarity
macro dataset_str(name::String)
    :(@datadep_str($name))
end

struct MaestroDataset <: Dataset end

trainfile(::MaestroDataset) = dataset"Maestro_LM/train.jld2"
devfile(::MaestroDataset) = dataset"Maestro_LM/validation.jld2"
testfile(::MaestroDataset) = dataset"Maestro_LM/test.jld2"

function reader(::Type{M}, ds::MaestroDataset; batchsize=32, shuffle=true) where {M <: Mode}
    @assert batchsize >= 1

    chan = Channel{Vector{Vector{Int}}}()
    filename = datafile(M, ds)
    num_examples = load(filename, "num_examples")
    load_order = shuffle ? randperm(num_examples) : collect(1:num_examples)

    loader = @async begin
        for batch_load_order in Iterators.partition(load_order, batchsize)
            length(batch_load_order) != batchsize && continue
            data = collect(load(filename, string.(batch_load_order)...))
            put!(chan, ifelse(batchsize == 1, [data], data))
        end
    end

    bind(chan, loader)
    chan
end

function register_datasets(configs)
    for (dataset, config) in pairs(configs)
        description = config["description"]
        checksum = config["checksum"]
        url = config["url"]
        if startswith(url, "https://docs.google.com")
            fetch_method = Transformers.Datasets.download_gdrive
        else
            fetch_method = DataDeps.fetch_http
        end
        dep = DataDep(dataset, description, url, checksum;
                      fetch_method=fetch_method, post_fetch_method=DataDeps.unpack)
        DataDeps.register(dep)
    end
end

function list_datasets()
    println.(keys(dataset_configs))
    return
end


text:
# Map moment names into functions that make raw file paths
raw_file_map() = Dict([
    # :fracDrop_qtM => raw_frac_drop_qual_gpa,
    :fracEnroll_qV => raw_entry_qual_parental,
    :fracEnter_gV => raw_entry_gpa_parental,
    :fracEnter_pV => raw_entry_gpa_parental,
    :fracEnter_gpM => raw_entry_gpa_parental,
    :fracGrad => raw_grad_rate_qual_gpa,
    :fracGrad_gV => raw_grad_rate_qual_gpa,
    :fracGrad_qV => raw_grad_rate_qual_gpa,
    :fracGrad_gpM => raw_frac_grad_qual_parental,
    :fracGrad_qgM => raw_grad_rate_qual_gpa,
    :mass_gpM => raw_mass_gpa_parental,
    :massEntry_qgM => raw_mass_entry_qual_gpa,
    :timeToDrop4y_gV => raw_time_to_drop_4y_qual_gpa,
    :timeToDrop4y_qV => raw_time_to_drop_4y_qual_gpa,
    :timeToDrop4y_qgM => raw_time_to_drop_4y_qual_gpa,
    :timeToGrad4y_qV => raw_time_to_grad_4y_qual_gpa,
    :timeToGrad4y_gV => raw_time_to_grad_4y_qual_gpa,
    :timeToGrad4y_pV => raw_time_to_grad_4y_qual_parental,
    :timeToGrad4y_qgM => raw_time_to_grad_4y_qual_gpa,
    :timeToGrad4y_qpM => raw_time_to_grad_4y_qual_parental,
    :workTime_pV => raw_work_hours_qual_parental,
    :workTime_qV => raw_work_hours_qual_parental,
    :workTime_qpM => raw_work_hours_qual_parental
]);

## ---------  Methods

# Mapping from characteristics to directory names
dSelfTranscript = Dict([:selfReport => "SelfReport",  :transcript => "Transcripts"]);

dMomentType = Dict([:mean => "Means",  :count => "Counts",
    :std => "StandardDeviations",
	:regression => "Regressions"]);

dGroup = Dict([:finance => "Financing",  :freshmen => "Fresh_Char",
	:hsGrads => "HS_Char",  :progress => "Progress",  :none => ""]);

data_settings(rf :: RawDataFile) = rf.ds;
set_moment_type(rf :: RawDataFile, momentType :: Symbol) =
    rf.momentType = momentType;

## ---------  Mapping moments to raw data files

"""
    $(SIGNATURES)

Retrieve raw file path from moment name. The optional `momentType` argument allows the user to load counts or std deviations instead of means for moments where those exist.
"""
function raw_file_path(ds :: DataSettings, mName :: Symbol; momentType = nothing)
    rf = raw_file(ds, mName; momentType);
    return data_file(rf)
end

# Make a RawDataFile from a moment name
function raw_file(ds :: DataSettings, mName :: Symbol; momentType = nothing)
    fileMap = raw_file_map();
    rf_fct = fileMap[mName];
    rf = rf_fct(ds; momentType = momentType);
    return rf
end

# File with counts for a given data moment
function count_file(ds :: DataSettings, mName :: Symbol)
    rf = raw_file(ds, mName);
    set_moment_type(rf, MtCount);
    return rf
end

function std_file(ds :: DataSettings, mName :: Symbol)
    rf = raw_file(ds, mName);
    set_moment_type(rf, MtStd);
    return rf
end

"""
	$(SIGNATURES)

Make a list of missing data files. Returns a String vector.
Only handles files in `raw_file_map`.
This is not perfect. It only checks for files containing means (not counts or std deviations). But useful for spotting that files have been unexpectedly renamed.
"""
function missing_file_list(ds)
    missList = Vector{String}();
    # Mapping of moments to files
    rawMap = raw_file_map();
    for mName in keys(rawMap)
        rf = raw_file(ds, mName; momentType = MtMean);
        @assert isa(rf, CollegeStratData.RawDataFile)
        fPath = raw_file_path(ds, mName; momentType = MtMean);
        if !isfile(fPath)
            push!(missList, fPath);
        end
    end
    return missList
end

"""
	$(SIGNATURES)

Generate a list of files that differ (in name only) between existing data (as stored in `DataCollegeStrat` somewhere in `.julia`) and new data (in `Dropbox`).
This should be run before importing the new data files.

Excludes files that refer to public or private colleges.
"""
function file_name_differences(ds :: DataSettings)
    newDir = dropbox_dir(ds);
    @assert isdir(newDir);

    fName = "file_name_differences.txt";
    logFn = joinpath(out_dir(ds), fName);
    open(logFn, "w") do io
        println(io, "\nFile name differences for $ds:");
        println(io, "  Comparing existing data dir with new data in Dropbox.");
        println(io, "  Old:  $(data_dir(ds))");
        println(io, "  New:  $newDir");
        dir_diff_report(data_dir(ds), newDir; io = io,
            exclude = ["_PRI", "_PUB"]);  # add io argument
    end

    # Mirror on screen
    for line in eachline(logFn)
        println(line);
    end
    println("Results written to $logFn");
    return nothing
end

# Base dir for new data files in dropbox. Only used for checking new files before import.
dropbox_base_dir() = "/Users/lutz/Dropbox/Dropout Policies/Data/empiricaltargets";

dropbox_dir(ds :: DataSettings) = joinpath(dropbox_base_dir(), data_sub_dir(ds));

"""
	$(SIGNATURES)

Compare all `dat` files in existing data directory with the corresponding `Dropbox` directory. Useful before new data files are imported.
"""
function compare_dirs(ds :: DataSettings)
    fName = "compare_dirs.txt";
    logFn = joinpath(out_dir(ds), fName);
    open(logFn, "w") do io
        compare_dirs(data_dir(ds), dropbox_dir(ds); io = io);
    end

    # Mirror on screen
    for line in eachline(logFn)
        println(line);
    end
    println("Results written to $logFn");
end

"""
	$(SIGNATURES)

Compare all `dat` files for two base directories. Report files where headers do not match. Useful for checking new data files before importing.
"""
function compare_dirs(dir1 :: AbstractString, dir2 :: AbstractString;
    io = stdout)

    @assert isdir(dir1)
    @assert isdir(dir2)
    fList = files_in_dir(dir1);
    if isempty(fList)
        println(io, "Empty directory");
        return nothing
    end

    allValid = true;
    for f ∈ fList
        _, fExt = splitext(f);
        if fExt ∈ (".dat", )
            newFile = joinpath(dir2, f);
            if !isfile(newFile)
                println(io, "Not found: $newFile");
            else
                isValid = compare_delim_files(joinpath(dir1, f), newFile; io = io);
                if !isValid
                    allValid = false;
                    println(io, "  in $f");
                end
            end
        end
    end
    return allValid
end

"""
	$(SIGNATURES)

Compare two delimited files. Report discrepancies in headers and dimensions.
Useful for checking new data files before importing.
Tries to never error. Just report.
"""
function compare_delim_files(fn1, fn2; io = stdout)
    isValid = true;
    if !isfile(fn1)
        println(io, "Not found: $fn1");
        isValid = false;
    elseif !isfile(fn2)
        println(io, "Not found: $fn2");
        isValid = false;
    end
    if isValid
        m1 = readdlm(fn1);
        m2 = readdlm(fn2);
        if !isequal(m1[1,:], m2[1,:])
            println(io, "Col header mismatch");
            println(io, "    Old: ", m1[1,:]);
            println(io, "    New: ", m2[1,:]);
            isValid = false;
        elseif !isequal(m1[:,1], m2[:, 1])
            println(io, "Row header mismatch");
            isValid = false;
        end
    end
    return isValid
end

## --------  Directories

# Raw data dir with characteristics of each moment
# raw_data_dir(rf :: RawDataFile) =
# 	raw_data_dir(rf.selfOrTranscript, rf.momentType, rf.group);

# function raw_data_dir(ds :: DataSettings, selfReportOrTranscript :: Symbol,
#     momentType :: Symbol, 	momentGroup :: Symbol)

# 	rawPath = joinpath(raw_data_base_dir(ds),
# 		data_sub_dir(selfReportOrTranscript, momentType, momentGroup));
# 	@assert isdir(rawPath)  "Directory not found: $rawPath"
# 	return rawPath
# end

"""
	$(SIGNATURES)

Path for a raw data file.
"""
data_file(rf :: RawDataFile) =
	joinpath(data_dir(data_settings(rf)), data_sub_dir(rf), file_name(rf));

# Subdir relative to `raw_data_base_dir` or `data_dir`
data_sub_dir(rf :: RawDataFile) =
	data_sub_dir(rf.selfOrTranscript, rf.momentType, rf.group);

data_sub_dir(selfReportOrTranscript :: Symbol,  momentType :: Symbol,
	momentGroup :: Symbol) =
    joinpath(dSelfTranscript[selfReportOrTranscript],  "dat_files",
		dMomentType[momentType],  dGroup[momentGroup])

function file_name(rf :: RawDataFile)
    fn = rf.rawFile;
    if rf.momentType == :count
        # `_N` appended to name (why??)
        fn2, fExt = splitext(fn);
        fn = fn2 * "_N" * fExt;
    elseif rf.momentType == :std
        # `_SD` appended to name (why??)
        fn2, fExt = splitext(fn);
        fn = fn2 * "_SD" * fExt;
    end
    return fn
end

# """
#     $(SIGNATURES)

# Copy raw data files from Dropbox to local dir.
# This only needs to be done when data files get updated.
# Not all moments have raw data files. Those are skipped.
# """
# function copy_raw_data_files(ds :: DataSettings; trialRun :: Bool = false)
#     println("\nCopying raw data files to local dir");
#     srcDir = raw_data_base_dir(ds);
#     tgDir = data_dir(ds);
#     rsync_dir(srcDir, tgDir;  trialRun = trialRun,  doDelete = false);
# end

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


text:
module BitArrayTests2

using Base.Test
using Base: findprevnot, findnextnot

tc{N}(r1::NTuple{N,Any}, r2::NTuple{N,Any}) = all(x->tc(x...), [zip(r1,r2)...])
tc{N}(r1::BitArray{N}, r2::Union{BitArray{N},Array{Bool,N}}) = true
tc(r1::RowVector{Bool,BitVector}, r2::Union{RowVector{Bool,BitVector},RowVector{Bool,Vector{Bool}}}) = true
tc{T}(r1::T, r2::T) = true
tc(r1,r2) = false

bitcheck(b::BitArray) = Base._check_bitarray_consistency(b)
bitcheck(x) = true

function check_bitop_call(ret_type, func, args...)
    r1 = func(args...)
    r2 = func(map(x->(isa(x, BitArray) ? Array(x) : x), args)...)
    ret_type ≢ nothing && !isa(r1, ret_type) && @show ret_type, r1
    ret_type ≢ nothing && @test isa(r1, ret_type)
    @test tc(r1, r2)
    @test isequal(r1, ret_type ≡ nothing ? r2 : convert(ret_type, r2))
    @test bitcheck(r1)
end
macro check_bit_operation(ex, ret_type)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop_call, esc(ret_type), map(esc, ex.args)...)
end
macro check_bit_operation(ex)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop_call, nothing, map(esc, ex.args)...)
end

let t0 = time()
    global timesofar
    function timesofar(str)
        return # no-op, comment to see timings
        t1 = time()
        println(str, ": ", t1-t0, " seconds")
        t0 = t1
    end
end

# vectors size
v1 = 260
# matrices size
n1, n2 = 17, 20
# arrays size
s1, s2, s3, s4 = 5, 8, 3, 7

allsizes = [((), BitArray{0}), ((v1,), BitVector),
            ((n1,n2), BitMatrix), ((s1,s2,s3,s4), BitArray{4})]

@testset "Binary arithmetic operators" begin
    @testset "Matrix{Bool}/Matrix{Bool}" begin
        b1 = bitrand(n1, n2)
        b2 = bitrand(n1, n2)
        @check_bit_operation broadcast(&, b1, b2)  BitMatrix
        @check_bit_operation broadcast(|, b1, b2)  BitMatrix
        @check_bit_operation broadcast(xor, b1, b2)  BitMatrix
        @check_bit_operation (+)(b1, b2)  Matrix{Int}
        @check_bit_operation (-)(b1, b2)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, b2) BitMatrix
        @check_bit_operation broadcast(/, b1, b2) Matrix{Float64}
        @check_bit_operation broadcast(^, b1, b2) BitMatrix
        @check_bit_operation (/)(b1,1) Matrix{Float64}

        b2 = trues(n1, n2)
        @check_bit_operation broadcast(div, b1, b2) BitMatrix
        @check_bit_operation broadcast(mod, b1, b2) BitMatrix
        @check_bit_operation broadcast(div, b1, Array(b2)) BitMatrix
        @check_bit_operation broadcast(mod, b1, Array(b2)) BitMatrix
        @check_bit_operation broadcast(div, Array(b1), b2) BitMatrix
        @check_bit_operation broadcast(mod, Array(b1), b2) BitMatrix

        b1 = bitrand(n1, n1)
        while abs(det(Array{Float64}(b1))) ≤ 1e-6
            b1 = bitrand(n1, n1)
        end
        b2 = bitrand(n1, n1)

        @check_bit_operation (*)(b1, b2) Matrix{Int}
        @check_bit_operation (/)(b1, b1) Matrix{Float64}
        @check_bit_operation (\)(b1, b1) Matrix{Float64}

        b0 = falses(0)
        @check_bit_operation broadcast(&, b0, b0)  BitVector
        @check_bit_operation broadcast(|, b0, b0)  BitVector
        @check_bit_operation broadcast(xor, b0, b0)  BitVector
        @check_bit_operation broadcast(*, b0, b0) BitVector
        @check_bit_operation (*)(b0, b0') BitMatrix
    end

    @testset "Matrix{Bool}/Matrix{Int}" begin
        b1 = bitrand(n1, n2)
        i2 = rand(1:10, n1, n2)
        @check_bit_operation broadcast(&, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(|, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(xor, b1, i2)  Matrix{Int}
        @check_bit_operation (+)(b1, i2)  Matrix{Int}
        @check_bit_operation (-)(b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, i2) Matrix{Int}
        @check_bit_operation broadcast(/, b1, i2) Matrix{Float64}
        @check_bit_operation broadcast(^, b1, i2) BitMatrix
        @check_bit_operation broadcast(div, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(mod, b1, i2)  Matrix{Int}
    end

    @testset "Matrix{Bool}/Matrix{Float64}" begin
        b1 = bitrand(n1, n2)
        f2 = 1.0 .+ rand(n1, n2)
        @check_bit_operation broadcast(*, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(/, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(^, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(mod, b1, f2)  Matrix{Float64}
    end

    @testset "Number/Matrix" begin
        b2 = bitrand(n1, n2)
        i1 = rand(1:10)
        u1 = UInt8(i1)
        f1 = Float64(i1)
        ci1 = complex(i1)
        cu1 = complex(u1)
        cf1 = complex(f1)

        @check_bit_operation broadcast(&, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(|, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(xor, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(+, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(-, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(*, i1, b2) Matrix{Int}

        @check_bit_operation broadcast(&, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(|, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(xor, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(+, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(-, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(*, u1, b2) Matrix{UInt8}

        for (x1,t1) = [(f1, Float64),
                       (ci1, Complex{Int}),
                       (cu1, Complex{UInt8}),
                       (cf1, Complex128)]
            @check_bit_operation broadcast(+, x1, b2)  Matrix{t1}
            @check_bit_operation broadcast(-, x1, b2)  Matrix{t1}
            @check_bit_operation broadcast(*, x1, b2) Matrix{t1}
        end

        b2 = trues(n1, n2)
        @check_bit_operation broadcast(/, true, b2)  Matrix{Float64}
        @check_bit_operation broadcast(div, true, b2)   BitMatrix
        @check_bit_operation broadcast(mod, true, b2)   BitMatrix
        @check_bit_operation broadcast(/, false, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, false, b2)  BitMatrix
        @check_bit_operation broadcast(mod, false, b2)  BitMatrix

        @check_bit_operation broadcast(/, i1, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(mod, i1, b2)  Matrix{Int}

        @check_bit_operation broadcast(/, u1, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(mod, u1, b2)  Matrix{UInt8}

        @check_bit_operation broadcast(/, f1, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, f1, b2)  Matrix{Float64}
        @check_bit_operation broadcast(mod, f1, b2)  Matrix{Float64}

        @check_bit_operation broadcast(/, ci1, b2) Matrix{Complex128}
        @check_bit_operation broadcast(/, cu1, b2) Matrix{Complex128}
        @check_bit_operation broadcast(/, cf1, b2) Matrix{Complex128}

        b2 = bitrand(n1, n2)
        @check_bit_operation broadcast(^, false, b2) BitMatrix
        @check_bit_operation broadcast(^, true, b2)  BitMatrix
        @check_bit_operation broadcast(^, 0x0, b2)   Matrix{UInt8}
        @check_bit_operation broadcast(^, 0x1, b2)   Matrix{UInt8}
        @check_bit_operation broadcast(^, -1, b2)    Matrix{Int}
        @check_bit_operation broadcast(^, 0, b2)     Matrix{Int}
        @check_bit_operation broadcast(^, 1, b2)     Matrix{Int}
        @check_bit_operation broadcast(^, 0.0, b2)   Matrix{Float64}
        @check_bit_operation broadcast(^, 1.0, b2)   Matrix{Float64}
        @check_bit_operation broadcast(^, 0.0im, b2) Matrix{Complex128}
        @check_bit_operation broadcast(^, 1.0im, b2) Matrix{Complex128}
        @check_bit_operation broadcast(^, 0im, b2)   Matrix{Complex{Int}}
        @check_bit_operation broadcast(^, 1im, b2)   Matrix{Complex{Int}}
        @check_bit_operation broadcast(^, 0x0im, b2) Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(^, 0x1im, b2) Matrix{Complex{UInt8}}
    end

    @testset "Matrix/Number" begin
        b1 = bitrand(n1, n2)
        i2 = rand(1:10)
        u2 = UInt8(i2)
        f2 = Float64(i2)
        ci2 = complex(i2)
        cu2 = complex(u2)
        cf2 = complex(f2)
        b2 = Array(bitrand(n1,n2))

        @check_bit_operation broadcast(&, b1, true)   BitMatrix
        @check_bit_operation broadcast(&, b1, false)  BitMatrix
        @check_bit_operation broadcast(&, true, b1)   BitMatrix
        @check_bit_operation broadcast(&, false, b1)  BitMatrix
        @check_bit_operation broadcast(|, b1, true)   BitMatrix
        @check_bit_operation broadcast(|, b1, false)  BitMatrix
        @check_bit_operation broadcast(|, true, b1)   BitMatrix
        @check_bit_operation broadcast(|, false, b1)  BitMatrix
        @check_bit_operation broadcast(xor, b1, true)   BitMatrix
        @check_bit_operation broadcast(xor, b1, false)  BitMatrix
        @check_bit_operation broadcast(xor, true, b1)   BitMatrix
        @check_bit_operation broadcast(xor, false, b1)  BitMatrix
        @check_bit_operation broadcast(+, b1, true)   Matrix{Int}
        @check_bit_operation broadcast(+, b1, false)  Matrix{Int}
        @check_bit_operation broadcast(-, b1, true)   Matrix{Int}
        @check_bit_operation broadcast(-, b1, false)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, true)  BitMatrix
        @check_bit_operation broadcast(*, b1, false) BitMatrix
        @check_bit_operation broadcast(*, true, b1)  BitMatrix
        @check_bit_operation broadcast(*, false, b1) BitMatrix
        @check_bit_operation broadcast(/, b1, true)  Matrix{Float64}
        @check_bit_operation broadcast(/, b1, false) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, true)   BitMatrix
        @check_bit_operation broadcast(mod,b1, true)    BitMatrix

        @check_bit_operation broadcast(&, b1, b2)  BitMatrix
        @check_bit_operation broadcast(|, b1, b2)  BitMatrix
        @check_bit_operation broadcast(xor, b1, b2)  BitMatrix
        @check_bit_operation broadcast(&, b2, b1)  BitMatrix
        @check_bit_operation broadcast(|, b2, b1)  BitMatrix
        @check_bit_operation broadcast(xor, b2, b1)  BitMatrix
        @check_bit_operation broadcast(&, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(|, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(xor, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(+, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(-, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, i2) Matrix{Int}
        @check_bit_operation broadcast(/, b1, i2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(mod, b1, i2)  Matrix{Int}

        @check_bit_operation broadcast(&, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(|, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(xor, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(+, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(-, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(*, b1, u2) Matrix{UInt8}
        @check_bit_operation broadcast(/, b1, u2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(mod, b1, u2)  Matrix{UInt8}

        @check_bit_operation broadcast(+, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(-, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(*, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(/, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(mod, b1, f2)  Matrix{Float64}

        @check_bit_operation broadcast(+, b1, ci2)  Matrix{Complex{Int}}
        @check_bit_operation broadcast(-, b1, ci2)  Matrix{Complex{Int}}
        @check_bit_operation broadcast(*, b1, ci2) Matrix{Complex{Int}}
        @check_bit_operation broadcast(/, b1, ci2) Matrix{Complex128}

        @check_bit_operation broadcast(+, b1, cu2)  Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(-, b1, cu2)  Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(*, b1, cu2) Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(/, b1, cu2) Matrix{Complex128}

        @check_bit_operation broadcast(+, b1, cf2)  Matrix{Complex128}
        @check_bit_operation broadcast(-, b1, cf2)  Matrix{Complex128}
        @check_bit_operation broadcast(*, b1, cf2) Matrix{Complex128}
        @check_bit_operation broadcast(/, b1, cf2) Matrix{Complex128}

        @check_bit_operation broadcast(^, b1, false) BitMatrix
        @check_bit_operation broadcast(^, b1, true)  BitMatrix
        @check_bit_operation broadcast(^, b1, 0x0)   BitMatrix
        @check_bit_operation broadcast(^, b1, 0x1)   BitMatrix
        @check_bit_operation broadcast(^, b1, 0)     BitMatrix
        @check_bit_operation broadcast(^, b1, 1)     BitMatrix
        @check_bit_operation broadcast(^, b1, -1.0)  Matrix{Float64}
        @check_bit_operation broadcast(^, b1, 0.0)   Matrix{Float64}
        @check_bit_operation broadcast(^, b1, 1.0)   Matrix{Float64}
        @check_bit_operation broadcast(^, b1, 0.0im) Matrix{Complex128}
        @check_bit_operation broadcast(^, b1, 0x0im) Matrix{Complex128}
        @check_bit_operation broadcast(^, b1, 0im)   Matrix{Complex128}
        @test_throws DomainError broadcast(^, b1, -1)

        b1 = trues(n1, n2)
        @check_bit_operation broadcast(^, b1, -1.0im) Matrix{Complex128}
        @check_bit_operation broadcast(^, b1, 1.0im)  Matrix{Complex128}
        @check_bit_operation broadcast(^, b1, -1im)   Matrix{Complex128}
        @check_bit_operation broadcast(^, b1, 1im)    Matrix{Complex128}
        @check_bit_operation broadcast(^, b1, 0x1im)  Matrix{Complex128}
    end
end

timesofar("binary arithmetic")

@testset "Binary comparison operators" begin
    b1 = bitrand(n1, n2)
    b2 = bitrand(n1, n2)
    @check_bit_operation broadcast(==, b1, b2) BitMatrix
    @check_bit_operation broadcast(!=, b1, b2) BitMatrix
    @check_bit_operation broadcast(<, b1, b2) BitMatrix
    @check_bit_operation broadcast(<=, b1, b2) BitMatrix
end

timesofar("binary comparison")

@testset "Data movement" begin
    b1 = bitrand(s1, s2, s3, s4)
    for d = 1:4
        j = rand(1:size(b1, d))
        #for j = 1 : size(b1, d)
            @check_bit_operation slicedim(b1, d, j) BitArray{3}
        #end
        @check_bit_operation flipdim(b1, d) BitArray{4}
    end
    @test_throws ArgumentError flipdim(b1, 5)

    b1 = bitrand(n1, n2)
    for k = 1:4
        @check_bit_operation rotl90(b1, k) BitMatrix
    end

    for m = 0:v1
        b1 = bitrand(m)
        @check_bit_operation reverse(b1) BitVector
    end

    b1 = bitrand(v1)
    for m = [rand(1:v1)-1, 0, 1, 63, 64, 65, 191, 192, 193, v1-1]
        @test isequal(b1 << m, [ b1[m+1:end]; falses(m) ])
        @test isequal(b1 >>> m, [ falses(m); b1[1:end-m] ])
        @test isequal(b1 << -m, b1 >> m)
        @test isequal(b1 >>> -m, b1 << m)
        @test isequal(rol(b1, m), [ b1[m+1:end]; b1[1:m] ])
        @test isequal(ror(b1, m), [ b1[end-m+1:end]; b1[1:end-m] ])
        @test isequal(ror(b1, m), rol(b1, -m))
        @test isequal(rol(b1, m), ror(b1, -m))
    end

    b = bitrand(v1)
    i = bitrand(v1)
    for m = [rand(1:v1), 63, 64, 65, 191, 192, 193, v1-1]
        j = rand(1:m)
        b1 = ror!(i, b, j)
        i1 = ror!(b, j)
        @test b1 == i1
        b2 = rol!(i1, b1, j)
        i2 = rol!(b1, j)
        @test b2 == i2
    end
end

timesofar("datamove")

@testset "countnz & find" begin
    for m = 0:v1, b1 in Any[bitrand(m), trues(m), falses(m)]
        @check_bit_operation countnz(b1) Int

        @check_bit_operation findfirst(b1) Int

        @check_bit_operation findfirst(b1, true)  Int
        @check_bit_operation findfirst(b1, false) Int
        @check_bit_operation findfirst(b1, 3)     Int

        @check_bit_operation findfirst(x->x, b1)     Int
        @check_bit_operation findfirst(x->!x, b1)    Int
        @check_bit_operation findfirst(x->true, b1)  Int
        @check_bit_operation findfirst(x->false, b1) Int

        @check_bit_operation find(b1) Vector{Int}
    end

    b1 = trues(v1)
    for i = 0:(v1-1)
        @test findfirst(b1 >> i) == i+1
        @test Base.findfirstnot(.~(b1 >> i)) == i+1
    end

    for i = 3:(v1-1), j = 2:i
        submask = b1 << (v1-j+1)
        @test findnext((b1 >> i) .| submask, j) == i+1
        @test findnextnot((.~(b1 >> i)) .⊻ submask, j) == i+1
    end

    b1 = bitrand(n1, n2)
    @check_bit_operation findnz(b1) Tuple{Vector{Int}, Vector{Int}, BitArray}
end

timesofar("nnz&find")

@testset "Findnext/findprev" begin
    b1 = trues(v1)
    b2 = falses(v1)
    for i = 1:v1
        @test findprev(b1, i) == findprev(b1, true, i) == findprev(identity, b1, i)
        @test findprevnot(b2, i) == findprev(!, b2, i) == i
    end

    odds = broadcast(isodd, 1:2000)
    evens = broadcast(iseven, 1:2000)

    for i = 1:2:2000
        @test findprev(odds,i)  == findprevnot(evens,i) == i
        @test findnext(odds,i)  == findnextnot(evens,i) == i
        @test findprev(evens,i) == findprevnot(odds,i)  == i-1
        @test findnext(evens,i) == findnextnot(odds,i)  == (i < 2000 ? i+1 : 0)
    end
    for i = 2:2:2000
        @test findprev(odds,i)  == findprevnot(evens,i) == i-1
        @test findprev(evens,i) == findprevnot(odds,i)  == i
        @test findnext(evens,i) == findnextnot(odds,i)  == i
        @test findnext(odds,i)  == findnextnot(evens,i) == (i < 2000 ? i+1 : 0)
    end

    elts = (1:64:(64*64+1)) .+ (0:64)
    n = maximum(elts)
    for c = [falses, trues]
        b1 = c(n)
        b1[elts] = .!b1[elts]
        b2 = .~b1
        i1 = Array(b1)
        for i = 1:n
            @test findprev(b1, i) == findprev(i1, i) == findprevnot(b2, i) == findprev(!, b2, i)
            @test findnext(b1, i) == findnext(i1, i) == findnextnot(b2, i) == findnext(!, b2, i)
        end
    end

    b1 = falses(1000)
    b1[77] = true
    b1[777] = true
    b2 = .~b1
    @test_throws BoundsError findprev(b1, 1001)
    @test_throws BoundsError findprevnot(b2, 1001)
    @test_throws BoundsError findprev(!, b2, 1001)
    @test_throws BoundsError findprev(identity, b1, 1001)
    @test_throws BoundsError findprev(x->false, b1, 1001)
    @test_throws BoundsError findprev(x->true, b1, 1001)
    @test findprev(b1, 1000) == findprevnot(b2, 1000) == findprev(!, b2, 1000) == 777
    @test findprev(b1, 777)  == findprevnot(b2, 777)  == findprev(!, b2, 777)  == 777
    @test findprev(b1, 776)  == findprevnot(b2, 776)  == findprev(!, b2, 776)  == 77
    @test findprev(b1, 77)   == findprevnot(b2, 77)   == findprev(!, b2, 77)   == 77
    @test findprev(b1, 76)   == findprevnot(b2, 76)   == findprev(!, b2, 76)   == 0
    @test findprev(b1, -1)   == findprevnot(b2, -1)   == findprev(!, b2, -1)   == 0
    @test findprev(identity, b1, -1) == findprev(x->false, b1, -1) == findprev(x->true, b1, -1) == 0
    @test_throws BoundsError findnext(b1, -1)
    @test_throws BoundsError findnextnot(b2, -1)
    @test_throws BoundsError findnext(!, b2, -1)
    @test_throws BoundsError findnext(identity, b1, -1)
    @test_throws BoundsError findnext(x->false, b1, -1)
    @test_throws BoundsError findnext(x->true, b1, -1)
    @test findnext(b1, 1)    == findnextnot(b2, 1)    == findnext(!, b2, 1)    == 77
    @test findnext(b1, 77)   == findnextnot(b2, 77)   == findnext(!, b2, 77)   == 77
    @test findnext(b1, 78)   == findnextnot(b2, 78)   == findnext(!, b2, 78)   == 777
    @test findnext(b1, 777)  == findnextnot(b2, 777)  == findnext(!, b2, 777)  == 777
    @test findnext(b1, 778)  == findnextnot(b2, 778)  == findnext(!, b2, 778)  == 0
    @test findnext(b1, 1001) == findnextnot(b2, 1001) == findnext(!, b2, 1001) == 0
    @test findnext(identity, b1, 1001) == findnext(x->false, b1, 1001) == findnext(x->true, b1, 1001) == 0

    @test findlast(b1) == Base.findlastnot(b2) == 777
    @test findfirst(b1) == Base.findfirstnot(b2) == 77

    b0 = BitVector(0)
    @test findprev(x->true, b0, -1) == 0
    @test_throws BoundsError findprev(x->true, b0, 1)
    @test_throws BoundsError findnext(x->true, b0, -1)
    @test findnext(x->true, b0, 1) == 0

    b1 = falses(10)
    @test findprev(x->true, b1, 5) == 5
    @test findnext(x->true, b1, 5) == 5
    @test findprev(x->true, b1, -1) == 0
    @test findnext(x->true, b1, 11) == 0
    @test findprev(x->false, b1, 5) == 0
    @test findnext(x->false, b1, 5) == 0
    @test findprev(x->false, b1, -1) == 0
    @test findnext(x->false, b1, 11) == 0
    @test_throws BoundsError findprev(x->true, b1, 11)
    @test_throws BoundsError findnext(x->true, b1, -1)

    for l = [1, 63, 64, 65, 127, 128, 129]
        f = falses(l)
        t = trues(l)
        @test findprev(f, l) == findprevnot(t, l) == 0
        @test findprev(t, l) == findprevnot(f, l) == l
        b1 = falses(l)
        b1[end] = true
        b2 = .~b1
        @test findprev(b1, l) == findprevnot(b2, l) == l
        @test findprevnot(b1, l) == findprev(b2, l) == l-1
        if l > 1
            b1 = falses(l)
            b1[end-1] = true
            b2 = .~b1
            @test findprev(b1, l) == findprevnot(b2, l) == l-1
            @test findprevnot(b1, l) == findprev(b2, l) == l
        end
    end
end

@testset "Reductions" begin
    b1 = bitrand(s1, s2, s3, s4)
    m1 = 1
    m2 = 3
    @check_bit_operation maximum(b1, (m1, m2)) BitArray{4}
    @check_bit_operation minimum(b1, (m1, m2)) BitArray{4}
    @check_bit_operation sum(b1, (m1, m2)) Array{Int,4}

    @check_bit_operation maximum(b1) Bool
    @check_bit_operation minimum(b1) Bool
    @check_bit_operation any(b1) Bool
    @check_bit_operation all(b1) Bool
    @check_bit_operation sum(b1) Int

    b0 = falses(0)
    @check_bit_operation any(b0) Bool
    @check_bit_operation all(b0) Bool
    @check_bit_operation sum(b0) Int
end

timesofar("reductions")

@testset "map over bitarrays" begin
    for l = [0, 1, 63, 64, 65, 127, 128, 129, 255, 256, 257, 6399, 6400, 6401]
        b1 = bitrand(l)
        b2 = bitrand(l)
        @test map(~, b1) == map(x->~x, b1) == broadcast(~, b1)
        @test map(identity, b1) == map(x->x, b1) == b1

        @test map(&, b1, b2) == map((x,y)->x&y, b1, b2) == broadcast(&, b1, b2)
        @test map(|, b1, b2) == map((x,y)->x|y, b1, b2) == broadcast(|, b1, b2)
        @test map(⊻, b1, b2) == map((x,y)->x⊻y, b1, b2) == broadcast(⊻, b1, b2) == broadcast(xor, b1, b2)

        @test map(^, b1, b2) == map((x,y)->x^y, b1, b2) == b1 .^ b2
        @test map(*, b1, b2) == map((x,y)->x*y, b1, b2) == b1 .* b2

        @test map(min, b1, b2) == map((x,y)->min(x,y), b1, b2) == min.(b1, b2)
        @test map(max, b1, b2) == map((x,y)->max(x,y), b1, b2) == max.(b1, b2)

        @test map(<, b1, b2)  == map((x,y)->x<y, b1, b2)  == (b1 .< b2)
        @test map(<=, b1, b2) == map((x,y)->x<=y, b1, b2) == (b1 .<= b2)
        @test map(==, b1, b2) == map((x,y)->x==y, b1, b2) == (b1 .== b2)
        @test map(>=, b1, b2) == map((x,y)->x>=y, b1, b2) == (b1 .>= b2)
        @test map(>, b1, b2)  == map((x,y)->x>y, b1, b2)  == (b1 .> b2)
        @test map(!=, b1, b2) == map((x,y)->x!=y, b1, b2) == (b1 .!= b2)

        @testset "map! for length $l" begin
            b = BitArray(l)
            @test map!(~, b, b1) == map!(x->~x, b, b1) == broadcast(~, b1) == b
            @test map!(!, b, b1) == map!(x->!x, b, b1) == broadcast(~, b1) == b
            @test map!(identity, b, b1) == map!(x->x, b, b1) == b1 == b
            @test map!(zero, b, b1) == map!(x->false, b, b1) == falses(l) == b
            @test map!(one, b, b1) == map!(x->true, b, b1) == trues(l) == b

            @test map!(&, b, b1, b2) == map!((x,y)->x&y, b, b1, b2) == broadcast(&, b1, b2) == b
            @test map!(|, b, b1, b2) == map!((x,y)->x|y, b, b1, b2) == broadcast(|, b1, b2) == b
            @test map!(⊻, b, b1, b2) == map!((x,y)->x⊻y, b, b1, b2) == broadcast(⊻, b1, b2) == broadcast(xor, b1, b2) == b

            @test map!(^, b, b1, b2) == map!((x,y)->x^y, b, b1, b2) == b1 .^ b2 == b
            @test map!(*, b, b1, b2) == map!((x,y)->x*y, b, b1, b2) == b1 .* b2 == b

            @test map!(min, b, b1, b2) == map!((x,y)->min(x,y), b, b1, b2) == min.(b1, b2) == b
            @test map!(max, b, b1, b2) == map!((x,y)->max(x,y), b, b1, b2) == max.(b1, b2) == b

            @test map!(<, b, b1, b2)  == map!((x,y)->x<y, b, b1, b2)  == (b1 .< b2)  == b
            @test map!(<=, b, b1, b2) == map!((x,y)->x<=y, b, b1, b2) == (b1 .<= b2) == b
            @test map!(==, b, b1, b2) == map!((x,y)->x==y, b, b1, b2) == (b1 .== b2) == b
            @test map!(>=, b, b1, b2) == map!((x,y)->x>=y, b, b1, b2) == (b1 .>= b2) == b
            @test map!(>, b, b1, b2)  == map!((x,y)->x>y, b, b1, b2)  == (b1 .> b2)  == b
            @test map!(!=, b, b1, b2) == map!((x,y)->x!=y, b, b1, b2) == (b1 .!= b2) == b
        end
    end

    @testset "Issue #17970" begin
        A17970 = [1,2,3] .== [3,2,1]
        B17970 = map(x -> x ? 1 : 2, A17970)
        @test B17970::Array{Int,1} == [2,1,2]
        C17970 = map(x -> x ? false : true, A17970)
        @test C17970::BitArray{1} == map(~, A17970)
    end
end

## Filter ##

# TODO

@testset "Transpose" begin
    b1 = bitrand(v1)
    @check_bit_operation transpose(b1) RowVector{Bool,BitVector}

    for m1 = 0:n1, m2 = 0:n2
        b1 = bitrand(m1, m2)
        @check_bit_operation transpose(b1) BitMatrix
    end
end

timesofar("transpose")

@testset "Permutedims" begin
    b1 = bitrand(s1, s2, s3, s4)
    p = randperm(4)
    @check_bit_operation permutedims(b1, p) BitArray{4}
    @check_bit_operation permutedims(b1, tuple(p...)) BitArray{4}
end

timesofar("permutedims")

@testset "Concatenation" begin
    b1 = bitrand(v1)
    b2 = bitrand(v1)
    @check_bit_operation hcat(b1, b2) BitMatrix
    for m = 1:(v1-1)
        @check_bit_operation vcat(b1[1:m], b1[m+1:end]) BitVector
    end
    @test_throws DimensionMismatch hcat(b1,trues(n1+1))
    @test_throws DimensionMismatch hcat(hcat(b1, b2),trues(n1+1))

    b1 = bitrand(n1, n2)
    b2 = bitrand(n1)
    b3 = bitrand(n1, n2)
    b4 = bitrand(1, n2)
    @check_bit_operation hcat(b1, b2, b3) BitMatrix
    @check_bit_operation vcat(b1, b4, b3) BitMatrix
    @test_throws DimensionMismatch vcat(b1, b4, trues(n1,n2+1))

    b1 = bitrand(s1, s2, s3, s4)
    b2 = bitrand(s1, s3, s3, s4)
    b3 = bitrand(s1, s2, s3, s1)
    @check_bit_operation cat(2, b1, b2) BitArray{4}
    @check_bit_operation cat(4, b1, b3) BitArray{4}
    @check_bit_operation cat(6, b1, b1) BitArray{6}

    b1 = bitrand(1, v1, 1)
    @check_bit_operation cat(2, 0, b1, 1, 1, b1) Array{Int,3}
    @check_bit_operation cat(2, 3, b1, 4, 5, b1) Array{Int,3}
    @check_bit_operation cat(2, false, b1, true, true, b1) BitArray{3}

    b1 = bitrand(n1, n2)
    for m1 = 1:(n1-1), m2 = 1:(n2-1)
        @test isequal([b1[1:m1,1:m2] b1[1:m1,m2+1:end]; b1[m1+1:end,1:m2] b1[m1+1:end,m2+1:end]], b1)
    end
end

timesofar("cat")

@testset "Linear algebra" begin
    b1 = bitrand(v1)
    b2 = bitrand(v1)
    @check_bit_operation dot(b1, b2) Int

    b1 = bitrand(n1, n2)
    for k = (-n1):n2
        @check_bit_operation tril(b1, k) BitMatrix
        @check_bit_operation triu(b1, k) BitMatrix
    end

    for sz = [(n1,n1), (n1,n2), (n2,n1)], (f,isf) = [(tril,istril), (triu,istriu)]
        b1 = bitrand(sz...)
        @check_bit_operation isf(b1) Bool
        b1 = f(bitrand(sz...))
        @check_bit_operation isf(b1) Bool
    end

    b1 = bitrand(n1,n1)
    b1 .|= b1.'
    @check_bit_operation issymmetric(b1) Bool
    @check_bit_operation ishermitian(b1) Bool

    b1 = bitrand(n1)
    b2 = bitrand(n2)
    @check_bit_operation kron(b1, b2) BitVector

    b1 = bitrand(s1, s2)
    b2 = bitrand(s3, s4)
    @check_bit_operation kron(b1, b2) BitMatrix

    b1 = bitrand(v1)
    @check_bit_operation diff(b1) Vector{Int}

    b1 = bitrand(n1, n2)
    @check_bit_operation diff(b1) Matrix{Int}

    b1 = bitrand(n1, n1)
    @check_bit_operation svd(b1)
    @check_bit_operation qr(b1)

    b1 = bitrand(v1)
    @check_bit_operation gradient(b1)
    @check_bit_operation gradient(b1, 1.0)

    b1 = bitrand(v1)
    @check_bit_operation diagm(b1) BitMatrix

    b1 = bitrand(n1, n2)
    @test_throws DimensionMismatch diagm(b1)

    b1 = bitrand(n1, n1)
    @check_bit_operation diag(b1)
end

timesofar("linalg")

@testset "findmax, findmin" begin
    b1 = trues(0)
    @test_throws ArgumentError findmax(b1)
    @test_throws ArgumentError findmin(b1)

    for b1 in [falses(v1), trues(v1),
               BitArray([1,0,1,1,0]),
               BitArray([0,0,1,1,0]),
               bitrand(v1)]
        @check_bit_operation findmin(b1)
        @check_bit_operation findmax(b1)
    end
end

@testset "I/O" begin
    b1 = bitrand(v1)
    fname = ""
    b1[v1 ÷ 2 + 1] = true
    b1[end] = true
    try
        fname = tempname()
        open(fname, "w") do f
            write(f, b1)
        end
        b2 = falses(v1)
        read!(fname, b2)
        @test bitcheck(b2)
        @test b1 == b2
        b2 = falses(v1 ÷ 10, 10)
        read!(fname, b2)
        @test bitcheck(b2)
        @test reshape(b1, v1 ÷ 10, 10) == b2
        b2 = falses(v1 + 65)
        @test bitcheck(b2)
        @test_throws EOFError read!(fname, b2)
        @test bitcheck(b2)
        b2 = falses(v1 ÷ 2)
        @test_throws DimensionMismatch read!(fname, b2)
        @test bitcheck(b2)
        b2 = falses(v1 - 1)
        @test_throws DimensionMismatch read!(fname, b2)
        @test bitcheck(b2)

        b1 = BitArray(0)
        open(fname, "w") do f
            write(f, b1)
        end
        b2 = BitArray(0)
        read!(fname, b2)
        @test b1 == b2
        @test bitcheck(b2)
        b2 = trues(1)
        @test_throws EOFError read!(fname, b2)
        @test bitcheck(b2)
    finally
         isfile(fname) && rm(fname)
    end
end

timesofar("I/O")

end # module

max_stars_repo_name: jamesonquinn/julia

max_issues_repo_name: jamesonquinn/julia

max_forks_repo_name: jamesonquinn/julia


text:
using Mimi

y = macroexpand(:(
    @defcomp foo begin
        index1 = Index()

        par1 = Parameter()
        par2::Bool = Parameter(index=[time,index1], description="description par 1")
        par3 = Parameter(index=[time])

        var1 = Variable()
        var2 = Variable(index=[time])
        var3 = Variable(index=[time,index1])
        var4::Bool = Variable(index=[3])
        var5 = Variable(index=[index1,4])
    end
    ))

println(y)


text:
using Gridap
using Test
n = 10
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(n,n))
labels = FaceLabels(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])
D = 2
order = 2
fespace1 = FESpace(
  reffe=:Lagrangian, conformity=:H1, valuetype=VectorValue{D,Float64},
  model=model, labels=labels, order=order, diritags=["diri0","diri1"])
fespace2 = FESpace(
  reffe=:PLagrangian, conformity=:L2, valuetype=Float64,
  model=model, order=order-1, constraint=:zeromean)
uD0(x) = VectorValue(0.0,0.0)
uD1(x) = VectorValue(1.0,0.0)

V = TestFESpace(fespace1)
Q = TestFESpace(fespace2)
Y = [V, Q]

U = TrialFESpace(fespace1,[uD0,uD1])
P = TrialFESpace(fespace2)
X = [U, P]
const Re = 10.0
@law conv(x,u,∇u) = Re*(∇u')*u
@law dconv(x,du,∇du,u,∇u) = conv(x,u,∇du)+conv(x,du,∇u)

function a(y,x)
  u, p = x
  v, q = y
  inner(∇(v),∇(u)) - inner(div(v),p) + inner(q,div(u))
end

c(v,u) = inner(v,conv(u,∇(u)))
dc(v,du,u) = inner(v,dconv(du,∇(du),u,∇(u)))

function res(x,y)
  u, p = x
  v, q = y
  a(y,x) + c(v,u)
end

function jac(x,y,dx)
  u, p = x
  v, q = y
  du, dp = dx
  a(y,dx)+ dc(v,du,u)
end
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=(order-1)*2)
t_Ω = NonLinearFETerm(res,jac,trian,quad)
op = NonLinearFEOperator(Y,X,t_Ω)
using LineSearches: BackTracking
nls = NLSolver(
  show_trace=false, method=:newton, linesearch=BackTracking())
solver = NonLinearFESolver(nls)
uh, ph = solve(solver,op)
#writevtk(trian,"ins-results",cellfields=["uh"=>uh,"ph"=>ph])

meanval = sum( integrate(ph,trian,quad))
@test meanval + 1 ≈ 1


text:
module TSMLextra

using Reexport

@reexport using TSML

include("system.jl")
using .System

include("datareader.jl")
using .DataReaders

include("datawriter.jl")
using .DataWriters

if LIB_CRT_AVAILABLE # from System module
    include("caret.jl")
    using .CaretLearners
end

if LIB_SKL_AVAILABLE # from System module
    include("scikitlearn.jl")
    using .SKLearners
end

end # module


text:
# use https://jacobadenbaum.github.io/TexTables.jl
# use https://github.com/magerton/TexTables.jl
using Revise
using ShaleDrillingLikelihood

TEX_TABLE_PATHS = [
    "E:/projects/royalty-rates-and-drilling/writeup/JMPEdited/tables/results-2020-01.tex",
    "E:/projects/ShaleDrillingResults/results-2020-01.tex"
]
MAIN_TEX_FILE = "../JMPEdited.tex"

FIRST_COST_YEAR = 2009
LAST_COST_YEAR = 2016

scriptdir = joinpath(dirname(pathof(ShaleDrillingLikelihood)), "..", "example")
include(joinpath(scriptdir, "display-coefs-fcts.jl"))

# JLD2FILE = "E:/projects/ShaleDrillingResults/2019-12/16481962/estimation-results-16481962-no-rigs-data_all_leases.jld2"

title_dataset = [
     # ("All"          , "E:/projects/ShaleDrillingResults/2019-12/16481962/estimation-results-16481962-no-rigs-data_all_leases.jld2",                  ),
     # ("First"        , "E:/projects/ShaleDrillingResults/2019-12/16627219/estimation-results-16627219-no-rigs-data_first_lease_only.jld2",            ),
     # ("First, restr" , "E:/projects/ShaleDrillingResults/2019-12/16482515/estimation-results-16482515-no-rigs-data_first_lease_only_all_leased.jld2", ),
     # ("Last"         , "E:/projects/ShaleDrillingResults/2019-12/16482550/estimation-results-16482550-no-rigs-data_last_lease_only.jld2",             ),
     # ("All, rigs"    , "E:/projects/ShaleDrillingResults/2019-12/16486124/estimation-results-16486124-WITH-rigs-data_all_leases.jld2",                ),
     # ("T1EV"         , "E:/projects/ShaleDrillingResults/2019-12/16627762/estimation-results-16627762-no-rigs-data_all_leases.jld2",                  ),
     ("Baseline", "E:/projects/ShaleDrillingResults/2020-01/17529657/estimation-results-17529657-no-rigs-data_all_leases.jld2", ),
     ("2008.25"  , "E:/projects/ShaleDrillingResults/2020-01/17608936/estimation-results-17608936-no-rigs-data_all_leases.jld2", ),
     ("2008.5 "  , "E:/projects/ShaleDrillingResults/2020-01/17493679/estimation-results-17493679-no-rigs-data_all_leases.jld2", ),
     ("2008.75"  , "E:/projects/ShaleDrillingResults/2020-01/17525858/estimation-results-17525858-no-rigs-data_all_leases.jld2", ),
     ("First"    , "E:/projects/ShaleDrillingResults/2020-01/17608948/estimation-results-17608948-no-rigs-data_first_lease_only.jld2", ),
     ("Restr"  , "E:/projects/ShaleDrillingResults/2020-01/17608940/estimation-results-17608940-no-rigs-data_first_lease_only_all_leased.jld2", ),
     ("Last"     , "E:/projects/ShaleDrillingResults/2020-01/17608951/estimation-results-17608951-no-rigs-data_last_lease_only.jld2", ),
     ("T1EV"     , "E:/projects/ShaleDrillingResults/2020-01/17608965/estimation-results-17608965-no-rigs-data_all_leases.jld2", ),
 ]

rcs = [regcol(t, nms_coef_se_sumstat(j)...) for (t,j) in title_dataset]

# idx_rho = 12
# idxs = vcat(13:21,1:12,22:25)
# thetarho = coef[idx_rho]
# se[idx_rho] *= _dρdθρ(thetarho)
# coef[idx_rho] = _ρ(thetarho)
# coef[idxs], se[idxs]

mynote = "Mean well costs are the mean per-well drilling cost over the period $FIRST_COST_YEAR--$LAST_COST_YEAR. " *
    "This is calculated as \$\\frac{\\sigma_\\epsilon}{1-\\tau_{k}} \\frac{cost(2,s_{it},z_{it})}{2}\$ " *
    "where \$\\sigma_\\epsilon\$ is computed from \\eqref{eq:sigma-epsilon}. " *
    "The effective marginal corporate income tax is $(100*MGL_CORP_INC_TAX)\\%, " *
    "and the marginal tax rate on capital investment is \$\\tau_{k} \\approx $(round(100*CAPITAL_MGL_TAX_RATE;digits=1))\\% \$"

tblcols = regtable(rcs...)
mt = make_table(tblcols, "Estimates for full model", "tab:all-estimates", mynote)
println(mt)

for t in TEX_TABLE_PATHS
    io = open(t, "w")
        write(io, "%!TEX root = " * MAIN_TEX_FILE * "\n\n\n\n")
        write(io, mt); write(io, "\n\n\n\n\n")
    close(io)
end


text:
# ---
# title: 672. Bulb Switcher II
# id: problem672
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/bulb-switcher-ii/description/>
# hidden: true
# ---
#
# There is a room with `n` lights which are turned on initially and 4 buttons on
# the wall. After performing exactly `m` unknown operations towards buttons, you
# need to return how many different kinds of status of the `n` lights could be.
#
# Suppose `n` lights are labeled as number [1, 2, 3 ..., n], function of these 4
# buttons are given below:
#
#   1. Flip all the lights.
#   2. Flip lights with even numbers.
#   3. Flip lights with odd numbers.
#   4. Flip lights with (3k + 1) numbers, k = 0, 1, 2, ...
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1, m = 1.
#     Output: 2
#     Explanation: Status can be: [on], [off]
#
#
#
#
# **Example 2:**
#
#
#
#     Input: n = 2, m = 1.
#     Output: 3
#     Explanation: Status can be: [on, off], [off, on], [off, off]
#
#
#
#
# **Example 3:**
#
#
#
#     Input: n = 3, m = 1.
#     Output: 4
#     Explanation: Status can be: [off, on, off], [on, off, on], [off, off, off], [off, on, on].
#
#
#
#
# **Note:** `n` and `m` both fit in range [0, 1000].
#
#
## @lc code=start
using LeetCode

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


text:
module astar

DIAGONALS = true
HEURISTIC = 1

mutable struct Node
    pos::Int64
    goal_position::Int64
    fCost::Float64
    gCost::Float64
    hCost::Float64
    isClosed::Bool
    isOpen::Bool
    parent
end

mutable struct Path
    pos::Array{Int64}
    costs::Array{Float64}
end

mutable struct Environment
    length::Int64
    height::Int64
    wallPercentage::Float64
    grid::Array{Int8}
end

function get_actions(env)
    actions = [1 1;
               env.length 1;
              -1 1;
              -env.length 1]

    if DIAGONALS
        d = [-env.length-1 sqrt(2);
             -env.length+1 sqrt(2);
              env.length-1 sqrt(2);
              env.length+1 sqrt(2)]
        actions = [actions;d]
    end
    return actions
end

function print_color(color, string)
    if color == "red"
        print("\033[1;31m")
        print(string)
        print("\033[0m")
    elseif color == "blue"
        print("\033[1;34m")
        print(string)
        print("\033[0m")
    elseif color == "yellow"
        print("\033[1;33m")
        print(string)
        print("\033[0m")
    elseif color == "green"
        print("\033[1;32m")
        print(string)
        print("\033[0m")
    elseif color == "cyan"
        print("\033[1;36m")
        print(string)
        print("\033[0m")
    else
        print(string)
    end
end

function generate_environment(length, height, wallPercentage)
    gd = zeros(Int8, height, length)

    for (index, _) in enumerate(gd)
        if index % length == 0 || index <= length || index > length*height - length || (index - 1) % length == 0
            gd[index] = 1
            continue
        end
        if index//height % length == 0
            gd[index] = 1
            continue
        end
        if rand() < wallPercentage
            gd[index] = 1
        end

    end

    env = Environment(length, height, wallPercentage, gd)
    return env
end

function get_random_pos(env)
    size = env.length * env.height
    n = rand(1:size)
    while env.grid[n] == 1
        n = rand(1:size)
    end
    return n
end

function get_start_goal(env, p_max_length)
    if p_max_length > 0.95
        p_max_length = 0.95
    end

    x = 0

    start = get_random_pos(env)
    goal = get_random_pos(env)
    size = env.length * env.height
    max_h = calculate_heuristic(env.length+2, size - env.length - 1, env)

    while calculate_heuristic(start, goal, env) <= max_h*p_max_length
        start = get_random_pos(env)
        goal = get_random_pos(env)
        x += 1
        if x > 1000
            print_color("red", "Could not find start and goal respecting constraints...\n")
            break
        end
    end

    return start, goal
end

function print_grid(env)
    for (index, val) in enumerate(env.grid)
        if index % env.length != 0
            if val == 1
                print_color("red", "█")
            elseif val == 0
                print(" ")
            elseif val == 2
                print_color("blue", "¤")
            elseif val == 3
                print_color("green", "x")
            elseif val == 4
                print_color("cyan", "o")
            end
        else
            if val == 1
                print_color("red", "█")
            elseif val == 0
                print(" ")
            elseif val == 2
                print_color("blue", "¤")
            elseif val == 3
                print_color("green", "x")
            elseif val == 4
                print_color("cyan", "o")
            end
            print("\n")
        end
    end
end

function add_elem_to_env(env, pos, val)
    env.grid[pos] = val
    return env
end

function add_path(env, path)
    for (index, val) in enumerate(path.pos)
        if val == path.pos[1] || val == last(path.pos)
            continue
        end
        env = add_elem_to_env(env, val, 2)
    end

    return env
end

function calculate_heuristic(current_pos, goal_pos, env)
    row_current = convert(Int64, floor(current_pos / env.length)) + 1
    col_current = current_pos % env.length
    row_target = convert(Int64, floor(goal_pos / env.length)) + 1
    col_target = goal_pos % env.length

    dy = abs(row_target - row_current)
    dx = abs(col_target - col_current)

    if HEURISTIC == 0
        return(sqrt(dy^2 + dx^2))
    elseif HEURISTIC == 1
        return(dx + dy)
    else
        print_color("red", "Reached unreachable statement.")
        exit()
    end
end

function find_lowest_cost_node(nodes)
    x = 0
    currentNode = nothing
    for (_, val) in nodes
        if val.isOpen
            x += 1
            currentNode = val
            break
        end
    end

    if !(x > 0)
        return 0, nothing
    end

    for (_, val) in nodes
        if val.isOpen && val != currentNode
            x += 1
            if val.fCost < currentNode.fCost
                currentNode = val
            end
        end
    end

    return x, currentNode
end

function get_neighbours(current_node, goal_pos, env)
    c_pos = current_node.pos
    size = env.length * env.height

    children = []

    if c_pos < 0 || c_pos >= env.length*env.height || env.grid[c_pos] == 1
        print_color("red", "Tried to get neighbours from unreachable tile, exiting...")
        exit()
    end

    for row in eachrow(get_actions(env))
        if !(c_pos + row[1] > 0 && c_pos + row[1] <= size)
            continue
        end

        if env.grid[Int(c_pos + row[1])] == 1
            continue
        end

        child = node_constructor(Int(c_pos + row[1]),
                                 goal_pos,
                                 current_node.gCost + row[2],
                                 false, false, current_node, env)

        children = push!(children, child)
    end

    return children
end

function node_constructor(position, goal_position, gcost, isclosed, isopen, parent, env)
    h = calculate_heuristic(position, goal_position, env)
    return(Node(position, goal_position, gcost + h, gcost, h, isclosed, isopen, parent))
end

function retrieve_path(last_node)
    current_node = last_node
    path = Path([], [])
    while current_node !== nothing
        path.pos = push!(path.pos, current_node.pos)
        path.costs = push!(path.costs, current_node.gCost)
        current_node = current_node.parent
    end
    return path
end

function search(start, goal, env)
    nodes = Dict{Int64, Node}()
    start_node = node_constructor(start, goal, 0, false, true, nothing, env)
    nodes[start] = start_node

    while true
        open_nodes_cnt, currentNode = find_lowest_cost_node(nodes)

        if open_nodes_cnt < 1
            print_color("red", "The graph has no open nodes left...\n")
            break
        end

        if currentNode.pos == goal
            print_color("yellow", "Path found!\n")
            return(retrieve_path(currentNode))
        end

        nodes[currentNode.pos].isOpen = false
        nodes[currentNode.pos].isClosed = true

        children = get_neighbours(currentNode, goal, env)

        for (_, child) in enumerate(children)
            if haskey(nodes, child.pos)
                if nodes[child.pos].isClosed
                    continue
                end

                if !nodes[child.pos].isOpen
                    child.isOpen = true
                    nodes[child.pos] = child
                else
                    if child.gCost < nodes[child.pos].gCost
                        child.isOpen = true
                        nodes[child.pos] = child
                    end
                end
            else
                child.isOpen = true
                nodes[child.pos] = child
            end
        end
    end
end

function main()
    println(PROGRAM_FILE, " starting...")
    length = 50
    height = 20
    wallP = .25
    env = generate_environment(length, height, wallP)

    start, goal = get_start_goal(env, 0.5)

    env = add_elem_to_env(env, start, 3)
    env = add_elem_to_env(env, goal, 3)

    print_grid(env)

    path = search(start, goal, env)

    if path !== nothing
        env = add_path(env, path)
        print_grid(env)
    end
end

main()

end #astar

max_stars_repo_name: guilyx/AStar

max_issues_repo_name: guilyx/AStar

max_forks_repo_name: guilyx/AStar


text:
module RTUtils

using ...Scruff
using ..Utils
using ..SFuncs

import ..Operators: Support

using ..MultiInterface

include("utils/range_utils.jl")
include("utils/expansions.jl")
include("utils/initializer.jl")

end


text:
using BinaryBuilder

# Collection of sources required to build Darknet
name = "Darknet"
version = v"2019.11.15"
sources = [
    "https://github.com/AlexeyAB/darknet/archive/71e835458904f782a905a06d28b4558d9e9830b4.zip" =>
    "d77017462ae49f9ce2540c3e47589e68b2ce565573bf7d2f011b560bc989fbfa",
    "./bundled",
]

script = raw"""
cd $WORKSPACE/srcdir/darknet-*

if [[ "${target}" == *-mingw* ]]; then
    # Fix case of some Windows headers
    atomic_patch -p1 ../patches/windows_headers_case.patch
    # Comment out a couple of definitions
    atomic_patch -p1 ../patches/getopt_windows.patch
    # Link against ws2_32 on Windows
    atomic_patch -p1 ../patches/windows_ldflags.patch
    # This makes it work for 64-bit Windows
    atomic_patch -p1 ../patches/gemmc_windows_64bit.patch
fi

# if [[ "${target}" = powerpc64le-* ]] || [[ "${target}" = arm* ]] || [[ "${target}" == aarch* ]] || [[ "${target}" == *-mingw* ]] || [[ "${target}" == *-apple* ]]; then
#     # Disable AVX on powerpc, arm, aarch, windows, apple
#     export AVXENABLE=0
# else
#     # Enable everywhere else (linux)
#     export AVXENABLE=1
# fi
export AVXENABLE=0

# Make sure to have the directories, before building
make obj backup results setchmod
make -j${nproc} libdarknet.${dlext} \
    LIBNAMESO="libdarknet.${dlext}" \
    LIBSO=1 \
    GPU=0 \
    CUDNN=0 \
    CUDNN_HALF=0 \
    OPENCV=0 \
    DEBUG=0 \
    OPENMP=0 \
    ZED_CAMERA=0 \
    AVX=${AVXENABLE}

mkdir -p "${libdir}"
cp libdarknet.${dlext} "${libdir}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libdarknet", :libdarknet),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
"""
A package to allow serializing and deserializing data from a Pluto notebook.
Provides `serialize` and `deserialize` functions which follow the same
interface as those provided by the `Serialization`, but do not extend it.
Importing both will result in name conflicts, use one or the other, not both.
"""
module PlutoSerialization

# Imports:

import Serialization

# Exports:

export deserialize
export serialize

# Public interface:

"""
    serialize(filename::AbstractString, object)
    serialize(io::IO, object)

Save the given `object` to `filename` or `io` using Julia's serialization, but
additionally handle Pluto's workspace modules such that the saved data can be
`deserialized` in a new session.
"""
function serialize end

function serialize(io::IO, @nospecialize(object))
    ps = PlutoSerializer(io)
    Serialization.writeheader(ps)
    Serialization.serialize(ps, object)
end
serialize(fn::AbstractString, @nospecialize(object)) = open(io -> serialize(io, object), fn, "w")

"""
    deserialize(filename::AbstractString)
    deserialize(io::IO)

Deserialize `filename` or `io` into a Julia object using Julia's serialization,
but correctly handle Pluto's workspace modules.
"""
function deserialize end

deserialize(io::IO) = Serialization.deserialize(PlutoSerializer(io))
deserialize(fn::AbstractString) = open(deserialize, fn)

# A mock top-level module to mark module references serialized in notebooks.
baremodule PlutoWorkspaces end
_parentmodule(m::Module) = m === PlutoWorkspaces ? PlutoWorkspaces : parentmodule(m)
_root_module_key(m::Module) = m === PlutoWorkspaces ? Base.PkgId("PlutoWorkspaces") : Base.root_module_key(m)
_root_module(pkgid::Base.PkgId) = Base.PkgId("PlutoWorkspaces") === pkgid ? PlutoWorkspaces : Base.root_module(pkgid)

# Implementation:

# Serialization and deserialization code for handling pluto notebook
# environments correctly. Notebooks incrementally add `workspaceN` modules when
# they re-evaluate cells, but when restarting a notebook `N` begins again from
# `1` rather than preserving the previous state (nothing wrong with this), but
# it means that using the default `serialize` and `deserialize` will likely
# store and load `workspaceN` modules that probably don't exist.
#
# To handle this, we implement an `AbstractSerializer` that checks through all
# available workspaces when deserializing objects to find the most recent
# reference to the object name, rather than expecting the specific named
# workspace to exist.

"""
Find all Pluto workspace modules and return an `id => mod` mapping for each.
"""
function get_ws()
    dict = Dict{Int,Module}()
    for name in names(Main)
        id = ws_number(name)
        isnothing(id) && continue
        if isdefined(Main, name)
            obj = getfield(Main, name)
            if isa(obj, Module) && isdefined(obj, :PlutoRunner)
                dict[id] = obj
            end
        end
    end
    return dict
end

"""
The `Int` representing a Pluto workspace, or `nothing` if not valid.
"""
function ws_number(name::Symbol)
    str = String(name)
    m = match(r"^workspace[#]?(\d+)$", str)
    return isnothing(m) ? nothing : tryparse(Int, m[1])
end
ws_number(m::Module) = ws_number(nameof(m))

"""
Is the `Module` a valid Pluto workspace?
"""
is_workspace(m::Module) = !isnothing(ws_number(m)) && isdefined(m, :PlutoRunner)

"""
Check through workspace modules from most recent and find the first module in
which the `name` is defined, that's *probably* the right one... otherwise just
give up and return the `default` module instead.
"""
function maybe_defined_in_workspace(workspaces::Dict, name::Symbol, default::Module)
    for id in length(workspaces):-1:1
        if haskey(workspaces, id)
            mod = workspaces[id]
            isdefined(mod, name) && return mod
        end
    end
    return isdefined(Main, name) ? Main : default
end

# Serialization implementation:
#
# Adapts some of the methods of `AbstractSerializer` to handle searching
# through all available workspace modules to find valid deserialized objects,
# but for the most part simply reuses the default `Serialization` methods.

mutable struct PlutoSerializer{I<:IO} <: Serialization.AbstractSerializer
    io::I
    counter::Int
    table::IdDict{Any,Any}
    pending_refs::Vector{Int}
    known_object_data::Dict{UInt64,Any}
    version::Int
    workspaces::Dict{Int,Module}
    PlutoSerializer{I}(io::I) where I<:IO = new(io, 0, IdDict(), Int[], Dict{UInt64,Any}(), Serialization.ser_version, get_ws())
end

PlutoSerializer(io::IO) = PlutoSerializer{typeof(io)}(io)

function Serialization.serialize_mod_names(s::PlutoSerializer, m::Module)
    # Swap out the workspace module for a custom toplevel one that we control
    # to "mark" it as a workspace.
    m = haskey(s.workspaces, ws_number(m)) ? PlutoWorkspaces : m
    p = _parentmodule(m)
    if p === m || m === Base
        key = _root_module_key(m)
        Serialization.serialize(s, key.uuid === nothing ? nothing : key.uuid.value)
        Serialization.serialize(s, Symbol(key.name))
    else
        Serialization.serialize_mod_names(s, p)
        Serialization.serialize(s, nameof(m))
    end
end

function Serialization.deserialize_module(s::PlutoSerializer)
    mkey = Serialization.deserialize(s)
    if isa(mkey, Tuple)
        # old version, TODO: remove
        if mkey === ()
            return Main
        end
        m = _root_module(mkey[1])
        for i = 2:length(mkey)
            m = getfield(m, mkey[i])::Module
        end
    else
        name = String(Serialization.deserialize(s)::Symbol)
        pkg = (mkey === nothing) ? Base.PkgId(name) : Base.PkgId(Base.UUID(mkey), name)
        m = _root_module(pkg)
        mname = Serialization.deserialize(s)
        while mname !== ()
            if m === PlutoWorkspaces || is_workspace(m)
                m = maybe_defined_in_workspace(s.workspaces, mname, m)
            end
            m = getfield(m, mname)::Module
            mname = Serialization.deserialize(s)
        end
    end
    return m === PlutoWorkspaces ? get(s.workspaces, length(s.workspaces), Main) : m
end

function Serialization.deserialize_datatype(s::PlutoSerializer, full::Bool)
    slot = s.counter
    s.counter += 1
    if full
        tname = Serialization.deserialize(s)::Core.TypeName
        ty = tname.wrapper
    else
        name = Serialization.deserialize(s)::Symbol
        mod = Serialization.deserialize(s)::Module
        if mod === PlutoWorkspaces || is_workspace(mod)
            mod = maybe_defined_in_workspace(s.workspaces, name, mod)
        end
        ty = getfield(mod, name)
    end
    if isa(ty,DataType) && isempty(ty.parameters)
        t = ty
    else
        np = Int(read(s.io, Int32)::Int32)
        if np == 0
            t = Base.unwrap_unionall(ty)
        elseif ty === Tuple
            # note np==0 has its own tag
            if np == 1
                t = Tuple{Serialization.deserialize(s)}
            elseif np == 2
                t = Tuple{Serialization.deserialize(s), Serialization.deserialize(s)}
            elseif np == 3
                t = Tuple{Serialization.deserialize(s), Serialization.deserialize(s), Serialization.deserialize(s)}
            elseif np == 4
                t = Tuple{Serialization.deserialize(s), Serialization.deserialize(s), Serialization.deserialize(s), Serialization.deserialize(s)}
            else
                t = Tuple{Any[Serialization.deserialize(s) for _ = 1:np]...}
            end
        else
            t = ty
            for _ = 1:np
                t = t{Serialization.deserialize(s)}
            end
        end
    end
    s.table[slot] = t
    return t
end

# Might need to adjust, but the idea here is that *all* the required
# definitions should already exist within the notebook, we shouldn't be
# deserializing anything that doesn't already have the whole type available.
# May need to re-evaluate this choice.
Serialization.should_send_whole_type(s::PlutoSerializer, t::DataType) = false

end


text:
@with_kw mutable struct MomentumParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    μ::T = 0.9
    ϵ::T = 1e-3
end

abstract type MomentumTrait end
abstract type Classical <: MomentumTrait end
abstract type Nesterov <: MomentumTrait end

struct Momentum{T <: MomentumTrait} <: AbstractBoosting
    params::MomentumParameters{Float64}
    ν::Vector{Float64}

    function (::Type{Momentum})(::Type{T}; kwargs...) where T <: MomentumTrait
        return new{T}(MomentumParameters(; kwargs...),Vector{Float64}())
    end
end

function initialize!(momentum::Momentum,x₀::Vector{Float64})
    resize!(momentum.ν,length(x₀))
    momentum.ν .= zeros(length(x₀))
end

params(momentum::Momentum) = momentum.params

function boost!(momentum::Momentum{Classical},wid::Integer,klocal::Integer,kglobal::Integer,gprev::AbstractVector,gcurr::AbstractVector)
    @unpack μ,ϵ = momentum.params
    momentum.ν .= μ*momentum.ν + ϵ*gprev
    gcurr .= momentum.ν
end

function boost!(momentum::Momentum{Nesterov},wid::Integer,klocal::Integer,kglobal::Integer,gprev::AbstractVector,gcurr::AbstractVector)
    @unpack μ,ϵ = momentum.params
    νprev = copy(momentum.ν)
    momentum.ν .= μ*νprev+ϵ*gprev
    gcurr .= μ^2*νprev+(1+μ)*ϵ*gprev
end

Momentum(; kwargs...) = Momentum(Classical; kwargs...)
Nesterov(; kwargs...) = Momentum(Nesterov; kwargs...)


text:
using BinaryBuilder, Pkg

name = "LibGD"
version = v"2.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/libgd/libgd/releases/download/gd-$(version)/libgd-$(version).tar.gz",
                  "32590e361a1ea6c93915d2448ab0041792c11bae7b18ee812514fe08b2c6a342")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libgd-*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
# For some reasons (something must be off in the configure script), on some
# platforms the build system tries to use iconv but without adding the `-liconv`
# flag.  Give a hint to make to use the right flag everywhere
make -j${nproc} LIBICONV="-liconv" LTLIBICONV="-liconv"
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("pngtogd2", :pngtogd2),
    ExecutableProduct("webpng", :webpng),
    ExecutableProduct("pngtogd", :pngtogd),
    ExecutableProduct("gdtopng", :gdtopng),
    ExecutableProduct("gdcmpgif", :gdcmpgif),
    ExecutableProduct("gd2topng", :gd2topng),
    ExecutableProduct("gdparttopng", :gdparttopng),
    ExecutableProduct("gd2copypal", :gd2copypal),
    ExecutableProduct("gd2togif", :gd2togif),
    ExecutableProduct("giftogd2", :giftogd2),
    LibraryProduct("libgd", :libgd),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8")),
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a")),
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f")),
    Dependency(PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828")),
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid = "c4d99508-4286-5418-9131-c86396af500b")),
    Dependency(PackageSpec(name="Xorg_libXpm_jll", uuid = "1a3ddb2d-74e3-57f3-a27b-e9b16291b4f2")),
    Dependency(PackageSpec(name="Libiconv_jll", uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531")),
    Dependency(PackageSpec(name="libwebp_jll", uuid = "c5f90fcd-3b7e-5836-afba-fc50a0988cb2")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# using CSV, DataFrames, Interact, Plots
#
#
# loadbutton = filepicker();
# columnbuttons = Observable{Any}(dom"div"());
# data = Observable{Any}(DataFrame);
# plt = Observable{Any}(plot());
# map!(CSV.read, data, loadbutton);
#
# function makebuttons(df)
#     buttons = button.(string.(names(df)))
#     for (btn, name) in zip(buttons, names(df))
#         map!(t -> histogram(df[name]), plt, btn)
#     end
#     dom"div"(hbox(buttons))
# end
#
# map!(makebuttons, columnbuttons, data);
#
# ui = dom"div"(loadbutton, columnbuttons, plt);
#
# using Blink
# w = Window()
# body!(w, ui)
#
#
# # function serverLauout(destinationPort)
# #     try
# #         WebIO.webio_serve(page("/",req -> creatGui()),destinationPort)
# #     catch e
# #         if isa(e,IOError)
# #             # sleep and then try later
# #             sleep(0.1)
# #             serverLauout(destinationPort)
# #         else
# #             throw(e)
# #         end
# #     end
# # end
# #
# # serverLauout(8000)

using Interact, DifferentialEquations, Latexify, Plots, WebIO
gr()

"""
A function which simulates and plots an ODE.
The system is first equilibrated for an input value of 1.
The plotted dynamics is the relaxation of the system after the input is changed to 10.
"""
function plotODE(ode, parameters, plotvars, initcond)
    tspan = (0.,100.)
    prob = ODEProblem(ode, initcond, tspan, parameters)
    sol = solve(prob)

    plot(
    plot(sol, vars=(:v, :w), xlabel="v", ylabel="w", title="Phase plot"), # make the phase plot
    plot(sol, vars=(:t, :v, :w), legend=:none, title="Phase(t)"),   # plot phase vs time
    plot(sol, title="FHN solution", xlabel="Time"), # plot the solution versus time
    layout=@layout [a b c]
    )

end

"""
Automatically generate sliders for all the ODE's parameters and map the results
to the plotODE function.
"""
function interactivePlot(ode, initcond)
    display(latexalign(ode))
    params = [slider(round.(exp10.(range(-1, stop=2, length=101)), sigdigits=3), label=latexify(p)) for p in ode.params]
    plotvars = ode.syms
    display(hbox(vbox(params...)))
    map((x...)->plotODE(ode, collect(x[1:end-1]), x[end], initcond),params..., plotvars)
end

"""
Define some ODEs
"""
ode1 = @ode_def NegativeFeedback begin
    dx = r_x * (e_x * input * y - x)
    dy = r_y * (e_y / x - y)
end input r_x e_x r_y e_y

ode2 = @ode_def IncoherentFeedForward begin
    dx = r_x * (e_x * input - x)
    dy = r_y * (e_y * input / x - y)
end input r_x e_x r_y e_y

ode3 = @reaction_network InducedDegradation begin
    (input*r_bind, r_unbind), X_free ↔ X_bound
    (p_free, d_free), 0 ↔ X_free
    d_bound, X_bound --> 0
end input r_bind r_unbind p_free d_free d_bound

max_stars_repo_name: kobied/interactTests

max_issues_repo_name: kobied/interactTests

max_forks_repo_name: kobied/interactTests


text:
using HDF5
using ProgressMeter

function embeds2h5(path::String)
    keys = String[]
    values = Float32[]
    lines = open(readlines, path)
    prog = Progress(length(lines))
    for line in lines
        items = split(line, " ")
        push!(keys, items[1])
        for i = 2:length(items)
            push!(values, parse(Float32,items[i]))
        end
        next!(prog)
    end
    value = reshape(values, length(values)÷length(keys), length(keys))

    outpath = path * ".h5"
    info("Key size: $(length(keys))")
    info("Value size: $(size(value))")
    h5write(outpath, "key", keys)
    h5write(outpath, "value", value)
    info("Finished.")
end

if isempty(ARGS)
    println("julia embeds2h5.jl <filename>")
else
    embeds2h5(ARGS[1])
end


text:
using UnivariateFunctions: create_quadratic_spline
using UnivariateFunctions: evaluate, years_between, PE_Function, Sum_Of_Functions
using UnivariateFunctions: change_base_of_PE_Function, derivative, indefinite_integral
using UnivariateFunctions: Piecewise_Function, right_integral, left_integral, evaluate_integral
using UnivariateFunctions: UnivariateFunction, Undefined_Function, years_from_global_base
using Dates

tol = 10*eps()
const global_base_date = Date(2000,1,1)
StartDate = Date(2018, 7, 21)

x = Array{Date}(undef, 1000)
for i in 1:1000
    x[i] = StartDate +Dates.Day(2* (i-1))
end

function ff(x::Date)
    days_between = years_from_global_base(x)
    return log(days_between) + sqrt(days_between)
end
y = ff.(x)

spline = create_quadratic_spline(x,y)
# Test if interpolating
all(abs.(evaluate.(Ref(spline), x) .- y) .< tol)

# Testing third derivatives
third_derivatives = evaluate.(Ref(derivative(derivative(derivative(spline)))), x)
maximum(third_derivatives) < tol

# Testing Integrals
function analytic_integral(lhs,rhs)
    lhs_in_days = years_from_global_base(lhs)
    rhs_in_days = years_from_global_base(rhs)
    return rhs_in_days*log(rhs_in_days) - rhs_in_days + (2/3)*rhs_in_days^(3/2) - (lhs_in_days*log(lhs_in_days) - lhs_in_days)- (2/3)*lhs_in_days^(3/2)
end

lhs = StartDate
rhs = StartDate + Dates.Month(16)
numerical_integral  = evaluate_integral(spline, lhs,rhs)
numerical_integral2 = evaluate(right_integral(spline, lhs), rhs)
numerical_integral3 = evaluate(left_integral(spline, rhs), lhs)

analytical = analytic_integral(lhs,rhs)
abs(  analytical - numerical_integral  ) < 0.0001
abs(  analytical - numerical_integral2  ) < 0.0001
abs(  analytical - numerical_integral3  ) < 0.0001


text:
include("reflection.jl")
include("brig.jl")


text:
using DifferentialEquations
using LinearAlgebra, DiffEqSensitivity, Optim
using Flux: flatten, params
using DiffEqFlux, Flux
using Plots
using Flux:train!
using NNlib

function dxdt_train(du, u, p, t)
    x, y, z = u
    du[1] = 10.0*(y - x)
    du[2] = x*(28.0-z) - y
    du[3] = x*y - (8/3)*z
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
dt = 0.1
prob = ODEProblem(dxdt_train,u0,tspan)
ts = collect(0.0:0.1:tspan[2])
prob_train = ODEProblem{true}(dxdt_train, u0, tspan, p=nothing)
data_train = Array(solve(prob_train, Tsit5(), saveat=ts))
A = [PolynomialBasis(2), PolynomialBasis(2)]
nn0 = TensorLayer(A, 1)
nn1 = TensorLayer(A, 1)
nn2 = TensorLayer(A, 1)
nn2([0.1,0.1])
f = x -> min(100one(x), x)

function dxdt_pred(du, u, p, t)
    x, y, z = u
    du[1] = 10*(y - x)
    du[2] = x*(p[1]-z) - y
    du[3] = x*y -  (8/3)*z
end

function dxdt_pred1(du, u, p, t)
    x, y, z = u
    du[1] = nn0([x,y],p[1:4])[1]
    du[2] = nn1([x,z],p[5:8])[1] - p[13]*y
    du[3] = nn2([x,y],p[9:12])[1] - p[14]*z
end

α = 27*ones(1)

prob_pred = ODEProblem{true}(dxdt_pred, u0, tspan, p=nothing)
function predict_adjoint(θ)
    x = Array(solve(prob_pred, Tsit5(), p=θ, saveat=ts))
end

function loss_adjoint(θ)
    x = predict_adjoint(θ)
    loss = log.(x) - log.(data_train)
    return loss
end

function cb(θ, l)
    @show l
    return false
end
res1 = DiffEqFlux.sciml_train(loss_adjoint, α, ADAM(0.02), cb=cb, maxiters=500)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u, BFGS(), cb=cb)#, maxiters=10000)
opt = res2.u
data_pred = predict_adjoint(opt)
plot(ts, data_train[1,:], label="X (ODE)")
plot!(ts, data_train[2,:], label="Y (ODE)")
plot!(ts, data_train[3,:], label="Z (ODE)")
plot!(ts, data_pred[1,:], label="X (NN)")
plot!(ts, data_pred[2,:],label="V (NN)")
plot!(ts, data_pred[3,:], label="Z (NN)")

using DiffEqFlux, DifferentialEquations, Plots

ann_node = FastChain(FastDense(2, 32, tanh), FastDense(32,2, tanh))

using Pkg, Dates
 Pkg.gc(collect_delay = Day(0))


text:
@testset "stats" begin
  res = stats("BTC-USD")
  @test isa(res, DataFrame)
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

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

equations = CompressibleEulerEquations2D(1.4)

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

Compressible Euler setup for a Mach 3 wind tunnel flow with a forward facing step.
Results in strong shock interactions as well as Kelvin-Helmholtz instabilities at later times.
See Section IV b on the paper below for details.

- Paul Woodward and Phillip Colella (1984)
  The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks.
  [DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6)
"""
@inline function initial_condition_mach3_flow(x, t, equations::CompressibleEulerEquations2D)
  # set the freestream flow parameters
  rho_freestream = 1.4
  v1 = 3.0
  v2 = 0.0
  p_freestream = 1.0

  prim = SVector(rho_freestream, v1, v2, p_freestream)
  return prim2cons(prim, equations)
end

initial_condition = initial_condition_mach3_flow

boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_mach3_flow)

# Outflow boundary condition.
# FIXME: For now only works for supersonic outflow where all values come from the internal state.
# The bones are here for the subsonic outflow as well. One simply needs to pass the reference pressure
# to set it from the outside, the rest comes from the internal solution.
# Once fixed this function should probably move to `compressible_euler_2d.jl`
# See the reference below for a discussion on inflow/outflow boundary conditions.
#
# - Jan-Reneé Carlson (2011)
#   Inflow/Outflow Boundary Conditions with Application to FUN3D.
#   [NASA TM 20110022658](https://ntrs.nasa.gov/citations/20110022658)
@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
                                            surface_flux_function, equations::CompressibleEulerEquations2D)
  # # This would be for the general case where we need to check the magnitude of the local Mach number
  # norm_ = norm(normal_direction)
  # # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
  # normal = normal_direction / norm_

  # # Rotate the internal solution state
  # u_local = Trixi.rotate_to_x(u_inner, normal, equations)

  # # Compute the primitive variables
  # rho_local, v_normal, v_tangent, p_local = cons2prim(u_local, equations)

  # # Compute local Mach number
  # a_local = sqrt( equations.gamma * p_local / rho_local )
  # Mach_local = abs( v_normal / a_local )
  # if Mach_local <= 1.0
  #   p_local = # Set to the external reference pressure value (somehow? maybe stored in `equations`)
  # end

  # # Create the `u_surface` solution state where the local pressure is possibly set from an external value
  # prim = SVector(rho_local, v_normal, v_tangent, p_local)
  # u_boundary = prim2cons(prim, equations)
  # u_surface = Trixi.rotate_from_x(u_boundary, normal, equations)

  # Compute the flux using the appropriate mixture of internal / external solution states
  # flux = Trixi.flux(u_surface, normal_direction, equations)

  # NOTE: Only for the supersonic outflow is this strategy valid
  # Calculate the boundary flux entirely from the internal solution state
  flux = Trixi.flux(u_inner, normal_direction, equations)

  return flux
end

boundary_conditions = Dict( :Bottom      => boundary_condition_slip_wall,
                            :Step_Front  => boundary_condition_slip_wall,
                            :Step_Top    => boundary_condition_slip_wall,
                            :Top         => boundary_condition_slip_wall,
                            :Right       => boundary_condition_outflow,
                            :Left        => boundary_condition_inflow   )

volume_flux = flux_ranocha
surface_flux = flux_lax_friedrichs

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

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

mesh = P4estMesh{2}(mesh_file)

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

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

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

summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval=analysis_interval)

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

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

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

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

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

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

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


text:
## Read data from the mangal.io database

# read the id and other metadata of all ecological networks archived on mangal.io
# id: network ID
# S: number of species in the network
# L: number of interactions in the network
# P: number of interactions of predation in the network
# H: number of interactions of herbivory in the network

number_of_networks = count(MangalNetwork)
count_per_page = 100
number_of_pages = convert(Int, ceil(number_of_networks/count_per_page))

mangal_networks = DataFrame(fill(Int64, 5),
                 [:id, :S, :L, :P, :H],
                 number_of_networks)

global cursor = 1
@progress "Paging networks" for page in 1:number_of_pages
    global cursor
    networks_in_page = Mangal.networks("count" => count_per_page, "page" => page-1)
    @progress "Counting items" for current_network in networks_in_page
        S = count(MangalNode, current_network)
        L = count(MangalInteraction, current_network)
        P = count(MangalInteraction, current_network, "type" => "predation")
        H = count(MangalInteraction, current_network, "type" => "herbivory")
        mangal_networks[cursor,:] .= (current_network.id, S, L, P, H)
        cursor = cursor + 1
    end
end

CSV.write(joinpath("data", "raw", "mangal_networks.csv"), mangal_networks, append=true)

max_stars_repo_name: FrancisBanville/Trophic-METE

max_issues_repo_name: FrancisBanville/Trophic-METE

max_forks_repo_name: FrancisBanville/Trophic-METE