text:
# Check if curl is installed
try
	print("Checking for curl...")
	cmd = `curl --version`
	(VERSION >= v"0.5")? readstring(cmd) : readall(`curl --version`)
	println("OK!")
catch
	error("Curl not found. Please install curl to use BCBData.")
end

"""
A Julia package to read Brazilian Central Bank (BCB) time series data.
"""
module BCBData

	using Base.Dates

	export ProxyConfig,
			readData

	"""
		Holds proxy configuration.
	"""
	type ProxyConfig
		host::String
		port::Int64
		user::String
		password::String

		function ProxyConfig()
			new("", 80, "", "")
		end #function
		function ProxyConfig(	host::String, port::Int64, user::String, password::String)
			new(host, port, user, password)
		end #function
	end # type

	"""
	    readData(code::Int64, startDate::Date, endDate::Date, proxy::ProxyConfig=ProxyConfig())

	Read data from Brazilian Central Bank (BCB) database.

	# Parameters
	* `code::Int64`: time series ID in BCB database. See https://www3.bcb.gov.br/sgspub/ for a list of series and IDs.
	* `startDate::Date`: start date.
	* `endDate::Date`: end date.
	* `proxy::ProxyConfig`: (optional) instance to hold proxy definitions. If not supplied, no proxy config is set.

	# Returns
	Returns two arrays, one with strings representind dates, the other with corresponding numeric values.
	"""
	function readData(code::Int64, startDate::Date, endDate::Date, proxy=ProxyConfig())

		str_start = Dates.format(startDate, "dd/mm/yyyy")
		str_end = Dates.format(endDate, "dd/mm/yyyy")

		if(proxy.host == "")
			str_proxy = ``
		else
			host =
			str_proxy = `--proxy $(proxy.host):$(proxy.port) --proxy-user $(proxy.user):$(proxy.password)`
		end #if

		cmd = `curl -s -S --header "SOAPAction:\"http://publico.ws.casosdeuso.sgs.pec.bcb.gov.br#getValoresSeriesXML\"" --header "Accept: text/xml" --header "Accept: multipart/*" --header "Content-Type: text/xml; charset=utf-8" --header "SOAPAction: \"http://publico.ws.casosdeuso.sgs.pec.bcb.gov.br#getValoresSeriesXML\"" -d "<?xml version=\"1.0\"?><SOAP-ENV:Envelope xmlns:SOAP-ENC=\"http://schemas.xmlsoap.org/soap/encoding/\" xmlns:SOAP-ENV=\"http://schemas.xmlsoap.org/soap/envelope/\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xmlns:xsd=\"http://www.w3.org/2001/XMLSchema\" SOAP-ENV:encodingStyle=\"http://schemas.xmlsoap.org/soap/encoding/\">  <SOAP-ENV:Body>    <getValoresSeriesXML xmlns=\"http://publico.ws.casosdeuso.sgs.pec.bcb.gov.br\">      <in0 xsi:type=\"SOAP-ENC:Array\" SOAP-ENC:arrayType=\"NA[1]\">        <item>$code</item>      </in0>      <in1 xsi:type=\"xsd:string\">$str_start</in1>      <in2 xsi:type=\"xsd:string\">$str_end</in2>    </getValoresSeriesXML>  </SOAP-ENV:Body></SOAP-ENV:Envelope>" $str_proxy https://www3.bcb.gov.br/wssgs/services/FachadaWSSGS`

		response = (VERSION >= v"0.5")? readstring(cmd) : readall(cmd)

		str_dates = String[]
		str_values = Float64[]

		pattern = r"DATA>(.*?)</DATA>\n\t\t\t<VALOR>(.*?)&"

		for m in eachmatch(pattern, response)
			push!(str_dates, m.captures[1])
			push!(str_values, parse(Float64, m.captures[2]))
		end #for

		return(str_dates, str_values)

	end #function

end #module


text:
module TestScikitLearnWrapper

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
nfcp = NumericFeatureClassification()

using FactCheck

importall Orchestra.Types
importall Orchestra.Transformers.ScikitLearnWrapper
using PyCall
@pyimport sklearn.neighbors as NN
@pyimport random as RAN

function skl_fit_and_transform!(learner::Learner, problem::MLProblem, seed=1)
  RAN.seed(seed)
  srand(seed)
  return fit_and_transform!(learner, problem, seed)
end

function backend_fit_and_transform!(sk_learner, seed=1)
  RAN.seed(seed)
  srand(seed)
  sk_learner[:fit](nfcp.train_instances, nfcp.train_labels)
  return collect(sk_learner[:predict](nfcp.test_instances))
end

function behavior_check(learner::Learner, sk_learner)
  # Predict with Orchestra learner
  orchestra_predictions = skl_fit_and_transform!(learner, nfcp)

  # Predict with original backend learner
  original_predictions = backend_fit_and_transform!(sk_learner)

  # Verify same predictions
  @fact orchestra_predictions => original_predictions
end

facts("scikit-learn learners") do
  context("SKLLearner gives same results as its backend") do
    learner_names = collect(keys(ScikitLearnWrapper.learner_dict))
    for learner_name in learner_names
      sk_learner = ScikitLearnWrapper.learner_dict[learner_name]()
      impl_options = Dict()

      fragile_learners = [
        "AdaBoostClassifier",
        "BaggingClassifier",
        "ExtraTreesClassifier",
        "GradientBoostingClassifier",
        "RandomForestClassifier",
        "SGDClassifier",
        "RadiusNeighborsClassifier",
        "RidgeClassifierCV"
      ]

      learner = SKLLearner({
        :learner => learner_name,
        :impl_options => impl_options
      })

      if !in(learner_name, fragile_learners)
        behavior_check(learner, sk_learner)
      else
        skl_fit_and_transform!(learner, nfcp)
        @fact 1 => 1
      end
    end
  end
end

end # module

ext: jl


text:
using Discretizers
using Distances
using LinearAlgebra

function distances!(
    r::AbstractArray{T,3} where T,
    t::Trajectory{TA,NonPeriodic} where {TA<:AbstractAtomNames},
    ur1::AbstractUnitRange,
    ur2::AbstractUnitRange
    )
    metric = Euclidean()
    for i in 1:length(t)
        pairwise!(view(r,:,:,i), metric, view(t,ur1,i), view(t,ur2,i),dims=2)
    end
    return r
end

function distances!(
    r::AbstractArray{T,3} where T,
    t::Trajectory{TA,TC} where {TA<:AbstractAtomNames, TC<:AbstractOrthorombicCell},
    ur1::AbstractUnitRange,
    ur2::AbstractUnitRange
    )
    metric = (PeriodicEuclidean ∘ celldiag)(t)
    for i in 1:length(t)
        pairwise!(view(r,:,:,i), metric, view(t,ur1,i), view(t,ur2,i),dims=2)
    end
    return r
end

function distances!(
    r::AbstractArray{T,3} where T,
    t::Trajectory{TA,VariableCell} where {TA<:AbstractAtomNames},
    ur1::AbstractUnitRange,
    ur2::AbstractUnitRange
    )
    get_metric(c::NonPeriodic) = Euclidean()
    get_metric(c::AbstractOrthorombicCell) = (PeriodicEuclidean ∘ celldiag)(c)
    for i in 1:length(t)
        pairwise!(view(r,:,:,i), get_metric(c.cell[i]), view(t,ur1,i), view(t,ur2,i),dims=2)
    end
    return r
end

function distances(t::AbstractTrajectory, i1::Integer, i2::Integer)
    out = Vector{Float64}(undef, length(t))
    distances!(out, t, i1, i2)
    return out
end

function distances!(
    r::AbstractVector,
    t::Trajectory{TA,NonPeriodic} where {TA<:AbstractAtomNames},
    i1::Integer,
    i2::Integer
    )
    colwise!(r, Euclidean(), view(t,i1,:), view(t,i2,:))
    return r
end

function distances!(
    r::AbstractVector,
    t::Trajectory{TA,TC} where {TA<:AbstractAtomNames, TC<:AbstractOrthorombicCell},
    i1::Integer,
    i2::Integer
    )
    metric = (PeriodicEuclidean ∘ celldiag)(t)
    colwise!(r, metric, view(t,i1,:), view(t,i2,:))
    return r
end

cellvolume(t::AbstractTrajectory) = volume(t.cell)

"""
    angle(t::AbstractTrajectory, i, j, k) -> Vector{Float64}

Computes angle for atoms i, j, k in degrees
"""
function Base.angle(t::AbstractTrajectory, i, j, k)
    @warn "angle does not understand periodicity"
    r1 = view(t,i,:) .- view(t,j,:)
    r2 = view(t,k,:) .- view(t,j,:)
    acos.(1 .- colwise(CosineDist(),r1,r2)) .* 180 ./ π
end

"""
    dihedral(t::AbstractTrajectory, i,j,k,m) -> Vector{Float64}
"""
function dihedral(t::AbstractTrajectory, i,j,k,m)
    out = zeros(length(t))
    @warn "dihedral does not understand periodicity"
    for n in 1:length(t)
        b1 = view(t,j,n) .- view(t,i,n)
        b2 = view(t,k,n) .- view(t,j,n)
        b3 = view(t,m,n) .- view(t,k,n)
        out[n]=atand((b1×b2)×(b2×b3)⋅b2 / norm(b2), (b1×b2)⋅(b2×b3))
    end
    return out
end

"""
    angletoframe(t::AbstractTrajectory, i, j, k, m; frame=:)

Computes angle betweem two vectors. First from atoms i to j second from k to m.
Only given frame is calculated for first given vector.
If frame=: (default) then calculate for whole trajectory.
"""
function angletoframe(t::AbstractTrajectory, i, j, k, m; frame=:)
    r1 = view(t,j,frame) .- view(t,i,frame)
    r2 = view(t,m,:) .- view(t,k,:)
    acosd.(1 .- colwise(CosineDist(),r1,r2))
end

"""
    sphericalview(t::AbstractTrajectory, i, j) -> Dict

Transfroms vector from atom i to j to spherical coordinates
"""
function sphericalview(t::AbstractTrajectory, i, j)
    function rdis(r)
         rr = sqrt.(sum(r.^2,dims=1)) #TODO fix for periodic
         reshape(rr, (length(rr),1))
    end
    r = view(t,j,:) .- view(t,i,:)
    rr = rdis(r)
    ϕ = atan.(view(r,2,:), view(r,1,:))
    θ = acos.( view(r,3,:) ./ rr  )
    return Dict("r"=>rr, "θ"=>θ.*180 ./ π, "ϕ"=>ϕ.*180 ./ π)
end

"""
    compute_rdf(t, ur1, ur2; mindis=undef, maxdis=9.0, nbins=100) -> Dict

Calculates radial distribution function for trajectory or file `t`
"""
function compute_rdf(
    t::Trajectory{TA,TC} where {TA<:AbstractAtomNames, TC<:AbstractOrthorombicCell},
    ur1::AbstractUnitRange,
    ur2::AbstractUnitRange;
    mindis=undef,
    maxdis=9.0,
    nbins=100
    )

    dis = distances(t, ur1, ur2)
    di = DiscretizeUniformWidth(nbins)

    # Collect data that are in defined range [mindis, maxdis]
    data = Dict{Int, Vector{Float64}}()
    for (i,j) in enumerate(ur1)
        tmp = @view dis[i,:,:]
        if mindis != undef
            q1 = tmp .< maxdis
            q2 = tmp .> mindis
            push!(data, j=> tmp[ q1 .& q2 ]  )
        else
            push!(data, j=> tmp[ tmp .< maxdis ]  )
        end
    end

    #Discretize data to form histograms
    ρ = length(ur2)*length(t) / cellvolume(t)  # Average density for ur2  in trajectry
    edges = Dict()
    counts = Dict()
    radius = Dict()
    for (key,val) in data
        if mindis != undef
            e = Array(LinRange(mindis, maxdis, nbins+1))
        else
            e = binedges(di, val)
        end
        r = 0.5 .* (e[1:end-1] .+ e[2:end])
        dr = diff(e)
        disc = LinearDiscretizer(e)
        count = get_discretization_counts(disc, val) ./ (4*π*ρ .* r.^2 .* dr)
        push!(edges, key=> e)
        push!(counts, key=> count)
        push!(radius, key=>r)
    end
    return Dict("r"=>radius, "rdf"=>counts, "edges"=>edges)
end

function compute_rdf(fname::AbstractString, ur1::AbstractUnitRange,
                     ur2::AbstractUnitRange; mindis=undef, maxdis=9.0, nbins=100)
    t = read_trajectory(fname)
    return compute_rdf(t, ur1, ur2, mindis=mindis, maxdis=maxdis, nbins=nbins)
end

"""
    compute_rdf(ur1, ur2, fnames...; mindis=0.0, maxdis=9.0, nbins=100) -> Dict()

Calculates radial distribution function from given input and adds up results
"""
function compute_rdf(ur1::AbstractUnitRange, ur2::AbstractUnitRange,
                        fnames...; mindis=0.0, maxdis=9.0, nbins=100)
    dtmp = pmap(
        x -> compute_rdf(x, ur1, ur2; mindis=mindis, maxdis=maxdis, nbins=nbins),
        fnames
    )

    # Sum up results
    data = []
    for x in dtmp
        if length(keys(x)) == 0
            continue
        end
        push!(data,x)
    end
    rdf = Dict()
    if length(data) >0
        for (k,v) in data[1]["rdf"]
            push!(rdf, k => deepcopy(v))
        end
    end
    if length(data) > 1
        for x in data[2:end]
            for (k,v) in x["rdf"]
                rdf[k] .+= v
            end
        end
    end
    for k in keys(rdf)
        rdf[k] ./= length(data)
    end
    rk = collect(keys(data[1]["r"]))[1]
    return Dict("r"=>data[1]["r"][rk] , "rdf" => rdf)
end

function get_close_atoms(traj::AbstractTrajectory, i; maxdis=5.0)
    t = subtrajectory(traj,1)
    d = distances(t,1:natoms(t),i)[:,:,1]
    m = minimum(d; dims=2)
    a = m .< maxdis
    return getatoms(traj, a[:,1])
end

function get_close_atoms(traj::AbstractTrajectory, i::Integer; maxdis=5.0)
    return get_close_atoms(traj, i:i; maxdis=maxdis)
end


text:
# This example was adapted from: Coding Challenge #115: Snake Game Redux
# by Daniel Shiffman
#
# https://thecodingtrain.com/CodingChallenges/115-snake-game-redux.html

module SnakeGame

using Lynx
using Luxor
using Utils

include("config.jl")
include("food.jl")
include("snake.jl")
include("main.jl")

run!(update, setup, await=true)

end # SnakeGame

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
#const C = -0.62772 -0.42193im
const C = 0.285 + 0.0im
const xmin = -2.
const xmax = 2.
const ymin = -2.
const ymax = 2.
const samples = 2000
@printf("%s", "Dump configuration
 -> C       = $(C)
 -> xrange  = [ $(xmin), $(xmax) ]
 -> yrange  = [ $(ymin), $(ymax) ]
 -> samples = $(samples)\n")

# This is LJR's modified version of Julia Set Generator

function gen_im!(file::AbstractString, xmin::Real, xmax::Real, ymin::Real, ymax::Real, samples::Int, C::Complex)
    y = ((ymin-ymax)/samples)
    x = ((xmax-xmin)/samples)
    f = open(file, "w")
    imgrange = ymax:y:ymin-y
    relrange = xmin:x:xmax-x
    write(f, "P2\n# Julia Set image\n$(samples) $(samples)\n255\n")
    @simd for im in imgrange
        for re in relrange
            z = complex(re, im)
            n = 255
            while abs(z) < 10 && n >= 5
                z = z*z + C
                n -= 5
            end
            write(f, "$n ")
        end
        write(f, "\n")
    end
    close(f)
end

gen_im!("julia.pgm", xmin, xmax, ymin, ymax, 200, C)
t = @elapsed gen_im!("julia.pgm", xmin, xmax, ymin, ymax, samples, C)
println("\n$(t)s\n")

# Her version 2 code

#C = -0.59+0.000im
#xmin = ymin = -1.5
#ymax = xmax = 1.5
#samples = 2500
#@printf("%s", "Dump configuration
# -> C       = $(C)
# -> xrange  = [ $(xmin), $(xmax) ]
# -> yrange  = [ $(ymin), $(ymax) ]
# -> samples = $(samples)\n")
#
#function gen_im!(file::String, Fun::Function, xmin::Real, xmax::Real, ymin::Real, ymax::Real, samples::Int, C::Complex)
#    y = ((ymin-ymax)/samples)
#    x = ((xmax-xmin)/samples)
#    realrange = xmin:x:xmax-x
#    imagrange = ymax:y:ymin-y
#    f = open(file, "w")
#    write(f, "P2\n# Julia Set image\n$(samples) $(samples)\n255\n")
#    A = Array(UInt8, samples)
#    for im in imagrange
#        counter = 0
#        @simd for re in realrange
#            z = complex(re, im)
#            n = 255
#            while abs(z) < 10 && n >= 5
#                z = Fun(z, C)
#                n -= 5
#            end
#            counter += 1
#            A[counter] = n
#        end
#        writedlm(f, A, ' ')
#        write(f, '\n')
#    end
#    close(f)
#end
#
#Fun(z::Complex, C::Complex) = muladd(z, z, C)
#
#gen_im!("julia.pgm", Fun, xmin, xmax, ymin, ymax, 200, C)
#t = @elapsed gen_im!("julia.pgm", Fun, xmin, xmax, ymin, ymax, samples, C)
#println("\n$(t)s\n")

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using TargetedLearning: LReg, CTMLEs, Parameters
import TargetedLearning.CTMLEs: CTMLE, Qmodel, Qchunk, FluctuationInfo

facts("Test CTMLEs") do
    srand(1234)
    n,p = 100, 5
    w = [ones(n) rand(n, p)]
    a = round(rand(n))
    y = round(rand(n))

    qfit = lreg([w a], y)
    logitQnA1 = linpred(qfit, [w ones(n)])
    logitQnA0 = linpred(qfit, [w zeros(n)])

    context("subsetting regimens") do
        @fact CTMLEs.subset(StaticRegimen(1.0), [1,2,3]) --> StaticRegimen(1.0)
        @fact CTMLEs.subset(DynamicRegimen([1.0, 0.0, 1.0]), [1, 3]) -->
            DynamicRegimen([1.0, 1.0])
        subset_ate = CTMLEs.subset(ATE(StaticRegimen(1.0), DynamicRegimen([1.0, 0.0, 1.0]), ),  [1, 3])
        @fact subset_ate.d1 --> StaticRegimen(1)
        @fact subset_ate.d0 --> DynamicRegimen([1.0, 1.0])
    end

    context("making chunk subsets") do
        qc = CTMLEs.make_chunk_subset(logitQnA1, logitQnA0, w, a, y, ATE(), 1:3, [0.0, 1.0], false)
        @fact qc --> is_a(Qchunk)
        @fact CTMLEs.nobs(qc) --> 3
        @fact qc.risk --> CTMLEs.risk(Qmodel(logitQnA1[1:3], logitQnA0[1:3]), a[1:3], y[1:3], ATE(), false)
        qc_pen = CTMLEs.make_chunk_subset(logitQnA1, logitQnA0, w, a, y, ATE(), 1:3, [0.0, 1.0], true)
        @fact qc_pen.risk --> Inf
    end

    context("FluctuationInfo") do
        g1 = LReg.LR(zeros(1), [1], false)
        g2 = LReg.LR(zeros(2), [1, 3], false)
        g1b = LReg.LR(zeros(2), [1, 5], false)

        f1 = FluctuationInfo(1, [g1], [true])
        @fact f1.covar_order --> Vector{Int}[[1]]
        @fact f1.added_each_fluc --> Vector{Int}[[1]]

        f2 = FluctuationInfo(2, [g1, g2], [true, false])
        @fact f2.covar_order -->  Vector{Int}[[1], [3]]
        @fact f2.added_each_fluc --> Vector{Int}[[1, 3]]

        f3 = FluctuationInfo(2, [g1b, g2], [true, true])
        @fact f3.covar_order -->  Vector{Int}[[1, 5], [3]]
        @fact f3.added_each_fluc -->  Vector{Int}[[1, 5], [3]]
    end

    context("pcor") do
        #test pcor
        srand(235)
        xx = rand(100, 3)
        yy = rand(100, 2)
        zz = rand(100)
        oz = [ones(100) zz]

        @fact CTMLEs.pcor(xx,yy,zz) --> roughly(cor(xx .- oz*(oz\xx), yy .- oz*(oz\yy)))
    end

    context("It runs...") do
        est=ctmle(logitQnA1, logitQnA0, w, a, y)
        @fact est --> is_a(CTMLE)
        est2 = ctmle(logitQnA1, logitQnA0, w, a, y, searchstrategy=CTMLEs.PreOrdered(CTMLEs.LogisticOrdering()))
        @fact est2 --> is_a(CTMLE)
        est3 = ctmle(logitQnA1, logitQnA0, w, a, y, searchstrategy=CTMLEs.PreOrdered(CTMLEs.PartialCorrOrdering()))
        @fact est3 --> is_a(CTMLE)
        est4 = ctmle(logitQnA1, logitQnA0, w, a, y, searchstrategy=CTMLEs.PreOrdered(CTMLEs.ManualOrdering([2,3,5,4,6])))
        @fact est4 --> is_a(CTMLE)
        @fact_throws ErrorException ctmle(logitQnA1, logitQnA0, w, a, y, searchstrategy=CTMLEs.PreOrdered(CTMLEs.ManualOrdering([2,3,4,6])))
        @fact_throws ErrorException ctmle(logitQnA1, logitQnA0, w, a, y, searchstrategy=CTMLEs.PreOrdered(CTMLEs.ManualOrdering([1,2,3,5,4,6])))
    end
    @pending "Actually test that it works" --> nothing
end


text:
using TabularDisplay
using Test

io = IOBuffer()

@testset "everything" begin

    displaytable(io, [string("randomstr", i) for i in 1:10];
        displaywidth=50)
    @test String(take!(io)) ==
        "randomstr1   randomstr5   randomstr9   \n" *
        "randomstr2   randomstr6   randomstr10  \n" *
        "randomstr3   randomstr7                \n" *
        "randomstr4   randomstr8                "

    displaytable(io, [string("randomstr", i) for i in 1:10];
        displaywidth=50, orientation=:row)
    @test String(take!(io)) ==
        "randomstr1   randomstr2   randomstr3   \n" *
        "randomstr4   randomstr5   randomstr6   \n" *
        "randomstr7   randomstr8   randomstr9   \n" *
        "randomstr10                            "

    displaytable(io, [string("randomstr", i) for i in 1:10];
        displaywidth=50, orientation=:row, index=true)
    @test String(take!(io)) ==
    "1:randomstr1    2:randomstr2    3:randomstr3    \n" *
    "4:randomstr4    5:randomstr5    6:randomstr6    \n" *
    "7:randomstr7    8:randomstr8    9:randomstr9    \n" *
    "10:randomstr10                                  "

end


text:
"""
    module Symata.SymataIO

This module provides "1d" text representations of Symata expressions. It is
for use at the REPL, and for other purposes, such as writing to files.
"""
module SymataIO

import Symata: Mxpr, SJSym, SSJSym, ListT, TimesT, PowerT,
       getsym, symname, mhead, margs,  getoptype, mtojsym, newargs, mapmargs,
       mxpr, mxprcf, Infinity, getkerneloptions, unicode_output, Qsym,
       CurrentContext, wrapout, using_unicode_output, comparison_translation,
       symnumerator, symdenominator, Formatting

import Symata: print_with_function_to_string

const infix_with_space = Dict( :&& => true , :|| => true, :| => true)

# A space, or maybe not.
opspc(sym) = haskey(infix_with_space,sym) ? " " : getkerneloptions(:compact_output) ? "" : " "
const binaryops_space = collect(keys(comparison_translation) )

## Some of these, at least "^:=" can't be parsed as quoted symbols, but can be
## read as command line input.... So, we use some strings.
const binaryops_space_strings = [ "=>",  "->" , ":>" , "^:=" , "⇒", "∈" ]

function binaryopspc(ss)
    ss in binaryops_space && return " "
    s = string(ss)
    s in binaryops_space_strings && return " "
    return ""
end

## TODO: pass formatting choices to the functions rather than hardcoding them.

# Julia-like syntax
const FUNCL = '('
const FUNCR = ')'
const LISTL = '['
const LISTR = ']'

Istring() =  using_unicode_output() ? "𝕚" : "I"

# We could probably replace all instances of mtojsym below with this.
# It depends on which positions unicode symbols may occur.
function outsym(s)
    haskey(unicode_output, s) && using_unicode_output() && return unicode_output[s]
    os = mtojsym(s)
    haskey(unicode_output, os) && using_unicode_output() && return unicode_output[os]
    os
end

mutable struct FullFormData
    lfunc::String
    rfunc::String
end

const juliafullformdata = FullFormData("(",")")
const mmafullformdata = FullFormData("[","]")

fullform(io::IO, mx::Union{Mxpr, Rational, Complex}) = fullform(io, mx, juliafullformdata)

function fullform(io::IO, mx::Mxpr, data::FullFormData)
    # Really fucking weird. I have to do wrapout here, but not in the method of wrapout(::Mxpr)
    fullform(io, wrapout(mhead(mx)), data)
    print(io, data.lfunc)
    if length(mx) > 0 fullform(io, mx[1], data) end
    for i in 2:length(mx)
        print(io,",")
        fullform(io, mx[i], data)
    end
    print(io, data.rfunc)
end

fullform(io::IO, x, data::FullFormData) = Base.show(io, x)
fullform(x, data::FullFormData) = fullform(stdout, x, data)

fullform(x) = fullform(stdout, x, juliafullformdata)

mmafullform(io::IO, x) = fullform(io, x, mmafullformdata)
mmafullform(x) = mmafullform(stdout, x)

function Base.show(io::IO, mx::Mxpr{:FullForm})
    fullform(io, mx[1], juliafullformdata)
end

symata_to_mma_fullform_string(x) = print_with_function_to_string(mmafullform, wrapout(x))

# This puts unecessary parens on prefix functions.
#needsparen(x::Mxpr) = length(x) > 1
needsparen(x::Mxpr) = (length(x) > 1  && getoptype(mhead(x)) == :infix)
needsparen(x::Rational{T}) where {T<:Integer} = true

needsparen(x::Real) =  x < 0
needsparen(x::Complex{T}) where {T<:Integer} = (real(x) != 0)
needsparen(x::Complex{T}) where {T<:Real} = true

needsparen(x) = false

# Mma displays gensyms with all the linenoise.
# So, we try disabling de_gensym
de_gensym(x) = x
# function de_gensym{T<:AbstractString}(str::T)
#     if str[1] == '#' && str[2] == '#'  # De-gensym local variables for display
#         return split(str,['#'],keepempty=false)[1]
#     else
#         return str
#     end
# end

#### Wrap output

abstract type AbstractWO end

needsparen(y::AbstractWO) = needsparen(y.x)

struct WOSymbol <: AbstractWO
    x::Symbol
end

struct WOBool <: AbstractWO
    x::Bool
end

struct WORational{T}  <: AbstractWO
    x::T
end

abstract type WOComplex <: AbstractWO end

struct WOComplexInteger{T}  <: WOComplex
    x::T
end
wrapout(x::Complex{T}) where {T<:Integer}=  WOComplexInteger(x)
Base.show(io::IO, ws::WOComplexInteger) = show_complexinteger(io, ws.x)

struct WOComplexReal{T}  <: WOComplex
    x::T
end
wrapout(x::Complex{T}) where {T<:Real} = WOComplexReal(x)
Base.show(io::IO,ws::WOComplexReal) = show_complexreal(io, ws.x)

struct WOComplexRational{T}  <: WOComplex
    x::T
end
wrapout(x::Complex{Rational{T}}) where {T<:Integer} = WOComplexRational(x)
Base.show(io::IO,ws::WOComplexRational) = show_complexrational(io, ws.x)

unwrapout(x::AbstractWO) = x.x
unwrapout(x) = x

needsparen(x::WORational) = true
needsparen(x::WOComplexRational) = true

function fullform(io::IO, x::Union{WORational, Rational}, data::FullFormData)
    print(io, "Rational")
    print(io, data.lfunc)
    print(io, numerator(unwrapout(x)), ", ", denominator(unwrapout(x)))
    print(io, data.rfunc)
end

function fullform(io::IO, x::WOComplex, data::FullFormData)
    print(io, "Complex")
    print(io, data.lfunc)
    fullform(io, real(x.x), data)
    print(io, ", ")
    fullform(io, imag(x.x), data)
    print(io, data.rfunc)
end

struct WOAbstractFloat{T}  <: AbstractWO
    x::T
end
wrapout(x::AbstractFloat) = WOAbstractFloat(x)
Base.show(io::IO, ws::WOAbstractFloat) = show_float(io,ws.x)
## Not performant. We do a dictionary lookup every time we print a float.
function show_float(io,x::Float64)
    _show_float(io,x,getkerneloptions(:float_format))
end

function show_float(io,x::BigFloat)
    _show_float(io,x,getkerneloptions(:bigfloat_format))
end

function _show_float(io,x,fmt)
    if fmt == ""
        Base.show(io,x)
    else
        s = Formatting.sprintf1(fmt,x)   # always print a `.` for floating point numbers
        s = occursin(r"\.",s) ? s : s * "."
        print(io,s)
    end
end

function show_float(io,x)
    Base.show(io,x)
end

# This breaks printing
#Base.string{T<:AbstractWO}(y::T) = string(y.x)

#wrapout(x) = x  # defined in wrapout.jl

function wrapout(mx::Mxpr)
    # nargs = newargs()
    # for x in margs(mx)
    #     push!(nargs,wrapout(x))
    # end
    mxprcf(mhead(mx), mapmargs(wrapout, margs(mx)))  # In v0.6, this chooses most efficient Array type. We want Any.
#    mxprcf(mhead(mx), nargs)
end

function wrapout(s::Symbol)
    WOSymbol(s)
end

function wrapout(s::Bool)
    WOBool(s)
end

function wrapout(x::Rational)
    WORational(x)
end

function Base.show(io::IO, ws::WOSymbol)
    show_symbol(io, ws.x)
end

#### Symbol and SSJSym

# NB: This comment is only relevant if we change the Symata symbol implementation.
# This may break. It will only work if the value of s
# is the symbol name in the symbol table that is associated
# with s; ie it is a 'free' symbol. SSJSym does not carry symbol name information.
# This is only stored as the key to the symbol table.

Base.show(io::IO, s::SSJSym) = show_symbol(io, symname(s))

# SJSym is just a symbol, so this translates everything
# We need this, because SSJSym symbols are sometimes stored as Symbol's (ie SJSym)
#Base.show(io::IO, s::SJSym) = show_symbol(io,s)

function show_symbol(io::IO, s::SJSym)
    if haskey(unicode_output, s) && using_unicode_output()
        Base.show_unquoted(io, unicode_output[s])
    else
        ss = string(symname(s))
        ss = de_gensym(ss) # remove gensym characters
        Base.show_unquoted(io, Symbol(ss))
    end
end

Base.show(io::IO, x::WORational) = show_rational(io, x.x)

function show_rational(io::IO, x::Rational)
    Base.show(io, symnumerator(x))
    print(io, "/")
    Base.show(io, symdenominator(x))
end

# Yes, it is REPL.REPLCompletions.latex_symbols
function Base.show(io::IO, x::Mxpr{:Subscript})
    if using_unicode_output()
        try
            s1 = string(x[1].x) # unwrap symbol
            for i in 2:length(x)
                y = typeof(x[i]) <: AbstractWO ? string(x[i].x) : string(x[i])
                s1 *=  REPL.REPLCompletions.latex_symbols["\\_" * y]
            end
            Base.show(io, wrapout(Symbol(s1)))
        catch
            show_prefix_function(io,x)
        end
    else
        show_prefix_function(io,x)
    end
end

# We display real part if it is 0.0
function show_complexreal(io::IO, z::Complex{T}) where T<:Real
#    Base.show(io,real(z))
    show_float(io,real(z))
    print(io," + ")
    show_float(io,imag(z))
#    Base.show(io,imag(z))
    print(io,Istring())
end

# Do not display real part if it is 0
function show_complexinteger(io::IO, z::Complex{T}) where T<:Integer
    if real(z) != 0
        Base.show(io,real(z))
        print(io," + ")
    end
    if imag(z) == 1
        print(io,Istring())
    else
        iz = imag(z)
        if iz == -1
            print(io,"-")
        else
            Base.show(io,iz)
        end
        print(io,Istring())
    end
end

function show_complexrational(io::IO, z::Complex{Rational{T}}) where T<:Integer
    if real(z) != 0
        show_rational(io,real(z))
        print(io," + ")
    end
    if imag(z) == 1
        print(io,Istring())
    else
        print(io,Istring(), " * ")
        show_rational(io,imag(z))
    end
end

function  Base.show(io::IO, mx::Mxpr{:DirectedInfinity})
    if length(mx) == 0
        show_symbol(io, :ComplexInfinity)
    elseif mx[1] == 1
        show_symbol(io, :Infinity)
    elseif mx[1] == -1
        Base.show(io, wrapout(mxprcf(:Times, -1, Infinity)))
    else
        print(io, "DirectedInfinity(")
        print(io, mx[1])
        print(io, ")")
    end
end

# This is causing stack overflow in some cases.
# Base.show{T<:BigFloat}(io::IO,x::T) = Base.showcompact(io,x)

Base.show(io::IO, v::WOBool) = show_bool(io,v.x)

function show_bool(io::IO, v::Bool)
    v ? Base.show_unquoted(io,:True) : Base.show_unquoted(io,:False)
end

# For Holdform, arguments are not evaluated, as in Hold.
# But, in addition, Holdform is not printed.
function Base.show(io::IO, s::Mxpr{:HoldForm})
    Base.show(io,s[1])
end

function Base.show(io::IO, s::Mxpr)
    if getoptype(mhead(s)) == :binary
        return show_binary(io,s)
    elseif getoptype(mhead(s)) == :infix
        return show_infix(io,s)
    end
    show_prefix_function(io,s)
end

function Base.show(io::IO, mx::Mxpr{:Comparison})
    args = margs(mx)
    for i in 1:length(args)-1
        Base.show(io, outsym(args[i]))
        print(io," ")
    end
    isempty(args) || Base.show(io,args[end])
end

function show_prefix_function(io::IO, mx::Mxpr)
    if ! isa(mx,ListT)
        mh = mhead(mx)
        wantparen::Bool = false
        if isa(mh,Mxpr)
            wantparen = getoptype(mhead(mh)) == :infix ? true : false
        end
        if wantparen
            print(io,"(")  # refactor this! --> maybeparen(wantparen,obj) or s.t.
        end
        Base.show(io, wrapout(mhead(mx)))
        if wantparen
            print(io,")")
        end
    end
    args = mx.args
    print(io,mhead(mx) == getsym(:List) ? LISTL : FUNCL)
    wantparens = true
    if mhead(mx) == :List wantparens = false end
    for i in 1:length(args)-1
        if needsparen(args[i]) && wantparens
            print(io,"(")
        end
        Base.show(io,args[i])
        if needsparen(args[i]) && wantparens
            print(io,")")
        end
        print(io,",")
    end
    isempty(args) || Base.show(io,args[end])
    print(io,mx.head == getsym(:List) ? LISTR : FUNCR)
end

function show_binary(io::IO, mx::Mxpr)
    if length(mx) != 2
        show_prefix_function(io,mx)
    else
        opstr = outsym(mhead(mx))
        lop = mx[1]
        if needsparen(lop)
            print(io,"(")
            Base.show(io,lop)
            print(io,")")
        else
            Base.show(io,lop)
        end
        print(io, binaryopspc(opstr), opstr , binaryopspc(opstr))
        rop = mx[2]
        if  needsparen(rop)
            print(io,"(")
            Base.show(io,rop)
            print(io,")")
        else
            Base.show(io,rop)
        end
    end
end

function Base.show(io::IO, mx::Mxpr{:Part})
    args = margs(mx)
    Base.show(io,args[1])
    print(io,"[")
    Base.show(io,args[2])
    for i in 3:length(args)
        print(io,",")
        Base.show(io,args[i])
    end
    print(io,"]")
end

# unary minus
function Base.show(io::IO, mx::Mxpr{:Minus})
    arg = mx.args[1]
    if isa(arg,Number) || isa(arg,SJSym)
        print(io,"-")
        Base.show(io,arg)
    else
        print(io,"-(")
        Base.show(io,arg)
        print(io,")")
    end
end

function Base.show(io::IO, mx::Mxpr{:Plus})
    args = margs(mx)
    if length(args) < 1
        error("show: can't show Plus with no args.")
    end
    Base.show(io,args[1])
    for i in 2:length(args)
        if isa(args[i],Mxpr{:Minus})
            print(io, " - ")
            Base.show(io,(args[i]).args[1])
        else
            if isa(args[i], TimesT) && typeof(args[i][1])  <:Union{AbstractFloat,Integer} && args[i][1] < 0
                show_infix(io, args[i], true)
            else
                print(io, " + ")
                Base.show(io, args[i])
            end
        end
    end
end

# This includes Times. No spaces surround the infix symbol

show_infix(io::IO, mx::Mxpr)  = show_infix(io,mx, false)
function show_infix(io::IO, mx::Mxpr, spaceminus::Bool)
    args = margs(mx)
    np = false
    sepsym = mtojsym(mhead(mx))
    spc = opspc(sepsym)
    # uncomment following to print space for multiplication rather than *. But, I want "InputForm" for now,
    # so we can copy output to input.
#    if mhead(mx) == :Times sepsym = " " end # not a sym. Maybe we should make them all strings
    startind = 1
    if isa(mx,TimesT) && ! isempty(args)
        a1 = args[1]
        if a1 == -1
            print(io, spaceminus ? " - " : "-")
            startind = 2
        elseif typeof(a1) <:Union{AbstractFloat,Integer}
            if a1 < 0 && spaceminus
                print(io, " - ", - a1)
            else
                print(io, a1)
                if length(args) > 1 && isa(args[2],Union{Number,PowerT})
                    print(io, spc, sepsym, spc)
                end
            end
            startind = 2
        end
    end
    for i in startind:length(args)-1
        arg = args[i]
        if needsparen(arg)
            np = true
            print(io,"(")
        else
            np = false
        end
        Base.show(io,arg)
        if np
            print(io,")")
        end
        print(io, spc, sepsym, spc)
    end
    if ! isempty(args)
        if needsparen(args[end])
            np = true
            print(io,"(")
        else
            np = false
        end
        Base.show(io, args[end])
        if np
            print(io,")")
        end
    end
end

function Base.show(io::IO, mx::Mxpr{:Blank})
    print(io,"_")
    if length(mx) > 0
        Base.show(io,mx[1])
    end
end

function Base.show(io::IO, mx::Mxpr{:BlankSequence})
    print(io,"__")
    if length(mx) > 0
        Base.show(io,mx[1])
    end
end

function Base.show(io::IO, mx::Mxpr{:BlankNullSequence})
    print(io,"___")
    if length(mx) > 0
        Base.show(io,mx[1])
    end
end

function Base.show(io::IO, mx::Mxpr{:Pattern})
    Base.show(io,mx[1])
    if isa(mx[2],Mxpr{:Blank})
        Base.show(io,mx[2])
    else
        print(io,"::(")
        Base.show(io,mx[2])
        print(io,")")
    end
end

function Base.show(io::IO, mx::Mxpr{:Element})
    if using_unicode_output()
        show_binary(io,mx)
    else
        show_prefix_function(io,mx)
    end
end

function Base.show(io::IO, qs::Qsym)
    if qs.context != CurrentContext.name
        Base.show(io,qs.context)
        print(io,".")
    end
    Base.show(io,qs.name)
end

end # module SymataIO


text:
using Colors

function colorize_field(images)
    num_colormap_entries = 100
    cmap = colormap("RdBu", num_colormap_entries)

    min, max = extrema(images)
    images_scaled_for_indexing = clamp.((images .- min) / (max - min) * num_colormap_entries .+ 1, 1, num_colormap_entries)
    images_indexed = Int.(round.(images_scaled_for_indexing))

    return cmap[images_indexed]
end


text:
struct TDHH3DSingleLayerNear{T}
    speed_of_light::T
end


text:
export nx_g, ny_g, nz_g, x_g, y_g, z_g, tic, toc

import MPI
using CUDA

macro nx_g()    esc(:( global_grid().nxyz_g[1] )) end
macro ny_g()    esc(:( global_grid().nxyz_g[2] )) end
macro nz_g()    esc(:( global_grid().nxyz_g[3] )) end
macro nx()      esc(:( global_grid().nxyz[1] )) end
macro ny()      esc(:( global_grid().nxyz[2] )) end
macro nz()      esc(:( global_grid().nxyz[3] )) end
macro coordx()  esc(:( global_grid().coords[1] )) end
macro coordy()  esc(:( global_grid().coords[2] )) end
macro coordz()  esc(:( global_grid().coords[3] )) end
macro olx()     esc(:( global_grid().overlaps[1] )) end
macro oly()     esc(:( global_grid().overlaps[2] )) end
macro olz()     esc(:( global_grid().overlaps[3] )) end
macro periodx() esc(:( convert(Bool, global_grid().periods[1]) )) end
macro periody() esc(:( convert(Bool, global_grid().periods[2]) )) end
macro periodz() esc(:( convert(Bool, global_grid().periods[3]) )) end

"""
    nx_g()

Return the size of the global grid in dimension x.
"""
nx_g()::GGInt = @nx_g()

"""
    ny_g()

Return the size of the global grid in dimension y.
"""
ny_g()::GGInt = @ny_g()

"""
    nz_g()

Return the size of the global grid in dimension z.
"""
nz_g()::GGInt = @nz_g()

"""
    nx_g(A)

Return the size of array `A` in the global grid in dimension x.
"""
nx_g(A::GGArray)::GGInt = @nx_g() + (size(A,1)-@nx())

"""
    ny_g(A)

Return the size of array `A` in the global grid in dimension y.
"""
ny_g(A::GGArray)::GGInt = @ny_g() + (size(A,2)-@ny())

"""
    nz_g(A)

Return the size of array `A` in the global grid in dimension z.
"""
nz_g(A::GGArray)::GGInt = @nz_g() + (size(A,3)-@nz())

"""
    x_g(ix, dx, A)

Return the global x-coordinate for the element `ix` in the local array `A` (`dx` is the space step between the elements).

# Examples
```jldoctest
julia> using ImplicitGlobalGrid

julia> lx=4; nx=3; ny=3; nz=3;

julia> init_global_grid(nx, ny, nz);
Global grid: 3x3x3 (nprocs: 1, dims: 1x1x1)

julia> dx = lx/(nx_g()-1)
2.0

julia> A  = zeros(nx,ny,nz);

julia> Vx = zeros(nx+1,ny,nz);

julia> [x_g(ix, dx, A) for ix=1:size(A, 1)]
3-element Array{Float64,1}:
 0.0
 2.0
 4.0

julia> [x_g(ix, dx, Vx) for ix=1:size(Vx, 1)]
4-element Array{Float64,1}:
 -1.0
  1.0
  3.0
  5.0
```
"""
function x_g(ix::Integer, dx::GGNumber, A::GGArray)::GGNumber
    x0 = 0.5*(@nx()-size(A,1))*dx;
    x  = (@coordx()*(@nx()-@olx()) + ix-1)*dx + x0;
    if @periodx()
        x = x - dx;                                     # The first cell of the global problem is a ghost cell; so, all must be shifted by dx to the left.
        if (x > (@nx_g()-1)*dx) x = x - @nx_g()*dx; end # It must not be (nx_g()-1)*dx as the distance between the local problems (1*dx) must also be taken into account!
        if (x < 0)              x = x + @nx_g()*dx; end # ...
    end
    return x
end

"""
    y_g(iy, dy, A)

Return the global y-coordinate for the element `iy` in the local array `A` (`dy` is the space step between the elements).

# Examples
```jldoctest
julia> using ImplicitGlobalGrid

julia> ly=4; nx=3; ny=3; nz=3;

julia> init_global_grid(nx, ny, nz);
Global grid: 3x3x3 (nprocs: 1, dims: 1x1x1)

julia> dy = ly/(ny_g()-1)
2.0

julia> A  = zeros(nx,ny,nz);

julia> Vy = zeros(nx,ny+1,nz);

julia> [y_g(iy, dy, A) for iy=1:size(A, 1)]
3-element Array{Float64,1}:
 0.0
 2.0
 4.0

julia> [y_g(iy, dy, Vy) for iy=1:size(Vy, 2)]
4-element Array{Float64,1}:
 -1.0
  1.0
  3.0
  5.0
```
"""
function y_g(iy::Integer, dy::GGNumber, A::GGArray)::GGNumber
    y0 = 0.5*(@ny()-size(A,2))*dy;
    y  = (@coordy()*(@ny()-@oly()) + iy-1)*dy + y0;
    if @periody()
        y = y - dy;
        if (y > (@ny_g()-1)*dy) y = y - @ny_g()*dy; end
        if (y < 0)              y = y + @ny_g()*dy; end
    end
    return y
end

"""
    z_g(iz, dz, A)

Return the global z-coordinate for the element `iz` in the local array `A` (`dz` is the space step between the elements).

# Examples
```jldoctest
julia> using ImplicitGlobalGrid

julia> lz=4; nx=3; ny=3; nz=3;

julia> init_global_grid(nx, ny, nz);
Global grid: 3x3x3 (nprocs: 1, dims: 1x1x1)

julia> dz = lz/(nz_g()-1)
2.0

julia> A  = zeros(nx,ny,nz);

julia> Vz = zeros(nx,ny,nz+1);

julia> [z_g(iz, dz, A) for iz=1:size(A, 1)]
3-element Array{Float64,1}:
 0.0
 2.0
 4.0

julia> [z_g(iz, dz, Vz) for iz=1:size(Vz, 3)]
4-element Array{Float64,1}:
 -1.0
  1.0
  3.0
  5.0
```
"""
function z_g(iz::Integer, dz::GGNumber, A::GGArray)::GGNumber
    z0 = 0.5*(@nz()-size(A,3))*dz;
    z  = (@coordz()*(@nz()-@olz()) + iz-1)*dz + z0;
    if @periodz()
        z = z - dz;
        if (z > (@nz_g()-1)*dz) z = z - @nz_g()*dz; end
        if (z < 0)              z = z + @nz_g()*dz; end
    end
    return z
end

# Timing tools.
@doc """
    tic()

Start chronometer once all processes have reached this point.

!!! warning
    The chronometer may currently add an overhead of multiple 10th of miliseconds at the first usage.

See also: [`toc`](@ref)
"""
tic

@doc """
    toc()

Return the elapsed time from chronometer (since the last call to `tic`) when all processes have reached this point.

!!! warning
    The chronometer may currently add an overhead of multiple 10th of miliseconds at the first usage.

See also: [`tic`](@ref)
"""
toc

let
    global tic, toc
    t0 = nothing

    tic() = ( check_initialized(); MPI.Barrier(comm()); t0 = time() )
    toc() = ( check_initialized(); MPI.Barrier(comm()); time() - t0 )
end


text:
using StaticArrays: SVector, StaticArray

Base.:+(x::StaticArray) = x

struct DualArray{T,N,V<:AbstractArray,D<:AbstractArray,E,M} <: AbstractArray{E,M}
    data::V
    partials::D
    # ndims(data) + 1 == ndims(partials) when chunking
    # ndims(data) == ndims(partials) when computing Jacobian-vector multiply
    function DualArray{T,N}(v, p) where {T,N}
        # only use SVector for now
        V, D = typeof(v), typeof(p)
        X = SVector{N,eltype(p)}
        E = Dual{T,eltype(v),X}
        M = ndims(v)
        return new{T,N,V,D,E,M}(v, p)
    end
end

isjvp(v, p) = ndims(v) === ndims(p)
isjvp(d::DualArray) = isjvp(data(d), allpartials(d))
function DualArray{T}(v, p) where {T}
    N = isjvp(v, p) ? 1 : size(p, 1)
    return DualArray{T,N}(v, p)
end
DualArray(a::AbstractArray, b::AbstractArray) = DualArray{typeof(dualtag())}(a, b)
data(d::DualArray) = d.data
allpartials(d::DualArray) = d.partials

npartials(d::DualArray) = npartials(typeof(d))
npartials(::Type{<:DualArray{T,N}}) where {T,N} = N

###
### Printing
###

function Base.print_array(io::IO, da::DualArray)
    _dispsize((w, h),) = (w, div(max(h-npartials(da)-1,0), (npartials(da) + 1)))
    sz = :displaysize => _dispsize(get(io, :displaysize, displaysize(io)))

    Base.printstyled(io, "Primals:\n", bold=false, color=2)
    prev_params = io isa IOContext ? io.dict : ()
    ioc = IOContext(io, prev_params..., sz)
    Base.print_array(ioc, data(da))
    Base.println(io)
    Base.printstyled(io,"Partials:\n", bold=false, color=3)
    Base.print_array(ioc, allpartials(da))
    return nothing
end

###
### Array interface
###

Base.size(d::DualArray) = size(data(d))
Base.IndexStyle(d::DualArray) = Base.IndexStyle(data(d))

# Having `Union{}` in `similar` usually means some people did unholy things
# like `promote_op`, and they should stop.
Base.similar(d::DualArray, ::Type{<:Union{}}, dims::Dims) = similar(d, eltype(d), dims)
Base.similar(d::DualArray, ::Type{<:Union{}}) = similar(d, eltype(d))
function Base.similar(d::DualArray{T,N}, ::Type{<:Dual{T,V,P}}, dims::Dims) where {T,N,V,P}
    new_data = similar(data(d), V, dims)
    new_partials = similar(allpartials(d), eltype(P), (N, reverse(dims)...))
    return DualArray{T,N}(new_data, new_partials)
end
function Base.similar(d::DualArray{T,N}, ::Type{<:Dual{T,V,P}}) where {T,N,V,P}
    new_data = similar(data(d), V)
    new_partials = similar(allpartials(d), eltype(P))
    return DualArray{T,N}(new_data, new_partials)
end

Base.eachindex(d::DualArray) = eachindex(data(d))
Base.copy(d::DualArray{T,N}) where {T,N} = DualArray{T,N}(copy(data(d)), copy(allpartials(d)))

# neccsarry to make array code go fast
Base.mightalias(::AbstractArray, ::DualArray) = false
Base.mightalias(::DualArray, ::AbstractArray) = false
Base.mightalias(x::DualArray, y::DualArray) = Base.mightalias(data(x), data(y))

partial_type(::Type{Dual{T,V,P}}) where {T,V,P} = P
Base.@propagate_inbounds function Base.getindex(d::DualArray{T}, i::Int...) where {T}
    ps = allpartials(d)
    P = partial_type(eltype(d))
    partials_ = isjvp(d) ? ps[i...] : ntuple(j->ps[j, i...], Val(npartials(d)))
    return Dual{T}(data(d)[i...], P(partials_))
end

Base.@propagate_inbounds function Base.setindex!(d::DualArray{T}, dual::Dual{T}, i::Int...) where {T}
    data(d)[i...] = value(dual)
    if isjvp(d)
        allpartials(d)[i...] = first(partials(dual))
    else
        allpartials(d)[:, i...] .= partials(dual)
    end
    return dual
end

Base.@propagate_inbounds function Base.setindex!(d::DualArray, x::Number, i::Int...)
    data(d)[i...] = x
    return x
end


text:
module UtilsModule

import Printf: @printf
using Distributed
import ..CoreModule: SRThreaded, SRSerial, SRDistributed

function debug(verbosity, string...)
    if verbosity > 0
        println(string...)
    end
end

function debug_inline(verbosity, string...)
    if verbosity > 0
        print(string...)
    end
end

pseudo_time = 0

function get_birth_order(; deterministic=false)::Int
    """deterministic gives a birth time with perfect resolution, but is not thread safe."""
    if deterministic
        global pseudo_time
        pseudo_time += 1
        return pseudo_time
    else
        resolution = 1e7
        return round(Int, resolution * time())
    end
end

function check_numeric(n)
    return tryparse(Float64, n) !== nothing
end

function is_anonymous_function(op)
    op_string = string(nameof(op))
    return length(op_string) > 1 && op_string[1] == '#' && check_numeric(op_string[2:2])
end

function recursive_merge(x::AbstractVector...)
    return cat(x...; dims=1)
end

function recursive_merge(x::AbstractDict...)
    return merge(recursive_merge, x...)
end

function recursive_merge(x...)
    return x[end]
end

isgood(x::T) where {T<:Number} = !(isnan(x) || !isfinite(x))
isgood(x) = true
isbad(x) = !isgood(x)

macro return_on_false(flag, retval)
    :(
        if !$(esc(flag))
            return ($(esc(retval)), false)
        end
    )
end

# Returns two arrays
macro return_on_false2(flag, retval, retval2)
    :(
        if !$(esc(flag))
            return ($(esc(retval)), $(esc(retval2)), false)
        end
    )
end

function next_worker(worker_assignment::Dict{Tuple{Int,Int},Int}, procs::Vector{Int})::Int
    job_counts = Dict(proc => 0 for proc in procs)
    for (key, value) in worker_assignment
        @assert haskey(job_counts, value)
        job_counts[value] += 1
    end
    least_busy_worker = reduce(
        (proc1, proc2) -> (job_counts[proc1] <= job_counts[proc2] ? proc1 : proc2), procs
    )
    return least_busy_worker
end

function next_worker(worker_assignment::Dict{Tuple{Int,Int},Int}, procs::Nothing)::Int
    return 0
end

macro sr_spawner(parallel, p, expr)
    quote
        if $(esc(parallel)) == SRSerial
            $(esc(expr))
        elseif $(esc(parallel)) == SRDistributed
            @spawnat($(esc(p)), $(esc(expr)))
        else
            Threads.@spawn($(esc(expr)))
        end
    end
end

# Fastest way to check for NaN in an array.
# (due to optimizations in sum())
is_bad_array(array) = !isfinite(sum(array))

end


text:
"""
myTSLS(y, D, instrument, control)

A simple two stage least squares implementation.
"""
struct myTSLS <: myEstimator
    β::Array{Float64} # coefficient
    y::Array{Float64} # response
	Z::Array{Float64} # combined first stage variables
    X::Array{Float64} # combined second stage variables
	FS::Array{Float64} # first stage coefficients

	# Define constructor function
	function myTSLS(y::Array{Float64}, D::Array{Float64},
                    instrument::Array{Float64}, control = nothing)
        # Add constant if no control is passed
        if isnothing(control) control = ones(length(y)) end

		# Define data matrices
		Z = hcat(control, instrument) # combined first stage variables
		X = hcat(D, control) # combined second stage variables

		# Calculate matrix products
		ZZ = Z' * Z
		DZ = X' * Z
		Zy = Z' * y
		FS = inv(ZZ) * DZ'

		# Calculate TSLS coefficient
		β = (DZ * FS)' \ (FS' * Zy)

		# Return output
		new(β, y, Z, X, FS)
	end #MYTSLS
end #MYTSLS

# Methods ======================================================================
"""
inference(fit::myLS; heteroskedastic, print_df)

A method to calculate standard errors of a myLS object.

To do: clustered se
"""
function inference(fit::myTSLS; heteroskedastic::Bool=false, cluster=nothing,
    print_df::Bool=true)
    # Obtain data parameters
    N = length(fit.y)
    Kx = size(fit.X, 2)
    Kz = size(fit.Z, 2)

	# Calculate matrix products
	PZ = fit.Z * fit.FS
	PZZPinv = inv(PZ' * PZ)

	# Calculate covariance matrix
    u = fit.y - predict(fit) # residuals
	if !heteroskedastic
		# Covariance under homoskedasticity
		covar = sum(u.^2) .* PZZPinv ./ (N - Kz)
	else
		# Covariance under heteroskedasticity
		PZuuZP = ((PZ .* (u.^2))' * PZ) .* (N / (N - Kz))
		covar = PZZPinv * PZuuZP * PZZPinv
	end

	# Get standard errors, t-statistics and p-values
    se = sqrt.(covar[diagind(covar)])
    t_stat = fit.β ./ se
    p_val = 2 * cdf.(Normal(), -abs.(t_stat))

    # Print estimates
    if print_df
        out_df = DataFrame(hcat(fit.β, se, t_stat, p_val), :auto)
        rename!(out_df, ["coef", "se", "t-stat", "p-val"])
        display(out_df)
    end

    # Organize and return output
    output = (β = fit.β, se = se, t = t_stat, p = p_val)
    return output
end #INFERENCE.MYTSLS


text:
"""
Temporary solution for gaining google cloud access using access token
obtained from `gcloud` command. It assumes that Google Cloud SDK is
installed and authorized with `gcloud auth login --update-adc` command.
"""
module PoorGCloudAuth

export GCSClient

import HTTP
import JSON
import Dates

include("utils.jl")

"""
    AccessToken()

Obtain a short-lived access token for google cloud api with
`gcloud auth print-access-token` command.
"""
struct AccessToken
    access_token::String
    expiry::Dates.DateTime
    function AccessToken(token::AbstractString)
        info = http_access_token_info(token)
        expiry = Dates.now() + Dates.Second(info["expires_in"])
        new(token, expiry)
    end
    function AccessToken()
        token = read_access_token()
        AccessToken(token)
    end
end

"""
Obtain a new access token if needed.
"""
reauth, reauth!

"""
Get expiration time of the access token.
"""
expiry

"""
Download data from google bucket.
"""
fetch_gcs

"""
    fetch_object_info(gsurl, cred)
    fetch_object_size(gsurl, cred)

Obtain all meta data or object size for a google bucket object.
"""
fetch_object_info, fetch_object_size

function reauth(token::AccessToken)
    # for now, we won't use refresh_token, instead just rerunning "gcloud auth print-access-token"
    AccessToken()
end

function expiry(token::AccessToken)
    token.expiry
end

function fetch_gcs(gsurl::AbstractString, access_token::AccessToken; range=nothing)
    fetch_gcs(gsurl, access_token.access_token; range=range)
end

function fetch_object_info(gsurl::AbstractString, access_token::AccessToken)
    fetch_object_info(gsurl, access_token.access_token)
end

"""
Client that holds google cloud credentials. It manages to create short-lived access token from refresh token.
"""
mutable struct GCSClient
    token::AccessToken
    expiry_limit::Dates.Second
    lock::ReentrantLock
    function GCSClient(token::AccessToken)
        new(token, Dates.Second(Dates.Minute(30)), ReentrantLock())
    end
end

function GCSClient()
    GCSClient(AccessToken())
end

function reauth!(client::GCSClient)
    lock(client.lock)
    try
        time_now = Dates.now()
        # obtain new access token if time passed limit
        if time_now > expiry(client) - client.expiry_limit
            client.token = reauth(client.token)
        end
    finally
        unlock(client.lock)
    end
    client
end

function expiry(client::GCSClient)
    expiry(client.token)
end

"""
Obtain access token from the client.
"""
function accesstoken(client::GCSClient)
    reauth!(client) # reauth if needed
    lock(client.lock) do
        client.token
    end
end

function fetch_gcs(gsurl::AbstractString, client::GCSClient; range=nothing)
    fetch_gcs(gsurl, accesstoken(client); range=range)
end

function fetch_object_info(gsurl::AbstractString, client::GCSClient)
    fetch_object_info(gsurl, accesstoken(client))
end

function fetch_object_size(gsurl::AbstractString, cred)
    parse(Int, fetch_object_info(gsurl, cred)["size"])
end

end


text:
"""
write to a trajectory file
"""
function flush_snapshot(io::IOStream, x::Vector{Float64}, temperature::Float64)
    @printf(io, "%f", temperature)
    for i in length(x)
        @printf(io, " %f", x[i])
    end
    @printf(io, "\n")
end

"""
propagate_mcmc

Markov Chain Monte Carlo test code
"""
function propagate_mcmc(pot_fun::Function, proposal_fun::Function, x_current::Vector{Float64}, temperature::Float64; nstep=1, io=nothing, outstep=1)::Vector{Float64}
    x = copy(x_current)
    x_proposal = similar(x)
    beta = 1 / temperature
    for istep = 1:nstep
        x_proposal .= proposal_fun(x)
        delta_pot = pot_fun(x_proposal) - pot_fun(x)
        if exp( -beta * delta_pot) > rand()
            x .= x_proposal
        end
        if isa(io, IOStream) && (mod(istep, outstep) == 0)
            flush_snapshot(io, x, temperature)
        end
    end
    return x
end

"""
propagate_md

Molecular dynamics test code
"""
function propagate_md(grad_fun::Function, x_current::Vector{Float64}, temperature::Float64; delta_t=0.01, nstep=1, io=nothing, outstep=1)
    x = copy(x_current)
    beta = 1 / temperature
    for istep = 1:nstep
        x_current .= x
        x .= x_current .- (delta_t .* grad_fun(x_current)) .+ (sqrt(2.0 .* delta_t ./ beta) .* randn(size(x_current)))
        if isa(io, IOStream) && (mod(istep, outstep) == 0)
            flush_snapshot(io, x, temperature)
        end
    end
    return x
end


text:
export CompositeShape

struct CompositeShape{K} <: AbstractShape{K}
    c::SFloat{K}  # center of ball
    svec::Vector{AbstractShape{K}}
    CompositeShape{K}(c,svec) where {K} = new(c,svec)
end

CompositeShape(c::AbsVecReal, s::AbstractShape...) = CompositeShape(c, [s...])
CompositeShape(c::AbsVecReal, svec::AbsVec{<:AbstractShape}) = CompositeShape{length(c)}(c,svec)

level(x::SReal{K}, s::CompositeShape{K}) where {K} = mapreduce(sᵢ->level(x,sᵢ), max, s.svec; init=-Inf)
center(s::CompositeShape) = s.c

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
using ConstrainedDynamics
using StaticArrays

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 1.0
x, z = .1, .1
b1 = Box(x, l1, l1, l1 * 1, color = RGBA(1., 1., 0.))
b2 = Box(x, l1, z, l1 * 1, color = RGBA(1., 0., 0.))

vert11 = [0;l1 / 2;l1 / 2]
vert12 = -vert11

vert21 = [0;l1 / 2;0.]
vert22 = -vert21

# Links
origin = Origin{Float64}()
link1 = Body(b1)
link2 = Body(b2)

# Constraints
# joint0to1 = EqualityConstraint(Floating(origin,link1))
# joint0to1 = EqualityConstraint(Fixed(origin,link1))
joint0to1 = EqualityConstraint(Revolute(origin, link1, ex))
joint1to2 = EqualityConstraint(Cylindrical(link1, link2, -ey; p1=vert12, p2=vert21))

links = [link1; link2]
constraints = [joint0to1;joint1to2]
shapes = [b1;b2]

function joint_force_control!(mechanism, k)
    F = SA[0.1; 0.0]
    setForce!(mechanism, joint1to2, F * 0)
    return
end

mech = Mechanism(origin, links, constraints, shapes = shapes, g = 0.)
setPosition!(link1,q = UnitQuaternion(RotX(pi / 2)))
setPosition!(link1,link2,p1 = vert12,p2 = vert21)
# setVelocity!(link1,ω=[2.;0;0])
# setVelocity!(link1,link2,p1=vert12,p2=vert21)

# setForce!([[0.1];nothing],joint1to2,mech)


text:
struct GeometricThinDisc{T} <: AbstractAccretionDisc{T}
    inner_radius::T
    outer_radius::T
    inclination::T
end

function in_nearby_region(d::GeometricThinDisc{T}, line_element) where {T}
    p = line_element[2]
    d.inner_radius < p[1] < d.outer_radius
end

function has_intersect(d::GeometricThinDisc{T}, line_element) where {T}
    u1, u2 = line_element
    sinα = sin(d.inclination)
    cosα = cos(d.inclination)

    s1 = u1[1] * (cosα * sin(u1[2]) * cos(u1[3]) + sinα * cos(u1[2]))
    s2 = u2[1] * (cosα * sin(u2[2]) * cos(u2[3]) + sinα * cos(u2[2]))

    s1 * s2 ≤ 0
end

export GeometricThinDisc


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "tippecanoe"
version = v"1.36.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mapbox/tippecanoe.git",
                  "97ace997ff2d709fe2e26e2dda909cf4afd458ac"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/tippecanoe
make -j${nproc} INCLUDES="-I${includedir} -I." LIBS="-L${libdir}"
mkdir -p $prefix
PREFIX=$prefix make install
"""

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

platforms = filter(
    # Windows outright unsupported;
    # FreeBSD probably needs `#include <sys/statfs.h>` wrapped in an #ifndef upstream
    p -> !(Sys.iswindows(p) || Sys.isfreebsd(p)),
    supported_platforms(),
)

# The products that we will ensure are always built
products = [
    ExecutableProduct("tippecanoe", :tippecanoe),
    ExecutableProduct("tippecanoe-enumerate", :tippecanoe_enumerate),
    ExecutableProduct("tippecanoe-decode", :tippecanoe_decode),
    ExecutableProduct("tippecanoe-json-tool", :tippecanoe_json_tool),
    ExecutableProduct("tile-join", :tile_join),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("SQLite_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    abstract type AbstractParticles end

Abstract supertype for storing particle data efficiently in arrays.

Particles generally have IDs, positions, and velocities as properties.
The properties are expected to be stored in a `Dict{Symbol,Any}` as vectors of length ``N``,
matrices of size ``m×N``, or as scalars (when the property is equal for all particles).

The inbuilt functionality of `AbstractParticles` includes accessing and setting properties
via the following syntax:
```julia
p.id = [1, 2, 3]
p[:id] === p.id
```

If the struct has additional fields, these can also be accessed by `p.field`, but not by
`p[:field]`. The latter syntax can only be used to access the particle properties.

# Property keys
- `:id`: vector of IDs
- `:pos`: ``2×N`` or ``3×N`` matrix with positions
- `:vel`: ``2×N`` or ``3×N`` matrix with velocities

Any arrays may be of the type `AbstractArray`, provided the arrays are 1-indexed.

# Methods
The methods `Base.keys`, `Base.values`, `Base.haskey`, `Base.empty`, `Base.empty!`, `Base.isempty`,
and `Base.copy!` are forwarded to the property `Dict`.

Concrete types of `AbstractParticles` should have the following methods implemented
(also see the implementation of [`Particles`](@ref)):
- `Base.copy`: returns new object containing a copy of the `Dict`
- `Base.empty`: returns an identical object, but with an empty `Dict`
- `Base.:(==)`: should call `Base.isequal` to prevent obtaining `missing` results for properties with missing values (relevant for [`AllParticles`](@ref))
- [`CosmoParticles.particle_name`](@ref): returns the name of the struct to be printed via `Base.show`
- `Base.propertynames`: implement this if there are additional struct fields
- `Base.show(io, mime, p)`
"""
abstract type AbstractParticles end

"""
    CosmoParticles.get_props(p::AbstractParticles)

Return the property `Dict` belonging to the particles.

This returns `p.props` by default if not overridden.

This is not exported and should not be used outside of the defining files for particle types.
"""
get_props(p::AbstractParticles) = p.props

function Base.getproperty(p::AP, sym::Symbol) where {AP<:AbstractParticles}
    if sym in fieldnames(AP)
        return getfield(p, sym)
    else
        return getindex(p, sym)
    end
end

function Base.setproperty!(p::AP, sym::Symbol, val) where {AP<:AbstractParticles}
    if sym in fieldnames(AP)
        setfield!(p, sym, val)
    else
        setindex!(p, val, sym)
    end
end

Base.getindex(p::AbstractParticles, sym::Symbol) = getindex(get_props(p), sym)
Base.setindex!(p::AbstractParticles, val, sym::Symbol) = (setindex!(get_props(p), val, sym); p)
Base.keys(p::AbstractParticles) = keys(get_props(p))
Base.haskey(p::AbstractParticles, key) = key in keys(p)
Base.values(p::AbstractParticles) = values(get_props(p))
Base.propertynames(p::AbstractParticles) = keys(p) |> collect
Base.empty!(p::AbstractParticles) = (empty!(get_props(p)); p)
Base.isempty(p::AbstractParticles) = isempty(get_props(p))
Base.copy!(dst::AbstractParticles, src::AbstractParticles) = (copy!(get_props(dst), get_props(src)); dst)

Base.getindex(p::AbstractParticles, ind::AbstractVector) = applyind(p, ind)
Base.deleteat!(p::AbstractParticles, ind::AbstractVector) = removeind!(p, ind)
deleteat(p::AbstractParticles, ind::AbstractVector) = removeind(p, ind)

# to implement:
# Base.copy(p::AbstractParticles) = AbstractParticles(copy(p.props))
# Base.empty(p::AbstractParticles) = AbstractParticles(empty(p.props))
# Base.:(==)(p1::AbstractParticles, p2::AbstractParticles) = isequal(p1.props, p2.props)

"""
    CosmoParticles.particle_name(p::AbstractParticles)

Returns the name of the particle type, which is used when printing an object of this type via
`Base.show`.

This is not exported.
"""
function particle_name(::AbstractParticles) end

"""
    CosmoParticles.particle_number(p::AbstractParticles)

Returns the number of particles by determining the length of one of the property arrays.

This is not exported.
"""
function particle_number(p::AbstractParticles)
    for val in values(p)
        if val isa AbstractArray
            return size(val, ndims(val))
        end
    end
    return 0
end

"""
    CosmoParticles.show_properties(io::IO, mime, p::AbstractParticles)

Prints the number of particles, the name of the type, and the property names.

Should be used internally when overriding `Base.show` for concrete implementations of
`AbstractParticles`.

This is not exported.
"""
show_properties(io::IO, mime::AbstractString, p::AbstractParticles) = show_properties(io, MIME(mime), p)

function show_properties(io::IO, ::MIME"text/plain", p::AbstractParticles)
    n = particle_number(p)
    name = particle_name(p)
    println(io, "$n $name")
    print(io, " ")
    props = join(keys(p) |> collect |> sort, " ")
    isempty(props) || print(io, props)
end


text:
"""
Solve least squares problem via steepest descent.
"""
struct SD <: AbstractMMAlg end

# Initialize data structures.
function __mm_init__(::SD, problem, ::Nothing)
    @unpack X, coeff = problem
    n, p, _ = probdims(problem)
    nparams = ifelse(problem.kernel isa Nothing, p, n)

    # residuals subroutine requires an object named Z; need to fix
    Z = nothing

    return (;
    projection=StructuredL0Projection(nparams), Z=Z,
    )
end

# Check for data structure allocations; otherwise initialize.
function __mm_init__(::SD, problem, extras)
    if :projection in keys(extras) && :Z in keys(extras) # TODO
        return extras
    else
        __mm_init__(SD(), problem, nothing)
    end
end

# Update data structures due to change in model subsets, k.
__mm_update_sparsity__(::SD, problem, ϵ, ρ, k, extras) = nothing

# Update data structures due to changing ρ.
__mm_update_rho__(::SD, problem, ϵ, ρ, k, extras) = nothing

# Update data structures due to changing λ.
__mm_update_lambda__(::SD, problem, ϵ, λ, extras) = nothing

# Apply one update.
function __mm_iterate__(::SD, problem, ϵ, ρ, k, extras)
    @unpack coeff, proj, grad, res = problem
    @unpack projection = extras
    X = get_design_matrix(problem)
    n, _, _ = probdims(problem)
    T = floattype(problem)
    a², b² = 1/n, ρ

    # Project and then evaluate gradient.
    apply_projection(projection, problem, k)
    __evaluate_residuals__(problem, ϵ, extras, true, true, false)
    __evaluate_gradient__(problem, ρ, extras)

    # Find optimal step size
    B, G, XG = coeff.all, grad.all, res.main.all
    mul!(XG, X, G)
    C1 = dot(G, G)
    C2 = dot(XG, XG)
    t = C1 / (a²*C2 + b²*C1 + eps())

    # Move in the direction of steepest descent.
    axpy!(-t, G, B)

    return nothing
end

# Apply one update in regularized problem.
function __reg_iterate__(::SD, problem, ϵ, λ, extras)
    @unpack coeff, grad, res = problem
    X = get_design_matrix(problem)
    n, _, _ = probdims(problem)
    a², b² = 1/n, λ

    # Evaluate the gradient using residuals.
    __evaluate_residuals__(problem, ϵ, extras, true, false, false)
    __evaluate_reg_gradient__(problem, λ, extras)

    # Find optimal step size
    B, G, XG = coeff.all, grad.all, res.main.all
    mul!(XG, X, G)
    C1 = dot(G, G)
    C2 = dot(XG, XG)
    t = C1 / (a²*C2 + b²*C1 + eps())

    # Move in the direction of steepest descent.
    axpy!(-t, G, B)

    return nothing
end


text:
const CACHE_INCLUSIVITY = let
    if !((Sys.ARCH === :x86_64) || (Sys.ARCH === :i686))
         (false,false,false,false)
    else
        # source: https://github.com/m-j-w/CpuId.jl/blob/401b638cb5a020557bce7daaf130963fb9c915f0/src/CpuInstructions.jl#L38
        # credit Markus J. Weber, copyright: https://github.com/m-j-w/CpuId.jl/blob/master/LICENSE.md
        function get_cache_edx(subleaf)
            Base.llvmcall(
                """
                    ; leaf = %0, subleaf = %1, %2 is some label
                    ; call 'cpuid' with arguments loaded into registers EAX = leaf, ECX = subleaf
                    %2 = tail call { i32, i32, i32, i32 } asm sideeffect "cpuid",
                        "={ax},={bx},={cx},={dx},{ax},{cx},~{dirflag},~{fpsr},~{flags}"
                        (i32 4, i32 %0) #2
                    ; retrieve the result values and return eax and edx contents
                    %3 = extractvalue { i32, i32, i32, i32 } %2, 0
                    %4 = extractvalue { i32, i32, i32, i32 } %2, 3
                    %5  = insertvalue [2 x i32] undef, i32 %3, 0
                    %6  = insertvalue [2 x i32]   %5 , i32 %4, 1
                    ; return the value
                    ret [2 x i32] %6
                """
                # llvmcall requires actual types, rather than the usual (...) tuple
                , Tuple{UInt32,UInt32}, Tuple{UInt32}, subleaf % UInt32
            )
        end
        # eax0, edx1 = get_cache_edx(0x00000000)
        t = (false,false,false,false)
        i = zero(UInt32)
        j = 0
        while (j < 4)
            eax, edx = get_cache_edx(i)
            i += one(UInt32)
            iszero(eax & 0x1f) && break
            iszero(eax & 0x01) && continue
            ci = ((edx & 0x00000002) != 0x00000000) & (eax & 0x1f != 0x00000000)
            t = Base.setindex(t, ci, (j += 1))
        end
        t
    end
end


text:
module TestUtils

using LinearAlgebra
using Random
using Statistics: mean
using Test

using StaticArrays: SVector, SMatrix

using Bifurcations.ArrayUtils: container_array_of, _lq!, canonicalize, nan_
using Bifurcations.Codim2: cast_container, as_reals, _ds_eigvec

@testset "container_array_of" begin
    @test container_array_of([1, 2, 3]) == Vector
    @test container_array_of([1 2; 3 4]) == Matrix
    @test container_array_of(SVector(1.0, 2.0, 3.0)) === SVector{3}
    @test container_array_of(SMatrix{2, 2}(1, 2, 3, 4)) === SMatrix{2, 2}
end

@testset "cast_container" begin
    # No cast if it's already of the correct type:
    for v in [ones(3),
              SVector(1.0, 2.0),
              SVector(1 + 2im, 3 + 4im)]
        @test cast_container(typeof(v), v) === v
    end

    v = [1.0, 2.0]
    @test cast_container(SVector{2, Float64}, v) === SVector{2, Float64}(v)

    v = [1.0 + 2.0im, 3.0 + 4.0im]
    @test cast_container(SVector{2, Float64}, v) === SVector{2, ComplexF64}(v)
end

@testset "_lq!" begin
    rng = MersenneTwister(0)
    for _ in 1:10, array_type in [Array, SMatrix{3, 3}]
        A = array_type(randn(rng, (3, 3)))

        L0, Q0 = lq(Array(A))

        L1, Q1 = _lq!(copy(A))

        @test L0 ≈ L1
        @test Q0 ≈ Q1
    end
end

@testset "canonicalize" begin
    x_list = Any[
        [1.0, 2],
        [1.0 + 2im, 3 + 4im],
        [1.0 + 0im, 0.0 + 1im],
    ]
    for x in [x1 for x0 in x_list for x1 in [x0, SVector(x0...)]]
        y = canonicalize(x)
        @test typeof(y) === typeof(x)
        @test norm(y) ≈ 1
        @test abs(real(y) ⋅ imag(y)) < 1e-7

        # Test idempotence
        z = canonicalize(y)
        @test typeof(z) === typeof(x)
        @test norm(z) ≈ 1
        @test abs(real(z) ⋅ imag(z)) < 1e-7
    end
end

nanmean(args...) = nan_(mean, args...)
ceil_to_int(x) = ceil(Int, x)

@testset "nan_(mean)" begin
    @testset "$v (-> finite)" for v in [
            10:30,
            collect(-3:-1) .+ 0.0,
            ]
        @test nanmean(v) == mean(v)
        @testset "$f" for f in [
                exp,
                sin,
                ceil_to_int,
                ]
            @test nanmean(f, v) == mean(f, v)
        end
    end
    @testset "$v (-> NaN)" for v in [
            Float64[],
            Int[],
            ]
        @test isnan(nanmean(v))
        # Note: `nanmean(f, v)` throws but that's OK; `mean(f, v)` does so too
    end
    for (x, y) in [
            ([1, NaN, 2, 3], 1:3),
            ([1, NaN, 2, 3], collect(1:3) .+ 0.0),
            ]
        @test nanmean(x) == mean(y)
    end

    A = [
        1   2 3   NaN
        NaN 2 NaN 6
        7   8 9 10
    ]
    @test nanmean(A, 1) == [4  4  6  8]
    @test nanmean(A, 2) == [2 4 8.5]'
end

@testset "as_reals" begin
    # No convert if it's already a real vector;
    for v in [ones(3),
              SVector(1.0, 2.0)]
        @test as_reals(v) === v
    end

    vc = [1.0 + 2.0im, 3.0 + 4.0im]
    vr = [1.0, 2.0, 3.0, 4.0]
    @test as_reals(vc) == vr
    @test as_reals(SVector{2, ComplexF64}(vc)) == SVector{4, Float64}(vr)
end

@testset "_ds_eigvec" begin
    @test _ds_eigvec(Float64, [0, 0, 1, 1, 2, 2]) == [1, 1]
    @test _ds_eigvec(ComplexF64, [0, 0, 1, 1, 1, 1, 3, 2, 2]) ==
        [1 + 1im, 1 + 1im]

    @test _ds_eigvec(Float64, SVector(0, 0, 1, 1, 2, 2)) ==
        SVector(1, 1)
    @test _ds_eigvec(ComplexF64, SVector(0, 0, 1, 1, 1, 1, 3, 2, 2)) ==
        SVector(1 + 1im, 1 + 1im)
end

end  # module


text:
"""
Die Funktion divide_jl ist eine julia Implementierung von divide.py.
Diese Funktion arbeitet nicht mit .txt Dateien sondern mit arrays.

### Argumente
- `data::Array`: enthält die Messzeiten (erster Parameter)
- `Δt::Real`: Länge der Intervalle in s (zweiter Parameter)

### Rückgabewert
- `counts::Array`: enthält Anzahl der Signale für jedes Interval

### Beispiel
```julia
counts = divide_jl(data, 10)
```
"""
function divide_jl(data::Array, Δt::Real)
    counts = []
    for x in data
        i = round(Int, x/Δt, RoundDown)
        while length(counts) <= i
            append!(counts, 0)
        end
        counts[i+1] += 1
    end
    pop!(counts)
    return counts
end
"""
Die Funktion binomial_jl ist eine julia Implementierung von binomial.py.
Diese Funktion arbeitet nicht mit .txt Dateien sondern mit arrays.

### Argumente
- `data::Array`: enthält die Messzeiten (erster Parameter)
- `Δt::Real`: Länge der Intervalle in s (zweiter Parameter)
- `n::Int`: Binningparameter (dritter Parameter)

### Rückgabewert
- `x::Array`: enthält die x-Werte der Binomialverteilung
- `y::Array`: enthält die y-Werte der Binomialverteilung

### Beispiel

```julia
x, y = binomial_jl(data, 0.2, 1)
```
"""
function binomial_jl(data::Array, Δt::Real, n::Int)
    counts = divide_jl(data, Δt)
    hist1 = []
    for i in 1:length(counts)
        while length(hist1) < counts[i]+1
            append!(hist1, 0)
        end
        hist1[counts[i]+1] += 1
    end

    hist2 = []
    for i in 1:length(hist1)
        j = round(Int, (i-1)/n, RoundDown)
        while length(hist2) <= j
            append!(hist2, 0)
        end
        hist2[j+1] += hist1[i]
    end

    return [i for i in 0:length(hist2)-1], hist2
end
"""
Die Funktion interval_jl ist eine julia Implementierung von interval.py.
Diese Funktion arbeitet nicht mit .txt Dateien sondern mit arrays.

### Argumente
- `data::Array`: enthält die Messzeiten (erster Parameter)
- `Δt::Real`: Länge der Intervalle in s (zweiter Parameter)
- `n::Int`: step, gibt an ob Zeit bis zum nächsten, übernächsten, ... betrachtet wird (dritter Parameter)

### Rückgabewert
- `x::Array`: enthält die x-Werte der Intervallverteilung
- `y::Array`: enthält die y-Werte der Intervallverteilung

### Beispiel

```julia
x, y = interval_jl(data, 0.001, 1)
```
"""
function interval_jl(data::Array, Δt::Real, n::Int)
    hist = []
    for i in 1:length(data)-n
        t = data[i+n]-data[i]
        k = round(Int, t/Δt, RoundDown)
        while length(hist) <= k
            append!(hist, 0)
        end
        hist[k+1] += 1
    end
    return [i for i in 0:length(hist)-1], hist
end


text:
# example from https://github.com/malmaud/TensorFlow.jl
using Distributions
using Printf

# Generate some synthetic data
x = randn(100, 50)
w = randn(50, 10)
y_prob = exp.(x*w)
y_prob ./= sum(y_prob,dims=2)

function draw(probs)
    y = zeros(size(probs))
    for i in 1:size(probs, 1)
        idx = rand(Categorical(probs[i, :]))
        y[i, idx] = 1
    end
    return y
end

y = draw(y_prob)

# Build the model
reset_default_graph()
sess = Session()

X = placeholder(Float64, shape=[nothing, 50])
Y_obs = placeholder(Float64, shape=[nothing, 10])

variable_scope("test_logisitic_model"; initializer=random_normal_initializer(0, .001), reuse=AUTO_REUSE) do
    global W = get_variable("W", shape=[50, 10], dtype=Float64)
    global B = get_variable("B", shape=[10], dtype=Float64)
end

Y=softmax(X*W + B)

Loss = -sum(log(Y).*Y_obs)
# adam = AdamOptimizer()
# minimize_op = minimize(adam, Loss)

# # Run training
# init(sess)
# for epoch in 1:1000
#     cur_loss, _ = run(sess, [Loss, minimize_op], feed_dict=Dict(X=>x, Y_obs=>y))
#     println(@sprintf("Current loss is %.2f.", cur_loss))
# end

function print_loss(loss_evaled)
    println(loss_evaled)
end

opt = ScipyOptimizerInterface(Loss, method="L-BFGS-B",options=Dict("maxiter"=> 100))
init(sess)
ScipyOptimizerMinimize(sess, opt, feed_dict=Dict(X=>x, Y_obs=>y), loss_callback=print_loss, fetches=[Loss])


text:
include("../base.jl")
include("fitnessFunctions.jl")

setRandomSeed(5)
setIndividualType(BinaryGenotype)
setStopCondition(maxIterations=500)
setEvaluator(FitnessFunction(countBinaryOnes))
setGenerator(randomBinaryGenerator, 50, popSize = 100)
setSelector(tournamentSelector, 20, samplingWithRep=true, nSelected=2.)
setCrossoverOperator(uniformCross, nChildren=2)
setMutationOperator(standardMutation, 0.1)
setReplacementOperator(replaceAllPopulation)

println("HOla")
println("HOla")
a=typeof(precompile(cross_, (CrossoverOperator, Array{Individual}, ExperimentInfo)))

runGenJ(verbose=false)

"""
genPopulation!(GenJ)

evaluate!(GenJ, GenJ._population) # TODO Tarda mas, revisar

setSelector(tournamentSelector, 20, samplingWithRep=true, nSelected=200)

selected = selectParents(GenJ)

@time offspring = cross(GenJ, selected) # TODO Tarda maxSymbols

offspring = mutate(GenJ, offspring)

evaluate!(GenJ, offspring)

replacePopulation!(GenJ, offspring) # TODO Tarda maxSymbols
"""

"""
setIndividualType(BinaryGenotype)
setAlgorithm(SPEA, 20)
setStopCondition(maxIterations=100)
setEvaluator(countBinaryOnes, 1, countBinaryZeros, 1, globalFitnessFunction=FitnessFunction(pareto))
setGenerator(randomBinaryGenerator, 50, popSize = 100)
setSelector(tournamentSelector, 4)
setCrossoverOperator(uniformCross)
setMutationOperator(standardMutation, 0.1)
setReplacementOperator(replaceAllPopulation)

@time runGenJ(verbose=false)

#a = GenJ._evaluator._compareFitness([10, 5], [7, 10], 1, 2)

population = genPopulation(GenJ._generator, GenJ._experimentInfo)
evaluate!(GenJ._evaluator, GenJ._experimentInfo, population, GenJ._stopCondition)
weights = getWeights(GenJ._evaluator)

conjunto1 = Array{Individual}(undef, 10)
conjunto2 = [population[2], population[3], population[1]]

launion=(union(conjunto1,conjunto2))

paretoLevel = pareto(population, weights)

"""

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using Test

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/Classic.jl")
    using .Classic
else
    # invoked during test
    using Classic
end

function test_SimpleDataPoint()
    # Constructor
    data = Float64.([1,2])
    p = SimpleDataPoint(data)
    @test length(data) == p.dimension
    @test data == p.original
    @test data == p.derived

    # (==) checks `derived` for equality
    x = SimpleDataPoint(Float64.([1,2]))
    y = SimpleDataPoint(Float64.([1,2]))
    @test x == y

    x = SimpleDataPoint(Float64.([1,2]))
    y = SimpleDataPoint(Float64.([1,2]))
    x.original[1] = 13
    @test x == y

    x = SimpleDataPoint(Float64.([1,2]))
    y = SimpleDataPoint(Float64.([1,2]))
    x.derived[1] = 13
    @test x != y

end

function test_set_original!()
    data = Float64.([1,2])
    p = SimpleDataPoint(data)
    @test data == p.original
    @test data == p.derived

    data = Float64.([7,13])
    set_original!(p, data)
    @test data == p.original
    @test data == p.derived
end

function test_distance()
    # distance
    p1 = SimpleDataPoint(Float64.([1,2]))
    p2 = SimpleDataPoint(Float64.([2,1]))
    @test sqrt(2) ≈ distance(p1,p2)

    # Error: Points with different dimension
    p3 = SimpleDataPoint(Float64.([2]))
    @test_throws DomainError distance(p1,p3)
end

function test_Cluster()
end

function test_KMeans()
    # Constructor
    pointvec = Vector{SimpleDataPoint}()
    push!(pointvec, SimpleDataPoint([1.0,4.0]))
    @test KMeans(13, pointvec) isa Any

    # Error
    @test_throws DomainError KMeans(0, pointvec)
end

function test_slice_derived()
    # Setup
    pointvec = Vector{SimpleDataPoint}()
    push!(pointvec, SimpleDataPoint([1.0,4.0]))
    push!(pointvec, SimpleDataPoint([2.0,5.0]))
    push!(pointvec, SimpleDataPoint([3.0,6.0]))
    pointvec[1].derived .= [7.0, 10.0]
    pointvec[2].derived .= [8.0, 11.0]
    pointvec[3].derived .= [9.0, 12.0]
    # Test
    data = slice_derived(pointvec, 1)
    @test Float64.([7,8,9]) == data
    data = slice_derived(pointvec, 2)
    @test Float64.([10,11,12]) == data
end

function test_slice_original()
    # Setup
    pointvec = Vector{SimpleDataPoint}()
    push!(pointvec, SimpleDataPoint([1.0,4.0]))
    push!(pointvec, SimpleDataPoint([2.0,5.0]))
    push!(pointvec, SimpleDataPoint([3.0,6.0]))
    pointvec[1].derived .= [7.0, 10.0]
    pointvec[2].derived .= [8.0, 11.0]
    pointvec[3].derived .= [9.0, 12.0]

    # Test
    data = slice_original(pointvec, 1)
    @test Float64.([1,2,3]) == data
    data = slice_original(pointvec, 2)
    @test Float64.([4,5,6]) == data
end

function test_kmeans(io::IO)
    pointvec = Vector{SimpleDataPoint}()
    push!(pointvec, SimpleDataPoint(Float64.([2,1,1])))
    push!(pointvec, SimpleDataPoint(Float64.([2,2,5])))
    push!(pointvec, SimpleDataPoint(Float64.([3,1.5,2.5])))

    kmeans = KMeans(2, pointvec)
    run!(kmeans, 100)
    for (ind,cluster) in enumerate(kmeans.clustervec)
        println(io, "Cluster $(ind): $(cluster.pointvec)")
    end
    # for cluster in kmeans.clustervec
    #     if 1 == length(cluster.pointvec)
    #         point = cluster.pointvec[1]
    #         @test point.original == Float64.([2,2,5])
    #     elseif 2 == length(cluster.pointvec)
    #         for point in cluster.pointvec
    #             @test point.original == Float64.([2,1,1]) || point.original == Float64.([3,1.5,2.5])
    #         end
    #     end
    # end
end

function test_Governor()
    p = Governor(-89.79113, 72, "Alabama")
    @test 2 == p.dimension
    @test [-89.79113, 72] == p.original
    @test [-89.79113, 72] == p.derived
    @test -89.79113       == p.longitude
    @test 72              == p.age
    @test "Alabama"       == p.state
    @test false == (p.original === p.derived)

    data = Float64.([7, 13])
    p = Governor(data)
    @test 2    == p.dimension
    @test data == p.original
    @test data == p.derived
end

function test_kmeans_governor(io::IO)
    pointvec = Vector{Governor}()
    begin # add data
        push!(pointvec, Governor(-86.79113, 72, "Alabama"))
        push!(pointvec, Governor(-152.404419, 66, "Alaska"))
        push!(pointvec, Governor(-111.431221, 53, "Arizona"))
        push!(pointvec, Governor(-92.373123, 66, "Arkansas"))
        push!(pointvec, Governor(-119.681564, 79, "California"))
        push!(pointvec, Governor(-105.311104, 65, "Colorado"))
        push!(pointvec, Governor(-72.755371, 61, "Connecticut"))
        push!(pointvec, Governor(-75.507141, 61, "Delaware"))
        push!(pointvec, Governor(-81.686783, 64, "Florida"))
        push!(pointvec, Governor(-83.643074, 74, "Georgia"))
        push!(pointvec, Governor(-157.498337, 60, "Hawaii"))
        push!(pointvec, Governor(-114.478828, 75, "Idaho"))
        push!(pointvec, Governor(-88.986137, 60, "Illinois"))
        push!(pointvec, Governor(-86.258278, 49, "Indiana"))
        push!(pointvec, Governor(-93.210526, 57, "Iowa"))
        push!(pointvec, Governor(-96.726486, 60, "Kansas"))
        push!(pointvec, Governor(-84.670067, 50, "Kentucky"))
        push!(pointvec, Governor(-91.867805, 50, "Louisiana"))
        push!(pointvec, Governor(-69.381927, 68, "Maine"))
        push!(pointvec, Governor(-76.802101, 61, "Maryland"))
        push!(pointvec, Governor(-71.530106, 60, "Massachusetts"))
        push!(pointvec, Governor(-84.536095, 58, "Michigan"))
        push!(pointvec, Governor(-93.900192, 70, "Minnesota"))
        push!(pointvec, Governor(-89.678696, 62, "Mississippi"))
        push!(pointvec, Governor(-92.288368, 43, "Missouri"))
        push!(pointvec, Governor(-110.454353, 51, "Montana"))
        push!(pointvec, Governor(-98.268082, 52, "Nebraska"))
        push!(pointvec, Governor(-117.055374, 53, "Nevada"))
        push!(pointvec, Governor(-71.563896, 42, "New Hampshire"))
        push!(pointvec, Governor(-74.521011, 54, "New Jersey"))
        push!(pointvec, Governor(-106.248482, 57, "New Mexico"))
        push!(pointvec, Governor(-74.948051, 59, "New York"))
        push!(pointvec, Governor(-79.806419, 60, "North Carolina"))
        push!(pointvec, Governor(-99.784012, 60, "North Dakota"))
        push!(pointvec, Governor(-82.764915, 65, "Ohio"))
        push!(pointvec, Governor(-96.928917, 62, "Oklahoma"))
        push!(pointvec, Governor(-122.070938, 56, "Oregon"))
        push!(pointvec, Governor(-77.209755, 68, "Pennsylvania"))
        push!(pointvec, Governor(-71.51178, 46, "Rhode Island"))
        push!(pointvec, Governor(-80.945007, 70, "South Carolina"))
        push!(pointvec, Governor(-99.438828, 64, "South Dakota"))
        push!(pointvec, Governor(-86.692345, 58, "Tennessee"))
        push!(pointvec, Governor(-97.563461, 59, "Texas"))
        push!(pointvec, Governor(-111.862434, 70, "Utah"))
        push!(pointvec, Governor(-72.710686, 58, "Vermont"))
        push!(pointvec, Governor(-78.169968, 60, "Virginia"))
        push!(pointvec, Governor(-121.490494, 66, "Washington"))
        push!(pointvec, Governor(-80.954453, 66, "West Virginia"))
        push!(pointvec, Governor(-89.616508, 49, "Wisconsin"))
        push!(pointvec, Governor(-107.30249, 55, "Wyoming"))
    end

    #println(io, pointvec)
    kmeans = KMeans(2, pointvec)
    run!(kmeans, 100)
    for (ind,cluster) in enumerate(kmeans.clustervec)
        println(io, "Cluster $(ind):")
        for point in cluster.pointvec
            println(io, "$(point.state)\t\t$(point.age)\t\t$(point.longitude)")
        end
    end
end

function test_Album()
    p = Album("Got to Be There", 1972, 35.45, 10.0)
    @test 2 == p.dimension
    @test [35.45, 10.0] == p.original
    @test [35.45, 10.0] == p.derived
    @test "Got to Be There" == p.name
    @test 1972 == p.year
    @test 35.45 == p.length
    @test 10.0 == p.tracks

    data = [35.45, 10.0]
    p = Album(data)
    @test 2 == p.dimension
    @test data == p.original
    @test data == p.derived
end

function test_kmeans_album(io::IO)
    pointvec = Vector{Album}()
    begin
        push!(pointvec, Album("Got to Be There", 1972, 35.45, 10.0))
        push!(pointvec, Album("Ben", 1972, 31.31, 10.0))
        push!(pointvec, Album("Music & Me", 1973, 32.09, 10.0))
        push!(pointvec, Album("Forever, Michael", 1975, 33.36, 10.0))
        push!(pointvec, Album("Off the Wall", 1979, 42.28, 10.0))
        push!(pointvec, Album("Thriller", 1982, 42.19, 9.0))
        push!(pointvec, Album("Bad", 1987, 48.16, 10.0))
        push!(pointvec, Album("Dangerous", 1991, 77.03, 14.0))
        push!(pointvec, Album("HIStory: Past, Present and Future, Book I", 1995, 148.58, 30.0))
        push!(pointvec, Album("Invincible", 2001, 77.05, 16.0))
    end

    kmeans = KMeans(2, pointvec)
    run!(kmeans, 100)
    for (ind,cluster) in enumerate(kmeans.clustervec)
        println(io, "Cluster $(ind):")
        for point in cluster.pointvec
            println(io, "$(point.name)\t\t$(point.length)\t\t$(point.tracks)")
        end
    end

end

io = stdout
#io = devnull

@testset "DataPoint" begin
    test_SimpleDataPoint()
    test_Cluster()
    test_KMeans()
    test_Governor()
    test_Album()

    test_distance()
    test_set_original!()
    test_slice_derived()
    test_slice_original()
    test_kmeans(io)
    test_kmeans_governor(io)
    test_kmeans_album(io)
end

nothing


text:
is_positive(cur_value::SymEngine.Basic) = true

is_positive(cur_value::Any) = ( cur_value > 0 )


text:
@testset "Two ligand model." begin
    @testset "LS: Test if bothLigands function works." begin
        pp = TAMode.Lsparam(fill(0.2, 9))
        ss = TAMode.getAutocrine(pp)
    end

    @testset "LS: Make sure code upholds mass conservation." begin
        tt = TAMode.Lsparam(fill(0.2, 9))
        firstV = TAMode.getAutocrine(tt)

        tt.kDeg = 0
        tt.expression = 0
        tt.curL = (100, 100)

        secondV = TAMode.runTAMinit([1.0e6], tt, firstV)

        @test dot(firstV, TAMode.totalLS) - dot(secondV, TAMode.totalLS) < 0.0001
    end

    @testset "LS: Ensure that system reaches detailed balance." begin
        rr = TAMode.Lsparam(fill(0.2, 9))
        rr.curL = (0.0, 0.0)

        autoC = TAMode.getAutocrine(rr)

        rr.expression = 0.0
        rr.kDeg = 0.0

        uLong = TAMode.runTAMinit([1.0e7], rr, autoC)

        # Get the Jacobian matrix
        du = zero(uLong)
        J = ForwardDiff.jacobian((y, x) -> TAMode.TAMreact(y, x, rr, 0.0), du, uLong)
        GK = J * diagm(vec(uLong))

        IDXmismatch = findall(GK - transpose(GK) .≠ 0.0)

        @test length(IDXmismatch) == 2
        @test_broken norm(GK - transpose(GK)) < 1.0e-5
    end

    @testset "LS: Make sure that surface total is preserved without trafficking." begin
        tt = TAMode.Lsparam(fill(0.2, 9))
        firstSurf = TAMode.getAutocrine(tt)

        tt.expression = 0.0
        tt.kRec = 0.0
        tt.internalize = 0.0
        tt.pYinternalize = 0.0
        tt.curL = (10.0, 10.0)

        secondSurf = TAMode.runTAMinit([100.0], tt, firstSurf)

        reduce = TAMode.surfaceLS .* TAMode.totalLS
        @test isapprox(dot(firstSurf, reduce), dot(secondSurf, reduce), rtol = 1.0e-5)
    end

    @testset "LS: Ensure that system reaches equilibrium." begin
        tt = TAMode.Lsparam(fill(0.2, 9))

        uLong = TAMode.getAutocrine(tt)
        dnorm = zeros(30)
        TAMode.TAMreact(dnorm, uLong, tt, 0.0)
        @test norm(dnorm) .< 0.001
    end

    @testset "LS: Test that if no Gas6 is present, we don’t see any." begin
        rr = TAMode.Lsparam(fill(0.2, 9))
        rr.curL = (0.0, 0.0)

        data = TAMode.runTAM(tps, rr, (0.0, 10.0))

        @test all(aboutZero.(data * TAMode.GasLS))
    end

    @testset "LS: Test that if no Protein S is present, we don’t see any." begin
        rr = TAMode.Lsparam(fill(0.2, 9))
        rr.curL = (0.0, 0.0)

        data = TAMode.runTAM(tps, rr, (10.0, 0.0))

        @test all(aboutZero.(data * TAMode.PROSLS))
    end

    @testset "LS: Test that if no ligand is present, we don’t see any pY." begin
        rr = TAMode.Lsparam(fill(0.2, 9))
        rr.curL = (0.0, 0.0)

        data = TAMode.runTAM(tps, rr, (0.0, 0.0))

        @test all(aboutZero.(data * TAMode.pYLS))
        @test all(aboutZero.(data * TAMode.PROSLS))
        @test all(aboutZero.(data * TAMode.GasLS))
    end

    @testset "LS: Test for symmetry wrt to both ligands." begin
        rr = TAMode.Lsparam(fill(0.2, 9))

        Gstim = TAMode.runTAM(tps, rr, (10.0, 0.0))
        Pstim = TAMode.runTAM(tps, rr, (0.0, 10.0))
        dataDiff = Gstim .- Pstim

        @test all(aboutZero.(dataDiff * TAMode.pYLS))
        @test all(aboutZero.(dataDiff * TAMode.totalLS))
        @test all(aboutZero.(dataDiff * TAMode.surfaceLS))
        @test all(aboutZero.(Gstim * TAMode.GasLS .- Pstim * TAMode.PROSLS))
    end
end


text:
module Keras

include("./wrapper.jl")
include("./utils.jl")

end # module


text:
# Test ability to support different memory layouts
module ExtraColor

using Color
import Color.Fractional

immutable BGR{T<:Fractional} <: AbstractRGB{T}
    b::T
    g::T
    r::T

    BGR(r::Number, g::Number, b::Number) = new(b, g, r)
end
BGR(r::Integer, g::Integer, b::Integer) = BGR{Float64}(r, g, b)
BGR(r::Fractional, g::Fractional, b::Fractional) = (T = promote_type(typeof(r), typeof(g), typeof(b)); BGR{T}(r, g, b))

immutable ARGB{T<:Fractional} <: AbstractAlphaColorValue{RGB{T}, T}
    alpha::T
    c::RGB{T}

    ARGB(x1::T, x2::T, x3::T, alpha::T) = new(alpha, RGB{T}(x1, x2, x3))
    ARGB(c::RGB{T}, alpha::T) = new(alpha, c)
end
ARGB{T<:Fractional}(c::RGB{T}, alpha::T = one(T)) = ARGB{T}(c, alpha)

end

using Color

const redbgr = ExtraColor.BGR(1, 0, 0)
const redrgb = RGB(1, 0, 0)
@test convert(HSV, redbgr) == convert(HSV, redrgb)

ac = rgba(redrgb)
ac_argb = ExtraColor.ARGB(redrgb)
@test convert(Uint32, convert(ARGB32, ac)) == convert(Uint32, convert(ARGB32, ac_argb))

if ENDIAN_BOM == 0x04030201
    # little-endian
    const redbgr8 = ExtraColor.BGR(0xffuf8, 0x00uf8, 0x00uf8)
    acv = AlphaColorValue(redbgr8, 0xffuf8)
    @test convert(Uint32, argb32(redbgr8)) == reinterpret(Uint32, typeof(acv)[acv])[1]
end


text:
using LoopVectorization, Test

const GAMMA = 1.4
@testset "Reject Unroll" begin
  function r_manual_inline!(r, u)
    @inbounds @fastmath begin
      m = size(u,2)
      n = size(u,3)
      fill!(r, 0.0)
      for j=2:n-1
        for i=1:m-1
          rl = u[1,i,j]
          rul = u[2,i,j]
          rvl = u[3,i,j]
          rEl = u[4,i,j]

          rr = u[1,i+1,j]
          rur = u[2,i+1,j]
          rvr = u[3,i+1,j]
          rEr = u[4,i+1,j]

          rr1  = 1.0/rr
          ur   = rur*rr1
          vr   = rvr*rr1
          Er   = rEr*rr1
          u2r  = ur*ur+vr*vr
          pr   = (GAMMA-1.0)*(rEr-0.5*rr*u2r)
          hr   = Er+pr*rr1
          unr  = ur

          rhr  = rEr+pr

          rl1  = 1.0/rl
          ul   = rul*rl1
          vl   = rvl*rl1
          El   = rEl*rl1
          u2l  = ul*ul+vl*vl
          pl   = (GAMMA-1.0)*(rEl-0.5*rl*u2l)
          hl   = El+pl*rl1
          unl  = ul

          rhl  = rEl+pl

          di    = sqrt(rr*rl1)
          d1    = 1.0/(di+1.0)
          ui    = (di*ur+ul)*d1
          vi    = (di*vr+vl)*d1
          hi    = (di*hr+hl)*d1
          ci2   = (GAMMA-1.0)*(hi-0.5*(ui*ui+vi*vi))
          ci    = sqrt(ci2)
          af    = 0.5*(ui*ui+vi*vi)
          uni   = ui

          dr    = rr-rl
          dru   = rur-rul
          drv   = rvr-rvl
          drE   = rEr-rEl

          rlam1 = abs(uni+ci)
          rlam2 = abs(uni-ci)
          rlam3 = abs(uni)

          s1    = 0.5*(rlam1+rlam2)
          s2    = 0.5*(rlam1-rlam2)
          al1x  = (GAMMA-1.0)*(af*dr-ui*dru-vi*drv+drE)
          al2x  = -uni*dr+dru
          cc1   = ((s1-rlam3)*al1x/ci2)+(s2*al2x/ci)
          cc2   = (s2*al1x/ci)+(s1-rlam3)*al2x

          f1 = 0.5*(rr*unr+rl*unl)               - 0.5*(rlam3*dr+cc1)
          f2 = 0.5*(rur*unr+rul*unl)+0.5*(pr+pl) - 0.5*(rlam3*dru+cc1*ui+cc2)
          f3 = 0.5*(rvr*unr+rvl*unl)             - 0.5*(rlam3*drv+cc1*vi)
          f4 = 0.5*(rhr*unr+rhl*unl)             - 0.5*(rlam3*drE+cc1*hi+cc2*uni)

          r[1,i,j] -= f1
          r[2,i,j] -= f2
          r[3,i,j] -= f3
          r[4,i,j] -= f4

          r[1,i+1,j] += f1
          r[2,i+1,j] += f2
          r[3,i+1,j] += f3
          r[4,i+1,j] += f4
        end
      end
    end
    r
  end

  function r_turbo!(r, u)
    @inbounds begin
      m = size(u,2)
      n = size(u,3)
      fill!(r, 0.0)
      @turbo for j=2:n-1
        for i=1:m-1
          rl = u[1,i,j]
          rul = u[2,i,j]
          rvl = u[3,i,j]
          rEl = u[4,i,j]

          rr = u[1,i+1,j]
          rur = u[2,i+1,j]
          rvr = u[3,i+1,j]
          rEr = u[4,i+1,j]

          rr1  = 1.0/rr
          ur   = rur*rr1
          vr   = rvr*rr1
          Er   = rEr*rr1
          u2r  = ur*ur+vr*vr
          pr   = (GAMMA-1.0)*(rEr-0.5*rr*u2r)
          hr   = Er+pr*rr1
          unr  = ur

          rhr  = rEr+pr

          rl1  = 1.0/rl
          ul   = rul*rl1
          vl   = rvl*rl1
          El   = rEl*rl1
          u2l  = ul*ul+vl*vl
          pl   = (GAMMA-1.0)*(rEl-0.5*rl*u2l)
          hl   = El+pl*rl1
          unl  = ul

          rhl  = rEl+pl

          di    = sqrt(rr*rl1)
          d1    = 1.0/(di+1.0)
          ui    = (di*ur+ul)*d1
          vi    = (di*vr+vl)*d1
          hi    = (di*hr+hl)*d1
          ci2   = (GAMMA-1.0)*(hi-0.5*(ui*ui+vi*vi))
          ci    = sqrt(ci2)
          af    = 0.5*(ui*ui+vi*vi)
          uni   = ui

          dr    = rr-rl
          dru   = rur-rul
          drv   = rvr-rvl
          drE   = rEr-rEl

          rlam1 = abs(uni+ci)
          rlam2 = abs(uni-ci)
          rlam3 = abs(uni)

          s1    = 0.5*(rlam1+rlam2)
          s2    = 0.5*(rlam1-rlam2)
          al1x  = (GAMMA-1.0)*(af*dr-ui*dru-vi*drv+drE)
          al2x  = -uni*dr+dru
          cc1   = ((s1-rlam3)*al1x/ci2)+(s2*al2x/ci)
          cc2   = (s2*al1x/ci)+(s1-rlam3)*al2x

          f1 = 0.5*(rr*unr+rl*unl)               - 0.5*(rlam3*dr+cc1)
          f2 = 0.5*(rur*unr+rul*unl)+0.5*(pr+pl) - 0.5*(rlam3*dru+cc1*ui+cc2)
          f3 = 0.5*(rvr*unr+rvl*unl)             - 0.5*(rlam3*drv+cc1*vi)
          f4 = 0.5*(rhr*unr+rhl*unl)             - 0.5*(rlam3*drE+cc1*hi+cc2*uni)

          r[1,i,j] -= f1
          r[2,i,j] -= f2
          r[3,i,j] -= f3
          r[4,i,j] -= f4

          r[1,i+1,j] += f1
          r[2,i+1,j] += f2
          r[3,i+1,j] += f3
          r[4,i+1,j] += f4
        end
      end
    end
    r
  end

  u = rand(4,128,128)
  u[[1,4],:,:] .= 1.0
  u[[2,3],:,:] .= 0.1
  @test r_manual_inline!(similar(u),u) ≈ r_turbo!(similar(u),u)

end


text:
"""
    RieszRepresenterCotangentVector(M::AbstractManifold, p, X)

Cotangent vector in Riesz representer form on manifold `M` at point `p` with Riesz
representer `X`.
"""
struct RieszRepresenterCotangentVector{TM<:AbstractManifold,TP,TX}
    manifold::TM
    p::TP
    X::TX
end

function allocate(ξ::RieszRepresenterCotangentVector)
    return RieszRepresenterCotangentVector(ξ.manifold, copy(ξ.p), allocate(ξ.X))
end

function (ξ::RieszRepresenterCotangentVector)(Y)
    return inner(ξ.manifold, ξ.p, ξ.X, Y)
end

@decorator_transparent_signature flat!(
    M::AbstractDecoratorManifold,
    ξ::CoTFVector,
    p,
    X::TFVector,
)

@doc raw"""
    flat(M::AbstractManifold, p, X)

Compute the flat isomorphism (one of the musical isomorphisms) of tangent vector `X`
from the vector space of type `M` at point `p` from the underlying [`AbstractManifold`](@ref).

The function can be used for example to transform vectors
from the tangent bundle to vectors from the cotangent bundle
$♭ : T\mathcal M → T^{*}\mathcal M$
"""
flat(M::AbstractManifold, p, X) = RieszRepresenterCotangentVector(M, p, X)
function flat(M::AbstractManifold, p, X::TFVector{<:Any,<:AbstractBasis})
    return CoTFVector(X.data, dual_basis(M, p, X.basis))
end

function flat!(::AbstractManifold, ξ::RieszRepresenterCotangentVector, p, X)
    # TODO: maybe assert that ξ.p is equal to p? Allowing for varying p in ξ leads to
    # issues with power manifold.
    copyto!(ξ.X, X)
    return ξ
end

function get_coordinates(
    M::AbstractManifold,
    p,
    ξ::RieszRepresenterCotangentVector,
    ::DefaultOrthonormalBasis{𝔽,CotangentSpaceType},
) where {𝔽}
    return get_coordinates(M, p, ξ.X, DefaultOrthonormalBasis{𝔽}())
end

function get_coordinates!(
    M::AbstractManifold,
    v,
    p,
    ξ::RieszRepresenterCotangentVector,
    ::DefaultOrthonormalBasis{𝔽,CotangentSpaceType},
) where {𝔽}
    get_coordinates!(M, v, p, ξ.X, DefaultOrthonormalBasis{𝔽}())
    return v
end

function get_vector(
    M::AbstractManifold,
    p,
    v,
    ::DefaultOrthonormalBasis{𝔽,CotangentSpaceType},
) where {𝔽}
    X = get_vector(M, p, v, DefaultOrthonormalBasis{𝔽}())
    return RieszRepresenterCotangentVector(M, p, X)
end

function get_vector!(
    M::AbstractManifold,
    ξr::RieszRepresenterCotangentVector,
    p,
    v,
    ::DefaultOrthonormalBasis{𝔽,CotangentSpaceType},
) where {𝔽}
    get_vector!(M, ξr.X, p, v, DefaultOrthonormalBasis{𝔽}())
    return ξr
end

@doc raw"""
    sharp(M::AbstractManifold, p, ξ)

Compute the sharp isomorphism (one of the musical isomorphisms) of vector `ξ`
from the vector space `M` at point `p` from the underlying [`AbstractManifold`](@ref).

The function can be used for example to transform vectors
from the cotangent bundle to vectors from the tangent bundle
$♯ : T^{*}\mathcal M → T\mathcal M$
"""
sharp(::AbstractManifold, p, ξ)

sharp(::AbstractManifold, p, ξ::RieszRepresenterCotangentVector) = ξ.X
function sharp(M::AbstractManifold, p, X::CoTFVector{<:Any,<:AbstractBasis})
    return TFVector(X.data, dual_basis(M, p, X.basis))
end

@decorator_transparent_signature sharp!(
    M::AbstractDecoratorManifold,
    X::TFVector,
    p,
    ξ::CoTFVector,
)

function sharp!(::AbstractManifold, X, p, ξ::RieszRepresenterCotangentVector)
    copyto!(X, ξ.X)
    return X
end


text:
module NLARX

include("nlarx_node_description.jl")
include("nlarx_rules_prototypes.jl")
include("nlarx_update_rules.jl")
include("util.jl")

end

max_stars_repo_name: wmkouw/NLARX

max_issues_repo_name: wmkouw/NLARX

max_forks_repo_name: wmkouw/NLARX