text:
using Test, FileIO, JLD2
using Knet.Ops20: RNN

@testset "JLD" begin

    #needed for load macro test: https://github.com/simonster/JLD2.jl/blob/cc56a4d6da116d6172a4ea89f4bec9d17154a0ba/test/loadsave.jl#L5L17
    fn = tempname()*".jld2"
    save(fn,"model",RNN(1,1))

    @eval begin
        function macro_load()
            @load $fn
            model2 = RNN(1,1)
            return typeof(model)==typeof(model2)
         end
    end

    #@test macro_load()

    function macro_save()
        model = RNN(1,1)
        @save fn model
        true
     end

    @test macro_save()
    @test macro_load()

    function fun_sl()
        model = RNN(1,1)
        save(fn,"model",model)
        model2 = load(fn,"model")
        typeof(model)==typeof(model2)
    end

    @test fun_sl()

    function collections_sl()
        model  = RNN(1,1)
        model2 = RNN(1,1)
        save(fn,"model",[model,model2])
        models = load(fn,"model")
        test1 = typeof(models[1])==typeof(model)
        save(fn,"model",(model,model2))
        models = load(fn,"model")
        test2 = typeof(models[1])==typeof(model)
        save(fn,"model",Dict("model"=>model,"model2"=>model2))
        models = load(fn,"model")
        test3 = typeof(models["model"])==typeof(model)
        save(fn,Dict("model"=>model,"model2"=>model2))
        models = load(fn)
        test3 = typeof(models["model"])==typeof(model)
        test1 && test2 && test3
    end

    @test collections_sl()

    rm(fn)
end


text:
using Vann
using BlackBoxOptim

################################################################################

# Test calibration of Gr4j

# Read data

filename = joinpath(Pkg.dir("Vann"), "data", "airgr", "test_data.txt")

data = readdlm(filename, ',', header = true)

prec  = data[1][1:1000,1]
epot  = data[1][1:1000,2]
q_obs = data[1][1:1000,3]

prec = transpose(prec)
epot = transpose(epot)

frac = zeros(Float64, 1)

tstep = 24.0

time = DateTime(2000,1,1)

# Select model

st_hydro = Gr4j(tstep, time)

# Run calibration

param_opt = run_model_calib(st_hydro, prec, epot, q_obs; warmup = 1, verbose = :verbose)

# Parameter values

param_test  = [257.238, 1.012, 88.235, 2.208]

# Test whether calibration gives correct answer

param_range = get_param_range(st_hydro)

println("Finished test calibration of Gr4j")
println(param_test)
println(param_opt)

for iparam in eachindex(param_test)

  range_param = param_range[iparam][2] - param_range[iparam][1]

  if !(param_test[iparam]-0.05*range_param < param_opt[iparam] < param_test[iparam]+0.05*range_param)
    error("Calibration of gr4j resulted in wrong parameter values")
  end

end

################################################################################

# Test calibration of TinBasic and Gr4j

# Load data

filename = joinpath(Pkg.dir("Vann"), "data", "atnasjo")

date, tair, prec, q_obs, frac = load_data(filename, "Q_ref.txt")

tstep = 24.0

time = date[1]

# Crop data

date = date[1:1000]
tair = tair[:,1:1000]
prec = prec[:,1:1000]
q_obs = q_obs[1:1000]

# Compute potential evapotranspiration

epot = epot_zero(date)

# Select model

st_snow = TinBasic(tstep, time, frac)
st_hydro = Gr4j(tstep, time)

# Run calibration

param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs; warmup = 1, verbose = :verbose)

param_opt = vcat(param_snow, param_hydro)

# Parameter values

param_test  = [0.0, 3.69, 1.02, 74.59, 0.81, 214.98, 1.24]

# Test whether calibration gives correct answer

param_range_snow  = get_param_range(st_snow)
param_range_hydro = get_param_range(st_hydro)

param_range = vcat(param_range_snow, param_range_hydro)

println("Finished test calibration of TinBasic and Gr4j")
println(param_test)
println(param_opt)

for iparam in eachindex(param_test)

  range_param = param_range[iparam][2] - param_range[iparam][1];

  if !(param_test[iparam]-0.25*range_param < param_opt[iparam] < param_test[iparam]+0.25*range_param)
    error("Calibration of gr4j + snow resulted in wrong parameter values")
  end

end

# ################################################################################

# # Test calibration of TinStandard and Gr4j

# # Load data

# filename = joinpath(dirname(@__FILE__), "../data/atnasjo")

# date, tair, prec, q_obs, frac = load_data(filename, "Q_ref.txt")

# # Compute potential evapotranspiration

# epot = epot_zero(date)

# # Select model

# st_snow = TinStandard(frac)
# st_hydro = Gr4j()

# # Run calibration

# res = run_model_calib(st_snow, st_hydro, date, tair, prec, epot, q_obs)

# println("Finished test calibration of TinStandard and Gr4j")
# println(res)

ext: jl

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
function _process_raw_csv_api_data(api_call_raw_data::String)::(Union{T, Nothing} where T<:Any)

    # check: do we have an empty string?
    # check: legit string?

    # need to check to see if legit data is coming back from the service -
    if (is_valid_json(api_call_raw_data) == true)
        api_data_dictionary = JSON.parse(api_call_raw_data)
        if (haskey(api_+data_dictionary,"Error Message") == true)

            # grab the error mesage -
            error_message = api_data_dictionary["Error Message"]

            # throw -
            return Result{AVKError}(AVKError(error_message))
        else

            # formulate an error message -
            error_message = "Error: CSV type returns JSON without error message"

            # throw -
            return Result{AVKError}(AVKError(error_message))
        end
    end

    # create a data table from the CSV data -
    tmp_data_table = CSV.read(IOBuffer(api_call_raw_data))

    # sort the table according to the timestamps -
    idx_sort = sortperm(tmp_data_table[:,1])

    # create a sorted data table -
    sorted_data_table = tmp_data_table[idx_sort,:]

    # return the sorted table -
    return Return{DataFrame}(sorted_data_table)
end

function _process_raw_json_api_data_sts_adjusted(api_call_raw_data::String, data_series_key::String)::(Union{T, Nothing} where T<:Any)
    # is the data coming back well formed, and does it contain valid data?
    check_result = _check_json_api_return_data(api_call_raw_data)
    if check_result != nothing
        return check_result
    end

    # if we get here, we have valid JSON. Build dictionary -
    api_data_dictionary = JSON.parse(api_call_raw_data)

    # grab the time series data -
    time_series_key = data_series_key
    if (haskey(api_data_dictionary, time_series_key) == false)

        # throw an error -
        error_message = "Error: Missing the Time series key = $(time_series_key)"

        # throw -
        return AVKError(error_message)
    end
    #all repeated up to here possibly look to make the above a function for just checking data
    # array of keys -- this is changed per function
    data_key_label_array = ["1. open", "2. high", "3. low", "4. close", "5. adjusted close", "6. volume", "7. dividend amount"]
    number_of_fields = length(data_key_label_array)

    #initialize arrays to hold these fields for all dates
    timestamp_array = Dates.Date[]
    open_price_array = Float64[]
    high_price_array = Float64[]
    low_price_array = Float64[]
    close_price_array = Float64[]
    adjusted_close_array = Float64[]
    volume_array = Int64[]
    dividend_amount_array = Float64[]

    # ok, we have the time series key, go through the data and load into the table -
    time_series_dictionary = api_data_dictionary[time_series_key]
    time_series_key_array = collect(keys(time_series_dictionary))
    for timestamp_value in time_series_key_array

        # get the local_dictionary -
        local_dictionary = time_series_dictionary[timestamp_value]

        # cache -
        push!(timestamp_array, Dates.Date(timestamp_value,"yyyy-mm-dd"))

        # add the price data -
        for key_index = 1:number_of_fields

            # grab key -
            local_key = data_key_label_array[key_index]
            value = local_dictionary[local_key]

            # populate the array's -
            if (key_index == 1)
                push!(open_price_array, parse(Float64, value))
            elseif (key_index == 2)
                push!(high_price_array, parse(Float64, value))
            elseif (key_index == 3)
                push!(low_price_array, parse(Float64, value))
            elseif (key_index == 4)
                push!(close_price_array, parse(Float64, value))
            elseif (key_index == 5)
                push!(adjusted_close_array, parse(Float64, value))
            elseif (key_index == 6)
                push!(volume_array, parse(Int64, value))
            else
                push!(dividend_amount_array, parse(Float64, value))
            end
        end
    end
    # we need to sort the timestamps, to make them in reverse order -
    idx_sort = sortperm(timestamp_array)

    # build the data frame -
    data_frame = DataFrame(timestamp=timestamp_array[idx_sort], open=open_price_array[idx_sort], high=high_price_array[idx_sort], low=low_price_array[idx_sort], close=close_price_array[idx_sort], adjusted_close = adjusted_close_array[idx_sort], volume=volume_array[idx_sort], dividend_amount = dividend_amount_array[idx_sort])

    # return the data back to the caller -
    return Result{DataFrame}(data_frame)
end

function _process_raw_json_api_data_sts_endpoint(api_call_raw_data::String, data_series_key::String)::(Union{T, Nothing} where T<:Any)
end

function _process_raw_json_api_data_sts(api_call_raw_data::String, data_series_key::String)::(Union{T, Nothing} where T<:Any)

    # is the data coming back well formed, and does it contain valid data?
    check_result = _check_json_api_return_data(api_call_raw_data)
    if check_result != nothing
        return check_result
    end

    # TODO - check for missing data series key -
    # ...

    # if we get here, we have valid JSON. Build dictionary -
    api_data_dictionary = JSON.parse(api_call_raw_data)

    # grab the time series data -
    time_series_key = data_series_key
    if (haskey(api_data_dictionary, time_series_key) == false)

        # throw an error -
        error_message = "Error: Missing the Time series key = $(time_series_key)"

        # throw -
        return AVKError(error_message)
    end

    # array of keys -
    data_key_label_array = ["1. open", "2. high", "3. low", "4. close", "5. volume"]
    number_of_fields = length(data_key_label_array)

    # initialize storage for the fields -
    timestamp_array = Dates.Date[]
    open_price_array = Float64[]
    high_price_array = Float64[]
    low_price_array = Float64[]
    close_price_array = Float64[]
    volume_array = Int64[]

    # ok, we have the time series key, go through the data and load into the table -
    time_series_dictionary = api_data_dictionary[time_series_key]
    time_series_key_array = collect(keys(time_series_dictionary))
    for timestamp_value in time_series_key_array

        # get the local_dictionary -
        local_dictionary = time_series_dictionary[timestamp_value]

        # cache -
        push!(timestamp_array, Dates.Date(timestamp_value,"yyyy-mm-dd"))

        # add the price data -
        for key_index = 1:number_of_fields

            # grab key -
            local_key = data_key_label_array[key_index]
            value = local_dictionary[local_key]

            # populate the array's -
            if (key_index == 1)
                push!(open_price_array, parse(Float64, value))
            elseif (key_index == 2)
                push!(high_price_array, parse(Float64, value))
            elseif (key_index == 3)
                push!(low_price_array, parse(Float64, value))
            elseif (key_index == 4)
                push!(close_price_array, parse(Float64, value))
            else
                push!(volume_array, parse(Int64, value))
            end
        end
    end

    # we need to sort the timestamps, to make them in reverse order -
    idx_sort = sortperm(timestamp_array)

    # build the data frame -
    data_frame = DataFrame(timestamp=timestamp_array[idx_sort], open=open_price_array[idx_sort], high=high_price_array[idx_sort], low=low_price_array[idx_sort], close=close_price_array[idx_sort], volume=volume_array[idx_sort])

    # return the data back to the caller -
    return Result{DataFrame}(data_frame)

end

max_stars_repo_name: teamalphavantage/AlphaVantageAPI

max_issues_repo_name: teamalphavantage/AlphaVantageAPI

max_forks_repo_name: teamalphavantage/AlphaVantageAPI


text:
using Revise
using PyPlot
using Infiltrator
using LinearAlgebra
using Bem2d
using PyCall
using FEniCS

include("Helpers.jl")

mu = 30e9
lambda = 30e9
eta = 1e19 # just made up a number!
VE_z = -16e3 # Depth below which we use a Maxwell rheology
nu = 0.25
g = 9.81
rho = 2700
siay = 3600 * 24 * 365.2425

# Element geometries and data structures
elsbox = Elements(Int(1e5))

# Box geometry
B = -50e3 # Bottom
R = 20e3 # Right
T = 0e3 # Top
L = -20e3; # Left
nside = 50 # number of elements per box side

addelsez!(elsbox, discretizedline(L, B, R, B, nside)..., "B")
addelsez!(elsbox, discretizedline(R, B, R, T, nside)..., "R")
addelsez!(elsbox, discretizedline(R, T, L, T, nside)..., "T")
addelsez!(elsbox, discretizedline(L, T, L, B, nside)..., "L")

# Fault geometry
nfault = 1
x1, y1, x2, y2 = discretizedline(-10e3, -10e3, -5e3, -5e3, nfault) # 45 degree dipping fault
addelsez!(elsbox, x1, y1, x2, y2, "F")

idx = getidxdict(elsbox);
BRTL_idx = [idx["B"] ; idx["R"] ; idx["T"] ; idx["L"]]
RTL_idx = [idx["R"] ; idx["T"] ; idx["L"]];

npts = 50
offset = 1
xgrid, ygrid = obsgrid(L + offset, B + offset, R - offset, T - offset, npts);

### Slip from fault only
T_TB_F, H_TB_F = PUTC(slip2dispstress, elsbox, BRTL_idx, idx["F"], mu, nu)
Fslip = [100; 100]; # y-direction slip only
Uslip = T_TB_F * Fslip
Tslip = H_TB_F * Fslip;

# Kernels and solve
T_B_BRTL, H_B_BRTL = PUTC(slip2dispstress, elsbox, idx["B"], BRTL_idx, mu, nu)
T_RTL_BRTL, H_RTL_BRTL = PUTC(slip2dispstress, elsbox, RTL_idx, BRTL_idx, mu, nu)
bcs_fault_only = zeros(8 * nside)
bcs_fault_only[1:2 * nside] = -Uslip[1:2 * nside] # Bottom
bcs_fault_only[2 * nside + 1:4 * nside] = -Tslip[2 * nside + 1:4 * nside] # Right
bcs_fault_only[4 * nside + 1:6 * nside] = -Tslip[4 * nside + 1:6 * nside] # Top
bcs_fault_only[6 * nside + 1:8 * nside] = -Tslip[6 * nside + 1:8 * nside] # Left
THbox = [T_B_BRTL ; H_RTL_BRTL];

THboxinv = inv(THbox);

Ueff_fault_only = THbox \ bcs_fault_only
UTB, STB = constdispstress(slip2dispstress, xgrid, ygrid, elsbox, BRTL_idx,
                           Ueff_fault_only[1:2:end], Ueff_fault_only[2:2:end], mu, nu)
UF, SF = constdispstress(slip2dispstress, xgrid, ygrid, elsbox, idx["F"],
                         Fslip[1:2:end], Fslip[2:2:end], mu, nu)
Ufaultonly = UTB .+ UF
Sfaultonly = STB .+ SF
twopanel(xgrid, ygrid, npts, Ufaultonly, Sfaultonly, idx, elsbox, figsize=(12, 6), ylim=[B,T]);

fenics_mesh = RectangleMesh(Point((L - 5000, T + 5000)), Point((R + 5000, B - 5000)), 60, 60)
scalar_fs = FunctionSpace(fenics_mesh, "P", 1)
vector_fs = VectorFunctionSpace(fenics_mesh, "P", 2)
tensor_fs = FunctionSpace(FEniCS.fenics.TensorFunctionSpace(fenics_mesh.pyobject, "P", 1))

include("Helpers.jl")

x_dofs = py"lambda x: x.tabulate_dof_coordinates()"(scalar_fs.pyobject);
_, stress_from_fault = constdispstress(slip2dispstress, x_dofs[:,1], x_dofs[:,2], elsbox, idx["F"],
                                       Fslip[1:2:end], Fslip[2:2:end], mu, nu)

# no viscous behavior in the surface layer, mu/eta below VE_z
# We can't use a threshold for eta because we divide by eta and dividing by zero is bad.
fenics_mu_over_eta = interpolate(Expression("0.0 + mu/eta*(x[1] <= VE_z)", degree=0, mu=mu, eta=eta, VE_z=VE_z), scalar_fs)
fenics_Ve = FeFunction(vector_fs)
fenics_disp = FeFunction(vector_fs)

v = TestFunction(vector_fs)
bc = DirichletBC(vector_fs, Constant((0, 0)), "on_boundary")
function fenics_lhs_assemble(V)
    u = TrialFunction(V)
    a = inner(fenics_stress(u, lambda, mu), fenics_strain(v)) * dx
    A = assemble(a)
    apply(bc, A)
    A
end
function fenics_rhs_assemble(V, f)
    L = FEniCS.dot(f, v) * dx
    b = assemble(L)
    apply(bc, b)
    b
end
lhs = fenics_lhs_assemble(vector_fs)
lu_solver = py"""
from fenics import *
def create_lu_solver(lhs):
    s = LUSolver(lhs, method='umfpack')
    s.parameters['symmetric'] = True
    return s
"""
lu_solver = py"create_lu_solver"(lhs.pyobject)
function fenics_lu_solve(lu_solver, u, rhs)
    py"lambda lu_solver, u, rhs: lu_solver.solve(u.vector(), rhs)"(lu_solver, u.pyobject, rhs.pyobject)
end

# Step 1: solve FEM+BEM problem
function solve_elastic_problem(dofs)
    vector(fenics_Ve).set_local(dofs)
    rhs = fenics_rhs_assemble(vector_fs, fenics_Ve)
    fenics_lu_solve(lu_solver, fenics_disp, rhs)

    fenics_sigma = fenics_stress(fenics_disp, lambda, mu)
    sxx_sxy, sxy_syy = fenics_tensor_components(vector_fs, fenics_sigma)

    bcs_B = -fenics_eval_u(fenics_disp, idx["B"], elsbox.xcenter, elsbox.ycenter)
    bcs_RTL_S = -fenics_eval_s(sxx_sxy, sxy_syy, RTL_idx, elsbox.xcenter, elsbox.ycenter)
    bcs_RTL = stress_to_trac(bcs_RTL_S, RTL_idx, elsbox.xnormal, elsbox.ynormal);

    bcs_ve_fem = zeros(2 * elsbox.endidx - 2)
    bcs_ve_fem[1:2 * nside] = bcs_B
    bcs_ve_fem[2 * nside + 1:end] = bcs_RTL;

    bcscombined = bcs_fault_only .+ bcs_ve_fem
    Ueffcombined = THbox \ bcscombined

    fenics_sigma, Ueffcombined
end

# Step 2: convert BEM stress to a fenics object
# Three sources of stress: u_eff, fault, FEM
function calc_total_stress(fenics_sigma, Ueffcombined)
    _, stress_from_ueff = constdispstress(slip2dispstress, x_dofs[:,1], x_dofs[:,2], elsbox, BRTL_idx,
                                          Ueffcombined[1:2:end], Ueffcombined[2:2:end], mu, nu)
    total_bem_stress = stress_from_ueff .+ stress_from_fault

    # set stress values to zero when they are outside the box domain. TODO: I don't think this is right solution even though it's easy.
    total_bem_stress[x_dofs[:,1] .<= L, :] .= 0
    total_bem_stress[x_dofs[:,1] .>= R, :] .= 0
    total_bem_stress[x_dofs[:,2] .<= B, :] .= 0
    total_bem_stress[x_dofs[:,2] .>= T, :] .= 0

    fenics_stress_dofs = zeros(4, size(x_dofs)[1])
    fenics_stress_dofs[1,:] = total_bem_stress[:,1]
    fenics_stress_dofs[2,:] = total_bem_stress[:,3]
    fenics_stress_dofs[3,:] = total_bem_stress[:,3]
    fenics_stress_dofs[4,:] = total_bem_stress[:,2];
    bem_sigma = FeFunction(tensor_fs)
    vector(bem_sigma).set_local(reshape(fenics_stress_dofs, (:)))
    total_sigma = bem_sigma + fenics_sigma
end

# Step 3: calculate the stress deviator and dVdt. Extract dofs from fenics function and return the array.
function stress_to_dVdt(total_sigma)
    # https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Stress_deviator_tensor
    stress_deviator = total_sigma - (1.0 / 3) * FEniCS.tr(total_sigma) * Identity(2);
    dVdt = project(fenics_mu_over_eta * div(stress_deviator), vector_fs)
    vector(dVdt).get_local()
end

function volume_disp_stress(Ueffcombined)
    UTB, STB = constdispstress(slip2dispstress, xgrid, ygrid, elsbox, BRTL_idx,
                               Ueffcombined[1:2:end], Ueffcombined[2:2:end], mu, nu)

    fenics_sigma = fenics_stress(fenics_disp, lambda, mu)
    sxx_sxy, sxy_syy = fenics_tensor_components(vector_fs, fenics_sigma)

    Ufem = transpose(reshape(fenics_eval_u(fenics_disp, 1:length(xgrid), xgrid, ygrid), (2, :)))
    Sfem = transpose(reshape(fenics_eval_s(sxx_sxy, sxy_syy, 1:length(xgrid), xgrid, ygrid), (3, :)));

    Utotal = UTB .+ UF .+ Ufem
    Stotal = STB .+ SF .+ Sfem
    Utotal, Stotal
end

function calc_dVdt(Ve_force_dofs)
    fenics_sigma, Ueffcombined = solve_elastic_problem(Ve_force_dofs)
    total_sigma = calc_total_stress(fenics_sigma, Ueffcombined)
    dVdt  = stress_to_dVdt(total_sigma)
    dVdt, Ueffcombined
end

function main()
    V = zeros(vector(fenics_Ve).size())
    dt = siay * 2
    t = 0
    Uhistory = []
    for i in 1:100
        println(t)
        dVdt, Ueffcombined = calc_dVdt(V)
        V += dt .* dVdt
        t += dt
        Utotal, Stotal = volume_disp_stress(Ueffcombined)
        append!(Uhistory, [Utotal])
    # @time twopanel(xgrid, ygrid, npts, Utotal, Stotal, idx, elsbox, figsize=(12,6), ylim=[B,T]);
    end
    Uhistory = reshape(hcat(Uhistory...), (:, 2, length(Uhistory)));

    # sanity check: velocity measured in mm/yr
    Vhistory = (Uhistory[:,:,2:end] - Uhistory[:,:,1:end - 1]) ./ dt * siay * 1000;

    maxV = dropdims(maximum(sqrt.(Vhistory[:,1,:].^2 + Vhistory[:,2,:].^2), dims=(1)), dims=(1));

    maxV

    plot(maxV)

    for i in 1:10:41
        figure()
        field = sqrt.(Vhistory[:, 1, i].^2 + Vhistory[:, 2, i].^2)
        ncontours = 10
        lowfield = minimum(field)
        highfield = maximum(field)
        ncontours = LinRange(lowfield, highfield, 11)

        xlim = [-20000 20000]
        scale = 1.0
        fieldmax = maximum(@.abs(field))
        contourf(reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
             reshape(field, npts, npts), ncontours,
             vmin=lowfield, vmax=highfield,
             cmap=get_cmap("magma"))
        clim(lowfield, highfield)
        colorbar(fraction=0.020, pad=0.05, extend="both", label=L"$||u||$ (m)")
        contour(reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
            reshape(field, npts, npts), ncontours,
            vmin=lowfield, vmax=highfield,
            linewidths=0.25, colors="w")
        gca().set_aspect("equal")
        gca().set_xlim([-20000, 20000])
        gca().set_ylim([-50000, 0])
        xlabel(L"$x$ (m)")
        ylabel(L"$y$ (m)")
    end

    plot(log10.(maxV))

end
main()


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_instances_by_tag("ClusterName", ec2_clustername)
hosts = ec2_hostnames(instances)

println("Creating addprocs_headnode.jl")
f = open("addprocs_headnode.jl", "w")
machines = ec2_addprocs(instances, ec2_sshkey_file;
                        use_public_dnsname=false,
                        num_workers=ec2_num_workers,
                        workers_per_instance=ec2_workers_per_instance,
                        machines_only=true)

for mset in machines
    print(f, "addprocs([")
    lenmset = length(mset)
    for (i,m) in enumerate(mset)
        print(f, "\"$m", i == lenmset ? "\"" : "\"," )
    end
    print(f, "]; sshflags=`-o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no -o LogLevel=ERROR`)\n\n")
end
println(f, "println(\"nworkers = \$(nworkers())\")")
flush(f)
close(f)

println("Copy files to head node into directory /home/ubuntu/run")

sshflags = `-o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no -o LogLevel=ERROR -i $(ec2_sshkey_file)`
run(`ssh $sshflags ubuntu@$(hosts[1][2]) "mkdir -p /home/ubuntu/run"`)
run(`scp $sshflags compute_headnode.jl config.jl addprocs_headnode.jl test_workers.jl ubuntu@$(hosts[1][2]):/home/ubuntu/run`)
if isfile("compute.jl")
    run(`scp $sshflags compute.jl ubuntu@$(hosts[1][2]):/home/ubuntu/run`)
end

println()
println("Use the following commands to connect, cd and execute on the head node")
println()
println("  ssh -i $(ec2_sshkey_file) ubuntu@$(hosts[1][2])")
println("  cd run")
println("  julia compute_headnode.jl")
println("OR")
println("  ssh -i $(ec2_sshkey_file) ubuntu@$(hosts[1][2])")
println("  cd run")
println("  julia")
println("  julia> include(\"compute_headnode.jl\")")
println()


text:
module TestSpectralRadius
using Test
using EcologicalNetworks

# All examples are from Figure 1 of SK&A 2013

A = [1 1 1 1; 1 1 1 1; 1 1 1 0; 1 1 0 0; 1 1 0 0; 1 1 0 0]
N = BipartiteNetwork(A.>0)
@test ρ(N; range=EcologicalNetworks.ρ_raw) ≈ 3.82 atol = 0.01

A = [0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0; 1 1 1 0; 1 0 0 1]
N = BipartiteNetwork(A.>0)
@test ρ(N; range=EcologicalNetworks.ρ_raw) ≈ 3.515 atol = 0.01

A = [1 1 1 1; 1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0; 1 0 0 0]
N = BipartiteNetwork(A.>0)
@test ρ(N; range=EcologicalNetworks.ρ_raw) ≈ 3.944 atol = 0.01

A = [1 0 1 1; 1 1 0 1; 1 1 1 0; 1 1 1 0; 1 1 1 0; 0 1 0 1]
N = BipartiteNetwork(A.>0)
@test ρ(N; range=EcologicalNetworks.ρ_raw) ≈ 3.595 atol = 0.01

end


text:
export NeuralNetworkApproximator, ActorCritic

using Flux
import Functors: functor
using MacroTools: @forward

"""
    NeuralNetworkApproximator(;kwargs)

Use a DNN model for value estimation.

# Keyword arguments

- `model`, a Flux based DNN model.
- `optimizer=nothing`
"""
Base.@kwdef struct NeuralNetworkApproximator{M,O} <: AbstractApproximator
    model::M
    optimizer::O = nothing
end

# some model may accept multiple inputs
(app::NeuralNetworkApproximator)(args...; kwargs...) = app.model(args...; kwargs...)

@forward NeuralNetworkApproximator.model Flux.testmode!,
Flux.trainmode!,
Flux.params,
device

functor(x::NeuralNetworkApproximator) =
    (model = x.model,), y -> NeuralNetworkApproximator(y.model, x.optimizer)

RLBase.update!(app::NeuralNetworkApproximator, gs) =
    Flux.Optimise.update!(app.optimizer, params(app), gs)

Base.copyto!(dest::NeuralNetworkApproximator, src::NeuralNetworkApproximator) =
    Flux.loadparams!(dest.model, params(src))

#####
# ActorCritic
#####

"""
    ActorCritic(;actor, critic, optimizer=ADAM())

The `actor` part must return logits (*Do not use softmax in the last layer!*), and the `critic` part must return a state value.
"""
Base.@kwdef struct ActorCritic{A,C,O} <: AbstractApproximator
    actor::A
    critic::C
    optimizer::O = ADAM()
end

functor(x::ActorCritic) =
    (actor = x.actor, critic = x.critic), y -> ActorCritic(y.actor, y.critic, x.optimizer)

RLBase.update!(app::ActorCritic, gs) = Flux.Optimise.update!(app.optimizer, params(app), gs)


text:
using CSTParser: @cst_str, headof, parentof, check_span, EXPR
jl_parse(s) = CSTParser.remlineinfo!(Meta.parse(s))

function check_parents(x::EXPR)
    if x.args isa Vector{EXPR}
        for a in x.args
            @test a.parent == x
            check_parents(a)
        end
    end
    if x.trivia isa Vector{EXPR}
        for a in x.trivia
            @test a.parent == x
        end
    end
end

function test_iter(ex)
    total = 0
    for x in ex
        @test x isa EXPR
        test_iter(x)
        total += x.fullspan
    end
    if length(ex) > 0
        @test total == ex.fullspan
    end
end

function test_expr(s, head, n, endswithtrivia = false)
    x = CSTParser.parse(s)
    head === nothing || @test headof(x) === head
    @test length(x) === n
    @test x.args === nothing || all(x === parentof(a) for a in x.args)
    @test x.trivia === nothing || all(x === parentof(a) for a in x.trivia)
    @test Expr(x) == jl_parse(s)
    @test isempty(check_span(x))
    check_parents(x)
    test_iter(x)
    @test endswithtrivia ? (x.fullspan-x.span) == (last(x.trivia).fullspan - last(x.trivia).span) : (x.fullspan-x.span) == (last(x.args).fullspan - last(x.args).span)
end

@testset ":local" begin
    test_expr("local a", :local, 2)
end

@testset ":global" begin
    test_expr("global a", :global, 2)
    test_expr("global a, b", :global, 4)
    test_expr("global a, b = 2", :global, 2)
    test_expr("global const a = 1", :const, 3)
    test_expr("global const a = 1, b", :const, 3)
end

@testset ":const" begin
    test_expr("const a = 1", :const, 2)
    test_expr("const global a = 1", :const, 2)
end

@testset ":return" begin
    test_expr("return a", :return, 2)
end

@testset ":abstract" begin
    test_expr("abstract type sig end", :abstract, 4, true)
end

@testset ":primitive" begin
    test_expr("primitive type sig spec end", :primitive, 5, true)
end

@testset ":call" begin
    test_expr("f()", :call, 3, true)
    test_expr("f(a, b)", :call, 6, true)
    test_expr("a + b", :call, 3, false)
end

@testset ":brackets" begin
    test_expr("(a)", :brackets, 3, true)
end

@testset ":begin" begin
    test_expr("begin end", :block, 2, true)
    test_expr("begin a end", :block, 3, true)
    test_expr("quote end", :quote, 3, true)
    test_expr("while cond end", :while, 4, true)
    test_expr("for i = I end", :for, 4, true)
    test_expr("module m end", :module, 5, true)
    test_expr("function f() end", :function, 4, true)
    test_expr("macro f() end", :macro, 4, true)
    test_expr("struct T end", :struct, 5, true)
    test_expr("mutable struct T end", :struct, 6, true)
end

@testset ":export" begin
    test_expr("export a", :export, 2, false)
    test_expr("export a, b", :export, 4, false)
end

@testset ":import" begin
    test_expr("import a", :import, 2, false)
    test_expr("import a, b", :import, 4, false)
    test_expr("import a.b", :import, 2, false)
    test_expr("import a.b, c", :import, 4, false)
    test_expr("import a:b", :import, 2, false)
    test_expr("import a:b.c", :import, 2, false)
    test_expr("import a:b.c, d", :import, 2, false)
end

@testset ":kw" begin
    test_expr("f(a=1)", :call, 4, false)
end

@testset ":tuple" begin
    test_expr("a,b ", :tuple, 3, false)
    test_expr("(a,b) ", :tuple, 5, true)
    test_expr("a,b,c ", :tuple, 5, false)
    test_expr("(a,b),(c) ", :tuple, 3, false)
end

@testset ":curly" begin
    test_expr("x{a}", :curly, 4, true)
    test_expr("x{a,b}", :curly, 6, true)
end

@testset "operators" begin
    test_expr("!a", :call, 2, false)
    test_expr("&a", nothing, 2, false)
    test_expr(":a", :quotenode, 2, false)
    test_expr("a + 1", :call, 3, false)
    test_expr(":(a + 1)", :quote, 2, false)
    test_expr("a ? b : c", :if, 5, false)
    test_expr("a:b", :call, 3, false)
    test_expr("a:b:c", :call, 5, false)
    test_expr("a = b", nothing, 3, false)
    test_expr("a += b", nothing, 3, false)
    test_expr("a < b", :call, 3, false)
    test_expr("a < b < c", :comparison, 5, false)
    test_expr("a^b", :call, 3, false)
    test_expr("!a^b", nothing, 2, false)
    test_expr("a.b", nothing, 3, false)
    test_expr("a.:b", nothing, 3, false)
    test_expr("a + b + c", :call, 5, false)
    test_expr("a where b", :where, 3, false)
    test_expr("a where {b }", :where, 5, true)
    test_expr("a where {b,c }  ", :where, 7, true)
    test_expr("a...", nothing, 2, false)
    @test let x = cst"a... "; x.fullspan - x.span == 1 end
    test_expr("a <: b", nothing, 3, false)

    # https://github.com/julia-vscode/CSTParser.jl/issues/278
    test_expr("*(a)*b*c", :call, 5, false)
    test_expr("+(a)+b+c", :call, 5, false)
    test_expr("(\na +\nb +\nc +\n d\n)", :brackets, 3, true)
end

@testset ":parameters" begin
    test_expr("f(a;b = 1)", nothing, 5, true)
end

@testset "lists" begin
    @testset ":vect" begin
        test_expr("[]", :vect, 2, true)
        test_expr("[a]", :vect, 3, true)
        test_expr("[a, b]", :vect, 5, true)
        test_expr("[a ]", :vect, 3, true)
    end

    @testset ":vcat" begin
        test_expr("[a\nb]", :vcat, 4, true)
        test_expr("[a;b]", :vcat, 4, true)
        test_expr("[a b\nc d]", :vcat, 4, true)
        test_expr("[a\nc d]", :vcat, 4, true)
        test_expr("[a;c d]", :vcat, 4, true)
    end

    @testset ":hcat" begin
        test_expr("[a b]", :hcat, 4, true)
    end
    @testset ":ref" begin
        test_expr("T[a]", :ref, 4, true)
        test_expr("T[a,b]", :ref, 6, true)
    end
    @testset ":typed_hcat" begin
        test_expr("T[a b]", :typed_hcat, 5, true)
    end
    @testset ":typed_vcat" begin
        test_expr("T[a;b]", :typed_vcat, 5, true)
    end
end

@testset ":let" begin
    test_expr("let\n end", :let, 4, true)
    test_expr("let x = 1 end", :let, 4, true)
    test_expr("let x = 1, y =1  end", :let, 4, true)
end

@testset ":try" begin
    test_expr("try catch end", :try, 6, true)
    test_expr("try a catch end", :try, 6, true)
    test_expr("try catch e end", :try, 6, true)
    test_expr("try a catch e end", :try, 6, true)
    test_expr("try a catch e b end", :try, 6, true)
    test_expr("try a catch e b end", :try, 6, true)
    test_expr("try finally end", :try, 8, true)
    test_expr("try finally a end", :try, 8, true)
    test_expr("try a catch e b finally c end", :try, 8, true)
end

@testset ":macrocall" begin
    test_expr("@m", :macrocall, 2, false)
    test_expr("@m a", :macrocall, 3, false)
    test_expr("@m a b", :macrocall, 4, false)
end

@testset ":if" begin
    test_expr("if c end", :if, 4, true)
    test_expr("if c a end", :if, 4, true)
    test_expr("if c a else end", :if, 6, true)
    test_expr("if c a else b end", :if, 6, true)
    test_expr("if c elseif c end", :if, 5, true)
    test_expr("if c a elseif c b end", :if, 5, true)
    test_expr("if c a elseif c b else d end", :if, 5, true)
end

@testset ":do" begin
    test_expr("f() do x end", :do, 4, true)
    test_expr("f() do x,y end", :do, 4, true)
end

@testset "strings" begin
    test_expr("a\"txt\"", :macrocall, 3, false)
    test_expr("a\"txt\"b", :macrocall, 4, false)
end

@testset ":string" begin
    test_expr("\"\$a\"", :string, 4, false)
    test_expr("\" \$a\"", :string, 4, false)
    test_expr("\" \$a \"", :string, 4, false)
end

@testset ":for" begin
    test_expr("for i = I end", :for, 4, true)
    test_expr("for i in I end", :for, 4, true)
    test_expr("for i ∈ I end", :for, 4, true)
    test_expr("for i ∈ I, j in J end", :for, 4, true)
end

@testset "docs" begin
    test_expr("\"doc\"\nT", :macrocall, 4, false)
end

@testset ":generator" begin
    test_expr("(arg for x in X)", :brackets, 3, true)
    test_expr("(arg for x in X if x)", :brackets, 3, true)
    test_expr("(arg for x in X, y in Y)", :brackets, 3, true)
    test_expr("(arg for x in X, y in Y if x)", :brackets, 3, true)
    test_expr("(arg for x in X for  y in Y)", :brackets, 3, true)
    test_expr("(arg for x in X for  y in Y if x)", :brackets, 3, true)
    test_expr("(arg for x in X for y in Y for z in Z)", :brackets, 3, true)
end

@testset ":cmd" begin
    let s = "``"
        x = CSTParser.parse(s)
        x1 = jl_parse(s)
        @test x1 == Expr(x)
    end
    let s = "`a`"
        x = CSTParser.parse(s)
        x1 = jl_parse(s)
        @test x1 == Expr(x)
    end
    let s = "`a \$a`"
        x = CSTParser.parse(s)
        x1 = jl_parse(s)
        @test x1 == Expr(x)
    end
    test_expr("a``", nothing, 3, false)
    test_expr("a`a`", nothing, 3, false)
end

@testset ":macrocall" begin
    test_expr("@m", :macrocall, 2, false)
    test_expr("@m a", :macrocall, 3, false)
    test_expr("@m(a)", :macrocall, 5, false)
    test_expr("@horner(r) + r", nothing, 3, false)
end

@testset "_str" begin
    test_expr("a\"txt\"", :macrocall, 3, false)
    test_expr("a.b\"txt\"", :macrocall, 3, false)
end


text:
using Documenter, Juno

makedocs(modules=Module[Juno],
         doctest=false, clean=true,
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         sitename="Juno Documentation",
         linkcheck = !("skiplinks" in ARGS),
         pages = Any[
         "Home" => "index.md",
         "Manual" => Any[
             "man/installation.md",
             "man/basic_usage.md",
             "man/debugging.md",
             "man/faq.md",
             "man/juno_frontend.md",
             "man/info_developer.md"
         ],
         "Developer Documentation" => Any[
             "devdocs/contributor_guide.md",
             "devdocs/dev_install.md",
             "devdocs/package_info.md",
             "devdocs/juno_ci.md",
             "devdocs/communication.md"
         ]])

deploydocs(
   repo = "github.com/JunoLab/JunoDocs.jl.git",
)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

abstract type MethodTableView; end

"""
    struct InternalMethodTable <: MethodTableView

A struct representing the state of the internal method table at a
particular world age.
"""
struct InternalMethodTable <: MethodTableView
    world::UInt
end

"""
    struct OverlayMethodTable <: MethodTableView

Overlays the internal method table such that specific queries can be redirected to an
external table, e.g., to override existing method.
"""
struct OverlayMethodTable <: MethodTableView
    world::UInt
    mt::Core.MethodTable
end

struct MethodLookupResult
    # Really Vector{Core.MethodMatch}, but it's easier to represent this as
    # and work with Vector{Any} on the C side.
    matches::Vector{Any}
    valid_worlds::WorldRange
    ambig::Bool
end
length(result::MethodLookupResult) = length(result.matches)
function iterate(result::MethodLookupResult, args...)
    r = iterate(result.matches, args...)
    r === nothing && return nothing
    match, state = r
    return (match::MethodMatch, state)
end
getindex(result::MethodLookupResult, idx::Int) = getindex(result.matches, idx)::MethodMatch

"""
    findall(sig::Type, view::MethodTableView; limit::Int=typemax(Int)) -> MethodLookupResult or missing

Find all methods in the given method table `view` that are applicable to the
given signature `sig`. If no applicable methods are found, an empty result is
returned. If the number of applicable methods exceeded the specified limit,
`missing` is returned.
"""
function findall(@nospecialize(sig::Type), table::InternalMethodTable; limit::Int=typemax(Int))
    return _findall(sig, nothing, table.world, limit)
end

function findall(@nospecialize(sig::Type), table::OverlayMethodTable; limit::Int=typemax(Int))
    result = _findall(sig, table.mt, table.world, limit)
    result === missing && return missing
    if !isempty(result)
        if all(match->match.fully_covers, result)
            # no need to fall back to the internal method table
            return result
        else
            # merge the match results with the internal method table
            fallback_result = _findall(sig, nothing, table.world, limit)
            return MethodLookupResult(
                vcat(result.matches, fallback_result.matches),
                WorldRange(min(result.valid_worlds.min_world, fallback_result.valid_worlds.min_world),
                           max(result.valid_worlds.max_world, fallback_result.valid_worlds.max_world)),
                result.ambig | fallback_result.ambig)
        end
    end
    # fall back to the internal method table
    return _findall(sig, nothing, table.world, limit)
end

function _findall(@nospecialize(sig::Type), mt::Union{Nothing,Core.MethodTable}, world::UInt, limit::Int)
    _min_val = RefValue{UInt}(typemin(UInt))
    _max_val = RefValue{UInt}(typemax(UInt))
    _ambig = RefValue{Int32}(0)
    ms = _methods_by_ftype(sig, mt, limit, world, false, _min_val, _max_val, _ambig)
    if ms === false
        return missing
    end
    return MethodLookupResult(ms::Vector{Any}, WorldRange(_min_val[], _max_val[]), _ambig[] != 0)
end

"""
    findsup(sig::Type, view::MethodTableView) -> Tuple{MethodMatch, WorldRange} or nothing

Find the (unique) method `m` such that `sig <: m.sig`, while being more
specific than any other method with the same property. In other words, find
the method which is the least upper bound (supremum) under the specificity/subtype
relation of the queried `signature`. If `sig` is concrete, this is equivalent to
asking for the method that will be called given arguments whose types match the
given signature. This query is also used to implement `invoke`.

Such a method `m` need not exist. It is possible that no method is an
upper bound of `sig`, or it is possible that among the upper bounds, there
is no least element. In both cases `nothing` is returned.
"""
function findsup(@nospecialize(sig::Type), table::InternalMethodTable)
    return _findsup(sig, nothing, table.world)
end

function findsup(@nospecialize(sig::Type), table::OverlayMethodTable)
    result = _findsup(sig, table.mt, table.world)
    result === nothing || return result
    return _findsup(sig, nothing, table.world) # fall back to the internal method table
end

function _findsup(@nospecialize(sig::Type), mt::Union{Nothing,Core.MethodTable}, world::UInt)
    min_valid = RefValue{UInt}(typemin(UInt))
    max_valid = RefValue{UInt}(typemax(UInt))
    result = ccall(:jl_gf_invoke_lookup_worlds, Any, (Any, Any, UInt, Ptr{Csize_t}, Ptr{Csize_t}),
                   sig, mt, world, min_valid, max_valid)::Union{MethodMatch, Nothing}
    return result === nothing ? result : (result, WorldRange(min_valid[], max_valid[]))
end

isoverlayed(::MethodTableView)     = error("unsatisfied MethodTableView interface")
isoverlayed(::InternalMethodTable) = false
isoverlayed(::OverlayMethodTable)  = true

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
using Patchwork
using Compat
using Base.Test

import Patchwork:
           Patch,
           TextNode,
           Elem,
           Overwrite,
           Insert,
           Delete,
           Reorder,
           DictDiff,
           are_equal

# p
p(x...) = Elem(:xhtml, :p, x)

# Structurally same
sameas(x, y) = false
sameas{T}(x::T, y::T) = x == y
sameas{T <: Patch}(x::T, y::T) =
    all([sameas(getfield(x, n), getfield(x, n)) for n in fieldnames(x)])
sameas(l::AbstractArray, r::AbstractArray) = all(map(sameas, l, r))

function sameas(l::Associative, r::Associative)
    if length(r) != length(r)
        return false
    end
    for (k, v) in l
        if !haskey(r, k) || !sameas(v, r[k])
            println("Mismatched dict key ", k, ": ", v, " <> ", r[k])
            return false
        end
    end
    true
end

@testset "Testing Diffs" begin
    @testset "testing are_equal" begin

        @test are_equal(1, 1) == true
        @test are_equal(1, 2) == false
        @test are_equal(1, 1.0) == true
        @test are_equal(1, 1.5) == false
        @test are_equal(:x, "x") == true
        @test are_equal("x", :x) == true
        @test are_equal(:x, :x) == true
        @test are_equal(:x, :y) == false
        a = [:x, :y]
        b = [:x, :y]
        @test are_equal(a, a) == true
        @test are_equal(a, b) == true
    end

    @testset "things that should return empty patches" begin
        e1 = p("a")
        e2 = p("b")
        a = @compat Dict(:x => 1)
        b = @compat Dict(:x => @compat Dict(:y => 1))
        @test isempty(diff(e1, e1))
        @test isempty(diff(p(e1), p(e1)))
        @test isempty(diff(e1, p("a")))
        @test isempty(diff(p(e1, e2), p(e1, e2)))
        @test isempty(diff(p(p("a"), e2), p(e1, p("b"))))
        @test isempty(diff(e1 & a, p("a") & a))
        @test isempty(diff(e1 & b, p("a") & b))
    end

    @testset "testing Overwrite" begin
        e1 = p("a")
        e2 = p("b")
        a = @compat Dict(1=>[Overwrite(TextNode("b"))])
        @test TextNode("a") == TextNode("a")
        @test sameas(diff(e1, e2), a)
    end
end


text:
push!(LOAD_PATH,"..")
using Documenter
using FastHenryHelper
# wait(Timer(1)) # why do I need this?

const PAGES = [
  "Home" => "index.md",
  "Installation" => "install.md",
  "Introduction" => "introduction.md",
  "Examples" => "examples.md",
  "Public API" => "public.md",
]

makedocs(
  modules = [FastHenryHelper],
  doctest   = "doctest" in ARGS,
  linkcheck = "linkcheck" in ARGS,
  clean   = true,
  format = "pdf" in ARGS ? :latex : :html,
  build = "site",
  sitename = "FastHenryHelper.jl",
  authors = "Chris Stook",
  pages = PAGES,
  html_prettyurls = "deploy" in ARGS,
)

if "deploy" in ARGS
  fake_travis = "C:/Users/Chris/fake_travis_FastHenryHelper.jl"
  if isfile(fake_travis)
    include(fake_travis)
  end
  deploydocs(
    repo = "github.com/cstook/FastHenryHelper.jl.git",
    target = "site",
    branch = "gh-pages",
    latest = "master",
    osname = "linux",
    julia  = "1.0",
    deps = nothing,
    make = nothing,
  )
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct FieldMask <: ProtoType
    paths::Base.Vector{AbstractString}
    FieldMask(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FieldMask

export FieldMask


text:
@testset "PIBasis"  begin

#---

using ACE, Random
using Printf, Test, LinearAlgebra, ACE.Testing, StaticArrays
using ACE: evaluate, evaluate_d, NaiveTotalDegree
using ACEbase.Testing: dirfdtest, fdtest, print_tf
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
      EuclideanVectorState, Product1pBasis

#---

@info("Basic test of PIBasis construction and evaluation")

D = NaiveTotalDegree()
maxdeg = 6
ord = 3
φ = ACE.Invariant()

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, D = D)

pibasis = PIBasis(B1p, ord, maxdeg; property = φ)

# generate a configuration
nX = 10
Xs = rand(EuclideanVectorState, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

AA = evaluate(pibasis, cfg)

println(@test(length(pibasis) == length(AA)))

spec = ACE.get_spec(pibasis)
println(@test all(length(b) > 0 for b in spec))

spec_naive = [
    [ ACE.get_spec(B1p, pibasis.spec.iAA2iA[iAA, t])
      for t = 1:pibasis.spec.orders[iAA] ]   for iAA = 1:length(pibasis)
    ]

println(@test spec == spec_naive)

# get inverse Aspec
inv_spec1 = Dict{Any, Int}()
for (i, b1) in enumerate(ACE.get_spec(B1p))
  inv_spec1[b1] = i
end

# a really naive implementation of PIBasis to check correctness
A = evaluate(B1p, cfg)
AA_naive =  [
      real(prod( A[ inv_spec1[ b1 ] ] for b1 in b )) for b in spec ]
println(@test( AA_naive ≈ AA ))

## Testing derivatives

@info("Derivatives of PIbasis")
tmpd = ACE.alloc_temp_d(pibasis, cfg)
AA1, dAA = ACE.evaluate_ed(pibasis, cfg)
println(@test AA1 ≈ AA)

for ntest = 1:30
  Us = randn(SVector{3, Float64}, length(Xs))
  c = randn(length(pibasis))
  F = t -> sum(c .* ACE.evaluate(pibasis, ACEConfig(Xs + t[1] * Us)))
  dF = t -> [ Us' * sum(c .* ACE.evaluate_ed(pibasis, ACEConfig(Xs + t[1] * Us))[2], dims=1)[:] ]
  print_tf(@test fdtest(F, dF, [0.0], verbose=false))
end
println()
##

end


text:
using PhantomJS
using Base.Test

# execjs test

execjs("""
  "use strict";
  console.log('using PhantomJS version ' +
    phantom.version.major + '.' +
    phantom.version.minor + '.' +
    phantom.version.patch);
  phantom.exit();
  """)

# renderhtml test

src = joinpath(dirname(@__FILE__), "../examples/Example Domain.html")

io = open(src)
ret = renderhtml(io, format="pdf")
close(io)

@test isa(ret, Vector{UInt8} )
@test length(ret) > 10000


text:
function discrete_returns(series::Vector{Float64}) :: Vector{Float64}
    n = length(series)
    r = Vector{Float64}(undef, n-1)

    @inbounds for i in 2:n
        r[i-1] = series[i]/series[i-1] - 1
    end

    return r
end

function discrete_returns(series::Array{Float64, 2}) :: Array{Float64, 2}
    (n, num_series) = size(series)
    r = Array{Float64, 2}(undef, n-1, num_series)

    @inbounds for c in 1:num_series, i in 2:n
        r[i-1, c] = series[i, c]/series[i-1, c] - 1
    end

    return r
end

function log_returns(series::Vector{Float64}) :: Vector{Float64}
    n = length(series)
    r = Vector{Float64}(undef, n-1)

    @inbounds for i in 2:n
        r[i-1] = log(series[i]/series[i-1])
    end

    return r
end

function log_returns(series::Array{Float64, 2}) :: Array{Float64, 2}
    (n, num_series) = size(series)
    r = Array{Float64, 2}(undef, n-1, num_series)

    @inbounds for c in 1:num_series, i in 2:n
        r[i-1, c] = log(series[i, c]/series[i-1, c])
    end

    return r
end


text:
#
println("De) Test of the AutoIonization module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-AutoIonization.sum")
setDefaults("method: continuum, asymptotic Coulomb")  ## setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")       ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

wa = Atomic.Computation("xx",  Nuclear.Model(12.); grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 600),
                        initialConfigs  =[Configuration("1s 2s^2 2p^6")],
                        finalConfigs    =[Configuration("1s^2 2s^2 2p^4"), Configuration("1s^2 2s 2p^5"), Configuration("1s^2 2p^6")], process = JAC.Auger,
                        processSettings = AutoIonization.Settings(true, true, true, Tuple{Int64,Int64}[(1,0)], 0., 1.0e6, 2, "Coulomb") )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
using Plots

#=
Continuation of Three. I thought of a different one

This is unfinished. I got stuck on converting the binary [hence the mess in the while loop].
I believe I'm going to have to manually make a binary conversion.
Maybe I'll return to it.
=#

(function ()
  println("Enter Binary Length Please")
  binLen = parse(UInt64, readline())

  flBinHalf = convert(UInt64, floor(((2^(binLen - 1))/2)))
  clBinHalf = convert(UInt64, ceil(((2^(binLen - 1))/2)))

  dnI = flBinHalf * 2
  upI = 1
  I = flBinHalf

  dnIArr = zeros(UInt64, flBinHalf)
  upIArr = zeros(UInt64, clBinHalf)

  while (dnI != clBinHalf && upI != flBinHalf)

     dnIArr[I] = length( string(dnI, base=2))
     dnI -= 1
      I -= 1

      upIArr[upI] = length( string(upI, base=2))
      upI += 1

  end

end)()


text:
include("copyFolder.jl")
using Base.Filesystem

ElementsSEB="/media/wang/ElementsSEB"
LinuxData_1T="/media/wang/LinuxData_1T"
Ex4= "/media/wang/ext4"
TwoTData = "/media/wang/2TData"
########################################

USB_64G ="/media/wang/64G-USB"
JXDocs = "/media/wang/JXdocs/USB_64G"
BackuRoot="backup"

BRIX8250Root=Base.Filesystem.joinpath(  JXDocs )
USB64GRoot=Base.Filesystem.joinpath(  USB_64G ,"" )

srcRoot= BRIX8250Root
tgrRoot=USB64GRoot
for folder in ["a","tem"]
        proCopyFolder( "$srcRoot/$folder","$tgrRoot/$folder")
end


text:
# Julia wrapper for header: /usr/local/cuda/include/cusolver.h

#helper functions
function cusolverDnCreate(handle)
  statuscheck(ccall( (:cusolverDnCreate, libcusolver), cusolverStatus_t, (Ptr{cusolverDnHandle_t},), handle))
end
function cusolverDnDestroy(handle)
  statuscheck(ccall( (:cusolverDnDestroy, libcusolver), cusolverStatus_t, (cusolverDnHandle_t,), handle))
end
function cusolverDnSetStream(handle, streamId)
  statuscheck(ccall( (:cusolverDnSetStream, libcusolver), cusolverStatus_t, (cusolverDnHandle_t, cudaStream_t), handle, streamId))
end
function cusolverDnGetStream(handle, streamId)
  statuscheck(ccall( (:cusolverDnGetStream, libcusolver), cusolverStatus_t, (cusolverDnHandle_t, Ptr{cudaStream_t}), handle, streamId))
end
function cusolverSpCreate(handle)
  statuscheck(ccall( (:cusolverSpCreate, libcusolver), cusolverStatus_t, (Ptr{cusolverSpHandle_t},), handle))
end
function cusolverSpDestroy(handle)
  statuscheck(ccall( (:cusolverSpDestroy, libcusolver), cusolverStatus_t, (cusolverSpHandle_t,), handle))
end
function cusolverSpSetStream(handle, streamId)
  statuscheck(ccall( (:cusolverSpSetStream, libcusolver), cusolverStatus_t, (cusolverSpHandle_t, cudaStream_t), handle, streamId))
end
function cusolverSpGetStream(handle, streamId)
  statuscheck(ccall( (:cusolverSpGetStream, libcusolver), cusolverStatus_t, (cusolverSpHandle_t, Ptr{cudaStream_t}), handle, streamId))
end
function cusolverSpCreateCsrqrInfo(info)
  statuscheck(ccall( (:cusolverSpCreateCsrqrInfo, libcusolver), cusolverStatus_t, (Ptr{csrqrInfo_t},), info))
end
function cusolverSpDestroyCsrqrInfo(info)
  statuscheck(ccall( (:cusolverDestroyCsrqrInfo, libcusolver), cusolverStatus_t, (csrqrInfo_t,), info))
end
function cusolverRfCreate(handle)
  statuscheck(ccall( (:cusolverRfCreate, libcusolver), cusolverStatus_t, (Ptr{cusolverRfHandle_t},), handle))
end
function cusolverRfDestroy(handle)
  statuscheck(ccall( (:cusolverRfDestroy, libcusolver), cusolverStatus_t, (cusolverRfHandle_t,), handle))
end
function cusolverRfSetStream(handle, streamId)
  statuscheck(ccall( (:cusolverRfSetStream, libcusolver), cusolverStatus_t, (cusolverRfHandle_t, cudaStream_t), handle, streamId))
end
function cusolverRfGetStream(handle, streamId)
  statuscheck(ccall( (:cusolverRfGetStream, libcusolver), cusolverStatus_t, (cusolverRfHandle_t, Ptr{cudaStream_t}), handle, streamId))
end


text:
"""
    selectSite!(df::DataFrame, Sites::Vector{String}, ColSite::String)

Filters the data to keep only the selected sites.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColSite::String`: The name of the column of `df` containing the sites.
- `Site::Vector{String}`: The names of the selected sites.
"""
function selectSite!(df::DataFrame, Sites::Vector{String}, ColSite::String)
    filter!(Symbol(ColSite) => x -> x in Sites, df)
end

"""
    selectSite(df::DataFrame, Sites::String, ColSite::Vector{String})

Creates a new dataframe containing only the data for the selected sites.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColSite::String`: The name of the column of `df` containing the sites.
- `Site::Vector{String}`: The names of the selected sites.
"""
function selectSite(df::DataFrame, Sites::Vector{String}, ColSite::String)
    filter(Symbol(ColSite) => x -> x in Sites, df)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
module Clp

if haskey(ENV, "JULIA_CLP_LIBRARY_PATH") || VERSION < v"1.3"
    deps_file = joinpath(dirname(@__DIR__), "deps", "deps.jl")
    if isfile(deps_file)
        include(deps_file)
    else
        error("Clp not properly installed. Please run import `Pkg; Pkg.build(\"Clp\")`.")
    end
else
    import Clp_jll: libClp
end

using CEnum

include("../lib/libClp.jl")

const _CLP_VERSION =
    VersionNumber("$(Clp_VersionMajor()).$(Clp_VersionMinor()).$(Clp_VersionRelease())")

if !(v"1.17.2" <= _CLP_VERSION <= v"1.17.6")
    error(
        "You have installed version $_CLP_VERSION of Clp, which is not " *
        "supported by Clp.jl. If the version change was breaking, changes " *
        "will need to be made to the Julia code. Please open an issue at " *
        "https://github.com/jump-dev/Clp.jl.",
    )
end

include("MOI_wrapper/MOI_wrapper.jl")

# TODO(odow): remove at Clp.jl v1.0.0.
function ClpSolver(args...; kwargs...)
    error(
        "`ClpSolver` is no longer supported. If you are using JuMP, upgrade " *
        "to the latest version and use `Clp.Optimizer` instead. If you are " *
        "using MathProgBase (e.g., via `lingprog`), you will need to upgrade " *
        "to MathOptInterface (https://github.com/jump-dev/MathOptInterface.jl).",
    )
end
export ClpSolver

# Clp exports all `Clp_xxx` symbols. If you don't want all of these symbols in
# your environment, then use `import Clp` instead of `using Clp`.

for sym in names(@__MODULE__, all = true)
    sym_string = string(sym)
    if startswith(sym_string, "Clp_")
        @eval export $sym
    end
end

end


text:
module MultiplicationTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: SPMultiplicationOutNGP, SPMultiplicationOutNPG, SPMultiplicationIn1GNP, SPMultiplicationAGPN,
                 SPMultiplicationOutNPP, SPMultiplicationIn1PNP, SPMultiplicationAPPN, SPMultiplicationOutNΓP,
                 SPMultiplicationOutNPΓ, SPMultiplicationIn1ΓNP, SPMultiplicationAΓPN

@testset "Multiplication node construction through * syntax" begin
    g = FactorGraph()
    a = constant(1.0)
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = a*x
    @test isa(z, Variable)
    @test isa(g.nodes[:multiplication_1], Multiplication)

    g = FactorGraph()
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = 1.0*x
    @test isa(z, Variable)
    @test isa(g.nodes[:multiplication_1], Multiplication)

    g = FactorGraph()
    @RV x ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV z = x*1.0
    @test isa(z, Variable)
    @test isa(g.nodes[:multiplication_1], Multiplication)

    g = FactorGraph()
    @RV x ~ Gamma(constant(1.0), constant(1.0))
    @RV z = x*2.0
    @test isa(z, Variable)
    @test isa(g.nodes[:multiplication_1], Multiplication)
end

#-------------
# Update rules
#-------------

@testset "SPMultiplicationOutNGP" begin
    @test SPMultiplicationOutNGP <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationOutNGP) == Message{GaussianMeanVariance}
    @test isApplicable(SPMultiplicationOutNGP, [Nothing, Message{Gaussian}, Message{PointMass}])

    @test ruleSPMultiplicationOutNGP(nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=3.0), Message(Univariate, PointMass, m=2.0)) == Message(Univariate, GaussianMeanVariance, m=2.0, v=12.0)
    @test ruleSPMultiplicationOutNGP(nothing, Message(Univariate, GaussianMeanVariance, m=1.0, v=3.0), Message(Multivariate, PointMass, m=[2.0])) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(12.0))
    @test ruleSPMultiplicationOutNGP(nothing, Message(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(3.0)), Message(MatrixVariate, PointMass, m=mat(2.0))) == Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(12.0))
end

@testset "SPMultiplicationOutNPG" begin
    @test SPMultiplicationOutNPG <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationOutNPG) == Message{GaussianMeanVariance}
    @test isApplicable(SPMultiplicationOutNPG, [Nothing, Message{PointMass}, Message{Gaussian}])

    @test ruleSPMultiplicationOutNPG(nothing, Message(Univariate, PointMass, m=2.0), Message(Univariate, GaussianMeanVariance, m=1.0, v=3.0)) == Message(Univariate, GaussianMeanVariance, m=2.0, v=12.0)
end

@testset "SPMultiplicationOutNPP" begin
    @test SPMultiplicationOutNPP <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationOutNPP) == Message{PointMass}
    @test isApplicable(SPMultiplicationOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPMultiplicationOutNPP(nothing, Message(Univariate, PointMass, m=2.0), Message(Univariate, PointMass, m=1.0)) == Message(Univariate, PointMass, m=2.0)
    @test ruleSPMultiplicationOutNPP(nothing, Message(Multivariate, PointMass, m=[2.0]), Message(MatrixVariate, PointMass, m=mat(1.0))) == Message(Multivariate, PointMass, m=[2.0])
end

@testset "SPMultiplicationIn1GNP" begin
    @test SPMultiplicationIn1GNP <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationIn1GNP) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(SPMultiplicationIn1GNP, [Message{Gaussian}, Nothing, Message{PointMass}])

    @test ruleSPMultiplicationIn1GNP(Message(Univariate, GaussianWeightedMeanPrecision, xi=1.0, w=3.0), nothing, Message(Univariate, PointMass, m=2.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=12.0)
    @test ruleSPMultiplicationIn1GNP(Message(Multivariate, GaussianWeightedMeanPrecision, xi=[1.0], w=mat(3.0)), nothing, Message(Multivariate, PointMass, m=[2.0])) == Message(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=12.0 + tiny)
    @test ruleSPMultiplicationIn1GNP(Message(Multivariate, GaussianWeightedMeanPrecision, xi=[1.0], w=mat(3.0)), nothing, Message(MatrixVariate, PointMass, m=mat(2.0))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(12.0 + tiny))
end

@testset "SPMultiplicationIn1PNP" begin
    @test SPMultiplicationIn1PNP <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationIn1PNP) == Message{PointMass}
    @test isApplicable(SPMultiplicationIn1PNP, [Message{PointMass}, Nothing, Message{PointMass}])

    @test ruleSPMultiplicationIn1PNP(Message(Univariate, PointMass, m=1.0), nothing, Message(Univariate, PointMass, m=2.0)) == Message(Univariate, PointMass, m=0.5)
    @test ruleSPMultiplicationIn1PNP(Message(Multivariate, PointMass, m=[1.0]), nothing, Message(MatrixVariate, PointMass, m=mat(2.0))) == Message(Multivariate, PointMass, m=[0.5])
end

@testset "SPMultiplicationAGPN" begin
    @test SPMultiplicationAGPN <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationAGPN) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(SPMultiplicationAGPN, [Message{Gaussian}, Message{PointMass}, Nothing])

    @test ruleSPMultiplicationAGPN(Message(Univariate, GaussianWeightedMeanPrecision, xi=1.0, w=3.0), Message(Univariate, PointMass, m=2.0), nothing) == Message(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=12.0)
end

@testset "SPMultiplicationAPPN" begin
    @test SPMultiplicationAPPN <: SumProductRule{Multiplication}
    @test outboundType(SPMultiplicationAPPN) == Message{PointMass}
    @test isApplicable(SPMultiplicationAPPN, [Message{PointMass}, Message{PointMass}, Nothing])

    @test ruleSPMultiplicationAPPN(Message(Univariate, PointMass, m=1.0), Message(Univariate, PointMass, m=2.0), nothing) == Message(Univariate, PointMass, m=0.5)
end

@testset "SPMultiplicationOutNPΓ" begin
    @test SPMultiplicationOutNPΓ <: SumProductRule{Multiplication}

    @test outboundType(SPMultiplicationOutNPΓ) == Message{Gamma}
    @test isApplicable(SPMultiplicationOutNPΓ, [Nothing, Message{PointMass}, Message{Gamma}])

    @test ruleSPMultiplicationOutNPΓ(nothing, Message(Univariate, PointMass, m=2.0), Message(Univariate, Gamma, a=1.0, b=2.0)) == Message(Gamma, a=1.0, b=1.0)
end

@testset "SPMultiplicationOutNΓP" begin
    @test SPMultiplicationOutNΓP <: SumProductRule{Multiplication}

    @test outboundType(SPMultiplicationOutNΓP) == Message{Gamma}
    @test isApplicable(SPMultiplicationOutNΓP, [Nothing, Message{Gamma}, Message{PointMass}])

    @test ruleSPMultiplicationOutNΓP(nothing, Message(Univariate, Gamma, a=1.0, b=2.0), Message(Univariate, PointMass, m=2.0)) == Message(Gamma, a=1.0, b=1.0)
end

@testset "SPMultiplicationIn1ΓNP" begin
    @test SPMultiplicationIn1ΓNP <: SumProductRule{Multiplication}

    @test outboundType(SPMultiplicationIn1ΓNP) == Message{Gamma}
    @test isApplicable(SPMultiplicationIn1ΓNP, [Message{Gamma}, Nothing, Message{PointMass}])

    @test ruleSPMultiplicationIn1ΓNP( Message(Univariate, Gamma, a=1.0, b=2.0), nothing, Message(Univariate, PointMass, m=2.0)) == Message(Gamma, a=1.0, b=4.0)
end

@testset "SPMultiplicationAΓPN" begin
    @test SPMultiplicationAΓPN <: SumProductRule{Multiplication}

    @test outboundType(SPMultiplicationAΓPN) == Message{Gamma}
    @test isApplicable(SPMultiplicationAΓPN, [Message{Gamma}, Message{PointMass}, Nothing])

    @test ruleSPMultiplicationAΓPN( Message(Univariate, Gamma, a=1.0, b=2.0), Message(Univariate, PointMass, m=2.0), nothing) == Message(Gamma, a=1.0, b=4.0)
end

@testset "messagePassingAlgorithm" begin
    fg = FactorGraph()

    @RV τ ~ Gamma(1, 1)
    @RV x ~ GaussianMeanPrecision(0.0, 10*τ)
    placeholder(x, :x)

    q = PosteriorFactorization(τ, ids=[:T])
    algo = messagePassingAlgorithm(free_energy=true)

    @test isa(algo, InferenceAlgorithm)
end

end # module


text:
using Emerald
using Test

include("general/types.jl")


text:
module SignalDecomposition

export decompose, Decomposition

"Supertype of all decomposition methods."
abstract type Decomposition end

"""
    decompose([t, ] s, method::Decomposition) → x, r
Decompose an 1D input signal or timeseries `s(t)` into components, `x, r`,
using the given `method`. `t` defaults to `1:length(s)`.

What are `x` and `r` really depends on your point of view and your application.
They can be structure `x` and noise `r` (i.e. noise reduction). They can be
seasonal/periodic `x` and residual component `r`. They can even be multiplier `x` and
input `r`.
"""
decompose(s::AbstractVector, method::Decomposition; kwargs...) =
decompose(1:length(s), s, method; kwargs...)

using Statistics
include("utils.jl")
include("linear/fourier.jl")
include("linear/lpv.jl")
include("product/matrixinversion.jl")
include("nonlinear/extremelysimple.jl")
include("nonlinear/projection.jl")
include("misc/anomaly.jl")

end # module


text:
mutable struct TransmissionOperations{R,I,T,P}

    # Parameters

    limits::Vector{Float64}

    # Variables

    flows::Array{VariableRef,3} # (i x t x p)

    # Expressions

    exports::Array{ExpressionRef,3} # (r x t x p)

    # Constraints

    maxflow_forward::Array{LessThanConstraintRef,3} # (i x t x p)
    maxflow_back::Array{LessThanConstraintRef,3}

    function TransmissionOperations{R,T,P}(
        limits::Vector{Float64}
    ) where {R,T,P}

        I = length(limits)
        @assert all(x -> x >= 0, limits)

        new{R,I,T,P}(limits)

    end

end

function TransmissionOperations{T,P}(
    interfaces::Interfaces{I,R}, transmissionpath::String
) where {I,R,T,P}

    transmissiondata =
        DataFrame!(CSV.File(joinpath(transmissionpath, "parameters.csv"),
                            types=scenarios_transmission_param_types))

    interfacelookup = Dict(zip(interfaces.name, 1:I))
    limits = zeros(Float64, I)

    for row in eachrow(transmissiondata)
        int_idx = interfacelookup[row.interface]
        limits[int_idx] = row.limit
    end

    return TransmissionOperations{R,T,P}(limits)

end

function setup!(
    tx::TransmissionOperations{R,I,T,P},
    interfaces::Interfaces{I,R},
    m::Model
) where {R,I,T,P}

    # Variables

    tx.flows =
        @variable(m, [i in 1:I, t in 1:T, p in 1:P])

    # Expressions

    tx.exports =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    sum(flowout(r, l, tx.flows[i,t,p])
                        for (i, l) in enumerate(interfaces.regions)))

    # Constraints

    tx.maxflow_forward =
        @constraint(m, [i in 1:I, t in 1:T, p in 1:P],
                    tx.flows[i,t,p] <= tx.limits[i])

    tx.maxflow_back =
        @constraint(m, [i in 1:I, t in 1:T, p in 1:P],
                    -tx.limits[i] <= tx.flows[i,t,p])

end

flowout(r::Int, label::Tuple{Int,Int}, x) =
    if r == first(label)
        x
    elseif r == last(label)
        -x
    else
        zero(x)
    end


text:
export move2zero, move2xplus, rotation, reflect_across, same_side
export in_up, up_in

# This code enables LFT's to act on HObjects.

(f::LFT)(p::HPoint) = HPoint(f(getz(p)))

(f::LFT)(L::HSegment) = HSegment(f(L.A), f(L.B))

(f::LFT)(T::HTriangle) = HTriangle(f(T.A), f(T.B), f(T.C))

function (f::LFT)(L::HLine)
    s = L.s
    t = L.t
    w = exp(s*im)
    z = exp(t*im)
    ww = f(w)
    zz = f(z)
    ss = angle(ww)
    tt = angle(zz)
    return HLine(ss,tt)
end

function (f::LFT)(R::HRay)
    z = exp(im*R.t)
    p = R.pt
    zz = f(z)
    tt = angle(zz)
    pp = f(p)
    return HRay(pp,tt)
end

function (f::LFT)(HC::Horocycle)
    p = HC.pt
    t = HC.t

    pp = f(p)
    z = exp(t*im)
    zz = f(z)
    tt = angle(zz)

    return Horocycle(pp,tt)
end

function (f::LFT)(X::HPolygon)
    pts = f.(X.plist)
    return HPolygon(pts)
end

(f::LFT)(H::HPlane) = HPlane()

# Various LFT creation functions.

"""
`in_up` is a `LFT` that maps the Poincaré disk to the upper half plane.
"""
const in_up = LFT(-im, -im, 1, -1)

"""
`up_in` is a `LFT` that maps the upper half plane to the Poincaré disk.
"""
const up_in = LFT(-1, im, -1, -im)

"""
`move2zero(P::Hpoint)`
returns a `LFT` that's an isometry of H^2 that maps `P` to the origin.
"""
function move2zero(z::Complex)::LFT
    # map to upper half plane and find x-displacement
    zz = in_up(z)
    x = real(zz)

    # move horizontally to place above origin
    f = LFT(1, -x, 0, 1)
    zz = f(zz)

    # move down to 0 + im
    y = imag(zz)
    g =  LFT( 1, -y*im, 1, y*im )

    return g*f*in_up
end

move2zero(P::HPoint) = move2zero(getz(P))

"""
`rotation(theta)` is an isometry of H^2 corresponding to a
rotation about the origin of the amount `theta`
"""
rotation(theta::Real)= LFT( exp(im*theta), 0, 0, 1)

"""
`move2xplus(P::HPoint)` returns an isometry of H^2 that maps `P` onto
the positive real axis.
"""
function move2xplus(z::Complex)::LFT
    if z == 0
        return LFT()
    end
    theta = angle(z)
    return rotation(-theta)
end

move2xplus(P::HPoint) = move2xplus(getz(P))

function move2xplus(a::Complex, b::Complex)
    f = move2zero(a)
    bb = f(b)
    theta = angle(bb)
    g = rotation(-theta)
    return g*f
end

"""
`move2xplus(A,B)` or `move2xplus(L::HSegment)`
gives an isometry `f` so that `f(A)` is 0 and `f(B)` is on the
positive real axis.
"""
move2xplus(A::HPoint, B::HPoint) = move2xplus(getz(A),getz(B))
move2xplus(L::HSegment) = move2xplus(endpoints(L)...)

"""
`move2xplus(L::HLine)` returns a linear fractional transformation
that maps points on `L` to the positive x-axis but is *not* an
isometry of the hyperbolic plane.
"""
function move2xplus(L::HLine)
    a = exp(im*L.s)
    b = getz(point_on_line(L))
    c = exp(im*L.t)
    f = LFT(a,b,c)
    return f
end

function move2xplus(R::HRay)
    A = get_vertex(R)
    w = getz(A)
    z = exp(im * R.t)
    return move2xplus(w,z)
end

"""
`reflect_across(X::HObject,L::HSegment/HLine)` returns the object
formed by refecting `X across the line segment/line `L`.
"""
function reflect_across(p::HPoint, L::Union{HLine,HSegment})
    f = move2xplus(L)
    z = getz(p)
    zz = f(z)'
    w = (inv(f))(zz)
    return HPoint(w)
end

function reflect_across(S::HSegment, L::Union{HLine,HSegment})
    A,B = endpoints(S)
    AA = reflect_across(A,L)
    BB = reflect_across(B,L)
    return HSegment(AA,BB)
end

function reflect_across(T::HTriangle, L::Union{HLine,HSegment})
    A,B,C = endpoints(T)
    AA = reflect_across(A,L)
    BB = reflect_across(B,L)
    CC = reflect_across(C,L)
    return HTriangle(AA,BB,CC)
end

function reflect_across(X::HPolygon, L::Union{HLine,HSegment})
    pts = [ reflect_across(p,L) for p in X.plist ]
    return HPolygon(pts)
end

function reflect_across(X::HLine, L::Union{HLine,HSegment})
    a = exp(im * L.s)
    b = getz(point_on_line(L))
    c = exp(im * L.t)
    f = LFT(a,-1,b,0,c,1)

    Y = f(X)
    s = Y.s
    t = Y.t
    Z = HLine(-s,-t)
    return (inv(f))(Z)
end

reflect_across(X::HPlane, L::Union{HLine,HSegment}) = HPlane()

"""
`same_side(P,Q,L)` determines if the points `P` and `Q`
lie in the same (closed) halfplane as determined by `L`.
If either point is on `L` then the result is `true`.
"""
function same_side(P::HPoint, Q::HPoint, L::HLine)::Bool

    if in(P,L) || in(Q,L)
        return true
    end

    f = move2xplus(L)
    a = imag(f(getz(P)))
    b = imag(f(getz(Q)))

    return sign(a) == sign(b)
end

same_side(P::HPoint, Q::HPoint, S::Union{HSegment,HRay}) = same_side(P,Q,HLine(S))

## UNARY MINUS

"""
`-X` where `X` is a hyperbolic object is a new object reflected through
the origin.
"""
(-)(H::HPlane) = HPlane()
(-)(P::HPoint) = HPoint(-getz(P))
(-)(P::HPolygon) = HPolygon((-).(P.plist))
(-)(C::HCircle) = HCircle( -(C.ctr), C.rad )

function (-)(L::HSegment)
    a,b = endpoints(L)
    return HSegment(-a,-b)
end

function (-)(R::HRay)
    p = -R.pt
    t = R.t + pi
    return HRay(p,t)
end

function (-)(T::HTriangle)
    a,b,c = endpoints(T)
    return HTriangle(-a,-b,-c)
end

function (-)(L::HLine)
    s = L.s
    t = L.t
    return HLine(pi+s,pi+t)
end

function (-)(HC::Horocycle)
    pt = -HC.pt
    t  = HC.t + pi
    return Horocycle(pt,t)
end

## ADJOINT
"""
`adjoint(X::HObject)` (that is, `X'`) returns a new `X` that is reflected
across the `x`-axis.
"""
adjoint(H::HPlane) = HPlane()
adjoint(P::HPoint) = HPoint(getz(P)')
adjoint(P::HPolygon) = HPolygon(adjoint.(P.plist))
adjoint(C::HCircle) = HCircle( (C.ctr)', C.rad )

function adjoint(L::HSegment)
    a,b = endpoints(L)
    return HSegment(a',b')
end

function adjoint(L::HLine)
    s = L.s
    t = L.t
    return HLine(-s,-t)
end

function adjoint(R::HRay)
    p = (R.pt)'
    t = -(R.t)
    return HRay(p,t)
end

function adjoint(HC::Horocycle)
    pt = (HC.pt)'
    t  = -HC.t
    return Horocycle(pt,t)
end

function adjoint(T::HTriangle)
    a,b,c = endpoints(T)
    return HTriangle(a',b',c')
end


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rastername, rasterpath, rasterurl

@testset "EarthEnv LandCover" begin
    using RasterDataSources: rastername, rasterurl, rasterpath

    @test rastername(EarthEnv{LandCover{:DISCover}}, 2) == "consensus_full_class_2.tif"
    landcover_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "EarthEnv", "LandCover")
    @test rasterpath(EarthEnv{LandCover{:DISCover}}) == joinpath(landcover_path, "with_DISCover")
    @test rasterpath(EarthEnv{LandCover}) == joinpath(landcover_path, "without_DISCover")
    @test rasterpath(EarthEnv{LandCover{:DISCover}}, 2) == joinpath(landcover_path, "with_DISCover", "consensus_full_class_2.tif")
    @test rasterurl(EarthEnv{LandCover{:DISCover}}, 2) ==
        URI(scheme="https", host="data.earthenv.org", path="/consensus_landcover/with_DISCover/consensus_full_class_2.tif")
    without_discover_path = joinpath(landcover_path, "without_DISCover", "Consensus_reduced_class_2.tif")
    with_discover_path = joinpath(landcover_path, "with_DISCover", "consensus_full_class_2.tif")
    @test getraster(EarthEnv{LandCover{:DISCover}}, 2) == with_discover_path
    @test isfile(with_discover_path)
    @test getraster(EarthEnv{LandCover}, :evergreen_broadleaf_trees) == without_discover_path
    @test getraster(EarthEnv{LandCover}, (2,)) == (evergreen_broadleaf_trees=without_discover_path,)
    @test getraster(EarthEnv{LandCover}, (:evergreen_broadleaf_trees,)) == (evergreen_broadleaf_trees=without_discover_path,)
    @test getraster(EarthEnv{LandCover}, 2:2) == (evergreen_broadleaf_trees=without_discover_path,)
    @test getraster(EarthEnv{LandCover}, [:evergreen_broadleaf_trees]) == (evergreen_broadleaf_trees=without_discover_path,)
    @test isfile(without_discover_path)
    for layer in RasterDataSources.layers(EarthEnv{LandCover})
        getraster(EarthEnv{LandCover{:DISCover}}, layer)
    end
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
function constructors(n, es)
  if n < 32
    floatdef = "f_to_p$(n)e$(es)(a)"
    doubledef = "f_to_p$(n)e$(es)((float) a)"
  else
    floatdef = "f_to_p$(n)e$(es)((double) a)"
    doubledef = "f_to_p$(n)e$(es)(a)"
  end

  """

  $(c(n,es))::$(c(n,es))(){ this->data = P$(n)ZER; }

  $(c(n,es))::$(c(n,es))(const float a){
    this->data = $floatdef.udata;
    if (errno) {
      throw std::domain_error("attempted to construct a posit from a NaN IEEE value");
    }
  }

  $(c(n,es))::$(c(n,es))(const double a){
    this->data = $doubledef.udata;
    if (errno) {
      throw std::domain_error("attempted to construct a posit from a NaN IEEE value");
    }
  }

  $(c(n,es))::$(c(n,es))(const $(c(n,es)) &a){ this->data = a.data; }

  $(c(n,es))::$(c(n,es))(const $(p(n,es)) a){ this->data = a.udata; }

  """
end

const CPP_ARITH_OPS = Dict(:add => :+, :sub => :-, :mul => :*, :div => :/)

function arith_operators(n,es)
  #start with the negation code.
  operator_code = ["""
  $(c(n,es)) $(c(n,es))::operator -() const{
    $(c(n,es)) res;
    res.data = -(this->data);
    return res;
  }
  """]

  for op in keys(CPP_ARITH_OPS)
    arith_opsymbol = CPP_ARITH_OPS[op]
    assn_opsymbol = Symbol(arith_opsymbol, :(=))
    #creates the reassignment binding then the arithmetic binding.
    push!(operator_code, """
    $(c(n,es)) &$(c(n,es))::operator $(assn_opsymbol)(const $(c(n,es)) rhs){
      $(p(n,es)) res;

      res = $(p(n,es,op))_e(*this, rhs);
      if (errno) {
        throw std::domain_error("NaN value obtained in operator $assn_opsymbol");
      }

      this->data = res.udata;
      return (*this);
    }

    $(c(n,es)) $(c(n,es))::operator $(arith_opsymbol)(const $(c(n,es)) rhs) const{
      $(c(n,es)) res;          //create a return value from the stack.

      res = $(p(n,es))($(p(n,es,op))_e(*this, rhs));

      if (errno) {
        throw std::domain_error("NaN value obtained in operator $arith_opsymbol");
      }

      return $(c(n,es))(res);
    }
    """)
  end

  join(operator_code, "\n")
end

const CPP_BOOL_OPS = Dict(:lt => "<", :gt => ">", :lte => "<=", :gte => ">=", :eq => "==")

function bool_operators(n,es)
  operator_code = []

  for op in keys(CPP_BOOL_OPS)
    bool_opsymbol = CPP_BOOL_OPS[op]
    #creates the reassignment binding then the arithmetic binding.
    push!(operator_code, """
    bool $(c(n,es))::operator $(bool_opsymbol)(const $(c(n,es)) rhs) const{
      $(p(n,es)) lhs_p, rhs_p;
      lhs_p.udata = this->data;
      rhs_p.udata = rhs.data;

      return $(p(n,es,op))(lhs_p, rhs_p);
    }

    """)
  end

  join(operator_code, "\n")
end

function conversion_operators(n,es)
  """
  $(c(n,es))::operator float() const{
    return (float) p$(n)e$(es)_to_f(($(p(n,es)))(*this));
  }

  $(c(n,es))::operator double() const{
    return (double) p$(n)e$(es)_to_f(($(p(n,es)))(*this));
  }

  $(c(n,es))::operator $(p(n,es))() const{
    $(p(n,es)) res;
    res.udata = this->data;
    return res;
  }
  """
end

################################################################################

const CPP_UNARY_FNS = [:mulinv, :log2, :exp2 , :sqrt, :log1p, :log, :log10, :exp, :sin, :cos, :atan]

function unary_functions(n, es)
  function_code = []

  for fn in CPP_UNARY_FNS
    #creates the reassignment binding then the arithmetic binding.
    push!(function_code, """

    $(c(n,es)) $(fn) (const $(c(n,es)) x){
      $(p(n,es)) res;

      res = $(p(n,es,fn))_e(x);

      if (errno) {
        throw std::domain_error("NaN value obtained in function $fn");
      }

      $(c(n,es)) res_c(res);
      return res_c;
    }
    """)
  end

  join(function_code, "\n")
end

################################################################################

const CPP_BINARY_FNS = [:pow, :atan2]

function binary_functions(n, es)
  function_code = []

  for fn in CPP_BINARY_FNS
    #creates the reassignment binding then the arithmetic binding.
    push!(function_code, """

    $(c(n,es)) $(fn) (const $(c(n,es)) a, const $(c(n,es)) b){
      $(p(n,es)) res;

      res = $(p(n,es,fn))_e(a, b);
      if (errno) {
        throw std::domain_error("NaN value obtained in function $fn");
      }

      $(c(n,es)) res_c(res);
      return res_c;
    }
    """)
  end

  join(function_code, "\n")
end

################################################################################

const CPP_TERNARY_FNS = [:fma, :fms, :nfma, :nfms]

function ternary_functions(n, es)
  function_code = []

  for fn in CPP_TERNARY_FNS
    #creates the reassignment binding then the arithmetic binding.
    push!(function_code, """

    $(c(n,es)) $(fn)(const $(c(n,es)) a, const $(c(n,es)) b, const $(c(n,es)) c){
      $(p(n,es)) res;

      res = $(p(n,es,fn))_e(a, b, c);

      if (errno) {
        throw std::domain_error("NaN value obtained in function $fn");
      }

      $(c(n,es)) res_c(res);
      return res_c;
    }
    """)
  end

  join(function_code, "\n")
end

doc"""
  adlibbed- file which fills in the implementation of the functions
"""
function generate_posit_class_cpp(io, n, es)
  write(io, """

  #include "include/posit_conv.h"
  #include "include/posit_ops.h"
  #include "include/posit_ops_errno.h"
  #include "include/$(c(n,es)).h"
  #include "errno.h"
  #include <stdexcept>

  $(constructors(n,es))

  $(arith_operators(n,es))

  $(bool_operators(n,es))

  $(conversion_operators(n,es))

  bool isequal(const $(c(n,es)) a, const $(c(n,es)) b){
    return a.data == b.data;
  }

  $(unary_functions(n,es))

  $(binary_functions(n,es))

  $(ternary_functions(n,es))
  """)

end


text:
using SimilarityMetrics

# Set

a, b = Set{Int}(), Set{Int}()

add!(a, 1); add!(a, 2)
add!(b, 2); add!(b, 3)
@assert matching(a, b) == length(intersect(a, b))

add!(a, 4)
@assert matching(a, b) == length(intersect(a, b))

add!(b, 1)
@assert matching(a, b) == length(intersect(a, b))

# IntSet

a, b = IntSet(), IntSet()

add!(a, 1); add!(a, 2)
add!(b, 2); add!(b, 3)
@assert matching(a, b) == length(intersect(a, b))

add!(a, 4)
@assert matching(a, b) == length(intersect(a, b))

add!(b, 1)
@assert matching(a, b) == length(intersect(a, b))

# AbstractVector

a, b = falses(4), falses(4)

a[1] = true; a[2] = true
b[2] = true; b[3] = true
@assert matching(a, b) == 1

a[4] = true
@assert matching(a, b) == 1

b[1] = true
@assert matching(a, b) == 2


text:
@doc raw"""
    LocalLevelCycle(y::Vector{Fl}) where Fl

The local level model with a cycle component is defined by:
```math
\begin{gather*}
    \begin{aligned}
        y_{t} &=  \mu_{t} + c_{t} + \varepsilon_{t} \quad &\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\
        \mu_{t+1} &= \mu_{t} + \eta_{t} \quad &\eta_{t} \sim \mathcal{N}(0, \sigma^2_{\eta})\\
        c_{t+1} &= c_{t} \cos(\lambda_c) + c_{t}^{*} \sin(\lambda_c)\ \quad & \tilde\omega_{t} \sim \mathcal{N}(0, \sigma^2_{\tilde\omega})\\
        c_{t+1}^{*} &= -c_{t} \sin(\lambda_c) + c_{t}^{*} \sin(\lambda_c) \quad &\tilde\omega^*_{t} \sim \mathcal{N}(0, \sigma^2_{\tilde\omega})\\
    \end{aligned}
\end{gather*}
```

# Example
```jldoctest
julia> model = LocalLevelCycle(rand(100))
LocalLevelCycle model
```

See more on TODO RJ_TEMPERATURE

# References
 * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods: Second Edition." Oxford University Press. pp. 48
"""
mutable struct LocalLevelCycle <: StateSpaceModel
    hyperparameters::HyperParameters
    system::LinearUnivariateTimeInvariant
    results::Results

    function LocalLevelCycle(y::Vector{Fl}) where Fl

        # Define system matrices
        Z = Fl[1, 1, 0]
        T = Fl[
            1 0 0
            0 0 0
            0 0 0
        ]
        R = Matrix{Fl}(I, 3, 3)
        d = zero(Fl)
        c = zeros(Fl, 3)
        H = one(Fl)
        Q = Matrix{Fl}(I, 3, 3)

        system = LinearUnivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q)

        # Define hyperparameters names
        names = ["sigma2_ε", "sigma2_η", "sigma2_ω", "λ_c"]
        hyperparameters = HyperParameters{Fl}(names)

        return new(hyperparameters, system, Results{Fl}())
    end
end

# Obligatory methods
function default_filter(model::LocalLevelCycle)
    Fl = typeof_model_elements(model)
    a1 = zeros(Fl, num_states(model))
    P1 = Fl(1e6) .* Matrix{Fl}(I, num_states(model), num_states(model))
    steadystate_tol = Fl(1e-5)
    return UnivariateKalmanFilter(a1, P1, num_states(model), steadystate_tol)
end

function initial_hyperparameters!(model::LocalLevelCycle)
    Fl = typeof_model_elements(model)
    observed_variance = var(model.system.y[findall(!isnan, model.system.y)])
    initial_hyperparameters = Dict{String,Fl}(
        "sigma2_ε" => observed_variance,
        "sigma2_η" => observed_variance,
        "sigma2_ω" => one(Fl),
        # Durbin and Koopman (2012) comment possible values
        # in their book pp. 48
        "λ_c" => Fl(2 * pi / 12),
    )
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return model
end

function constrain_hyperparameters!(model::LocalLevelCycle)
    Fl = typeof_model_elements(model)
    constrain_variance!(model, "sigma2_ε")
    constrain_variance!(model, "sigma2_η")
    constrain_variance!(model, "sigma2_ω")
    # Durbin and Koopman (2012) comment possible values in their book pp. 48
    constrain_box!(model, "λ_c", Fl(2 * pi / 100), Fl(2 * pi / 1.5))
    return model
end

function unconstrain_hyperparameters!(model::LocalLevelCycle)
    Fl = typeof_model_elements(model)
    unconstrain_variance!(model, "sigma2_ε")
    unconstrain_variance!(model, "sigma2_η")
    unconstrain_variance!(model, "sigma2_ω")
    # Durbin and Koopman (2012) comment possible values in their book pp. 48
    unconstrain_box!(model, "λ_c", Fl(2 * pi / 100), Fl(2 * pi / 1.5))
    return model
end

function fill_model_system!(model::LocalLevelCycle)
    model.system.H = get_constrained_value(model, "sigma2_ε")
    model.system.Q[1] = get_constrained_value(model, "sigma2_η")
    model.system.Q[5] = get_constrained_value(model, "sigma2_ω")
    model.system.Q[9] = get_constrained_value(model, "sigma2_ω")
    model.system.T[5] = cos(get_constrained_value(model, "λ_c"))
    model.system.T[6] = -sin(get_constrained_value(model, "λ_c"))
    model.system.T[8] = sin(get_constrained_value(model, "λ_c"))
    model.system.T[9] = cos(get_constrained_value(model, "λ_c"))
    return model
end

function reinstantiate(::LocalLevelCycle, y::Vector{Fl}) where Fl
    return LocalLevelCycle(y)
end

has_exogenous(::LocalLevelCycle) = false


text:
"""
    NonlinearSciMLData

A package to automatically download various datasets related to physical experiments on
nonlinear systems.
"""
module NonlinearSciMLData

# Functions
export paper_doi, data_doi, describe, data_path

# Datasets
export NTMD_MSSP_2017

# Dependencies
using DataDeps: DataDeps, register, DataDep, unpack
using MAT: MAT

"""
    paper_doi(dataset)

Return the DOI of the paper associated with the data.
"""
function paper_doi end

"""
    data_doi(dataset)

Return the DOI of the data.
"""
function data_doi end

"""
    describe(dataset)

Return a short description of the dataset.
"""
function describe end

"""
    experiments(dataset)

Return a `Dict{Symbol,String}` of experiments within the dataset.
"""
function experiments end

"""
    experiment(dataset, label)

Return an array of measured inputs and outputs for a particular experiment. The array might
contain nested arrays to indicate where observations have been made under the same
conditions.
"""
function experiment end

"""
    data_path(dataset)

Return the directory containing the dataset.
"""
function data_path(dataset)
    datadep = string(@__MODULE__) * "." * string(nameof(dataset))
    return DataDeps.@datadep_str datadep
end

"""
    init_data()

Register all the data dependencies with DataDeps.jl.
"""
function init_data()
    for dataset in DATASETS
        register_data(dataset)
    end
end

"""
    register_data(dataset)

Register a particular data dependency.
"""
function register_data end

# Package initialisation
function __init__()
    return init_data()
end

# Actual datasets
include("datasets/NTMD_MSSP_2017.jl")
using ._NTMD_MSSP_2017

const DATASETS = [NTMD_MSSP_2017]

end

# Other datasets
# https://data.bris.ac.uk/data/dataset/15xkehnayh6st11gbc9dmeqg6v
# https://data.bris.ac.uk/data/dataset/1pfqddbxy8ry91n8j5c0cbeqxe
# https://data.bris.ac.uk/data/dataset/1jlp5gsztp3xr2hngpf0q72uii
# https://data.bris.ac.uk/data/dataset/1ttza2mrigjbt2mmbp6raqrxen


text:
"""
    createtimeset()

# 処理概要
時刻のデータセットを作成

# 引数
* なし

# 結果（戻り値）
* 時間の組み合わせ（文字列）
"""
function createtimeset()
  timemaster=Matrix{String}(undef,1,1441)
  HH=["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20","21","22","23"]
  MM=["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36","37","38","39","40","41","42","43","44","45","46","47","48","49","50","51","52","53","54","55","56","57","58","59"]

  timemaster[1,1441]=string("24:00")

  for i in 1:size(HH)[1]
      for j in 1:size(MM)[1]
          timemaster[1,60*(i-1)+j]=string(HH[i],":",MM[j])
      end
  end

  timessymbol=map(Symbol,timemaster)
  timesdate=map(dt->DateTime(dt,"HH:MM"),timemaster)

  return timemaster,timessymbol,timesdate
end

"""
    ContinuousOperation(x...)

# 処理概要
処理時間が連続しているか確認

# 引数
* `x`:

# 結果（戻り値）
* 連続している時間数
"""
function ContinuousOperation(x...)
  sigma=0
  for i in 1:length(x)-1
      sigma=sigma+x[i]*x[i+1]
  end
  return sigma
end


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

"""
    LinearPath(domain)

A linear path on spatial `domain`.
"""
struct LinearPath{D<:AbstractDomain} <: AbstractPath{D}
  domain::D
end

Base.iterate(p::LinearPath, state=1) = state > npoints(p.domain) ? nothing : (state, state + 1)

# ------------
# IO methods
# ------------
function Base.show(io::IO, path::LinearPath)
  print(io, "LinearPath")
end

function Base.show(io::IO, ::MIME"text/plain", path::LinearPath)
  println(io, path)
end


text:
using MultiQuad, Test, Unitful

@testset "quad" begin
    rtol = 1e-10

    xmin = 0
    xmax = 4
    func1(x) = x^2 * exp(-x)
    res = 2 - 26 / exp(4)
    int, err = quad(func1, xmin, xmax, rtol = rtol)
    @test isapprox(int, res, atol = err)
    @test typeof(quad(
        func1,
        xmin,
        xmax,
        rtol = rtol,
        method = :vegas,
    )[1]) == Float64
    @test typeof(quad(
        func1,
        xmin,
        xmax,
        rtol = rtol,
        method = :suave,
    )[1]) == Float64

    xmin2 = 0 * u"m"
    xmax2 = 4 * u"m"
    func2(x) = x^2
    res2 = 64 / 3
    int2, err2 = quad(func2, xmin2, xmax2, rtol = rtol)
    @test isapprox(int2, res2 * u"m^3", atol = err2)
    @test typeof(quad(
        func2,
        xmin2,
        xmax2,
        rtol = rtol,
        method = :vegas,
    )[1]) == typeof(1.0 * u"m^3")
    @test typeof(quad(
        func2,
        xmin2,
        xmax2,
        rtol = rtol,
        method = :suave,
    )[1]) == typeof(1.0 * u"m^3")

    @test_throws ErrorException quad(func1, 0, 4, method = :none)
end

@testset "dblquad" begin
    rtol = 1e-10

    xmin = 0
    xmax = 4
    ymin(x) = x^2
    ymax(x) = x^3

    func1(y, x) = y * x^3
    res1 = 241664 / 5
    int1, err1 = dblquad(func1, xmin, xmax, ymin, ymax, rtol = rtol)
    @test isapprox(int1, res1, atol = err1)
    @test typeof(dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :vegas,
    )[1]) == Float64
    @test typeof(dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :suave,
    )[1]) == Float64
    @test typeof(dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :divonne,
    )[1]) == Float64

    xmin2 = 0 * u"m"
    xmax2 = 4 * u"m"
    ymin2(x) = x^2 * u"m"
    ymax2(x) = x^3
    func2(y, x) = y * x^3

    int2, err2 = dblquad(func2, xmin2, xmax2, ymin2, ymax2, rtol = rtol)

    @test isapprox(int2, res1 * u"m^10", atol = err2)
    @test typeof(dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :vegas,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :suave,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(dblquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        rtol = rtol,
        method = :divonne,
    )[1]) == typeof(1.0 * u"m^10")

    @test_throws ErrorException dblquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        rtol = rtol,
        method = :none,
    )
end

@testset "tblquad" begin
    rtol = 1e-10

    xmin = 0
    xmax = 4
    ymin(x) = x^2
    ymax(x) = x^3
    zmin(x, y) = 3
    zmax(x, y) = 4

    func1(z, y, x) = y * x^3 * sin(z)
    res = 241664 / 5 * (cos(3) - cos(4))
    int1, err1 = tplquad(func1, xmin, xmax, ymin, ymax, zmin, zmax, rtol = rtol)
    @test isapprox(int1, res, atol = err1)
    @test typeof(tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :vegas,
    )[1]) == Float64
    @test typeof(tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :suave,
    )[1]) == Float64
    @test typeof(tplquad(
        func1,
        xmin,
        xmax,
        ymin,
        ymax,
        zmin,
        zmax,
        rtol = rtol,
        method = :divonne,
    )[1]) == Float64

    xmin2 = 0 * u"m"
    xmax2 = 4 * u"m"
    ymin2(x) = x^2 * u"m"
    ymax2(x) = x^3
    zmin2(x, y) = 3
    zmax2(x, y) = 4

    func2(z, y, x) = y * x^3 * sin(z)
    int2, err2 = tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
    )
    @test isapprox(int2, res * u"m^10", atol = err2)
    @test typeof(tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :vegas,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :suave,
    )[1]) == typeof(1.0 * u"m^10")
    @test typeof(tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :divonne,
    )[1]) == typeof(1.0 * u"m^10")

    @test_throws ErrorException tplquad(
        func2,
        xmin2,
        xmax2,
        ymin2,
        ymax2,
        zmin2,
        zmax2,
        rtol = rtol,
        method = :none,
    )[1]
end


text:
# Renyi entanglement entropy of Free fermions in 1D.

push!(LOAD_PATH, joinpath(dirname(@__FILE__), "src"))
using FreeFermionsOEE

using ArgParse
using Printf
using LinearAlgebra

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "L"
        help = "number of sites"
        arg_type = Int
        required = true
    "N"
        help = "number of particles"
        arg_type = Int
        required = true
    "--out"
        metavar = "FILE"
        help = "path to output file"
        required = true
end

add_arg_group(s, "BH parameters")
@add_arg_table s begin
    "--l-min"
        metavar = "l"
        help = "minimum lA"
        arg_type = Int
        default = 1
    "--l-max"
        metavar = "l"
        help = "maximum lA"
        arg_type = Int
        default = 1
    "--l-step"
        metavar = "l"
        help = "lA step"
        arg_type = Int
        default = 1
    "--l-logstep"
        metavar = "Delta"
        help = " log(lA) step"
        arg_type = Float64
        default = 0.1
    "--l-num"
        metavar = "n"
        help = "number of lA"
        arg_type = Int
        default = 1
    "--l-log"
        help = "use logarithmic scale for l"
        action = :store_true
    "--precision"
        metavar = "precision"
        help = "precision"
        arg_type = Int
        default = 54
    "--alpha"
        metavar = "alpha"
        help = "alpha value"
        arg_type = Float64
        default = 2.0
    "--eigs"
        help = "report the spectrum of the projected correlation matrix."
        action = :store_true
    "--prob"
        help = "report the particle-number probability distribution."
        action = :store_true
end

function main()

   c = parsed_args = parse_args(ARGS, s, as_symbols=true)
   # Number of sites
   L = c[:L]
   # Number of particles
   N = c[:N]
   # Output file
   output = c[:out]
   # minimum lA
   lA_min = c[:l_min]
   #precision
   precision=c[:precision]
   alpha=c[:alpha]
   eigenvalues= c[:eigs]
   probabilities= c[:prob]

   if c[:l_log]
      l_cut=trunc(Int,1.0/(exp(c[:l_logstep])-1.0))+1
      if lA_min<l_cut
         loglnum=trunc(Int,log(c[:l_max]*1.0/l_cut)/c[:l_logstep])
         l_range = zeros(Int, l_cut-lA_min+1+ loglnum)
         for i=1: l_cut-lA_min+1
            l_range[i]=i+ lA_min-1
         end
         for i=1:loglnum
            l_range[i+l_cut-lA_min+1]= trunc(Int,l_cut*exp(i*c[:l_logstep]))
         end
         lA_max=l_range[l_cut+ loglnum-lA_min]
      else
         loglnum=trunc(Int,log(c[:l_max]*1.0/lA_min)/c[:l_logstep])
         l_range = zeros(Int, loglnum+1)
         for i=1: loglnum+1
            l_range[i]= trunc(Int,lA_min*exp((i-1)*c[:l_logstep]))
         end
         lA_max=l_range[loglnum+1]
      end
   else
      lA_max=(lA_min+c[:l_step]*(c[:l_num]-1))
      l_range = lA_min:c[:l_step]: lA_max
   end
#_______________________
   open(output, "w") do f
      write(f, "# L=$(L), N=$(N), precision=$(precision), alpha=$(alpha) \n")
      write(f, @sprintf "#%7s%24s%24s%24s%24s%24s%24s%24s%24s%24s%24s\n" "l"    "S1"    "S1_OP"    "S_alpha"    "S_alpha_OP(1)"    "sigma2_n" "H_alpha" "navge_alpha" "S_alpha_OP(3)" "S_alpha_OP(4)" "S_alpha_OP(5)")

   for l in l_range
      # ======================= Correlation matrix =============================
      Corr_Matrix = zeros(Float64, l, l)
      phase_factor= Float64
      even_odd=N % 2
      phase_factor=2*pi/L
      for i=1: l
         for j=i+1: l
            Corr_Matrix[i,j]=sin(N*phase_factor/2*(i-j))/sin(phase_factor/2*(i-j))/L
            Corr_Matrix[j,i]=Corr_Matrix[i,j]
         end
         Corr_Matrix[i,i]=N/L
      end
	# ==================================End===================================

        # Create C_A
        nu=eigvals!(Symmetric(Corr_Matrix))
        nu_increasing= zeros(Float64, l)
        nu_increasing[l:-1:1]=nu[1:l]
        navge_alpha=0.0
        for k=1:l
         if nu[k]>=0 && nu[k]<=1
             navge_alpha+=nu[k]^alpha/(nu[k]^alpha+(1-nu[k])^alpha)
         elseif nu[k]>1
             navge_alpha+=1.0
         end
        end

 # Calculate spatial & operational entanglement entropy
        s1_spatial, s1_operational, salpha_spatial, salpha_operational, sigma2_n, LogSumPNalpha, H4, H5, H6 = operational_entanglement(L, N , l, precision, alpha, nu_increasing, eigenvalues, probabilities)
        write(f, @sprintf "%8s%24.12E%24.12E%24.12E%24.12E%24.12E%24.12E%24.12E%24.12E%24.12E%24.12E\n" "$(l)"  s1_spatial s1_operational salpha_spatial salpha_operational sigma2_n LogSumPNalpha navge_alpha  H4 H5 H6)

        flush(f)
   end
   end
end

main()