text:
module Pathogen

  # Dependencies
  using Distributed,
        DataFrames,
        Distributions,
        RecipesBase,
        Logging,
        StatsBase,
        Statistics,
        ProgressMeter,
        LinearAlgebra,
        OnlineStats,
        PhyloModels,
        Random

  # Methods for functions not in Base
  import StatsBase.mean,
         StatsBase.mode,
         ProgressMeter.next!

  # Core
  include("core/IndividualLevelModel.jl")
  include("core/DiseaseStateSequence.jl")
  include("core/DiseaseState.jl")
  include("core/Population.jl")
  include("core/Events/Event.jl")
  include("core/Events/Events.jl")
  include("core/Events/EventRates.jl")
  include("core/Events/EventObservations.jl")
  include("core/Risks/AbstractRisk.jl")
  include("core/Risks/RiskFunctions.jl")
  include("core/Risks/RiskParameters.jl")
  include("core/Transmissions/Transmission.jl")
  include("core/Transmissions/AbstractTransmissionNetwork.jl")
  include("core/Transmissions/TransmissionNetwork.jl")
  include("core/Transmissions/TransmissionRates.jl")
  include("core/Transmissions/TransmissionRateCache.jl")
  include("core/initialize.jl")
  include("core/update!.jl")

  # Simulation
  include("simulation/Simulation.jl")
  include("simulation/generate/Event.jl")
  include("simulation/generate/Tree.jl")
  include("simulation/generate/EventObservations.jl")
  include("simulation/generate/Transmission.jl")
  include("simulation/update!.jl")
  include("simulation/simulate!.jl")

  # Inference
  include("inference/TransmissionNetworkDistribution.jl")
  include("inference/RiskPriors.jl")
  include("inference/EventExtents.jl")
  include("inference/MarkovChain.jl")
  include("inference/MCMC.jl")
  include("inference/_accept.jl")
  include("inference/_bounds.jl")
  include("inference/generate.jl")
  include("inference/initialize.jl")
  include("inference/update!.jl")
  include("inference/loglikelihood.jl")
  include("inference/logprior.jl")
  include("inference/start!.jl")
  include("inference/iterate!.jl")
  include("inference/summary.jl")

  # Visualization
  include("visualization/epidemic_curve.jl")
  include("visualization/epidemic_curve_distribution.jl")
  include("visualization/observations.jl")
  include("visualization/population.jl")
  include("visualization/transmission_network.jl")
  include("visualization/transmission_network_distribution.jl")
  include("visualization/degree_distribution.jl")
  include("visualization/trace.jl")

  # re-export PhyloModels, DataFrames, Distributions
  for pkg in [PhyloModels, DataFrames, Distributions]
    for name in names(pkg)
      @eval export $(name)
    end
  end

  export
    ILM, TransmissionNetworkILM, TNILM, PhylodynamicILM, PhyloILM,
    SEIR, SEI, SIR, SI, DiseaseStateSequence,
    DiseaseState, DiseaseStates,
    State_S, State_E, State_I, State_R,
    RiskFunctions, RiskParameters, RiskPriors,
    Population,
    AbstractTransmissionNetwork,
    TransmissionNetwork,
    TransmissionNetworkDistribution, TNDistribution,
    TransmissionNetworkPrior, TNPrior,
    TransmissionNetworkPosterior, TNPosterior,
    EndogenousTransmission, ExogenousTransmission, NoTransmission,
    Simulation,
    next!, simulate!, generate,
    Events, EventObservations, EventExtents,
    observe,
    MCMC, start!, iterate!,
    summary, mean, mode, innerjoin, leftantijoin
end


text:
@testset "test: newick parsing" begin
global net
@testset "readTopology White Symbol Tests" begin
    global n1, n2, n3
	#Test newlines, spaces, tabs, and carriage returns
	n1 = readTopology("(A,((B,#H1),(C,(D)#H1)));")
	n2 = readTopology("(A\n,\t (\n\r (B , #H1 ),( C ,(D )#H1 ) \n\n\n\n\n) );")
	@test writeTopology(n1) == writeTopology(n2)

	#Test spaces in troublesome locations
	n3 = readTopology("( A, ( ( B , #H 1 ) ,( C , ( D )#H1 ) ) ) ;")
	@test writeTopology(n1) == writeTopology(n3)

	#Test spaces in names
	resultant = readTopology("('Homo sapiens', ((B,#H1),(C,(D)#H1)));")
	expected = readTopology("('Homosapiens',((B,#H1),(C,(D)#H1)));")
	@test writeTopology(resultant) == writeTopology(expected)

	#Test other white symbols in names
	resultant = readTopology("('H\tom\ro\nsa   p \n\n\n\n\n\n\n\ni\t\t\t\t\t\t\t\ten s   ', ((B,#H1),(C,(D)#H1)));")
	expected = readTopology("('Homosapiens',((B,#H1),(C,(D)#H1)));")
	@test writeTopology(resultant) == writeTopology(expected)
end
@testset "isMajor and gamma consistency" begin
	net = readTopology("((((B)#H1)#H2,((D,C,#H2:::0.8),(#H1,A))));");
	@test writeTopology(net, round=true, digits=8) == "(#H2:::0.2,((D,C,((B)#H1)#H2:::0.8),(#H1,A)));"
	net = readTopologyLevel1("(E,((B)#H1:::.5,((D,C),(#H1:::.5,A))));");
	@test writeTopology(net) == "(D:1.0,C:1.0,((#H1:1.0::0.5,A:1.0):1.0,((B:1.0)#H1:1.0::0.5,E:1.0):1.0):1.0);"
    originalstdout = stdout
	redirect_stdout(open("/dev/null", "w")) # not portable to Windows
	@test_logs PhyloNetworks.printEverything(net)
	redirect_stdout(originalstdout)
end
@testset "internal nodes, writemulti" begin
	@test writeTopology(readTopology("(a,b):0.5;")) == "(a,b);"
	@test writeTopology(readTopology("((a,(b)#H1)i1,(#H1,c))r;")) == "((a,(b)#H1)i1,(#H1,c))r;"
	@test writeTopology(readTopology("((a,(b)#H1)i1,(#H1,c))r;"), internallabel=false) == "((a,(b)#H1),(#H1,c));"
	n11 = (@test_logs readTopology("((a,(b)#H1)i1,(#H1,c)i2)root:0.5;"));
	n12 = (@test_logs readTopology("(((a,(b)#H2)i1,(#H2,c)i2)root:0.5);")); # root edge was deleted
	# writeTopology(net) == "((a,(b)#H1)i1,(#H1,c)i2);"
	# readTopology("((((a,(b)#H1)i1,(#H1,c)i2)root:0.5));"); still has 1 (of the 2) root edges
	# writeMultiTopology([n1,n2], stdout)
	n11.root = 2 # below the hybrid node: will trigger RootMismatch and message below
	originalstdout = stdout
	redirect_stdout(open("/dev/null", "w"))
	writeMultiTopology([n11,n12], "test_relaxedreading.net")
	rm("test_relaxedreading.net")
	redirect_stdout(originalstdout)
end
@testset "write for hybrid-lambda" begin
	net = readTopology("((a:1,(bH:1)#H1:1::0.8)H2:5,(#H1:0::0.2,c:1):1);");
	@test hybridlambdaformat(net) == "((a:1.0,(bH:1.0)H1#0.8:1.0)H2:5.0,(H1#0.8:0.0,c:1.0)I1:1.0)I2;"
	net = readTopology("((#H1:::0.2,c),(I3,(b)#H1:::0.8):5);")
	@test hybridlambdaformat(net) == "((H1#0.2,c)I4,(I3,(b)H1#0.2)I5:5.0)I6;"
end
end

ext: jl


text:
#TODO fully implement this!!!
# TODO: vectorize this
# TODO: make vectorized version
"""
    Will add a Met to a given idx. See Met struct for details
"""
function set_rxn!(metnet::MetNet, rxni::Int, rxn::Rxn)
    metnet.rxns[rxni] = rxn.id
    metnet.lb[rxni] = rxn.lb
    metnet.ub[rxni] = rxn.ub
    metnet.c[rxni] = rxn.c
    for (met, s) in zip(rxn.mets, rxn.S)
        meti = metindex(metnet, met)
        metnet.S[meti, rxni] = s
    end
    return metnet
end

set_rxn!(metnet::MetNet, rxn::Rxn) =
    set_rxn!(metnet, findempty(metnet, :rxns; check = true), rxn)

set_rxn!(metnet::MetNet, rxni::Int, id::String) = set_rxn!(metnet, rxni, Rxn(id))
set_rxn!(metnet::MetNet, id::String) =
    set_rxn!(metnet, findempty(metnet, :rxns; check = true), id)


text:
module Fezzik

include("blacklist.jl")
include("import_modules.jl")
include("bootstrap.jl")
include("startup.jl")
include("compile_macro.jl")

end # module


text:
######## Internal Simulation Object Structs ########
abstract type AbstractStage end

mutable struct _Stage{M<:AbstractOperationModel} <: AbstractStage
    key::Int64
    reference::ModelReference
    op_model::Type{M}
    sys::PSY.System
    canonical::CanonicalModel
    optimizer::JuMP.OptimizerFactory
    executions::Int64
    execution_count::Int64
    chronology_ref::Dict{Int64, <:Chronology}
    ini_cond_chron::Union{<:Chronology, Nothing}
    cache::Dict{Type{<:AbstractCache}, AbstractCache}

    function _Stage(key::Int64,
                    reference::ModelReference,
                    op_model::Type{M},
                    sys::PSY.System,
                    canonical::CanonicalModel,
                    optimizer::JuMP.OptimizerFactory,
                    executions::Int64,
                    chronology_ref::Dict{Int64, <:Chronology},
                    cache::Vector{<:AbstractCache}) where M <: AbstractOperationModel

    ini_cond_chron = get(chronology_ref, 0, nothing)
    if !isempty(get_initial_conditions(canonical))
        if isnothing(ini_cond_chron)
            @warn("Initial Conditions chronology set for Stage $(key) which contains Initial conditions")
        end
    end

    pop!(chronology_ref, 0, nothing)

    cache_dict = Dict{Type{<:AbstractCache}, AbstractCache}()
    for c in cache
        cache_dict[typeof(c)] = c
    end

    new{M}(key,
           reference,
           op_model,
           sys,
           canonical,
           optimizer,
           executions,
           0,
           chronology_ref,
           ini_cond_chron,
           cache_dict)

    end

end

######## Exposed Structs to define a Simulation Object ########

mutable struct Stage <: AbstractStage
    op_model::Type{<:AbstractOperationModel}
    model::ModelReference
    execution_count::Int64
    sys::PSY.System
    optimizer::JuMP.OptimizerFactory
    chronology_ref::Dict{Int64, <:Chronology}
    cache::Vector{<:AbstractCache}

    function Stage(::Type{M},
                   model::ModelReference,
                   execution_count::Int64,
                   sys::PSY.System,
                   optimizer::JuMP.OptimizerFactory,
                   chronology_ref=Dict{Int, <:Chronology}(),
                   cache::Vector{<:AbstractCache}=Vector{AbstractCache}()) where M<:AbstractOperationModel

        new(M,
            model,
            execution_count,
            sys,
            optimizer,
            chronology_ref,
            cache)
    end

end

function Stage(::Type{M},
                model::ModelReference,
                execution_count::Int64,
                sys::PSY.System,
                optimizer::JuMP.OptimizerFactory,
                chronology_ref::Dict{Int64, <:Chronology},
                cache::Union{Nothing, AbstractCache}=nothing) where M<:AbstractOperationModel

    cacheinput = isnothing(cache) ? Vector{AbstractCache}() : [cache]
    return Stage(M, model, execution_count, sys, optimizer, chronology_ref, cacheinput)

end

function Stage(model::ModelReference,
               execution_count::Int64,
               sys::PSY.System,
               optimizer::JuMP.OptimizerFactory,
               chronology_ref::Dict{Int64, <:Chronology},
               cache::Union{Nothing, AbstractCache}=nothing)

    return Stage(DefaultOpModel, model, execution_count, sys, optimizer, chronology_ref, cache)

end

get_execution_count(s::S) where S <: AbstractStage = s.execution_count
get_sys(s::S) where S <: AbstractStage = s.sys
get_chronology_ref(s::S) where S <: AbstractStage = s.chronology_ref

get_model_ref(s::Stage) = s.model


text:
@testset "SBitstream" begin
    x = SBitstream(0.5)
    y = SBitstream(0.1)
    w = SBitstream(-0.1)
    X = SBitstream(0.5 * rand(2, 2))
    Y = SBitstream(-0.5 * rand(2, 2))
    v = SBitstream([0.5 * rand(), -0.5 * rand()])

    @testset "FP Values" begin
        @test x.value == 0.5 && y.value == 0.1
        @test float(x) == x.value
        @testset for op in (+, -, *)
            @test op(x, y).value == op(x.value, y.value)
            @test op(x, w).value == op(x.value, w.value)
        end
        @test (y / x).value == y.value / x.value
        @test_logs (:warn, "SBitstream can only be ∈ [-1, 1] (saturation occurring).") (x / y).value == 1
        @test (x ÷ 2).value == x.value / 2
        @test (w ÷ 2).value == w.value / 2
        @test_throws ErrorException x ÷ 0.2
        @test sqrt(x).value == sqrt(x.value)
        @test decorrelate(x).value == x.value
        @test float.(X * Y) == float.(X) * float.(Y)
        @test float.(X * v) == float.(X) * float.(v)
        @test norm(v).value == norm(float.(v))
    end

    @testset "Sample Generation" begin
        xbits = generate(x, 100000)
        @test mean(map(xbit -> pos(xbit) - neg(xbit), xbits)) ≈ x.value atol = 0.01
    end

    @testset "Bit-Level Ops" begin
        T = 30000
        @testset for op in (+, -, *, /, ÷)
            if op == ÷
                @test begin
                    z = 0
                    for t in 1:T
                        bit = pop!(op(x, 2))
                        z += pos(bit) - neg(bit)
                    end

                    z / T
                end ≈ op(x, 2).value atol = 0.01
                @test begin
                    z = 0
                    for t in 1:T
                        bit = pop!(op(w, 2))
                        z += pos(bit) - neg(bit)
                    end

                    z / T
                end ≈ op(w, 2).value atol = 0.01
            elseif op == /
                @test begin
                    z = 0
                    for t in 1:T
                        bit = pop!(op(y, x))
                        z += pos(bit) - neg(bit)
                    end

                    z / T
                end ≈ op(y, x).value atol = 0.01
            else
                @test begin
                    z = 0
                    for t in 1:T
                        bit = pop!(op(x, y))
                        z += pos(bit) - neg(bit)
                    end

                    z / T
                end ≈ op(x, y).value atol = 0.01
                @test begin
                    z = 0
                    for t in 1:T
                        bit = pop!(op(y, x))
                        z += pos(bit) - neg(bit)
                    end

                    z / T
                end ≈ op(y, x).value atol = 0.01
                @test begin
                    z = 0
                    for t in 1:T
                        bit = pop!(op(x, w))
                        z += pos(bit) - neg(bit)
                    end

                    z / T
                end ≈ op(x, w).value atol = 0.01
            end
        end
        @testset for op in (sqrt, decorrelate)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(op(x))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ op(x).value atol = 0.01
        end
        @testset "op = * (matrix, matrix)" begin
            @test begin
                Z = zeros(2, 2)
                for t in 1:T
                    bit = pop!.(X * Y)
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(X * Y) atol = 0.05
            @test begin
                Z = zeros(2, 2)
                for t in 1:T
                    bit = pop!.(Y * X)
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(Y * X) atol = 0.05
        end
        @testset "op = * (matrix, vector)" begin
            @test begin
                T = 50000
                Z = zeros(2)
                for t in 1:T
                    bit = pop!.(X * v)
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(X * v) atol = 0.05
            @test_broken begin
                T = 100000
                Z = zeros(1, 2)
                for t in 1:T
                    bit = pop!.(transpose(v) * X)
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(transpose(v) * X) atol = 0.05
        end
        @testset "op = norm" begin
            @test begin
                T = 50000
                z = 0
                for t in 1:T
                    bit = pop!(norm(v))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ norm(v).value atol = 0.01
        end
    end
end


text:
abstract type Call end

method_expr(f, Ts::Type{<:Tuple}) =
  :($f($([:(::$T) for T in Ts.parameters]...)))

function loc(c::Call)
  meth = method(c)
  "$(meth.file):$(meth.line)"
end

struct DynamicCall{F,A} <: Call
  f::F
  a::A
  DynamicCall{F,A}(f,a...) where {F,A} = new(f, a)
end

DynamicCall(f, a...) = DynamicCall{typeof(f),typeof(a)}(f, a...)

argtypes(c::DynamicCall) = Base.typesof(c.a...)
types(c::DynamicCall) = (typeof(c.f), argtypes(c).parameters...)
method(c::DynamicCall) = which(c.f, argtypes(c))
method_expr(c::DynamicCall) = method_expr(c.f, argtypes(c))

function code(c::DynamicCall; optimize = false)
  codeinfo = @static if VERSION >= v"1.1.0"
    code_typed(c.f, argtypes(c), optimize = optimize, debuginfo=:source)
  else
    code_typed(c.f, argtypes(c), optimize = optimize)
  end
  @assert length(codeinfo) == 1
  codeinfo = codeinfo[1]
  return codeinfo
end

# struct StaticCall <: Call
#   method_instance::MethodInstance
# end
#
# argtypes(c::StaticCall) = Tuple{c.method_instance.specTypes.parameters[2:end]...}
# types(c::StaticCall) = c.method_instance.specTypes
# method(c::StaticCall) = c.method_instance.def
#
# method_expr(c::StaticCall) = method_expr(method(c).name, argtypes(c))
#
# function code(c::StaticCall; optimize = false)
#   # TODO static call graph can only be computed with optimize=true, so analyzing with optimized=false will skip inlined methods
#   codeinfo = get_code_info(c.method_instance, optimize=true)
#   linearize!(codeinfo[1])
#   return codeinfo
# end

function eachline(f, code, line = -1)
  for (i, l) in enumerate(code.code)
    ind = code.codelocs[i]
    1 <= ind <=  length(code.linetable) ?
      line = code.linetable[ind].line :
      line = -1
    f(line, l)
  end
end

exprtype(code, x) = typeof(x)
exprtype(code, x::Core.TypedSlot) = x.typ
exprtype(code, x::QuoteNode) = typeof(x.value)
exprtype(code, x::Core.SSAValue) = code.ssavaluetypes[x.id+1]
# exprtype(code, x::Core.SlotNumber) = code.slottypes[x.id]

rebuild(code, x) = x
rebuild(code, x::QuoteNode) = x.value
rebuild(code, x::Expr) = Expr(x.head, rebuild.(Ref(code), x.args)...)
rebuild(code, x::Core.SlotNumber) = code.slotnames[x.id]

struct Warning
  call::Call
  line::Int
  message::String
  stack::Vector{Call}
end

getmod(w::Warning) = method(w.call).module

Warning(call, message) = Warning(call, -1, message)
Warning(call, line, message) = Warning(call, line, message, Call[])

function warning_printer(modules=[])
  (w) -> begin
    meth = method(w.call)
    if isempty(modules) || getmod(w) in modules
      @safe_warn w.message _file=String(meth.file) _line=w.line _module=nothing
    end
  end
end

# global variables

function globals(warn, call)
  c = code(call)[1]
  eachline(c) do line, ex
    ex isa Expr || return
    for ref in ex.args
      ref isa GlobalRef || continue
      isconst(ref.mod, ref.name) && continue
      getfield(ref.mod, ref.name) isa Module && continue
      warn(call, line, "uses global variable $(ref.mod).$(ref.name)")
    end
  end
end

# fields

function fields(warn, call)
  c = code(call)[1]
  eachline(c) do line, x
    (isexpr(x, :(=)) && isexpr(x.args[2], :call) &&
     rebuild(c, x.args[2].args[1]) == GlobalRef(Core,:getfield)) ||
      return
    x, field = x.args[2].args[2:3]
    x, x_expr, field = exprtype(c, x), rebuild(c, x), rebuild(c, field)
    (isconcretetype(x) && !(x.name.wrapper == Type) && field isa Symbol) ||
      return
    isconcretetype(fieldtype(x, field)) || warn(call, line, "field $x_expr.$field::$(fieldtype(x, field)), $x_expr::$x")
  end
end

# local variables

function assignments(code, l = -1)
  assigns = Dict()
  idx = 0
  eachline(code, l) do line, ex
    idx += 1
    (isexpr(ex, :(=)) && isexpr(ex.args[1], Core.SlotNumber)) || return
    typ = Core.Compiler.widenconst(code.ssavaluetypes[idx])
    push!(get!(assigns, ex.args[1], []), (line, typ))
  end
  return assigns
end

function locals(warn, call)
  c = code(call)[1]
  as = assignments(c)
  for (x, as) in as
    (length(unique(map(x->x[2],as))) == 1 && (isconcretetype(as[1][2]) || istype(as[1][2]))) && continue
    var = c.slotnames[x.id]
    startswith(string(var), '#') && continue
    for (l, t) in as
      warn(call, l, "$(var) is assigned as $(t)")
    end
  end
end

# dynamic dispatch

rebuild(code, x::Core.SSAValue) = rebuild(code, code.code[x.id])

function dispatch(warn, call)
  c = code(call, optimize = true)[1]
  eachline(c) do line, ex
    isexpr(ex, :(=)) && (ex = ex.args[2])
    isexpr(ex, :call) || return
    callex = rebuild(c, ex)
    f = callex.args[1]
    (f isa GlobalRef && isprimitive(getfield(f.mod, f.name)) || isprimitive(f)) && return
    warn(call, line, string("dynamic dispatch to ", callex))
  end
end

# return type

function issmallunion(t)
  ts = Base.uniontypes(t)
  length(ts) == 1 && isconcretetype(first(ts)) && return true
  length(ts) > 2 && return false
  (Missing in ts || Nothing in ts) && return true
  return false
end

istype(::Type{T}) where T = true
istype(::Union) = false
istype(x) = false

function rettype(warn, call)
  c, out = code(call)

  if out == Any || !(issmallunion(out) || isconcretetype(out) || istype(out))
    warn(call, method(call).line, "$(call.f) returns $out")
  end
end

# overall analysis

function analyse(warn, call)
  globals(warn, call)
  locals(warn, call)
  fields(warn, call)
  dispatch(warn, call)
  rettype(warn, call)
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_gthr #
#############

function test_gthr!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_y = CUSPARSE.gthr!(d_x,d_y,'O')
    h_x = to_host(d_x)
    nx = SparseVector(m,x.nzind,y[x.nzind])
    @test_approx_eq(h_x,nx)
end

function test_gthr(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_z = CUSPARSE.gthr(d_x,d_y,'O')
    h_z = to_host(d_z)
    nx = SparseVector(m,x.nzind,y[x.nzind])
    @test_approx_eq(h_z,nx)
end

##############
# test_gthrz #
##############

function test_gthrz!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_x,d_y = CUSPARSE.gthrz!(d_x,d_y,'O')
    h_x = to_host(d_x)
    h_y = to_host(d_y)
    nx = SparseVector(m,x.nzind,y[x.nzind])
    @test_approx_eq(h_x,nx)
    y[x.nzind] = zero(elty)
    @test_approx_eq(h_y,y)
end

function test_gthrz(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_z,d_w = CUSPARSE.gthrz(d_x,d_y,'O')
    h_w = to_host(d_w)
    h_z = to_host(d_z)
    nx = SparseVector(m,x.nzind,y[x.nzind])
    @test_approx_eq(h_z,nx)
    y[x.nzind] = zero(elty)
    @test_approx_eq(h_w,y)
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_gthr!(elty)
    println("gthr! took ", toq(), " for ", elty)
    tic()
    test_gthr(elty)
    println("gthr took ", toq(), " for ", elty)
    tic()
    test_gthrz!(elty)
    println("gthrz! took ", toq(), " for ", elty)
    tic()
    test_gthrz(elty)
    println("gthrz took ", toq(), " for ", elty)
end


text:
# function to calculate the volume of a sphere
function sphere_vol(r)
    # julia allows Unicode names (in UTF-8 encoding)
    # so either "pi" or the symbol π can be used
    return 4 / 3 * pi * r^3
end

# functions can also be defined more succinctly
quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a

# calculates x for 0 = a*x^2+b*x+c, arguments types can be defined in function definitions
function quadratic2(a::Float64, b::Float64, c::Float64)
    # unlike other languages 2a is equivalent to 2*a
    # a^2 is used instead of a**2 or pow(a,2)
    sqr_term = sqrt(b^2 - 4a * c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)
    # multiple values can be returned from a function using tuples
    # if the return keyword is omitted, the last term is returned
    r1, r2
end

unit_volume = sphere_vol(1)

max_stars_repo_name: Sushant-Padha/julia-in-python

max_issues_repo_name: Sushant-Padha/julia-in-python

max_forks_repo_name: Sushant-Padha/julia-in-python


text:
"""
Abstract model type. All models are subtypes of this one, see *e.g.* [`AbstractAModel`](@ref)
"""
abstract type AbstractModel end

"""
Abstract structure to hold all models for object components (*e.g.* leaves, metamers...).
"""
abstract type AbstractComponentModel <: AbstractModel end

"""
Assimilation (photosynthesis) abstract model. All photosynthesis models must be a subtype of
this.
"""
abstract type AbstractAModel <: AbstractModel end

"""
Stomatal conductance abstract model. All stomatal conductance models must be a subtype of
this.

An AbstractGsModel subtype struct must implement at least a g0 field.
"""
abstract type AbstractGsModel <: AbstractModel end

"""
Light interception abstract struct. All light interception models must be a subtype of this.
"""
abstract type AbstractInterceptionModel <: AbstractModel end

"""
Energy balance abstract struct. All energy balance models must be a subtype of this.
"""
abstract type AbstractEnergyModel <: AbstractModel end


text:
module Implicit3DPlotting

export plot_implicit_surface,
       plot_implicit_surface!,
       plot_implicit_curve,
       plot_implicit_curve!,
       GLMakiePlottingLibrary,
       WGLMakiePlottingLibrary

import GLMakie: xlims!, ylims!, zlims!, wireframe!, linesegments!, mesh!, Scene, cam3d!, Point3f0, scatter!, scatter
import GLMakie as GLMakiePlottingLibrary
import WGLMakie as WGLMakiePlottingLibrary
import Meshing: MarchingCubes, MarchingTetrahedra
import GLMakie.GeometryBasics: Mesh, Rect, Vec, decompose, TriangleFace, Point
import Polyhedra: vrep, intersect, polyhedron

"""
Adds an implicitly defined surface to the scene.
"""
function plot_implicit_surface!(
    scene,
    f;
    x_min = -3.0,
    xmin = x_min,
    x_max = 3.0,
    xmax = x_max,
    y_min = xmin,
    ymin = y_min,
    y_max = xmax,
    ymax = y_max,
    z_min = xmin,
    zmin = z_min,
    z_max = xmax,
    zmax = z_max,
    xlims = (xmin, xmax),
    ylims = (ymin, ymax),
    zlims = (zmin, zmax),
    color = :steelblue,
    transparency = true,
    samples=(35,35,35),
    shading = true,
    wireframe=false,
    MarchingModeIsCubes=true,
    WGLMode = false,
)
    if WGLMode
        WGLMakiePlottingLibrary.activate!()
    else
        GLMakiePlottingLibrary.activate!()
    end

    try
        f([1,1,1])
    catch
        throw("The specified function does not take 3 arguments.")
    end
    implicit_mesh = Mesh(f,Rect(Vec(xlims[1], ylims[1],zlims[1]),
                            Vec(xlims[2]-xlims[1], ylims[2]-ylims[1],zlims[2]-zlims[1])), samples=samples, MarchingModeIsCubes ? MarchingCubes() : MarchingTetrahedra())
    if wireframe
        if WGLMode
            WGLMakiePlottingLibrary.wireframe!(scene, implicit_mesh, shading=shading, color=color, transparency=transparency)
        else
            wireframe!(scene, implicit_mesh, shading=shading, color=color, transparency=transparency)
        end
    else
        vertices = decompose(Point{3, Float64}, implicit_mesh)
        triangles = decompose(TriangleFace{Int}, implicit_mesh)
        if WGLMode
            WGLMakiePlottingLibrary.mesh!(scene, vertices, triangles, shading=shading, color=color, transparency=transparency)
        else
            mesh!(scene, vertices, triangles, shading=shading, color=color, transparency=transparency)
        end
    end

    return(scene)
end

"""
Plots an implicitly defined surface.
"""
function plot_implicit_surface(
    f;
    show_axis = true,
    resolution=(800,800),
    scale_plot=false,
    WGLMode=false,
    in_line=false,
    kwargs...
)
    if WGLMode
        WGLMakiePlottingLibrary.activate!()
        WGLMakiePlottingLibrary.AbstractPlotting.inline!(in_line)
        scene = WGLMakiePlottingLibrary.Scene(resolution=resolution, scale_plot=scale_plot, camera=WGLMakiePlottingLibrary.cam3d!, show_axis=show_axis)
    else
        GLMakiePlottingLibrary.activate!()
        GLMakiePlottingLibrary.AbstractPlotting.inline!(in_line)
        scene = GLMakiePlottingLibrary.Scene(resolution=resolution, scale_plot=scale_plot, camera=cam3d!, show_axis=show_axis)
    end
    plot_implicit_surface!(scene, f; WGLMode=WGLMode, kwargs...)
    return(scene)
end

"""
Adds a space curve, implicitly defined by two equations, to a scene.
"""
function plot_implicit_curve!(
    scene,
    f,
    g;
    x_min = -2.0,
    xmin = x_min,
    x_max = 2.0,
    xmax = x_max,
    y_min = xmin,
    ymin = y_min,
    y_max = xmax,
    ymax = y_max,
    z_min = xmin,
    zmin = z_min,
    z_max = xmax,
    zmax = z_max,
    xlims = (xmin, xmax),
    ylims = (ymin, ymax),
    zlims = (zmin, zmax),
    color = :steelblue,
    samples=(30,30,30),
    linewidth=1.5,
    MarchingModeIsCubes=true,
    WGLMode = false,
    kwargs...
)
    if WGLMode
        WGLMakiePlottingLibrary.activate!()
    else
        GLMakiePlottingLibrary.activate!()
    end

    try
        f([1,1,1])
        g([1,1,1])
    catch
        throw("One of the specified function does not take 3 arguments.")
    end
    f_implicit_mesh = Mesh(f, Rect(Vec(xlims[1], ylims[1],zlims[1]),
                            Vec(xlims[2]-xlims[1], ylims[2]-ylims[1],zlims[2]-zlims[1])), samples=samples, MarchingModeIsCubes ? MarchingCubes() : MarchingTetrahedra())
    g_implicit_mesh = Mesh(g, Rect(Vec(xlims[1], ylims[1],zlims[1]),
                            Vec(xlims[2]-xlims[1], ylims[2]-ylims[1],zlims[2]-zlims[1])), samples=samples, MarchingModeIsCubes ? MarchingCubes() : MarchingTetrahedra())

    lines=[]
    for triangle_f in f_implicit_mesh, triangle_g in g_implicit_mesh
        if(min(sum((triangle_f[1]-triangle_g[1]).^2), sum((triangle_f[2]-triangle_g[2]).^2), sum((triangle_f[3]-triangle_g[3]).^2))
                < max(sum((triangle_f[1]-triangle_f[2]).^2), sum((triangle_f[2]-triangle_f[3]).^2), sum((triangle_g[3]-triangle_g[1]).^2), sum((triangle_g[1]-triangle_g[2]).^2), sum((triangle_g[2]-triangle_g[3]).^2), sum((triangle_g[3]-triangle_g[1]).^2)))
            intersection=intersectTwoTriangles(triangle_f, triangle_g)
            if(length(intersection)>=2)
                push!(lines, [Point3f0(intersection[1]), Point3f0(intersection[2])])
            end
        end
    end
    if WGLMode
        # 3*linewidth, as the lines seem to be drawn way thinner than in GLMakie
        try
            foreach(line->WGLMakiePlottingLibrary.linesegments!(scene, line; color=color, linewidth=3*linewidth, kwargs...), lines)
            WGLMakiePlottingLibrary.xlims!(scene, (xlims[1],xlims[2]))
            WGLMakiePlottingLibrary.ylims!(scene, (ylims[1],ylims[2]))
            WGLMakiePlottingLibrary.zlims!(scene, (zlims[1],zlims[2]))
        catch e
            println("No curve in WebGL-Mode detected! Check for the relative generality of the implicit surfaces!")
        end
    else
        try
            foreach(line->linesegments!(scene, line; color=color, linewidth=linewidth, kwargs...), lines)
            GLMakiePlottingLibrary.xlims!(scene, (xlims[1],xlims[2]))
            GLMakiePlottingLibrary.ylims!(scene, (ylims[1],ylims[2]))
            GLMakiePlottingLibrary.zlims!(scene, (zlims[1],zlims[2]))
        catch e
            println("No curve in OpenGL-Mode detected! Check for the relative generality of the implicit surfaces!")
        end
    end
    return(scene)
end

"""
Plots a space curve, implicitly defined by two equations.
"""
function plot_implicit_curve(
    f,
    g;
    resolution=(800,800),
    scale_plot=false,
    show_axis=true,
    in_line=false,
    WGLMode=false,
    kwargs...
)
    if WGLMode
        WGLMakiePlottingLibrary.activate!()
        WGLMakiePlottingLibrary.AbstractPlotting.inline!(in_line)
        scene = WGLMakiePlottingLibrary.Scene(resolution=resolution, scale_plot=scale_plot, camera=WGLMakiePlottingLibrary.cam3d!, show_axis=show_axis)
        plot_implicit_curve!(scene, f, g; WGLMode=WGLMode, kwargs...)
    else
        GLMakiePlottingLibrary.activate!()
        GLMakiePlottingLibrary.AbstractPlotting.inline!(in_line)
        scene = GLMakiePlottingLibrary.Scene(resolution=resolution, scale_plot=scale_plot, camera=cam3d!, show_axis=show_axis)
        plot_implicit_curve!(scene, f, g; WGLMode=WGLMode, kwargs...)
    end
    return(scene)
end

"""
Calculates the intersection of two TriangleFaces.
"""
function intersectTwoTriangles(triangle1, triangle2)
    points1 = [[point[1],point[2],point[3]] for point in triangle1]
    points2 = [[point[1],point[2],point[3]] for point in triangle2]
    polyhedron1 = polyhedron(vrep(points1))
    polyhedron2 = polyhedron(vrep(points2))
    intersection = intersect(polyhedron1, polyhedron2)
    return(vrep(intersection).points.points)
end

end


text:
hi = 61
car = carmichael(hi)

curp = tcnt = 0
print("Carmichael 3 (p×q×r) pseudoprimes, up to p = $hi:")
for j in sortperm(1:size(car)[2], by=x->(car[1,x], car[2,x], car[3,x]))
    p, q, r = car[:, j]
    c = prod(car[:, j])
    if p != curp
        curp = p
        @printf("\n\np = %d\n  ", p)
        tcnt = 0
    end
    if tcnt == 4
        print("\n  ")
        tcnt = 1
    else
        tcnt += 1
    end
    @printf("p× %d × %d = %d  ", q, r, c)
end
println("\n\n", size(car)[2], " results in total.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function GillIntegrator(arg0::jdouble)
    return GillIntegrator((jdouble,), arg0)
end

function get_a(obj::GillIntegrator)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::GillIntegrator)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::GillIntegrator)
    return jcall(obj, "getC", Vector{jdouble}, ())
end


text:
using Documenter, BufferedStreams

makedocs(
    format = :html,
    sitename = "BufferedStreams",
    pages = [
        "Home" => "index.md",
        "Input Streams" => "inputstreams.md",
        "Output Streams" => "outputstreams.md"
    ],
    authors = "D. C. Jones, Ben J. Ward"
)

deploydocs(
    repo = "github.com/BioJulia/BufferedStreams.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
# Autogenerated wrapper script for pylon_julia_wrapper_jll for x86_64-linux-gnu-cxx11-julia_version+1.6.0
export pylon_julia_wrapper

using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("pylon_julia_wrapper")
JLLWrappers.@declare_library_product(pylon_julia_wrapper, "libpylon_julia_wrapper.so")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        pylon_julia_wrapper,
        "lib/libpylon_julia_wrapper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2022-Jan-04: add a short cut function to read the dimension names
#     2022-Jan-28: fix documentation
#
#######################################################################################################################################################################################################
"""

    dimname_nc(file::String)

Return all the names of the dimensions, given
- `file` Path of the netcdf dataset

---
# Examples
```julia
dims = dimname_nc("test.nc");
```
"""
function dimname_nc(file::String)
    _dset = Dataset(file, "r");
    _dims = keys(_dset.dim);
    close(_dset);

    return _dims
end

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2022-Feb-03: abstractize the function to read the variable names recursively
#
#######################################################################################################################################################################################################
"""
This function returns the variable names in the netcdf dataset, and the supported methids are

$(METHODLIST)

"""
function varname_nc end

#######################################################################################################################################################################################################
#
# Changes to the method
# General
#     2022-Jan-04: add a short cut function to read the variable names excluding dimension names
#     2022-Jan-28: fix documentation
#     2022-Feb-03: remove dimension name control
#
#######################################################################################################################################################################################################
"""

    varname_nc(file::String)

Return all the names of the variables (excluding the dimensions), given
- `file` Path of the netcdf dataset

---
# Examples
```julia
vars = varname_nc("test.nc");
```
"""
varname_nc(file::String) = (
    _dset = Dataset(file, "r");
    _vars = varname_nc(_dset);
    close(_dset);

    return _vars
);

#######################################################################################################################################################################################################
#
# Changes to the method
# General
#     2022-Feb-03: add recursive variable query feature
#
#######################################################################################################################################################################################################
"""

    varname_nc(ds::Dataset)

Return all the names of the variables (excluding the dimensions), given
- `ds` NCDatasets.Dataset type dataset

---
# Examples
```julia
dset = Dataset("test.nc");
vars = varname_nc(dset);
close(dset);
```
"""
varname_nc(ds::Dataset) = (
    _vars = String[];

    # read the variables from dataset directly
    _vars = [_vars...; keys(ds)...];

    # loop through the groups
    for _group in keys(ds.group)
        _group_vars = varname_nc(ds.group[_group]);
        _vars = [_vars...; _group_vars...];
    end;

    return _vars
);

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: add a short cut function to read the size of a dataset
#     2022-Jan-28: move the function from size.jl to info.jl
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add recursive variable query feature
#
#######################################################################################################################################################################################################
"""

    size_nc(file::String, var_name::String)

Return the dimensions and size of a NetCDF dataset, given
- `file` Path of the netcdf dataset
- `var_name` Variable name

---
# Examples
```julia
ndims,sizes = size_nc("test.nc", "test");
```
"""
function size_nc(file::String, var_name::String)
    _dset = Dataset(file, "r");

    _fvar = find_variable(_dset, var_name);
    if _fvar === nothing
        close(_dset)
        return error("$(var_name) does not exist in $(file)!");
    end;

    _ndim = ndims(_fvar);
    _size = size(_fvar);
    close(_dset);

    return _ndim, _size
end


text:
"""
    DoesNotExist()

This differential indicates that the derivative Does Not Exist (D.N.E).
This is not the cast that it is not implemented, but rather that it mathematically
is not defined.
"""
struct DoesNotExist <: AbstractDifferential end

function extern(x::DoesNotExist)
    throw(ArgumentError("Derivative does not exit. Cannot be converted to an external type."))
end

Base.Broadcast.broadcastable(::DoesNotExist) = Ref(DoesNotExist())

Base.iterate(x::DoesNotExist) = (x, nothing)
Base.iterate(::DoesNotExist, ::Any) = nothing


text:
#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Loaded Cross Band Pass Filter
# This example is originally from Fig. 7.9 of B. Munk, *Frequency Selective Surfaces,
# Theory and Design,* John Wiley and Sons, 2000.  The same case was analyzed in L. Li,
# D. H. Werner et al, "A Model-Based Parameter Estimation Technique for
# Wide-Band Interpolation of Periodic Moment Method Impedance Matrices With Application to
# Genetic Algorithm Optimization of Frequency Selective Surfaces", *IEEE Trans. AP-S*,
# vol. 54, no. 3, March 2006, pp. 908-924, Fig. 6.  Unfortunately, in neither reference
# are the dimensions of the loaded cross stated, except for the square unit cell
# period of 8.61 mm.  I estimated the dimensions from the sketch in Fig. 6 of the second
# reference.  To provide a reliable comparison, I enlisted my colleague
# [Mike Maybell](https://www.linkedin.com/in/mike-maybell-308b77ba),
# principal of Planet Earth Communications, who generously offered to
# analyze the filter using
# [CST Microwave Studio](https://www.3ds.com/products-services/simulia/products/cst-studio-suite/),
# a rigorous commercial finite volume electromagnetic solver.

# Two identical loaded cross slot-type elements are separated by a 6 mm layer of dielectric
# constant 1.9.  Outboard of each sheet is a 1.1 cm layer of dielectric constant 1.3.
# The closely spaced sheets are a good test of the generalized scattering formulation
# implemented in PSSFSS.  The sheet geometry is shown below.  Remember that the entire
# sheet is metalized *except* for the region of the triangulation.

using Plots, PSSFSS
sheet = loadedcross(class='M', w=0.023, L1=0.8, L2=0.14,
            s1=[0.861,0.0], s2=[0.0,0.861], ntri=600, units=cm)
plot(sheet, unitcell=true)
#md savefig("bpf1.png"); nothing  # hide
#-
#md # ![](bpf1.png)

steering = (ϕ=0, θ=0)
strata = [  Layer()
            Layer(ϵᵣ=1.3, width=1.1cm)
            sheet
            Layer(ϵᵣ=1.9, width=0.6cm)
            sheet
            Layer(ϵᵣ=1.3, width=1.1cm)
            Layer()  ]
flist = 1:0.1:20
results = analyze(strata, flist, steering, resultfile=devnull,
                  logfile=devnull, showprogress=false)
data = extract_result(results, @outputs FGHz s21db(v,v) s11db(v,v))
using DelimitedFiles
dat = readdlm("../src/assets/MaybellLoadedCrossResults.csv", ',', skipstart=1)
p = plot(xlabel="Frequency (GHz)", ylabel="Reflection Coefficient (dB)",
         legend=:left, title="Loaded Cross Band-Pass Filter", xtick=0:2:20, ytick=-30:5:0,
         xlim=(-0.1,20.1), ylim=(-35,0.1))
plot!(p, data[:,1], data[:,3], label="PSSFSS", color=:red)
plot!(p, dat[:,1], dat[:,2], label="CST", color=:blue)
#md savefig("bpf2.png"); nothing  # hide
#-
#md # ![](bpf2.png)
#-
p2 = plot(xlabel="Frequency (GHz)", ylabel="Transmission Coefficient (dB)",
          legend=:bottom, title="Loaded Cross Band-Pass Filter", xtick=0:2:20, ytick=-80:10:0,
         xlim=(-0.1,20.1), ylim=(-80,0.1))
plot!(p2, data[:,1], data[:,2], label="PSSFSS", color=:red)
plot!(p2, dat[:,1], dat[:,4], label="CST", color=:blue)
#md savefig("bpf3.png"); nothing  # hide
#-
#md # ![](bpf3.png)

# This analysis takes about 90 seconds for 191 frequencies on my machine.  Note that
# rather than including two separate invocations of the `loadedcross` function when
# defining the strata, I referenced the same sheet object in the two different locations.
# This allows PSSFSS to recognize that the triangulations are identical, and to exploit
# this fact in making the analysis more efficient.  In fact, if both sheets had been embedded
# in similar dielectric claddings (in the same order), then the GSM (generalized scattering matrix)
# computed for the sheet in its first location could be reused without additional computation for its
# second location.  In this case, though, only the spatial integrals are re-used.  For an oblique
# incidence case, computing the spatial integrals is often the most expensive part of the analysis,
# so the savings from reusing the same sheet definition can be substantial.

# ### Conclusion
# Very good agreement is obtained versus CST over a large dynamic range.


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

using Markdown
using InteractiveUtils

# ╔═╡ f7637354-2aa7-11eb-12d6-c11516eca78b
begin
	using Images
	philip = load("philip.jpg")
end

# ╔═╡ b5d696fa-2aa7-11eb-39ab-cdecffa866ba
url = "https://i.imgur.com/VGPeJ6s.jpg"

# ╔═╡ ecded2d4-2aa7-11eb-0754-f99bb18d3a81
download(url, "philip.jpg")

# ╔═╡ ce58f136-2aa8-11eb-3023-e5eca91fa4da
begin
	(h, w) = size(philip)
	head = philip[(h ÷ 2):h, (w ÷ 10):(9w ÷ 10)]
end

# ╔═╡ 3c39b564-2aa9-11eb-10a5-f16a3356d99b
head

# ╔═╡ 51f03d56-2aa9-11eb-0be6-f139e7839284
[
	head reverse(head, dims=2)
	reverse(head, dims=1) reverse(reverse(head, dims=1), dims=2)
]

# ╔═╡ 5bc99070-2aa9-11eb-051c-1f54b361cd46
new_phil = copy(philip)

# ╔═╡ ac80c394-2aa9-11eb-0003-65caf5a6ba15
red = RGB(1, 0, 0)

# ╔═╡ b73b3b16-2aa9-11eb-1118-35f4b79eb3b3
for i in 1:100
	for j in 1:300
		new_phil[i, j] = red
	end
end

# ╔═╡ d6bd88a4-2aa9-11eb-1009-cb82a3b4856f
new_phil

# ╔═╡ da49cfd2-2aa9-11eb-2b69-bb102cf3f718
begin
	new_phil2 = copy(new_phil)
	new_phil2[100:200, 1:100] .=  RGB(0, 1, 0)
	new_phil2
end

# ╔═╡ 15cb0a44-2aaa-11eb-0512-4b8fc99949ca
function redify(color)
	return RGB(color.r, 0, 0)
end

# ╔═╡ 3007c0aa-2aaa-11eb-1aa2-3d38a0360fa2
begin
	color = RGB(0.5, 0.5, 0.2)
	[color, redify(color)]
end

# ╔═╡ 57301754-2aaa-11eb-22f5-91a2e3031bea
redify.(head)

# ╔═╡ 83f7ac7a-2aaa-11eb-0c55-6be1f6b21921
edge_detect(head)

# ╔═╡ 21844fe2-2aac-11eb-24b4-b7176e263ccb

# ╔═╡ Cell order:
# ╠═b5d696fa-2aa7-11eb-39ab-cdecffa866ba
# ╠═ecded2d4-2aa7-11eb-0754-f99bb18d3a81
# ╠═f7637354-2aa7-11eb-12d6-c11516eca78b
# ╠═ce58f136-2aa8-11eb-3023-e5eca91fa4da
# ╠═3c39b564-2aa9-11eb-10a5-f16a3356d99b
# ╠═51f03d56-2aa9-11eb-0be6-f139e7839284
# ╠═5bc99070-2aa9-11eb-051c-1f54b361cd46
# ╠═ac80c394-2aa9-11eb-0003-65caf5a6ba15
# ╠═b73b3b16-2aa9-11eb-1118-35f4b79eb3b3
# ╠═d6bd88a4-2aa9-11eb-1009-cb82a3b4856f
# ╠═da49cfd2-2aa9-11eb-2b69-bb102cf3f718
# ╠═15cb0a44-2aaa-11eb-0512-4b8fc99949ca
# ╠═3007c0aa-2aaa-11eb-1aa2-3d38a0360fa2
# ╠═57301754-2aaa-11eb-22f5-91a2e3031bea
# ╠═83f7ac7a-2aaa-11eb-0c55-6be1f6b21921
# ╠═21844fe2-2aac-11eb-24b4-b7176e263ccb

max_stars_repo_name: zuckberj/comp-thinking-lectures

max_issues_repo_name: zuckberj/comp-thinking-lectures

max_forks_repo_name: zuckberj/comp-thinking-lectures


text:
module SeDuMi

using SparseArrays
using MATLAB

export sedumi

# The fields should be integer values but the type should be `Float64` to work
# with SeDuMi
mutable struct Cone
    f::Float64 # number of free primal variables / linear dual equality constraints
    l::Float64 # length of LP cone
    q::Vector{Float64} # list of second-order cone dimension
    r::Vector{Float64} # list of rotated second-order cone dimension
    s::Vector{Float64} # list of semidefinite constraints
end
Cone(f::Real, l::Real) = Cone(f, l, Float64[], Float64[], Float64[])
dimension(K::Cone) = K.f + K.l + sum(K.q) + sum(K.r) + sum(K.s.^2)

to_vec(x::Vector) = x
to_vec(x::Float64) = [x]
function to_vec(x::AbstractMatrix)
    @assert isempty(x) || isone(size(x, 2))
    return vec(x)
end

# Solve the primal/dual pair
# min c'x,      max b'y
# s.t. Ax = b,   c - A'x ∈ K
#       x ∈ K
#
# Note, if `A` is square then SeDuMi assumes that `A'` is passed instead,
# see https://github.com/sqlp/sedumi/issues/42
function sedumi(A::Union{Matrix{Float64}, SparseMatrixCSC{Float64}},
                b::Vector{Float64}, c::Vector{Float64},
                K::Cone = Cone(0, size(A, 2));
                kws...)
    pars = Dict{String, Any}(string(key) => value for (key, value) in kws)
    @assert size(A, 1) == length(b)
    @assert size(A, 2) == length(c)
    @assert size(A, 2) == dimension(K)
    # There are 3 output arguments: x, y and info so we use `3` above
    x, y, info = mxcall(:sedumi, 3, A, b, c, K, pars)
    return to_vec(x), to_vec(y), info
end

include("MOI_wrapper.jl")

end # module


text:
export DefaultArray

# To solve this problem, you'll have to improve the `struct` definition and modify/add constructors
struct DefaultArray{T, N} <: AbstractArray{T,N}
    parentarray::Array{T,N}
    defaultvalue::T
end
# DefaultArray(parentarray::AbstractArray{T,N} , defaultvalue) where {N, T} = DefaultArray{eltype(parentarray), ndims(parentarray)}(parentarray, defaultvalue)
DefaultArray(parentarray::Array{T,N}, defaultvalue) where {T, N} = DefaultArray{T, N}(parentarray, convert(T,defaultvalue))
DefaultArray(parentarray::Array{T,N}, defaultvalue::Nothing) where {T, N} = DefaultArray{Union{Nothing, T}, N}(parentarray, convert(Union{Nothing,T},defaultvalue))
# DefaultArray(parentarray::AbstractArray{T, N}, defaultvalue) = DefaultArray{T, N}(parentarray, defaultvalue)
@inline Base.getindex(a::DefaultArray{T,N}, i::Int) where {T, N}   = checkbounds(Bool, a, i) ? a.parentarray[i] : a.defaultvalue
@inline Base.getindex(a::DefaultArray{T,N}, i::Vararg{Int,N}) where {N, T} = checkbounds(Bool, a, i...) ? a.parentarray[i...] : a.defaultvalue
Base.size(a::DefaultArray) = size(a.parentarray)


text:
# nanmean.jl
#
#   Mean, ignoring NaNs
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# 02.04.2020
# Last Edit: 08.04.2021

nanmean(x::Array{T}) where T <: Number = mean(filter(!isnan,x))

nanmean(x::Array{Union{T,Missing}}) where T <: Number = mean(filter(!isnan,filter(!ismissing,x)))

function nanmean(x::Matrix{T},y::Integer) where T Union{S, Union{S, Missing}} where S <: Number
	return mapslices(nanmean,x,dims=y)
end


text:
### Auxiliary functions used as scores for penalising deviation of observed from target acceptance rate

## logistic_rate_score

# logistic_rate_score uses the logistic function to scale the acceptance rate by a factor ranging from 0 to 2
# In other words, it allows to nearly eliminate or double the rate depending on its observed value
# k gives the curve's steepness. For larger k, the curve becomes more steep

logistic_rate_score(x::Real, k::Real=7.) = logistic(x, 2., k, 0., 0.)

## erf_rate_score

# erf_rate_score uses the error function (erf) to scale the acceptance rate by a factor ranging from 0 to 2
# In other words, it allows to nearly eliminate or double the rate depending on its observed value
# k gives the curve's steepness. For larger k, the curve becomes more steep

erf_rate_score(x::Real, k::Real=3.) = erf(k*x)+1

### AcceptanceRateMCTuner

# AcceptanceRateMCTuner tunes empirically on the basis of the discrepancy between observed and target acceptance rate
# This discrepancy is pernalised via a score function set by the user
# The default score function is a stretched logistic map

immutable AcceptanceRateMCTuner <: MCTuner
  targetrate::Real # Target acceptance rate
  score::Function # Score function for penalising discrepancy between observed and target acceptance rate
  period::Integer # Tuning period over which acceptance rate is computed
  verbose::Bool # If verbose=false then the tuner is silent, else it is verbose

  function AcceptanceRateMCTuner(targetrate::Real, score::Function, period::Integer, verbose::Bool)
    @assert 0 < targetrate < 1 "Target acceptance rate should be between 0 and 1"
    @assert period > 0 "Tuning period should be positive"
    new(targetrate, score, period, verbose)
  end
end

AcceptanceRateMCTuner(
  targetrate::Real;
  score::Function=logistic_rate_score,
  period::Integer=100,
  verbose::Bool=false
) =
  AcceptanceRateMCTuner(targetrate, score, period, verbose)

tune!(tune::BasicMCTune, tuner::AcceptanceRateMCTuner) = (tune.step *= tuner.score(tune.rate-tuner.targetrate))

Base.show(io::IO, tuner::AcceptanceRateMCTuner) =
  print(io, "AcceptanceRateMCTuner: target rate = $(tuner.targetrate), period = $(tuner.period), verbose = $(tuner.verbose)")

Base.writemime(io::IO, ::MIME"text/plain", tuner::AcceptanceRateMCTuner) = show(io, tuner)


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module CalculatedABC

using Interpolations
using Plots

include("ABCcurve.jl")
include("ABCanalysis.jl")
include("Gini_coeff.jl")
include("utils.jl")

export
    ABCcurve, ABCanalysis, ABCanalysis4Curve, ABCanalysis4Data,
    gini_coeff,
    remove_small_yields

end # module


text:
module CommunalHelperAdventureHelper

using ..Ahorn, Maple
using Ahorn.CommunalHelper

# Stolen from AdventureHelper
@mapdef Entity "CommunalHelper/AdventureHelper/CustomCrystalHeart" CustomCrystalHeart(
    x::Integer,
    y::Integer,
    color::String="00a81f",
    path::String="",
    nodes::Array{Tuple{Integer,Integer},1}=Tuple{Integer,Integer}[],
)

function getPlacements()
    # Check if there's an adventurehelper folder/zip in the mods folder
    # This is bad and can easily fail, but whatever
    if CommunalHelper.detectMod("adventurehelper")
        return Ahorn.PlacementDict(
            "Crystal Heart (Communal Helper, Adventure Helper)" => Ahorn.EntityPlacement(
                CustomCrystalHeart,
                "point",
                Dict{String,Any}(),
                function (entity)
                    entity.data["nodes"] = [
                        (Int(entity.data["x"]) - 20, Int(entity.data["y"]) - 20),
                        (Int(entity.data["x"]), Int(entity.data["y"]) - 20),
                        (Int(entity.data["x"]) + 20, Int(entity.data["y"]) - 20),
                    ]
                end,
            ),
        )
    else
        @warn "AdventureHelper not detected: CommunalHelper+AdventureHelper plugins not loaded."
        return Ahorn.PlacementDict()
    end
end

placements = getPlacements()

Ahorn.nodeLimits(entity::CustomCrystalHeart) = 1, -1

const shardSprites = Tuple{String,String}[(
        "collectables/CommunalHelper/heartGemShard/shard_outline0$i",
        "collectables/CommunalHelper/heartGemShard/shard0$i",
    ) for i in 0:2
]

function Ahorn.selection(entity::CustomCrystalHeart)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle("collectables/heartGem/3/00", x, y)]

    for i in 1:length(nodes)
        nx, ny = nodes[i]

        push!(res, Ahorn.getSpriteRectangle(shardSprites[mod1(i, 3)][1], nx, ny))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::CustomCrystalHeart)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())

    for i in 1:length(nodes)
        nx, ny = nodes[i]

        Ahorn.drawLines(
            ctx,
            Tuple{Number,Number}[(x, y), (nx, ny)],
            Ahorn.colors.selection_selected_fc,
        )
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CustomCrystalHeart, room::Maple.Room)
    x, y = Ahorn.position(entity)
    path = get(entity.data, "path", "")
    tint = CommunalHelper.hexToRGBA(get(entity.data, "color", "663931"))
    sprite = String

    if path == "heartgem0"
        sprite = "collectables/heartGem/0/00"
    elseif path == "heartgem1"
        sprite = "collectables/heartGem/1/00"
    elseif path == "heartgem2"
        sprite = "collectables/heartGem/2/00"
    elseif path == "heartgem3"
        sprite = "collectables/heartGem/3/00"
    elseif path == ""
        sprite = "collectables/AdventureHelper/RecolorHeart_Outline/00"
        tint = CommunalHelper.hexToRGBA(get(entity.data, "color", "663931"))
        Ahorn.drawSprite(
            ctx,
            "collectables/AdventureHelper/RecolorHeart/00",
            x,
            y,
            tint=tint,
        )
    else
        sprite = "collectables/heartGem/3/00"
    end

    nodes = get(entity.data, "nodes", ())
    for i in 1:length(nodes)
        nx, ny = nodes[i]
        sprIdx = mod1(i, 3)
        Ahorn.drawSprite(ctx, shardSprites[sprIdx][1], nx, ny)
        Ahorn.drawSprite(ctx, shardSprites[sprIdx][2], nx, ny, tint=tint)
    end

    Ahorn.drawSprite(ctx, sprite, x, y)
end

end


text:
#securityId :
# year : "Year"
#data : [{ priority:'1', source: {name:'quandl', id:'quandldataid' }, columns:[x,y,,....z], data:[[]]},
#           { priority:'2', source: {name:'quandl', id:'quandldataid2' }, columns:[x,y,,....z], data:[[]]}

#Before this is executed the we need to generate the Unique ID

function createsecurityid(exchange :: UTF8String ,Ticker :: UTF8String ,securityType :: UTF8String)

end

"""
Just some helper functions starts
"""
function getsecurityid(symbol ::Any , frequency :: Int , exchange ::AbstractString , dataType ::UTF8String)

    mongoQuery = Dict("ticker" => symbol , "dataSources.frequency" => frequency ,
                        "exchange" => exchange , "dataSources.dataType" => dataType)
    doc = find(securityCollection , mongoQuery)
    return get(doc, "securityID", "NULL")
end

function getdataarray(securityID :: AbstractString, year::Int ; priority :: Int = 1)
    dataDoc = find(dataCollection , query("securityID" => securityID ,
                                          "Year" => year,
                                          "data.priority" => priority
                                        )
                )
    dataColoumns = dataDoc["data"]["column"]
    actualData = dataDoc["data"]["column"]["data"]
    return dataColoumns , actualData
end
"""
Just some helper functions ends
"""

#function to gethistory on the basis of two dates

function gethistory(sDate :: AbstractString ,#    Here start date  is greater than end date --- Think of reverse nature of backtester
                    eDate:: AbstractString,
                    dataType::UTF8String,
                    symbols::Array{Any,1},
                    frequency::Int;
                    priority::Int = 1,
                    exchange::AbstractString = "NSE"
                    )

    for symbol in symbols

        securityID = getsecurityid(symbol , frequency , exchange , dataType)
        if securityID == "NULL"
            println("Data not found")
        else
            dataColoumns,actualData = get(securityID , Int(Dates.Year(DateTime(sDate))) , priority)
            if( addData(securityID,dataColoumns , actualData ,date , horizon) == -1 )
                #create Error over here
            else
                println("Data for the symbol " * symbol * " successfully stored")
            end
        end
    end

end

function gethistory(sDate :: AbstractString ,
                    horizon:: Int,
                    dataType::UTF8String,
                    symbols::Array{Any , 1},
                    frequency::Int;
                    priority::Int = 1
                    )
    for symbol in symbols

        securityID = getsecurityid(symbol , frequency , exchange , dataType)
        if securityID == "NULL"
            println("Data not found")
        else
            dataColoumns,actualData = get(securityID , Int(Dates.Year(DateTime(sDate))) , priority)
            if( addData(securityID,dataColoumns , actualData ,date , horizon) == -1 )
                #create Error over here
            else
                println("Data for the symbol " * symbol * " successfully stored")
            end
        end
    end
end

function gethistorybasedonsymbolids(sDate :: AbstractString ,
                                    eDate:: AbstractString,
                                    dataType::UTF8String,
                                    symbolIds::Array{Any , 1},
                                    frequency::Int;
                                    priority::Int = 1
                                    )
    for symbolId in symbolIds

        dataColoumns,actualData = get(securityID , Int(Dates.Year(DateTime(sDate))) , priority)
        if( addData(securityID,dataColoumns , actualData ,sdate , eDate) == -1 )
            #create Error over here
        else
            println("Data for the symbol " * symbol * " successfully stored")
        end
    end
end

function gethistorybasedonsymbolids(sDate :: AbstractString ,
                                    horizon:: Int,
                                    dataType::UTF8String,
                                    symbolIds::Array{Any , 1},
                                    frequency::Int;
                                    priority::Int = 1
                                    )
    for symbolId in symbolIds

        dataColoumns,actualData = get(securityID , Int(Dates.Year(DateTime(sDate))) , priority)
        if( addData( securityID,dataColoumns , actualData ,sDate , horizon) == -1 )
            #create Error over here
        else
            println("Data for the symbol " * symbol * " successfully stored")
        end
    end

end

function getindexofdate(dataColoumns :: Array{Any , 1})

    for i = 1: length(dataColoumns)
        if(dataColoumns[i] == "Date")
            return i
        end
    end
    return -1
end

function addData(securityID :: AbstractString ,
                 dataColoumns :: Array{Any,1},
                 actualData :: Array{Any , 1},
                 sDate :: AbstractString,
                 eDate:: AbstractString)
    sDateObject = DateTime(sDate)
    eDateObject = DateTime(eDate)
    yearDiff = Int(Dates.Year(sDateObject) - Dates.Year(eDateObject))

    dateIndex = getindexofdate(dataColoumns)

    dataToReturn = Array{Any ,1}()
    if yearDiff == 0
        length = length(actualData)
        for i = 1 : length
            presentDate = actualData[i][dateIndex]
            presentDate = DateTime(presentDate)
            if (presentDate <= sDate && presentDate >= eDate)
                append!(dataToReturn , actualData[i])
            end
        end
    else
        #need to fetch multiple year --- difficult implemenation --- need to figure out
        append!(dataToReturn , actualData)
        dateObject = sDateObject - Dates.Year(1)
        yearDiff = yearDiff - 1;
        while yearDiff > 0
            dataColoumnNew , actualDataNew = getdataarray(securityID, Int(Dates.Year(DateTime(dateObject))) , 1)
            append!(dataToReturn , actualDataNew)
            yearDiff = yearDiff - 1;
            dateObject = dateObject - Dates.Year(1)
        end
        yearTobeFeatchedInt = Int(Dates.Year(dateObject)) - 1
        lastDate = string(yearTobeFeatchedInt)*"-12-31"
        dataColoumn , actualData = getdataarray( securityID, yearTobeFeatchedInt)
        dateIndex = getindexofdate(dataColoumn)
        for i = 1 : length(actualData)
            presentDateObject = DateTime(actualData[i][dateIndex])
            if (presentDateObject >= eDate)
                append!(dataToReturn , actualData[i])
            end
        end
    end
    return  dataToReturn
end

function addData(securityID :: AbstractString ,
                 dataColoumns :: Array{Any,1},
                 actualData :: Array{Any , 1},
                 sDate :: AbstractString,
                 horizon:: Int)

    dateTimeObjectSDate = DateTime(sDate)
    dateIndex = getindexofdate(dataColoumns) #can be replaced by lambda function

    """
    Filling the data over here
    """

    dataToReturn = Array{Any , 1}()
    if(dateIndex == -1)
        return -1
    end
    length = length(dataColoumns)
    count = 0
    while count <= horizon
        i = 1
        for i = 1 : length
            presentDateObject = DateTime(actualData[i][dateIndex])
            if(presentDateObject <= dateTimeObjectSDate)
                append!(dataToReturn , actualData[i])
                count = count + 1
            end
        end
        if(count <= horizon)
            yearTobeFeatched = Int(Dates.Year(dateTimeObjectSDate)) - 1

            coloumns , actualData = getdataarray(securityID , Int(Dates.Year(yearTobeFeatched)))
            dateIndex = getindexofdate(coloumns)
            length = length(actualData)
        end
    end
    return dataToReturn
end

max_stars_repo_name: AnkitAggarwalPEC/Quandl-and-Julia

max_issues_repo_name: AnkitAggarwalPEC/Quandl-and-Julia

max_forks_repo_name: AnkitAggarwalPEC/Quandl-and-Julia


text:
 #returns the specification symbol in multicomponent
#returns :singlecomponent if there are no specifications
#throws error if more than one specification is found

const AMOUNT_CONST =Dict{Int,Any}(
     0000 => OneMol()
    ,0001 => MaterialAmount{MOLAR}()
    ,0002 => MaterialAmount{MASS}()
    ,0099 => SingleComponent()
    ,1100 => MaterialCompounds{MOLAR,FRACTION}()
    ,1200 => MaterialCompounds{MASS,FRACTION}()
    ,2100 => MaterialCompounds{MOLAR,TOTAL_AMOUNT}()
    ,2200 => MaterialCompounds{MASS,TOTAL_AMOUNT}()
    ,3100 => HumiditySpec{WetBulbTemperature}()
    ,3200 => HumiditySpec{HumidityRatio}()
    ,3300 => HumiditySpec{MolarHumidity}()
    ,3400 => HumiditySpec{MassHumidity}()
    ,3500 => HumiditySpec{RelativeHumidity}()
    ,3600 => HumiditySpec{HumidityDewPoint}()
)

#the following functions are essencial in the ordering
#do not change the relative order!!

#equilibria properties are always first
_stateorder(x::PhaseFractions) = 1
_stateorder(x::TwoPhaseEquilibrium) = 2
_stateorder(x::VaporFraction) = 3

#pressure is always first in the scalar properties
_stateorder(x::Pressure) = 4
_stateorder(x::Helmholtz) = 5
_stateorder(x::Gibbs) = 6
_stateorder(x::InternalEnergy) = 7
_stateorder(x::VolumeAmount) = 8
_stateorder(x::Enthalpy) = 9
_stateorder(x::Entropy) = 10
#temperature always for last of the scalar properties
_stateorder(x::Temperature) = 11

#From here, it doesnt matter the order, as those are treated specially

_stateorder(x::HumiditySpec) = 1000

_stateorder(x::MaterialCompounds) = 10002

_stateorder(x::MaterialAmount ) = 10001
_stateorder(x::PhaseTag) = 10000
_stateorder(x::Options) = 10004
_stateorder(x::MolecularWeight) = 10005
_stateorder(x::VariableSpec) = 10006

#using specs here allows to delete all reduce_check_mass(::Type)
function _static_specs(x::Type{ThermodynamicState{S,C}}) where {S,C}
    @nospecialize S,C
    tvec = S.parameters
    res1 = [tvec[i].parameters[1]() for i in 1:length(tvec)]
    if C === Tuple{}
        return res1
    else
        tvec2 = C.parameters
        res2 = [tvec2.parameters[i]() for i in 1:length(tvec)]
        append!(res1,res2)
        return res1
    end
end

@generated function static_amount_type(x::T) where T<: ThermodynamicState
    sps = _static_specs(x)
    val = sum(_reduce_check_mass(i) for i in sps)
    modval = mod(val,MATERIAL_SINGLE_MAX)
    if val in (2100,2200)
        _1 = val
        _2 = val
    else
        _1 = max(val-modval,99)
        _2 = modval
    end
    res =  (AMOUNT_CONST[_1],AMOUNT_CONST[_2])

    return :($res)
end

function _state_from_type(x::Type{T}) where T<: ThermodynamicState
    @nospecialize T
    state_specs =  _static_specs(T)
    has_compound = false
    for spec in state_specs
        if spec isa MaterialCompounds
            has_compound = true
            comp_spec = spec
            break
        end
    end

    #double iteration, used in moist air to create an invalid
    #state where compound and humidity spec are present
    #in this case, the resulting type is MultiPT instead of HumPT
    for spec in state_specs
        if (spec isa HumiditySpec) & !has_compound
            has_compound = true
            comp_spec = spec
            break
        end
    end

    if !has_compound
        comp_spec = SingleComponent()
    end
    sorted_specs = sort(state_specs,by= _stateorder)
    res = (sorted_specs[1],sorted_specs[2],comp_spec)
    return res
end

__unwrap(x::Type{T}) where T = x.parameters[1]
__unwrap(x) = x
function _dynamic_get_spec_idx(@nospecialize(x::Type{ThermodynamicState{S,C}}),@nospecialize(sp::Type{S2}))::Int where {S,C,S2}
    state_specs =  _static_specs(x)
    S3 = __unwrap(sp)
    for (i,spec) in pairs(state_specs)
        if typeof(spec) <: S3
            return i
        end
    end
    return 0
end

@generated function _static_get_spec_idx(@nospecialize(x::ThermodynamicState{S,C}),@nospecialize(sp::Type{S2}))::Int where {S,C,S2}
    res = _dynamic_get_spec_idx(x,sp)
    return :($res)
end
"""
    state_type(x::ThermodynamicState)

Returns a sorted tuple of thermodynamic specifications of the input state,
 the tuple has 3 elements, corresponding to two general specifications and a last amount specification,
  indicating if the state is single or multicomponent.

"""
@generated function state_type(x::ThermodynamicState)
    res = _state_from_type(x)
    return :($res)
end

"""
    get_spec(val::AbstractSpec,st::ThermodynamicState)
    get_spec(val::Type{AbstractSpec},st::ThermodynamicState)

Returns the first specification of the type `val` present in the state `st`. returns `nothing` if not found.

## Example

```julia-repl
julia> a1 = state(n = [1.1,2.2,3.3],T=273.15,mol_h = 101.342)
ThermodynamicState with 3 properties:
  Molar amounts : [1.1, 2.2, 3.3][mol]
  Temperature : 273.15[K]
  Molar enthalpy : 101.342[J mol^-1]

julia> get_spec(Enthalpy,a1)
  spec(mol_h = 101.342[J mol^-1])

julia> get_spec(VolumeAmount,a1) == nothing
  true
```

"""
function get_spec end
get_spec(val::T,st::ThermodynamicState) where T<:AbstractSpec = get_spec(T,st)

function get_spec(val::Type{T},st::ThermodynamicState) where T<: AbstractSpec
    i = _static_get_spec_idx(st,T)
    i === 0 && return nothing
    return st.specs[i]
    #return nothing
end

function throw_get_spec(val,st::ThermodynamicState)::Spec
    res = get_spec(val,st)
    res === nothing && throw(error("$st does not have $val"))
    return res
end

"""
    has_spec(val::AbstractSpec,st::ThermodynamicState)::Bool
    has_spec(val::Type{AbstractSpec},st::ThermodynamicState)::Bool

checks if the specification of the type `val` present in the state `st`.

## Example

```julia-repl
julia> a1 = state(n = [1.1,2.2,3.3],T=273.15,mol_h = 101.342)
ThermodynamicState with 3 properties:
  Molar amounts : [1.1, 2.2, 3.3][mol]
  Temperature : 273.15[K]
  Molar enthalpy : 101.342[J mol^-1]

julia> has_spec(Enthalpy,a1)
  true

julia> has_spec(VolumeAmount,a1) == nothing
  false
```

"""
function has_spec end
function has_spec(val::Type{T},st::ThermodynamicState) where T<: AbstractSpec
    static_has_spec(st,val)
end

"""
    amount_type(st::ThermodynamicState)

given a state, returns a tuple with information about how is the amount of matter specified.

The first element of the tuple is related to the number of components present in the state. if it has only one component. it will return `SingleComponent`. if a multicomponent specification (`MaterialCompounds`)
is defined, it will return that type instead.

The second element is related to the type of material specified (molar or mass). if there isn't any specification, it will assume one mol, via the type `OneMol`. It will return a `MaterialAmount`
if it is specified.

## Examples

```julia-repl
julia> tp = state(t=1,p=2)
ThermodynamicState with 2 properties:
  Temperature : 1[K]
  Pressure : 2[Pa]

julia> amount_type(tp)
(SingleComponent(), OneMol())

julia> tpn = state(t=1,p=2,mass=3)
ThermodynamicState with 3 properties:
  Temperature : 1[K]
  Pressure : 2[Pa]
  Mass : 3[kg]

julia> amount_type(tpn)
(SingleComponent(), MaterialAmount{MASS}())

julia> tpxm = state(t=1,p=2,mass=3,xn = [0.1,0.9])
ThermodynamicState with 4 properties:
  Temperature : 1[K]
  Pressure : 2[Pa]
  Mass : 3[kg]
  Molar fraction : [0.1, 0.9]

julia> amount_type(tpxm)
(MaterialCompounds{MOLAR, FRACTION}(), MaterialAmount{MASS}())

julia> tpm = state(t=1,p=2,m = [1,2])
ThermodynamicState with 3 properties:
  Temperature : 1[K]
  Pressure : 2[Pa]
  Mass amounts : [1, 2][kg]

julia> amount_type(tpm)
(MaterialCompounds{MASS, TOTAL_AMOUNT}(), MaterialCompounds{MASS, TOTAL_AMOUNT}())
```
"""
function amount_type(st::T) where T<:ThermodynamicState
    return static_amount_type(st)
end

@generated function static_has_spec(st::ThermodynamicState,val::Type{T}) where T<:AbstractSpec
    idx =_static_get_spec_idx(st,T)
    res = !iszero(idx)
    return :($res)
end

reduce_numtype(::T) where T <: Spec = reduce_numtype(T)
reduce_numtype(::Type{Spec{V,T}}) where {V,T<:Unitful.Quantity{W}} where W = W
reduce_numtype(::Type{Spec{V,T}}) where {V,T<:AbstractVector{<:Unitful.Quantity{W}}} where W = W
reduce_numtype(::Type{Spec{V,T}}) where {V,T<:AbstractVector{<:W}} where W = W
reduce_numtype(::Type{Spec{V,T}}) where {V,T} = Bool
reduce_numtype(::Type{Spec{V,T}}) where {V,T<:Number} = T
reduce_numtype(::Type{T}) where T<:Number = T
reduce_numtype(A,B) = promote_type(reduce_numtype(A),reduce_numtype(B))
function numtype(st::T) where T <: ThermodynamicState
    res = foldl(reduce_numtype,st.specs)
    return res
end
_one(st::T) where T <: ThermodynamicState = one(numtype(st))


text:
module CompressedCellFields

using Gridap
import Gridap: gradient

const CompressedCellFieldLike{D,T,N,R<:FieldLike{D,T,N}} = CompressedCellValue{R}

function gradient(f::CompressedCellFieldLike)
  broadcast(gradient,f)
end

end # module


text:
const Region5_tup = (
    no = [-0.131_799_836_742_01E2,
    0.685_408_416_344_34E1,
   -0.248_051_489_334_66E-1,
    0.369_015_349_803_33,
   -0.311_613_182_139_25E1,
   -0.329_616_265_389_17],

    Jo = [0,
    1,
    -3,
    -2,
    -1,
    2],

    nr = [0.157_364_048_552_59E-2,
    0.901_537_616_739_44E-3,
    -0.502_700_776_776_48E-2,
    0.224_400_374_094_85E-5,
    -0.411_632_754_534_71E-5,
    0.379_194_548_229_55E-7],

    Ir = [1,
    1,
    1,
    2,
    2,
    3],

    Jr = [1,
    2,
    3,
    3,
    9,
    7],
)
function Region5(Output::Symbol, P, T)
    Pstar = 1.0     #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T

    γo    =  log(π) + sum([no[i]*(τ^Jo[i]) for i=1:6])
    γo_π  =  1/π
    γo_ππ = -1/(π^2)
    γo_τ  =          sum([no[i]*Jo[i]*(τ^(Jo[i]-1)) for i=1:6])
    γo_ττ =          sum([no[i]*Jo[i]*(Jo[i]-1)*(τ^(Jo[i]-2)) for i=1:6])
    γo_πτ = 0

    γr    = sum([nr[i]*(π^Ir[i])*(τ^Jr[i]) for i=1:6])
    γr_π  = sum([nr[i]*Ir[i]*(π^(Ir[i]-1))*(τ^Jr[i]) for i=1:6])
    γr_ππ = sum([nr[i]*Ir[i]*(Ir[i]-1)*(π^(Ir[i]-2))*(τ^Jr[i]) for i=1:6])
    γr_τ  = sum([nr[i]*(π^Ir[i])*Jr[i]*(τ^(Jr[i]-1)) for i=1:6])
    γr_ττ = sum([nr[i]*(π^Ir[i])*Jr[i]*(Jr[i]-1)*(τ^(Jr[i]-2)) for i=1:6])
    γr_πτ = sum([nr[i]*Ir[i]*(π^(Ir[i]-1))*Jr[i]*(τ^(Jr[i]-1)) for i=1:6])

    if Output == :SpecificG            #kJ/kg
        return IF97_R*T*(γo + γr)
    elseif Output == :SpecificF            #kJ/kg
        return IF97_R*T*((γo + γr) - π*(γo_π + γr_π)/1000)
    elseif Output == :SpecificV        #m3/kg
        return IF97_R*T*π*(γo_π + γr_π)/P/1000
    elseif Output == :SpecificU        #kJ/kg
        return IF97_R*T*(τ*(γo_τ + γr_τ) - π*(γo_π + γr_π))
    elseif Output == :SpecificS        #kJ/kgK
        return IF97_R*(τ*(γo_τ + γr_τ) - (γo + γr))
    elseif Output == :SpecificH        #kJ/kg
        return IF97_R*T*τ*(γo_τ + γr_τ)
    elseif Output == :SpecificCP       #kJ/kgK
        return -IF97_R*τ^2*(γo_ττ + γr_ττ)
    elseif Output == :SpecificCV       #kJ/kgK
        return IF97_R*((-τ^2)*(γo_ττ + γr_ττ) - (1 + π*γr_π - τ*π*γr_πτ)^2/(1-(π^2)*γr_ππ))
    elseif Output == :SpeedOfSound     #m/s
        return √(1000*IF97_R*T*(1 + 2*π*γr_π + (π^2)*(γr_π^2))/((1 - (π^2)*(γr_ππ))+((1+π*γr_π-τ*π*γr_πτ)^2)/(τ^2*(γo_ττ + γr_ττ))))
    else
        throw(DomainError(Output, "Unknown value requested."))
    end
end

function mass_gibbs_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T
    γo = log(π) + sum(no[i]*τ^Jo[i] for i=1:9)
    γr = sum(nr[i]*π^Ir[i]*(τ-0.5)^Jr[i] for i = 1:6)
    res = IF97_R*T*(γo + γr)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_helmholtz_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr
    π = P / Pstar
    τ = Tstar / T
    γo    = log(π) + sum(no[i]*τ^Jo[i] for i=1:9)
    γr    = sum(nr[i]*π^Ir[i]*(τ-0.5)^Jr[i] for i = 1:6)
    γo_π  = 1/π
    γr_π  = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*(τ-0.5)^Jr[i] for i=1:6)
    res =  IF97_R*T*(γo + γr - π*(γo_π + γr_π)/1000)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_volume_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T

    γo_π  = 1/π
    γr_π  = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*(τ-0.5)^Jr[i] for i=1:6)

    res = IF97_R*T*π*(γo_π + γr_π)/P/1000
    return res
end

function mass_internal_energy_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T
    γo_τ  = sum(no[i]*Jo[i]*τ^(Jo[i]-1) for i=1:9)
    γr_τ  = sum(nr[i]*π^Ir[i]*Jr[i]*(τ-0.5)^(Jr[i]-1) for i=1:6)
    γo_π  = 1/π
    γr_π  = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*(τ-0.5)^Jr[i] for i=1:6)

    res = IF97_R*T*(τ*(γo_τ + γr_τ) - π*(γo_π + γr_π))

    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_entropy_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T
    γo_τ  = sum(no[i]*Jo[i]*τ^(Jo[i]-1) for i=1:9)
    γr_τ  = sum(nr[i]*π^Ir[i]*Jr[i]*(τ-0.5)^(Jr[i]-1) for i=1:6)

    γo    = log(π) + sum(no[i]*τ^Jo[i] for i=1:9)
    γr    = sum(nr[i]*π^Ir[i]*(τ-0.5)^Jr[i] for i = 1:6)

    res = IF97_R*(τ*(γo_τ + γr_τ) - (γo + γr))

    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)
end

function mass_enthalpy_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T
    γo_τ  = sum(no[i]*Jo[i]*τ^(Jo[i]-1) for i=1:9)
    γr_τ  = sum(nr[i]*π^Ir[i]*Jr[i]*(τ-0.5)^(Jr[i]-1) for i=1:6)

    res = IF97_R*T*τ*(γo_τ + γr_τ)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_cp_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
       Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T
    γo_ττ = sum(no[i]*Jo[i]*(Jo[i]-1)*τ^(Jo[i]-2) for i=1:9)
    γr_ττ = sum(nr[i]*π^Ir[i]*Jr[i]*(Jr[i]-1)*(τ-0.5)^(Jr[i]-2) for i=1:6)

    res = -IF97_R*τ^2*(γo_ττ + γr_ττ)

    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)

end

function mass_cv_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T

    γo_ττ = sum(no[i]*Jo[i]*(Jo[i]-1)*τ^(Jo[i]-2) for i=1:9)
    γr_ττ = sum(nr[i]*π^Ir[i]*Jr[i]*(Jr[i]-1)*(τ-0.5)^(Jr[i]-2) for i=1:6)
    γr_π  = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*(τ-0.5)^Jr[i] for i=1:6)
    γr_πτ = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*Jr[i]*(τ-0.5)^(Jr[i]-1) for i = 1:6)
    γr_ππ = sum(nr[i]*Ir[i]*(Ir[i]-1)*π^(Ir[i]-2)*(τ-0.5)^Jr[i] for i=1:6)

    res = IF97_R*(-τ^2*(γo_ττ+γr_ττ) - (1+π*γr_π-τ*π*γr_πτ)^2/(1-π^2*γr_ππ))

    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)
end

function sound_speed_impl(mt::SinglePT,model::IF97Region{:r5},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 1.0    #MPa
    Tstar = 1000.0  #K
    no = Region5_tup.no
    Jo = Region5_tup.Jo
    nr = Region5_tup.nr
    Ir = Region5_tup.Ir
    Jr = Region5_tup.Jr

    π = P / Pstar
    τ = Tstar / T

    γr_π  = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*(τ-0.5)^Jr[i] for i=1:6)
    γr_ππ = sum(nr[i]*Ir[i]*(Ir[i]-1)*π^(Ir[i]-2)*(τ-0.5)^Jr[i] for i=1:6)
    γr_πτ = sum(nr[i]*Ir[i]*π^(Ir[i]-1)*Jr[i]*(τ-0.5)^(Jr[i]-1) for i = 1:6)
    γo_ττ = sum(no[i]*Jo[i]*(Jo[i]-1)*τ^(Jo[i]-2) for i=1:9)
    γr_ττ = sum(nr[i]*π^Ir[i]*Jr[i]*(Jr[i]-1)*(τ-0.5)^(Jr[i]-2) for i=1:6)

    res =  sqrt((1000*IF97_R*T)*(1+2*π*γr_π+π^2*γr_π^2)/
    ((1-π^2*γr_ππ)+(1+π*γr_π-τ*π*γr_πτ)^2/τ^2/(γo_ττ + γr_ττ)))
    return res
end


text:
module GfxBase

export Vec4, Vec3, Vec2, Ray, distance

using Images
using StaticArrays
import LinearAlgebra.cross
import LinearAlgebra.normalize

export RGB

# some type aliases:
const Vec4 = SVector{4, Float32} # 4-vector of floats
const Vec3 = SVector{3, Float32} # 3-vector of floats
const Vec2 = SVector{2, Float32} # 2-vector of floats

struct Ray
    origin::Vec3
    direction::Vec3
end

function distance(a::Vec3, b::Vec3)
    sqrt((b[1] - a[1])^2 + (b[2] - a[2])^2 + (b[3] - a[3])^2)
end

end # module GfxBase

max_stars_repo_name: aedans/raytracer

max_issues_repo_name: aedans/raytracer

max_forks_repo_name: aedans/raytracer


text:
left(rect::Rect2D) = minimum(rect)[1]
right(rect::Rect2D) = maximum(rect)[1]
bottom(rect::Rect2D) = minimum(rect)[2]
top(rect::Rect2D) = maximum(rect)[2]

bottomleft(bbox::Rect2D) = Point(left(bbox), bottom(bbox))
topleft(bbox::Rect2D) = Point(left(bbox), top(bbox))
bottomright(bbox::Rect2D) = Point(right(bbox), bottom(bbox))
topright(bbox::Rect2D) = Point(right(bbox), top(bbox))

topline(bbox::Rect2D) = (topleft(bbox), topright(bbox))
bottomline(bbox::Rect2D) = (bottomleft(bbox), bottomright(bbox))
leftline(bbox::Rect2D) = (bottomleft(bbox), topleft(bbox))
rightline(bbox::Rect2D) = (bottomright(bbox), topright(bbox))

function shrinkbymargin(rect, margin)
    return IRect(minimum(rect) .+ margin, (widths(rect) .- 2 .* margin))
end

function limits(r::Rect{N, T}) where {N, T}
    mini, maxi = extrema(r)
    return ntuple(i -> (mini[i], maxi[i]), N)
end

function limits(r::Rect, dim::Integer)
    return (minimum(r)[dim], maximum(r)[dim])
end

xlimits(r::Rect) = limits(r, 1)
ylimits(r::Rect) = limits(r, 2)

function enlarge(bbox::Rect2D, l, r, b, t)
    BBox(left(bbox) - l, right(bbox) + r, bottom(bbox) - b, top(bbox) + t)
end

function center(bbox::Rect2D)
    Point2((right(bbox) + left(bbox)) / 2, (top(bbox) + bottom(bbox)) / 2)
end

"""
Converts a point in fractions of rect dimensions into real coordinates.
"""
function fractionpoint(bbox::FRect2D, point::T) where T <: Point2
    T(left(bbox) + point[1] * width(bbox), bottom(bbox) + point[2] * height(bbox))
end

function anglepoint(center::Point2, angle::Real, radius::Real)
    Ref(center) .+ Ref(Point2(cos(angle), sin(angle))) .* radius
end


text:
function notdecreasing(s)
    s[1] <= s[2] <= s[3] <= s[4] <= s[5] <= s[6]
end

function adjacentsame(s)
    for i=1:5
        if s[i] == s[i+1]
            return true
        end
    end
    false
end

function validpassword(password)
    s = string(password)
    notdecreasing(s) && adjacentsame(s)
end

puzzle_range = 168630:718098

println(sum(map(validpassword, puzzle_range)))

function singleduplet(s)
    for i=1:5
        if s[i] == s[i+1]
            if (i > 1) && (s[i-1] == s[i])
                continue
            end
            if (i < 5) && (s[i+1] == s[i+2])
                continue
            end
            return true
        end
    end
    false
end

function validpassword2(password)
    s = string(password)
    notdecreasing(s) && adjacentsame(s) && singleduplet(s)
end

println(sum(map(validpassword2, puzzle_range)))

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
#!/usr/bin/env julia

# Lets add workers programatically, and see if they distribute to different cores.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--cpus"
    help = "Which CPUS to add"
    arg_type = Array{Int, 1}
end

parsed_args = parse_args(ARGS, s)

cpus = parsed_args["cpus"]
p = length(cpus)

using ClusterManagers

addprocs(LocalAffinityManager(;affinities=cpus))

nchild = nworkers()
println("Workers allocated ", nchild)
println("Total workers ", workers())

println("Total processes ", nprocs())

@everywhere function print_cpus()
    run(pipeline(`ps -eo pid,psr,pcpu,cmd`, `grep julia`))
end

println("Before map")
print_cpus()

@everywhere function pcpus(mat)
    println("In map")
    print_cpus()
    svd(mat)
end

M = [rand(1000,1000) for i=1:p]

pmap(pcpus, M)

println("After map.")
print_cpus()


text:
#####################
# Timing and retiming
#####################

function normalize_start(x::TimeSeries, y...)
  mutate_time(x->normalize_start(x, y...), x)
end
normalize_start(time, start = time[1]) = map(x -> x - start, time)

# intersect general time series
const default_num_samples = 100
function intersect_intervals(series; num_samples = default_num_samples)
  lower = maximum(get_time(x)[1] for x in series if length(get_time(x)) > 0)
  upper = minimum(get_time(x)[end] for x in series if length(get_time(x)) > 0)

  range(lower, stop=upper, length=num_samples)
end

# intersect regular time series (such as produced using published iterations or
# with a auto-generated linear index)
function intersect_regular(series::AbstractArray{<:AbstractArray})
  sort(reduce(intersect, series[2:end], init=series[1]))
end

# Interpolate a deconstructed time series at a single point
# ts: sorted timestamps for data
# return tuple of value and time for lower bound on time
#
# This is a special internal version that also takes and returns a lower bound
# for efficient linear indexing
function interpolate_(sample_time::Real, ts, data, lower_bound = 0)
  if sample_time <= ts[1]
    return (data[1,:], 0)
  elseif sample_time >= ts[end]
    return (data[length(ts),:], length(ts))
  end

  for ii = lower_bound+1:length(ts)
    if sample_time < ts[ii]
      dl = data[ii-1,:]
      tl = ts[ii-1]

      du = data[ii,:]
      tu = ts[ii]

      value = dl + (du - dl) * (sample_time - tl) / (tu - tl)

      return (value, ii - 1)
    end
  end

  # Control flow should not be able to get here
  throw(ArgumentError("Interpolate escaped bounds"))
end

function interpolate_(sample_time::Real, series::TimeSeries, x...)
  interpolate_(sample_time, get_time(series), get_data(series), x...)
end

interpolate(sample_time::Real, x...) = interpolate_(sample_time, x...)[1]

# hack to get type of array to allocate after division
function division_output_type(x::Type)
  types = Base.return_types(/, Tuple{x, Float64})
  promote_type(types...)
end

# Interpolate a deconstructed time series at a series of points
# sample_times should be sorted
function interpolate(sample_times::AbstractArray{<:Real,1}, ts, data)
  assert_sorted(sample_times)

  output_dimension = (length(sample_times), size(data)[2:end]...)
  output_element_type = division_output_type(eltype(data))

  ret = similar(data, output_element_type, output_dimension)
  lower_bound = 0

  for ii = 1:size(ret, 1)
    (ret[ii,:], lower_bound) = interpolate_(sample_times[ii], ts, data,
                                          lower_bound)
  end

  ret
end

# Interpolate a TimeSeries at one or more points
interpolate(sample_time::Real, x::TimeSeries) =
  interpolate(sample_time, get_time(x), get_data(x))
function interpolate(sample_times::AbstractArray, x::TimeSeries)
  data = interpolate(sample_times, get_time(x), get_data(x))

  TimeSeries(sample_times, data)
end

# Access a time series at an exact time
#
# This is a special internal version that also takes and returns a lower bound
# for efficient linear indexing
function get_at_time_(time::Real, x::TimeSeries, lower_bound = 1)
  for ii = lower_bound:size(x, 1)
    if get_time(x)[ii] == time
      return (get_data(x)[ii,:], ii)
    end
  end

  error("Attempt to get at invalid time: Try interpolate instead?")
end

# Access a time series at a given time or array of times
get_at_time(time::Real, x::TimeSeries) = get_at_time_(time, x)[1]
function get_at_time(times::AbstractArray{<:Real,1}, x::TimeSeries)
  assert_sorted(times)

  output_dimension = (length(times), size(x)[2:end]...)
  ret = similar(get_data(x), output_dimension)
  lower_bound = 1

  for ii = 1:size(ret, 1)
    (ret[ii,:], lower_bound) = get_at_time_(times[ii], x, lower_bound)
  end

  TimeSeries(times, ret)
end

function get_at_nearest_time(time::Real, x::TimeSeries)
  ii = searchsortedlast(get_time(x), time)
  jj = searchsortedfirst(get_time(x), time)

  output_index = 0
  best_error = Inf

  for index in (ii, jj)
    if 0 < index <= size(x,1)
      error = abs(time - get_time(x)[index])
      if error < best_error
        best_error = error
        output_index = index
      end
    end
  end

  get_data(x)[output_index]
end

# Intersect intervals (and resample) for a set of time series and return a
# concatenated time series
function intersect_interpolate(series::AbstractArray{<:TimeSeries}; x...)
  # compute the interval
  interval = intersect_intervals(series; x...)

  cat_dim = maximum(ndims, series) + 1

  # interpolate and concatenate data
  cat(dims=cat_dim, map(x->interpolate(interval, x), series)...)
end

# Exactly intersect intervals  for a set of time series and return a
# concatenated time series
function intersect_series(series::AbstractArray{<:TimeSeries})
  time = intersect_regular(map(get_time, series))

  cat_dim = maximum(ndims, series) + 1

  cat(dims=cat_dim, map(x->get_at_time(time, x), series)...)
end

# Returns a view into x at given indices
function select_indices(inds, x::TimeSeries)
  TimeSeries(selectdim(get_time(x), 1, inds), selectdim(get_data(x), 1, inds))
end

# Returns a view into x for times where f is true
function select_times(f, x::TimeSeries)
  inds = [ind for (ind, t) in enumerate(get_time(x)) if f(t)]
  select_indices(inds, x)
end

##########
# Calculus
##########

# The following code relates to numeric integration and differentiation of time
# series data

# differentiation algorithms
struct Windowed end; struct Discrete end

# default algorithm
function differentiate(x::TimeSeries, method=Windowed(); kws...)
  differentiate(x, method=method; kws...)
end

# concrete implementations

# Compute the derivative over a window. Interpolate for smoothness
function differentiate_windowed(x::TimeSeries, lower_time, upper_time,
                                lower_lb = 0, upper_lb = 0)
  lower_time = max(lower_time, get_time(x)[1])
  upper_time = min(upper_time, get_time(x)[end])

  (lower_value, lower_lb) = interpolate_(lower_time, x, lower_lb)
  (upper_value, upper_lb) = interpolate_(upper_time, x, upper_lb)

  derivative = (upper_value - lower_value) / (upper_time - lower_time)

  (derivative, lower_lb, upper_lb)
end

# Differentiate the TimeSeries using a fixed time window
# I could potentially implement this instead by interpolating arrays rather than
# value. That would be slightly slower but also cleaner.
function differentiate(x::TimeSeries, method::Windowed;
                       time=get_time(x), window=(time[end]-time[1])/(length(time)-1))
  assert_sorted(time)

  output_dimension = (length(time), size(x)[2:end]...)
  ret = TimeSeries(time, similar(get_data(x), output_dimension))

  lower_lb = 0
  upper_lb = 0

  hw = window / 2

  for (ii, t) = enumerate(time)
    (ret[ii,:], lower_lb, upper_lb) = differentiate_windowed(x, t-hw, t+hw,
                                                             lower_lb, upper_lb)
  end

  ret
end

############
# Statistics
############

# Compute standard error in a given dimension.
# Default dimension is two as samples are typically in the second dimension for
# TimeSeries.
function standard_error(x, dim = 2)
  std(x, dims=dim) ./ sqrt(size(x, dim))
end

#######
# Other
#######

# Returns first point to cros y = value
function first_horizontal_intersection(x::TimeSeries{N,1}, value) where {N}
  for ii = 1:size(x, 1)-1
    y1 = x[ii]
    y2 = x[ii+1]

    # Test for crossing
    if y1 <= value <= y2 || y2 <= value <= y1
      x1 = get_time(x)[ii]
      x2 = get_time(x)[ii+1]

      return x1 + (x2 - x1) * (value - y1) / (y2 - y1)
    end
  end

  # If it doesn't reach the value, return infinity
  # (Assuming a monotonic input which reaches at infinity)
  return Inf
end

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
abstract type CalciumScore end

"""
    function Integrated(
        vol,
        I=sum(vol),
        N=length(vol)
        )

# Arguments
- vol: region of interest
- I: integrated intensity of the entire `vol`
- N: total number of voxels contained in the entire `vol`

# Citation
'Accurate quantification of vessel cross-sectional area using CT
angiography: a simulation study' [Molloi, Johnson, Ding, Lipinski]
(DOI: 10.1007/s10554-016-1007-9)
"""
struct Integrated{T1<:AbstractArray,T2<:AbstractFloat} <: CalciumScore
    vol::T1
    I::T2
    N::T2
end

function Integrated(
    vol::AbstractArray,
    I::AbstractFloat=Float64.(sum(vol)),
    N::AbstractFloat=Float64.(length(vol)),
)
    return Integrated(I, N)
end

"""
    score(vol, S_Bkg, S_Obj, algorithm::Integrated)
    score(vol, S_Bkg, S_Obj, size, algorithm::Integrated)
    score(vol, S_Bkg, S_Obj, size, ρ, algorithm::Integrated)

Given a volume `vol` of interest, use the integrated intensity
approach to find the total number of object voxels/volume/mass
contained in `vol`

# Arguments
- vol: region of interest
- S_Bkg: pure background signal intensity in the `vol`
- S_Obj: pure object signal intensity in the `vol`
- size: size of the voxels in the given `vol`
- ρ: density of the given object of interest contained within
    the `vol`
- algorithm: the `Integrated` algorithm which accounts for
    noise and partial volume effect by integrating the
    the intensity of the entire volume `vol`

# Citation
'Accurate quantification of vessel cross-sectional area using CT
angiography: a simulation study' [Molloi, Johnson, Ding, Lipinski]
(DOI: 10.1007/s10554-016-1007-9)
"""
function score(S_Bkg, S_Obj, algorithm::Integrated)
    I = algorithm.I
    N = algorithm.N

    N_Obj = (I - (N * S_Bkg)) / (S_Obj - S_Bkg)
    return N_Obj
end

function score(S_Bkg, S_Obj, size, algorithm::Integrated)
    I = algorithm.I
    N = algorithm.N

    N_Obj = (I - (N * S_Bkg)) / (S_Obj - S_Bkg)
    Vol_Obj = N_Obj * size[1] * size[2] * size[3]
    return Vol_Obj
end

function score(S_Bkg, S_Obj, size, ρ, algorithm::Integrated)
    I = algorithm.I
    N = algorithm.N

    N_Obj = (I - (N * S_Bkg)) / (S_Obj - S_Bkg)
    Vol_Obj = N_Obj * size[1] * size[2] * size[3]
    Mass_Obj = Vol_Obj * ρ
    return Mass_Obj
end


text:
# lets see if we can interpret the data

using Caesar, RoME, IncrementalInference
using CloudGraphs, Neo4j
# using Distributions
using JSON

using KernelDensityEstimate

using Base.Test

## Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true)
# session = addrdict["session"]
# Nparticles = parse(Int,addrdict["num particles"])

# interactive operation
session = "SESSTRUT45"
Nparticles = 100
include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

fg = Caesar.initfg(sessionname=session, cloudgraph=cloudGraph)

updatenewverts!(fg, N=Nparticles)

ls(fg,:x754)

writeGraphPdf(fg)
@async run(`evince fg.pdf`)

setDBAllReady!(fg)

# totally reset to frontend built state in DB
resetentireremotesession(conn,session)

fg.fIDs

sortedvd = getnewvertdict(fg.cg.neo4j.connection, fg.sessionname)

@show collect(keys(sortedvd))

populatenewvariablenodes!(fg, sortedvd, N=N)

populatenewfactornodes!(fg, sortedvd)

ls(fg)

sortedvd[200050]

loadtx = transaction(conn)
query = "match (n:$(session))-[:DEPENDENCE]-(f:NEWDATA:$(session):FACTOR) where n.ready=1 or f.ready=1 return distinct n, f"
cph = loadtx(query, submit=true)

#
# sortedvd[110240]
#
elem = sortedvd[110592][:frtend]
#
# elem["range"]
#
rn = recoverConstraintType(elem)

# pty = convert(PackedPose2DPoint2DBearingRangeDensity, rn)
# upt = convert(Pose2DPoint2DBearingRangeDensity, pty)

# @test norm(getPoints(rn.bearing)-getPoints(upt.bearing)) < 1e-10
# @test norm(getPoints(rn.range)-getPoints(upt.range)) < 1e-10

# plotKDE(upt.bearing)
plotKDE(rn.bearing)

vc = startdefaultvisualization(draworigin=true)

drawposepoints!(vc, fg, :x8)

getVal(fg, :l200050, api=dlapi)

visualizeallposes!(vc, fg, drawlandms=false, drawtype=:fit)

# for lb in [:x1,:x2,:x3,:x4,:x5,:x6,:x7]
#   @show getKDEMax(getVertKDE(fg, lb))
# end
plotKDE(marginal(getKDE(getVert(fg,:l200050,api=dlapi)),[1;2]))

#