text:
#=
Low-level wrapper for PCRE2 library

Copyright 2018 Gandalf Software, Inc., Scott P. Jones, and contributors to pcre.jl and pcre2.h
(Based in part on julia/base/pcre.jl, and on pcre2.h (copyright University of Cambridge))
Licensed under MIT License, see LICENSE.md
=#

__precompile__()
module PCRE2

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("PCRE2 not installed properly, run Pkg.build(\"PCRE2\"), restart Julia and try again")
end
include(depsjl_path)

PCRE_LOCK = nothing

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
    @static if isdefined(Base.PCRE, :PCRE_COMPILE_LOCK)
        global PCRE_LOCK = Base.PCRE.PCRE_COMPILE_LOCK
    else
        global PCRE_LOCK = Threads.SpinLock()
    end
end

const CodeUnitTypes = Union{UInt8, UInt16, UInt32}

import Base.GC: @preserve
evr(str, rep, sub) = Core.eval(@__MODULE__, Meta.parse(replace(str, rep => sub)))
const _ncodeunits = ncodeunits
create_vector(T, len)  = Vector{T}(undef, len)

import Base: RefValue

include("pcre2_h.jl")
include("pcre2_api.jl")

# Create Julia bindings to PCRE2 C API

for siz in (8,16,32), (nam, ret, sig) in funclist
    l = SubString("a,b,c,d,e,f,g,h,i,j,k,l,m", 1, length(sig)*2-1)
    #parms = string(["$('a'+i-1)::$(sig[i]), " for i=1:length(sig)]...)[1:end-1]
    sub = "UInt$siz"
    str = "$nam(::Type{$sub},$l)=ccall((:pcre2_$(nam)_$siz, libpcre2_$siz),$ret,$sig,$l)"
    evr(str, "PCRE2._UCHAR", sub)
end

const UNSET = ~Csize_t(0)  # Indicates that an output vector element is unset

struct PCRE2_Error <: Exception
    errmsg::String
    errno::Int32
    erroff::Int32
    PCRE2_Error(errmsg, errno, erroff = -1) = new(errmsg, errno%Int32, erroff%Int32)
end

Base.show(io::IO, exc::PCRE2_Error) =
    print(io, "PCRE2: ", exc.errmsg,
          exc.errno < 0 ? "" : "error: $(err_message(exc.errno))",
          exc.erroff < 0 ? "" : "at offset $(exc.erroff)")

pcre_error(msg::AbstractString) = throw(PCRE2_Error(msg, -1))
pcre_error(errno::Integer)      = throw(PCRE2_Error("", errno))
jit_error(errno::Integer)       = throw(PCRE2_Error("JIT ", errno))
compile_error(errno, erroff)    = throw(PCRE2_Error("compilation ", errno, erroff))

function info_error(errno)
    errno == Base.PCRE.ERROR_NULL && pcre_error("NULL regex object")
    errno == Base.PCRE.ERROR_BADMAGIC && pcre_error("Invalid regex object")
    errno == Base.PCRE.ERROR_BADOPTION && pcre_error("Invalid option flags")
    pcre_error(errno)
end

"""PCRE2 pattern_info call wrapper"""
function info(::Type{T}, regex::CodeP, what::INFO, ::Type{S}) where {S,T<:CodeUnitTypes}
    buf = RefValue{S}()
    ret = pattern_info(T, regex, what, buf)
    ret == 0 ? buf[] : info_error(ret)
end

get_ovec(::Type{T}, md) where {T<:CodeUnitTypes} =
    unsafe_wrap(Array, get_ovector_pointer(T, md), 2 * get_ovector_count(T, md))

"""Wrapper for compile() function, throw error on error return with error info"""
function compile(pattern::T, options::Integer) where {T<:AbstractString}
    errno = RefValue{Cint}(0)
    erroff = RefValue{Csize_t}(0)
    re_ptr = compile(codeunit(T), pattern, _ncodeunits(pattern), options, errno, erroff, C_NULL)
    re_ptr == C_NULL ? compile_error(errno[], erroff[]) : re_ptr
end

"""Wrapper for jit_compile() function, throw error on error return with error info"""
jit_compile(::Type{T}, regex::CodeP) where {T<:CodeUnitTypes} =
    ((errno = jit_compile(T, regex, JIT_COMPLETE)) == 0 || jit_error(errno) ; nothing)

function err_message(errno)
    buffer = create_vector(UInt8, 256)
    get_error_message(UInt8, errno, buffer, sizeof(buffer))
    @preserve buffer unsafe_string(pointer(buffer))
end

function sub_length_bynumber(::Type{T}, match_data, num) where {T<:CodeUnitTypes}
    s = RefValue{Csize_t}()
    rc = substring_length_bynumber(T, match_data, num, s)
    rc < 0 ? pcre_error(rc) : convert(Int, s[])
end

function sub_copy_bynumber(::Type{T}, match_data, num, buf, siz::Integer) where {T<:CodeUnitTypes}
    s = RefValue{Csize_t}(siz)
    rc = substring_copy_bynumber(T, match_data, num, buf, s)
    rc < 0 ? pcre_error(rc) : convert(Int, s[])
end

function capture_names(T, re)
    name_count = info(T, re, INFO_NAMECOUNT, UInt32)
    name_entry_size = info(T, re, INFO_NAMEENTRYSIZE, UInt32)
    nametable_ptr = info(T, re, INFO_NAMETABLE, Ptr{UInt8})
    names = Dict{Int, String}()
    for i=1:name_count
        offset = (i-1)*name_entry_size + 1
        # The capture group index corresponding to name 'i' is stored as a
        # big-endian 16-bit value.
        high_byte = UInt16(unsafe_load(nametable_ptr, offset))
        low_byte = UInt16(unsafe_load(nametable_ptr, offset+1))
        idx = (high_byte << 8) | low_byte
        # The capture group name is a null-terminated string located directly
        # after the index.
        names[idx] = unsafe_string(nametable_ptr+offset+1)
    end
    names
end

end # module PCRE2


text:
src_dir = "Queller"
output_dir = "QuellerCLI"
build_env = "julia_build_env"

using Pkg
Pkg.activate(build_env)
Pkg.instantiate()

# Create build dir
mkpath(output_dir)

# Build program
using PackageCompiler
create_app(src_dir, output_dir, executables=["Queller"=>"main"], force=true)

# Create a shortcuts to program
shortcut = joinpath(output_dir, "QuellerCLI")

if Sys.iswindows()
	exe = joinpath("bin", "Queller.exe")
	shortcut *= ".lnk"

	create_sc = "\$s=(New-Object -COM WScript.Shell).CreateShortcut('$shortcut');"
	set_sc_target = "\$s.TargetPath='%windir%\\explorer.exe';"
	set_sc_argument = "\$s.Arguments='$exe';"
	set_sc_icon = "\$s.IconLocation='Winver.exe,0';"
	save_sc = "\$s.Save()"
	cmd = `powershell "$create_sc $set_sc_target $set_sc_argument $set_sc_icon $save_sc"`

	run(cmd)
else
	symlink("bin/Queller", shortcut)
end

ext: jl

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# Test conversion to and from unix
@test Dates.unix2datetime(Dates.datetime2unix(DateTime(2000,1,1))) == DateTime(2000,1,1)
@test Dates.value(Dates.DateTime(1970)) == Dates.UNIXEPOCH

# Tests from here: http://en.wikipedia.org/wiki/Unix_time
@test string(Dates.unix2datetime(1095379198.75)) == string("2004-09-16T23:59:58.75")
@test string(Dates.unix2datetime(1095379199.00)) == string("2004-09-16T23:59:59")
@test string(Dates.unix2datetime(1095379199.25)) == string("2004-09-16T23:59:59.25")
@test string(Dates.unix2datetime(1095379199.50)) == string("2004-09-16T23:59:59.5")
@test string(Dates.unix2datetime(1095379199.75)) == string("2004-09-16T23:59:59.75")
@test string(Dates.unix2datetime(1095379200.00)) == string("2004-09-17T00:00:00")
@test string(Dates.unix2datetime(1095379200.25)) == string("2004-09-17T00:00:00.25")
@test string(Dates.unix2datetime(1095379200.50)) == string("2004-09-17T00:00:00.5")
@test string(Dates.unix2datetime(1095379200.75)) == string("2004-09-17T00:00:00.75")
@test string(Dates.unix2datetime(1095379201.00)) == string("2004-09-17T00:00:01")
@test string(Dates.unix2datetime(1095379201.25)) == string("2004-09-17T00:00:01.25")
@test string(Dates.unix2datetime(915148798.75)) == string("1998-12-31T23:59:58.75")
@test string(Dates.unix2datetime(915148799.00)) == string("1998-12-31T23:59:59")
@test string(Dates.unix2datetime(915148799.25)) == string("1998-12-31T23:59:59.25")
@test string(Dates.unix2datetime(915148799.50)) == string("1998-12-31T23:59:59.5")
@test string(Dates.unix2datetime(915148799.75)) == string("1998-12-31T23:59:59.75")
@test string(Dates.unix2datetime(915148800.00)) == string("1999-01-01T00:00:00")
@test string(Dates.unix2datetime(915148800.25)) == string("1999-01-01T00:00:00.25")
@test string(Dates.unix2datetime(915148800.50)) == string("1999-01-01T00:00:00.5")
@test string(Dates.unix2datetime(915148800.75)) == string("1999-01-01T00:00:00.75")
@test string(Dates.unix2datetime(915148801.00)) == string("1999-01-01T00:00:01")
@test string(Dates.unix2datetime(915148801.25)) == string("1999-01-01T00:00:01.25")

@test Date(Dates.rata2datetime(734869)) == Dates.Date(2013,1,1)
@test Dates.datetime2rata(Dates.rata2datetime(734869)) == 734869

# Tests from here: http://mysite.verizon.net/aesir_research/date/back.htm#JDN
@test Dates.julian2datetime(1721119.5) == Dates.DateTime(0,3,1)
@test Dates.julian2datetime(1721424.5) == Dates.DateTime(0,12,31)
@test Dates.julian2datetime(1721425.5) == Dates.DateTime(1,1,1)
@test Dates.julian2datetime(2299149.5) == Dates.DateTime(1582,10,4)
@test Dates.julian2datetime(2415020.5) == Dates.DateTime(1900,1,1)
@test Dates.julian2datetime(2415385.5) == Dates.DateTime(1901,1,1)
@test Dates.julian2datetime(2440587.5) == Dates.DateTime(1970,1,1)
@test Dates.julian2datetime(2444239.5) == Dates.DateTime(1980,1,1)
@test Dates.julian2datetime(2452695.625) == Dates.DateTime(2003,2,25,3)
@test Dates.datetime2julian(Dates.DateTime(2013,12,3,21)) == 2456630.375

@test typeof(Dates.now()) <: Dates.DateTime
@test typeof(Dates.today()) <: Dates.Date
@test typeof(Dates.now(Dates.UTC)) <: Dates.DateTime


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkWatcherListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{NetworkWatcher} } # spec name: value

    function NetworkWatcherListResult(;value=nothing)
        o = new()
        validate_property(NetworkWatcherListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type NetworkWatcherListResult

const _property_map_NetworkWatcherListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_NetworkWatcherListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{NetworkWatcher}")
Base.propertynames(::Type{ NetworkWatcherListResult }) = collect(keys(_property_map_NetworkWatcherListResult))
Swagger.property_type(::Type{ NetworkWatcherListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkWatcherListResult[name]))}
Swagger.field_name(::Type{ NetworkWatcherListResult }, property_name::Symbol) =  _property_map_NetworkWatcherListResult[property_name]

function check_required(o::NetworkWatcherListResult)
    true
end

function validate_property(::Type{ NetworkWatcherListResult }, name::Symbol, val)
end


text:
"""
    MITgcm_namelist(groups,params)

Data structure representing a MITgcm _namelist_ file, such as `data.pkg`, which contains

```
    groups :: Array{Symbol,1} = Array{Symbol,1}(undef, 0)
    params :: Array{OrderedDict{Symbol,Any},1} = Array{OrderedDict{Symbol,Any},1}(undef, 0)
```

with model parameters (`params`) being organized in `groups` as done in the files.

```
using MITgcmTools
fil=joinpath(MITgcm_path[1],"verification","advect_xy","run","data")
nml=read_namelist(fil)
MITgcm_namelist(nml.groups,nml.params)
MITgcm_namelist(groups=nml.groups,params=nml.params)
MITgcm_namelist(groups=nml.groups)
```
"""
Base.@kwdef struct MITgcm_namelist
    groups :: Array{Symbol,1} = Array{Symbol,1}(undef, 0)
    params :: Array{OrderedDict{Symbol,Any},1} = Array{OrderedDict{Symbol,Any},1}(undef, 0)
end

import Base:read,write
read(fil::AbstractString,nml::MITgcm_namelist) = read_namelist(fil)
write(fil::AbstractString,nml::MITgcm_namelist) = write_namelist(fil,nml)

"""
    MITgcm_config()

Concrete type of `AbstractModelConfig` for `MITgcm` (as part of the `ClimateModels.jl` interface for `MITgcm`)
which contains

```
    model :: String = "MITgcm"
    configuration :: String = ""
    options :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    outputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    channel :: Channel{Any} = Channel{Any}(10)
    folder :: String = tempdir()
    ID :: UUID = UUIDs.uuid4()
```

and can be constructed using keywords as follows

```
unknown_config=MITgcm_config(configuration="unknown")
```
"""
Base.@kwdef struct MITgcm_config <: AbstractModelConfig
    model :: String = "MITgcm"
    configuration :: String = ""
    options :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    outputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    channel :: Channel{Any} = Channel{Any}(10)
    folder :: String = tempdir()
    ID :: UUID = UUIDs.uuid4()
end


text:
immutable URange{T<:Real} <: AbstractUnitRange{T}
    start::T
    stop::T
    URange(start, stop) = new(start, urange_last(start,stop))
end
URange{T<:Real}(start::T, stop::T) = URange{T}(start, stop)

urange_last(::Bool, stop::Bool) = stop
urange_last{T<:Integer}(start::T, stop::T) =
    ifelse(stop >= start, stop, convert(T,start-one(stop-start)))
urange_last{T}(start::T, stop::T) =
    ifelse(stop >= start, convert(T,start+floor(stop-start)),
                          convert(T,start-one(stop-start)))

Base.start{T}(r::URange{T}) = oftype(r.start + one(T), r.start)

Base.intersect{T1<:Integer,T2<:Integer}(r::URange{T1}, s::URange{T2}) = URange(max(first(r),first(s)), min(last(r),last(s)))

@inline function Base.getindex{R,S<:Integer}(r::URange{R}, s::AbstractUnitRange{S})
    @boundscheck checkbounds(r, s)
    f = first(r)
    strt = f + first(s) - 1
    URange{R}(strt, strt+length(s)-1)
end

Base.promote_rule{T1,T2}(::Type{URange{T1}},::Type{URange{T2}}) =
    URange{promote_type(T1,T2)}
Base.convert{T<:Real}(::Type{URange{T}}, r::URange{T}) = r
Base.convert{T<:Real}(::Type{URange{T}}, r::URange) = URange{T}(r.start, r.stop)

Base.promote_rule{T1,UR<:AbstractUnitRange}(::Type{URange{T1}}, ::Type{UR}) =
    URange{promote_type(T1,eltype(UR))}
Base.promote_rule{T1,T2}(::Type{UnitRange{T2}}, ::Type{URange{T1}}) =
    URange{promote_type(T1,T2)}
Base.convert{T<:Real}(::Type{URange{T}}, r::AbstractUnitRange) = URange{T}(first(r), last(r))

let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    Base.start{T<:smallint}(r::URange{T}) = convert(Int, r.start)
end

Base.show(io::IO, r::URange) = print(io, typeof(r).name, '(', repr(first(r)), ',', repr(last(r)), ')')


text:
#######################################################################
# ParametricHamiltonian
#######################################################################
struct ParametricHamiltonian{P,N,M<:NTuple{N,ElementModifier},D<:NTuple{N,Any},H<:Hamiltonian}
    baseh::H
    h::H
    modifiers::M                   # N modifiers
    ptrdata::D                     # P is an NTuple{N,Vector{Vector{ptrdata}}}, one per harmonic
    allptrs::Vector{Vector{Int}}   # ptrdata may be a nzval ptr, a (ptr,r) or a (ptr, r, dr)
    parameters::NTuple{P,NameType} # allptrs are modified ptrs in each harmonic
end

function Base.show(io::IO, ::MIME"text/plain", pham::ParametricHamiltonian{N}) where {N}
    i = get(io, :indent, "")
    print(io, i, "Parametric")
    show(io, pham.h)
    print(io, i, "\n", "$i  Parameters       : $(parameters(pham))")
end

"""
    parametric(h::Hamiltonian, modifiers::ElementModifier...)

Builds a `ParametricHamiltonian` that can be used to efficiently apply `modifiers` to `h`.
`modifiers` can be any number of `@onsite!(args -> body; kw...)` and `@hopping!(args -> body;
kw...)` transformations, each with a set of parameters `ps` given as keyword arguments of
functions `f = (...; ps...) -> body`. The resulting `ph::ParamtricHamiltonian` can be used
to produced the modified Hamiltonian simply by calling it with those same parameters as
keyword arguments.

Note 1: for sparse `h`, `parametric` only modifies existing onsites and hoppings in `h`,
so be sure to add zero onsites and/or hoppings to `h` if they are originally not present but
you need to apply modifiers to them.

Note 2: `optimize!(h)` is called prior to building the parametric Hamiltonian. This can lead
to extra zero onsites and hoppings being stored in sparse `h`s.

    h |> parametric(modifiers::ElementModifier...)

Function form of `parametric`, equivalent to `parametric(h, modifiers...)`.

# Examples

```jldoctest
julia> ph = LatticePresets.honeycomb() |> hamiltonian(onsite(0) + hopping(1, range = 1/√3)) |>
       unitcell(10) |> parametric(@onsite!((o; μ) -> o - μ))
ParametricHamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 200 × 200
  Orbitals         : ((:a,), (:a,))
  Element type     : scalar (Complex{Float64})
  Onsites          : 0
  Hoppings         : 600
  Coordination     : 3.0
  Parameters       : (:μ,)

julia> ph(μ = 2)
Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 200 × 200
  Orbitals         : ((:a,), (:a,))
  Element type     : scalar (Complex{Float64})
  Onsites          : 200
  Hoppings         : 600
  Coordination     : 3.0
```
# See also
    `@onsite!`, `@hopping!`
"""
function parametric(h::Hamiltonian, ts::ElementModifier...)
    ts´ = resolve(ts, h.lattice)
    optimize!(h)  # to avoid ptrs getting out of sync if optimize! later
    allptrs = Vector{Int}[Int[] for _ in h.harmonics]
    ptrdata = parametric_ptrdata_tuple!(allptrs, h, ts´)
    foreach(sort!, allptrs)
    foreach(unique!, allptrs)
    params = parameters(ts...)
    return ParametricHamiltonian(h, copy(h), ts´, ptrdata, allptrs, params)
end

parametric(ts::ElementModifier...) = h -> parametric(h, ts...)

parametric_ptrdata_tuple!(allptrs, h, ts´::Tuple) = _parametric_ptrdata!(allptrs, h, ts´...)
_parametric_ptrdata!(allptrs, h, t, ts...) =
    (parametric_ptrdata!(allptrs, h, t), _parametric_ptrdata!(allptrs, h, ts...)...)
_parametric_ptrdata!(allptrs, h) = ()

function parametric_ptrdata!(allptrs, h::Hamiltonian{LA,L,M,<:AbstractSparseMatrix}, t::ElementModifier) where {LA,L,M}
    harmonic_ptrdata = empty_ptrdata(h, t)
    lat = h.lattice
    selector = t.selector
    for (har, ptrdata, allptrs_har) in zip(h.harmonics, harmonic_ptrdata, allptrs)
        matrix = har.h
        dn = har.dn
        rows = rowvals(matrix)
        for col in 1:size(matrix, 2), ptr in nzrange(matrix, col)
            row = rows[ptr]
            selected = selector(lat, (row, col), (dn, zero(dn)))
            if selected
                push!(ptrdata, ptrdatum(t, lat, ptr, (row, col), dn))
                push!(allptrs_har, ptr)
            end
        end
    end
    return harmonic_ptrdata
end

# Uniform case, one vector of nzval ptr per harmonic
empty_ptrdata(h, t::UniformModifier)  = Vector{Int}[Int[] for _ in h.harmonics]

# Non-uniform case, one vector of (ptr, r, dr) per harmonic
function empty_ptrdata(h, t::OnsiteModifier)
    S = positiontype(h.lattice)
    return Vector{Tuple{Int,S}}[Tuple{Int,S}[] for _ in h.harmonics]
end

function empty_ptrdata(h, t::HoppingModifier)
    S = positiontype(h.lattice)
    return Vector{Tuple{Int,S,S}}[Tuple{Int,S,S}[] for _ in h.harmonics]
end

# Uniform case
ptrdatum(t::UniformModifier, lat, ptr, (row, col), dn) = ptr

# Non-uniform case
ptrdatum(t::OnsiteModifier, lat, ptr, (row, col), dn) = (ptr, sites(lat)[col])

function ptrdatum(t::HoppingModifier, lat, ptr, (row, col), dn)
    r, dr = _rdr(sites(lat)[col], sites(lat)[row] + bravais(lat) * dn)
    return (ptr, r, dr)
end

function (ph::ParametricHamiltonian)(; kw...)
    checkconsistency(ph, false) # only weak check for performance
    reset_harmonic!.(ph.h.harmonics, ph.baseh.harmonics, ph.allptrs)
    applymodifier_ptrdata!.(Ref(ph.h), ph.modifiers, ph.ptrdata, Ref(values(kw)))
    return ph.h
end

function reset_harmonic!(har, basehar, ptrs)
    nz = nonzeros(har.h)
    onz = nonzeros(basehar.h)
    @simd for ptr in ptrs
        nz[ptr] = onz[ptr]
    end
    return har
end

function applymodifier_ptrdata!(h, modifier, ptrdata, kw)
    for (har, hardata) in zip(h.harmonics, ptrdata)
        nz = nonzeros(har.h)
        @simd for data in hardata  # @simd is valid because ptrs are not repeated
            ptr = first(data)
            args = modifier_args(nz, data)
            val = modifier(args...; kw...)
            nz[ptr] = val
        end
    end
    return h
end

# A negative ptr corresponds to a forced-hermitian element
modifier_args(nz, ptr::Int) = (nz[ptr],)
modifier_args(nz, (ptr, r)::Tuple{Int,SVector}) = (nz[ptr], r)
modifier_args(nz, (ptr, r, dr)::Tuple{Int,SVector,SVector}) = (nz[ptr], r, dr)

function checkconsistency(ph::ParametricHamiltonian, fullcheck = true)
    isconsistent = true
    length(ph.baseh.harmonics) == length(ph.h.harmonics) || (isconsitent = false)
    if fullcheck && isconsistent
        for (basehar, har) in zip(ph.baseh.harmonics, ph.h.harmonics)
            length(nonzeros(har.h)) == length(nonzeros(basehar.h)) || (isconsistent = false; break)
            rowvals(har.h) == rowvals(basehar.h) || (isconsistent = false; break)
            getcolptr(har.h) == getcolptr(basehar.h) || (isconsistent = false; break)
        end
    end
    isconsistent ||
        throw(error("ParametricHamiltonian is not internally consistent, it may have been modified after creation"))
    return nothing
end

# ParametricHamiltonian's are already optimized upon creation
optimize!(ph::ParametricHamiltonian) = ph

"""
    parameters(ph::ParametricHamiltonian)

Return the names of the parameter that `ph` depends on
"""
parameters(ph::ParametricHamiltonian) = ph.parameters

matrixtype(ph::ParametricHamiltonian) = matrixtype(parent(ph))

blockeltype(ph::ParametricHamiltonian) = blockeltype(parent(ph))

bravais(ph::ParametricHamiltonian) = bravais(ph.hamiltonian.lattice)

Base.parent(ph::ParametricHamiltonian) = ph.h

Base.copy(ph::ParametricHamiltonian) =
    ParametricHamiltonian(copy(ph.baseh), copy(ph.h), ph.modifiers, copy(h.ptrdata), copy(h.allptrs))

Base.size(ph::ParametricHamiltonian, n...) = size(ph.h, n...)

Base.eltype(ph::ParametricHamiltonian) = eltype(ph.h)

#######################################################################
# bloch! for parametric
#######################################################################

bloch(ph::ParametricHamiltonian, args...) = bloch!(similarmatrix(ph), ph, args...)

bloch!(matrix, ph::ParametricHamiltonian, pϕs = (), axis = 0) =
    bloch!(matrix, h_phases(ph, toSVector(pϕs))..., axis)

@inline function h_phases(ph::ParametricHamiltonian, pϕs)
    pnames = parameters(ph)
    ps, ϕs = extract_parameters_phases(pnames, pϕs)
    h = ph(; ps...)
    return (h, ϕs)
end

extract_parameters_phases(pnames::NTuple{N,NameType}, ϕs::SVector{M}) where {N,M} =
    (NamedTuple{pnames}(ntuple(i->ϕs[i], Val(N))), ntuple(i->ϕs[i+N], Val(M-N)))


text:
using Distributed

# Number of ITensorNetworkMaps in the sum
nmats = 2

addprocs(nmats)

# Set up for the calculation. Define everything
# in the `@everywhere` block that you need
# on all of the processes to perform the
# calculation.
@everywhere begin
  using ITensors
  using ITensors.ITensorNetworkMaps
  using LinearAlgebra
  using KrylovKit

  BLAS.set_num_threads(1)
  ITensors.Strided.disable_threads()
  ITensors.disable_threaded_blocksparse()

  include(joinpath("src", "DistributedSums.jl"))

  # Function that defines the tensor network
  # that will map the input vector (i.e.
  # a transfer matrix or effective Hamiltonian).
  function itensor_map(i, j)
    A1 = randomITensor(i', dag(i), dag(l))
    A2 = randomITensor(j', dag(j), l)
    return ITensorNetworkMap([A1, A2])
  end

  d = 2
  χ = 3
  i = Index(d, "i")
  j = Index(d, "j")
  l = Index(χ, "l")

  v0 = randomITensor(i, j)
end

# Defines an object representing a sum of
# linear maps (generalized matrices) built
# from contractions of ITensors that are
# distributed across `nmats` workers
# (each linear map lives on a worker).
A = DistributedSum(_ -> itensor_map(i, j), nmats)

# Apply the linear map, which internally
# is performed in parallel where the linear
# map is created.
Av0 = A(v0)

@show norm(Av0 - sum(fetch.(A))(v0))

# Solve an eigenvector equation in
# parallel.
λ⃗, v⃗ = eigsolve(A, v0, 1, :LM)
λ = λ⃗[1]
v = v⃗[1]

@show norm(A(v) - λ * v)

max_stars_repo_name: mtfishman/ITensorParallel


text:
# create_sysimage(
#            [:CSV, :DataFrames, :FluxDash, :Flux, :Random],
#                sysimage_path="sysnew_fluxdash.so",
#                    precompile_execution_file="FluxDash\\test\\nn_train.jl"
#                    )
run(`julia --sysimage sysnew_fluxdash.so --project run.jl`)

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
module Emporos

using MacroTools: MacroTools

export IBKR
export @wrapaccessors, gen_cpp_wrapaccessors

const PACKAGE_ROOT = @__DIR__

include("macros.jl")
include("util.jl")

include("./IBKR/IBKR.jl")
using .IBKR

end


text:
using Swagger
using Test

# set the current julia executable to be used henceforth
("JULIA" in keys(ENV)) || (ENV["JULIA"] = joinpath(Sys.BINDIR, "julia"))

const gencmd = joinpath(dirname(@__FILE__()), "petstore", "generate.sh")
println("Generating petstore using $gencmd")
run(`$gencmd`)

println("Running petstore tests...")
include("petstore/runtests.jl")


text:
using PSCOPF

using Test
using Dates

@testset verbose=true "test_compute_firmness" begin

    #=
       10h30        12h           14h
    TS-Delta        ECH            TS
           |          |             |
            <---------------------->
                    Delta
                    3h30
    =#

    @testset "ech_late" begin
        ech = Dates.DateTime("2015-01-01T12:00:00")
        ts = Dates.DateTime("2015-01-01T14:00:00")
        delta = Dates.Minute(210) #3h30

        next_ech = Dates.DateTime("2015-01-01T05:00:00")
        @test_throws ErrorException PSCOPF.compute_firmness(ech, next_ech, ts, delta)

        next_ech = Dates.DateTime("2015-01-01T13:00:00")
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.DECIDED

        next_ech = nothing
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.DECIDED
    end

    #=
       10h30                      14h
    TS-Delta                       TS
         ECH
           |                        |
            <---------------------->
                    Delta
                    3h30
    =#
    @testset "ech_last_moment" begin
        ech = Dates.DateTime("2015-01-01T10:30:00")
        ts = Dates.DateTime("2015-01-01T14:00:00")
        delta = Dates.Minute(210) #3h30

        next_ech = Dates.DateTime("2015-01-01T05:00:00")
        @test_throws ErrorException PSCOPF.compute_firmness(ech, next_ech, ts, delta)

        next_ech = Dates.DateTime("2015-01-01T10:30:00")
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.FREE

        next_ech = Dates.DateTime("2015-01-01T13:00:00")
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.TO_DECIDE

        next_ech = nothing
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.TO_DECIDE
    end

    #=
      7h             10h30                      14h
     ECH          TS-Delta                       TS
       |                 |                        |
                          <---------------------->
                                    Delta
                                     3h30
    =#
    @testset "ech_early" begin
        ech = Dates.DateTime("2015-01-01T07:00:00")
        ts = Dates.DateTime("2015-01-01T14:00:00")
        delta = Dates.Minute(210) #3h30

        next_ech = Dates.DateTime("2015-01-01T05:00:00")
        @test_throws ErrorException PSCOPF.compute_firmness(ech, next_ech, ts, delta)

        next_ech = Dates.DateTime("2015-01-01T10:30:00")
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.FREE

        next_ech = Dates.DateTime("2015-01-01T13:00:00")
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.TO_DECIDE

        next_ech = nothing
        @test PSCOPF.compute_firmness(ech, next_ech, ts, delta) == PSCOPF.TO_DECIDE
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
`gen`

look for a list of genotypes from the MouseID
"""
function gen(str; dir = joinpath(dirname(@__DIR__), "genotypes"))
    genotype = "missing"
    for file in readdir(dir)
        if endswith(file, ".csv")
            df = FileIO.load(joinpath(dir, file)) |> DataFrame
            n = names(df)[1]
            if str in df[!,n]
                genotype = string(n)
            end
        end
    end
    genotype
end

"""
`pharm`
"""
function pharm(str; dir = joinpath(dirname(@__DIR__), "pharmacology"))
    drug = "missing"
    for file in readdir(dir)
        if endswith(file, ".csv")
            df = FileIO.load(joinpath(dir, file)) |> DataFrame
            n = names(df)[1]
            if str in df[!,n]
                drug = String(n)
            end
        end
    end
    return drug
end


text:
@symbol_func function K_F(cur_reactor::AbstractReactor)
  cur_K = 278.3

  cur_K *= cur_reactor.f_D ^ 2

  cur_K *= cur_reactor.epsilon ^ 2

  cur_K *= cur_reactor.kappa_95

  cur_K *= g(cur_reactor)

  cur_K
end


text:
using FunctionObserving
using Test

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


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %%
?propertynames

# %%
module O

Base.@kwdef mutable struct Foo
    pub::String = "meow"
    priv::Int = 3
end

function Base.getproperty(foo::Foo, name::Symbol)
    name === :multiple && return getfield(foo, :pub)^getfield(foo, :priv)
    hasproperty(foo, name) && return getfield(foo, name)
    error("type Foo has no public property $name")
end

Base.show(io::IO, foo::Foo) = print(io, "Foo(pub = ", repr(foo.pub), ')')

Base.propertynames(foo::Foo, private::Bool=false) =
    private ? (fieldnames(typeof(foo))..., :multiple) : (:pub, :multiple)

end

# %%
foo = O.Foo()

# %%
propertynames(foo)

# %%
@which propertynames(foo)

# %%
methodswith(O.Foo)

# %%
propertynames(foo, true)

# %%
foo.multiple

# %%
foo.pub = "bow"
foo

# %%
foo.multiple

# %%
foo.priv

# %%
foo.pub

# %%
f(foo::O.Foo) = foo.multiple

# %%
f(foo)

# %%
@code_typed f(foo)

# %%
g(foo::O.Foo) = foo.pub

# %%
g(foo)

# %%
@code_typed g(foo)

# %%
foo.priv = 10

# %%
foo.multiple

# %%

max_stars_repo_path: 0005/propertynames example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/propertynames example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/propertynames example.jl

max_forks_repo_name: genkuroki/public


text:
struct TimeFunctionContainer
    times::Vector{Function}
    dtimes::Vector{Function}
    ddtimes::Vector{Function}

    function TimeFunctionContainer(network, f)
        times = Vector{Function}()
        dtimes = Vector{Function}()
        ddtimes = Vector{Function}()
        for l in links(network)
            a, b, c = f(props(l))
            push!(times, a)
            push!(dtimes, b)
            push!(ddtimes, c)
        end

        new(times, dtimes, ddtimes)
    end
end

function times(tfc::TimeFunctionContainer, x::Array{<:Real,1}, ids::Union{Nothing,Array{<:Integer,1}}=nothing)
    t = try
        funcmap(tfc.times[ids], x)
    catch
        funcmap(tfc.times[ids], x[ids])
    end
    t
end

function dtimes(tfc::TimeFunctionContainer, x::Array{<:Real,1}, ids::Union{Nothing,Array{<:Integer,1}}=nothing)
    dt = try
        funcmap(tfc.dtimes[ids], x)
    catch
        funcmap(tfc.dtimes[ids], x[ids])
    end
    dt
end

function ddtimes(tfc::TimeFunctionContainer, x::Array{<:Real,1}, ids::Union{Nothing,Array{<:Integer,1}}=nothing)
    ddt = try
        funcmap(tfc.ddtimes[ids], x)
    catch DimensionMismatch
        funcmap(tfc.ddtimes[ids], x[ids])
    end
    ddt
end


text:
using CubicEoSDatabase

"""
    load(MBWREoSComponent; name::AbstractString[, component_dbs])

Return `BrusilovskyEoSComponent` named `name` with parameters from `component_dbs`.

# Arguments
- the first argument is `BrusilovskyEoSComponent`

# Keywords
- `name::AbstractString`: the component name
- `component_dbs`: an iterable of `ComponentDatabase` objects (by default,
(Data.martinez(), Data.brusilovsky_comp()))

# Returns
- `BrusilovskyEoSComponent`: the component generated by merging the EoS properties for the
given name from databases
"""
function load(
    ::Type{<:MBWREoSComponent};
    name::AbstractString,
    component_dbs = (ComponentDatabase("data/nist.csv"),)
)
    comp_properties = foldl(
        (dict, db) -> merge!(dict, getentry(db, name)),
        component_dbs;
        init = Dict{Symbol, Any}()
    )

    return __load_mbwr_comp__(; comp_properties...)
end

function __load_mbwr_comp__(
    ;
    name::AbstractString,
    molecular_mass,
    number_carbons,
    critical_temperature,
    critical_pressure,
    critical_compressibility,
    acentric_factor,
    extra_kw...  # unneeded keywords to construct component
)
    __load_mbwr_comp__(
        name,
        molecular_mass,
        number_carbons,
        critical_temperature,
        critical_pressure,
        critical_compressibility,
        acentric_factor,
    )
end

function __load_mbwr_comp__(
    name::AbstractString,
    molecular_mass,
    number_carbons,
    critical_temperature,
    critical_pressure,
    critical_compressibility,
    acentric,
)

    return MBWREoSComponent(
        ;
        name = name,
        critical_pressure = critical_pressure,
        critical_temperature = critical_temperature,
        acentric_factor = acentric,
        Zc = critical_compressibility,
        molar_mass = molecular_mass,
        carbon_number = number_carbons
    )
end

"""
    load(BrusilovskyEoSMixture; names[, component_dbs])

Return `BrusilovskyEoSMixture` by loading parameters for individual substances named `names`
from `component_dbs` and their interaction parameters from `mix_eos_db`.

# Arguments
- the first argument is `BrusilovskyEoSMixture`

# Keywords
- `names`: an iterable of `AbstractString`s containing the names of components to load
- `component_dbs`: an iterable of `ComponentDatabase` objects (by default,
`(Data.martinez(), Data.brusilovsky_comp())`)
- `mix_eos_db::MixtureDatabase`: the database with the cross-component interaction
parameters (by default, `Data.brusilovsky_mix()`)

# Returns
- `BrusilovskyEoSMixture`: the EoS for the mixture of the specified substances
"""
function load(
    ::Type{<:MBWREoSMixture};
    names,
    component_dbs = (Data.martinez(), Data.brusilovsky_comp()),
    mix_eos_db::MixtureDatabase = Data.brusilovsky_mix_adjusted(),
)
    components = [
        load(
            BrusilovskyEoSComponent;
            name = name,
            component_dbs = component_dbs
        ) for name in names
    ]
    corrections = getmatrix(mix_eos_db, names)
    return MBWREoSMixture(;
        components = components,
        corrections...
    )
end


text:
##
using FixedPointNumbers
using FixedPointDecimals
##
x = FixedDecimal{Int,2}(5)
## Emulate the Mario Kart system to rotate a dinosaur
MKFP = Fixed{Int16,8} # use 16 bits for a signed integer and 8 bits for fraction
val = MKFP(5.25)

## Get some data off the internet and plot it
data = [24,18, 24,17, 22,16, 20,11, 19,6, 19,2, 17,2, 17,6,
16,5, 15,2, 13,2, 14,5, 14,6, 12,6, 12,2, 10,2, 10,4,
9,2, 7,2, 9,6, 7,6, 4,4, 2,3, 0,2, 1,3, 3,5, 5,9,
9,11, 17,11, 21,17, 23,18, 24,18]
P = reshape(data,2,div(length(data),2))
x = P[1,:]
y = P[2,:]
using Plots
scatter(x,y,legend=false)

## Plot it in Mario Kart numbers and connect the dots
x = MKFP.(x)
y = MKFP.(y)
plot!(x,y,legend=false) # no change!
##
deg = 63
theta = deg/180*pi
R(theta) = [cos(theta) -sin(theta); sin(theta) cos(theta)]
Pp = R(theta)*P # True
plot(Pp[1,:],Pp[2,:],linestyle = :solid,marker = :circle,color=2,lab="Float64")
# Rotate in fixed Point
Pmk = MKFP.(R(MKFP(theta)))*MKFP.(P)
plot!(Pmk[1,:],Pmk[2,:],linestyle = :solid,marker = :circle,color=4,lab="Mario Kart")
##
theme(:dark)
plot!(size=(350,350))
savefig("dinofig-rotated.pdf")
##

## Save original fig
theme(:dark)
plot(x,y,legend=false,marker=:circle,color=2,size=(350,350)) # no change!
savefig("dinofig.pdf")


text:
using Base.Test
using DataFrames

let
	unsorted_dv = DataVector[2, 1, NA]
	sorted_dv = DataVector[NA, 1, 2]

	@assert isequal(sort(unsorted_dv), sorted_dv)
	@assert isequal(sortperm(unsorted_dv), [3, 2, 1])
	# TODO: Make this work
	# tiedrank(dv)

	@assert first(unsorted_dv) == 2
	@assert isna(last(unsorted_dv))
end


text:
"""
     Credential(; api_key)

A credential object for Airtable.
If the api_key or api_token are not provided,
they will be read from the `AIRTABLE_KEY` environment variable.
Go to [Airtable account settings](https://airtable.com/account)
to aquire your credentials.

```jldoctest; setup = :(using Airtable)
# after running `export AIRTABLE_KEY=<api key>` in the shell
julia> key = Airtable.Credential()
Airtable.Credential(<secrets>)
```
"""
struct Credential
    api_key::String
end

function Credential(; api_key=Base.get(ENV, "AIRTABLE_KEY", nothing))
    isnothing(api_key) && throw(ArgumentError("Environment does not have `$AIRTABLE_KEY` set. Must past api key directly"))
    return Credential(api_key)
end

Base.show(io::IO, ::Credential) = println(io, "Airtable.Credential(<secrets>)")

@testset "Credentials" begin
    key = Airtable.Credential()
    @test key isa Airtable.Credential
    @test key.api_key isa String
end


text:
export getEdgeMassMatrix, getdEdgeMassMatrix

# Edge mass matrix for isotropic coefficient
function getEdgeMassMatrix(mesh::TetraMesh, Sigma::Array{Float64,1})

	nc = mesh.nc
  ne = mesh.ne
  nn = mesh.nn

	M = zeros(6,6,nc)
  N = round(Int128,nn)
  globalEdgeNumbering = arrayToUniqueInt(mesh.edges,N)
  localEdgeNumbering  = zeros(Int64,nc,6)

	for i=1:nc

		Ti = mesh.Tetras[i,:]
    Xi = mesh.Points[Ti,:]

		localEdges = [
      Ti[1] Ti[2]
      Ti[1] Ti[3]
      Ti[1] Ti[4]
      Ti[2] Ti[3]
      Ti[2] Ti[4]
      Ti[3] Ti[4]]
		localEdgeNumbering[i,:] = arrayToUniqueInt(localEdges,N)

		M[:,:,i] = Sigma[i] * getLocalEdgeMassMatrix(Xi)

	end

	edgeInd = reshape(indexin(localEdgeNumbering,globalEdgeNumbering),nc,6)

	ii = zeros(Int64,36*nc)
	jj = zeros(Int64,36*nc)

	cnt = 1
	for i=1:nc
  	for i2=1:6
	  	for i1=1:6
				ii[cnt] = edgeInd[i,i1]
				jj[cnt] = edgeInd[i,i2]
				cnt    += 1
			end
		end
	end
  mm = vec(M)

	return sparse(ii,jj,mm,ne,ne)

end

# Edge mass matrix for anisotropic coefficient (symmetric tensor)
function getEdgeMassMatrix(mesh::TetraMesh, Sigma::Array{Float64,2})

	nc = mesh.nc
  ne = mesh.ne
  nn = mesh.nn

	M = zeros(6,6,nc)
  N = round(Int128,nn)
  globalEdgeNumbering = arrayToUniqueInt(mesh.edges,N)
  localEdgeNumbering  = zeros(Int64,nc,6)

	for i=1:nc

		Ti = mesh.Tetras[i,:]
    Xi = mesh.Points[Ti,:]

		localEdges = [
      Ti[1] Ti[2]
      Ti[1] Ti[3]
      Ti[1] Ti[4]
      Ti[2] Ti[3]
      Ti[2] Ti[4]
      Ti[3] Ti[4]]
		localEdgeNumbering[i,:] = arrayToUniqueInt(localEdges,N)

		M[:,:,i] = getLocalEdgeMassMatrix(Xi,Sigma[i,:])

	end

	edgeInd = reshape(indexin(localEdgeNumbering,globalEdgeNumbering),nc,6)

	ii = zeros(Int64,36*nc)
	jj = zeros(Int64,36*nc)

	cnt = 1
	for i=1:nc
  	for i2=1:6
	  	for i1=1:6
				ii[cnt] = edgeInd[i,i1]
				jj[cnt] = edgeInd[i,i2]
				cnt    += 1
			end
		end
	end
  mm = vec(M)

	return sparse(ii,jj,mm,ne,ne)

end

# Derivative of edge mass matrix times vector product with
# isotropic coefficient
function getdEdgeMassMatrix(mesh::TetraMesh, Sigma::Array{Float64,1}, v::Array{Float64,1})

  nc = mesh.nc
  ne = mesh.ne
  nn = mesh.nn

  N                   = round(Int128,nn)
  globalEdgeNumbering = arrayToUniqueInt(mesh.edges,N)
  localEdgeNumbering  = zeros(Int64,nc,6)

  for i=1:nc

		Ti = mesh.Tetras[i,:]
		localEdges = [
      Ti[1] Ti[2]
      Ti[1] Ti[3]
      Ti[1] Ti[4]
      Ti[2] Ti[3]
      Ti[2] Ti[4]
      Ti[3] Ti[4]]
		localEdgeNumbering[i,:] = arrayToUniqueInt(localEdges,N)

  end

  edgeInd  = reshape(indexin(localEdgeNumbering,globalEdgeNumbering),nc,6)
  sigmaInd = collect(1:nc)

  dMv = zeros(6,nc)

	for i=1:nc

		Ti = mesh.Tetras[i,:]
    Xi = mesh.Points[Ti,:]

    vi = v[edgeInd[i,:]]
    Mi = getLocalEdgeMassMatrix(Xi)
		dMv[:,i] = Mi * vi

	end

	ii = zeros(Int64,6*nc)
	jj = zeros(Int64,6*nc)

	cnt = 1
	for i=1:nc
  	for i1=1:6 # edges
			ii[cnt] = edgeInd[i,i1]
			jj[cnt] = sigmaInd[i]
			cnt    += 1
		end
	end
  mm = vec(dMv)

	return sparse(ii,jj,mm,ne,nc)

end

# Derivative of edge mass matrix times vector product with
# anisotropic coefficient (symmetric tensor)
function getdEdgeMassMatrix(mesh::TetraMesh, Sigma::Array{Float64,2}, v::Array{Float64,1})

  nc = mesh.nc
  ne = mesh.ne
  nn = mesh.nn

  N                   = round(Int128,nn)
  globalEdgeNumbering = arrayToUniqueInt(mesh.edges,N)
  localEdgeNumbering  = zeros(Int64,nc,6)

  for i=1:nc

		Ti = mesh.Tetras[i,:]
		localEdges = [
      Ti[1] Ti[2]
      Ti[1] Ti[3]
      Ti[1] Ti[4]
      Ti[2] Ti[3]
      Ti[2] Ti[4]
      Ti[3] Ti[4]]
		localEdgeNumbering[i,:] = arrayToUniqueInt(localEdges,N)

  end

  edgeInd  = reshape(indexin(localEdgeNumbering,globalEdgeNumbering),nc,6)
  sigmaInd = reshape(collect(1:6*nc),nc,6)

  dMv = zeros(6,6,nc)

	for i=1:nc

		Ti = mesh.Tetras[i,:]
    Xi = mesh.Points[Ti,:]

    vi = v[edgeInd[i,:]]
    dMv[:,:,i] = getdLocalEdgeMassMatrix(Xi,vi)

	end

	ii = zeros(Int64,36*nc)
	jj = zeros(Int64,36*nc)

  cnt = 1
  for i=1:nc
    for i2=1:6 # sigma components
      for i1=1:6 # edges
        ii[cnt] = edgeInd[i,i1]
        jj[cnt] = sigmaInd[i,i2]
        cnt    += 1
      end
    end
  end
  mm = vec(dMv)

	return sparse(ii,jj,mm,ne,6*nc)

end

# Element edge mass matrix with unit coefficient
function getLocalEdgeMassMatrix(X::Array{Float64})
#              * X4
#             /|\
#            / | \
#           /  |  \
#          /   |   \
#         /    |    \
#     t6 /     |     \t5
#       /      |      \
#      /       |t3     \
#     /        |        \
#    /         |         \
# X3/_________ |t4_______ \ X2
#    .         |         .
#      .       |       .
#        .     |     . t1
#     t2   .   |   .
#            . * .
#             X1

	psi = [ zeros(3) for i = 1:4, j = 1:6 ]
	s = (1 - 1 / sqrt(5)) / 4
	t = 1 - 3 * s

	b1 = X[1,:]
	b2 = X[2,:]
	b3 = X[3,:]
	b4 = X[4,:]
	B  = [b2-b1 b3-b1 b4-b1]
	Binv = inv3X3(B)

	gradLambda1 = Binv' * [-1,-1,-1]
	gradLambda2 = Binv' * [ 1, 0, 0]
	gradLambda3 = Binv' * [ 0, 1, 0]
	gradLambda4 = Binv' * [ 0, 0, 1]

  L        = norm(b2 - b1)
	psi[1,1] = L * (t * gradLambda2 - s * gradLambda1) #psi_12 = |b2-b1|*psi_12   basis functions
	psi[2,1] = L * (s * gradLambda2 - t * gradLambda1)
	psi[3,1] = L * (s * gradLambda2 - s * gradLambda1)
	psi[4,1] = L * (s * gradLambda2 - s * gradLambda1)

  L        = norm(b3 - b1)
	psi[1,2] = L * (t * gradLambda3 - s * gradLambda1) #psi_13 = |b3-b1|*psi_13   basis functions
	psi[2,2] = L * (s * gradLambda3 - s * gradLambda1)
	psi[3,2] = L * (s * gradLambda3 - t * gradLambda1)
	psi[4,2] = L * (s * gradLambda3 - s * gradLambda1)

  L        = norm(b4 - b1)
	psi[1,3] = L * (t * gradLambda4 - s * gradLambda1) #psi_14 = |b4-b1|*psi_14   basis functions
	psi[2,3] = L * (s * gradLambda4 - s * gradLambda1)
	psi[3,3] = L * (s * gradLambda4 - s * gradLambda1)
	psi[4,3] = L * (s * gradLambda4 - t * gradLambda1)

  L        = norm(b3 - b2)
	psi[1,4] = L * (s * gradLambda3 - s * gradLambda2) #psi_23 = |b3-b2|*psi_23   basis functions
	psi[2,4] = L * (t * gradLambda3 - s * gradLambda2)
	psi[3,4] = L * (s * gradLambda3 - t * gradLambda2)
	psi[4,4] = L * (s * gradLambda3 - s * gradLambda2)

  L        = norm(b4 - b2)
	psi[1,5] = L * (s * gradLambda4 - s * gradLambda2) #psi_24 = |b4-b2|*psi_24   basis functions
	psi[2,5] = L * (t * gradLambda4 - s * gradLambda2)
	psi[3,5] = L * (s * gradLambda4 - s * gradLambda2)
	psi[4,5] = L * (s * gradLambda4 - t * gradLambda2)

  L        = norm(b4 - b3)
	psi[1,6] = L * (s * gradLambda4 - s * gradLambda3) #psi_34 = |b4-b3|*psi_34   basis functions
	psi[2,6] = L * (s * gradLambda4 - s * gradLambda3)
	psi[3,6] = L * (t * gradLambda4 - s * gradLambda3)
	psi[4,6] = L * (s * gradLambda4 - t * gradLambda3)

	M = zeros(6,6)
  for q=1:4
    for k=1:6
      vk = psi[q,k]
      for j=1:6
        M[j,k] += dot(psi[q,j], vk)
      end
    end
  end

	V = getCellVolume(X[1,:], X[2,:], X[3,:], X[4,:])
	return (V/4)*M

end

# Element edge mass matrix with symmetric tensor coefficient
function getLocalEdgeMassMatrix(X::Array{Float64}, sigma::Array{Float64})
#              * X4
#             /|\
#            / | \
#           /  |  \
#          /   |   \
#         /    |    \
#     t6 /     |     \t5
#       /      |      \
#      /       |t3     \
#     /        |        \
#    /         |         \
# X3/_________ |t4_______ \ X2
#    .         |         .
#      .       |       .
#        .     |     . t1
#     t2   .   |   .
#            . * .
#             X1

	psi = [ zeros(3) for i = 1:4, j = 1:6 ]
	s = (1 - 1 / sqrt(5)) / 4
	t = 1 - 3 * s

	b1 = X[1,:]
	b2 = X[2,:]
	b3 = X[3,:]
	b4 = X[4,:]
	B  = [b2-b1 b3-b1 b4-b1]
	Binv = inv3X3(B)

	gradLambda1 = Binv' * [-1,-1,-1]
	gradLambda2 = Binv' * [ 1, 0, 0]
	gradLambda3 = Binv' * [ 0, 1, 0]
	gradLambda4 = Binv' * [ 0, 0, 1]

  L        = norm(b2 - b1)
	psi[1,1] = L * (t * gradLambda2 - s * gradLambda1) #psi_12 = |b2-b1|*psi_12   basis functions
	psi[2,1] = L * (s * gradLambda2 - t * gradLambda1)
	psi[3,1] = L * (s * gradLambda2 - s * gradLambda1)
	psi[4,1] = L * (s * gradLambda2 - s * gradLambda1)

  L        = norm(b3 - b1)
	psi[1,2] = L * (t * gradLambda3 - s * gradLambda1) #psi_13 = |b3-b1|*psi_13   basis functions
	psi[2,2] = L * (s * gradLambda3 - s * gradLambda1)
	psi[3,2] = L * (s * gradLambda3 - t * gradLambda1)
	psi[4,2] = L * (s * gradLambda3 - s * gradLambda1)

  L        = norm(b4 - b1)
	psi[1,3] = L * (t * gradLambda4 - s * gradLambda1) #psi_14 = |b4-b1|*psi_14   basis functions
	psi[2,3] = L * (s * gradLambda4 - s * gradLambda1)
	psi[3,3] = L * (s * gradLambda4 - s * gradLambda1)
	psi[4,3] = L * (s * gradLambda4 - t * gradLambda1)

  L        = norm(b3 - b2)
	psi[1,4] = L * (s * gradLambda3 - s * gradLambda2) #psi_23 = |b3-b2|*psi_23   basis functions
	psi[2,4] = L * (t * gradLambda3 - s * gradLambda2)
	psi[3,4] = L * (s * gradLambda3 - t * gradLambda2)
	psi[4,4] = L * (s * gradLambda3 - s * gradLambda2)

  L        = norm(b4 - b2)
	psi[1,5] = L * (s * gradLambda4 - s * gradLambda2) #psi_24 = |b4-b2|*psi_24   basis functions
	psi[2,5] = L * (t * gradLambda4 - s * gradLambda2)
	psi[3,5] = L * (s * gradLambda4 - s * gradLambda2)
	psi[4,5] = L * (s * gradLambda4 - t * gradLambda2)

  L        = norm(b4 - b3)
	psi[1,6] = L * (s * gradLambda4 - s * gradLambda3) #psi_34 = |b4-b3|*psi_34   basis functions
	psi[2,6] = L * (s * gradLambda4 - s * gradLambda3)
	psi[3,6] = L * (t * gradLambda4 - s * gradLambda3)
	psi[4,6] = L * (s * gradLambda4 - t * gradLambda3)

	Sigma = [
    sigma[1] sigma[4] sigma[5]
    sigma[4] sigma[2] sigma[6]
    sigma[5] sigma[6] sigma[3] ]

	M = zeros(6,6)
  for q=1:4
    for k=1:6
      vk = Sigma * psi[q,k]
      for j=1:6
        M[j,k] += dot(psi[q,j], vk)
      end
    end
  end

	V = getCellVolume(X[1,:], X[2,:], X[3,:], X[4,:])
	return (V/4)*M

end

# Derivative of element edge mass matrix - vector product
# with symmetric tensor coefficient
function getdLocalEdgeMassMatrix(X::Array{Float64}, v::Array{Float64,1})
#              * X4
#             /|\
#            / | \
#           /  |  \
#          /   |   \
#         /    |    \
#     t6 /     |     \t5
#       /      |      \
#      /       |t3     \
#     /        |        \
#    /         |         \
# X3/_________ |t4_______ \ X2
#    .         |         .
#      .       |       .
#        .     |     . t1
#     t2   .   |   .
#            . * .
#             X1

	psi = [ zeros(3) for i = 1:4, j = 1:6 ]
	s = (1 - 1 / sqrt(5)) / 4
	t = 1 - 3 * s

	b1 = X[1,:]
	b2 = X[2,:]
	b3 = X[3,:]
	b4 = X[4,:]
	B  = [b2-b1 b3-b1 b4-b1]
	Binv = inv3X3(B)

	gradLambda1 = Binv' * [-1,-1,-1]
	gradLambda2 = Binv' * [ 1, 0, 0]
	gradLambda3 = Binv' * [ 0, 1, 0]
	gradLambda4 = Binv' * [ 0, 0, 1]

  L        = norm(b2 - b1)
	psi[1,1] = L * (t * gradLambda2 - s * gradLambda1) #psi_12 = |b2-b1|*psi_12   basis functions
	psi[2,1] = L * (s * gradLambda2 - t * gradLambda1)
	psi[3,1] = L * (s * gradLambda2 - s * gradLambda1)
	psi[4,1] = L * (s * gradLambda2 - s * gradLambda1)

  L        = norm(b3 - b1)
	psi[1,2] = L * (t * gradLambda3 - s * gradLambda1) #psi_13 = |b3-b1|*psi_13   basis functions
	psi[2,2] = L * (s * gradLambda3 - s * gradLambda1)
	psi[3,2] = L * (s * gradLambda3 - t * gradLambda1)
	psi[4,2] = L * (s * gradLambda3 - s * gradLambda1)

  L        = norm(b4 - b1)
	psi[1,3] = L * (t * gradLambda4 - s * gradLambda1) #psi_14 = |b4-b1|*psi_14   basis functions
	psi[2,3] = L * (s * gradLambda4 - s * gradLambda1)
	psi[3,3] = L * (s * gradLambda4 - s * gradLambda1)
	psi[4,3] = L * (s * gradLambda4 - t * gradLambda1)

  L        = norm(b3 - b2)
	psi[1,4] = L * (s * gradLambda3 - s * gradLambda2) #psi_23 = |b3-b2|*psi_23   basis functions
	psi[2,4] = L * (t * gradLambda3 - s * gradLambda2)
	psi[3,4] = L * (s * gradLambda3 - t * gradLambda2)
	psi[4,4] = L * (s * gradLambda3 - s * gradLambda2)

  L        = norm(b4 - b2)
	psi[1,5] = L * (s * gradLambda4 - s * gradLambda2) #psi_24 = |b4-b2|*psi_24   basis functions
	psi[2,5] = L * (t * gradLambda4 - s * gradLambda2)
	psi[3,5] = L * (s * gradLambda4 - s * gradLambda2)
	psi[4,5] = L * (s * gradLambda4 - t * gradLambda2)

  L        = norm(b4 - b3)
	psi[1,6] = L * (s * gradLambda4 - s * gradLambda3) #psi_34 = |b4-b3|*psi_34   basis functions
	psi[2,6] = L * (s * gradLambda4 - s * gradLambda3)
	psi[3,6] = L * (t * gradLambda4 - s * gradLambda3)
	psi[4,6] = L * (s * gradLambda4 - t * gradLambda3)

  dSigma = [ zeros(3,3) for m = 1:6 ]
  dSigma[1][1,1] = 1.0
  dSigma[2][2,2] = 1.0
  dSigma[3][3,3] = 1.0
  dSigma[4][1,2] = dSigma[4][2,1] = 1.0
  dSigma[5][1,3] = dSigma[5][3,1] = 1.0
  dSigma[6][2,3] = dSigma[6][3,2] = 1.0

	dMv = zeros(6,6)
  for q = 1:4
    for m = 1:6
      for k = 1:6
        u = (dSigma[m] * psi[q,k]) * v[k]
        for j = 1:6
          dMv[j,m] += dot(psi[q,j], u)
        end
      end
    end
  end

	V = getCellVolume(X[1,:], X[2,:], X[3,:], X[4,:])
	return (V/4)*dMv

end

function inv3X3(A)
  Ainv = 1/det(A) * [
	  A[2,2]*A[3,3]-A[2,3]*A[3,2]  A[1,3]*A[3,2]-A[1,2]*A[3,3]  A[1,2]*A[2,3]-A[1,3]*A[2,2]
    A[2,3]*A[3,1]-A[2,1]*A[3,3]  A[1,1]*A[3,3]-A[1,3]*A[3,1]  A[1,3]*A[2,1]-A[1,1]*A[2,3]
	  A[2,1]*A[3,2]-A[2,2]*A[3,1]  A[1,2]*A[3,1]-A[1,1]*A[3,2]  A[1,1]*A[2,2]-A[1,2]*A[2,1] ]
	return Ainv
end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
module DiscretisedFluidQueues

import Jacobi, LinearAlgebra, SparseArrays, StaticArrays
import Base: *, +, size, show, getindex, setindex!, sum

# Types
export Model, DiscretisedFluidQueue, BoundedFluidQueue, Phase, PhaseSet # Queues are <:Model
export AbstractMatrixExponential, ConcentratedMatrixExponential, MatrixExponential
export Mesh, DGMesh, FRAPMesh, FVMesh # are <:Mesh
export Generator, FullGenerator, LazyGenerator # are <:Generator
export SFMDistribution, AbstractIntegrationMethod, AutoQuadrature, Quadrature, TrapezoidRule
export ExplicitRungeKuttaScheme, ForwardEuler, Heuns, StableRK3, StableRK4 # Euler, RungeKutta4 <: TimeIntegrationScheme
export Simulation

# Functions
export augment_model, membership, N₋, N₊, n_phases, phases, rates # Model methods
export cell_nodes, Δ, n_bases_per_cell, n_bases_per_phase, n_intervals, total_n_bases # Mesh methods
export interior_point_mass, left_point_mass, right_point_mass, integrate_time # SFMDistribution methods
export simulate, fixed_time, n_jumps, first_exit_x # Simulation methods
export build_lazy_generator, build_full_generator, static_generator
export cme_params, pdf, ccdf, cdf, build_me, cell_probs
export normalised_closing_operator_cdf, normalised_closing_operator_pdf
export naive_normalised_closing_operator_cdf, naive_normalised_closing_operator_pdf
export unnormalised_closing_operator_cdf, unnormalised_closing_operator_pdf
export limit, Limiter, NoLimiter, GeneralisedMUSCL

# model
include("SFM.jl")

include("abstract_mesh.jl") # things which apply to all meshes
include("discretised_fluid_queue.jl")
include("lazy_generators.jl")
include("full_generators.jl")

# auxillary functions
include("ME_tools.jl") # used in FRAPApproximation.jl
include("polynomials.jl") # used in discontinuous_Galerkin.jl

include("discontinuous_Galerkin.jl")
include("finite_volume_method.jl")
include("FRAP_approximation.jl")

include("distributions.jl")
include("extend_Base_operations.jl")
include("limiters.jl")

include("time_integration.jl")

include("simulate.jl")

end


text:
# https://juliacollections.github.io/IterTools.jl/stable/

using IterTools
using Test

function learn_distinct()
    data = [13, 5, 17, 3, 2, 2, 7, 5, 11, 13, 23]
    @test [13, 5, 17, 3, 2, 7, 11, 23] == collect(distinct(data))
end

@testset "Base" begin
    learn_distinct()
end;


text:
"""
```julia
function PulseDiff(; name, Val=1, dt=0.1)
```

# States
- `val(t)`
  Output value of the source

# Connectors
- `d`
  Output [`DigitalPin`](@ref)
"""
function PulseDiff(; name, Val=1, dt=0.1)
    @named d = DigitalPin()
    @variables val(t)
    D = ModelingToolkit.Difference(t; dt=dt)

    eqs = [
        D(val) ~ Val
        val ~ d.val
    ]

    ODESystem(eqs, t, [val], [], systems=[d], defaults=Dict(Val=>0), name=name)
end

"""
```julia
function Set(; name)
```

Source that outputs a constant signal of `1`.

# Connectors
- `d`
  Output [`DigitalPin`](@ref)
"""
function Set(; name)
    @named d = DigitalPin()

    eqs = [
        d.val ~ 1
    ]
    ODESystem(eqs, t, [], [], systems=[d],  name=name)
end

"""
```julia
function Reset(; name)
```

Source that outputs a constant signal of `1`

# Connectors
- `d`
  Output [`DigitalPin`](@ref)
"""
function Reset(; name)
    @named d = DigitalPin()

    eqs = [
        d.val ~ 0
    ]
    ODESystem(eqs, t, [], [], systems=[d], name=name)
end

"""
```julia
function Pulse(; name, duty_cycle=0.5, T=1.0)
```

Pulse output with specified `duty_cycle` and time period (`T`)

# Connectors
- `d`
  Output [`DigitalPin`](@ref)
"""
function Pulse(; name, duty_cycle=0.5, T=1.0)
    @named d = DigitalPin()

    eqs = [
        d.val ~ IfElse.ifelse(t%T > duty_cycle*T, 1, 0)
    ]
    ODESystem(eqs, t, [], [], systems=[d], name=name)
end


text:
const res_data_key = "data";

load_res_data(repo::ResRepoTree) =
    FileIO.load(repo.data_file)[res_data_key];

save_res_data!(repo::ResRepoTree, data) =
    FileIO.save(repo.data_file, res_data_key, data);


text:
println("Starting tests")

using Base.Test
using Evapotranspiration
using Dates

@test_approx_eq_eps convert(Radian, Degree(360)) 2π 0.0001
@test promote_type(Radian, Degree) == Radian

Evap = Evapotranspiration
latitude  = Degree(45.281258)
longitude = Degree(-75.575136)

debug=true
Tmax = Celcius(19.46)
Tmin = Celcius( 5.9 )
Tdew = Celcius( 6.25)
RHmean = RH(0.71)
latitude = Degree(45.1)

@test_approx_eq_eps convert(Kelvin, Tmax) 292.6 0.1

# Number of days through the year, J, days (1 = January 1)
J = Day(134)

# Lz longitude of the centre of the local time zone
Lz = Degree(75)

# Lm longitude of the measurement site [degrees west of Greenwich]
Lm = Degree(75.76498)

latent_heat_of_vaporaization = 2.54 #  MJ/kg

timestamp = 1368594000
d1 = unix2datetime(timestamp)
@test dayofyear(d1) == 135
d2 = Date(2013,2,27)
@test dayofyear(d2) == 58
d3 = Date(2012,12,31)
@test_approx_eq_eps dayofyear(d3) 366 1 # leap year

# Example 2:
@test_approx_eq_eps Evap.P(Meter(1800), Celcius(20)) 81.8 0.1
@test_approx_eq_eps Evap.γ(Meter(1800), Celcius(20)) 0.054 0.01
@test_approx_eq_eps Evap.γ(Meter(100)) 0.0642 0.001

@test_approx_eq_eps Evap.e°(Celcius(20)) 2.3383 0.0001
# Example 3:
@test_approx_eq_eps Evap.e°(Celcius(24.5)) 3.075 0.001
@test_approx_eq_eps Evap.e°(Celcius(15.0)) 1.705 0.001
@test_approx_eq_eps Evap.Es(Celcius(15), Celcius(24.5)) 2.39 0.01  # kPa
@test_approx_eq_eps Evap.Es(Celcius(17.5), Celcius(29)) 3.003 0.001  # kPa

## Saturation vapour pressure slope (Δ) in kPa/C:
@test_approx_eq_eps Evap.Δ(Celcius(28)) 0.2201 0.001 # kPa/C
@test_approx_eq_eps Evap.Δ(Celcius(26), Celcius(30)) 0.22008 0.00001 # kPa/C

## Actual vapour pressure (Ea) derived from dewpoint temperature
@test_approx_eq_eps Evap.Ea(Celcius(20)) 2.338 0.001 # kPa

# daily mean RH comes from forecast.io
# This is preferred over Ea(Tdew)
@test_approx_eq_eps Evap.Ea(Celcius(20), Celcius(20), RH(0.5)) 1.169 0.001 # kPa
@test_approx_eq_eps Evap.Ea(Celcius(18), Celcius(25), ave(RH(0.54),RH(0.82))) 1.7788 0.0001

# This is preferred over other methods
@test_approx_eq_eps Evap.Ea(Celcius(18), Celcius(25), RH(0.54), RH(0.82)) 1.7015 0.0001 # kPa

# Example 5
Tmin = Celcius(18); RHmax = RH(0.82)
Tmax = Celcius(25); RHmin = RH(0.54)
@test_approx_eq_eps Evap.e°(Tmin) 2.064 0.001
@test_approx_eq_eps Evap.e°(Tmax) 3.168 0.001
@test_approx_eq_eps Evap.Ea(Tmin, Tmax, RHmin, RHmax) 1.70 0.01
@test_approx_eq_eps Evap.Ea(Tmin, Tmax, ave(RHmin,RHmax)) 1.78 0.01
# Example 6
esa = Evap.Esa(Tmin, Tmax, RHmin, RHmax)
@test_approx_eq_eps esa 0.91 0.01

## Vapor pressure deficit, Es - Ea, kPa:
@test_approx_eq_eps Evap.Esa(Celcius(18), Celcius(25), Tdew=Celcius(17.8)) 0.5777 0.0001  # kPa
@test_approx_eq_eps Evap.Esa(Celcius(18), Celcius(25), RH(0.68)) 0.83708 0.00001  # kPa

Radian(2) == Radian(1)

# Example 8:
## Extraterrestrial radiation for daily periods, Ra, in MJ/(m^2*day)
sept3 = Date(2013,9,3)
latitude = Degree(-20) # degrees,
loc = Location(latitude, Lz, Lm)
ϕ = Evap.latitude(loc)
@test_approx_eq_eps ϕ -0.35 0.01
J = Day(dayofyear(sept3))
@test J == Day(246)
dᵣ = Evap.dᵣ(J)
@test_approx_eq_eps dᵣ 0.985 0.001
δ  = Evap.δ(J)
@test_approx_eq_eps δ  0.120 0.001
ωs = Evap.ωs(ϕ, δ)
@test_approx_eq_eps ωs 1.527 0.001
Ra = Evap.Ra(Day, J, loc) # MJ/(m^2*day)
@test_approx_eq_eps Ra 32.2 0.1 # MJ/(m^2*day)
Ra_mm_day = Evap.MJ2mm(Ra)
@test_approx_eq_eps Ra_mm_day 13.1 0.1 # mm/day

@test_approx_eq_eps Evap.ave(Hour, 0,1) 0.5 0.01
@test_approx_eq_eps Evap.ave(Hour, 24, 1) 24.5 0.01
@test_approx_eq_eps Evap.duration(Hour, 0,1) 1 0.01
@test_approx_eq_eps Evap.duration(Hour, 24, 1) 1 0.01
@test_approx_eq_eps Evap.duration(Hour, 24, 4) 4 0.01
@test_approx_eq_eps Evap.duration(Hour, 24, 0) 0 0.01
for hr in 0:6
  @test Evap.issunset(Hour, hr, J, loc) == true
end
for hr in 7:17
  @test Evap.issunset(Hour, hr, J, loc) == false
end
for hr in 18:24
  @test Evap.issunset(Hour, hr, J, loc) == true
end

# Example 9:
N = Evap.daylighthours(J, loc)
@test_approx_eq_eps N 11.7 0.1 # max daylight hours

t = 13.5
ωs = Evap.ωs(Hour, t, J, loc)
@test_approx_eq_eps ωs 0.385 0.001
# Extraterrestrial radiation for hourly periods, Ra, in mm/time_period)
hr1 = 13
hr2 = hr1+1
J = Day(134)  # dayofyear
ra_24hr = 0
step = 1
for hr1 in 0:step:24
  hr2 = (hr1+step)%24
  ra = Evap.MJ2mm(Evap.Ra(Hour, J, hr1, hr2, loc))
  ra_24hr += ra
end
@test_approx_eq_eps ra_24hr Evap.MJ2mm(Evap.Ra(Day, J, loc)) 0.2

# Example 10:
# The amount of radiation that penetrates the atmosphere, Rs, units same as Ra
latitude = Degree(-(22 + 54/60))
loc = Location(latitude, Lz, Lm)
J = Day(135)    # day of year
Ra = Evap.Ra(Day, J, loc)
@test_approx_eq_eps Ra 25.1 0.1
N = Evap.daylighthours(J, loc)
@test_approx_eq_eps N 10.9 0.1
# in May (31 days), 220 hours of sunshine were recorded
n = 220/31
@test_approx_eq_eps n 7.1 0.1
cloud_cover = Evap.cloudcover(n, J, loc)
@test_approx_eq_eps cloud_cover 0.348 0.001
z  = Meter(100)
as = 0.25
bs = 0.5
loc = Location(latitude, Lz, Lm, z, as, bs)
Rs = Evap.Rs(cloud_cover, J, Ra, loc)
@test_approx_eq_eps Rs 14.5 0.1
@test_approx_eq_eps Evap.MJ2mm(Rs) 5.9 0.1

Tmax = 25.1
Tmin = 19.1
Tdew = 18.1
ra = Evap.Ra(Day, J, loc)
cloud_cover = (1 - 7.1/10.9)
@test_approx_eq_eps Evap.MJ2mm(Evap.Rs(cloud_cover, J, ra, loc)) 5.90 0.01

# Example 11:
# Rio de Janeiro
latitude = Degree(-(22+54/60))
loc = Location(latitude, Lz, Lm, z, as, bs)
hrs_sunshine_per_month=220
n = hrs_sunshine_per_month/31
N = Evap.daylighthours(J, loc)
cloud_cover =  1 - n/N
@test_approx_eq_eps cloud_cover 0.349 0.001
Tmin = Celcius(19.1)
Tmax = Celcius(25.1)
Ea = kPa(2.1)  # kPa
Rso = Evap.Rso(Ra,as,bs)
@test_approx_eq_eps Rso 18.8 0.1
# Rs/Rso: Ratio of solar radiation that reaches the earth (compared to a clear sky day)
Rsso = Evap.Rsso(cloud_cover, J, Ra, loc)
@test_approx_eq_eps Rsso 0.77 0.01
## Net long wave radiation, Rnl, in mm/day
Rnl = Evap.Rnl(Rsso, Tmin, Tmax, Ea)
@test_approx_eq_eps Rnl 3.5 0.1

# Example 12:
## Net solar radiation, MJ/(m^2*day)
albedo = 0.23
loc = Location(latitude, Lz, Lm, z, as, bs, albedo)
Rns = Evap.Rns(Rs, loc)
@test_approx_eq_eps Rns 11.1 0.1
Rn = Evap.Rn(Rns, Rnl)
@test_approx_eq_eps Rn 7.6 0.1
@test_approx_eq_eps Evap.MJ2mm(Rn) 3.1 0.1
ea = kPa(2.1)
# Net radiation, Rn, in mm/day:
net_ra = Evap.net_radiation(Day, cloud_cover, J, Tmin, Tmax, RHmin, RHmax, loc)
@test_approx_eq_eps net_ra 7.2 0.1

hr1 = 12
hr2 = hr1+1
net_ra_hr = Evap.net_radiation(Hour, cloud_cover, J, Tmin, Tmax, RHmin, RHmax, hr1, hr2, loc)
@test_approx_eq_eps net_ra_hr -2.363 0.001

## For a day depth is between 0.1 and 0.2, 1 month is about 1, 4 months can be 2.0.
shf = Evap.soil_heat_flux(Day, 30, Celcius(14.1), Celcius(18.8))
@test_approx_eq_eps shf 0.08451 0.00001
@test_approx_eq_eps Evap.soil_heat_flux(Day, 1, Tmin, Tmax)  0.6048 0.0001
@test_approx_eq_eps Evap.soil_heat_flux(Day, 1, Tmax, Tmin) -0.6048 0.0001
# Example 13:
# Soil is warming
T1 = Celcius(14.1) # March
T2 = Celcius(16.1) # April
T3 = Celcius(18.8) # May
@test_approx_eq_eps Evap.soil_heat_flux(Month, 2, T1, T3) 0.33 0.01
@test_approx_eq_eps Evap.soil_heat_flux(Month, 1, T1, T2) 0.28 0.01

# Example 14:
## Windspeed height compensation (for weather stations)
# Find equiv windspeed at 2m
uz = 3.2  # wind speed at 10 m in m/s
z  = 10
u2 = Evap.u2(uz, z)
@test_approx_eq_eps u2 2.4 0.01

## Brussels, 6 July located at 50°48'N and at 100 m above sea level:
latitude = Degree(50 + 48/60)
z = Meter(100)
loc = Location(latitude, Lz, Lm, z, as, bs, albedo)
Tmin = Celcius(12.3)
Tmax = Celcius(21.5)
RHmin = RH(0.63)
RHmax = RH(0.84)
RHmean = ave(RHmin,RHmax)
P_kPa = kPa(100.1)
windspeed_10m_kmph = 10
cloud_cover = 0.43
J = Day(187)
eto_day = Evap.ETo(Day, 1, Tmin, Tmax, RHmin, RHmax, P_kPa, windspeed_10m_kmph, cloud_cover, J::Day, loc)
@test_approx_eq_eps eto_day 3.7 0.01  # mm/day

# Example 17:
# Given the monthly average climatic data of April of Bangkok (Thailand)
# located at 13°44'N and at an elevation of 2 m
latitude = Degree(13+44/60)
z  = Meter(2)  # elevation above sea level
# Lz longitude of the centre of the local time zone
Lz = Degree(100+29/60)
# Lm longitude of the measurement site [degrees west of Greenwich]
Lm = Lz
loc = Location(latitude, Lz, Lm, z, as, bs, albedo)
dt = Date(2015, 4, 15)
J = Day(dayofyear(dt))
Tmin = Celcius(25.6)
Tmax = Celcius(34.8)
u2 = 2
T2 = Celcius(30.2)  # mean monthly average temp, April
T1 = Celcius(29.2)  # mean monthly average temp, March
P_kPa = Evap.P(z, ave(T1,T2))
@test_approx_eq_eps P_kPa 101.3 0.1
n  = 8.5 # hrs/day of sunshine
cloud_cover = Evap.cloudcover(n, J, loc)
# Need to figure out what to do about RH values: TODO
eto_day = Evap.ETo(Day, 1, Tmin, Tmax, RHmin, RHmax, P_kPa, windspeed_10m_kmph, cloud_cover, J::Day, loc)

# Todo:

#function print_evaporation(forecast_data) {
#    vec = forecast_data
#    totalDewPoint = 0
#    totalRH = 0
#    numDewPoints = 0
#    for (i=0; i < vec.length; i++) {
#        v = vec[i]
#	v.date = Date(v.time*1000)
#	v.dom = v.date.getDate()
#	v.hr = v.date.getHours()
#        ave_temp_c = NaN
#        if ("temperature" in v) {
#            ave_temp_c = v.temperature
#        } else if (("temperatureMax" in v) && ("temperatureMin" in v)) {
#            ave_temp_c = v.temperatureMin + 0.75*(v.temperatureMax - v.temperatureMin)
#        }
#        dew_temp_c = v.dewPoint
#        solar_radiation = 6
#        sunshine = 1-v.cloudCover; # account for night having no sun.
#        wind_speed_kmph = v.windSpeed
#        E = evaporation_mm_per_day(wind_speed_kmph, ave_temp_c, dew_temp_c)
#        H = degree_days_mm_per_day(ave_temp_c, dew_temp_c, 6, sunshine)
#                U = potential_evaporation(ave_temp_c, dew_temp_c, wind_speed_kmph, solar_radiation, sunshine)
#        precipIntensity = NaN
#        if ("percipIntensityMax" in v) {
#            precipIntensity = v.percipIntensityMax
#        }
#        if ("precipIntensity" in v) {
#            precipIntensity = v.precipIntensity
#        }
#        percip_mm_day = v.precipIntensity; # * vec[i].precipProbability
#        if (v.dom == 18) {
#		console.log(v.date+': '+v.humidity+', '+v.dewPoint)
#		totalDewPoint += v.dewPoint
#		totalRH += v.humidity
#		numDewPoints += 1
#	}
#        if (typeof percip_mm_day == 'undefined') {
#            percip_mm_day = 0
#        }
#        AM = percip_mm_day - U
#        t = Date(1000*vec[i].time)
#        #console.log(vec[i].time+": "+t + ": E="+ E.toFixed(2) + ", H=" + H.toFixed(2)+", ETo="+U.toFixed(2)+", percip="+percip_mm_day.toFixed(2)+", AM="+AM.toFixed(2))
#
#	# New method:
#	Tmin = v.temperatureMin
#	Tmax = v.temperatureMax
#	Tdew = v.dewPoint
#	RHmean = v.humidity
#	P_kPa = v.pressure
#	# TODO: scale down wind because of shelter from trees, maybe 0.25
#	windspeed_10m_kmph = v.windSpeed
#	cloud_cover = v.cloudCover
#	dayofyear = dayofyear(v.time)
#        eto_day = ETo_day(Tmin, Tmax, RHmean, P_kPa, windspeed_10m_kmph, cloud_cover, latitude, dayofyear)
#	eto_diff = U - eto_day
#        #console.log("ETo day="+eto_day.toFixed(2)+" [mm/day]   diff="+eto_diff)
#    }
#    console.log("Total dew: "+ totalDewPoint+" RH: "+totalRH)
#    aveDewPoint = totalDewPoint/numDewPoints
#    aveRH = totalRH/numDewPoints
#    console.log("Ave: Dew Point: "+ aveDewPoint+", RH: "+ aveRH)
#}
#
#
#function print_daily_forecast(data) {
#  daily = data.daily.data
#  for (i=0; i<daily.length; i++) {
#      di = daily[i]
#      console.log(Date(1000*di.time) + ":")
#      console.log("  - sunrise: " + Date(1000*di.sunriseTime))
#      console.log("  - sunset:  " + Date(1000*di.sunsetTime))
#      console.log("  - TmaxTime: " + Date(1000*di.temperatureMaxTime))
#      console.log("  - TminTime:  " + Date(1000*di.temperatureMinTime))
#      console.log(di)
#  }
#}
#
#forecast.get(latitude, longitude, options, function (err, res, data) {
#  if (err) throw err
#  #console.log('res: ' + util.inspect(res))
#  #console.log('data: ' + util.inspect(data, true, 5))
#  #console.log('daily data: ' + util.inspect(data, true, 5))
#  #print_daily_forecast(data)
#  print_evaporation(data.daily.data)
#  print_evaporation(data.hourly.data)
#})
#
#

println("Done")


text:
module EachindexArrays

using Swizzles
using Swizzles.ShallowArrays
using Swizzles.ArrayifiedArrays
using Swizzles.WrapperArrays

using Base.Cartesian

using Swizzles: assign!, increment!

export EachindexArray, CartesianTiledIndices, Tile

struct EachindexArray{T, N, Arg, Indices} <: ShallowArray{T, N, Arg}
    arg::Arg
    indices::Indices
end

@inline EachindexArray(arg::AbstractArray, indices) = EachindexArray{eltype(arg), ndims(arg), typeof(arg), typeof(indices)}(arg, indices)

@inline Base.parent(arr::EachindexArray) = arr.arg
@inline WrapperArrays.iswrapper(arr::EachindexArray) = true
@inline WrapperArrays.adopt(arr::EachindexArray, arg) = EachindexArray(arg, arr.indices)

@inline function Base.eachindex(arr::EachindexArray)
    return arr.indices
end

@inline function Base.eachindex(arr::EachindexArray, args::AbstractArray...)
    return arr.indices
end

struct CartesianTiledIndices{N, Indices <: CartesianIndices{N}, Chunks <: Tuple{Vararg{Any, N}}} <: ShallowArray{CartesianIndex{N}, N, Indices}
    indices::Indices
    chunks::Chunks
    function CartesianTiledIndices{N, Indices, Chunks}(indices, chunks) where {N, Indices, Chunks}
        return new{N, Indices, Chunks}(indices, chunks)
    end
end

@inline Base.parent(arr::CartesianTiledIndices) = arr.indices
@inline WrapperArrays.iswrapper(arr::CartesianTiledIndices) = true

@inline function CartesianTiledIndices(indices::Indices, chunks::Chunks) where {Indices, Chunks}
    CartesianTiledIndices{ndims(indices), Indices, Chunks}(indices, chunks)
end

@generated function Swizzles.assign!(dst, index, src, drive::CartesianTiledIndices{N, Indices, chunks}) where {N, Indices, chunks}
    return quote
        Base.@_propagate_inbounds_meta
        @nloops $N i n -> 1:drive.chunks[n]:length(drive.indices.indices[n]) begin
            tile = @ntuple $N n -> drive.indices.indices[n][i_n:min(i_n + drive.chunks[n] - 1, length(drive.indices.indices[n]))]
            assign!(dst, index, src, CartesianIndices(tile))
        end
    end
end

@generated function Swizzles.increment!(op::Op, dst, index, src, drive::CartesianTiledIndices{N, Indices, chunks}) where {Op, N, Indices, chunks}
    return quote
        Base.@_propagate_inbounds_meta
        @nloops $N i n -> 1:drive.chunks[n]:length(drive.indices.indices[n]) begin
            tile = @ntuple $N n -> drive.indices.indices[n][i_n:min(i_n + drive.chunks[n] - 1, length(drive.indices.indices[n]))]
            increment!(op, dst, index, src, CartesianIndices(tile))
        end
    end
end

struct Tile{_Chunks} <: Swizzles.Intercept
    _chunks::_Chunks
end

@inline Tile(_chunks...) = Tile{typeof(_chunks)}(_chunks)

@inline (ctr::Tile)(arg) = ctr(arrayify(arg))
@inline function(ctr::Tile{<:Tuple{Vararg{Any, _N}}})(arg::Arg) where {_N, Arg <: AbstractArray}
    if @generated
        chunks = ntuple(n -> n <= _N ? :(ctr._chunks[$n]) : 1, ndims(arg))
        return :(return EachindexArray(arg, CartesianTiledIndices(CartesianIndices(arg), ($(chunks...),))))
    else
        chunks = ntuple(n -> n <= _N ? ctr._chunks[n] : 1, ndims(arg))
        return EachindexArray(arg, CartesianTiledIndices(CartesianIndices(arg), chunks))
    end
end

end


text:
module CellArrayApplyTests

using Test
using Gridap
using Gridap.CachedArrays
using ..CellValuesMocks
using TensorValues
using Gridap.Kernels: IntegrateNumberKernel
using Gridap.Kernels: IntegrateArrayKernel

l = 10

a = [1,2,3]
v = TestIterCellValue(a,l)
w = apply(-,v,broadcast=true)
o = [ CachedArray(-vi) for vi in v ]
test_iter_cell_array(w,o)

oo = collect(w)
@test oo == o

a = [1,2,3]
b = [3,2,1]
u = TestIterCellValue(a,l)
v = TestIterCellValue(b,l)
w = apply(-,u,v,broadcast=true)
o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
test_iter_cell_array(w,o)

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)

ax = [rand(1,3,4), 1.0      , [v1,v2,v3], [v1,v2,v3], 1         ]
bx = [rand(2,3,1), rand(2,3), [v2,v3,v1], v1        , [v2,v3,v1]]

for (a,b) in zip(ax,bx)
  u = TestIterCellValue(a,l)
  v = TestIterCellValue(b,l)
  o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
  w = apply(-,u,v,broadcast=true)
  test_iter_cell_array(w,o)
end

a = [1,2,3]
v = TestIndexCellValue(a,l)
w = apply(-,v,broadcast=true)
o = [ CachedArray(-vi) for vi in v ]
test_index_cell_array(w,o)

a = [1,2,3]
b = [3,2,1]
u = TestIndexCellValue(a,l)
v = TestIndexCellValue(b,l)
w = apply(-,u,v,broadcast=true)
o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
test_index_cell_array(w,o)

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)

ax = [rand(1,3,4), 1.0      , [v1,v2,v3], [v1,v2,v3], 1         ]
bx = [rand(2,3,1), rand(2,3), [v2,v3,v1], v1        , [v2,v3,v1]]

for (a,b) in zip(ax,bx)
  u = TestIndexCellValue(a,l)
  v = TestIndexCellValue(b,l)
  o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
  w = apply(-,u,v,broadcast=true)
  test_index_cell_array(w,o)
end

a = [1,2,3]
b = [3,2,1]
c = [3,2,4]
u = TestIndexCellValue(a,l)
v = TestIndexCellValue(b,l)
w = TestIndexCellValue(c,l)
k = IntegrateNumberKernel{Int}()
z = apply(k,u,v,w)
r = fill(sum(a.*b.*c),l)
test_index_cell_number(z,r)

a = rand(3,2,4)
b = rand(4)
c = rand(4)
u = TestIndexCellValue(a,l)
v = TestIndexCellValue(b,l)
w = TestIndexCellValue(c,l)
k = IntegrateArrayKernel()
z = apply(k,u,v,w)
ri =  reshape( sum(a .* reshape(b,(1,1,4)) .* reshape(c,(1,1,4)), dims=3), (3,2))
r = fill(collect(ri),l)
test_index_cell_array(z,r)

end # module