text:
"""
 Holds the tableau of a weak implicit Runge-Kutta method.

   According to Wang, Hong, Xu, "Construction of Symplectic Runge-Kutta Methods for Stochastic Hamiltonian Systems",
   Commun. Comput. Phys. 21(1), 2017
"""
struct TableauWIRK{T} <: AbstractTableauIRK{T}
    name::Symbol
    s::Int
    qdrift0::CoefficientsRK{T}
    qdrift1::CoefficientsRK{T}
    qdiff0::CoefficientsRK{T}
    qdiff1::CoefficientsRK{T}
    qdiff3::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift and qdiff are understood as the classical orders of these methods.
    #
    # qdrift0, qdrift1 correspond to A0, A1 in the paper
    # qdiff0, qdiff1, qdiff3 correspond to B0, B1, B3
    # qdrift0.b = alpha
    # qdiff0.b  = beta
    # qdrift0.c = c0
    # qdrift1.c = c1

    function TableauWIRK{T}(name, s, qdrift0, qdrift1, qdiff0, qdiff1, qdiff3) where {T}
        # THE COMMENTED OUT PART WAS FOR TableauFIRK. MAY IMPLEMENT SOMETHING
        # SIMILAR FOR TableauSFIRK LATER.

        # if (q.s > 1 && istrilstrict(q.a)) || (q.s==1 && q.a[1,1] == 0)
        #     warn("Initializing TableauFIRK with explicit tableau ", q.name, ".\n",
        #          "You might want to use TableauERK instead.")
        # elseif q.s > 1 && istril(q.a)
        #     warn("Initializing TableauFIRK with diagonally implicit tableau ", q.name, ".\n",
        #          "You might want to use TableauDIRK instead.")
        # end

        @assert s == qdrift0.s == qdrift1.s == qdiff0.s == qdiff1.s == qdiff3.s

        new(name, s, qdrift0, qdrift1, qdiff0, qdiff1, qdiff3)
    end
end

function TableauWIRK(name::Symbol, qdrift0::CoefficientsRK{T}, qdrift1::CoefficientsRK{T}, qdiff0::CoefficientsRK{T}, qdiff1::CoefficientsRK{T}, qdiff3::CoefficientsRK{T}) where {T}
    TableauWIRK{T}(name, qdrift0.s, qdrift0, qdrift1, qdiff0, qdiff1, qdiff3)
end

function TableauWIRK(name::Symbol, A0::Matrix{T}, A1::Matrix{T},
                                   B0::Matrix{T}, B1::Matrix{T}, B3::Matrix{T},
                                    α::Vector{T}, β1::Vector{T},
                                   c0::Vector{T}, c1::Vector{T} ) where {T}
    TableauWIRK(name, CoefficientsRK(name, 0, A0, α, c0), CoefficientsRK(name, 0, A1, α, c1),
                         CoefficientsRK(name, 0, B0, β1, c0), CoefficientsRK(name, 0, B1, β1, c1), CoefficientsRK(name, 0, B3, β1, c1))
end

"Parameters for right-hand side function of weak implicit Runge-Kutta methods."
mutable struct ParametersWIRK{DT, TT, ET <: SDE{DT,TT}, D, M, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauWIRK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}

    t::TT
    q::Vector{DT}
end

function ParametersWIRK(equ::ET, tab::TableauWIRK{TT}, Δt::TT, ΔW::Vector{DT}, ΔZ::Vector{DT}) where {DT, TT, ET <: SDE{DT,TT}}
    @assert equ.m == length(ΔW) == length(ΔZ)
    ParametersWIRK{DT, TT, ET, equ.d, equ.m, tab.s}(equ, tab, Δt, ΔW, ΔZ, 0, zeros(DT, equ.d))
end

struct IntegratorCacheWIRK{DT}
    # Structure for holding the internal stages Q0, and Q1 the values of the drift vector
    # and the diffusion matrix evaluated at the internal stages V=v(Q0), B=B(Q1),
    # and the increments Y = Δt*a_drift*v(Q) + a_diff*B(Q)*ΔW
    Q0::Vector{Vector{DT}}     # Q0[i][k]   - the k-th component of the internal stage Q^(0)_i
    Q1::Vector{Matrix{DT}}     # Q1[i][k,l] - the k-th component of the internal stage Q^(l)_i
    V ::Vector{Vector{DT}}     # V [i][k]   - the k-th component of the drift vector v(Q0[i][:])
    B ::Vector{Matrix{DT}}     # B [i][:,:] - the diffusion matrix at the stage i such that the l-th column B[i][:,l] is evaluated at Q^(l)_i
    Y0::Vector{Vector{DT}}     # Y0[i][:]   - the increment of the internal stage Q0[i][:]
    Y1::Vector{Matrix{DT}}     # Y1[i][:,l] - the increment of the internal stage Q1[i][:,l]

    v::Vector{DT}
    b::Matrix{DT}
    y::Vector{DT}

    Δy::Vector{DT}

    function IntegratorCacheWIRK{DT}(D, M, S) where {DT}
        # create internal stage vectors
        Q0 = create_internal_stage_vector(DT, D, S)
        Q1 = create_internal_stage_vector(DT, D, M, S)
        V  = create_internal_stage_vector(DT, D, S)
        B  = create_internal_stage_vector(DT, D, M, S)
        Y0 = create_internal_stage_vector(DT, D, S)
        Y1 = create_internal_stage_vector(DT, D, M, S)

        # create velocity and update vector
        v = zeros(DT,D)
        b = zeros(DT,D,M)
        y = zeros(DT,D)

        # create temporary vectors
        Δy  = zeros(DT,M)

        new(Q0, Q1, V, B, Y0, Y1, v, b, y, Δy)
    end
end

"Stochastic implicit Runge-Kutta integrator."
struct IntegratorWIRK{DT, TT, PT <: ParametersWIRK{DT,TT},
                              ST <: NonlinearSolver{DT}, N} <: StochasticIntegrator{DT,TT}
    params::PT
    solver::ST
    cache::IntegratorCacheWIRK{DT}
end

function IntegratorWIRK(equation::SDE{DT,TT,VT,BT,N}, tableau::TableauWIRK{TT}, Δt::TT) where {DT,TT,VT,BT,N}
    D = equation.d
    M = equation.m
    NS= max(equation.ns,equation.ni)
    S = tableau.s

    # create params
    params = ParametersWIRK(equation, tableau, Δt, zeros(DT,M), zeros(DT,M))

    # create solver
    solver = create_nonlinear_solver(DT, D*S*(M+1), params)

    # create cache for internal stage vectors and update vectors
    cache = IntegratorCacheWIRK{DT}(D, M, S)

    # create integrator
    IntegratorWIRK{DT, TT, typeof(params), typeof(solver), N}(params, solver, cache)
end

equation(integrator::IntegratorWIRK) = integrator.params.equ
timestep(integrator::IntegratorWIRK) = integrator.params.Δt
tableau(integrator::IntegratorWIRK) = integrator.params.tab
noisedims(integrator::IntegratorWIRK) = integrator.params.equ.m
Base.eltype(integrator::IntegratorWIRK{DT, TT, PT, ST, N}) where {DT, TT, PT, ST, N} = DT

function update_params!(int::IntegratorWIRK, sol::AtomicSolutionSDE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.ΔW .= sol.ΔW
    int.params.ΔZ .= sol.ΔZ
end

"""
Unpacks the data stored in
x = (Y0[1][1], Y0[1][2]], ... Y0[1][D], ... Y0[S][D], Y1[1][1,1], Y1[1][2,1], ... Y1[1][D,1], Y1[1][1,2], Y1[1][2,2], ... Y1[1][D,2], ... Y1[S][D,M]  )
into Y0::Vector{Vector} and Y1::Vector{Matrix}, calculates the internal stages Q0 and Q1, the values of the RHS
of the SDE ( v(Q0) and B(Q1) ), and assigns them to V and B.
Unlike for FIRK, here Y = Δt a v(Q) + ̃a B(Q) ΔW
"""
@generated function compute_stages!(x::Vector{ST}, Q0::Vector{Vector{ST}}, Q1::Vector{Matrix{ST}}, V::Vector{Vector{ST}},
                                                   B::Vector{Matrix{ST}}, Y0::Vector{Vector{ST}}, Y1::Vector{Matrix{ST}},
                                                   params::ParametersWIRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    tQ::Vector{ST} = zeros(ST,D)
    tV::Vector{ST} = zeros(ST,D)
    tB::Vector{ST} = zeros(ST,D)    # This was tB::Matrix{ST} = zeros(ST,D,M) for SIRK, but here we calculate B column by column.
                                    # tB is actually superfluous, could just use tV; keeping tB for clarity though

    quote
        local tᵢ::TT

        @assert S == length(Q0) == length(Q1) == length(V) == length(B)

        # copy x to Y0 and Y1, and calculate Q0 and Q1
        for i in eachindex(Q0)

            @assert D == length(Q0[i]) == size(Q1[i],1) == length(V[i]) == size(B[i],1)
            @assert M == size(Q1[i],2) == size(B[i],2)

            for k in eachindex(Q0[i])
                Y0[i][k] = x[D*(i-1)+k]
                Q0[i][k] = params.q[k] + Y0[i][k]
            end
        end

        for i in eachindex(Q1)
            for l in 1:size(Q1[i],2)
                for k in 1:size(Q1[i],1)
                    Y1[i][k,l] = x[ D*S + D*M*(i-1) + D*(l-1) + k ]
                    Q1[i][k,l] = params.q[k] + Y1[i][k,l]
                end
            end
        end

        # compute V = v(Q0) and B=B(Q1)
        for i in 1:S
            # time point with the coefficient c0
            tᵢ = params.t + params.Δt * params.tab.qdrift0.c[i]
            # calculates v(t,tQ) and assigns to the i-th column of V
            params.equ.v(tᵢ, Q0[i], V[i])

            # time point with the coefficient c1
            tᵢ = params.t + params.Δt * params.tab.qdrift1.c[i]

            for l=1:M
                # copies Q1[i][:,l] to the vector tQ
                simd_copy_xy_first!($tQ, Q1[i], l)
                # calculates the l-th column of B(t,tQ) and assigns to the vector tB
                params.equ.B(tᵢ, $tQ, $tB, col=l)
                # copies the vector tB to B[i][:,l]
                simd_copy_yx_first!($tB, B[i], l)
            end
        end
    end
end

"Compute stages of weak implicit Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersWIRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    cache = IntegratorCacheWIRK{ST}(D, M, S)

    quote
        compute_stages!(x, $cache.Q0, $cache.Q1, $cache.V, $cache.B, $cache.Y0, $cache.Y1, params)

        local y1::ST
        local y2::ST

        # compute b = - (Y-AV)
        # the terms corresponding to Y0
        for i in 1:S
            for k in 1:D
                y1 = 0
                y2 = 0
                for j in 1:S
                    y1 += params.tab.qdrift0.a[i,j] * $cache.V[j][k] * params.Δt
                    y2 += params.tab.qdrift0.â[i,j] * $cache.V[j][k] * params.Δt
                    for l in 1:M
                        y1 += params.tab.qdiff0.a[i,j] * $cache.B[j][k,l] * params.ΔW[l]
                        y2 += params.tab.qdiff0.â[i,j] * $cache.B[j][k,l] * params.ΔW[l]
                    end
                end
                b[D*(i-1)+k] = - $cache.Y0[i][k] + (y1 + y2)
            end
        end

        # compute b = - (Y-AV)
        # the terms corresponding to Y1
        for i in 1:S
            for l in 1:M
                for k in 1:D
                    y1 = 0
                    y2 = 0
                    for j in 1:S
                        # The drift terms
                        y1 += params.tab.qdrift1.a[i,j] * $cache.V[j][k] * params.Δt
                        y2 += params.tab.qdrift1.â[i,j] * $cache.V[j][k] * params.Δt

                        # The noise terms, calculated either with the B1 or B3 tableau
                        for noise_idx in 1:M
                            if noise_idx==l
                                y1 += params.tab.qdiff1.a[i,j] * $cache.B[j][k,noise_idx] * params.ΔW[noise_idx]
                                y2 += params.tab.qdiff1.â[i,j] * $cache.B[j][k,noise_idx] * params.ΔW[noise_idx]
                            else
                                y1 += params.tab.qdiff3.a[i,j] * $cache.B[j][k,noise_idx] * params.ΔW[noise_idx]
                                y2 += params.tab.qdiff3.â[i,j] * $cache.B[j][k,noise_idx] * params.ΔW[noise_idx]
                            end
                        end
                    end
                    b[D*S + D*M*(i-1) + D*(l-1) + k] = - $cache.Y1[i][k,l] + (y1 + y2)
                end
            end
        end
    end
end

"""
This function computes initial guesses for Y and assigns them to int.solver.x
For WIRK we are NOT IMPLEMENTING an InitialGuess.
"""
function initial_guess!(int::IntegratorWIRK{DT,TT}, sol::AtomicSolutionSDE{DT,TT}) where {DT,TT}
    # The simplest initial guess for Y is 0
    int.solver.x .= 0

    # Using an explicit integrator to predict the next step's value (like in SIRK)
    # does not seem to be a good idea here, because the integrators are convergent
    # in the weak sense only, and there is no guarantee that the explicit integrator
    # will produce anything close to the desired solution...
end

"""
Integrate SDE with a stochastic implicit Runge-Kutta integrator.
  Integrating the m-th sample path
"""
function integrate_step!(int::IntegratorWIRK{DT,TT}, sol::AtomicSolutionSDE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int, sol)

    # compute initial guess and assign to int.solver.x
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute the drift vector field and the diffusion matrix at internal stages
    compute_stages!(int.solver.x, int.cache.Q0, int.cache.Q1, int.cache.V, int.cache.B, int.cache.Y0, int.cache.Y1, int.params)

    # compute final update (same update function as for SIRK)
    update_solution!(sol.q, int.cache.V, int.cache.B, int.params.tab.qdrift0.b, int.params.tab.qdrift0.b̂, int.params.tab.qdiff0.b, int.params.tab.qdiff0.b̂, int.params.Δt, int.params.ΔW, int.cache.Δy)
end


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

# ----------------------------------------------------------------------------------- #
# Function: generate_problem_dictionary
# Description: Holds simulation and model parameters as key => value pairs in a Julia Dict()
# Generated on: 2018-03-15T00:00:56.939
#
# Output arguments:
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model and simulation parameters as key => value pairs
# ----------------------------------------------------------------------------------- #
function generate_problem_dictionary()

	# Load the stoichiometric network from disk -
	path_to_network_file = joinpath(_PATH_TO_CONFIG,"Network.dat")
	stoichiometric_matrix = readdlm(path_to_network_file);

	# What is the system dimension? -
	(number_of_species,number_of_reactions) = size(stoichiometric_matrix)

	# Setup the flux bounds array -
	flux_bounds_array = zeros(number_of_reactions,2)
	# TODO: update the flux_bounds_array for each reaction in your network, col 1 => lower bound, col 2 => upper bound, each row is a reaction
	#flux in units of micromole/gDW-sec
	flux_bounds_array[1,1]=0
	flux_bounds_array[1,2]=2.03 #Vmax for v1
	flux_bounds_array[2,1]=0
	flux_bounds_array[2,2]=0.345 #Vmax for v2
	flux_bounds_array[3,1]=0
	flux_bounds_array[3,2]=2.49 #Vmax for v3
	flux_bounds_array[4,1]=0
	flux_bounds_array[4,2]=0.881 #Vmax for v4
	flux_bounds_array[5,1]=-0.137 #Vmax for v5 (reverse rxn)
	flux_bounds_array[5,2]=0.137 #Vmax for v5 (forward rxn)

	flux_bounds_array[6,1]=0
	flux_bounds_array[6,2]=2.78 #upper bound for b1, in micromole/gDW-sec
	flux_bounds_array[7,1]=0
	flux_bounds_array[7,2]=2.78 #upper bound for b2, in micromole/gDW-sec
	flux_bounds_array[8,1]=0
	flux_bounds_array[8,2]=2.78 #upper bound for b3
	flux_bounds_array[9,1]=0
	flux_bounds_array[9,2]=2.78 #upper bound for b4
# additional b rxns:
	flux_bounds_array[10,1]=0
	flux_bounds_array[10,2]=2.78 #upper bound for b5
	flux_bounds_array[11,1]=0
	flux_bounds_array[11,2]=2.78 #upper bound for b6
	flux_bounds_array[12,1]=0
	flux_bounds_array[12,2]=2.78 #upper bound for b7
	flux_bounds_array[13,1]=0
	flux_bounds_array[13,2]=2.78 #upper bound for b8
	flux_bounds_array[14,1]=0
	flux_bounds_array[14,2]=2.78 #upper bound for b9
	flux_bounds_array[15,1]=0
	flux_bounds_array[15,2]=2.78 #upper bound for b10
	flux_bounds_array[16,1]=0
	flux_bounds_array[16,2]=2.78 #upper bound for b11
	flux_bounds_array[17,1]=0
	flux_bounds_array[17,2]=2.78 #upper bound for b12
	flux_bounds_array[18,1]=0
	flux_bounds_array[18,2]=2.78 #upper bound for b13
	flux_bounds_array[19,1]=0
	flux_bounds_array[19,2]=2.78 #upper bound for b14
	flux_bounds_array[20,1]=0
	flux_bounds_array[20,2]=2.78 #upper bound for b15
	flux_bounds_array[21,1]=0
	flux_bounds_array[21,2]=2.78 #upper bound for b16

	# Setup default species bounds array -
	species_bounds_array = zeros(number_of_species,2)
	# TODO: NOTE - we by default assume Sv = 0, so species_bounds_array should be a M x 2 array of zeros
	# TODO: however, if you formulate the problem differently you *may* need to change this
	# so...I don't need to change anything here, right?

	# Min/Max flag - default is minimum -
	is_minimum_flag = true

	# Setup the objective coefficient array -
	objective_coefficient_array = zeros(number_of_reactions)
	# TODO: update me to maximize Urea production (Urea leaving the virtual box)
	# what is the physical significance of this? Did I code it corectly?
	# TODO: if is_minimum_flag = true => put a -1 in the index for Urea export (reaction b4)

	if is_minimum_flag == true
		objective_coefficient_array[9]= -1.0
	end

# this is leaving all other coefficients as zero...is that what I want to be doing??

	# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{String,Any}()
	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
	data_dictionary["objective_coefficient_array"] = objective_coefficient_array
	data_dictionary["flux_bounds_array"] = flux_bounds_array;
	data_dictionary["species_bounds_array"] = species_bounds_array
	data_dictionary["is_minimum_flag"] = is_minimum_flag
	data_dictionary["number_of_species"] = number_of_species
	data_dictionary["number_of_reactions"] = number_of_reactions
	# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
	return data_dictionary
end

ext: jl


text:
using IncrementalInference
using Test

##

@testset "test Sphere1D" begin

##

fg = initfg()
getSolverParams(fg).useMsgLikelihoods = true

addVariable!.(fg, [Symbol("x$i") for i=0:4], Sphere1)
addFactor!(fg, [:x0], PriorSphere1(Normal(0.0,0.1)))
map(i->addFactor!(fg, [Symbol("x$i"),Symbol("x$(i+1)")], Sphere1Sphere1(Normal(1.0, 0.1))), 0:3)

solveTree!(fg);

##

sppes = map(var->getPPE(var).suggested[1], sortDFG(getVariables(fg),by=getLabel))

gt = rem2pi.(collect(0:4), RoundNearest)

@show sppes
@test all(isapprox.(sppes, gt, atol=0.35))

# test packing converters also
d = "/tmp/caesar/random/testfg"
saveDFG(d,fg)
lfg = loadDFG(d)
Base.rm(d*".tar.gz")
# check loaded fg for all variable and factors
@test issetequal(ls(fg), ls(lfg))
@test issetequal(lsf(fg), lsf(lfg))

##

end


text:
function blockranges(nblocks, total_len)
    rem = total_len % nblocks
    main_len = div(total_len, nblocks)

    starts=Int[1]
    ends=Int[]
    for ii in 1:nblocks
        len = main_len
        if rem>0
            len+=1
            rem-=1
        end
        push!(ends, starts[end]+len-1)
        push!(starts, ends[end] + 1)
    end
    @assert ends[end] == total_len
    starts[1:end-1], ends
end

function _threadedsort(data::Vector)
    starts, ends = blockranges(nthreads(), length(data))

    # Sort each block
    @threads for (ss, ee) in collect(zip(starts, ends))
        @inbounds sort!(@view(data[ss:ee]))
    end

    # Go through each sorted block taking out the smallest item and putting it in the new array
    # This code could maybe be optimised. see https://stackoverflow.com/a/22057372/179081
    # ret = similar(data) # main bit of allocation right here. avoiding it seems expensive.
    # # Need to not overwrite data we haven't read yet
    # @inbounds for ii in eachindex(ret)
    #     minblock_id = 1
    #     ret[ii]=data[starts[1]]
    #     @inbounds for blockid in 2:endof(starts) # findmin allocates a lot for some reason, so do the find by hand. (maybe use findmin! ?)
    #         ele = data[starts[blockid]]
    #         if ret[ii] > ele
    #             ret[ii] = ele
    #             minblock_id = blockid
    #         end
    #     end
    #     starts[minblock_id]+=1 # move the start point forward
    #     if starts[minblock_id] > ends[minblock_id]
    #         deleteat!(starts, minblock_id)
    #         deleteat!(ends, minblock_id)
    #     end
    # end
    # return ret
    # data.=ret  # copy back into orignal as we said we would do it inplace
    # return data
end

# parallel sorting algorithms based on https://stackoverflow.com/questions/47235390/how-can-i-parallelize-sorting/47235391#47235391
# the problem with this sort is that it only performs better if the number of unique elements is large
# see below for benchmarking code
threadedsort(a) = _threadedsort(copy(a))

if false
# a = rand(Int, 100_000_000);
# tic()
# @time threadedsort(a);
# toc()

# @time sort(a, alg=RadixSort);

    a = rand(1:40_000_000, 100_000_000);
    @time threadedsort(a);
    @time sort(a, alg=RadixSort);
    @time sort(a);
end


text:
Threads.nthreads()

a = zeros(12)

Threads.@threads for i = 1:12
    a[i] = Threads.threadid()
end
a

@macroexpand Threads.@threads for i = 1:12
    a[i] = Threads.threadid()
end

a = rand(10 ^ 6)

@timev begin
    for i = 1:length(a)
        a[i] = Threads.threadid()
    end
end

@timev begin
    Threads.@threads for i = 1:length(a)
        a[i] = Threads.threadid()
    end
end

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


text:
__precompile__()

module Weber
using TimedSound
using Unitful
using Juno
using Base.Iterators
using MacroTools
using Lazy: @>, @>>, @_
export resize_cache!, @>, @>>, @_

export sound, tone, noise, silence, harmonic_complex, irn, audible, highpass,
  lowpass, bandpass, bandstop, ramp, rampon, rampoff, fadeto, attenuate, mix,
  mult, envelope, duration, nchannels, nsamples, audiofn, leftright, left,
  right, play, setup_sound , playable, resample, stop, samplerate,
  current_sound_latency, pause_sounds, resume_sounds, ..,
  ms, s, kHz, Hz, samples

try
  @assert Sys.WORD_SIZE == 64
catch
  error("Weber can only be run as a 64-bit program. Please use a 64-bit ",
        "implementation of Julia.")
end

old = pwd()
try
  cd(Pkg.dir("Weber"))

  suffix = (success(`git diff-index HEAD --quiet`) ? "" : "-dirty")
  if !isempty(suffix)
    warn("Source files in $(Pkg.dir("Weber")) have been modified "*
                  "without being committed to git. Your experiment will not "*
                  "be reproduceable.")
  end
  global const version =
    convert(VersionNumber,chomp(readstring(`git describe --match 'v*' --tags`))*suffix)
catch
  try
    global const version = Pkg.installed("Weber")
    if !isempty(version.build)
      warn("Source files do not correspond to an official release "*
                    "of Weber. Your experiment will not be reproducable. "*
                    "Consider installing git and adding it to your PATH to "*
                    "record a more precise version number.")
    end
  catch
    warn("The Weber version number could not be determined. "*
         "Your experiment will not be reproducable. "*
         "It is recommended that you install Weber via Pkg.add(\"Weber\").")
  end
finally
  cd(old)
end

# load binary library dependencies
depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsjl)
  include(depsjl)
else
  error("Weber not properly installed. "*
        "Please run\nPkg.build(\"Weber\")")
end

# setup error reporting functions (these are the only calls to SDL that occur
# all that often, so they're the only calls I've wrapped directly).
SDL_GetError() = unsafe_string(ccall((:SDL_GetError,weber_SDL2),Cstring,()))
Mix_GetError = SDL_GetError
TTF_GetError = SDL_GetError

# this is a simple function for accessing aribtrary offsets in memory as any
# bitstype you want... this is used to read from a c union by determining the
# offset of various fields in the data using offsetof(struct,field) in c and
# then using that offset to access the memory in julia. SDL's
# core event type (SDL_Event) is a c union.
function at{T}(x::Ptr{Void},::Type{T},offset)
  unsafe_wrap(Array,reinterpret(Ptr{T},x + offset),1)[1]
end

import FileIO: load, save
export load, save

const sdl_is_setup = Array{Bool}()
sdl_is_setup[] = false

include(joinpath(@__DIR__,"timing.jl"))
include(joinpath(@__DIR__,"video.jl"))

function resize_cache!(size)
  resize!(image_cache,size)
  resize!(convert_cache,size)
end

include(joinpath(@__DIR__,"types.jl"))
include(joinpath(@__DIR__,"sound_hooks.jl"))
include(joinpath(@__DIR__,"event.jl"))
include(joinpath(@__DIR__,"trial.jl"))
include(joinpath(@__DIR__,"experiment.jl"))

include(joinpath(@__DIR__,"primitives.jl"))
include(joinpath(@__DIR__,"helpers.jl"))
include(joinpath(@__DIR__,"adaptive.jl"))

include(joinpath(@__DIR__,"precompile.jl"))

include(joinpath(@__DIR__,"extension_macro.jl"))
include(joinpath(@__DIR__,"extensions.jl"))

const localunits = Unitful.basefactors
const localpromotion = Unitful.promotion
function __init__()
  _precompile_()
  merge!(Unitful.basefactors,localunits)
  merge!(Unitful.promotion, localpromotion)
end

end


text:
const CANDLE_NAMES = ["unix_time", "low", "high", "open", "close", "volume"]
const GRANULARITY = [60, 300, 900, 3600, 21600, 86400]

function _candles_url(id::String, startDate::DateTime, endDate::DateTime, granularity::Int64)
  ENDPOINT * "products/$(id)/candles?start=$(startDate)&&end=$(endDate)&&granularity=$(granularity)"
end

function candles(id::String, startDate::DateTime, endDate::DateTime, granularity::Int64)

  url = _candles_url(id, startDate, endDate, granularity)

  res = HTTP.get(url, status_execption=false)
  rawData = JSON.parse(String(res.body))
  if(!isa(rawData, Array))
    error(rawData["message"])
  end

  df = vcat(DataFrame.(Dict.(zip.([CANDLE_NAMES], rawData)))...)
  df.time = unix2datetime.(df.unix_time)
  df
end


text:
using Base.MathConstants

"""
cf is a closed form approximation for the Lognormal characteristic function.
See the following for details: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.720.1797&rep=rep1&type=pdf
"""
function cf(d::LogNormal, θ)
    θ = -θ*im; σ = d.σ; μ = d.μ
    a = lambertw(θ*e^μ*σ^2)^2
    b = 2*lambertw(θ*e^μ*σ^2)
    d = sqrt(1+lambertw(θ*e^μ*σ^2))
    return exp(-(a + b)/(2*σ^2))/d
end

mutable struct Power <: ContinuousUnivariateDistribution
    L::Float64
    U::Float64
    d::Float64
end

function cf(dist::Power, t)
    L = dist.L; U = dist.U; d=dist.d
    Γ = inc_gamma_lower
    x1 = (-t+0im)^(1/d)*exp(1im*π/(2*d))*gamma(-1/d)*Γ(exp(1im*π/d),L^(-d)*t*exp(1im*π/2))/(d*gamma(1-1/d))
    x2 = (-t+0im)^(1/d)*exp(1im*π/(2*d))*gamma(-1/d)*Γ(exp(1im*π/d),U^(-d)*t*exp(1im*π/2))/(d*gamma(1-1/d))
    return x1 + x2
end


text:
using DCEGM
using DelimitedFiles
using Test

@testset "Running DCEGM tests" begin

    include("basics.jl")
    include("linetests.jl")
    include("Envtests.jl")
    include("F_test.jl")

end


text:
## 18.jl
# 各行を3コラム目の数値の降順にソート

f = open("./dat/hightemp.txt")
tsv = [split(chomp(line), "\t") for line in readlines(f)]
sort_tsv = sort(tsv, by= x -> x[3]) #　昇順
for x in reverse(sort_tsv) #　降順
  println(join(x," "))
end
close(f)


text:
using BenchmarkTools
using Distributed
addprocs()
@everywhere using TrajectoryOptimization
@everywhere using SharedArrays

# Set up problem
model, obj0 = Dynamics.cartpole_analytical
n,m = model.n, model.m

obj = copy(obj0)
obj.x0 = [0;0;0;0.]
obj.xf = [0.5;pi;0;0]
obj.tf = 2.
u_bnd = 50
x_bnd = [0.6,Inf,Inf,Inf]
obj_con = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd, x_min=-x_bnd, x_max=x_bnd)
obj_con = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd)
obj_con = to_static(obj_con)
obj = to_static(obj)
dt = 0.1

# Initialize trajectory
solver = Solver(model,obj,dt=dt,integration=:rk3_foh)
n,m,N = get_sizes(solver)
U0 = ones(1,N)*1
X0 = line_trajectory(obj.x0, obj.xf, N)

solver = Solver(model,obj,dt=dt,integration=:rk3)
res = UnconstrainedVectorResults(n,m,N)
reS = UnconstrainedStaticResults(n,m,N)

rollout!(res,solver)
rollout!(reS,solver)
J_prev2 = cost(solver, res)
J_prev3 = cost(solver, reS)

@btime TrajectoryOptimization.update_jacobians!(res, solver)
@btime TrajectoryOptimization.update_jacobians!(reS, solver)

X,U = res.X,res.U
fx,fu = res.fx, res.fu

fxu = [zeros(n,n+m) for i = 1:N-1]
fx_view = [view(F,1:n,1:n) for F in fxu]
fu_view = [view(F,1:n,n.+(1:m)) for F in fxu]

function calc_dyn_jacob(fx,fu,X,U)
    N = length(X)
    for k = 1:N-1
        res.fx[k], res.fu[k] = solver.Fd(res.X[k], res.U[k])
    end
    return nothing
end

function calc_dyn_jacob_map(fx,fu,X,U)
    N = solver.N
    jacob = map((x,u)->solver.Fd(x,u),res.X,res.U)
    for k = 1:N-1
        res.fx[k],res.fu[k] = jacob[k]
    end
end

function calc_dyn_jacob_map!(fxu,X,U)
    N = solver.N
    map!((x,u)->Fd(x,u),fxu,res.X,res.U)
    fx_view = [view(F,1:n,1:n) for F in fxu]
    fu_view = [view(F,1:n,n.+(1:m)) for F in fxu]
    return fx_view,fu_view
end

Xrand = rand(n,N)
Urand = rand(m,N)

copyto!(res.X,Xrand)
copyto!(res.U,Urand)

k = rand(1:N-1)
calc_dyn_jacob(fx,fu,X,U)
A,B = fx[k],fu[k]
calc_dyn_jacob_map(fx,fu,X,U)
fx[k] == A
fu[k] == B
fx_view,fu_view = calc_dyn_jacob_map!(fxu,X,U)
fx_view[k] == fx[k]
fu_view[k] == fu[k]
Fd(X[k],U[k]) == [A B]

@btime calc_dyn_jacob(fx,fu,res.X,res.U)
@btime calc_dyn_jacob_map(fx,fu,res.X,res.U)
@btime calc_dyn_jacob_map!(fxu,res.X,res.U)

###################################
#         PARALLEL STUFF          #
###################################
N = 1000

@everywhere Fd,Fc = TrajectoryOptimization.generate_dynamics_jacobians($model,$dt,TrajectoryOptimization.rk3,:zoh)
X = SharedArray{Float64,2}((n,N))
U = SharedArray{Float64,2}((m,N))
Fxu = SharedArray{Float64,3}((n,n+m,N))

X .= Xrand
U .= Urand

# Split the timesteps betwen processors
@everywhere function split_timesteps(N::Int,idx::Int=myid()-1)
    if idx == 0
        return 1:0
    end
    nchunks = length(workers())
    split = [round(Int,s) for s in range(0,stop=N,length=nchunks+1)]
    return split[idx]+1:split[idx+1]
end

@everywhere function jacobian_chunk!(Fxu::SharedArray,X::SharedArray,U::SharedArray,inds::UnitRange)
    n,N = size(X)
    m = size(U,1)
    for k in inds
        Fxu[:,:,k] = Fd(X[:,k],U[:,k])
    end
end

@everywhere function calc_jacobian_chunks(Fxu::SharedArray,X::SharedArray,U::SharedArray)
    N = size(X,2)
    @sync begin
        for w in workers()
            @async remotecall_wait(jacobian_chunk!,w,Fxu,X,U,split_timesteps(N,w-1))
        end
    end
end
inds = split_timesteps(N,1)
jacobian_chunk!(Fxu,X,U,inds)
@time calc_jacobian_chunks(Fxu,X,U)
Fxu[1:n,1:n,k] ≈ fx[k]

TrajectoryOptimization.to_array(fx) ≈ Fxu[1:n,1:n,1:end-1]

@btime calc_jacobian_chunks(Fxu,X,U)

function calculate_jacobians_parallel!(res::ConstrainedIterResults, solver::Solver)::Nothing #TODO change to inplace '!' notation throughout the code
    N = solver.N
    for k = 1:N-1
        if solver.control_integration == :foh
            res.fx[k], res.fu[k], res.fv[k] = solver.Fd(res.X[k], res.U[k], res.U[k+1])
            res.Ac[k], res.Bc[k] = solver.Fc(res.X[k], res.U[k])
        else
            res.fx[k], res.fu[k] = solver.Fd(res.X[k], res.U[k])
        end
        solver.c_jacobian(res.Cx[k], res.Cu[k], res.X[k],res.U[k])
    end

    if solver.control_integration == :foh
        res.Ac[N], res.Bc[N] = solver.Fc(res.X[N], res.U[N])
        solver.c_jacobian(res.Cx[N], res.Cu[N], res.X[N],res.U[N])
    end

    solver.c_jacobian(res.Cx_N, res.X[N])
    return nothing
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of alignments.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Column alignments" begin
    # Left
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""
    result = pretty_table(String, data;
                          alignment = :l,
                          backend = :html,
                          standalone = false)

    @test result == expected

    # Center
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          alignment = :c,
                          backend = :html,
                          standalone = false)

    @test result == expected

    # Per column configuration
    # ==========================================================================

    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""
    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          backend = :html,
                          standalone = false)

    @test result == expected
end

@testset "Cell alignment override" begin
    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""
    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          backend = :html,
                          cell_alignment = Dict((3,1) => :r,
                                                (3,2) => :l,
                                                (1,4) => :l,
                                                (3,4) => :c,
                                                (4,4) => :c,
                                                (6,4) => :l ),
                          standalone = false)

    @test result == expected
end

@testset "Header alignment" begin
    header = ["A" "B" "C" "D"
              "a" "b" "c" "d"]

    expected = """

    A
    B
    C
    D

    a
    b
    c
    d

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""
    result = pretty_table(String, data, header;
                          alignment = [:l,:r,:c,:r],
                          backend = :html,
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ),
                          header_alignment = [:l,:c,:r,:r],
                          standalone = false)

    @test result == expected

    expected = """

    A
    B
    C
    D

    a
    b
    c
    d

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data, header;
                          alignment = [:l,:r,:c,:r],
                          backend = :html,
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ),
                          header_alignment = [:l,:c,:r,:r],
                          header_cell_alignment = Dict( (1,1) => :r,
                                                        (2,2) => :l,
                                                        (2,4) => :c),
                          standalone = false)

    @test result == expected
end


text:
using ModelingToolkit, LinearAlgebra, SparseArrays

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const _DD = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 8
const X = reshape([i for i in 1:N for j in 1:N],N,N)
const Y = reshape([j for i in 1:N for j in 1:N],N,N)
const α₁ = 1.0.*(X.>=4*N/5)

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the initial condition as normal arrays
@variables du[1:N,1:N,1:3] u[1:N,1:N,1:3] MyA[1:N,1:N] AMx[1:N,1:N] DA[1:N,1:N]

# Define the discretized PDE as an ODE function
function f(du,u,p,t)
   A = @view  u[:,:,1]
   B = @view  u[:,:,2]
   C = @view  u[:,:,3]
  dA = @view du[:,:,1]
  dB = @view du[:,:,2]
  dC = @view du[:,:,3]
  mul!(MyA,My,A)
  mul!(AMx,A,Mx)
  @. DA = _DD*(MyA + AMx)
  @. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
  @. dB = α₂ - β₂*B - r₁*A*B + r₂*C
  @. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

f(du,u,nothing,0.0)

multithreadedf = eval(ModelingToolkit.build_function(du,u,parallel=ModelingToolkit.MultithreadedForm())[2])

MyA = zeros(N,N);
AMx = zeros(N,N);
DA = zeros(N,N);
# Loop to catch syncronization issues
for i in 1:100
   _du = rand(N,N,3)
   _u = rand(N,N,3)
   multithreadedf(_du,_u)
   _du2 = copy(_du)
   f(_du2,_u,nothing,0.0)
   @test _du ≈ _du2
end

#=
jac = sparse(ModelingToolkit.jacobian(vec(du),vec(u)))
fjac = eval(ModelingToolkit.build_function(jac,u,parallel=ModelingToolkit.SerialForm())[2])
multithreadedfjac = eval(ModelingToolkit.build_function(jac,u,parallel=ModelingToolkit.MultithreadedForm())[2])

u = rand(N,N,3)
J = similar(jac,Float64)
fjac(J,u)

J2 = similar(jac,Float64)
multithreadedfjac(J2,u)
@test J ≈ J2

using FiniteDiff
J3 = Array(similar(jac,Float64))
FiniteDiff.finite_difference_jacobian!(J2,(du,u)->f!(du,u,nothing,nothing),u)
maximum(J2 .- Array(J)) < 1e-5
=#

using Distributed
addprocs(4)
distributedf = eval(ModelingToolkit.build_function(du,u,parallel=ModelingToolkit.DistributedForm())[2])

using Dagger
daggerf = eval(ModelingToolkit.build_function(du,u,parallel=ModelingToolkit.DaggerForm())[2])

jac = sparse(ModelingToolkit.jacobian(vec(du),vec(u),simplify=false))
serialjac = eval(ModelingToolkit.build_function(vec(jac),u)[2])
multithreadedjac = eval(ModelingToolkit.build_function(vec(jac),u,parallel=ModelingToolkit.MultithreadedForm())[2])
distributedjac = eval(ModelingToolkit.build_function(vec(jac),u,parallel=ModelingToolkit.DistributedForm())[2])
daggerjac = eval(ModelingToolkit.build_function(vec(jac),u,parallel=ModelingToolkit.DaggerForm())[2])

MyA = zeros(N,N)
AMx = zeros(N,N)
DA  = zeros(N,N)
_du = rand(N,N,3)
_u = rand(N,N,3)

f(_du,_u,nothing,0.0)
multithreadedf(_du,_u)
#distributedf(_du,_u)
#daggerf(_du,_u)

#=
using BenchmarkTools
@btime f(_du,_u,nothing,0.0)
@btime multithreadedf(_du,_u)
@btime distributedf(_du,_u)
@btime daggerf(_du,_u)

_jac = similar(jac,Float64)
@btime serialjac(_jac,_u)
@btime multithreadedjac(_jac,_u)
@btime distributedjac(_jac,_u)
@btime daggerjac(_jac,_u)
=#


text:
"""
Wraps the mutation operator, so that it could
be used as crossover operator.
"""
struct MutationWrapper{OP<:MutationOperator} <: CrossoverOperator{1,1}
    inner::OP

    MutationWrapper(mutation::OP) where {OP<:MutationOperator} =
        new{OP}(mutation)
end

function apply!(wrapper::MutationWrapper,
                target::Individual, targetIndex::Int, pop, parentIndices)
    @assert length(parentIndices) == 1
    apply!(wrapper.inner, copyto!(target, viewer(pop, parentIndices[1])), targetIndex)
end

Base.show(io::IO, xover::MutationWrapper) = print(io, "MutationWrapper(", xover.inner, ")")

# convert mutation to crossover by wrapping it in `MutationWrapper`
Base.convert(::Type{CrossoverOperator}, op::MutationOperator) = MutationWrapper(op)


text:
export shcat, shcat!

@with_kw struct ShcatOptions
    show_ends = false
    number = false
    number_nonblank = false
end

function parse_options(::Type{ShcatOptions}, options)
    ShcatOptions(
        show_ends = get_option(options, :E, :show_ends, default = false).second,
        number = get_option(options, :n, :number, default = false).second,
        number_nonblank = get_option(
            options,
            :b,
            :number_nonblank,
            default = false,
        ).second,
    )
end

"""
    shcat(options..., files...)

Return lines of files concatenated into a single array.

With no files, or when file is "-", read standard input.

    :b, :number_nonblank    number nonempty output lines as `n => line`,
                            overrides :n
    :E, :show_ends          display \$ at end of each line
    :n, :number             number all output lines as `n => line`
      :help     display this help and exit

    FOLLOWING NOT IMPLEMENTED
    :A, :show_all           equivalent to :v, :E, :T
    :e                      equivalent to :v, :E
    :s, :squeeze_blank      suppress repeated empty output lines
    :t                      equivalent to :v, :T
    :T, :show_tabs          display TAB characters as ^I
    :u                      (ignored)
    :v, :show_nonprinting   use ^ and M- notation, except for LFD and TAB
      :version  output version information and exit

# Examples

    shcat("f", "-", "g") # Output f's contents, then standard input, then g's
                         # contents.
    shcat!()        # Copy standard input to standard output.

"""
function shcat(args...)
    options_seq, pos_args = parse_args(args)
    do_print_usage(options_seq, @doc shcat) && return
    options = parse_options(ShcatOptions, options_seq)
    _shcat(options, pos_args)
end

function _shcat(options, args)
    number_all = options.number

    lines = []
    if length(args) == 0
        args = ["-"]
    end
    for fn in args
        fn == "-" && (fn = stdin)

        file_lines = readlines(fn)
        if options.number_nonblank
            n = 1
            numbered_lines = Array{Any}(undef, length(file_lines))
            for i in eachindex(file_lines)
                if length(file_lines[i]) > 0
                    numbered_lines[i] = n => file_lines[i]
                    n += 1
                else
                    numbered_lines[i] = ""
                end
            end
            file_lines = numbered_lines
        elseif number_all
            file_lines = [i => line for (i, line) in enumerate(file_lines)]
        end
        append!(lines, file_lines)
    end

    lines
end

function shcat!(args...)
    options_seq, pos_args = parse_args(args)
    do_print_usage(options_seq, @doc shcat) && return
    options = parse_options(ShcatOptions, options_seq)
    lines = _shcat(options, pos_args)
    nlines = length(lines)
    nspaces = (nlines |> log10 |> ceil |> Int) - 1
    n = 0
    ndigits = 1
    for line in lines
        if line isa Pair
            n += 1
            if n % 10^ndigits == 0
                nspaces -= 1
                ndigits += 1
            end
            print(" "^nspaces, line.first, " ")
            line = line.second
        end
        print(line)
        options.show_ends && print("\$")
        println()
    end
end


text:
Base.@kwdef struct Word
    text::String
    style::Crayon = Crayon()
end

Base.@kwdef struct Paragraph
    block::Block
    words::Vector{Word}
    scroll::Int = 1
    number_of_lines::Ref{Int}
    function Paragraph(block, words, scroll)
        if scroll < 1
            scroll = 1
        end
        new(block, words, scroll, Ref{Int}(0))
    end
end

make_words(sentence::String, style = Crayon()) = [Word(word, style) for word in split(sentence)]

function draw(p::Paragraph, area::Rect, buf::Buffer)
    draw(p.block, area, buf);
    text_area = inner(p.block, area)

    height(text_area) < 1 && return

    lines = split(wrap(join([word.text for word in p.words], " "), width = Int(width(text_area))), "\n")
    p.number_of_lines[] = length(lines)
    word_count = 1
    for (i, line) in enumerate(lines)
        words = split(line)
        if i < p.scroll || i - p.scroll > height(text_area)
            for _ in words
                word_count += 1
            end
            continue
        end
        count = left(text_area)
        for (j, str) in enumerate(words)
            set(buf, count, top(text_area) + i - p.scroll, String(str), p.words[word_count].style)
            count += length(str)
            word_count += 1
            if j != length(words)
                set(buf, count, top(text_area) + i - p.scroll, " ", p.words[word_count].style)
                count += 1
            end
        end
    end

end


text:
using AtomEnergyLevels, Libxc, Printf

function libxcfun(xc::Symbol)
    fxc = Functional(xc)

    if getproperty(fxc, :kind) ≠ :exchange_correlation
        throw(DomainError(getproperty(fxc, :identifier), "is not an exchange-correlation functional"))
    end

    function xc_lda(ρ, vxc = similar(ρ), exc = similar(ρ))
        if !(length(ρ) == length(vxc) == length(exc))
            throw(DimensionMismatch("Dimensions of the arrays ρ, vxc, and exc doesn't match"))
        end
        evaluate!(fxc, rho = ρ, vrho = vxc, zk = exc)
        return vxc, exc
    end

    if getproperty(fxc, :family) == :lda
        return XC_Functional(xc_type_lda, xc_lda)
    else
        return nothing
    end
end

function libxcfun(x::Symbol, c::Symbol)
    fx, fc = Functional(x), Functional(c)

    if getproperty(fx, :kind) ≠ :exchange || getproperty(fc, :kind) ≠ :correlation
        throw(DomainError(getproperty(fxc, :identifier), "is not an exchange-correlation functional"))
    end

    function xc_lda(ρ, vxc = similar(ρ), exc = similar(ρ))
        if !(length(ρ) == length(vxc) == length(exc))
            throw(DimensionMismatch("Dimensions of the arrays ρ, vxc, and exc doesn't match"))
        end
        exch = evaluate(fx, rho = ρ)
        corr = evaluate(fc, rho = ρ)
        @inbounds for i in eachindex(ρ)
            vxc[i] = exch.vrho[i] + corr.vrho[i]
            exc[i] = exch.zk[i]   + corr.zk[i]
        end
        return vxc, exc
    end

    if getproperty(fx, :family) == :lda && getproperty(fc, :family) == :lda
        return XC_Functional(xc_type_lda, xc_lda)
    else
        return nothing
    end
end

lda(2, xc = libxcfun(:lda_x, :lda_c_vwn))

function zoo(z)
    function is_lda_xc(x::Functional)
        id = getproperty(x, :identifier)
        occursin(r"_\dd_", string(id)) && return false

        getproperty(x, :family) == :lda &&
        getproperty(x, :kind)   == :exchange_correlation &&

        intersect([:vxc, :exc], getproperty(x, :flags)) == [:vxc, :exc]
    end

    function is_lda_c(x::Functional)
        id = getproperty(x, :identifier)
        occursin(r"_\dd_", string(id)) && return false

        getproperty(x, :family) == :lda &&
        getproperty(x, :kind)   == :correlation &&

        intersect([:vxc, :exc], getproperty(x, :flags)) == [:vxc, :exc]
    end

    xc_funcs = [] # filter(is_lda_xc, Functional.(available_functionals()))
     c_funcs = filter(is_lda_c, Functional.(available_functionals()))

    ids_xc = map(λ -> (getproperty(λ, :identifier),), xc_funcs)
    ids_c  = map(λ -> (:lda_x, getproperty(λ, :identifier)), c_funcs)

    energies = []

    for xc in vcat(ids_xc, ids_c)
        @printf("Testing %s\n", string(xc))
        try
            push!(energies, lda(z, xc = libxcfun(xc...)).energy.total)
        catch
            push!(energies, NaN)
            continue
        end
    end

    @printf("== Total energies of the %s atom in various LDA functionals ==\n", keys(atom)[z])
    for (xc, en) in zip(vcat(xc_funcs, c_funcs), energies)
        if getproperty(xc, :kind) == :correlation
            abbr = string.("lda_x + ", getproperty(xc, :identifier))
        else
            abbr = string.("        ", getproperty(xc, :identifier))
        end
        name = getproperty(xc, :name)
        year = last(split(getproperty(xc, :references)[1].reference))
        @printf("%-30s\t%12.6f\t%s %s\n", abbr, en, name, year)
    end
end

zoo(2)
#=
== Total energies of the He atom in various LDA functionals ==
        lda_xc_teter93             -2.833976    Teter 93 (1996)
        lda_xc_zlp                 -2.913789    Zhao, Levy & Parr, Eq. (20) (1993)
        lda_xc_ksdt                -2.833100    Karasiev, Sjostrom, Dufty & Trickey (2014)
        lda_xc_lp_a                -2.869133    Lee-Parr reparametrization A (1990)
        lda_xc_lp_b                -2.903489    Lee-Parr reparametrization B (1990)
        lda_xc_gdsmfb              -2.832794    Groth, Dornheim, Sjostrom, Malone, Foulkes, Bonitz (2017)
        lda_xc_bn05                -2.244445    Baer and Neuhauser, gamma=1 (2005)
lda_x + lda_c_wigner               -2.818849    Wigner (1938)
lda_x + lda_c_rpa                  -2.837328    Random Phase Approximation (RPA) (1957)
lda_x + lda_c_hl                   -2.839923    Hedin & Lundqvist (1971)
lda_x + lda_c_gl                   -2.860137    Gunnarson & Lundqvist (1976)
lda_x + lda_c_xalpha               -3.170112    Slater's Xalpha (1951)
lda_x + lda_c_vwn                  -2.834836    Vosko, Wilk & Nusair (VWN5) (1980)
lda_x + lda_c_vwn_rpa              -2.872169    Vosko, Wilk & Nusair (VWN5_RPA) (1980)
lda_x + lda_c_pz                   -2.834291    Perdew & Zunger (1981)
lda_x + lda_c_pz_mod               -2.834313    Perdew & Zunger (Modified) parts
lda_x + lda_c_ob_pz                -2.830831    Ortiz & Ballone (PZ parametrization) (1994)
lda_x + lda_c_pw                   -2.834455    Perdew & Wang (1992)
lda_x + lda_c_pw_mod               -2.834455    Perdew & Wang (modified) (http://dft.rutgers.edu/pubs/PBE.asc)
lda_x + lda_c_ob_pw                -2.831420    Ortiz & Ballone (PW parametrization) (1994)
lda_x + lda_c_vbh                  -2.870439    von Barth & Hedin (1972)
lda_x + lda_c_ml1                  -2.779743    Modified LSD (version 1) of Proynov and Salahub (1994)
lda_x + lda_c_ml2                  -2.751033    Modified LSD (version 2) of Proynov and Salahub (1994)
lda_x + lda_c_gombas               -2.837973    Gombas (1965)
lda_x + lda_c_pw_rpa               -2.871024    Perdew & Wang (fit to the RPA energy) (1992)
lda_x + lda_c_rc04                 -2.816964    Ragot-Cortona (2004)
lda_x + lda_c_vwn_1                -2.834836    Vosko, Wilk & Nusair (VWN1) (1980)
lda_x + lda_c_vwn_2                -2.834836    Vosko, Wilk & Nusair (VWN2) (1980)
lda_x + lda_c_vwn_3                -2.834836    Vosko, Wilk & Nusair (VWN3) (1980)
lda_x + lda_c_vwn_4                -2.834836    Vosko, Wilk & Nusair (VWN4) (1980)
lda_x + lda_c_chachiyo             -2.831427    Chachiyo simple 2 parameter correlation (2016)
lda_x + lda_c_lp96                       NaN    Liu-Parr correlation (1996)
lda_x + lda_c_chachiyo_mod         -2.831427    Chachiyo simple 2 parameter correlation with modified spin scaling [cond-mat.mtrl-sci]
lda_x + lda_c_karasiev_mod         -2.833235    Karasiev reparameterization of Chachiyo [cond-mat.mtrl-sci]
lda_x + lda_c_mcweeny              -2.850426    McWeeny 76 3--31
lda_x + lda_c_br78                 -2.764565    Brual & Rothstein 78 (1978)
lda_x + lda_c_pk09                 -2.834520    Proynov and Kong 2009 (2009)
lda_x + lda_c_ow_lyp               -2.780839    Wigner with corresponding LYP parameters (1995)
lda_x + lda_c_ow                   -2.784822    Optimized Wigner (1995)
lda_x + lda_c_gk72                 -2.839369    Gordon and Kim 1972 (1972)
lda_x + lda_c_karasiev             -2.833235    Karasiev reparameterization of Chachiyo (2016)
lda_x + lda_c_upw92                -2.834749    Ruggeri, Rios, and Alavi unrestricted fit (2018)
lda_x + lda_c_rpw92                -2.834455    Ruggeri, Rios, and Alavi restricted fit (2018)
=#


text:
# This file holds different types of investment functions
@inline function quadratic_investment(p::S, χ::S, δ::S) where {S <: Real}
    return (p .- 1.) .^ 2 ./ (2. .* χ) + (p .- 1.) ./ χ .+ δ
end

@inline function growth_quadratic_investment(p::S, χ::S, δ::S) where {S <: Real}
    return (p - 1.) / χ  + δ
end

@inline function derivative_quadratic_investment(p::S, χ::S) where {S <: Real}
    return p ./ χ + 1. ./ χ
end

@inline function derivative_quadratic_investment_li2020(p::S, χ::S) where {S <: Real}
    return p ./ χ
end


text:
unsafe_load(pointer.(strs) + 3)

a = "def"
b  = unsafe_string(pointer(a))

pointer(a)
pointer(b)

a = UInt8[64, 65, C_NULL, 62, 53]
pp= pointer(a)
print(pp)
@time as = String(a)
pointer(as)

unsafe_wrap(a, Ptr{String})

a = "def"
unsafe_load(pointer(a)+3)

x = "id".*string.(rand(UInt8, 1_000_000))

io = open("c:/data/io.bin", "w")
ncw = write.(Ref(io), x)
close(io)

io = open("c:/data/io.bin", "r")
buffer = Vector{UInt8}(undef, sum(ncw))
readbytes!(io, buffer, sum(ncw))
close(io)

@time aa = test(buffer)

function test(array)
    start = 1
    strings = String[]
    GC.@preserve array begin
        ptr = pointer(array) - 1
        for i in eachindex(array)
            @inbounds char = array[i]
            if char == UInt8(',')
                len = i - start
                str = unsafe_string(ptr + start, len)
                push!(strings, str)
                start = i + 1
            end
        end
    end
    strings
end


text:
@testset "product" begin

    # Test the multiplication of a GP by a constant.
    let
        rng, N, N′, D = MersenneTwister(123456), 10, 11, 2
        X, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        g1, c, c′ = GP(ConstantMean(1.0), EQ(), GPC()), -4.3, 2.1
        g2, g2′ = c * g1, g1 * c′
        g3, g3′ = c * g2, g2′ * c′
        g4, g4′ = c * g3, g3′ * c′

        @test mean_vec(g2(X)) == c .* mean_vec(g1(X))
        @test mean_vec(g3(X)) == c .* mean_vec(g2(X))
        @test mean_vec(g4(X)) == c .* mean_vec(g3(X))
        @test mean_vec(g2′(X)) == mean_vec(g1(X)) .* c′
        @test mean_vec(g3′(X)) == mean_vec(g2′(X)) .* c′
        @test mean_vec(g4′(X)) == mean_vec(g3′(X)) .* c′

        @test cov(g2(X)) ≈ c^2 .* cov(g1(X))
        @test cov(g3(X)) ≈ c^2 .* cov(g2(X))
        @test cov(g4(X)) ≈ c^2 .* cov(g3(X))
        @test cov(g2′(X)) ≈ c′^2 .* cov(g1(X))
        @test cov(g3′(X)) ≈ c′^2 .* cov(g2′(X))
        @test cov(g4′(X)) ≈ c′^2 .* cov(g3′(X))

        @test xcov(g2(X), g1(X′)) ≈ c .* xcov(g1(X), g1(X′))
        @test xcov(g3(X), g1(X′)) ≈ c^2 .* xcov(g1(X), g1(X′))
        @test xcov(g4(X), g1(X′)) ≈ c^3 .* xcov(g1(X), g1(X′))
        @test xcov(g2′(X), g1(X′)) ≈ c′ .* xcov(g1(X), g1(X′))
        @test xcov(g3′(X), g1(X′)) ≈ c′^2 .* xcov(g1(X), g1(X′))
        @test xcov(g4′(X), g1(X′)) ≈ c′^3 .* xcov(g1(X), g1(X′))

        @test xcov(g2(X), g2′(X′)) ≈ (c * c′) .* xcov(g1(X), g1(X′))
        @test xcov(g3(X), g3′(X′)) ≈ (c * c′)^2 .* xcov(g1(X), g1(X′))
        @test xcov(g4(X), g4′(X′)) ≈ (c * c′)^3 .* xcov(g1(X), g1(X′))

        @test xcov(g2(X), g3′(X′)) ≈ (c * c′^2) .* xcov(g1(X), g1(X′))
        @test xcov(g3(X), g2′(X′)) ≈ (c^2 * c′) .* xcov(g1(X), g1(X′))
        @test xcov(g2(X), g4′(X′)) ≈ (c * c′^3) .* xcov(g1(X), g1(X′))
        @test xcov(g4(X), g2′(X′)) ≈ (c^3 * c′) .* xcov(g1(X), g1(X′))
    end

    # Test the multiplication of a GP by a known function.
    let
        rng, N, N′, D = MersenneTwister(123456), 10, 11, 2
        X, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        g1, f, f′ = GP(ConstantMean(1.0), EQ(), GPC()), x->sum(sin, x), x->sum(cos, x)
        g2, g2′ = f * g1, g1 * f′
        g3, g3′ = f * g2, g2′ * f′
        g4, g4′ = f * g3, g3′ * f′

        @test mean_vec(g2(X)) == map(f, X) .* mean_vec(g1(X))
        @test mean_vec(g3(X)) == map(f, X) .* mean_vec(g2(X))
        @test mean_vec(g4(X)) == map(f, X) .* mean_vec(g3(X))
        @test mean_vec(g2′(X)) == mean_vec(g1(X)) .* map(f′, X)
        @test mean_vec(g3′(X)) == mean_vec(g2′(X)) .* map(f′, X)
        @test mean_vec(g4′(X)) == mean_vec(g3′(X)) .* map(f′, X)

        fX, f′X = map(f, X), map(f′, X)
        @test cov(g2(X)) ≈ fX .* cov(g1(X)) .* fX'
        @test cov(g3(X)) ≈ fX .* cov(g2(X)) .* fX'
        @test cov(g4(X)) ≈ fX .* cov(g3(X)) .* fX'
        @test cov(g2′(X)) ≈ f′X .* cov(g1(X)) .* f′X'
        @test cov(g3′(X)) ≈ f′X .* cov(g2′(X)) .* f′X'
        @test cov(g4′(X)) ≈ f′X .* cov(g3′(X)) .* f′X'

        fX′, f′X′ = map(f, X′), map(f′, X′)
        @test xcov(g2(X), g1(X′)) ≈ fX .* xcov(g1(X), g1(X′))
        @test xcov(g3(X), g1(X′)) ≈ fX.^2 .* xcov(g1(X), g1(X′))
        @test xcov(g4(X), g1(X′)) ≈ fX.^3 .* xcov(g1(X), g1(X′))
        @test xcov(g2′(X), g1(X′)) ≈ f′X .* xcov(g1(X), g1(X′))
        @test xcov(g3′(X), g1(X′)) ≈ (f′X.^2) .* xcov(g1(X), g1(X′))
        @test xcov(g4′(X), g1(X′)) ≈ (f′X.^3) .* xcov(g1(X), g1(X′))

        @test xcov(g1(X′), g2(X)) == xcov(g2(X), g1(X′))'
        @test xcov(g1(X′), g3(X)) == xcov(g3(X), g1(X′))'
        @test xcov(g1(X′), g4(X)) == xcov(g4(X), g1(X′))'

        @test xcov(g2(X), g2′(X′)) ≈ fX .* xcov(g1(X), g1(X′)) .* f′X′'
        @test xcov(g3(X), g3′(X′)) ≈ fX.^2 .* xcov(g1(X), g1(X′)) .* (f′X′.^2)'
        @test xcov(g4(X), g4′(X′)) ≈ fX.^3 .* xcov(g1(X), g1(X′)) .* (f′X′.^3)'

        @test xcov(g2(X), g3′(X′)) ≈ fX .* xcov(g1(X), g1(X′)) .* (f′X′').^2
        @test xcov(g3(X), g2′(X′)) ≈ fX.^2 .* xcov(g1(X), g1(X′)) .* (f′X′')
        @test xcov(g2(X), g4′(X′)) ≈ fX .* xcov(g1(X), g1(X′)) .* (f′X′').^3
        @test xcov(g4(X), g2′(X′)) ≈ fX.^3 .* xcov(g1(X), g1(X′)) .* (f′X′')
    end
end


text:
function anglevariables(m::JuMP.AbstractModel, system_formulation::Type{S}, sys::PowerSystems.PowerSystem) where {S <: Union{PM.AbstractDCPForm, PM.AbstractACPForm}}

    on_set = [d.name for d in sys.buses if d.available == true]

    time_range = 1:sys.time_periods

    theta = @variable(m, theta[on_set,time_range])

end

function voltagevariables(m::JuMP.AbstractModel, system_formulation::Type{S}, sys::PowerSystems.PowerSystem) where {S <: PM.AbstractACPForm}

    on_set = [d.name for d in sys.buses if d.available == true]

    time_range = 1:sys.time_periods

    vm = @variable(m, vm[on_set,time_range])
end

function voltagevariables(m::JuMP.AbstractModel, system_formulation::Type{S}, sys::PowerSystems.PowerSystem) where {S <: PM.AbstractACRForm}

    on_set = [d.name for d in sys.buses if d.available == true]

    time_range = 1:sys.time_periods

    vm_re = @variable(m, vm_re[on_set,time_range])

    vm_im = @variable(m, vm_re[on_set,time_range])

end


text:
module TestSubDataFrame
    using Test, DataFrames

    @testset "copy - SubDataFrame" begin
        df = DataFrame(x = 1:10, y = 1.0:10.0)
        sdf = view(df, 1:2, 1:1)
        @test sdf isa SubDataFrame
        @test copy(sdf) isa DataFrame
        @test sdf == copy(sdf)
    end

    @testset "view -- DataFrame" begin
        df = DataFrame(x = 1:10, y = 1.0:10.0)
        @test view(df, 1, :) == DataFrameRow(df, 1)
        @test view(df, UInt(1), :) == DataFrameRow(df, 1)
        @test view(df, BigInt(1), :) == DataFrameRow(df, 1)
        @test view(df, UInt(1):UInt(1), :) == SubDataFrame(df, 1:1)
        @test view(df, BigInt(1):BigInt(1), :) == SubDataFrame(df, 1:1)
        @test view(df, 1:2, :) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), :) == first(df, 2)
        @test view(df, [1, 2], :) == first(df, 2)
        @test view(df, 1, :x) == view(df[:x], 1)
        @test view(df, 1, :x) isa SubArray
        @test size(view(df, 1, :x)) == ()
        @test view(df, 1:2, :x) == df[:x][1:2]
        @test view(df, 1:2, :x) isa SubArray
        @test view(df, vcat(trues(2), falses(8)), :x) == view(df[:x], vcat(trues(2), falses(8)))
        @test view(df, [1, 2], :x) == view(df[:x], [1, 2])
        @test view(df, 1, 1) == view(df[1], 1)
        @test view(df, 1, 1) isa SubArray
        @test size(view(df, 1, 1)) == ()
        @test view(df, 1:2, 1) == df[1][1:2]
        @test view(df, 1:2, 1) isa SubArray
        @test view(df, vcat(trues(2), falses(8)), 1) == view(df[1], vcat(trues(2), falses(8)))
        @test view(df, [1, 2], 1) == view(df[1], [1,2])
        @test view(df, 1:2, 1) == df[1][1:2]
        @test view(df, 1:2, 1) isa SubArray
        @test view(df, 1, [:x, :y]) == DataFrameRow(df[[:x, :y]], 1)
        @test view(df, 1:2, [:x, :y]) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [:x, :y]) == first(df, 2)
        @test view(df, [1, 2], [:x, :y]) == first(df, 2)
        @test view(df, 1, [1, 2]) == DataFrameRow(df[1:2], 1)
        @test view(df, 1:2, [1, 2]) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [1, 2]) == first(df, 2)
        @test view(df, [1, 2], [1, 2]) == first(df, 2)
        @test view(df, 1, trues(2)) == DataFrameRow(df[trues(2)], 1)
        @test view(df, 1:2, trues(2)) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), trues(2)) == first(df, 2)
        @test view(df, [1, 2], trues(2)) == first(df, 2)
        @test view(df, Integer[1, 2], :) == first(df, 2)
        @test view(df, UInt[1, 2], :) == first(df, 2)
        @test view(df, BigInt[1, 2], :) == first(df, 2)
        @test view(df, Union{Int, Missing}[1, 2], :) == first(df, 2)
        @test view(df, Union{Integer, Missing}[1, 2], :) == first(df, 2)
        @test view(df, Union{UInt, Missing}[1, 2], :) == first(df, 2)
        @test view(df, Union{BigInt, Missing}[1, 2], :) == first(df, 2)
        @test view(df, :) == df
        @test view(df, :, :) == df
        @test view(df, 1, :) == DataFrameRow(df, 1)
        @test view(df, :, 1) == df[:, 1]
        @test view(df, :, 1) isa SubArray
        @test_throws ArgumentError view(df, [missing, 1])
        @test_throws ArgumentError view(df, [missing, 1], :)
    end

    @testset "view -- SubDataFrame" begin
        df = view(DataFrame(x = 1:10, y = 1.0:10.0), 1:10, :)
        @test view(df, 1, :) == DataFrameRow(df, 1)
        @test view(df, UInt(1), :) == DataFrameRow(df, 1)
        @test view(df, BigInt(1), :) == DataFrameRow(df, 1)
        @test view(df, 1:2, :) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), :) == first(df, 2)
        @test view(df, [1, 2], :) == first(df, 2)
        @test view(df, 1, :x) == view(df[:x], 1)
        @test view(df, 1, 1) isa SubArray
        @test size(view(df, 1, 1)) == ()
        @test view(df, 1:2, :x) == view(df[:x], 1:2)
        @test view(df, vcat(trues(2), falses(8)), :x) == view(df[:x], vcat(trues(2), falses(8)))
        @test view(df, [1, 2], :x) == view(df[:x], [1, 2])
        @test view(df, 1, 1) == view(df[1], 1)
        @test view(df, 1, 1) isa SubArray
        @test size(view(df, 1, 1)) == ()
        @test view(df, 1:2, 1) == view(df[:x], 1:2)
        @test view(df, vcat(trues(2), falses(8)), 1) == view(df[:x], vcat(trues(2), falses(8)))
        @test view(df, [1, 2], 1) == view(df[:x], [1, 2])
        @test view(df, 1, [:x, :y]) == DataFrameRow(df[[:x,:y]], 1)
        @test view(df, 1:2, [:x, :y]) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [:x, :y]) == first(df, 2)
        @test view(df, [1, 2], [:x, :y]) == first(df, 2)
        @test view(df, 1, [1, 2]) == DataFrameRow(df[1:2], 1)
        @test view(df, 1:2, [1, 2]) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), [1, 2]) == first(df, 2)
        @test view(df, [1, 2], [1, 2]) == first(df, 2)
        @test view(df, 1, trues(2)) == DataFrameRow(df[trues(2)], 1)
        @test view(df, 1:2, trues(2)) == first(df, 2)
        @test view(df, vcat(trues(2), falses(8)), trues(2)) == first(df, 2)
        @test view(df, [1, 2], trues(2)) == first(df, 2)
        @test view(df, Integer[1, 2], :) == first(df, 2)
        @test view(df, UInt[1, 2], :) == first(df, 2)
        @test view(df, BigInt[1, 2], :) == first(df, 2)
        @test view(df, Union{Int, Missing}[1, 2], :) == first(df, 2)
        @test view(df, Union{Integer, Missing}[1, 2], :) == first(df, 2)
        @test view(df, Union{UInt, Missing}[1, 2], :) == first(df, 2)
        @test view(df, Union{BigInt, Missing}[1, 2], :) == first(df, 2)
        @test view(df, :) == df
        @test view(df, :, :) == df
        @test view(df, 1, :) == DataFrameRow(df, 1)
        @test view(df, :, 1) == df[:, 1]
        @test view(df, :, 1) isa SubArray
        @test_throws ArgumentError view(df, [missing, 1])
        @test_throws ArgumentError view(df, [missing, 1], :)
    end

    @testset "getproperty, setproperty! and propertynames" begin
        x = collect(1:10)
        y = collect(1.0:10.0)
        df = view(DataFrame(x = x, y = y), 2:6, :)

        @test Base.propertynames(df) == names(df)

        @test df.x == 2:6
        @test df.y == 2:6
        @test_throws KeyError df.z

        df.x = 1:5
        @test df.x == 1:5
        @test x == [1; 1:5; 7:10]
        df.y = 1
        @test df.y == [1, 1, 1, 1, 1]
        @test y == [1; 1; 1; 1; 1; 1; 7:10]
        @test_throws ErrorException df.z = 1:5
        @test_throws ErrorException df.z = 1
    end

    @testset "dump" begin
        y = 1.0:10.0
        df = view(DataFrame(y=y), 2:6, :)
        @test sprint(dump, df) == """
                                  SubDataFrame{UnitRange{$Int}}  5 observations of 1 variables
                                    y: [2.0, 3.0, 4.0, 5.0, 6.0]\n
                                  """
    end

    @testset "deleterows!" begin
        y = 1.0:10.0
        df = view(DataFrame(y=y), 2:6, :)
        @test_throws ArgumentError deleterows!(df, 1)
    end

    @testset "parent" begin
        df = DataFrame(a=Union{Int, Missing}[1, 2, 3, 1, 2, 2],
                       b=[2.0, missing, 1.2, 2.0, missing, missing],
                       c=["A", "B", "C", "A", "B", missing])
        @test parent(view(df, [4, 2], :)) === df
        @test parentindices(view(df, [4, 2], :)) == ([4,2], Base.OneTo(3))
        @test parent(view(df, [4, 2], 1:3)) !== df
        @test parentindices(view(df, [4, 2], 1:3)) == ([4, 2], Base.OneTo(3))
    end
end


text:
# # Housing data

# In this example, we create a linear regression model that predicts housing data.
# It replicates the housing data example from the [Knet.jl readme](https://github.com/denizyuret/Knet.jl).
# Although we could have reused more of Flux (see the MNIST example), the library's abstractions are very
# lightweight and don't force you into any particular strategy.

# A linear model can be created as a neural network with a single layer.
# The number of inputs is the same as the features that the data has.
# Each input is connected to a single output with no activation function.
# Then, the output of the model is a linear function that predicts unseen data.

# ![singleneuron](img/singleneuron.svg)

# Source: [Dive into Deep Learning](http://d2l.ai/chapter_linear-networks/linear-regression.html#from-linear-regression-to-deep-networks)

# To run this example, we need the following packages:

using Flux
using Flux: gradient
using Flux.Optimise: update!
using DelimitedFiles, Statistics
using Parameters: @with_kw

# We set default values for the learning rate (for the training routine) and the percentage of
# the data that we use when testing the model:

@with_kw mutable struct Hyperparams
    ## Learning rate
    lr::Float64 = 0.1
    ## Train Test split ratio, define percentage of data to be used as Test data
    split_ratio::Float64 = 0.1
end

# ## Data

# We create the function `get_processed_data` to load the housing data, normalize it,
# and finally split it into train and test datasets:

function get_processed_data(args)
    isfile("housing.data") ||
        download("https://raw.githubusercontent.com/MikeInnes/notebooks/master/housing.data",
            "housing.data")

    rawdata = readdlm("housing.data")'

    ## The last feature is our target -- the price of the house.
    split_ratio = args.split_ratio ## For the train test split

    x = rawdata[1:13,:]
    y = rawdata[14:14,:]

    ## Normalise the data
    x = (x .- mean(x, dims = 2)) ./ std(x, dims = 2)

    ## Split into train and test sets
    split_index = floor(Int,size(x,2)*split_ratio)
    x_train = x[:,1:split_index]
    y_train = y[:,1:split_index]
    x_test = x[:,split_index+1:size(x,2)]
    y_test = y[:,split_index+1:size(x,2)]

    train_data = (x_train, y_train)
    test_data = (x_test, y_test)

    return train_data,test_data
end

# This function performs the following tasks:

# 1. Downloads the housing data. The original size of the data is 505 rows and 14 columns.
# 2. Loads the data as a 14x505 matrix. This is the shape that Flux expects.
# 3. Splits the data into features and a target. Notice that the 14th row corresponds to the target for each example.
# 4. Normalizes the data. For more information on normalizing data, see [How to Use StandardScaler and MinMaxScaler Transforms in Python](https://machinelearningmastery.com/standardscaler-and-minmaxscaler-transforms-in-python/).
# 5. Splits the data into train and test datasets.

# ## Model
# We use a struct to define the model’s parameters.
# It contains an array for holding the weights *W* and a vector for the bias term *b*:

mutable struct model
    W::AbstractArray
    b::AbstractVector
end

# Also, we create the function `predict` to compute the model’s output:

predict(x, m) = m.W*x .+ m.b

# Notice that the function `predict` takes as an argument the model struct we defined above.

# ## Loss function

# The most commonly used loss function for Linear Regression is Mean Squared Error (MSE).
# We define the MSE function as:

meansquarederror(ŷ, y) = sum((ŷ .- y).^2)/size(y, 2)

# **Note:** An implementation of the MSE function is also available in
# [Flux](https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.mse).

# ## Train function
# Finally, we define the `train` function so that the model learns the best parameters (*W* and *b*):

function train(; kws...)
    ## Initialize the Hyperparamters
    args = Hyperparams(; kws...)

    ## Load the data
    (x_train,y_train),(x_test,y_test) = get_processed_data(args)

    ## The model
    m = model((randn(1,13)),[0.])

    loss(x, y) = meansquarederror(predict(x, m), y)

    ## Training
    η = args.lr
    θ = params(m.W, m.b)

    for i = 1:500
        g = gradient(() -> loss(x_train, y_train), θ)
        for x in θ
            update!(x, g[x]*η)
        end
        if i%100==0
            @show loss(x_train, y_train)
        end
    end

    ## Predict the RMSE on the test set
    err = meansquarederror(predict(x_test, m),y_test)
    println(err)
end

# The function above initializes the model’s parameters *W* and *b* randomly.
# Then, it sets the learning rate η and θ as a
# [params object](https://fluxml.ai/Flux.jl/stable/training/training/#Flux.params)
# that points to  W and b. Also, it sets a
# [custom training loop](https://fluxml.ai/Flux.jl/stable/training/training/#Custom-Training-loops)
# which is the [Gradient descent algorithm](https://en.wikipedia.org/wiki/Gradient_descent).
# Finally, it computes the MSE for the test set.

# ## Run the example
# We call the `train` function to run the Housing data example:

cd(@__DIR__)
train()


text:
import KnetLayers: IndexedDict
@testset "data" begin
    data = "this is an example data with example words"
    vocab = IndexedDict{String}()
    append!(vocab,split(data))
    vocab["example"]
    vocab[4]
    vocab[[1,2,3,4]]
    vocab[["example","this"]]
    push!(vocab,"a-new-word")
    vocab["a-new-word"]
    length(vocab)
    @test true

    #https://github.com/denizyuret/Knet.jl/pull/374
    #x = arrtype(rand(10))
    #indices = rand([0,1],10)
    #l2 = BCELoss()
    #l2(x,indices)
    #@test true

    #x = arrtype(randn(10))
    #indices = rand([-1,1],10)
    #l3 = LogisticLoss()
    #l3(x,indices)
    #@test true

end


text:
df_gens = CSV.read(".//inputs//data_generators.csv", header=true, DataFrame);
df_peakers = CSV.read(".//inputs//data_peakers.csv", header=true, DataFrame);
df_fuel_gens = CSV.read(".//inputs//data_fuel_price.csv", header=true, DataFrame);
df_fuel_peakers = CSV.read(".//inputs//data_fuel_price_peakers.csv", header=true, DataFrame);
# Reads the schedules into a dataframe
df_gen_sch = CSV.read(".//inputs//UnitSchedule_(DEC-DEP 2019).csv",
                  dateformat="mm-dd-yyyy THH:MM:SS.000Z", header=true, copycols=true, DataFrame);

df_nox_winter = CSV.read(".//inputs//NOX_Breakpoints_Winter_Results.csv",
                  header=true, copycols=true, DataFrame);

df_so2_winter = CSV.read(".//inputs//SO2_Breakpoints_Winter_Results.csv",
            header=true, copycols=true, DataFrame);

# Source: https://www.eia.gov/electricity/annual/
df_eia_co2 = CSV.read(".//inputs//eia_co2_emiss_per_mmbtu.csv", header=true, DataFrame);

select!(df_gens, Not(:Cogen));

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
struct GraphConv{A<:AbstractArray{Float32},B<:AbstractArray{Float32},C<:AbstractArray{Float32}}
    weight::A # Fout x Fin x k
    bias::B # Fout
    L::C #Graph operator V x V
    k::Integer
    #in_channel::Integer
    #out_channel::Integer
end

function GraphConv(L::AbstractArray,ch::Pair{<:Integer,<:Integer}, k::Integer;
    init = glorot_normal, bias::Bool=false)
    b = bias ? init(ch[2]) : zeros(Float32,ch[2])
    GraphConv(init(ch[2], ch[1], k), b, L, k)
end

# see https://fluxml.ai/Flux.jl/stable/models/advanced/
# https://discourse.julialang.org/t/implementing-a-custom-layer/38632
@functor GraphConv #(weight, bias,)
Flux.trainable(c::GraphConv) = (c.weight,c.bias,)

function (c::GraphConv)(X::AbstractArray)
    fin, b, v = size(X) #Features_in, Batch, Vertices
    fout = size(c.weight,1)
    Y = view(c.weight,:,:,1) * reshape(X, fin, b*v) #Fout x B*V
    Z_prev = reshape(X, fin*b, v) #B*Fin x V
    Z = Z_prev * c.L #B*Fin x V
    Y += view(c.weight,:,:,2) * reshape(Z, fin, b*v)
    for k = 3:c.k
        Z, Z_prev = 2*Z*c.L - Z_prev, Z #B*Fin x V
        Y += view(c.weight,:,:,k) * reshape(Z, fin, b*v)
    end
    out = reshape(Y,fout, b, v) .+ c.bias
    relu.(out)
end

struct GraphMaxPool
    k::Integer
end

function (m::GraphMaxPool)(X::AbstractArray{T}) where {T<:Float32}
    #size(X) Feat x B x V
    permutedims(maxpool(permutedims(X,[3,1,2]),(m.k,)) , [2,3,1])#Feat x B x V/k
end

function graphflatten(x::AbstractArray)
    #size(x) Feat x B x V
    x = permutedims(x, [1,3,2])# Feat x V x B
    return reshape(x, :, size(x)[end])# Feat*V x B
end


text:
"""
测试温度流求导
"""

using Test
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation

@testset "温度流" begin
    model = common_square_lattice(0.)
    Γ4 = TFGamma4(model, 4.0, 8, 40)
    Γ4.V .+= 1.0
    bubb_pp, bubb_fs, bubb_ex = all_bubble_tf_mt(Γ4, 1.0)
    @time dl = dl_tf_mt(Γ4, bubb_pp, bubb_fs, bubb_ex)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 2], -0.0127333, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 3], -0.0038054, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 4], 0.0143136, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 2, 2, 5], 0.0441989, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 2], 0.0181191, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 3], 0.0038054, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 4], -0.0089279, atol=1e-6)
    @test isapprox(dl[1, 1, 1, 1, 4, 2, 5], 0., atol=1e-6)
end


text:
using JuLIPMaterials: Vec2, Mat2, Ten22, Ten42

using Einsum, StaticArrays

export GreenFunction2D, IsoGreenFcn2D, grad

# ========== Implementation of the 2D Green's Function =============

struct GreenFunction2D{T} <: AbstractGreenFunction{T}
   Nquad::Int
   C::Ten42{T}
   remove_singularity::Bool
end

# 2D rotation through π/2
Rot = [0 -1; 1  0];

"""
`GreenFunction2D`
construct a green's function type
"""
function GreenFunction(C::Ten42; Nquad = nothing, remove_singularity = true)
   if Nquad == nothing
      error("still need to implement auto-tuning, please provide Nquad")
   end
   return GreenFunction2D(Nquad, C, remove_singularity)
end

# GreenFunction{T}(C::Array{T, 4}; kwargs...) = GreenFunction(Ten42{T}(C); kwargs...)
#
# GreenFunction(at::AbstractAtoms; kwargs...) =
#       GreenFunction(calculator(at), at; kwargs...)
#
# GreenFunction(calc::AbstractCalculator, at::AbstractAtoms; kwargs...) =
#       GreenFunction(elastic_moduli(calc, at); kwargs...)

function (G::GreenFunction2D{T})(x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 2, 2)
   end
   return eval_green(Vec2(x), G.C, G.Nquad)
end

function grad(G::GreenFunction2D{T}, x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SArray zeros(T, 2, 2, 2)
   end
   return grad_green(Vec2(x), G.C, G.Nquad)
end

"eval_green(x::Vec2, ℂ::Ten42, Nquad::Int)"
function eval_green(x::Vec2{TT}, ℂ::Ten42, Nquad::Int) where TT
   # allocate
   G = @SMatrix zeros(TT, 2, 2)
   zz = @MMatrix zeros(TT, 2, 2)
   # Initialise tensors.
   x̂ = x/norm(x)
   # two vectors orthogonal to x.
   x⟂ = Rot*x̂

   if iseven(Nquad)
      Nquad += 1
   end
   # Integrate
   for ω in range(0.0, step=pi/Nquad, length=Nquad)
      z = cos(ω) * x̂ + sin(ω) * x⟂
      @einsum zz[i,j] = z[α] * ℂ[i,α,j,β] * z[β]
      # Perform integration
      G += inv(zz)*log(dot(x,z))
   end
   # Normalise appropriately
   return G / (2*pi^2*Nquad)
end

struct IsoGreenFcn2D{T} <: AbstractGreenFunction{T}
   λ::T
   μ::T
   remove_singularity::Bool
end

IsoGreenFcn2D(λ, μ; remove_singularity = true) =
   IsoGreenFcn2D(λ, μ, remove_singularity)

function (G::IsoGreenFcn2D{T})(x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 2, 2)
   end
   return eval_greeniso(Vec2(x), G.λ, G.μ)
end

function grad(G::IsoGreenFcn2D{T}, x) where T
   if G.remove_singularity && norm(x) < 1e-10
      return @SArray zeros(T, 2, 2, 2)
   end
   return grad_greeniso(Vec2{T}(x), G.λ, G.μ)
end

"isotropic CLE Green's function"
function eval_greeniso(x::Vec2{T}, λ::Real, μ::Real) where T
   Id = one(SMatrix{2,2,T})
   return ( (λ+μ) * ( x*x'/dot(x,x)- 0.5 * Id ) -
            (λ+3*μ) * log(norm(x)) ) / (4.0*π*μ * (λ+2*μ))
end