text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Getting started with data and plotting

# In many cases we might need to read data available in an external file rather
# than type it into Julia ourselves.

# This tutorial is concerned with reading tabular data into Julia. We'll cover
# basic plotting along the way.

# ## Where to get help

# * Plots.jl documentation: http://docs.juliaplots.org/latest/
# * CSV.jl documentation: http://csv.juliadata.org/stable
# * DataFrames.jl documentation: https://dataframes.juliadata.org/stable/

## We need this constant to point to where the data file are.
const DATA_DIR = joinpath(@__DIR__, "data");

# !!! note
#     There are multiple ways to read the same kind of data into Julia. This
#     tutorial focuses on DataFrames.jl because it provides the ecosystem to
#     work with most of the required file types in a straightforward manner.

# ### DataFrames.jl

# The `DataFrames` package provides a set of tools for working with tabular
# data. It is available through the Julia package system.
# ```julia
# using Pkg
# Pkg.add("DataFrames")
# ```

import DataFrames

# ### Plots.jl

# The `Plots` package provides a set of tools for plotting. It is available
# through the Julia package system.
# ```julia
# using Pkg
# Pkg.add("Plots")
# ```

import Plots

# ### What is a DataFrame?

# A DataFrame is a data structure like a table or spreadsheet. You can use it
# for storing and exploring a set of related data values. Think of it as a
# smarter array for holding tabular data.

# ## Reading Tabular Data into a DataFrame

# We will begin by reading data from different file formats into a DataFrame
# object.

# ### CSV files

# CSV and other delimited text files can be read by the CSV.jl package.

# ```julia
# Pkg.add("CSV")
# ```

import CSV

# To read a CSV file into a DataFrame, we use the `CSV.read` function.

csv_df = CSV.read(joinpath(DATA_DIR, "StarWars.csv"), DataFrames.DataFrame)

# Let's try plotting some of this data

Plots.scatter(
    csv_df.Weight,
    csv_df.Height,
    xlabel = "Weight",
    ylabel = "Height",
)

# That doesn't look right. What happened? If you look at the dataframe above, it
# read `Weight` in as a `String` column because there are "NA" fields. Let's
# correct that, by telling CSV to consider "NA" as `missing`.

csv_df = CSV.read(
    joinpath(DATA_DIR, "StarWars.csv"),
    DataFrames.DataFrame,
    missingstring="NA",
)

# Then let's re-plot our data

Plots.scatter(
    csv_df.Weight,
    csv_df.Height,
    title = "Height vs Weight of StarWars characters",
    xlabel = "Weight",
    ylabel = "Height",
    label = false,
    ylims = (0, 3),
)

# Better! Read the [CSV documentation](https://csv.juliadata.org/stable/) for
# other parsing options.

# DataFrames.jl supports manipulation using functions similar to pandas. For
# example, split the dataframe into groups based on eye-color:

by_eyecolor = DataFrames.groupby(csv_df, :Eyecolor)

# Then recombine into a single dataframe based on a function operating over the
# split dataframes:

eyecolor_count = DataFrames.combine(by_eyecolor) do df
    return DataFrames.nrow(df)
end

# We can rename columns:

DataFrames.rename!(eyecolor_count, :x1 => :count)

# Drop some missing rows:

DataFrames.dropmissing!(eyecolor_count, :Eyecolor)

# Then we can visualize the data:

sort!(eyecolor_count, :count, rev = true)
Plots.bar(
    eyecolor_count.Eyecolor,
    eyecolor_count.count,
    xlabel = "Eyecolor",
    ylabel = "Number of characters",
    label = false,
)

# ### Other Delimited Files

# We can also use the `CSV.jl` package to read any other delimited text file
# format.

# By default, CSV.File will try to detect a file's delimiter from the first 10
# lines of the file.

# Candidate delimiters include `','`, `'\t'`, `' '`, `'|'`, `';'`, and `':'`. If
# it can't auto-detect the delimiter, it will assume `','`.

# Let's take the example of space separated data.

ss_df = CSV.read(joinpath(DATA_DIR, "Cereal.txt"), DataFrames.DataFrame)

# We can also specify the delimiter by passing the `delim` argument.

delim_df = CSV.read(
    joinpath(DATA_DIR, "Soccer.txt"),
    DataFrames.DataFrame,
    delim = "::",
)

# ## Working with DataFrames

# Now that we have read the required data into a DataFrame, let us look at some
# basic operations we can perform on it.

# ### Querying Basic Information

# The `size` function gets us the dimensions of the DataFrame.

DataFrames.size(ss_df)

# We can also use the `nrow` and `ncol` functions to get the number of rows and
# columns respectively.

DataFrames.nrow(ss_df), DataFrames.ncol(ss_df)

# The `describe` function gives basic summary statistics of data in a DataFrame.

DataFrames.describe(ss_df)

# Names of every column can be obtained by the `names` function.

DataFrames.names(ss_df)

# Corresponding data types are obtained using the broadcasted `eltype` function.

eltype.(ss_df)

# ### Accessing the Data

# Similar to regular arrays, we use numerical indexing to access elements of a
# DataFrame.

csv_df[1, 1]

# The following are different ways to access a column.

csv_df[!, 1]

#-

csv_df[!, :Name]

#-

csv_df.Name

#-

csv_df[:, 1] # Note that this creates a copy.

# The following are different ways to access a row.

csv_df[1:1, :]

#-

csv_df[1, :] # This produces a DataFrameRow.

# We can change the values just as we normally assign values.

# Assign a range to scalar.

csv_df[1:3, :Height] .= 1.83

# Vector to equal length vector.

csv_df[4:6, :Height] = [1.8, 1.6, 1.8]

#-

csv_df

# !!! tip
#     There are a lot more things which can be done with a DataFrame. Read the
#     [docs](https://juliadata.github.io/DataFrames.jl/stable/) for more
#     information.

# For information on dplyr-type syntax:
 # * Read: https://dataframes.juliadata.org/stable/man/querying_frameworks/
 # * Check out DataFramesMeta.jl: https://github.com/JuliaData/DataFramesMeta.jl

max_stars_repo_name: odow/jump-training-materials

max_issues_repo_name: odow/jump-training-materials

max_forks_repo_name: odow/jump-training-materials


text:
infile = "./src/Day8/input.txt"
input = readlines(infile)
insplit = split.(input, " ")
opcodes = first.(insplit)
vals = parse.(Int,last.(insplit))

mutable struct Line
    opcode::Symbol
    val::Int64
    count::Int64
end

lines = Dict{Int,Line}()
for i in 1:length(opcodes)
    l = Line(Symbol(opcodes[i]), vals[i], 0)
    lines[i] = l
end

function run(lines)
    currentno = 1
    maxrun = 0
    acc = 0
    while maxrun <= 1
        if !haskey(lines, currentno)
            return acc, true
        end
        currentline = lines[currentno]
        #println("Reached $currentno: $currentline")
        maxrun = currentline.count
        if currentline.count == 1
            # about to run a second time
            return acc, false
        end
        if currentline.opcode == :jmp
            currentno += currentline.val
            currentline.count += 1
        elseif currentline.opcode == :acc
            acc += currentline.val
            currentline.count += 1
            currentno += 1
        elseif currentline.opcode == :nop
            currentline.count += 1
            currentno += 1
        else
            error("Was on line $currentno and died")
        end
    end
    return 0, false
end

l = deepcopy(lines)

run(l)

function fix(lines)
    maybechange = findall(x -> x.opcode ∈ (:jmp, :nop), lines)
    for ix in maybechange
        fixed = deepcopy(lines)
        fixline = fixed[ix]
        if fixline.opcode == :jmp
            fixline.opcode = :nop
        else
            fixline.opcode = :jmp
        end
        val, status = run(fixed)
        if status
            return val, ix
        end
    end
    error("Nope")
end

l = deepcopy(lines)
fix(l)

ext: jl


text:
add_one(args...) = 1.0

""" Instantiate ModelForcing. """
function test_forcing(fld)
    kwarg = Dict(fld => add_one)
    forcing = ModelForcing(; kwarg...)
    f = getfield(forcing, fld)
    return f() == 1.0
end

""" Take one time step with three forcing functions on u, v, w. """
function time_step_with_forcing_functions(arch)
    @inline Fu(i, j, k, grid, clock, state) = @inbounds ifelse(k == grid.Nz, -state.velocities.u[i, j, k] / 60, 0)
    @inline Fv(i, j, k, grid, clock, state) = @inbounds ifelse(k == grid.Nz, -state.velocities.v[i, j, k] / 60, 0)
    @inline Fw(i, j, k, grid, clock, state) = @inbounds ifelse(k == grid.Nz, -state.velocities.w[i, j, k] / 60, 0)

    forcing = ModelForcing(u=Fu, v=Fv, w=Fw)

    grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=forcing)
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with ParameterizedForcing forcing functions. """
function time_step_with_parameterized_forcing(arch)
    @inline Fu_func(i, j, k, grid, clock, state, params) = @inbounds ifelse(k == grid.Nz, -state.velocities.u[i, j, k] / params.τ, 0)
    @inline Fv_func(i, j, k, grid, clock, state, params) = @inbounds ifelse(k == grid.Nz, -state.velocities.v[i, j, k] / params.τ, 0)
    @inline Fw_func(i, j, k, grid, clock, state, params) = @inbounds ifelse(k == grid.Nz, -state.velocities.w[i, j, k] / params.τ, 0)

    Fu = ParameterizedForcing(Fu_func, (τ=60,))
    Fv = ParameterizedForcing(Fv_func, (τ=60,))
    Fw = ParameterizedForcing(Fw_func, (τ=60,))

    forcing = ModelForcing(u=Fu, v=Fv, w=Fw)

    grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=forcing)
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with forcing functions containing sin and exp functions. """
function time_step_with_forcing_functions_sin_exp(arch)
    @inline Fu(i, j, k, grid, clock, state) = @inbounds sin(grid.xC[i])
    @inline FT(i, j, k, grid, clock, state) = @inbounds exp(-state.tracers.T[i, j, k])

    forcing = ModelForcing(u=Fu, T=FT)

    grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=forcing)
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a SimpleForcing forcing function. """
function time_step_with_simple_forcing(arch)
    u_forcing = SimpleForcing((x, y, z, t) -> sin(x))
    grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=ModelForcing(u=u_forcing))
    time_step!(model, 1, euler=true)
    return true
end

""" Take one time step with a SimpleForcing forcing function with parameters. """
function time_step_with_simple_forcing_parameters(arch)
    u_forcing = SimpleForcing((x, y, z, t, p) -> sin(p.ω * x), parameters=(ω=π,))
    grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=ModelForcing(u=u_forcing))
    time_step!(model, 1, euler=true)
    return true
end

@testset "Forcing" begin
    @info "Testing forcings..."

    @testset "Forcing function initialization" begin
        @info "  Testing forcing function initialization..."
        for fld in (:u, :v, :w, :T, :S)
            @test test_forcing(fld)
        end
    end

    for arch in archs
        @testset "Forcing function time stepping [$(typeof(arch))]" begin
            @info "  Testing forcing function time stepping [$(typeof(arch))]..."
            @test time_step_with_forcing_functions(arch)
            @test time_step_with_parameterized_forcing(arch)
            @test time_step_with_forcing_functions_sin_exp(arch)
            @test time_step_with_simple_forcing(arch)
            @test time_step_with_simple_forcing_parameters(arch)
        end
    end
end


text:
using Revise
using PortfolioBedtest
import PortfolioBedtest: order, execute
using Dates
using Statistics
using Setfield
using Convex, SCS # It can be simplifed, but existing solutions are fine

########################################
# Auxiliary functions
########################################
function get_optim_weights(w, risk, ret, λ)
    p = minimize(λ * risk - (1 - λ) * ret,
                 sum(w) == 1,
                 0.0 <= w,
                 w <= 1.0)
    solve!(p, () -> SCS.Optimizer(verbose = false))

    return sqrt(evaluate(risk)), evaluate(w)
end

function find_weights(M, r, risk0, tol = 1e-8)
    w = Variable(length(r))
    risk = quadform(w, M)
    ret = dot(w, r)
    l1 = 0.01
    risk1, _ = get_optim_weights(w, risk, ret, l1)
    l2 = 0.99
    risk2, _ = get_optim_weights(w, risk, ret, l2)
    l3 = (l1 + l2)/2
    risk3, _ = get_optim_weights(w, risk, ret, l3)
    while abs(l1 - l2) > tol
        if risk1 >= risk0 >= risk3
            l2 = l3
        else
            l1 = l3
        end
        l3 = (l1 + l2) / 2
        risk3, _ = get_optim_weights(w, risk, ret, l3)
    end
    risk3, w3 = get_optim_weights(w, risk, ret, l3)

    return risk3, w3
end

function build(assts, rng, risk = 0.1)
    # not very efficient, but ok
    wnd = assts[rng]
    wnd2 = wnd[2:end] ./ wnd[1:end - 1]
    ret = exp.((252 * mean(map.(log, wnd2))).event) .- 1
    # It can be done in row-wise format, but optimization is not important here
    # So we transform data to column-wise matrix
    M = 252 .* cov(hcat(map(i -> getindex.(getfield.(wnd2, :event), i), 1:length(ret))...))
    _, w = find_weights(M, ret, risk)

    return w
end

function create_signal(assts, risk = 0.1, period = 36)
    indx = filter_eom_indices(assts)
    w = build(assts, indx[1]:indx[period+1], risk)
    w2 = tuple(w...)
    res = [TimedEvent(assts[indx[period + 1]].ts, w2)]
    for i in period+2:length(indx)
        w = build(assts, indx[i - period]:indx[i], risk)
        w2 = tuple(w...)
        push!(res, TimedEvent(assts[indx[i]].ts, w2))
    end
    Signal(res)
end

########################################
# Simulation auxiliary structures
########################################
struct RatioStrategy end

struct Order{T} <: AbstractOrder
    shares::T
end

function order(strat::RatioStrategy, broker, signal, prices)
    a = assets(broker, prices)
    ntuple(length(signal)) do i
        floor(Int, a * signal[i] / prices[i])
    end |> Order
end

function execute(broker, ord::Order, prices)
    a = assets(broker, prices)
    for i in 1:length(ord.shares)
        a -= ord.shares[i] * prices[i]
    end
    Broker(a, ord.shares)
end

########################################
# Calculations
########################################

# ETF
spyts = yahoo("SPY") |> x -> map(x -> TimedEvent(x.Date, getproperty(x, Symbol("Adj Close"))), x);
aggts = yahoo("AGG") |> x -> map(x -> TimedEvent(x.Date, getproperty(x, Symbol("Adj Close"))), x);
gldts = yahoo("GLD") |> x -> map(x -> TimedEvent(x.Date, getproperty(x, Symbol("Adj Close"))), x);

assts = glue(spyts, aggts, gldts) do v1, v2, k
    k == 0 && return (v1..., v2...)
    k == 1 && return (v1..., map(_ -> NaN, v2)...)
    return (map(_ -> NaN, v1)..., v2...)
end

assts = filter(x -> !any(isnan.(x.event)), assts)
signal = create_signal(assts)
broker = Broker(1_000_000.0, ntuple(_ -> 0, length(assts[1].event)))
strat = RatioStrategy()
recorder = Recorder{Date, typeof((; shares = broker.shares, total = 0.0))}([])

sim(broker, strat, assts, signal, recorder)

########################################
# Volatility 7.5%
########################################
vol75_signal = create_signal(assts, 0.075)
vol75_recorder = Recorder{Date, typeof((; shares = broker.shares, total = 0.0))}([])

sim(broker, strat, assts, vol75_signal, vol75_recorder)

########################################
# 60/40 comparison
########################################

sfsignal = Signal(mutate.(_ -> (0.6, 0.4, 0.0), signal.signal))
sfrecorder = Recorder{Date, typeof((; shares = broker.shares, total = 0.0))}([])

sim(broker, strat, assts, sfsignal, sfrecorder)

########################################
# Pure SPY
########################################
spysignal = Signal(mutate.(_ -> (1.0, 0.0, 0.0), signal.signal))
spyrecorder = Recorder{Date, typeof((; shares = broker.shares, total = 0.0))}([])

sim(broker, strat, assts, spysignal, spyrecorder)

########################################
# Plotting
########################################
using Plots
tsts = getfield.(signal.signal, :ts)
plot(tsts, getfield.(mutate.(x -> x[1], signal.signal), :event), legend = :outerright, label = "SPY")
plot!(tsts, getfield.(mutate.(x -> x[2], signal.signal), :event), label = "AGG")
plot!(tsts, getfield.(mutate.(x -> x[3], signal.signal), :event), label = "GLD")

# Equities
equity_curve = mutate.(x -> x.total, recorder.data)
plot(getfield.(equity_curve, :ts), getfield.(equity_curve, :event), legend = :topleft, label = "10% volatilty rebalance", linewidth = 0.5)

sfequity_curve = mutate.(x -> x.total, sfrecorder.data)
plot!(getfield.(sfequity_curve, :ts), getfield.(sfequity_curve, :event), label = "60/40 SPY/AGG", linewidth = 0.5)

spyequity_curve = mutate.(x -> x.total, spyrecorder.data)
plot!(getfield.(spyequity_curve, :ts), getfield.(spyequity_curve, :event), label = "100% SPY", linewidth = 0.5)

vol75_equity_curve = mutate.(x -> x.total, vol75_recorder.data)
plot!(getfield.(vol75_equity_curve, :ts), getfield.(vol75_equity_curve, :event))

########################################
# Volatility spectre
########################################
function build_vol_plot(assts, rng)
    isfirst = true
    broker = Broker(1_000_000.0, ntuple(_ -> 0, length(assts[1].event)))
    strat = RatioStrategy()
    recorder = Recorder{Date, typeof((; shares = broker.shares, total = 0.0))}([])
    local g
    for vol in rng
        signal = create_signal(assts, vol)
        sim(broker, strat, assts, signal, recorder)
        equity_curve = mutate.(x -> x.total, recorder.data)
        g = if isfirst
            isfirst = false
            plot(getfield.(equity_curve, :ts), getfield.(equity_curve, :event), legend = :topleft, label = "Risk: $vol")
        else
            plot!(getfield.(equity_curve, :ts), getfield.(equity_curve, :event), label = "Risk: $vol")
        end
    end

    return g
end

build_vol_plot(assts, 0.06:0.02:0.12)

########################################
# Sliding window spectre
########################################

function build_win_plot(assts, rng, risk = 0.1)
    isfirst = true
    broker = Broker(1_000_000.0, ntuple(_ -> 0, length(assts[1].event)))
    strat = RatioStrategy()
    recorder = Recorder{Date, typeof((; shares = broker.shares, total = 0.0))}([])
    local g
    for wnd in rng
        signal = create_signal(assts, risk, wnd)
        sim(broker, strat, assts, signal, recorder)
        equity_curve = mutate.(x -> x.total, recorder.data)
        g = if isfirst
            isfirst = false
            plot(getfield.(equity_curve, :ts), getfield.(equity_curve, :event), legend = :topleft, label = "Window: $wnd", linewidth = 0.5)
        else
            plot!(getfield.(equity_curve, :ts), getfield.(equity_curve, :event), label = "Window: $wnd", linewidth = 0.5)
        end
    end

    return g
end

build_win_plot(assts, 12:6:48)
build_win_plot(assts, 20:2:28)
build_win_plot(assts, 20:2:28, 0.12)


text:
"""
$(TYPEDEF)
"""
struct NonlinearSolution{T,N,uType,R,P,A,O} <: AbstractNonlinearSolution{T,N}
  u::uType
  resid::R
  prob::P
  alg::A
  retcode::Symbol
  original::O
end

const SteadyStateSolution = NonlinearSolution

function build_solution(prob::AbstractNonlinearProblem,
                        alg,u,resid;calculate_error = true,
                        retcode = :Default,
                        original = nothing,
                        kwargs...)

  T = eltype(eltype(u))
  N = length((size(prob.u0)...,))

  NonlinearSolution{T,N,typeof(u),typeof(resid),
                    typeof(prob),typeof(alg),typeof(original)}(
                    u,resid,prob,alg,retcode,original)
end

function sensitivity_solution(sol::AbstractNonlinearProblem,u)
  T = eltype(eltype(u))
  N = length((size(sol.prob.u0)...,))

  NonlinearSolution{T,N,typeof(u),typeof(sol.resid),
                    typeof(sol.prob),typeof(sol.alg),
                    typeof(sol.original)}(
                    u,sol.resid,sol.prob,sol.alg,sol.retcode,sol.original)
end


text:
# Used for objectives and solvers where no gradient is available/exists
mutable struct NonDifferentiable{TF,TX} <: AbstractObjective
    f
    F::TF
    x_f::TX
    f_calls::Vector{Int}
end

# These could be the same if we could restrict g below not to be an AbstractArray
function NonDifferentiable(f, x::AbstractArray, F::Real = real(zero(eltype(x))); inplace = true)
    xnans = x_of_nans(x)
    NonDifferentiable{typeof(F),typeof(xnans)}(f, F, xnans, [0,])
end
function NonDifferentiable(f, x::AbstractArray, F::AbstractArray; inplace = true)
    f = !inplace && (F isa AbstractArray) ? f!_from_f(f, F, inplace) : f
    xnans = x_of_nans(x)
    NonDifferentiable{typeof(F),typeof(xnans)}(f, F, xnans, [0,])
end

# this is the g referred to above!
NonDifferentiable(f, g,        x::AbstractArray, F::Union{AbstractArray, Real} = real(zero(eltype(x)))) = NonDifferentiable(f, x, F)
NonDifferentiable(f, g, h,     x::TX, F) where TX  = NonDifferentiable(f, x, F)
NonDifferentiable(f, g, fg, h, x::TX, F) where TX  = NonDifferentiable(f, x, F)


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    eltypes = (N0f8, N0f16, Float32, Float64)        # eltypes of parametric colors
    pctypes = (Gray, RGB, AGray, GrayA, ARGB, RGBA)  # parametric colors
    cctypes = (Gray24, AGray32, RGB24, ARGB32)       # non-parametric colors
    realtypes = (Float16, Float32, Float64, Int)     # types for mixed Normed/Real operations
    # Constructors
    for R in realtypes
        for T in eltypes
            @assert precompile(Tuple{Type{Gray{T}},R})
            @assert precompile(Tuple{Type{AGray{T}},R})
            @assert precompile(Tuple{Type{GrayA{T}},R})
            @assert precompile(Tuple{Type{AGray{T}},R,R})
            @assert precompile(Tuple{Type{GrayA{T}},R,R})
            @assert precompile(Tuple{Type{RGB{T}},R,R,R})
            @assert precompile(Tuple{Type{RGBA{T}},R,R,R})
            @assert precompile(Tuple{Type{RGBA{T}},R,R,R,R})
            @assert precompile(Tuple{Type{ARGB{T}},R,R,R})
            @assert precompile(Tuple{Type{ARGB{T}},R,R,R,R})
        end
        @assert precompile(Tuple{Type{Gray},R})
        @assert precompile(Tuple{Type{AGray},R})
        @assert precompile(Tuple{Type{GrayA},R})
        @assert precompile(Tuple{Type{AGray},R,R})
        @assert precompile(Tuple{Type{GrayA},R,R})
        @assert precompile(Tuple{Type{RGB},R,R,R})
        @assert precompile(Tuple{Type{RGBA},R,R,R})
        @assert precompile(Tuple{Type{RGBA},R,R,R,R})
        @assert precompile(Tuple{Type{ARGB},R,R,R})
        @assert precompile(Tuple{Type{ARGB},R,R,R,R})
        @assert precompile(Tuple{Type{Gray24},R})
        @assert precompile(Tuple{Type{AGray32},R})
        @assert precompile(Tuple{Type{AGray32},R,R})
        @assert precompile(Tuple{Type{RGB24},R,R,R})
        @assert precompile(Tuple{Type{ARGB32},R,R,R})
        @assert precompile(Tuple{Type{ARGB32},R,R,R,R})
    end
    @assert precompile(Tuple{Type{Gray24},UInt32,Type{Val{true}}})
    @assert precompile(Tuple{Type{AGray32},UInt32,Type{Val{true}}})
    @assert precompile(Tuple{Type{RGB24},UInt32,Type{Val{true}}})
    @assert precompile(Tuple{Type{ARGB32},UInt32,Type{Val{true}}})
    # LCHab is used by distinguishable_colors
    for R in (Float32, Float64)
        @assert precompile(Tuple{Type{LCHab},R,R,R})
    end
    # convert
    for T1 in eltypes, T2 in eltypes
        @assert precompile(Tuple{typeof(convert),Type{Gray{T1}},Gray{T2}})
        @assert precompile(Tuple{typeof(convert),Type{RGB{T1}},RGB{T2}})
        @assert precompile(Tuple{typeof(convert),Type{RGB{T1}},Gray{T2}})
        @assert precompile(Tuple{typeof(convert),Type{RGBA{T1}},RGB{T2}})
        @assert precompile(Tuple{typeof(convert),Type{RGBA{T1}},RGBA{T2}})
        @assert precompile(Tuple{typeof(convert),Type{RGBA{T1}},ARGB{T2}})
        @assert precompile(Tuple{typeof(convert),Type{RGBA{T1}},Gray{T2}})
        @assert precompile(Tuple{typeof(convert),Type{ARGB{T1}},RGB{T2}})
        @assert precompile(Tuple{typeof(convert),Type{ARGB{T1}},RGBA{T2}})
        @assert precompile(Tuple{typeof(convert),Type{ARGB{T1}},ARGB{T2}})
        @assert precompile(Tuple{typeof(convert),Type{ARGB{T1}},Gray{T2}})
    end
    # traits
    for T in eltypes, C in pctypes
        @assert precompile(Tuple{typeof(alpha),C{T}})
        C <: AbstractGray && @assert precompile(Tuple{typeof(gray),C{T}})
        C <: AbstractRGB || continue
        @assert precompile(Tuple{typeof(red),C{T}})
        @assert precompile(Tuple{typeof(green),C{T}})
        @assert precompile(Tuple{typeof(blue),C{T}})
    end
    for C in cctypes
        @assert precompile(Tuple{typeof(alpha),C})
        C <: AbstractGray && @assert precompile(Tuple{typeof(gray),C})
        C <: AbstractRGB || continue
        @assert precompile(Tuple{typeof(red),C})
        @assert precompile(Tuple{typeof(green),C})
        @assert precompile(Tuple{typeof(blue),C})
    end
    # ccolor typically gets compiled as part of other things
    # hash
    for T in eltypes, C in pctypes
        @assert precompile(Tuple{typeof(hash),Type{C{T}},UInt})
    end
    for C in cctypes
        @assert precompile(Tuple{typeof(hash),Type{C},UInt})
    end
    # rand
    for T in eltypes, C in pctypes
        @assert precompile(Tuple{typeof(rand),Type{C{T}},Tuple{Int}})
        @assert precompile(Tuple{typeof(rand),Type{C{T}},Tuple{Int,Int}})
        @assert precompile(Tuple{typeof(rand),Type{C{T}},Int})
        @assert precompile(Tuple{typeof(rand),Type{C{T}},Int,Int})
    end
    for C in cctypes
        @assert precompile(Tuple{typeof(rand),Type{C},Tuple{Int}})
        @assert precompile(Tuple{typeof(rand),Type{C},Tuple{Int,Int}})
    end
    # mapc, reducec, and mapreducec are not really precompilable,
    # since they will be specialized for f
end


text:
using FinEconResearchUtils
using Test

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


text:
"""
    nextoffset = flatten!(v, domain, val, offset=1)
Flatten `val`, taking values in `domain`, and transform to Cartesian coordinates. Insert the transformed value into `v` starting at `offset`. Returns the next offset to use.
"""
flatten!(v, domain, val) = flatten!(v, domain, val, 1)

# recursively apply `flatten`
function flatten!(v, ::Type{T}, val::T, offset) where {T}
  map(zip(fieldnames(T), domains(T))) do (fieldname, domain)
    offset = flatten!(v, domain, getfield(val, fieldname), offset)
  end
  return offset
end

function flatten!(v, ::Type{T}, val::T, offset) where {T<:Real}
  v[offset] = val
  return offset+1
end

"""
    v = flatten(domain, val)
Flatten `val`, taking values in `domain`, and transform to Cartesian coordinates. Returns a vector `v` of length `dimension(domain)`.
"""
function flatten(domain, val)
  v = zeros(Float64, dimension(domain))
  flatten!(v, domain, val)
  v
end

"""
    val, nextoffset = unflatten(v, domain, offset)
Transform from `v`, starting at `offset`, stored in Cartesian coordinates to value in `domain`. Returns the transformed value and the next offset.
"""
function unflatten end

"""
    val = unflatten(v, domain)
Transform from `v` stored in Cartesian coordinates to value in `domain`. Returns the transformed value.
"""
unflatten(v, domain) = first(unflatten(v, domain, 1))

function unflatten(v, ::Type{T}, offset) where {T}
  vals = map(domains(T)) do domain
    val, offset = unflatten(v, domain, offset)
    val
  end
  return T(vals...), offset
end

function unflatten(v, ::Type{T}, offset) where {T<:Real}
  v[offset], offset+1
end


text:
@testset "HW5" begin
    using DMUStudent.HW5: mc
    using CommonRLInterface
    using Compose

    @test mc isa AbstractEnv

    A = actions(mc)
    done = false
    rsum = 0.0
    reset!(mc)
    s = observe(gw)
    t = 1
    γ = 0.95
    max_t = 200
    while t <= max_t && !terminated(mc)
        a = rand(A)

        r = act!(mc, a)
        sp = observe(mc)

        rsum += γ^t*r

        s = sp
        t += 1
    end
    @test -5.0 <= rsum <= 1000.0

    # rendering
    rndr = render(mc)
    @test rndr isa Context
    sprint((io, x)->show(io, MIME("text/html"), x), rndr)

    @test HW5.evaluate(s->0.0).score < 0.0
    @test HW5.evaluate(s->0.0, "zachary.sunberg@colorado.edu").score < 0.0
    @test HW5.evaluate(s->rand(actions(mc)), "zachary.sunberg@colorado.edu").score < 0.0

    # test whether it has enough energy
    disc = 1.0
    rsum = 0.0
    mc.s = [1.56, 0.0]
    while !terminated(mc) && disc > 0.001
        r = act!(mc, -1.0)
        rsum += disc*r
        disc *= 0.99
    end
    @test rsum > 200.0

    disc = 1.0
    rsum = 0.0
    mc.s = [1.4, 0.0]
    while !terminated(mc) && disc > 0.001
        r = act!(mc, -1.0)
        rsum += disc*r
        disc *= 0.99
    end
    @test rsum < 0.0
end


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init()

using ClimateMachine.Atmos
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.DGmethods.NumericalFluxes
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.MoistThermodynamics
using ClimateMachine.VariableTemplates

using Distributions
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra

using CLIMAParameters
using CLIMAParameters.Planet: cp_d, MSLP, grav, LH_v0
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.DGmethods:
    vars_state_conservative,
    vars_state_auxiliary,
    vars_integrals,
    vars_reverse_integrals,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ClimateMachine.DGmethods: boundary_state!
import ClimateMachine.Atmos: flux_second_order!

"""
  Initial Condition for Taylor-Green vortex (LES)
@article{taylorGreen1937,
author = {Taylor, G. I. and Green, A. E.},
title = {Mechanisms of production of small eddies from large ones},
journal = {Proc. Roy. Soc. A},
volume = {158}
number = {895},
year = {1937},
doi={doi.org/10.1098/rspa.1937.0036},
}

@article{rafeiEtAl2018,
author = {Rafei, M.E. and K\"on\"oszy, L. and Rana, Z.},
title = {Investigation of Numerical Dissipation in Classicaland Implicit Large Eddy Simulations,
journal = {Aerospace},
year = {2018},
}

@article{bullJameson2014,
author = {Bull, J.R. and Jameson, A.},
title = {Simulation of the Compressible {Taylor Green} Vortex
using High-Order Flux Reconstruction Schemes},
journal = {AIAA Aviation 7th AIAA theoretical fluid mechanics conference},
year = {2014},
}

"""
function init_greenvortex!(bl, state, aux, (x, y, z), t)
    # Problem float-type
    FT = eltype(state)

    # Unpack constant parameters
    R_gas::FT = R_d(bl.param_set)
    c_p::FT = cp_d(bl.param_set)
    c_v::FT = cv_d(bl.param_set)
    p0::FT = MSLP(bl.param_set)
    _grav::FT = grav(bl.param_set)
    γ::FT = c_p / c_v

    # Compute perturbed thermodynamic state:
    ρ = FT(1.178)      # density
    ρu = SVector(FT(0), FT(0), FT(0))  # momentum
    #State (prognostic) variable assignment
    e_pot = FT(0)# potential energy
    Pinf = 101325
    Uzero = FT(100)
    p = Pinf + (ρ * Uzero / 16) * (2 + cos(z)) * (cos(x) + cos(y))
    u = Uzero * sin(x) * cos(y) * cos(z)
    v = -Uzero * cos(x) * sin(y) * cos(z)
    e_kin = 0.5 * (u^2 + v^2)
    T = p / (ρ * R_gas)
    e_int = internal_energy(bl.param_set, T, PhasePartition(FT(0)))
    ρe_tot = ρ * (e_kin + e_pot + e_int)
    # Assign State Variables
    state.ρ = ρ
    state.ρu = SVector(FT(ρ * u), FT(ρ * v), FT(0))
    state.ρe = ρe_tot
end

"""
    config_greenvortex(FT, N, resolution, xmax, ymax, zmax)

Arguments
- FT = Floating-point type. Currently `Float64` or `Float32`
- N  = DG Polynomial order
- resolution = 3-component Tuple (Δx, Δy, Δz) with effective resolutions for Cartesian directions
- xmax, ymax, zmax = Domain maximum extents. Assumes (0,0,0) to be the domain minimum extents unless otherwise specified.

Returns
- `config` = Object using the constructor for the `AtmosLESConfiguration`
"""
function config_greenvortex(
    FT,
    N,
    resolution,
    xmax,
    ymax,
    zmax,
    xmin,
    ymin,
    zmin,
)

    ode_solver = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    _C_smag = FT(C_smag(param_set))
    model = AtmosModel{FT}(
        AtmosLESConfigType,                 # Flow in a box, requires the AtmosLESConfigType
        orientation = NoOrientation(),
        param_set;                          # Parameter set corresponding to earth parameters
        turbulence = Vreman(_C_smag),       # Turbulence closure model
        moisture = DryModel(),
        source = (),
        init_state_conservative = init_greenvortex!,             # Apply the initial condition
    )

    # Finally,  we pass a `Problem Name` string, the mesh information, and the model type to  the [`AtmosLESConfiguration`] object.
    config = ClimateMachine.AtmosLESConfiguration(
        "GreenVortex",       # Problem title [String]
        N,                       # Polynomial order [Int]
        resolution,              # (Δx, Δy, Δz) effective resolution [m]
        xmax,                    # Domain maximum size [m]
        ymax,                    # Domain maximum size [m]
        zmax,                    # Domain maximum size [m]
        param_set,               # Parameter set.
        init_greenvortex!,      # Function specifying initial condition
        boundary = ((0, 0), (0, 0), (0, 0)),
        periodicity = (true, true, true),
        xmin = xmin,
        ymin = ymin,
        zmin = zmin,
        solver_type = ode_solver,# Time-integrator type
        model = model,           # Model type
    )
    return config
end
# Here we define the diagnostic configuration specific to this problem.
function config_diagnostics(driver_config)
    interval = "10000steps"
    dgngrp = setup_atmos_default_diagnostics(interval, driver_config.name)
    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

function main()
    ClimateMachine.init()

    FT = Float64
    N = 4
    Ncellsx = 64
    Ncellsy = 64
    Ncellsz = 64
    Δx = FT(2 * pi / Ncellsx)
    Δy = Δx
    Δz = Δx
    resolution = (Δx, Δy, Δz)
    xmax = FT(pi)
    ymax = FT(pi)
    zmax = FT(pi)
    xmin = FT(-pi)
    ymin = FT(-pi)
    zmin = FT(-pi)
    t0 = FT(0)
    timeend = FT(0.1)
    CFL = FT(1.8)

    driver_config = config_greenvortex(
        FT,
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        xmin,
        ymin,
        zmin,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFL,
    )
    dgn_config = config_diagnostics(driver_config)

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do (init = false)
        Filters.apply!(solver_config.Q, 6, solver_config.dg.grid, TMARFilter())
        nothing
    end
    # information.
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        check_euclidean_distance = true,
    )

end

main()


text:
using Makie, GLMakie, Observables, Colors, Yao
AbstractPlotting.inline!(true)

function bloch_arrow(θ, ϕ)
    x = sin(2θ) * cos(ϕ)
    y = sin(2θ) * sin(ϕ)
    z = cos(2θ)
    return (x, y, z)
end

function bloch_arrow(r::ArrayReg)
    @assert nqubits(r) == 1 "invalid quantum state, expect only one qubit"
    st = statevec(r)
    global_phase = exp(im * angle(st[1]))
    st = st ./ global_phase
    θ = acos(real(st[1])) # st[1] is real
    ϕ = iszero(θ) ? zero(θ) : angle(st[2] / sin(θ))
    return bloch_arrow(θ, ϕ)
end

function plot_arrow!(scene, r::ArrayReg; kwargs...)
    x, y, z = bloch_arrow(r)
    lines!(scene, [0.0, x], [0.0, y], [0.0, z]; kwargs...)
    return scene
end

function plot_block!(scene, x::AbstractBlock; linewidth=5.0, kwargs...)
    r = ArrayReg(bit"0")
    plot_arrow!(scene, r; color=:red, linewidth=linewidth, kwargs...)
    rand_color = RGB(rand(0:255)/255, rand(0:255)/255, rand(0:255)/255)
    plot_arrow!(scene, apply!(r, x); color=rand_color ,linewidth=linewidth, kwargs...)
    return scene
end

function bloch_sphere!(scene)
    wireframe!(scene, Sphere(Point3f0(0), 1f0), color=:cyan3, linewidth=0.1, thickness = 0.6f0, transparency = true)
    return scene
end

function plotapply(x::AbstractBlock)
    AbstractPlotting.inline!(false)
    scene = Scene(resolution=(1000, 1000), center=false)
    bloch_sphere!(scene)
    plot_block!(scene, x)
    rotate_cam!(scene, 0.5, 0.0, 0.0)
    return scene
end

function animate_rot()
    scene = Scene(resolution=(1000, 1000))
    bloch_sphere!(scene)
    scene.center = false

    time = Node(0.0)
    scene = lift()
    bloch_arrow!(scene, x)
    rotate_cam!(scene, 0.5, 0.0, 0.0)
end

# Pauli Gates
function animate_rot(blk::AbstractBlock)
    AbstractPlotting.inline!(false)
    screen = GLMakie.Screen()
    GLMakie.GLFW.set_visibility!(GLMakie.to_native(screen), AbstractPlotting.use_display[])

    AbstractPlotting.inline!(true)
    scene = Scene(resolution=(1000, 1000), center=false)

    bloch_sphere!(scene)

    r = ArrayReg(bit"0")
    plot_arrow!(scene, r; color=:red, linewidth=5.0)
    time = Node(0.0)

    lifted = lift(time) do t
        x, y, z = bloch_arrow(apply!(r, rot(blk, t)))
        [0.0, x], [0.0, y], [0.0, z]
    end

    lines!(scene, lifted; color=:blue, linewidth=5.0)
    rotate_cam!(scene, 0.1, 0.0, 0.0)

    for _ in 1:100
        sleep(0.1)
        push!(time, to_value(time) + 0.001)
        display(screen, scene)
    end
    return scene
end

# 1. plot pauli gates
plotapply(X)

# 2. plot rotate animation
scene = animate_rot(chain(X, H, X))

# 3. make video
scene = Scene(resolution=(1000, 1000), center=false)
bloch_sphere!(scene)
r = ArrayReg(bit"0")
plot_arrow!(scene, r; color=:red, linewidth=5.0)

T = Node(0.0)
lifted = lift(T) do t
    x, y, z = bloch_arrow(apply!(r, rot(X, t)))
    [0.0, x], [0.0, y], [0.0, z]
end

lines!(scene, lifted; color=:blue, linewidth=5.0)
rotate_cam!(scene, 0.1, 0.0, 0.0)

record(scene, "bloch.gif", 1:100) do i
    push!(T, 0.001 * i)
end


text:
# modified from Amanda's legacy parsing code; see more notes in wntr_dict_parser.jl

"""
Convert dictionary of water network to WaterSystems structure
"""
function dict_to_struct(data::Dict{String,Any})
    junctions = junction_to_struct(data["Junction"])
    j_dict = Dict{String, Junction}(junction.name => junction for junction in junctions)
    arcs = link_to_struct(data["Link"], j_dict)
    a_dict = Dict{String, Arc}(arc.name => arc for arc in arcs)
    res = res_to_struct(data["Reservoir"], j_dict)
    demands = demand_to_struct(data["Demand"], j_dict)
    tanks = tank_to_struct(data["Tank"], j_dict)
    patterns = pattern_to_struct(data["Pattern"])
    curves = curve_to_struct(data["Curve"])
    pipes = pipe_to_struct(data["Pipe"], a_dict)
    pumps = pump_to_struct(data["Pump"], a_dict, curves, patterns)
    valves = valve_to_struct(data["Valve"], a_dict)

    return junctions, arcs, res, demands, tanks, patterns, curves, pipes, pumps, valves
end

"""
Create array of junctions using WaterSystems.Junction type
"""
function junction_to_struct(j_dict::Dict{String, Any})
    junctions = Vector{Junction}(undef, length(j_dict))
    for (i, (name, junction)) in enumerate(j_dict)
        junctions[i] = Junction(name, junction["elevation"], junction["press_head"],
                                junction["coordinates"], nothing)
    end
    return junctions
end

"""
Create array of arcs using WaterSystems.Arc type
"""
function link_to_struct(l_dict::Dict{String, Any}, j_dict::Dict{String,Junction})
    arcs = Vector{Arc}(undef, length(l_dict))
    for (i, (name,link)) in enumerate(l_dict)
        from_name = link["connectionpoints"][:from]["name"]
        to_name = link["connectionpoints"][:to]["name"]
        arcs[i] = Arc(name, j_dict[from_name], j_dict[to_name])
    end
    return arcs
end

"""
Create array of reservoirs using WaterSystems.Reservoir type.
"""
function res_to_struct(r_vec::Vector{Any}, j_dict::Dict{String,Junction})
    reservoirs = Vector{Reservoir}(undef,length(r_vec))
    for (i, reservoir) in enumerate(r_vec)
        name = reservoir["name"]
        junction = j_dict[name]
        # assign the reservoir head pattern to the junction
        junction.pattern_name = reservoir["pattern_name"]
        reservoirs[i] = Reservoir(name, true, junction)
    end
    return reservoirs
end

"""
Create array of tanks using WaterSystems.Tank types. Only cylindrical tanks are
currently supported.
"""
function tank_to_struct(t_vec::Vector{Any}, j_dict::Dict{String,Junction})
    tanks = Vector{Tank}(undef, length(t_vec))
    for (i, tank) in enumerate(t_vec)
        name = tank["name"]
        tanks[i] = CylindricalTank(name, true, j_dict[name], tank["diameter"], tank["level"],
                                   tank["levellimits"])
    end
    return tanks
end

"""
Create array of demands using WaterSystems.WaterDemand types. Only static demands are
currently supported.
"""
function demand_to_struct(d_vec::Vector{Any}, j_dict::Dict{String,Junction})
    # Only populating base_demand and the pattern name; the actual forecast will be created
    # and added separately, following the paradigm of PowerSystems.jl. Note that the
    # base_demand and pattern values can vary wildly with no consistent scaling. Would it be
    # useful (considering the computational cost) to rescale all of the base demands and the
    # patterns? See "demand_dict!" in wntr_dict_parser.jl, JJS 12/26/19
    demands = Vector{WaterDemand}(undef, length(d_vec))
    for (i,  demand) in enumerate(d_vec)
        name = demand["name"]
        base_demand = demand["base_demand"]
        pattern_name = demand["pattern_name"]
        demands[i] = StaticDemand(name, true, j_dict[name], base_demand, pattern_name)
    end
    return demands
end

"""
Create array of patterns using WaterSystems.Pattern type.
"""
function pattern_to_struct(p_vec::Vector{Any})
    patterns = Vector{Pattern}(undef, length(p_vec))
    for (i,curve) in enumerate(p_vec)
        patterns[i] = Pattern(curve["name"], curve["multipliers"])
    end
    return patterns
end

"""
Create array of curves using WaterSystems.Curve type.
"""
function curve_to_struct(c_vec::Vector{Any})
    curves = Vector{Curve}(undef, length(c_vec))
    for (i,curve) in enumerate(c_vec)
        curves[i] = Curve(curve["name"], curve["type"], curve["points"])
    end
    return curves
end

"""
Create array of pipes using WaterSystems.Pipe types.
"""
function pipe_to_struct(pi_vec::Vector{Any}, a_dict::Dict{String,Arc})
    pipes = Vector{Pipe}(undef, length(pi_vec))
    # define arbitrary flow limits -- just remove flow limits from pipes? JJS 1/28/19
    #unid_flowlimits = CVPipe(nothing).flowlimits # neat trick to get defaults, not using
    #bidi_flowlimits = GatePipe(nothing).flowlimits
    unid_flowlimits = (min = small, max = large)
    bidi_flowlimits = (min = -large, max = large)
    for (i, pipe) in enumerate(pi_vec)
        name = pipe["name"]
        if pipe["cv"] == true
            pipes[i] = CVPipe(name, a_dict[name], true, pipe["diameter"], pipe["length"],
                              pipe["roughness"], unid_flowlimits, nothing, nothing)
        elseif pipe["control_pipe"] == true
            # add initial status (open/closed)? would need to add the field, JJS 12/26/19
            pipes[i] = GatePipe(name, a_dict[name], true, pipe["diameter"], pipe["length"],
                               pipe["roughness"], bidi_flowlimits, nothing, nothing, nothing)
        else
            pipes[i] = OpenPipe(name, a_dict[name], true, pipe["diameter"], pipe["length"],
                                pipe["roughness"], bidi_flowlimits, nothing, nothing)
        end
    end
    return pipes
end

"""
Create array of pumps using WaterSystem.Pump type.
"""
function pump_to_struct(pu_vec::Vector{Any}, a_dict::Dict{String,Arc}, c_vec::Vector{Curve},
                        pa_vec::Vector{Pattern})
    pumps = Vector{Pump}(undef, length(pu_vec))
    c_dict = Dict{String, Curve}(curve.name => curve for curve in c_vec)
    pa_dict = Dict{String, Pattern}(patt.name => patt for patt in pa_vec)
    for (i, pump) in enumerate(pu_vec)
        name = pump["name"]
        # creat EPANETPumpParams object for the pump
        if pump["efficiency"] isa String
            efficiency = c_dict[pump["efficiency"]]
        else
            efficiency = pump["efficiency"]
        end
        epanet_params = EPANETPumpParams(pump["type"], pump["power"],
                                        c_dict[pump["head_curve_name"]],
                                        efficiency)
        # calculate normalized params and create the object
         # in utils/PumpCoefs.jl
        Qbep, etabep, Gbep, Pbep = norm_pump_params(pump, c_dict)
        flowlimits = (min = 0.5*Qbep, max = 1.5*Qbep)
        norm_coefs = NormPumpParams(Qbep, etabep, Gbep, Pbep)
        # create the pump parameters object, including base price and pattern
        pump_params = PumpParams(epanet_params, norm_coefs, pump["price"],
                                 pump["price_pattern"])
        # create the pump object
        pumps[i] = Pump(name, a_dict[name], true, pump_params, flowlimits, nothing, nothing,
                        nothing)
    end
    return pumps
end

#ignoring valves for now, JJS 12/5/19
function valve_to_struct(v_vec::Vector{Any}, a_dict::Dict{String,Arc})
    if size(v_vec)[1] != 0
        @warn("There appears to be pressure/flow control valves in this network, but these valves are not currently handled by WaterSystems.jl")
    end
    return nothing
end
# function valve_to_struct(data::Vector{Any}, junctions::Array{Junction,1})
#     valves = Array{PR where {PR <: PressureReducingValve},1}(undef, length(data))
#     for (ix, v) in enumerate(data)
#         #push!(valves, PressureReducingValve(...))
#         j_from = v["connectionpoints"].from
#         j_to = v["connectionpoints"].to
#         junction_to = [junc for junc in junctions if junc.name == j_to["name"]][1]
#         junction_from = [junc for junc in junctions if junc.name == j_from["name"]][1]
#         if v["valvetype"] == "PRV"
#             valves[ix] = PressureReducingValve(v["name"], (from = junction_from, to = junction_to) ,v["status"], v["diameter"], v["pressure_drop"])
#         elseif v["valvetype"] == "PSV"
#             valves[ix] = PressureSustainingValve(v["name"], (from = junction_from, to = junction_to) ,v["status"], v["diameter"], v["pressure_drop"])
#         elseif v["valvetype"] == "PBV"
#             valves[ix] = PressureBreakerValve(v["name"], (from = junction_from, to = junction_to) ,v["status"], v["diameter"], v["pressure_drop"])
#         elseif v["valvetype"] == "FCV"
#             valves[ix] = FlowControlValve(v["name"], (from = junction_from, to = junction_to) ,v["status"], v["diameter"], v["pressure_drop"])
#         elseif v["valvetype"] == "TCV"
#             valves[ix] = ThrottleControlValve(v["name"], (from = junction_from, to = junction_to) ,v["status"], v["diameter"], v["pressure_drop"])
#         else
#             valves[ix] = GeneralPurposeValve(v["name"], (from = junction_from, to = junction_to) ,v["status"], v["diameter"], v["pressure_drop"])
#         end
#     end
#     return valves
# end


text:
"""
  IsingJ

  ver 0.2.2 Released Sep 2020
  This code is distributed under the constitution
  of MIT LICENSE.

  module Initial

  2019/04/24  Released by PearCandy                 """

module Smoothing_mag
    using Printf
    export output_T_vs_mag

 """  plot_T_vs_mag(_temp,_mag)
      output T vs Magnetization
      _temp: temperature
      _mag : absolute values of magnetization  """

    function output_T_vs_mag(low_t,dt,_mag)
        open("t_mag.dat","w") do io
            println(io,"# temperature vs ave.magnetization")

            for i in 1:length(_mag)
                _temp = low_t + dt*float(i)

                # smoothing of magnetization
                if(3<i<length(_mag)-3)
                    mag_ave=(_mag[i-2]+_mag[i-1]+_mag[i]+_mag[i+1]+_mag[i+2])/5
                else
                    mag_ave=_mag[i]
                end
                # -- output --
                @printf(io, " %+06f",     _temp  )
                @printf(io, " %+06f\n",   mag_ave)
            end
        end
    end

end #module


text:
@testset "Deprecations" begin
    # RodriguesVec
    @test Base.isdeprecated(Rotations, :RodriguesVec)
    @test Base.isdeprecated(Rotations, :Quat)
    @test Base.isdeprecated(Rotations, :SPQuat)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module TableTransforms

using Tables
using ScientificTypes
using Distributions: Normal
using Transducers: tcollect
using LinearAlgebra
using Statistics

import Distributions: ContinuousUnivariateDistribution
import Distributions: quantile, cdf

include("assertions.jl")
include("distributions.jl")
include("transforms.jl")

export
  # interface
  Transform,
  Stateless,
  Colwise,
  assertions,
  isrevertible,
  apply, revert, reapply,
  colapply, colrevert,

  # built-in
  Select,
  Reject,
  Rename,
  Identity,
  Center,
  Scale,
  MinMax,
  Interquartile,
  ZScore,
  Quantile,
  Functional,
  EigenAnalysis,
  PCA, DRS, SDS,
  Sequential,
  Parallel,
  →, ⊔

end


text:
using MultivariateStats
using Base.Test

srand(34568)

## prepare data

d = 5
ns = [10, 15, 20]
nc = length(ns)
n = sum(ns)
Xs = Matrix{Float64}[]
ys = Vector{Int}[]
Ss = Matrix{Float64}[]
cmeans = zeros(d, nc)

for k = 1:nc
	R = qr(randn(d, d))[1]
	nk = ns[k]

	Xk = R * scale(2 * rand(d) + 0.5, randn(d, nk)) .+ randn(d)
	yk = fill(k, nk)
	uk = vec(mean(Xk, 2))
	Zk = Xk .- uk
	Sk = Zk * Zk'

	push!(Xs, Xk)
	push!(ys, yk)
	push!(Ss, Sk)
	cmeans[:,k] = uk
end

X = hcat(Xs...)
y = vcat(ys...)
mv = vec(mean(X, 2))
Sw = zeros(d, d)
for k = 1:nc
	Sw += Ss[k]
end

Sb = zeros(d, d)
for k = 1:nc
	dv = cmeans[:,k] - mv
	Sb += ns[k] * (dv * dv')
end

Z = X .- mv
Sa = Z * Z'
@test_approx_eq Sw + Sb Sa

@assert size(X) == (d, n)
@assert size(y) == (n,)

## Stats

S = multiclass_lda_stats(nc, X, y)

@test S.dim == d
@test S.nclasses == nc

@test classweights(S) == ns
@test_approx_eq classmeans(S) cmeans
@test_approx_eq mean(S) mv

@test_approx_eq withclass_scatter(S) Sw
@test_approx_eq betweenclass_scatter(S) Sb

## Solve

emax = maximum(eigvals(Sw))
Sw_c = copy(Sw)
Sb_c = copy(Sb)

lambda = 1.0e-3
Sw_r = Sw + (lambda * emax) * eye(d)

P1 = mclda_solve(Sb, Sw, :gevd, nc-1, lambda)
@test_approx_eq P1' * Sw_r * P1 eye(nc-1)
U = Sb * P1
V = Sw_r * P1
# test whether U is proportional to V,
# which indicates P is the generalized eigenvectors
@test_approx_eq U scale(V, vec(mean(U ./ V, 1)))

P2 = mclda_solve(Sb, Sw, :whiten, nc-1, lambda)
@test_approx_eq P2' * Sw_r * P2 eye(nc-1)

@test_approx_eq P1 P2

## LDA

M = fit(MulticlassLDA, nc, X, y; method=:gevd, regcoef=lambda)
@test indim(M) == d
@test outdim(M) == nc - 1
@test_approx_eq projection(M) P1
@test_approx_eq M.pmeans M.proj'cmeans
@test_approx_eq transform(M, X) M.proj'X

M = fit(MulticlassLDA, nc, X, y; method=:whiten, regcoef=lambda)
@test indim(M) == d
@test outdim(M) == nc - 1
# @test_approx_eq projection(M) P2  # signs may change
@test_approx_eq M.pmeans M.proj'cmeans
@test_approx_eq transform(M, X) M.proj'X


text:
# Get the promoted type of the Index objects in a collection
# of Index (Tuple, Vector, ITensor, etc.)
indtype(i::Index) = typeof(i)
indtype(T::Type{<:Index}) = T
indtype(is::Tuple{Vararg{<:Index}}) = eltype(is)
indtype(is::Vector{<:Index}) = eltype(is)
indtype(A::ITensor...) = indtype(inds.(A))

indtype(tn1, tn2) = promote_type(indtype(tn1), indtype(tn2))
indtype(tn) = mapreduce(indtype, promote_type, tn)

# More general siteind that allows specifying
# the space
function _siteind(site_tag, n::Int; space)
  return addtags(Index(space, "Site,n=$n"), site_tag)
end

_siteinds(site_tag, N::Int; space) = __siteinds(site_tag, N, space)

function __siteinds(site_tag, N::Int, space::Vector)
  return [_siteind(site_tag, n; space=space[n]) for n in 1:N]
end

function __siteinds(site_tag, N::Int, space)
  return [_siteind(site_tag, n; space=space) for n in 1:N]
end

infsiteinds(s::Vector{<:Index}) = CelledVector(addtags(s, celltags(1)))
function infsiteinds(site_tag, N::Int; space=nothing, kwargs...)
  if !isnothing(space)
    s = _siteinds(site_tag, N; space=space)
  else
    # TODO: add a shift option
    s = siteinds(site_tag, N; kwargs...)
  end
  return infsiteinds(s)
end

function ITensors.linkinds(ψ::InfiniteMPS)
  N = nsites(ψ)
  return CelledVector([linkinds(ψ, (n, n + 1)) for n in 1:N])
end

function InfMPS(s::Vector, f::Function)
  return InfMPS(infsiteinds(s), f)
end

function indval(iv::Pair)
  return ind(iv) => val(iv)
end

function zero_qn(i::Index)
  return zero(qn(first(space(i))))
end

function insert_linkinds!(A; left_dir=ITensors.Out)
  # TODO: use `celllength` here
  N = nsites(A)
  l = CelledVector{indtype(A)}(undef, N)
  n = N
  s = siteind(A, 1)
  dim = if hasqns(s)
    kwargs = (; dir=left_dir)
    qn_ln = zero_qn(s)
    [qn_ln => 1] #Default to 0 on the right
  else
    kwargs = ()
    1
  end
  l[N] = Index(dim, default_link_tags("l", n, 1); kwargs...)
  for n in 1:(N - 1)
    # TODO: is this correct?
    dim = if hasqns(s)
      qn_ln = flux(A[n]) * left_dir + qn_ln#Fixed a bug on flux conservation
      [qn_ln => 1]
    else
      1
    end
    l[n] = Index(dim, default_link_tags("l", n, 1); kwargs...)
  end
  for n in 1:N
    A[n] = A[n] * onehot(l[n - 1] => 1) * onehot(dag(l[n]) => 1)
  end

  @assert all(i -> flux(i) == zero_qn(s), A) "Flux not invariant under one unit cell translation, not implemented"

  return A
end

function UniformMPS(s::CelledVector, f::Function; left_dir=ITensors.Out)
  sᶜ¹ = s[Cell(1)]
  A = InfiniteMPS([ITensor(sⁿ) for sⁿ in sᶜ¹])
  N = length(sᶜ¹)
  for n in 1:N
    Aⁿ = A[n]
    Aⁿ[indval(s[n] => f(n))] = 1.0
    A[n] = Aⁿ
  end
  insert_linkinds!(A; left_dir=left_dir)
  return A
end

function InfMPS(s::CelledVector, f::Function)
  # TODO: rename cell_length
  N = length(s)
  ψL = UniformMPS(s, f; left_dir=ITensors.Out)
  ψR = UniformMPS(s, f; left_dir=ITensors.In)
  ψC = InfiniteMPS(N)
  l = linkinds(ψL)
  r = linkinds(ψR)
  for n in 1:N
    ψCₙ = ITensor(dag(l[n])..., r[n]...)
    ψCₙ[l[n]... => 1, r[n]... => 1] = 1.0
    ψC[n] = ψCₙ
  end
  return ψ = InfiniteCanonicalMPS(ψL, ψC, ψR)
end

function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n)
  s = siteinds(only, ψ.AL)
  return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[]
end

function ITensors.expect(ψ::InfiniteCanonicalMPS, h::MPO)
  l = linkinds(ITensorInfiniteMPS.only, ψ.AL)
  r = linkinds(ITensorInfiniteMPS.only, ψ.AR)
  s = siteinds(ITensorInfiniteMPS.only, ψ)
  δˢ(n) = ITensorInfiniteMPS.δ(dag(s[n]), prime(s[n]))
  δˡ(n) = ITensorInfiniteMPS.δ(l[n], prime(dag(l[n])))
  δʳ(n) = ITensorInfiniteMPS.δ(dag(r[n]), prime(r[n]))
  ψ′ = prime(dag(ψ))

  ns = ITensorInfiniteMPS.findsites(ψ, h)
  nrange = ns[end] - ns[1] + 1
  idx = 2
  temp_O = δˡ(ns[1] - 1) * ψ.AL[ns[1]] * ψ′.AL[ns[1]] * h[1]
  for n in (ns[1] + 1):(ns[1] + nrange - 1)
    if n == ns[idx]
      temp_O = temp_O * ψ.AL[n] * ψ′.AL[n] * h[idx]
      idx += 1
    else
      temp_O = temp_O * ψ.AL[n] * δˢ(n) * ψ′.AL[n]
    end
  end
  temp_O = temp_O * ψ.C[ns[end]] * δʳ(ns[end]) * ψ′.C[ns[end]]
  return temp_O[]
end

function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteSum)
  return [expect(ψ, h[j]) for j in 1:nsites(ψ)]
end


text:
module TestC

import ..TestB.c as y
import package_example.TestD:d as z

function print_vars()
    @show y
    @show z
end

end


text:
#==============================================================================#
# ImportExport.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/importexport-2010-06-01.normal.json
#==============================================================================#

__precompile__()

module ImportExport

using AWSCore

"""
    using AWSSDK.ImportExport.cancel_job
    cancel_job([::AWSConfig], arguments::Dict)
    cancel_job([::AWSConfig]; JobId=, <keyword arguments>)

    using AWSCore.Services.importexport
    importexport([::AWSConfig], "POST", "/?Operation=CancelJob", "CancelJob", arguments::Dict)
    importexport([::AWSConfig], "POST", "/?Operation=CancelJob", "CancelJob", JobId=, <keyword arguments>)

# CancelJob Operation

This operation cancels a specified job. Only the job owner can cancel it. The operation fails if the job has already started or is complete.

# Arguments

## `JobId = ::String` -- *Required*

## `APIVersion = ::String`

# Returns

`CancelJobOutput`

# Exceptions

`InvalidJobIdException`, `ExpiredJobIdException`, `CanceledJobIdException`, `UnableToCancelJobIdException`, `InvalidAccessKeyIdException` or `InvalidVersionException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/importexport-2010-06-01/CancelJob)
"""
@inline cancel_job(aws::AWSConfig=default_aws_config(); args...) = cancel_job(aws, args)

@inline cancel_job(aws::AWSConfig, args) = AWSCore.Services.importexport(aws, "POST", "/?Operation=CancelJob", "CancelJob", args)

@inline cancel_job(args) = cancel_job(default_aws_config(), args)

"""
    using AWSSDK.ImportExport.create_job
    create_job([::AWSConfig], arguments::Dict)
    create_job([::AWSConfig]; JobType=, Manifest=, ValidateOnly=, <keyword arguments>)

    using AWSCore.Services.importexport
    importexport([::AWSConfig], "POST", "/?Operation=CreateJob", "CreateJob", arguments::Dict)
    importexport([::AWSConfig], "POST", "/?Operation=CreateJob", "CreateJob", JobType=, Manifest=, ValidateOnly=, <keyword arguments>)

# CreateJob Operation

This operation initiates the process of scheduling an upload or download of your data. You include in the request a manifest that describes the data transfer specifics. The response to the request includes a job ID, which you can use in other operations, a signature that you use to identify your storage device, and the address where you should ship your storage device.

# Arguments

## `JobType = "Import" or "Export"` -- *Required*

## `Manifest = ::String` -- *Required*

## `ManifestAddendum = ::String`

## `ValidateOnly = ::Bool` -- *Required*

## `APIVersion = ::String`

# Returns

`CreateJobOutput`

# Exceptions

`MissingParameterException`, `InvalidParameterException`, `InvalidAccessKeyIdException`, `InvalidAddressException`, `InvalidManifestFieldException`, `MissingManifestFieldException`, `NoSuchBucketException`, `MissingCustomsException`, `InvalidCustomsException`, `InvalidFileSystemException`, `MultipleRegionsException`, `BucketPermissionException`, `MalformedManifestException`, `CreateJobQuotaExceededException`, `InvalidJobIdException` or `InvalidVersionException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/importexport-2010-06-01/CreateJob)
"""
@inline create_job(aws::AWSConfig=default_aws_config(); args...) = create_job(aws, args)

@inline create_job(aws::AWSConfig, args) = AWSCore.Services.importexport(aws, "POST", "/?Operation=CreateJob", "CreateJob", args)

@inline create_job(args) = create_job(default_aws_config(), args)

"""
    using AWSSDK.ImportExport.get_shipping_label
    get_shipping_label([::AWSConfig], arguments::Dict)
    get_shipping_label([::AWSConfig]; jobIds=, <keyword arguments>)

    using AWSCore.Services.importexport
    importexport([::AWSConfig], "POST", "/?Operation=GetShippingLabel", "GetShippingLabel", arguments::Dict)
    importexport([::AWSConfig], "POST", "/?Operation=GetShippingLabel", "GetShippingLabel", jobIds=, <keyword arguments>)

# GetShippingLabel Operation

This operation generates a pre-paid UPS shipping label that you will use to ship your device to AWS for processing.

# Arguments

## `jobIds = [::String, ...]` -- *Required*

## `name = ::String`

## `company = ::String`

## `phoneNumber = ::String`

## `country = ::String`

## `stateOrProvince = ::String`

## `city = ::String`

## `postalCode = ::String`

## `street1 = ::String`

## `street2 = ::String`

## `street3 = ::String`

## `APIVersion = ::String`

# Returns

`GetShippingLabelOutput`

# Exceptions

`InvalidJobIdException`, `ExpiredJobIdException`, `CanceledJobIdException`, `InvalidAccessKeyIdException`, `InvalidAddressException`, `InvalidVersionException` or `InvalidParameterException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/importexport-2010-06-01/GetShippingLabel)
"""
@inline get_shipping_label(aws::AWSConfig=default_aws_config(); args...) = get_shipping_label(aws, args)

@inline get_shipping_label(aws::AWSConfig, args) = AWSCore.Services.importexport(aws, "POST", "/?Operation=GetShippingLabel", "GetShippingLabel", args)

@inline get_shipping_label(args) = get_shipping_label(default_aws_config(), args)

"""
    using AWSSDK.ImportExport.get_status
    get_status([::AWSConfig], arguments::Dict)
    get_status([::AWSConfig]; JobId=, <keyword arguments>)

    using AWSCore.Services.importexport
    importexport([::AWSConfig], "POST", "/?Operation=GetStatus", "GetStatus", arguments::Dict)
    importexport([::AWSConfig], "POST", "/?Operation=GetStatus", "GetStatus", JobId=, <keyword arguments>)

# GetStatus Operation

This operation returns information about a job, including where the job is in the processing pipeline, the status of the results, and the signature value associated with the job. You can only return information about jobs you own.

# Arguments

## `JobId = ::String` -- *Required*

## `APIVersion = ::String`

# Returns

`GetStatusOutput`

# Exceptions

`InvalidJobIdException`, `ExpiredJobIdException`, `CanceledJobIdException`, `InvalidAccessKeyIdException` or `InvalidVersionException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/importexport-2010-06-01/GetStatus)
"""
@inline get_status(aws::AWSConfig=default_aws_config(); args...) = get_status(aws, args)

@inline get_status(aws::AWSConfig, args) = AWSCore.Services.importexport(aws, "POST", "/?Operation=GetStatus", "GetStatus", args)

@inline get_status(args) = get_status(default_aws_config(), args)

"""
    using AWSSDK.ImportExport.list_jobs
    list_jobs([::AWSConfig], arguments::Dict)
    list_jobs([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.importexport
    importexport([::AWSConfig], "POST", "/?Operation=ListJobs", "ListJobs", arguments::Dict)
    importexport([::AWSConfig], "POST", "/?Operation=ListJobs", "ListJobs", <keyword arguments>)

# ListJobs Operation

This operation returns the jobs associated with the requester. AWS Import/Export lists the jobs in reverse chronological order based on the date of creation. For example if Job Test1 was created 2009Dec30 and Test2 was created 2010Feb05, the ListJobs operation would return Test2 followed by Test1.

# Arguments

## `MaxJobs = ::Int`

## `Marker = ::String`

## `APIVersion = ::String`

# Returns

`ListJobsOutput`

# Exceptions

`InvalidParameterException`, `InvalidAccessKeyIdException` or `InvalidVersionException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/importexport-2010-06-01/ListJobs)
"""
@inline list_jobs(aws::AWSConfig=default_aws_config(); args...) = list_jobs(aws, args)

@inline list_jobs(aws::AWSConfig, args) = AWSCore.Services.importexport(aws, "POST", "/?Operation=ListJobs", "ListJobs", args)

@inline list_jobs(args) = list_jobs(default_aws_config(), args)

"""
    using AWSSDK.ImportExport.update_job
    update_job([::AWSConfig], arguments::Dict)
    update_job([::AWSConfig]; JobId=, Manifest=, JobType=, ValidateOnly=, <keyword arguments>)

    using AWSCore.Services.importexport
    importexport([::AWSConfig], "POST", "/?Operation=UpdateJob", "UpdateJob", arguments::Dict)
    importexport([::AWSConfig], "POST", "/?Operation=UpdateJob", "UpdateJob", JobId=, Manifest=, JobType=, ValidateOnly=, <keyword arguments>)

# UpdateJob Operation

You use this operation to change the parameters specified in the original manifest file by supplying a new manifest file. The manifest file attached to this request replaces the original manifest file. You can only use the operation after a CreateJob request but before the data transfer starts and you can only use it on jobs you own.

# Arguments

## `JobId = ::String` -- *Required*

## `Manifest = ::String` -- *Required*

## `JobType = "Import" or "Export"` -- *Required*

## `ValidateOnly = ::Bool` -- *Required*

## `APIVersion = ::String`

# Returns

`UpdateJobOutput`

# Exceptions

`MissingParameterException`, `InvalidParameterException`, `InvalidAccessKeyIdException`, `InvalidAddressException`, `InvalidManifestFieldException`, `InvalidJobIdException`, `MissingManifestFieldException`, `NoSuchBucketException`, `ExpiredJobIdException`, `CanceledJobIdException`, `MissingCustomsException`, `InvalidCustomsException`, `InvalidFileSystemException`, `MultipleRegionsException`, `BucketPermissionException`, `MalformedManifestException`, `UnableToUpdateJobIdException` or `InvalidVersionException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/importexport-2010-06-01/UpdateJob)
"""
@inline update_job(aws::AWSConfig=default_aws_config(); args...) = update_job(aws, args)

@inline update_job(aws::AWSConfig, args) = AWSCore.Services.importexport(aws, "POST", "/?Operation=UpdateJob", "UpdateJob", args)

@inline update_job(args) = update_job(default_aws_config(), args)

end # module ImportExport

#==============================================================================#
# End of file
#==============================================================================#


text:
function __init__()
    # very simple warmup
    # XXX this doesn't look like it help time to first plot
    # because the slowest part is actually opening and activating the
    # JUNO panel...

    # t = @elapsed begin
    #     continuous_preview(false)
    #     f = Figure("_")
    #     plot([1, 2], [1, 2], ls="--", color="blue")
    #     plot!([1, 2], lw=0.05, label="blah")
    #     scatter([1, 2], mcol="red")
    #     if GP_ENV["HAS_BACKEND"]
    #         savefig(f, "_", res=100, format="png", path=GP_ENV["TMP_PATH"],_noout=true)
    #     end
    #     destroy(f)
    # end; println("Init done in $(round(t, digits=2))s.")

    # check if in IJulia, and if that's the case disable the continuous_preview by default
    isdefined(Main, :IJulia) && Main.IJulia.inited && continuous_preview(false)
    isdefined(Main, :Atom) && Main.Atom.PlotPaneEnabled.x && continuous_preview(true)
    return nothing
end


text:
@testset "flatten" begin
    @testset "Reals" begin
        test_flatten_interface(1.0)

        @testset "Integers" begin
            test_flatten_interface(1)
            @test isempty(first(flatten(1)))
        end
    end

    @testset "AbstractArrays" begin
        test_flatten_interface(randn(10))
        test_flatten_interface(randn(5, 4))
        test_flatten_interface([randn(5) for _ in 1:3])

        # Prevent regression of https://github.com/invenia/ParameterHandling.jl/issues/31
        @testset for v in [[1, 2, 3], sparse([1, 0, 3])]
            test_flatten_interface(v)
            @test length(first(flatten(v))) == 0
        end
    end

    @testset "SparseMatrixCSC" begin
        test_flatten_interface(sprand(10, 10, 0.5))
    end

    @testset "Tuple" begin
        test_flatten_interface((1.0, 2.0); check_inferred=tuple_infers)

        test_flatten_interface((1.0, (2.0, 3.0), randn(5)); check_inferred=tuple_infers)
    end

    @testset "NamedTuple" begin
        test_flatten_interface(
            (a=1.0, b=(2.0, 3.0), c=(e=5.0,)); check_inferred=tuple_infers
        )
    end

    @testset "Dict" begin
        test_flatten_interface(Dict(:a => 4.0, :b => 5.0); check_inferred=false)
    end
end


text:
module TestUsesThreads
import PerformanceTestTools
if Threads.nthreads() > 1
    include("__test_uses_threads.jl")
else
    PerformanceTestTools.@include_foreach(
        "__test_uses_threads.jl",
        [["JULIA_NUM_THREADS" => "2"]],
    )
end
end  # module


text:
using Documenter
using ModelicaScriptingTools

makedocs(
    sitename = "InaMo",
    format = Documenter.HTML(),
    pages = [
        "index.md",
        "Models" => [
            "models/cells.md",
            "models/membrane.md",
            "models/channels.md",
            "models/concentrations.md",
            "models/pumps.md",
            "models/experimental-methods.md",
            "models/examples.md"
        ],
        "Modelica features" => [
            "extends.md",
            "outer.md",
            "replaceable.md",
            "acausal.md",
        ],
        "unit-tests.md",
        "ci.md",
        "reconstruction.md"
    ],
    workdir = joinpath(@__DIR__, "build")
)

deploydocs(
    repo = "github.com/CSchoel/inamo.git",
    devbranch = "main",
    versions = ["v^", "v#.#", "stable" => "v^", "dev" => "dev"]
)

max_stars_repo_name: CSchoel/inamo

max_issues_repo_name: CSchoel/inamo

max_forks_repo_name: CSchoel/inamo


text:
# A Julia implementation of a Vigenere cipher

function vigenere()
    #Get plaintext and key
    println("Please enter the text: ")
    text = readline()
    println("Please enter the key: ")
    key = readline()
    #Key must be of equal length to plaintext
    while length(key) != length(text)
        println("Error: key must be the same length as the text: ", length(text))
        key = readline()
    end

    #The original cipher operated in uppercase, so I feel it is appropriate here
    text = uppercase(text)
    key = uppercase(key)

    #=Now we have a suitable text and key, time to encrypt. Index points to the current position in the Key
    A caesar cipher is carried out each time, where A corresponds to the current position in the Key
    Julia is 1 indexed! =#
    index = 1
    for char in text
        distance = key[index] - 'A'
        newchar = char + distance
        if(newchar > 'Z')
            newchar -= 26
        end
        if(newchar < 'A')
            newchar += 26
        end
        index += 1
        print(newchar)
    end

end
vigenere()


text:
using MOBC, vOptGeneric, GLPKMathProgInterface, PyPlot
using Base.Test

m = vModel(solver = GLPKSolverMIP())

nbint = 80
nbcont = 100
range_p = 50:100
range_w = 50:100

p1 = rand(range_p, nbint)
p2, p2cont = rand(range_p, nbint), rand(range_p, nbcont)
w, wcont = rand(range_w, nbint), rand(range_w, nbcont)
c = (sum(w) + sum(wcont)) ÷ 2

@variable(m, x[1:length(p1)], Bin)
@variable(m, 0 <= y[1:length(p2cont)] <= 1)
@addobjective(m, Max, dot(x, p1))
@addobjective(m, Max, dot(x, p2) + dot(y, p2cont))
@constraint(m, dot(x, w) + dot(y, wcont)<= c)

YN, XE = MOBC.solve_BC(m, 5000)

solve(m, method=:dichotomy)
voptYN=getY_N(m)

sleep(2)
clf()
plot(map(x->x[1], voptYN), map(x->x[2], voptYN), "bo", markersize="3")
plot(map(x->x[1], YN), map(x->x[2], YN), "rx", markersize="3")


text:
# Optimized Product Quantization. Adapted from Mohammad Norouzi's code.

using Distances
using Distributions # for sampling

include("kmeans.jl")
include("../pq/PQ.jl")
include("../utils.jl")

function quantize_opq{T <: AbstractFloat}(
  X::Matrix{T}, # d-by-n matrix of data points to quantize
  R::Matrix{T}, # d-by-d matrix. Learned rotation for X
  C::Vector{Matrix{T}},  # m-long array. Each entry is a d-by-h codebooks
  V::Bool=false)

  # Apply rotation and quantize as in PQ
  return quantize_pq( R'*X, C, V )
end

function train_opq{T <: AbstractFloat}(
  X::Matrix{T},         # d-by-n matrix of data points to train on.
  m::Integer,           # number of codebooks
  h::Integer,           # number of entries per codebook
  niter::Integer,       # number of optimization iterations
  init::AbstractString, # how to initialize the optimization
  V::Bool=false )    # wheter to print progress

  d, n = size( X )

  C = Vector{Matrix{T}}(m) # codebooks

  obj = zeros( Float32, niter+1 )

  # Number of bits in the final codes.
  nbits = log2(h) * m
  CB    = zeros(T, size(X))

  if init == "natural" # Initialize R with identity
    R = eye(T, d)
  elseif init == "random"
    R, _, _ = svd( randn( T, d, d ))
  else
    error("Intialization $init unknown")
  end

  RX = R' * X # Rotate the data

  subdims = splitarray( 1:d, m )

  # Initialization sampling RX
  for i = 1:m
    perm = sample(1:n, h, replace=false)
    C[i] = RX[ subdims[i], perm ]
  end

  costs     = zeros(Float32, n)
  counts    = zeros(Int, h)
  cweights  = zeros(Float32, h)
  to_update = zeros(Bool, h)
  unused    = Int[]

  # Initialize the codes -- B
  B = zeros( Int64, n, m )
  for i=1:m
    dmat = pairwise( SqEuclidean(), C[i], RX[subdims[i],:] )
    dmat = convert(Array{Float32}, dmat)
    update_assignments!( dmat, true, view( B,:,i ), costs, counts, to_update, unused  )

    CB[subdims[i],:] = C[i][:, B[:,i]]
  end

  for iter=0:niter

    obj[iter+1] = sum( (R*CB - X).^2 ) ./ n
    @printf("%3d %e   \n", iter, obj[iter+1])

    # update R
    U, S, V = svd(X * CB', thin=true)
    R = U * V'

    # update R*X
    RX = R' * X

    for i=1:m
      # update C
      C[i] = update_centers!( RX[subdims[i],:], B[:,i], cweights, h )

      # update B
      dmat = pairwise( SqEuclidean(), C[i], RX[subdims[i],:] )
      dmat = convert(Array{Float32}, dmat)
      update_assignments!( dmat, false, view(B,:,i), costs, counts, to_update, unused  )

      # update D*B
      CB[ subdims[i], : ] = C[i][:, B[:, i]]
    end # for i=1:m
  end # for iter=0:niter

  C = convert( Vector{Matrix{T}}, C )
  return C, B', R, obj
end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
module ModifiedLatinHypercubeSampling
# Hess, Stephane, Kenneth E. Train and John W. Polak, (2006), "On the Use of a Modified Latin Hypercube
# Sampling (MLHS) Method in the Estimation of a Mixed Logit Model for Vehicle Choice,"
# Transportation Research Part B, Vol. 40, pp. 147 - 163.

using Random
export mlhs

# R draws for D dimensions
# R is the number of draws
# D is the number of dimensions

function mlhs(R, D)

    # ϕ from Equation (14)
    ϕ = (0:(R-1)) ./ R
    ϕ_grid = repeat(ϕ, 1, D)

    # Obtain random value
    Random.seed!(19900729)
    randval = rand(D)
    while 0.0 in randval || 1.0 in randval
        randval = rand(D)
    end
    x = transpose(randval / R)

    # Ψ from Equation (15)
    Ψ = ϕ .+ x

    # Shuffle sequences in each dimension
    Ψ_shuffled = zeros(R, D)

    for i in 1:D
        Ψ_shuffled[:,i] = Ψ[shuffle!(collect(1:R)),i]
    end
    Ψ_shuffled
end

end # Module


text:
"A module for function composition."
module Composite

import Base.call, Base.%, Base.<, Base.>, Base.^

export AbstractCallable, Callable,
    AbstractCompositionMethod,  SequentialComposition, ConcurrentComposition, CompositeFunction,
    call, ∘, %, <, >, ^, dup

abstract AbstractCallable
typealias Callable Union{AbstractCallable,Function}

abstract AbstractCompositionMethod
abstract SequentialComposition <: AbstractCompositionMethod
abstract ConcurrentComposition <: AbstractCompositionMethod

dup( x ) = ( x, x )

immutable CompositeFunction{ First <: Callable, Second <: Callable, Method <: AbstractCompositionMethod } <: AbstractCallable
    first::First
    second::Second
end

call{ First, Second, Xs }( f::CompositeFunction{First,Second,SequentialComposition }, xs::Xs ... ) = f.first( f.second( xs ... ) ... )
call{ First, Second, X, Xs }( f::CompositeFunction{ First,Second,ConcurrentComposition}, x::X, xs::Xs ... ) = tuple( f.first( x ), f.second( xs ... ) ... )

comp( f, g ) = CompositeFunction{typeof(f),typeof(g),SequentialComposition}( f, g )
comp( f, g, h, hs ... ) = comp( comp( f,  g ), h, hs ... )
∘{ F <: Callable, G <: Callable }( f::F, g::G, hs ... ) = comp( f, g, hs ... )

split( f, g ) = CompositeFunction{typeof(f),typeof(g),ConcurrentComposition}( f, g )
split( f, g, h, hs ... ) = split( split( f, g ), h, hs ... )
%{F<:Callable,G<:Callable}( f::F, g::G, hs ... ) = split( f, g, hs ... )

fan( f, g ) = (f % g ) ∘ dup
fan( f, g , h, hs ... ) = fan( f, fan( g, h, hs ... ))
^{F<:Callable,G<:Callable}( f::F, g::G, hs ... ) = fan( f, g, hs ... )

abstract AbstractPacking
abstract Left <: AbstractPacking
abstract Right <: AbstractPacking

immutable PartialApplication{F <: Callable, Packing <: AbstractPacking} <: AbstractCallable
    callable::F
    args::Tuple
end

call{F}( f::PartialApplication{F,Left}, xs ... ) = f.callable( f.args ..., xs ... )
call{F}( f::PartialApplication{F,Right}, xs ... ) = f.callable( xs ..., f.args ... )

partial{F <: Callable }( f::F, xs ... ) = PartialApplication{F,Left}( f, tuple( xs ... ))
rpartial{F <: Callable }( f::F, xs ... ) = PartialApplication{F,Right}( f, tuple( xs ... ))

<{F <: Callable}( f::F, xs ... ) = partial( f, xs ... )
>{F <: Callable}( f::F, xs ... ) = rpartial( f, xs ... )

end # module

max_stars_repo_name: sabjohnso/Composite

max_issues_repo_name: sabjohnso/Composite

max_forks_repo_name: sabjohnso/Composite


text:
using RecipesBase
using KernelDensity
using Random
using PlotUtils
using Colors

#=
!!!! IMPORTANT !!!!
Each code path of the @series macro must return the data for the series.
data = ... should be the last statement
=#
@recipe function abc_output_recipe(abco::ABCOutput;
      runs=nothing, params=nothing, population_colors=nothing)
      if params === nothing
            params = [i for i in 1:abco.n_params]
      end
      if runs === nothing
            runs = [i for i in 1:length(abco.population)]
      elseif isa(runs, Int)
            runs = [i for i in max(length(abco.population) - runs + 1, 1):length(abco.population)]
      end
      legend --> false
      layout := length(params) ^ 2
      for (i, par1) in enumerate(params)
            for (j, par2) in enumerate(params)
                  subplot := (i - 1) * length(params) + j
                  if i == j
                        @series begin
                              seriestype := :histogram
                              ylabel := "Accepted"
                              xlabel := "Parameter $(par1)"
                              bins --> 50
                              if isa(abco.population, Vector)
                                    data = abco.population[end][:, par1]
                              else
                                    data = abco.population[:, par1]
                              end # if Vector
                        end  # @series
                  elseif j < i
                        seriestype := :scatter
                        markerstrokecolor --> false
                        if isa(abco.population, Vector)
                              for r in runs
                                    @series begin
                                          xlabel --> "Parameter $(par2)"
                                          ylabel --> "Parameter $(par1)"
                                          if population_colors !== nothing
                                                idx = r % length(population_colors)
                                                if idx == 0
                                                      idx = length(population_colors)
                                                end
                                                seriescolor := population_colors[idx]
                                          end
                                          pop = abco.population[r]
                                          data = (pop[:, par2], pop[:, par1])
                                    end # @series
                              end # for r
                        else
                              @series begin
                                    data = (abco.population[:, par2], abco.population[:, par1])
                              end # @series
                        end # if Vector
                  else
                        @series begin
                              ylabel := ""
                              xlabel := ""
                              grid := false
                              xaxis := false
                              yaxis := false
                              data = []
                        end # @series
                  end # if/elseif
            end # for j, par2
      end # for i, par1
end # @recipe

function scale_outer_intervals(intervals, scale=1.0)
    int_min = min([int[1] for int in intervals]...)
    int_max = max([int[2] for int in intervals]...)
    int_scaled_half = (int_max - int_min) * scale / 2.0
    int_mid = (int_min + int_max) / 2.0
    int_mid - int_scaled_half, int_mid + int_scaled_half
end

my_linspace(bounds::Tuple{Float64, Float64}, length::Int64) = range(bounds[1], stop=bounds[2], length=length)

@recipe function emulation_vs_simulation(emu_out::EmulatedABCSMCOutput, sim_out::SimulatedABCSMCOutput, true_params::AbstractVector{<:Real}=Float64[];
            pop_idx=0,
            pop_color = "#007731",
            simulation_color = "#08519c",
            emulation_color = "#800013",
            contour_levels=8)
      if pop_idx == 0
            pop_idx = min(length(sim_out.population), length(emu_out.population))
      end
      kernel_bandwidth_scale = 0.09
      bounds_scale = 1.2
      sim_count = 30
      params = 1:emu_out.n_params
      legend --> false
      layout := length(params) ^ 2
      for i in params
            for j in params
                  subplot := (i - 1) * length(params) + j
                  if j < i
                        x_data_emu = emu_out.population[pop_idx][:,j]
                        y_data_emu = emu_out.population[pop_idx][:,i]
                        x_data_sim = sim_out.population[pop_idx][:,j]
                        y_data_sim = sim_out.population[pop_idx][:,i]
                        sim_size = size(sim_out.population[pop_idx], 1)
                        if sim_size > sim_count
                            idx = randperm(sim_size)[1:sim_count]
                            x_data_sim = x_data_sim[idx]
                            y_data_sim = y_data_sim[idx]
                        end
                        x_extr_emu = extrema(x_data_emu)
                        y_extr_emu = extrema(y_data_emu)
                        x_extr_sim = extrema(x_data_sim)
                        y_extr_sim = extrema(y_data_sim)

                        bandwidth = (-kernel_bandwidth_scale * -(x_extr_emu...), -kernel_bandwidth_scale * -(y_extr_emu...))
                        kde_joint = kde((x_data_emu, y_data_emu), bandwidth=bandwidth)
                        contour_x = my_linspace(scale_outer_intervals([x_extr_emu], bounds_scale), 100)
                        contour_y = my_linspace(scale_outer_intervals([y_extr_emu], bounds_scale), 100)
                        contour_z = pdf(kde_joint, contour_x, contour_y)
                        @series begin
                              seriestype := :contour
                              color := ColorGradient([:white, emulation_color])
                              levels := contour_levels
                              fill := true
                              (contour_x, contour_y, contour_z)
                        end
                        @series begin
                              seriestype := :scatter
                              markersize := 4
                              markershape := :x
                              color := simulation_color
                              (x_data_sim, y_data_sim)
                        end
                        if length(true_params) > 0
                              seriestype := :path
                              linestyle := :dash
                              color := :black
                              markershape := :none
                              @series begin
                                    ([true_params[j], true_params[j], minimum(contour_x)],
                                    [minimum(contour_y), true_params[i], true_params[i]])
                              end
                        end
                  elseif j > i
                        pop_colors = range(parse(Colorant, "white"), parse(Colorant, pop_color), length=pop_idx+1)
                        seriestype := :scatter
                        markershape := :circle
                        markersize := 2
                        for k in 1:pop_idx
                              @series begin
                                    markercolor := pop_colors[k+1]
                                    markerstrokecolor := pop_colors[k+1]
                                    (emu_out.population[k][:,j], emu_out.population[k][:,i])
                              end
                        end
                        if length(true_params) > 0
                              seriestype := :path
                              linestyle := :dash
                              color := :black
                              markershape := :none
                              @series begin
                                    ([true_params[j], true_params[j], minimum(emu_out.population[1][:,j])],
                                    [minimum(emu_out.population[1][:,i]), true_params[i], true_params[i]])
                              end
                        end
                  else # i == j
                        seriestype := :path
                        markershape := :none
                        emu_data = emu_out.population[pop_idx][:,i]
                        sim_data = sim_out.population[pop_idx][:,i]
                        extr_emu = extrema(emu_data)
                        extr_sim = extrema(sim_data)
                        kde_emu = kde(emu_data, bandwidth=-kernel_bandwidth_scale * -(extr_emu...))
                        kde_sim = kde(sim_data, bandwidth=-kernel_bandwidth_scale * -(extr_sim...))
                        x_bounds = scale_outer_intervals([extr_emu, extr_sim], bounds_scale)
                        x_plot = my_linspace(x_bounds, 100)
                        y_emu_plot = pdf(kde_emu,x_plot)
                        y_sim_plot = pdf(kde_sim,x_plot)
                        @series begin
                              linestyle := :solid
                              color := simulation_color
                              (x_plot, y_sim_plot)
                        end
                        @series begin
                              linestyle := :solid
                              color := emulation_color
                              (x_plot, y_emu_plot)
                        end
                        if length(true_params) > 0
                              true_param = true_params[i]
                              max_pdf = max(pdf(kde_emu, true_param), pdf(kde_sim, true_param))
                              @series begin
                                    linestyle := :dash
                                    color := :black
                                    ([true_param, true_param], [0, max_pdf])
                              end
                        end
                  end # if/else/elseif i/j
            end # for j
      end # for i
end

# Plot recipe for mdoel selection output
@recipe function modelselection_plotrecipe(::Type{ModelSelectionOutput}, mso::ModelSelectionOutput)
    seriestype := :line
    xlabel --> "Population"
    ylabel --> "Number of accepted particles"
    labels --> [string("Model ", m) for m in 1:mso.M]
    data = [[mso.n_accepted[i][j] for i in 1:size(mso.n_accepted,1)] for j in 1:mso.M]
end


text:
#=
   AR(1) in Turing.jl

   From Stan Users Guide, section 3.1 (page 51) "Autoregressive Models"
   (Though this is a port of my WebPPL model ar1.wppl)

   From the R model
   """
  Result:
          mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
   alpha  0.73    0.00 0.10  0.54  0.67  0.73  0.80  0.92 14499    1
   beta   0.16    0.00 0.10 -0.04  0.09  0.16  0.23  0.36 14601    1
   sigma  0.35    0.00 0.03  0.31  0.34  0.35  0.37  0.41 19812    1
   lp__  53.17    0.01 1.25 49.90 52.60 53.50 54.09 54.61 12256    1
   """

  Using Normal, Normal, and Normal (for alpha, beta, sigma)
  """
  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

       alpha    0.7332    0.0943     0.0009    0.0015   4013.6785    1.0001      388.6211
        beta    0.1554    0.1002     0.0010    0.0016   3916.5955    1.0001      379.2211
       sigma    0.3520    0.0259     0.0003    0.0004   5418.8333    1.0002      524.6740

  Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       alpha    0.5492    0.6695    0.7329    0.7962    0.9212
        beta   -0.0423    0.0878    0.1546    0.2228    0.3491
       sigma    0.3060    0.3339    0.3502    0.3683    0.4076
  """

  Using Normal, Normal, and Truncated(Gamma(2,2), 0,Inf) (for alpha, beta, sigma)
  """
  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

       alpha    0.7279    0.0977     0.0010    0.0015   3948.5885    1.0000      757.1598
        beta    0.1612    0.1038     0.0010    0.0016   4001.2383    1.0001      767.2557
       sigma    0.3532    0.0258     0.0003    0.0003   5222.2375    0.9999     1001.3878

  Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       alpha    0.5348    0.6627    0.7290    0.7928    0.9226
        beta   -0.0428    0.0929    0.1602    0.2296    0.3664
       sigma    0.3064    0.3351    0.3517    0.3693    0.4073

   """

   Cf ~/stan/ar1.R and ~/stan/ar1.stan
      ~/webppl/ar1.wppl

=#

using Turing
include("jl_utils.jl")

@model function ar1(ys)
    n = length(ys)

    # Uninformative priors
    # alpha ~ Uniform(-2,2)
    # beta ~ Uniform(-2,2)
    # sigma ~ Uniform(0,2)

    # Informative priors
    alpha ~ Normal(2,2)
    beta ~ Normal(2,2)
    sigma ~ Truncated(Gamma(2,2), 0, Inf)

    for i in 1:n
        if i == 1
            ys[i] ~ Normal(alpha + beta, sigma)
        else
            ys[i] ~ Normal(alpha + beta * ys[i-1], sigma)
        end
    end

end

ys = [0.705429, 1.43062, 0.618161, 0.315107, 1.09993, 1.49022, 0.690016, 0.587519, 0.882984,
      1.0278, 0.998615, 0.878366,  1.17405,  0.532718, 0.486417, 1.13685, 1.32453, 1.3661,
      0.914368, 1.07217, 1.1929, 0.418664, 0.889512, 1.47218, 1.13471, 0.410168, 0.639765,
      0.664874, 1.12101, 1.22703, -0.0931769, 0.4275, 0.901126, 1.01896, 1.27746, 1.17844,
      0.554775, 0.325423, 0.494777, 1.05813, 1.10177, 1.11225, 1.34575, 0.527594, 0.323462,
      0.435063, 0.739342, 1.05661, 1.42723, 0.810924, 0.0114801, 0.698537, 1.13063, 1.5286,
      0.968813, 0.360574, 0.959312, 1.2296, 0.994434, 0.59919, 0.565326, 0.855878, 0.892557,
      0.831705, 1.31114, 1.26013, 0.448281, 0.807847, 0.746235, 1.19471, 1.23253, 0.724155,
      1.1464, 0.969122, 0.431289, 1.03716, 0.798294, 0.94466, 1.29938, 1.03269, 0.273438,
      0.589431, 1.2741, 1.21863, 0.845632, 0.880577, 1.26184, 0.57157, 0.684231, 0.854955,
      0.664501, 0.968114, 0.472076, 0.532901, 1.4686, 1.0264, 0.27994, 0.592303, 0.828514,
      0.625841]

model = ar1(ys)
num_chains = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 10_000)
# chs = sample(model, IS(), 10_000)
# chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chs)
# chs = sample(model, SGLD(), 10_000)
chs = sample(model,NUTS(), 10_000)

display(chs)
# display(plot(chs))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank