text:
module Schedule
    abstract type AbstractSchedule end

    abstract type PeriodicSchedule <: AbstractSchedule end

end

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


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

using Markdown
using InteractiveUtils

# ╔═╡ d7ff772c-9d33-42e7-9df1-3e02ceee79c5
begin
	using PlutoUI
	PlutoUI.TableOfContents()
end

# ╔═╡ 7553b887-a16d-4f3b-a91b-043f89e963d7
md"""
# Control flow

- [Julia docs | control flow](https://docs.julialang.org/en/v1/manual/control-flow/)
- [Julia docs | functions](https://docs.julialang.org/en/v1/manual/functions/)

Julia programs are able to run nonlinearly by controlling its execution flow.
"""

# ╔═╡ 383e7736-dd64-44ce-ac21-60a9a5a7c0db
md"""
## Conditional statements

Compound expressions are already covered in the previous chapter: [intro to Julia](intro-to-julia.html)

### `if`/`else` block:

```julia
result = if cond1
	run_if_cond1_is_true()
elseif cond2
	run_if_cond2_is_true_and_cond1_is_false()
else
	run_if_neither_is_true()
end
```

- The `elseif` and `else` blocks are optional. `if` and `end` are mandatory.
- `if` blocks return a value. Capturing the value is optional.
- `if` blocks are "leaky", i.e. they do not introduce a local scope. (The same as Python)
- Only boolean (true / false) could be evaluated in `if` blocks. Using other types (e.g. Int) will generate an error.

### Ternary operator

`cond ? tvalue : fvalue` is eqialent to:

```julia
if cond
	tvalue
else
	fvalue
end
```

### ifelse function

All the arguments are evaluated first in `ifelse(cond, tvalue, fvalue)`.

### short cicuit evaluation

`&&` (logical and) and `||` (logical or) operators support [short circuit evaluation](
Short-circuit evaluation - Wikipedia
https://en.wikipedia.org › wiki › Short-circuit_evaluation).

- In the expression `a && b`, the subexpression `b` is only evaluated if `a` evaluates to true.
- In the expression `a || b`, the subexpression `b` is only evaluated if `a` evaluates to false.
"""

# ╔═╡ 0335796f-da7d-4acc-be19-a0b45cd6042a
true && 1  # && evaluates and returns the second argument if the first is true

# ╔═╡ 2ee995ac-d9e8-4ded-b6b9-5fd35d564a24
false && 1  # && otherwise returns false

# ╔═╡ ce2cbdb6-637d-4f07-bb91-580621485b66
score = 10

# ╔═╡ 648e6be1-5e8b-4406-bee5-708a93ea20f3
response = if 80 < score <= 100
  "Good"
elseif 60 < score <= 80
  "Okay"
else
  "Uh-Oh"
end

# ╔═╡ 2732609d-86f3-4594-8f02-4f8066d41720
md"""
## Loops

Loops are repeated evaluations.

- `while` loops are often controlled by evaluated conditions.
- `for` loops are often controlled by sequence length.

- `break`: exit the loop immediately.
- `continue`: move on to the next item / evaluation immediately.
"""

# ╔═╡ af291389-2e9c-41fc-880b-28fe1fe70e49
## Hailstone sequence (3n+1 problem) in a while loop
let n = 27
    step = 0
    while n != 1
        if iseven(n)
            n ÷= 2
        else
            n = 3n + 1
        end
        step += 1
    end
    step
end

# ╔═╡ 4bb1e57d-ab16-4c97-ad4b-3db6a76b7e11
let sum = 0, n = 100
	for i in 1:n
		sum += i
	end
	sum
end

# ╔═╡ 0117590e-213b-4a70-a9c8-f44075c58854
with_terminal() do
	for x in 1:9
		if x == 5
			continue # jump to line #2
		elseif x >=8
			break    # jump to line #9
		end
		println(x, "^2 = ", x^2)
	end
end

# ╔═╡ d52dfdce-cd96-4d1b-8aa6-f1137758221b
with_terminal() do
	# Use enumerate(seq) to get a pair of one idx and one element
	for (i, x) in enumerate(10:-1:1)
		println("xs[", i, "] = ", x)
	end
end

# ╔═╡ 81194437-5b15-4ab9-af47-33108dde7751
md"""
Multiple nested for loops can be combined into a single outer loop, forming the cartesian product of its iterables.
"""

# ╔═╡ 21c53be2-3d98-4fb5-87b5-dd53bac7dd76
with_terminal() do
	# Multiple nested for loops can be combined into a single outer loop, forming the cartesian product of its iterables
	for i = 'x':'z', j = '1':'3'
		println(i, j)
	end
end

# ╔═╡ bd4a9500-5b35-4f0a-bd1c-4e42a650308b
md"""
# Functions

> In Julia, a function is an object that maps a tuple of argument values to a return value.
> [Julia docs](https://docs.julialang.org/en/v1/manual/functions/)

Functions facilitate:
- Code reuse and encapsulation.
- Specializations ([Methods](https://docs.julialang.org/en/v1/manual/methods/))

- Functions are first-class objects and can be passed into a higher-order function.
- The arguments are "passed-by-sharing". Modifications to mutable argument values (such as `Arrays`) will be reflected to the caller. (Similar to Python)
- Functions that will update the arguments are named with a bang `!` by convention. (e.g. sort(arr) vs sort!(arr))
- Often only the scalar version of a function is required; for vector element-wise operations, there are broadcast (dot) syntax.
- You can write multiple function with the same name provided they have different parameter lists. Julia will choose the most apporpriate one for you.
"""

# ╔═╡ 0d1252a0-3ce8-4b96-ba08-a6f1895cc0dc
md"""
## Standard form
"""

# ╔═╡ 5d26cdd5-9a6d-4307-89b7-6629997b9132
"Mechaelis-Menton function"  # Function documentations
function mm(x, k)            # function name and parameter list
  result = x / (x +k)
  return result              # return statement is optional
end

# ╔═╡ ce413514-bf5c-4dc8-b2af-42e9f663fa87
md"""
## One-liner form
"""

# ╔═╡ b87e79d8-91b0-4ecb-85b8-e2e8d01da426
# One-liner function
mm(x) = mm(x, one(x))

# ╔═╡ 5d2c7cdc-d34a-4e97-82c4-ac58a9196fd8
mm(0.5, 2)

# ╔═╡ d137d92d-f629-41cf-824c-516391f26f11
mm(1)

# ╔═╡ 65de3953-6a86-4e3f-ba00-b6c352678284
md"""
## Anonymous functions

Often used with high-order functions e.g. `map()`
"""

# ╔═╡ f14ab8cb-db2e-47f7-adaa-d932ed25bb1b
map(x->x^2, 1:3)

# ╔═╡ 26ba95b5-1d1c-4bed-b6dd-e17eba8136e3
md"""
Use [`do` block](https://docs.julialang.org/en/v1/manual/functions/#Do-Block-Syntax-for-Function-Arguments) for long anonymous functions.

`do` block are also useful in opening files.

```julia
open("outfile", "w") do io
    write(io, data)
end
```
"""

# ╔═╡ e2801350-24b7-4453-af27-ee63afbb4087
val = rand(-6:6, 10)

# ╔═╡ 2f07f81c-346a-4f17-872c-f174adf6ff1d
map(val) do x
    if x < 0 && iseven(x)
        return 0
    elseif x == 0
        return 1
    else
        return x
    end
end

# ╔═╡ ea04a83f-74a1-46b0-b686-f125e0972ca2
md"""
## Optional arguments

[Optional (positional) arguments](https://docs.julialang.org/en/v1/manual/functions/#Optional-Arguments) are listed after mandatory ones.

```julia
function func(a, b, c=1)
# content
end
```

And they are called with `func(a, b)` or `func(a, b, 3)`

## Keyword arguments

[Keyword arguments](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments) are listed after `;`

```julia
function plot(x, y; style="solid", width=1, color="black")
    ###
end
```

And they are called with `plot(x, y, width=2)` or `plot(x, y; width=2)`
"""

# ╔═╡ d2e2e668-71f5-444b-98bf-5ea45df81c84
args_kwargs(args...; kwargs...) = (args, kwargs)  # mind the semicolon ;

# ╔═╡ 91e0412a-24ef-469a-9a67-3dcfa155d62a
args_kwargs(1,true, :IAmASymbol, b=4, c="IAmAString")

# ╔═╡ 33a48540-cb59-11eb-0602-8da0ecdcbe4c
md"""
# See also

- [Compositing functions and pipes](https://docs.julialang.org/en/v1/manual/functions/#Function-composition-and-piping)
- [Variable argument (vararg) functions](https://docs.julialang.org/en/v1/manual/functions/#Varargs-Functions)
- [Argument destructuring](https://docs.julialang.org/en/v1/manual/functions/#Argument-destructuring)
- [Scope of variables](https://docs.julialang.org/en/v1/manual/variables-and-scoping/)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═d7ff772c-9d33-42e7-9df1-3e02ceee79c5
# ╠═7553b887-a16d-4f3b-a91b-043f89e963d7
# ╠═383e7736-dd64-44ce-ac21-60a9a5a7c0db
# ╠═0335796f-da7d-4acc-be19-a0b45cd6042a
# ╠═2ee995ac-d9e8-4ded-b6b9-5fd35d564a24
# ╠═ce2cbdb6-637d-4f07-bb91-580621485b66
# ╠═648e6be1-5e8b-4406-bee5-708a93ea20f3
# ╠═2732609d-86f3-4594-8f02-4f8066d41720
# ╠═af291389-2e9c-41fc-880b-28fe1fe70e49
# ╠═4bb1e57d-ab16-4c97-ad4b-3db6a76b7e11
# ╠═0117590e-213b-4a70-a9c8-f44075c58854
# ╠═d52dfdce-cd96-4d1b-8aa6-f1137758221b
# ╠═81194437-5b15-4ab9-af47-33108dde7751
# ╠═21c53be2-3d98-4fb5-87b5-dd53bac7dd76
# ╠═bd4a9500-5b35-4f0a-bd1c-4e42a650308b
# ╠═0d1252a0-3ce8-4b96-ba08-a6f1895cc0dc
# ╠═5d26cdd5-9a6d-4307-89b7-6629997b9132
# ╠═ce413514-bf5c-4dc8-b2af-42e9f663fa87
# ╠═b87e79d8-91b0-4ecb-85b8-e2e8d01da426
# ╠═5d2c7cdc-d34a-4e97-82c4-ac58a9196fd8
# ╠═d137d92d-f629-41cf-824c-516391f26f11
# ╠═65de3953-6a86-4e3f-ba00-b6c352678284
# ╠═f14ab8cb-db2e-47f7-adaa-d932ed25bb1b
# ╠═26ba95b5-1d1c-4bed-b6dd-e17eba8136e3
# ╠═e2801350-24b7-4453-af27-ee63afbb4087
# ╠═2f07f81c-346a-4f17-872c-f174adf6ff1d
# ╠═ea04a83f-74a1-46b0-b686-f125e0972ca2
# ╠═d2e2e668-71f5-444b-98bf-5ea45df81c84
# ╠═91e0412a-24ef-469a-9a67-3dcfa155d62a
# ╠═33a48540-cb59-11eb-0602-8da0ecdcbe4c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
##############################################################

"""
  function fasta2matrix(filename::AbstractString; max_gap_fraction = 1)

  Reads an aligned fasta file and returns a matrix with amino acids converted
  to numbers following the convention A --> 1, C --> 2, ... '-' --> 21.

    "filename": full path of the fasta file.
    "max_gap_fraction": to be removed.

"""

function fasta2matrix(filename::AbstractString; max_gap_fraction = 1)
    f = FastaReader(filename)

    max_gap_fraction = Float64(max_gap_fraction)

    # pass 1

    seqs = Int[]
    inds = Int[]
    fseqlen = 0

    for (name, seq) in f
        ngaps = 0
        if f.num_parsed == 1
            ls = length(seq)
            resize!(inds, ls)
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    fseqlen += 1
                    inds[fseqlen] = i
                    c == '-' && (ngaps += 1)
                end
            end
        else
            ls = length(seq)
            ls == length(inds) || error("inputs are not aligned")
            tstfseqlen = 0
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    tstfseqlen += 1
                    inds[tstfseqlen] == i || error("inconsistent inputs")
                    c == '-' && (ngaps += 1)
                end
            end
            tstfseqlen == fseqlen || error("inconsistent inputs")
        end
        ngaps / fseqlen <= max_gap_fraction && push!(seqs, f.num_parsed)
    end

    length(seqs) > 0 || error("Out of $(f.num_parsed) sequences, none passed the filter (max_gap_fraction=$max_gap_fraction)")

    # pass 2
    Z = Array{Int16}(undef, fseqlen, length(seqs))
    seqid = 1
    for (name, seq) in f
        seqs[end] < f.num_parsed && break
        seqs[seqid] == f.num_parsed || continue
        for i = 1:fseqlen
            c = seq[inds[i]]
            Z[i, seqid] = letter2num(c)
        end
        seqid += 1
    end
    @assert seqid == length(seqs) + 1

    close(f)

    return Z'
end

##############################################################
"""
    letter2num(c::Union{Char,UInt8})

    Takes as input a char (representing an amino acid)
    and returns its conventional number representation.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

let alphabet = [ 1, 21, 2, 3, 4, 5, 6, 7, 8, 21, 9, 10, 11, 12, 21, 13, 14, 15, 16, 17, 21, 18, 19, 21, 20]
               # A, B,  C, D, E, F, G, H, I, J,  K, L,  M,  N,  O,  P,  Q,  R,  S,  T,  U,  V,  W,  X,  Y
    global letter2num
    function letter2num(c::Union{Char,UInt8})
        i = UInt8(c) - 0x40
        1 <= i <= 25 && return alphabet[i]
        return 21
     end
end

##############################################################
"""
    num2letter(i::Integer)

    Takes as input an integer (representing an amino acid)
    and returns its conventional character representation.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

let alphabet = ["A", "C", "D", "E", "F", "G", "H", "I",  "K", "L",  "M",  "N", "P",  "Q",  "R",
"S",  "T", "V",  "W",  "Y"]
    global num2letter
    function num2letter(i :: Integer)
        1 <= i <= 20 && return alphabet[i]
        return "-"
    end
end

##############################################################
"""
    vec2string(v)

    Takes as input a vector of integers (representing an amino acids)
    and returns a list of characters, the corresponding amino acids.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

function vec2string(v)
    s = ""
    for i in v
        s = s*num2letter(i)
    end
    return s
end

##############################################################
"""
    Takes as input a string (representing amino acids)
    and returns a vector of numbers, the corresponding number representation.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

function string2vec(s)
    v = Vector{Int8}(undef, length(s))
    for (i, l) in enumerate(s)
        v[i] = letter2num(l)
    end
    return v
end

max_stars_repo_name: sabrinacotogno/KitMSA

max_issues_repo_name: sabrinacotogno/KitMSA

max_forks_repo_name: sabrinacotogno/KitMSA


text:
# This module extends the julia markdown parser to improve compatibility with Jupyter, Pandoc etc.
module WeaveMarkdown

using ..Weave: isnothing, take2string!
using Markdown
import Markdown: @trigger, @breaking, Code, MD, withstream, startswith, LaTeX

# HACK: that this definition is type-piracy. It is required since `Markdown`
# does not have a built in system for contextual rendering by users. `io` here
# should always be either `IOBuffer` or `IOContext` since it is reached via
# `sprint` in all user-facing code paths in `Markdown`.
function Markdown.latex(io::Union{IOBuffer,IOContext}, tex::Markdown.LaTeX)
    math_envs = ["align", "equation", "eqnarray"]
    use_dollars =
        !any([occursin("\\begin{$me", tex.formula) for me in math_envs])
    use_dollars && write(io, "\\[")
    write(io, string("\n", tex.formula, "\n"))
    use_dollars && write(io, "\\]\n")
end

mutable struct Comment
    text::String
end

@breaking true -> function dollarmath(stream::IO, block::MD)
    withstream(stream) do
        str = Markdown.startswith(stream, r"^\$\$$"m)
        isempty(str) && return false
        trailing = strip(readline(stream))
        buffer = IOBuffer()
        while !eof(stream)
            line_start = position(stream)
            estr = Markdown.startswith(stream, r"^\$\$$"m)
            if !isempty(estr)
                estr = Markdown.startswith(stream, r"^\$\$$"m)
                if isempty(estr)
                    push!(block, LaTeX(take2string!(buffer) |> chomp))
                end
                return true
            else
                seek(stream, line_start)
            end
            write(buffer, readline(stream, keep = true))
        end
        return false
    end
end

@breaking true -> function topcomment(stream::IO, block::MD)
    buffer = IOBuffer()
    withstream(stream) do
        str = Markdown.startswith(stream, r"^<!--")
        isempty(str) && return false
        while !eof(stream)
            line = readline(stream, keep = true)
            write(buffer, line)
            if occursin(r"-->$", line)
                s = replace(take2string!(buffer) |> chomp, r"-->$" => "")
                push!(block, Comment(s))
                return true
            end
        end
        return false
    end
end

@trigger '<' -> function comment(stream::IO, md::MD)
    withstream(stream) do
        Markdown.startswith(stream, "<!--") || return
        text = Markdown.readuntil(stream, "-->")
        isnothing(text) && return
        return Comment(text)
    end
end

# Create own flavor and copy all the features from julia flavor
Markdown.@flavor weavemd [dollarmath, comment, topcomment]
weavemd.breaking = [weavemd.breaking; Markdown.julia.breaking]
weavemd.regular = [weavemd.regular; Markdown.julia.regular]
for key in keys(Markdown.julia.inner)
    if haskey(weavemd.inner, key)
        weavemd.inner[key] = [weavemd.inner[key]; Markdown.julia.inner[key]]
    else
        weavemd.inner[key] = Markdown.julia.inner[key]
    end
end

include("html.jl")
include("latex.jl")

end  # module


text:
export DryModel, EquilMoist, NonEquilMoist

#### Moisture component in atmosphere model
abstract type MoistureModel end

vars_state(::MoistureModel, ::AbstractStateType, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::MoistureModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_moisture!(
    ::MoistureModel,
    ::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    ::MoistureModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end
function flux_second_order!(
    ::MoistureModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_t,
) end
function flux_first_order!(::MoistureModel, _...) end
function compute_gradient_argument!(
    ::MoistureModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

internal_energy(atmos::AtmosModel, state::Vars, aux::Vars) =
    internal_energy(atmos.moisture, atmos.orientation, state, aux)

@inline function internal_energy(
    moist::MoistureModel,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    Thermodynamics.internal_energy(
        state.ρ,
        state.ρe,
        state.ρu,
        gravitational_potential(orientation, aux),
    )
end

"""
    DryModel

Assumes the moisture components is in the dry limit.
"""
struct DryModel <: MoistureModel end

vars_state(::DryModel, ::Auxiliary, FT) = @vars(θ_v::FT, air_T::FT)
@inline function atmos_nodal_update_auxiliary_state!(
    moist::DryModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ts = new_thermo_state(atmos, state, aux)
    aux.moisture.θ_v = virtual_pottemp(ts)
    aux.moisture.air_T = air_temperature(ts)
    nothing
end

"""
    EquilMoist

Assumes the moisture components are computed via thermodynamic equilibrium.
"""
struct EquilMoist{FT} <: MoistureModel
    maxiter::Int
    tolerance::FT
end
EquilMoist{FT}(;
    maxiter::IT = 3,
    tolerance::FT = FT(1e-1),
) where {FT <: AbstractFloat, IT <: Int} = EquilMoist{FT}(maxiter, tolerance)

vars_state(::EquilMoist, ::Prognostic, FT) = @vars(ρq_tot::FT)
vars_state(::EquilMoist, ::Gradient, FT) = @vars(q_tot::FT)
vars_state(::EquilMoist, ::GradientFlux, FT) = @vars(∇q_tot::SVector{3, FT})
vars_state(::EquilMoist, ::Auxiliary, FT) =
    @vars(temperature::FT, θ_v::FT, q_liq::FT, q_ice::FT)

@inline function atmos_nodal_update_auxiliary_state!(
    moist::EquilMoist,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ts = new_thermo_state(atmos, state, aux)
    aux.moisture.temperature = air_temperature(ts)
    aux.moisture.θ_v = virtual_pottemp(ts)
    aux.moisture.q_liq = PhasePartition(ts).liq
    aux.moisture.q_ice = PhasePartition(ts).ice
    nothing
end

function compute_gradient_argument!(
    moist::EquilMoist,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.moisture.q_tot = state.moisture.ρq_tot * ρinv
end

function compute_gradient_flux!(
    moist::EquilMoist,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    # diffusive flux of q_tot
    diffusive.moisture.∇q_tot = ∇transform.moisture.q_tot
end

function flux_moisture!(
    moist::EquilMoist,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    u = state.ρu / ρ
    flux.moisture.ρq_tot += u * state.moisture.ρq_tot
end

function flux_second_order!(
    moist::EquilMoist,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_t,
)
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    flux_second_order!(moist, flux, state, d_q_tot)
end
#TODO: Consider whether to not pass ρ and ρu (not state), foc BCs reasons
function flux_second_order!(moist::EquilMoist, flux::Grad, state::Vars, d_q_tot)
    flux.ρ += d_q_tot * state.ρ
    flux.ρu += d_q_tot .* state.ρu'
    flux.moisture.ρq_tot += d_q_tot * state.ρ
end

"""
    NonEquilMoist

Does not assume that the moisture components are in equilibrium.
"""
struct NonEquilMoist <: MoistureModel end

vars_state(::NonEquilMoist, ::Prognostic, FT) =
    @vars(ρq_tot::FT, ρq_liq::FT, ρq_ice::FT)
vars_state(::NonEquilMoist, ::Gradient, FT) =
    @vars(q_tot::FT, q_liq::FT, q_ice::FT)
vars_state(::NonEquilMoist, ::GradientFlux, FT) = @vars(
    ∇q_tot::SVector{3, FT},
    ∇q_liq::SVector{3, FT},
    ∇q_ice::SVector{3, FT}
)
vars_state(::NonEquilMoist, ::Auxiliary, FT) = @vars(temperature::FT, θ_v::FT)

@inline function atmos_nodal_update_auxiliary_state!(
    moist::NonEquilMoist,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ts = new_thermo_state(atmos, state, aux)
    aux.moisture.temperature = air_temperature(ts)
    aux.moisture.θ_v = virtual_pottemp(ts)
    nothing
end

function compute_gradient_argument!(
    moist::NonEquilMoist,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.moisture.q_tot = state.moisture.ρq_tot * ρinv
    transform.moisture.q_liq = state.moisture.ρq_liq * ρinv
    transform.moisture.q_ice = state.moisture.ρq_ice * ρinv
end

function compute_gradient_flux!(
    moist::NonEquilMoist,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    # diffusive fluxes of moisture variables
    diffusive.moisture.∇q_tot = ∇transform.moisture.q_tot
    diffusive.moisture.∇q_liq = ∇transform.moisture.q_liq
    diffusive.moisture.∇q_ice = ∇transform.moisture.q_ice
end

function flux_moisture!(
    moist::NonEquilMoist,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    u = state.ρu / ρ
    flux.moisture.ρq_tot += u * state.moisture.ρq_tot
    flux.moisture.ρq_liq += u * state.moisture.ρq_liq
    flux.moisture.ρq_ice += u * state.moisture.ρq_ice
end

function flux_second_order!(
    moist::NonEquilMoist,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_t,
)
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    d_q_liq = (-D_t) .* diffusive.moisture.∇q_liq
    d_q_ice = (-D_t) .* diffusive.moisture.∇q_ice

    flux_second_order!(moist, flux, state, d_q_tot, d_q_liq, d_q_ice)
end

function flux_second_order!(
    moist::NonEquilMoist,
    flux::Grad,
    state::Vars,
    d_q_tot,
    d_q_liq,
    d_q_ice,
)
    flux.ρ += d_q_tot * state.ρ
    flux.ρu += d_q_tot .* state.ρu'
    flux.moisture.ρq_tot += d_q_tot * state.ρ
    flux.moisture.ρq_liq += d_q_liq * state.ρ
    flux.moisture.ρq_ice += d_q_ice * state.ρ
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module BioBenchmarks

export
    @benchmark,
    datafile

# Benchmark Runner
# ----------------

"""
    run(revision;
        force=false,
        stop_on_error=false,
        benchmark_directory=Pkg.dir("BioBenchmarks", "benchmarks"),
        benchmark_scripts=readdir(benchmark_directory))

Run benchmark scripts of `revision` and return the filepath of the results.

# Arguments
* `revision`: Git revision of Bio.jl (branch, tag, or hash).
* `force`: if `true`, run benchmarks even if the result file exists.
* `stop_one_error`: if `true`, immediately stop benchmarks when a script fails.
* `benchmark_directory`: directory of benchmark scripts.
* `benchmark_scripts`: benchmark script names in `benchmark_directory`.
"""
function run(revision;
             force=false,
             stop_on_error=false,
             benchmark_directory=Pkg.dir("BioBenchmarks", "benchmarks"),
             benchmark_scripts=readdir(benchmark_directory))
    suffix = cd(Pkg.dir("Bio")) do
        Base.run(`git checkout $revision`)
        gitrevision()
    end

    outputfile = "results_$(suffix).tsv"
    if isfile(outputfile) && !force
        info("use the existing benchmark results of $revision")
        return outputfile
    end

    out = open(outputfile, "w")
    print_header(out)
    flush(out)  # this seems to be necessary (but don't know why)
    for script in benchmark_scripts
        info("running $script")
        try
            Base.run(pipeline(
                `julia $(joinpath(benchmark_directory, script))`,
                out))
        catch ex
            if stop_on_error
                rethrow()
            else
                warn(ex)
            end
        end
    end
    close(out)

    return outputfile
end

function gitrevision()
    commit = chomp(readall(`git rev-parse HEAD`))
    branch = chomp(readall(`git rev-parse --abbrev-ref HEAD`))
    return string(branch, '_', commit[1:7])
end

macro benchmark(ex)
    quote
        results = ResultSet(splitext(basename(@__FILE__))[1])
        total_elapsedtime = UInt64(0)
        $(esc(ex))  # JIT compiling
        for _ in 1:5
            if total_elapsedtime > 60 * 10^9  # 60 sec
                break
            end
            gcstats = Base.gc_num()
            elapsedtime = Base.time_ns()
            $(esc(ex))
            elapsedtime = Base.time_ns() - elapsedtime
            total_elapsedtime += elapsedtime
            push!(results, Result(elapsedtime, Base.GC_Diff(Base.gc_num(), gcstats)))
        end
        BioBenchmarks.print_results(results)
    end
end

# Benchmark Result
# ----------------

immutable Result
    elapsed::UInt64
    diff::Base.GC_Diff

    Result(elapsed, diff) = new(elapsed, diff)
end

immutable ResultSet
    name::Symbol
    results::Vector{Result}

    ResultSet(name) = new(name, Result[])
end

function Base.push!(s::ResultSet, r::Result)
    push!(s.results, r)
    return s
end

function print_header(out)
    println(out, "name\ttrial\telapsed\tgc_time\tallocated")
end

function print_results(set::ResultSet)
    for (i, r) in enumerate(set.results)
        println(
            set.name,          '\t',
            i,                 '\t',
            r.elapsed,         '\t',
            r.diff.total_time, '\t',
            r.diff.allocd
        )
    end
end

# Utils for Data Files
# --------------------

function datafile(filename)
    datadir = Pkg.dir("BioBenchmarks", "data")
    if !isdir(datadir)
        mkdir(datadir)
    end

    filepath = joinpath(datadir, filename)
    if isfile(filepath)
        return filepath
    end

    if ismatch(r"^chr(\d{1,2}|[XYM])\.fa$", filename)
        # human chromosome
        recipe = pipeline(
            `curl http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/$filename.gz`,
            `gzip -cd`,
            filepath
        )
    elseif ismatch(r"^\w{6}\.fasta$", filename)
        # protein sequence of UniProt
        recipe = pipeline(
            `curl http://www.uniprot.org/uniprot/$filename`,
            filepath
        )
    else
        error("don't know how to get $filename")
    end

    Base.run(recipe)

    return filepath
end

end # module


text:
module RandomUseful

using LinearAlgebra

include("fourier_diff.jl")
include("cheb_diff.jl")
include("delaydifferentialequation.jl")

end # module


text:
const _js_loaded = [false]

function ijulia()

    # if we're in IJulia call setup the notebook js interop
    if isdefined(Main, :IJulia) && Main.IJulia.inited
        # borrowed from https://github.com/plotly/plotly.py/blob/2594076e29584ede2d09f2aa40a8a195b3f3fc66/plotly/offline/offline.py#L64-L71
        # and https://github.com/JuliaLang/Interact.jl/blob/cc5f4cfd34687000bc6bc70f0513eaded1a7c950/src/IJulia/setup.jl#L15
        if !_js_loaded[1]

            const _pages_js = readstring(joinpath(dirname(@__FILE__),"pages.js"))

            display("text/html", """

pages.js loaded.

            """)
            _js_loaded[1] = true
        end
    end

end


text:
module MixedBlockArrays

using Reexport

@reexport using BlockArrays
@reexport using FillArrays

import Base.Cartesian: @nloops, @ntuple
import BlockArrays: block, blockindex

export MixedBlockArray

struct MixedBlockArray{T, N, BS<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractBlockArray{T, N}
    blocks::AbstractArray{AbstractArray{T,N},N}
    axes::BS

    MixedBlockArray(blocks::R, blocked_axes::BS) where {T, N, R<:AbstractArray{AbstractArray{T,N},N}, BS<:NTuple{N,AbstractUnitRange{Int}}} =
        new{T, N, BS}(blocks, blocked_axes)
end

@generated function MixedBlockArray{T}(axes::NTuple{N,AbstractUnitRange{Int}}) where {T, N, R<:AbstractArray{AbstractArray{T,N},N}}
    return quote
        n_blocks = map(blocklength, axes)
        blocks = Array{AbstractArray{T,N},N}(undef, n_blocks)
        @nloops $N block dim->blockaxes(axes[dim],1) begin
            block_indices = @ntuple $N block
            indices = getindex.(axes, block_indices)
            block_size = length.(indices)
            blocks[Int.(block_indices)...] = Zeros{T}(block_size...)
        end
        return MixedBlockArray(blocks, axes)
    end
end

BlockArrays.getindex(A::MixedBlockArray, b1::Block{1}, b2::Block{1}, btail::Vararg{Block{1}}) = getindex(A, Block(b1.n..., b2.n...), btail...)
BlockArrays.getindex(A::MixedBlockArray, I::Block) = A.blocks[I.n...]

function BlockArrays.setindex!(A::MixedBlockArray, v,b1::Block{1}, b2::Block{1},
    btail::Vararg{Block{1}})

    setindex!(A, v, Block(b1.n..., b2.n...), btail...)
end
BlockArrays.setindex!(A::MixedBlockArray, v, I::Block) = setindex!(A.blocks, v, I.n...)

Base.axes(A::MixedBlockArray) = A.axes
Base.size(A::MixedBlockArray) = length.(axes(A))
function Base.getindex(A::MixedBlockArray, I::Vararg{Int})
    Bi = findblockindex.(axes(A), I)
    B, i = block.(Bi), blockindex.(Bi)
    bl = A[B...]
    return bl[i...]
end

function Base.setindex!(A::MixedBlockArray, v, I::Vararg{Int})
    Bi = findblockindex.(axes(A), I)
    B, i = block.(Bi), blockindex.(Bi)
    bl = A[B...]
    setindex!(bl, v, i...)
end

end


text:
Base.cconvert(::Type{Ptr{TSTreeCursor}}, cursor::TSTreeCursor) =
    Base.cconvert(Ptr{TSTreeCursor}, Ref(cursor))

function cursor_new(node::TSNode)
    ts_tree_cursor_new(node)
end

function cursor_current_node(cursor::TSTreeCursor)
    ts_tree_cursor_current_node(cursor)
end

function cursor_current_field_name(cursor::TSTreeCursor)
    s = ts_tree_cursor_current_field_name(cursor)
    s == C_NULL && return nothing
    return unsafe_string(s)
end

function cursor_current_field_id(cursor::TSTreeCursor)
    ts_tree_cursor_current_field_id(cursor)
end

function cursor_goto_parent!(cursor::TSTreeCursor)
    parent_exists = ts_tree_cursor_goto_parent(cursor)
    parent_exists && return cursor
end

function cursor_goto_next!(cursor::TSTreeCursor)
    has_sibling = ts_tree_cursor_goto_next_sibling(cursor)
    has_sibling && return cursor
end

function cursor_goto_child!(cursor::TSTreeCursor)
    has_child = ts_tree_cursor_goto_first_child(cursor)
    has_child && return cursor
end

function cursor_copy(cursor::TSTreeCursor)
    new_cursor = ts_tree_cursor_copy(cursor)
    return new_cursor
end

function cursor_reset(cursor::TSTreeCursor, node::TSNode)
    ncursor::TSTreeCursor = ts_tree_cursor_reset(cursor, node)
    return ncursor
end

function node_type(cursor::TSTreeCursor)
    node_type(cursor_current_node(cursor))
end


text:
run(`julia configure_pycall.jl`)

using PyCall

const PIP = joinpath(split(PyCall.PYTHONHOME, ":")[end], "bin/pip")

function pipinstall(pkg)
    try
        pyimport(pkg)
    catch
        println("`$pkg` is not available in the current PyCall.jl.")
        println("Installing using pip ...")
        run(`$PIP install $pkg`)
    end
end

pipinstall("matplotlib")
pipinstall("wandb")


text:
function simulate!(s::Simulator{T};
	clock_time_noise=1.0e-3,
	verbose=false) where T
    status = false

    p = s.policy
	w = s.dist
	traj = s.traj

	N = length(traj.u)

	# state
	t = 1
	x = [traj.q[t+1]; (traj.q[t+1] - traj.q[t]) ./ s.h; traj.γ[t]]

	# last_computed_control = p.ref_traj.u[1] / N_sample
	clock_time = 0.0
    for t = 1:N
		# println("timestep $t / $N, time $time")
		policy_time = @elapsed traj.u[t] = exec_policy(p, x, clock_time)
		traj.u[t] .*= s.h # transform force from policy to impulse for RoboDojo simulator
		s.opts.record && (s.stats.policy_time[t] = policy_time)

        # disturbances
        traj.w[t] .= disturbances(w, traj.q[t+1], t)

        # step
        status = RoboDojo.step!(s, t, verbose=verbose)
		x .= [traj.q[t+2]; (traj.q[t+2] - traj.q[t+1]) ./ s.h; traj.γ[t]]
		clock_time += s.h + clock_time_noise * rand(1)[1]

        !status && break
    end

    return status
end

function exec_policy(p::CIMPC{T,NQ,NU,NW,NC}, x::Vector{T}, t::T) where {T,NQ,NU,NW,NC}
	# reset
	if t ≈ 0.0
		p.next_time_update = 0.0
		p.q0 .= p.ref_traj.q[1]
		p.altitude .= 0.0
		p.buffer_time = 0.0
		set_trajectory!(p.traj, p.ref_traj)
		set_implicit_trajectory!(p.im_traj, p.im_traj_cache)
		reset_window!(p.window)
	end

    if t >= p.next_time_update
		(p.opts.altitude_update && t > 0.0) && (update_altitude!(p.altitude, p.ϕ, p.s,
									x, NQ, NC, p.N_sample,
									threshold = p.opts.altitude_impact_threshold,
									verbose = p.opts.altitude_verbose))

		set_altitude!(p.im_traj, p.altitude)
		update!(p.im_traj, p.traj, p.s, p.altitude, p.κ[1], p.traj.H)

		# visualize
		# p.opts.live_plotting && live_plotting(p.s.model, p.traj, traj, p.newton, p.q0, traj.q[t+1], t)

		# shift trajectory
		rot_n_stride!(p.traj, p.traj_cache, p.stride, p.window)

		# update
		update_window!(p.window, p.ref_traj.H)

		# reset update time
		p.next_time_update = (t - t % p.traj.h) + p.traj.h
    end

	policy_time = @elapsed begin
		if p.buffer_time <= 0.0
			# @show t
			# optimize
			q1 = x[1:NQ]
			q0 = x[1:NQ] - x[NQ .+ (1:NQ)] .* p.traj.h
			newton_solve!(p.newton, p.s, q0, q1,
				p.window, p.im_traj, p.traj, warm_start = t > 0.0)

			update!(p.im_traj, p.traj, p.s, p.altitude, p.κ[1], p.traj.H)
		end
	end

	# update buffer time
	(p.buffer_time <= 0.0) && (p.buffer_time = policy_time)
	p.buffer_time -= p.traj.h

	# scale control
	p.u .= p.newton.traj.u[1]
	p.u ./= p.traj.h

	return p.u
end


text:
module GeneralizedStateSpace

include("utils.jl")
include("state_space_function.jl")

export StateSpaceFunction

end


text:
# -------- Async Wrapper ----------------------------------------------------- #

"""
Asynchronous model wrapper that allows a model to be called on a batch of games
in parallel when the single calls take place in an async context. Note that
an Async model always returns CPU arrays, even if the model it is based on works
on the GPU.
"""
mutable struct Async{G} <: Model{G, false}
  model          :: Model{G}
  channel        :: Channel
  thread

  max_batchsize  :: Int
  buffersize     :: Int     # Must not be smaller than max_batchsize!
end

"""
    Async(model; max_batchsize, buffersize)

Wraps `model` to an asynchronous model with maximal batchsize `max_batchsize`
for evaluation in parallel and a buffer of size `buffersize` for queuing.
"""
function Async( model :: Model{G};
                max_batchsize = 50,
                buffersize = 2*max_batchsize ) where {G <: Game}

  # Make sure that the buffer is larger than the maximal allowed batchsize
  @assert buffersize >= max_batchsize "buffersize must be larger than max_batchsize"

  # Open the channel in which the inputs are dumped
  channel = Channel(buffersize)

  # Start the worker thread in the background
  thread = @async worker_thread(channel, model, max_batchsize)

  # Create the instance
  amodel = Async{G}(model, channel, thread, max_batchsize, buffersize)

  # Register finalizer
  finalizer(m -> close(m.channel), amodel)

  # Return the instance
  amodel

end

function (m :: Async{G})( game :: G
                        , use_features = false
                        ) where {G <: Game}

  @assert !use_features "Features cannot be used in Async."

  out_channel = Channel(1)
  bind(out_channel, m.thread)
  put!(m.channel, (game, out_channel))
  take!(out_channel)

end

function (m :: Async{G})( games :: Vector{G}
                        , use_features = false
                        ) where {G <: Game}

  @assert !use_features "Features should never be used in Async. "
  @warn "Calling Async model in batched mode is not recommended." maxlog=1

  outputs = asyncmap(x -> m(x, use_features), games, ntasks = m.buffersize)
  cat_outputs(outputs)

end

function switch_model(m :: Async{G}, model :: Model{G}) where {G <: Game}
  Async( model
       , max_batchsize = m.max_batchsize
       , buffersize = m.buffersize )
end

swap(m :: Async) = @warn "Async cannot be swapped. Remaining on CPU."
Base.copy(m :: Async) = switch_model(m, copy(m.model))

ntasks(m :: Async) = m.buffersize
base_model(m :: Async) = m.model
training_model(m :: Async) = m.model
worker_model(m :: Async) = m.model

# Async networks cannot calculate features. To get the features of the
# network on which it is based, access them via training_model(...)
features(m :: Async) = Feature[]

worker_model_to_cpu(m :: Async) = switch_model(m, to_cpu(m.model))
worker_model_to_gpu(m :: Async) = switch_model(m, to_gpu(m.model))

# -------- Async Worker ------------------------------------------------------ #

closed_and_empty(channel) = try fetch(channel); false catch _ true end

function worker_thread(channel, model, max_batchsize)

  try

    while !closed_and_empty(channel)

      inputs = Vector()
      # If we arrive here, there is at least one thing to be done.
      while isready(channel) && length(inputs) < max_batchsize
        push!(inputs, take!(channel))
        yield()
      end

      v, p, f = model(first.(inputs))
      v = v |> to_cpu
      p = p |> to_cpu
      f = f |> to_cpu

      for i = 1:length(inputs)
        put!(inputs[i][2], (v[i], p[:,i], f[:,i]))
      end

    end

  catch err

    error("Worker died: $err")

  end

end

max_stars_repo_name: roSievers/jtac


text:
module  SummaryInfo

using PopGen
using DataFrames
using Test

cats = @nancycats;

@testset "allele_avg" begin
    @test length(allele_avg(cats,  rounding = true)) == 2
    @test typeof(allele_avg(cats,  rounding = true)) <: NamedTuple
    @test length(allele_avg(cats,  rounding = false)) == 2
    @test typeof(allele_avg(cats,  rounding = false)) <: NamedTuple
end

@testset "richness" begin
    @test typeof(richness(cats, by = "locus")) == DataFrame
    @test size(richness(cats, by = "locus")) == (9, 2)
    @test typeof(richness(cats, by = "population")) == DataFrame
    @test size(richness(cats, by = "population")) == (153, 3)
end

@testset "summary F/D/etc. stats" begin
    smry_glob = summary(cats)
    smry_loc = summary(cats, by = "locus")
    @test typeof(smry_glob) == DataFrame
    @test typeof(smry_loc) == DataFrame
    @test size(smry_glob) == (1, 10)
    @test size(smry_loc) ==(9, 11)
end

end # module


text:
#=

@authors : Spencer Lyon <spencer.lyon@nyu.edu>, John Stachurski

=#

using QuantEcon
using Distributions

"""
The Lucas asset pricing model --- parameters and grid data
"""
struct LucasTree{TF<:AbstractFloat}
    γ::TF                 # coefficient of risk aversion
    β::TF                 # Discount factor in (0, 1)
    α::TF                 # Correlation coefficient in the shock process
    σ::TF                 # Volatility of shock process
    ϕ::Distribution       # Distribution for shock process
    grid::Vector{TF}      # Grid of points on which to evaluate prices
    shocks::Vector{TF}    # Draws of the shock
    h::Vector{TF}         # The h function represented as a vector
end

"""
Constructor for the Lucas asset pricing model
"""
function LucasTree(;γ::AbstractFloat=2.0,
                β::AbstractFloat=0.95,
                α::AbstractFloat=0.9,
                σ::AbstractFloat=0.1,
                grid_size::Integer=100)

    ϕ = LogNormal(0.0, σ)
    shocks = rand(ϕ, 500)

    # == build a grid with mass around stationary distribution == #
    ssd = σ / sqrt(1 - α^2)
    grid_min, grid_max = exp(-4 * ssd), exp(4 * ssd)
    grid = collect(linspace(grid_min, grid_max, grid_size))

    # == set h(y) = β * int u'(G(y,z)) G(y,z) ϕ(dz) == #
    h = similar(grid)
    for (i, y) in enumerate(grid)
        h[i] = β * mean((y^α .* shocks).^(1 - γ))
    end

    return LucasTree(γ,
                     β,
                     α,
                     σ,
                     ϕ,
                     grid,
                     shocks,
                     h)
end

"""
The approximate Lucas operator, which computes and returns updated function
Tf on the grid points.
"""
function lucas_operator(lt::LucasTree, f::Vector)

    # == unpack names == #
    grid, α, β, h = lt.grid, lt.α, lt.β, lt.h
    z = lt.shocks

    Af = LinInterp(grid, f)

    Tf = [h[i] + β * mean(Af.(grid[i]^α.*z)) for i in 1:length(grid)]
    return Tf
end

"""
Compute the equilibrium price function associated with Lucas tree `lt`
"""
function solve_lucas_model(lt::LucasTree;
                           tol::AbstractFloat=1e-6,
                           max_iter::Integer=500)

    # == simplify notation == #
    grid, γ = lt.grid, lt.γ

    i = 0
    f = zeros(grid)  # Initial guess of f
    error = tol + 1

    while (error > tol) && (i < max_iter)
        f_new = lucas_operator(lt, f)
        error = maximum(abs, f_new - f)
        f = f_new
        i += 1
    end

    # p(y) = f(y) * y ^ γ
    price = f .* grid.^γ

    return price
end


text:
#
# Parameter explanations
# MvNormal: Distribution of phenotypes is multivariate normal
# q: number of cross validation folds
# init_beta: get good initial estimates for genetic and nongenetic regression coefficients
# min_iter: iterate at least 10 times before checking for convergence
#
using MendelIHT, Random
Random.seed!(2021)                        # seed for reproducibility
plinkfile = "ukb.plink.filtered"          # plink files without .bed/bim/fam
phenotypes = "normalized.bmi.sbp.dbp.txt" # comma-seprated phenotype file (no header line)
covariates = "iht.covariates.txt"         # comma-separated covariates file (no header line)

# cross validate k = 100, 200, ..., 1000
path = 100:100:1000
@time mses = cross_validate(plinkfile, MvNormal, path=path, q=3,
    covariates=covariates, phenotypes=phenotypes, init_beta=true, min_iter=10,
    cv_summaryfile="cviht.summary.roughpath1.txt")

# if best k was 200 in previous step, now cross validate k = 110, 120, ..., 290
k_rough_guess = path[argmin(mses)]
path = (k_rough_guess - 90):10:(k_rough_guess + 90)
@time mses = cross_validate(plinkfile, MvNormal, path=path, q=3,
    covariates=covariates, phenotypes=phenotypes, init_beta=true, min_iter=10,
    cv_summaryfile="cviht.summary.roughpath2.txt")

# if best k was 190 in previous step, now cross validate k = 181, 182, ..., 209
k_rough_guess = path[argmin(mses)]
dense_path = (k_rough_guess - 9):(k_rough_guess + 9)
@time mses_new = cross_validate(plinkfile, MvNormal, path=dense_path,
    covariates=covariates, q=3, phenotypes=phenotypes, min_iter=10, init_beta=true)

# now run IHT on best k
@time iht_result = iht(plinkfile, dense_path[argmin(mses_new)], MvNormal,
    covariates=covariates, phenotypes=phenotypes, min_iter=10, init_beta=true)


text:
# const epoch_dt_format2 = Dates.DateFormat("yy mm dd HH MM SS.sssssss")
function parse_2epoch_time(ephs::String, toUTC::Millisecond)

    dtstr = SubString(ephs, 2, 26)
    y = parse(Int64, dtstr[1:2])
    if y < 70 y += 2000
    else y += 1900
    end

    m = parse_Int_or_Empty(dtstr[4:5])
    d = parse_Int_or_Empty(dtstr[7:8])
    H = parse_Int_or_Empty(dtstr[10:11])
    M = parse_Int_or_Empty(dtstr[13:14])
    S = parse_Int_or_Empty(dtstr[16:17])
    # s = parse(Int64, dtstr[19:25])
    # TODO only Millisecond three digit maybe later use TimeDates Nanosecond.
    s = parse_Int_or_Empty(dtstr[19:21])
    this_dt = Dates.DateTime(y, m, d, H, M, S, s)
    # @show "  $ephs"
    # @show "$this_dt"
    this_dt += toUTC
    this_dt
end

function add_epoch_line2!(sods::Dict{SatID,SatelliteObsData},
                f::IOStream, satid::SatID, this_dt::DateTime,
                systemNumObs::D_SystemNumObs,
                recommend_length::Int64, )

    obsNum, obsIds = get(systemNumObs, satid.ss, (0, Array{ObsID,1}()))
    if obsNum == 0
        ss = collect(keys(systemNumObs))[1]
        obsNum, _ = systemNumObs[ss]
        nl = ceil(Int64, obsNum / 5.0)
        # @show obsNum, nl
        for i in 1:nl
            # @show i,
            readline(f)
        end
        return
        throw(error("no such system $satid"))
    end

    if haskey(sods, satid)
        this_sod = sods[satid]
        this_dts = this_sod.dts
        this_obs = this_sod.obs
    else
        this_dts  = Array{DateTime,1}()
        sizehint!(this_dts, recommend_length)
        this_obs = Array{Array{RinexDatum,1},1}(undef, obsNum)
        for ion in 1:obsNum
            this_obs[ion] = Array{RinexDatum,1}()
            sizehint!(this_obs[ion], recommend_length)
        end
        sods[satid] = SatelliteObsData(satid, obsIds,
                                        this_dts, this_obs)
    end
    si = 81
    line = ""
    for ion in 1:obsNum
        # @show ion, si
        if si > 80
            line = rpad(readline(f), 80)
            si = 1
        end
        push!(this_obs[ion],
                RinexDatum(SubString(line, si, si + 15)), )
        si += 16
    end
    push!(this_dts, this_dt)
    return
end

function add_this_epoch_comments2!(comments::Dict{DateTime,String},
                dt_before::DateTime, f::IOStream, numSVs::Int64, )
    this_comments = Array{String,1}(undef, numSVs)
    for ii_el in 1:numSVs
        this_comments[ii_el]  = readline(f)
    end
    if haskey(comments, dt_before)
        comments[dt_before] = string(comments[dt_before],"\n",
                                join(this_comments, "\n"))
    else
        comments[dt_before] = join(this_comments, "\n")
    end
end

function read_one_epoch2(f::IOStream,
    epochs::Dict{DateTime,EpochHeader},
    comments::Dict{DateTime,String},
    add_this_epoch_line2::Function,
    toUTC::Millisecond,
    ephs_dr_set::Set{Int64},
    dt_before::DateTime,)::DateTime
    ephs = readline(f)

    if ephs[1] != ' ' || ephs[4] != ' ' || ephs[7] != ' '
        throw(error("epoch line seems wrong! \n|$ephs|"))
    end

    epf = parse(Int64, ephs[29])

    numSVs = parse(Int64, ephs[30:32])

    if numSVs == 0
        @warn "0 satellite? \n|$ephs|"
    end

    if epf in ephs_dr_set # 0,1,6 data record
        this_dt = parse_2epoch_time(ephs, toUTC)

        if haskey(epochs, this_dt)
            @warn "duplicate epoch? $ephs"
            delete!(epochs, this_dt)
        end
        len_ephs = length(ephs)
        if len_ephs > 68
            index_end = len_ephs > 80 ? 80 : len_ephs
            rco = parse(Float64, ephs[69:index_end])
        else
            rco = 0.0
        end

        epochs[this_dt] = EpochHeader(this_dt, epf, numSVs, rco)

        # ephs_arr = Array{String,1}()

        satIndex = Array{SatID,1}(undef, numSVs)
        start_index = 30
        for ndx in 1:numSVs
            if start_index > 65
                ephs = readline(f)
                start_index = 30
            end
            start_index += 3
            satIndex[ndx] = SatID(ephs[start_index:start_index + 2])
        end
        for ii_el in 1:numSVs
            add_this_epoch_line2(f, satIndex[ii_el], this_dt)
        end
    else
        add_this_epoch_comments2!(comments, dt_before, f, numSVs)
        this_dt = dt_before
    end
    this_dt
end

function read_2obsdata(f::IOStream, sysNumObs::D_SystemNumObs,
    toUTC::Millisecond=Millisecond(0),
    recommend_length::Int64=4800,
    ephs_dr_set::Set{Int64}=Set([0, 1, 6]))
    # rod = RinexObsData()
    # @unpack_RinexObsData rod
    # epochs, comments, obs, flag, dt

    epochs  = Dict{DateTime,EpochHeader}()
    comments = Dict{DateTime,String}()

    sods = Dict{SatID,SatelliteObsData}()

    add_this_epoch_line2(f::IOStream,
        satid::SatID, this_dt::DateTime) =
                    add_epoch_line2!(sods, f, satid, this_dt,
                                sysNumObs, recommend_length, )

    dt_before = DateTime(0, 1, 1, 0, 0, 0)
    while !eof(f)
        dt_before = read_one_epoch2(f,
                epochs, comments,
                add_this_epoch_line2,
                toUTC, ephs_dr_set, dt_before)
    end
    # sys = Dict{SatSys,Array{SatID,1}}()
    # # , sys length of each sat
    # init_sys!(sys, dts)
    # RinexObsData(epochs, comments, obs, flags, dts, sys)

    # @bp
    # @show epochs
    # @show comments
    # @show sods
    RinexObsData(epochs, comments, sods)
end


text:
"Check the number of clusters in the sample."
nclusters(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GaussianCopulaVCModel, Poisson_Bernoulli_VCModel, GLMCopulaARModel, NBCopulaARModel, GaussianCopulaARModel, GLMCopulaCSModel, NBCopulaCSModel, GaussianCopulaCSModel}) = length(model.data)

"Check minimum cluster size."
minclustersize(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GaussianCopulaVCModel, GLMCopulaARModel, NBCopulaARModel, GaussianCopulaARModel, GLMCopulaCSModel, NBCopulaCSModel, GaussianCopulaCSModel}) = minimum(model.data[i].n for i in 1:nclusters(model))

"Check maximum cluster size."
maxclustersize(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GaussianCopulaVCModel, GLMCopulaARModel, NBCopulaARModel, GaussianCopulaARModel, GLMCopulaCSModel, NBCopulaCSModel, GaussianCopulaCSModel}) = maximum(model.data[i].n for i in 1:nclusters(model))

"Check the number of variance components."
nvc(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GaussianCopulaVCModel, Poisson_Bernoulli_VCModel}) = model.m

"Check the number of fixed effects to estimate."
nfixedeffects(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GaussianCopulaVCModel, Poisson_Bernoulli_VCModel, GLMCopulaARModel, NBCopulaARModel, GaussianCopulaARModel, GLMCopulaCSModel, NBCopulaCSModel, GaussianCopulaCSModel}) = model.p

"Check the base distribution."
basedistribution(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GLMCopulaARModel, NBCopulaARModel, GLMCopulaCSModel, NBCopulaCSModel}) = Base.typename(typeof(unique(model.d)[1])).wrapper

"Check the link function."
linkfunction(model::Union{GLMCopulaVCModel, NBCopulaVCModel, GLMCopulaARModel, NBCopulaARModel, GLMCopulaCSModel, NBCopulaCSModel}) = Base.typename(typeof(unique(model.link)[1])).wrapper

# gaussian base
"Check the base distribution."
basedistribution(model::Union{GaussianCopulaVCModel, GaussianCopulaARModel, GaussianCopulaCSModel}) = Base.typename(typeof(Normal())).wrapper

"Check the link function."
linkfunction(model::Union{GaussianCopulaVCModel, GaussianCopulaARModel, GaussianCopulaCSModel}) = Base.typename(typeof(IdentityLink())).wrapper

# better printing; customize how a type is summarized in a REPL
function Base.show(io::IO, x::Union{GLMCopulaVCModel, NBCopulaVCModel, GaussianCopulaVCModel})
    print(io, "Quasi-Copula Variance Component Model\n")
    print(io, "  * base distribution: $(basedistribution(x))\n")
    print(io, "  * link function: $(linkfunction(x))\n")
    print(io, "  * number of clusters: $(nclusters(x))\n")
    print(io, "  * cluster size min, max: $(minclustersize(x)), $(maxclustersize(x))\n")
	print(io, "  * number of variance components: $(nvc(x))\n")
    print(io, "  * number of fixed effects: $(nfixedeffects(x))\n")
    if x.penalized == true
        print(io, "  * L2 ridge penalty on variance components: $(x.penalized)")
    end
end

function Base.show(io::IO, x::Union{GLMCopulaARModel, NBCopulaARModel, GaussianCopulaARModel})
    print(io, "Quasi-Copula Autoregressive AR(1) Model\n")
    print(io, "  * base distribution: $(basedistribution(x))\n")
    print(io, "  * link function: $(linkfunction(x))\n")
    print(io, "  * number of clusters: $(nclusters(x))\n")
    print(io, "  * cluster size min, max: $(minclustersize(x)), $(maxclustersize(x))\n")
    print(io, "  * number of fixed effects: $(nfixedeffects(x))\n")
    if x.penalized == true
        print(io, "  * L2 ridge penalty on AR(1) variance parameter: $(x.penalized)")
    end
end

function Base.show(io::IO, x::Union{GLMCopulaCSModel, NBCopulaCSModel, GaussianCopulaCSModel})
    print(io, "Quasi-Copula Compound Symmetric CS Model\n")
    print(io, "  * base distribution: $(basedistribution(x))\n")
    print(io, "  * link function: $(linkfunction(x))\n")
    print(io, "  * number of clusters: $(nclusters(x))\n")
    print(io, "  * cluster size min, max: $(minclustersize(x)), $(maxclustersize(x))\n")
    print(io, "  * number of fixed effects: $(nfixedeffects(x))\n")
    if x.penalized == true
        print(io, "  * L2 ridge penalty on CS variance parameter: $(x.penalized)")
    end
end

# bivariate mixed vc
"Check cluster size."
clustersize(model::Poisson_Bernoulli_VCModel) = unique(model.data[i].n for i in 1:nclusters(model))[1]

"Check all the base distributions."
basedistributions(model::Poisson_Bernoulli_VCModel) = [Base.typename.(typeof.(unique(model.vecd)[1][i])).wrapper for i in 1:length(unique(model.vecd[1]))]

"Check all the link functions."
linkfunctions(model::Poisson_Bernoulli_VCModel) = [Base.typename.(typeof.(unique(model.veclink)[1][i])).wrapper for i in 1:length(unique(model.veclink[1]))]

# better printing; customize how a type is summarized in a REPL
function Base.show(io::IO, x::Poisson_Bernoulli_VCModel)
    print(io, "Bivariate Mixed Quasi-Copula Variance Component Model\n")
    print(io, "  * base Distributions for (y_1, y_2): ($(basedistributions(x)[1]), $(basedistributions(x)[2]))\n")
    print(io, "  * link functions for (y_1, y_2): ($(linkfunctions(x)[1]), $(linkfunctions(x)[2]))\n")
    print(io, "  * number of clusters: $(nclusters(x))\n")
    print(io, "  * cluster size: $(clustersize(x))\n")
	print(io, "  * number of variance components: $(nvc(x))\n")
    print(io, "  * number of fixed effects: $(nfixedeffects(x))")
end


text:
Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y

instantiate(x::ModelingToolkit.Variable{Real}) = 1.0
instantiate(x::ModelingToolkit.Variable) = oneunit(1*ModelingToolkit.vartype(x))
function instantiate(x::Num)
    x = value(x)
    if operation(x) isa Sym
        return instantiate(operation(x))
    elseif operation(x) isa Differential
        instantiate(arguments(x)[1])/instantiate(arguments(x)[1].args[1])
    else
        operation(x)(instantiate.(arguments(x))...)
    end
end

function validate(eq::ModelingToolkit.Equation)
    try
        return typeof(instantiate(eq.lhs)) == typeof(instantiate(eq.rhs))
    catch
        return false
    end
end

function validate(sys::AbstractODESystem)
    all(validate.(equations(sys)))
end


text:
Base.:(==)(x::A, y::B) where {A<:Symbolic,B<:Symbolic} = false
Base.promote(x::A, y::B) where {A<:Symbolic,B<:Symbolic} =
    (SymExpr(:identity, [x]), SymExpr(:identity, [y]))

macro new_number(type_name, comp_params, compare)
    Base_isequal = Expr(:call, Expr(Symbol("."), :Base, :(:(==))))
    signature,args = if comp_params.head == :where
        Expr(:where, Base_isequal, comp_params.args[2:end]...),comp_params.args[1].args
    else
        Base_isequal,comp_params.args
    end
    append!(Base_isequal.args, args)
    def = Expr(:(=), signature, compare)

    quote
        Base.:(==)(x::$(esc(type_name)), y::N) where {N<:Union{Real,Complex}} = false
        Base.:(==)(x::N, y::$(esc(type_name))) where {N<:Union{Real,Complex}} = false
        Base.:(==)(x::$(esc(type_name)), y::Sym) = false
        Base.:(==)(x::Sym, y::$(esc(type_name))) = false
        Base.:(==)(x::$(esc(type_name)), y::SymExpr) = false
        Base.:(==)(x::SymExpr, y::$(esc(type_name))) = false

        Base.promote(::Type{<:$(esc(type_name))}) = SymExpr

        Base.promote(x::TT, y::SymExpr) where {TT<:$(esc(type_name))} =
            (SymExpr(:identity, [x]), y)
        Base.promote(x::SymExpr, y::TT) where {TT<:$(esc(type_name))} =
            (x, SymExpr(:identity, [y]))

        $(esc(def))
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PCL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PCL")
JLLWrappers.@generate_main_file("PCL", UUID("19da7e43-5f78-5205-a494-8c2c52b4cca1"))
end  # module PCL_jll


text:
function read_mesh(name)
    """
        reading from a OFF file in 3 dimensions, returning X0 (coordinates) and F (faces)
    """
    file = open(name)

    #check type of file
    file_type = strip(readline(file))
    if file_type != "OFF"
        throw(ArgumentError("Wrong type of file, only reads OFF files"))
    end

    #number of vertices/faces/edges:
    (n_verts, n_faces, n_edges) = tuple([parse(Int,s) for s in split(strip(readline(file)),' ')]...)

    #vertices
    X0 = zeros(3,n_verts)
    for i in 1:n_verts
        X0[:,i]= [parse(Float64,s) for s in split(strip(readline(file)),' ')]
    end

    #faces
    F = zeros(Int,3,n_faces)
    for i in 1:n_faces
        F[:,i] = [parse(Int,s) for s in split(strip(readline(file)),' ')[2:end]]
    end

    return X0, F+1
end


text:
lpsuffix(per::Day)= "_LP$(per.value)d"
finiteor0(x) = isfinite(x) ? x : 0.0
finiteor0(::Missing) = 0.0
function lpfilter(df;
    focalcols::Vector{Symbol},
    lpfilterperiods::Vector{<:Union{Day, Nothing}} = PARAM[:analysislpfilterperiods],
    buf=PARAM[:analysislpfilterbuffer],
    op=:mean
  )

  @assert issorted(df, :date)
  #these are the columns we will filter
  #originalgrowthcols::Vector{Symbol} = growthcols(df)
  @assert allunique(lpfilterperiods)

  #used for the start heuristic
  averageinterval = mean((df.date[2:end] .- df.date[1:(end-1)]) .|> (d)->d.value)
  for per ∈ lpfilterperiods
    (per === nothing) && continue #indicates the unfiltered columns

    #create the new cols
    lpcols = (s->Symbol(s,lpsuffix(per))).(focalcols)
    for lpcol ∈ lpcols
      df[!, lpcol] = missings(Float64, nrow(df))
    end

    startrow = (per.value / averageinterval) |> round |> Int
    for (i,r) ∈ enumerate(eachrow(df))
      #we don't want to start the series too early, but allow for some tolerance via buf
      #note buf needs to calibrated based on 1) the size of each period and 2) noting
      #that the dates are at the end of each period
      (i < startrow) && continue
      window = view(df,
        (df.date .≥ (r.date - per + buf)) .& (df.date .≤ r.date), focalcols)

      if op === :mean
        for (Fold, Fnew) ∈ zip(focalcols, lpcols)
          r[Fnew] = mean(skipmissing(window[!, Fold])) |> finiteormissing
        end

        @assert (finiteor0.(Vector(r[lpcols])) .≈
          finiteor0.(vec(sum(x->coalesce(x,0.0), Matrix(window),dims=1)) ./
          vec(sum(x->!ismissing(x), Matrix(window), dims=1)))) |> all
      elseif op === :sum
        for (Fold, Fnew) ∈ zip(focalcols, lpcols)
          r[Fnew] = sum(skipmissing(window[!, Fold])) |> finiteormissing
        end

        @assert (finiteor0.(Vector(r[lpcols])) .≈
          finiteor0.(vec(sum(x->coalesce(x,0.0), Matrix(window),dims=1)))) |> all
      else
        throw("unrecognized op $op")
      end
    end

    #adjust for growth rows having a missing first row
    for (Fold, Fnew) ∈ zip(focalcols, lpcols)
      if df[1,Fold] === missing #heuristic- if this is true, assume a growth row
        df[startrow, Fnew] = missing
      end
    end

  end

  df |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:analysismeasurefilename])_testlp.csv")
  return df
end

#helper function to match smoothing levels bewteen measures and benchmarks
function matchcolsonlpsuffix(col, allcols)::Vector{Symbol}
  local strfocalcols
  strallcols = string.(allcols)
  if occursin(r"_LP[0-9]+d$", string(col))
    suffix = match(r"LP[0-9]+d$", string(col)).match
    strfocalcols = allcols[(s->occursin(Regex("$(suffix)\$"),s)).(strallcols)]
  else #no smoothing case
    strfocalcols = setdiff(allcols, allcols[(s->occursin(r"LP[0-9]+d$",s)).(strallcols)])
  end

  return Symbol.(strfocalcols)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# To make composition possible, we implement [`compose`](#), which
# defaults to returning a [`Sequence`](#).

"""
    Sequence(transforms...)

`Transform` that applies multiple `transformations`
after each other.

You should not use this explicitly. Instead use [`compose`](#).
"""
struct Sequence{T<:Tuple where N} <: Transform
    transforms::T
end

Sequence(tfms...) = Sequence{typeof(tfms)}(tfms)

getrandstate(seq::Sequence) = getrandstate.(seq.transforms)

compose(tfm1::Transform, tfm2::Transform) = Sequence(tfm1, tfm2)
compose(seq::Sequence, tfm::Transform) = Sequence(seq.transforms..., tfm)
compose(seq::Sequence, ::Identity) = seq

function apply(seq::Sequence, items::Tuple; randstate = getrandstate(seq))
    for (tfm, r) in zip(seq.transforms, randstate)
        items = apply(tfm, items; randstate = r)
    end
    return items
end

apply(seq::Sequence, item::Item; randstate = getrandstate(seq)) =
    apply(seq, (item,); randstate = randstate) |> only

makebuffer(pipeline::Sequence, item::Item) = only.(makebuffer(pipeline, (item,)))
function makebuffer(pipeline::Sequence, items::Tuple)
    buffers = []
    for tfm in pipeline.transforms
        push!(buffers, makebuffer(tfm, items))
        items = apply(tfm, items)
    end
    return Tuple(buffers)
end

function apply!(buffers::Tuple, pipeline::Sequence, items::Tuple; randstate = getrandstate(pipeline))
    @assert length(buffers) == length(pipeline.transforms)
    @assert length(buffers[1]) == length(items)
    for (tfm, buffer, r) in zip(pipeline.transforms, buffers, randstate)
        items = apply!(buffer, tfm, items; randstate = r)
    end
    return items
end

function apply!(buffers, pipeline::Sequence, item::AbstractItem; randstate = getrandstate(pipeline))
    @assert length(buffers) == length(pipeline.transforms)
    for (tfm, buffer, r) in zip(pipeline.transforms, buffers, randstate)
        item = apply!(buffer, tfm, item; randstate = r)
    end
    return item
end


text:
using SchumakerSpline
using Base.Test

# Run tests

tic()
println("Second Derivative Test")
@time @test include("SecondDerivativeTest.jl")
toc()


text:
using Base.Test
using JuMP
using MIPVerify

@testset "sequential.jl" begin

    @testset "Sequential" begin
        nnparams = Sequential([
            Conv2d(rand(1, 4, 4, 16), rand(16), 2),
            ReLU(),
            Flatten([3, 4, 1, 2]),
            Linear(rand(1000, 50), rand(50)),
            MaskedReLU(rand(MersenneTwister(0), [-1, 0, 1], 50)),
            Linear(rand(50, 10), rand(10))
        ], "testnet")
        io = IOBuffer()
        Base.show(io, nnparams)
        @test String(take!(io)) == """
        sequential net testnet
          (1) Conv2d(4, 16, kernel_size=(1, 4), stride=(2, 2), padding=same)
          (2) ReLU()
          (3) Flatten(): flattens 4 dimensional input, with dimensions permuted according to the order [3, 4, 1, 2]
          (4) Linear(1000 -> 50)
          (5) MaskedReLU with expected input size (50,). (18 zeroed, 17 as-is, 15 rectified).
          (6) Linear(50 -> 10)
        """
    end

end


text:
module TestRDA
    using Base.Test
    using DataFrames

    # R code generating test .rdas
        # df = data.frame(num=c(1.1, 2.2))
        # save(df, file='minimal.rda')

        # df['int'] = c(1L, 2L)
        # df['logi'] = c(TRUE, FALSE)
        # df['chr'] = c('ab', 'c')
        # df['factor'] = factor(df$chr)
        # #utf=c('Ж', '∰')) R handles it, read_rda doesn't.
        # save(df, file='types.rda')

        # df[2, ] = NA
        # df['chr'] = NULL  # NA characters breaking read_rda
        # save(df, file='NAs.rda')

        # names(df) = c('end', '!', '1', '%_B*\tC*')
        # save(df, file='names.rda')

    testdir = dirname(@__FILE__)

    df = DataFrame(num = [1.1, 2.2])
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/minimal.rda")["df"]), df)

    df[:int] = Int32[1, 2]
    df[:logi] = [true, false]
    df[:chr] = ["ab", "c"]
    df[:factor] = pool(df[:chr])
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/types.rda")["df"]), df)

    df[2, :] = NA
    df = df[:, [:num, :int, :logi, :factor]]  # (NA) chr breaks read_rda
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/NAs.rda")["df"]), df)

    rda_names = names(DataFrame(read_rda("$testdir/data/RDA/names.rda")["df"]))
    @test rda_names == [:_end, :x!, :x1, :_B_C_]
end


text:
using RDatasets, JDF, DataFrames

a = dataset("datasets", "iris");

first(a, 2)

@time jdffile = JDF.save("iris.jdf", a)
@time a2 = DataFrame(JDF.load("iris.jdf"))

all(names(a2) .== names(a)) # true
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)])) #true

a2_selected = DataFrame(JDF.load("iris.jdf", cols = [:Species, :SepalLength, :PetalWidth]))

jdf"path/to/JDF.jdf"

path_to_JDF = "path/to/JDF.jdf"
JDFFile(path_to_JDF)

afile = JDFFile("iris.jdf")

afile[!, :Species] # load Species column
afile[!, [:Species, :PetalLength]] # load Species and PetalLength column

afile[:, :Species] # load Species column
afile[:, [:Species, :PetalLength]] # load Species and PetalLength column

@view(afile[!, :Species]) # load Species column
@view(afile[!, [:Species, :PetalLength]]) # load Species and PetalLength column

using Tables
ajdf = JDFFile("iris.jdf")
Tables.columnaccess(ajdf)

Tables.columns(ajdf)

Tables.schema(ajdf)

getproperty(Tables.columns(ajdf), :Species)

jdffile = jdf"iris.jdf"
for col in eachcol(jdffile)
  # do something to col
  # where `col` is the content of one column of iris.jdf
end

jdffile = jdf"iris.jdf"
for (name, col) in zip(names(jdffile), eachcol(jdffile))
  # `name::Symbol` is the name of the column
  #  `col` is the content of one column of iris.jdf
end

using JDF, DataFrames
df = DataFrame(a = 1:3, b = 1:3)
savejdf(df, "plsdel.jdf")

names(jdf"plsdel.jdf") # [:a, :b]

nrow(jdf"plsdel.jdf") # 3

ncol(jdf"plsdel.jdf") # 2

size(jdf"plsdel.jdf") # (3, 2)

size(jdf"plsdel.jdf", 1) # 2

size(jdf"plsdel.jdf", 2) # 3

# clean up
rm("plsdel.jdf", force = true, recursive = true)

@time jdffile = ssavejdf("iris.jdf", a)
@time jdffile = sloadjdf("iris.jdf")

type_compress!(df)

type_compress!(df, compress_float = true)


text:
import Requires

function __init__()::Nothing
    Requires.@require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("pyplot.jl")

    Requires.@require Winston="bd07be1c-e76f-5ff0-9c0b-f51ef45303c6" include("winston.jl")

    return nothing
end


text:
VERSION >= v"0.5-" && __precompile__(true)

module Temporal
using Base.Dates
using Requests

export
    TS, ts, size, overlaps,
    ojoin, ijoin, ljoin, rjoin, merge, hcat, vcat, head, tail,
    nanrows, nancols, dropnan, fillnan, fillnan!, ffill!, bfill!, linterp!,
    numfun, arrfun, op,
    ones, zeros, trues, falses, isnan, countnz, sign,
    sum, mean, maximum, minimum, prod, cumsum, cumprod, diff, lag, nans,
    mondays, tuesdays, wednesdays, thursdays, fridays, saturdays, sundays,
    bow, eow, bom, eom, boq, eoq, boy, eoy, collapse, apply,
    tsread, tswrite, yahoo, quandl, quandl_auth, quandl_meta, quandl_search,
    acf

include("ts.jl")
include("indexing.jl")
include("combine.jl")
include("collapse.jl")
include("operations.jl")
include("models.jl")
include("slice.jl")
include("io.jl")

end # module


text:
using ITensors
using LinearAlgebra
using TBLIS

let
  χ = 500
  d = 3

  nthreads = 4

  l = Index(χ, "l")
  r = Index(χ, "r")
  s1 = Index(d, "s1")
  s2 = Index(d, "s2")
  h = Index(d, "h")

  indsA = (l, h, l')
  indsB = (s1, r, s2, l)
  indsC = (h, l', s1, r, s2)

  A = randomITensor(indsA)
  B = randomITensor(indsB)
  C = randomITensor(indsC)

  println("Contracting tensors A * B -> C, χ = $χ, d = $d")
  @show inds(A)
  @show inds(B)
  @show inds(C)

  #
  # Use BLAS
  #

  ITensors.disable_tblis()
  BLAS.set_num_threads(nthreads)

  println()
  println("Using BLAS with $nthreads threads")

  C_blas = copy(C)
  time_blas = @belapsed $C_blas .= $A .* $B samples = 100

  C_blas = copy(C)
  allocated_blas = @ballocated $C_blas .= $A .* $B samples = 100

  println()
  println("Time (BLAS) = ", time_blas, " seconds")
  println("Allocated (BLAS) = ", allocated_blas / 1024^2 , " MiB")

  #
  # Use TBLIS
  #

  ITensors.enable_tblis()
  TBLIS.set_num_threads(4)

  println()
  println("Using TBLIS with $(TBLIS.get_num_threads()) threads")

  C_tblis = copy(C)
  time_tblis = @belapsed $C_tblis .= $A .* $B samples = 100

  C_tblis = copy(C)
  allocated_tblis = @ballocated $C_tblis .= $A .* $B samples = 100

  println()
  println("Time (TBLIS) = ", time_tblis, " seconds")
  println("Allocated (TBLIS) = ", allocated_tblis * 1e-6, " MiB")

  println()
  @show norm(C_blas - C_tblis)

  println()
  println("Time (TBLIS) / Time (BLAS) = ", time_tblis / time_blas)
end


text:
# Part 1 - parse passports
passports = strip.(split(read("day04_input.txt", String), "\n\n"))
passports = map(p -> Dict(split.(split(p), ':')), passports)
required = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]
println("Part 1: There are $(sum(map(p -> issubset(required, keys(p)), passports))) valid passports")

# Part 2
function number_valid(p, key, re, mini, maxi)
    m = match(re, p[key])
    if isnothing(m) return false end
    mini <= parse(Int32, first(m.captures)) <= maxi
end

validators = [
    p -> number_valid(p, "byr", r"^(\d{4})$", 1920, 2002),
    p -> number_valid(p, "iyr", r"^(\d{4})$", 2010, 2020),
    p -> number_valid(p, "eyr", r"^(\d{4})$", 2020, 2030),
    p -> number_valid(p, "hgt", r"^(\d{3})cm$", 150, 193) || number_valid(p, "hgt", r"^(\d{2})in$", 59, 76),
    p -> !isnothing(match(r"^#[0-9a-f]{6}$", p["hcl"])),
    p -> p["ecl"] in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"],
    p -> !isnothing(match(r"^\d{9}$", p["pid"]))
]
is_valid(p) = issubset(required, keys(p)) && all(map(v -> v(p), validators))
println("Part 2: There are $(sum(is_valid.(passports))) valid passports")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
# Multivariate t-distribution

## Generic multivariate t-distribution class

abstract AbstractMvTDist <: ContinuousMultivariateDistribution

immutable GenericMvTDist{Cov<:AbstractPDMat} <: AbstractMvTDist
    df::Float64 # non-integer degrees of freedom allowed
    dim::Int
    zeromean::Bool
    μ::Vector{Float64}
    Σ::Cov

    function GenericMvTDist{Cov}(df::Float64, dim::Int, zmean::Bool, μ::Vector{Float64}, Σ::Cov)
      df > zero(df) || error("df must be positive")
      new(float64(df), dim, zmean, μ, Σ)
    end
end

function GenericMvTDist{Cov<:AbstractPDMat}(df::Float64, μ::Vector{Float64}, Σ::Cov, zmean::Bool)
    d = length(μ)
    dim(Σ) == d || throw(ArgumentError("The dimensions of μ and Σ are inconsistent."))
    GenericMvTDist{Cov}(df, d, zmean, μ, Σ)
end

function GenericMvTDist{Cov<:AbstractPDMat}(df::Float64, μ::Vector{Float64}, Σ::Cov)
    d = length(μ)
    dim(Σ) == d || throw(ArgumentError("The dimensions of μ and Σ are inconsistent."))
    GenericMvTDist{Cov}(df, d, allzeros(μ), μ, Σ)
end

function GenericMvTDist{Cov<:AbstractPDMat}(df::Float64, Σ::Cov)
    d = dim(Σ)
    GenericMvTDist{Cov}(df, d, true, zeros(d), Σ)
end

## Construction of multivariate normal with specific covariance type

typealias IsoTDist  GenericMvTDist{ScalMat}
typealias DiagTDist GenericMvTDist{PDiagMat}
typealias MvTDist GenericMvTDist{PDMat}

MvTDist(df::Float64, μ::Vector{Float64}, C::PDMat) = GenericMvTDist(df, μ, C)
MvTDist(df::Float64, C::PDMat) = GenericMvTDist(df, C)
MvTDist(df::Float64, μ::Vector{Float64}, Σ::Matrix{Float64}) = GenericMvTDist(df, μ, PDMat(Σ))
MvTDist(df::Float64, Σ::Matrix{Float64}) = GenericMvTDist(df, PDMat(Σ))

DiagTDist(df::Float64, μ::Vector{Float64}, C::PDiagMat) = GenericMvTDist(df, μ, C)
DiagTDist(df::Float64, C::PDiagMat) = GenericMvTDist(df, C)
DiagTDist(df::Float64, μ::Vector{Float64}, σ::Vector{Float64}) = GenericMvTDist(df, μ, PDiagMat(abs2(σ)))

IsoTDist(df::Float64, μ::Vector{Float64}, C::ScalMat) = GenericMvTDist(df, μ, C)
IsoTDist(df::Float64, C::ScalMat) = GenericMvTDist(df, C)
IsoTDist(df::Float64, μ::Vector{Float64}, σ::Real) = GenericMvTDist(df, μ, ScalMat(length(μ), abs2(float64(σ))))
IsoTDist(df::Float64, d::Int, σ::Real) = GenericMvTDist(df, ScalMat(d, abs2(float64(σ))))

## convenient function to construct distributions of proper type based on arguments

mvtdist(df::Float64, μ::Vector{Float64}, C::AbstractPDMat) = GenericMvTDist(df, μ, C)
mvtdist(df::Float64, C::AbstractPDMat) = GenericMvTDist(df, C)

mvtdist(df::Float64, μ::Vector{Float64}, σ::Real) = IsoTDist(df, μ, float64(σ))
mvtdist(df::Float64, d::Int, σ::Float64) = IsoTDist(d, σ)
mvtdist(df::Float64, μ::Vector{Float64}, σ::Vector{Float64}) = DiagTDist(df, μ, σ)
mvtdist(df::Float64, μ::Vector{Float64}, Σ::Matrix{Float64}) = MvTDist(df, μ, Σ)
mvtdist(df::Float64, Σ::Matrix{Float64}) = MvTDist(df, Σ)

# Basic statistics

dim(d::GenericMvTDist) = d.dim

mean(d::GenericMvTDist) = d.df>1 ? d.μ : NaN
mode(d::GenericMvTDist) = d.μ
modes(d::GenericMvTDist) = [mode(d)]

var(d::GenericMvTDist) = d.df>2 ? (d.df/(d.df-2))*diag(d.Σ) : Float64[NaN for i = 1:d.dim]
scale(d::GenericMvTDist) = full(d.Σ)
cov(d::GenericMvTDist) = d.df>2 ? (d.df/(d.df-2))*full(d.Σ) : NaN*ones(d.dim, d.dim)
invscale(d::GenericMvTDist) = full(inv(d.Σ))
invcov(d::GenericMvTDist) = d.df>2 ? ((d.df-2)/d.df)*full(inv(d.Σ)) : NaN*ones(d.dim, d.dim)
logdet_cov(d::GenericMvTDist) = d.df>2 ? logdet((d.df/(d.df-2))*d.Σ) : NaN

# For entropy calculations see "Multivariate t Distributions and their Applications", S. Kotz & S. Nadarajah
function entropy(d::GenericMvTDist)
    hdf, hdim = 0.5*d.df, 0.5*d.dim
    shdfhdim = hdf+hdim
    0.5*logdet(d.Σ)+hdim*log(d.df*pi)+lbeta(hdim, hdf)-lgamma(hdim)+shdfhdim*(digamma(shdfhdim)-digamma(hdf))
end

# evaluation (for GenericMvTDist)

function sqmahal(d::GenericMvTDist, x::Vector{Float64})
    z::Vector{Float64} = d.zeromean ? x : x - d.μ
    invquad(d.Σ, z)
end

function sqmahal!(r::Array{Float64}, d::GenericMvTDist, x::Matrix{Float64})
    if !(size(x, 1) == dim(d) && size(x, 2) == length(r))
        throw(ArgumentError("Inconsistent argument dimensions."))
    end
    z::Matrix{Float64} = d.zeromean ? x : bsubtract(x, d.μ, 1)
    invquad!(r, d.Σ, z)
end

# generic PDF evaluation (appliable to AbstractMvTDist)

insupport{T<:Real}(d::AbstractMvTDist, x::Vector{T}) = dim(d) == length(x) && allfinite(x)
insupport{T<:Real}(d::AbstractMvTDist, x::Matrix{T}) = dim(d) == size(x, 1) && allfinite(x)
insupport{G<:AbstractMvTDist,T<:Real}(::Type{G}, x::Vector{T}) = allfinite(x)
insupport{G<:AbstractMvTDist,T<:Real}(::Type{G}, x::Matrix{T}) = allfinite(x)

sqmahal(d::AbstractMvTDist, x::Matrix{Float64}) = sqmahal!(Array(Float64, size(x, 2)), d, x)

function logpdf(d::AbstractMvTDist, x::Vector{Float64})
  hdf, hdim = 0.5*d.df, 0.5*d.dim
  shdfhdim = hdf+hdim
  lgamma(shdfhdim)-lgamma(hdf)-hdim*log(d.df)-hdim*log(pi)-0.5*logdet(d.Σ)-shdfhdim*log(1+sqmahal(d, x)/d.df)
end

function logpdf!(r::Array{Float64}, d::AbstractMvTDist, x::Matrix{Float64})
  sqmahal!(r, d, x)
  hdf, hdim = 0.5*d.df, 0.5*d.dim
  shdfhdim = hdf+hdim
  for i = 1:size(x, 2)
    r[i] = lgamma(shdfhdim)-lgamma(hdf)-hdim*log(d.df)-hdim*log(pi)-0.5*logdet(d.Σ)-shdfhdim*log(1+r[i]/d.df)
  end
  r
end

function gradloglik(d::GenericMvTDist, x::Vector{Float64})
  z::Vector{Float64} = d.zeromean ? x : x - d.μ
  prz = invscale(d)*z
  -((d.df + d.dim) / (d.df + dot(z, prz))) * prz
end

# Sampling (for GenericMvTDist)

function rand!(d::GenericMvTDist, x::Vector{Float64})
  normdim = d.dim
  normd = GenericMvNormal{typeof(d.Σ)}(normdim, true, zeros(normdim), d.Σ)
  chisqd = Chisq(d.df)
  y = Array(Float64, d.dim)
  unwhiten!(normd.Σ, randn!(x))
  rand!(chisqd, y)
  y = sqrt(y/(d.df))
  divide!(x, y)
  if !d.zeromean
    add!(x, d.μ)
  end
  x
end

function rand!(d::GenericMvTDist, x::Matrix{Float64})
  normdim = d.dim
  normd = GenericMvNormal{typeof(d.Σ)}(normdim, true, zeros(normdim), d.Σ)
  chisqd = Chisq(d.df)
  y = Array(Float64, d.dim)
  unwhiten!(normd.Σ, randn!(x))
  rand!(chisqd, y)
  y = sqrt(y/(d.df))
  bdivide!(x, y, 1)
  if !d.zeromean
    badd!(x, d.μ, 1)
  end
  x
end


text:
function uri2filepath(uri::AbstractString)
    uri_path = normpath(URIParser.unescape(URIParser.URI(uri).path))

    if Sys.iswindows()
        if uri_path[1] == '\\' || uri_path[1] == '/'
            uri_path = uri_path[2:end]
        end
    end
    return uri_path
end

function filepath2uri(file::String)
    if Sys.iswindows()
        file = normpath(file)
        file = replace(file, "\\" => "/")
        file = URIParser.escape(file)
        file = replace(file, "%2F" => "/")
        return string("file:///", file)
    else
        file = normpath(file)
        file = URIParser.escape(file)
        file = replace(file, "%2F" => "/")
        return string("file://", file)
    end
end

function joinuriwithpath(uri::AbstractString, path::AbstractString)
    left_file_path = uri2filepath(uri)
    combined_path = joinpath(left_file_path, path)
    return filepath2uri(combined_path)
end

function should_file_be_linted(uri, server)
    !server.runlinter && return false

    uri_path = uri2filepath(uri)

    if length(server.workspaceFolders)==0
        return false
    else
        return any(i->startswith(uri_path, i), server.workspaceFolders)
    end
end

CompletionItemKind(t) = t in [:String, :AbstractString] ? 1 :
                                t == :Function ? 3 :
                                t == :DataType ? 7 :
                                t == :Module ? 9 : 6

SymbolKind(t) = t in [:String, :AbstractString] ? 15 :
                        t == :Function ? 12 :
                        t == :DataType ? 5 :
                        t == :Module ? 2 :
                        t == :Bool ? 17 : 13

# Find location of default datatype constructor
const DefaultTypeConstructorLoc= let def = first(methods(Int))
    Base.find_source_file(string(def.file)), def.line
end

function is_ignored(uri, server)
    fpath = uri2filepath(uri)
    fpath in server.ignorelist && return true
    for ig in server.ignorelist
        if !endswith(ig, ".jl")
            if startswith(fpath, ig)
                return true
            end
        end
    end
    return false
end

is_ignored(uri::URI2, server) = is_ignored(uri._uri, server)

# function toggle_file_lint(doc, server)
#     if doc._runlinter
#         doc._runlinter = false
#         empty!(doc.diagnostics)
#     else
#         doc._runlinter = true
#         # L = lint(doc, server)
#         # doc.diagnostics = L.diagnostics
#     end
#     publish_diagnostics(doc, server)
# end

function remove_workspace_files(root, server)
    for (uri, doc) in server.documents
        fpath = uri2filepath(uri._uri)
        doc._open_in_editor && continue
        if startswith(fpath, fpath)
            for folder in server.workspaceFolders
                if startswith(fpath, folder)
                    continue
                end
                delete!(server.documents, uri)
            end
        end
    end
end

function find_root(doc::Document, server)
    path = uri2filepath(doc._uri)
    for (uri1,f) in server.documents
        for incl in f.code.state.includes
            if path == incl.file
                if doc.code.index != incl.index
                    doc.code.index = incl.index
                    doc.code.nb = incl.pos
                end

                return find_root(f, server)
            end
        end
    end
    return doc
end

function Base.getindex(server::LanguageServerInstance, r::Regex)
    out = []
    for (uri,doc) in server.documents
        occursin(r, uri._uri) && push!(out, doc)
    end
    return out
end

function _offset_unitrange(r::UnitRange{Int}, first = true)
    return r.start-1:r.stop
end

function get_toks(doc, offset)
    ts = CSTParser.Tokenize.tokenize(doc._content)
    ppt = CSTParser.Tokens.RawToken(CSTParser.Tokens.ERROR, (0,0), (0,0), 1, 0, CSTParser.Tokens.NO_ERR, false)
    pt = CSTParser.Tokens.RawToken(CSTParser.Tokens.ERROR, (0,0), (0,0), 1, 0, CSTParser.Tokens.NO_ERR, false)
    t = CSTParser.Tokenize.Lexers.next_token(ts)
    if offset > length(doc._content)
        offset = sizeof(doc._content) - 1
    end

    while t.kind != CSTParser.Tokenize.Tokens.ENDMARKER
        if t.startbyte < offset <= t.endbyte + 1
            break
        end
        ppt = pt
        pt = t
        t = CSTParser.Tokenize.Lexers.next_token(ts)
    end
    return ppt, pt, t
end

function isvalidjlfile(path)
    isfile(path) &&
    endswith(path, ".jl") &&
    validchars(path)
end

function validchars(path)
    io = open(path)
    while !eof(io)
        c = read(io, Char)
        Base.ismalformed(c) && return false
    end
    close(io)
    return true
end

function goto_loc(x::EXPR, offset::Int, pos = 0)
    if x.args isa Vector{EXPR}
        for a in x.args
            if pos <= offset < pos + a.fullspan
                return goto_loc(a, offset, pos)
            else
                pos += a.fullspan
            end
        end
    end
    if offset == x.fullspan && !(x.parent isa EXPR)
        return last(x.args)
    end
    return x
end

function get_expr(x, offset, pos = 0)
    if pos > offset
        return nothing
    end
    if x.args !== nothing
        for a in x.args
            if pos < offset <= (pos + a.fullspan)
                return get_expr(a, offset, pos)
            end
            pos += a.fullspan
        end
    elseif (pos < offset <= (pos + x.fullspan)) || pos == 0
        return x
    end
end

function get_identifier(x, offset, pos = 0)
    if pos > offset
        return nothing
    end
    if x.args !== nothing
        for a in x.args
            if pos <= offset <= (pos + a.span)
                return get_identifier(a, offset, pos)
            end
            pos += a.fullspan
        end
    elseif x.typ === CSTParser.IDENTIFIER && (pos <= offset <= (pos + x.span)) || pos == 0
        return x
    end
end

@static if Sys.iswindows() && VERSION < v"1.3"
    if VERSION < v"1.2"
        function _splitdir_nodrive(a::String, b::String)
            m = match(r"^(.*?)([/\\]+)([^/\\]*)$",b)
            m === nothing && return (a,b)
            a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
            a, String(m.captures[3])
        end
    else
        _splitdir_nodrive = Base._splitdir_nodrive
    end
    function _dirname(path::String)
        m = match(r"^([^\\]+:|\\\\[^\\]+\\[^\\]+|\\\\\?\\UNC\\[^\\]+\\[^\\]+|\\\\\?\\[^\\]+:|)(.*)$"s, path)
        a, b = String(m.captures[1]), String(m.captures[2])
        _splitdir_nodrive(a,b)[1]
    end
else
    _dirname = dirname
end


text:
mutable struct AtomicMalloc
    blocksize::Int
    ptrs::Vector{Cptr}
    blockid::Int
    offset::Int

    function AtomicMalloc(blocksize::Int=1024*1000*1000)
        m = new(blocksize, Cptr[], 0, 0)
        # finalizer(m, dispose)
        m
    end
end

const ATOMIC_MALLOC = AtomicMalloc()

function (m::AtomicMalloc)(::Type{T}, dims::Dims{N}) where {T,N}
    bytesize = prod(dims) * sizeof(T)
    @assert 0 < bytesize <= m.blocksize

    if m.blockid == 0 || bytesize + m.offset > m.blocksize
        if m.blockid == length(m.ptrs)
            push!(m.ptrs, memalloc(m.blocksize))
        end
        m.blockid += 1
        ptr = m.ptrs[m.blockid]
        m.offset = 0
    else
        ptr = m.ptrs[m.blockid]
    end
    m.offset += bytesize
    CuArray(Ptr{T}(ptr), dims, getdevice())
end

function reset(m::AtomicMalloc)
    m.blockid = 1
    m.offset = 0
end

function dispose(m::AtomicMalloc)
    setdevice(m.device) do
        foreach(memfree, m.ptrs)
    end
end


text:
module Remarkable

include("Common.jl")
include("Tags.jl")
include("Articles.jl")
include("Talks.jl")
include("Remark/Remark.jl")
include("StaticSites/StaticSites.jl")

end  # module


text:
module EMLpower_plant
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	include("power_plant/CalorEspecifico.jl")
	include("power_plant/CoefGlobal_area.jl")
	include("power_plant/Dif_Pres.jl")
	include("power_plant/Dif_Temp.jl")
	include("power_plant/Eficiencia.jl")
	include("power_plant/EnergiaInterna.jl")
	include("power_plant/Entalpia.jl")
	include("power_plant/EntalpiaMol.jl")
	include("power_plant/Entropia.jl")
	include("power_plant/Fracao.jl")
	include("power_plant/Potencia.jl")
	include("power_plant/Pressao.jl")
	include("power_plant/MassaEspecifica.jl")
	include("power_plant/NoType.jl")
	include("power_plant/Temperatura.jl")
	include("power_plant/VazaoMassica.jl")
	include("power_plant/VolumeEspecifico.jl")
	include("power_plant/Pot_sinal.jl")
	include("power_plant/positive.jl")
	include("power_plant/mol.jl")
	include("power_plant/fraction.jl")
	include("power_plant/molweight.jl")
	include("power_plant/Volume.jl")
	include("power_plant/Reacao_mol.jl")
	include("power_plant/VazaoMolar.jl")
	include("power_plant/Corrente.jl")
	include("power_plant/CorrenteZ.jl")
	include("power_plant/CorrenteVap.jl")
	include("power_plant/Fonte2.jl")
	include("power_plant/Sumidouro.jl")
	include("power_plant/SumidouroQ.jl")
	include("power_plant/Turbina.jl")
	include("power_plant/Turbina_sangra.jl")
	include("power_plant/Condensador.jl")
	include("power_plant/Condensador_2alim.jl")
	include("power_plant/Tanque2.jl")
	include("power_plant/Tanque.jl")
	include("power_plant/Tanque3.jl")
	include("power_plant/Tanque4.jl")
	include("power_plant/Tanque4perdas.jl")
	include("power_plant/Desaerador5.jl")
	include("power_plant/Trocador.jl")
	include("power_plant/Torre.jl")
	include("power_plant/Bomba.jl")
	include("power_plant/Gerador_Vapor.jl")
	include("power_plant/Gerador_VaporMod.jl")
	include("power_plant/Gerador_Vapor_Simples.jl")
	include("power_plant/Gerador_Eletrico.jl")
	include("power_plant/Gerador_Eletrico_Simples.jl")
	include("power_plant/Splitter.jl")
	include("power_plant/Splitter4.jl")
	include("power_plant/Flash.jl")
	include("power_plant/ETA_CICLO.jl")
	include("power_plant/Equil_vap.jl")
end


text:
using PyCall

rclpy = pyimport("rclpy")
rosNode = pyimport("rclpy.node")
rosString = pyimport("std_msgs.msg")

@pydef mutable struct MinimialSubscriber <: rosNode.Node
    function __init__(self)
        rosNode.Node.__init__(self, "minimal_subscriber")
        function listener_callback(self, msg)
            self.get_logger().info("I heard: $(msg.data)")
        end
        self.subscription = self.create_subscription(rosString.String, "topic", msg -> listener_callback(self, msg), 10)
        self.subscription
    end
end

function main(args=nothing)
    rclpy.init(args=args)
    minimal_subscriber = MinimialSubscriber()
    rclpy.spin(minimal_subscriber)
    minimal_subscriber.destroy_node()
    rclpy.shutdown()
end