text:
module Allocs

# Most of this file was copied from the PProf.jl package, and then adapted to
# export a profile of the heap profile data from this package.
# This code is pretty hacky, and I could probably do a better job re-using
# logic from the PProf package, but :shrug:.

import Profile  # For Profile.Allocs structures

# Import the PProf generated protobuf types from the PProf package:
import PProf
using PProf.perftools.profiles: ValueType, Sample, Function, Location, Line, Label
using PProf: _enter!, _escape_name_for_pprof
const PProfile = PProf.perftools.profiles.Profile
using Base.StackTraces: StackFrame

using PProf.ProtoBuf
using PProf.OrderedCollections

using ProgressMeter

"""
PProf.Allocs.pprof([alloc_profile]; kwargs...)

The `kwargs` are the same as [`PProf.pprof`](@ref), except:
- `frame_for_type = true`: If true, add a frame to the FlameGraph for the Type:
  of every allocation. Note that this tends to make the Graph view harder to
  read, because it's over-aggregated, so we recommend filtering out the `Type:`
  nodes in the PProf web UI.
"""
function pprof(alloc_profile::Profile.Allocs.AllocResults = Profile.Allocs.fetch()
               ;
               web::Bool = true,
               webhost::AbstractString = "localhost",
               webport::Integer = 62261,  # Use a different port than PProf (chosen via rand(33333:99999))
               out::AbstractString = "alloc-profile.pb.gz",
               from_c::Bool = true,
               drop_frames::Union{Nothing, AbstractString} = nothing,
               keep_frames::Union{Nothing, AbstractString} = nothing,
               ui_relative_percentages::Bool = true,
               full_signatures::Bool = true,
               # Allocs-specific arguments:
               frame_for_type::Bool = true,
            )
    period = UInt64(0x1)

    @assert !isempty(basename(out)) "`out=` must specify a file path to write to. Got unexpected: '$out'"
    if !endswith(out, ".pb.gz")
        out = "$out.pb.gz"
        @info "Writing output to $out"
    end

    string_table = OrderedDict{AbstractString, Int64}()
    enter!(string) = _enter!(string_table, string)
    enter!(::Nothing) = _enter!(string_table, "nothing")
    ValueType!(_type, unit) = ValueType(_type = enter!(_type), unit = enter!(unit))

    # Setup:
    enter!("")  # NOTE: pprof requires first entry to be ""

    funcs_map  = Dict{String, UInt64}()
    functions = Vector{Function}()

    locs_map  = Dict{StackFrame, UInt64}()
    locations = Vector{Location}()

    sample_type = [
        ValueType!("allocs", "count"), # Mandatory
        ValueType!("size", "bytes")
    ]

    prof = PProfile(
        sample = [], location = [], _function = [],
        mapping = [], string_table = [],
        sample_type = sample_type,
        # We default to allocs, since the Profile.Allocs code currently uniformly samples
        # accross allocations, so allocs is a representative profile, while size is not.
        default_sample_type = 1, # allocs
        period = period, period_type = ValueType!("heap", "bytes")
    )

    if drop_frames !== nothing
        prof.drop_frames = enter!(drop_frames)
    end
    if keep_frames !== nothing
        prof.keep_frames = enter!(keep_frames)
    end

    function maybe_add_location(frame::StackFrame)::UInt64
        return get!(locs_map, frame) do
            loc_id = UInt64(length(locations) + 1)

            # Extract info from the location frame
            (function_name, file_name, line_number) =
                string(frame.func), string(frame.file), frame.line

            # Decode the IP into information about this stack frame
            function_id = get!(funcs_map, function_name) do
                func_id = UInt64(length(functions) + 1)

                # Store the function in our functions dict
                funcProto = Function()
                funcProto.id = func_id
                file = function_name
                simple_name = _escape_name_for_pprof(function_name)
                local full_name_with_args
                if frame.linfo !== nothing && frame.linfo isa Core.MethodInstance
                    linfo = frame.linfo::Core.MethodInstance
                    meth = linfo.def
                    file = string(meth.file)
                    io = IOBuffer()
                    Base.show_tuple_as_call(io, meth.name, linfo.specTypes)
                    name = String(take!(io))
                    full_name_with_args = _escape_name_for_pprof(name)
                    funcProto.start_line = convert(Int64, meth.line)
                else
                    # frame.linfo either nothing or CodeInfo, either way fallback
                    file = string(frame.file)
                    full_name_with_args = _escape_name_for_pprof(string(frame.func))
                    funcProto.start_line = convert(Int64, frame.line) # TODO: Get start_line properly
                end
                isempty(simple_name) && (simple_name = "[unknown function]")
                isempty(full_name_with_args) && (full_name_with_args = "[unknown function]")
                # WEIRD TRICK: By entering a separate copy of the string (with a
                # different string id) for the name and system_name, pprof will use
                # the supplied `name` *verbatim*, without pruning off the arguments.
                # So even when full_signatures == false, we want to generate two `enter!` ids.
                funcProto.system_name = enter!(simple_name)
                if full_signatures
                    funcProto.name = enter!(full_name_with_args)
                else
                    funcProto.name = enter!(simple_name)
                end
                file = Base.find_source_file(file_name)
                file = file !== nothing ? file : file_name
                funcProto.filename   = enter!(file)
                push!(functions, funcProto)

                return func_id
            end

            locationProto = Location(;id = loc_id,
                                line=[Line(function_id = function_id, line = line_number)])
            push!(locations, locationProto)

            return loc_id
        end
    end

    type_name_cache = Dict{Any,String}()

    function get_type_name(type::Any)
        return get!(type_name_cache, type) do
            return "Alloc: $(type)"
        end
    end

    function construct_location_for_type(typename)
        # TODO: Lol something less hacky than this:
        return maybe_add_location(StackFrame(get_type_name(typename), "nothing", 0))
    end

    # convert the sample.stack to vector of location_ids
    @showprogress "Analyzing $(length(alloc_profile.allocs)) allocation samples..." for sample in alloc_profile.allocs
        # for each location in the sample.stack, if it's the first time seeing it,
        # we also enter that location into the locations table
        location_ids = UInt64[
            maybe_add_location(frame)
            for frame in sample.stacktrace if (!frame.from_c || from_c)
        ]

        if frame_for_type
            # Add location_id for the type:
            pushfirst!(location_ids, construct_location_for_type(sample.type))
        end

        # report the value: allocs = 1 (count)
        # report the value: size (bytes)
        value = [
            1,                   # allocs
            sample.size,         # bytes
        ]

        labels = Label[
            Label(key = enter!("bytes"), num = sample.size, num_unit = enter!("bytes")),
            Label(key = enter!("type"), str = enter!(_escape_name_for_pprof(string(sample.type))))
        ]

        push!(prof.sample, Sample(;location_id = location_ids, value = value, label = labels))
    end

    # Build Profile
    prof.string_table = collect(keys(string_table))
    # If from_c=false funcs and locs should NOT contain C functions
    prof._function = functions
    prof.location  = locations

    # Write to disk
    open(out, "w") do io
        writeproto(io, prof)
    end

    if web
        PProf.refresh(webhost = webhost, webport = webport, file = out,
                      ui_relative_percentages = ui_relative_percentages,
        )
    end

    out
end

end  # module Allocs


text:
## A small driver file to run the portfolio tests and combine
#run like this
#julia -p 4 -L testHarness_Paper.jl tests_port.jl
#ARGS[1] is numRun per batch

using Distributed

const spath = "./Results/portExpOrig"
const param_path = "./Results/param_portExp_mtn2.csv"

n_grid = [2^i for i = 5:17]
numRuns = parse(Int, ARGS[1])

##VG
##original experimentswere run with GammaMin, Gammamax = 1, 20.
##  The "Big" Experiment is 5, 40 step = .01
## The Big100 Experiment is 10, 100 step = .5
start_time = time_ns()
a = @spawn test_ReadData(spath, numRuns, n_grid, 8675309, param_path, Gamma_min=1, Gamma_max=20., Gamma_step =.01)
b = @spawn test_ReadData(spath, numRuns, n_grid, 5164174290, param_path, Gamma_min=1, Gamma_max=20., Gamma_step=.01)
c = @spawn test_ReadData(spath, numRuns, n_grid, 123456, param_path, Gamma_min=1, Gamma_max=20., Gamma_step =.01)
d = @spawn test_ReadData(spath, numRuns, n_grid, 5167462266, param_path, Gamma_min=1, Gamma_max=20., Gamma_step =.01)

######
file_a = fetch(a)
file_b = fetch(b)
file_c = fetch(c)
file_d = fetch(d)

time_stamp = (time_ns() - start_time) * 1e-9

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)

data_t[:, 1] .+= numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] .+= 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] .+= 3numRuns
data = vcat(data, data_t)

#strip the name of file_a to make the numbers better
println("Filea \t", file_a)
println("spath \t", spath)

#indx = search(file_a, spath)[end] + 1
indx =  findfirst(spath, file_a)[end] + 1

f = open("$(spath)_$(file_a[indx:end])_full_$(4numRuns).csv", "w")
writedlm(f,  header, ',')
writedlm(f,  data, ',')
close(f)

println("Num Paths: \t $(numRuns) \t Time:", time_stamp )

ext: jl

max_stars_repo_path: Normal Model/testPortExp.jl

max_issues_repo_path: Normal Model/testPortExp.jl

max_forks_repo_path: Normal Model/testPortExp.jl


text:
function passTable(simulation, semesters=-1)
    frame = DataFrame()

    # Make Keys
    frame[:COUSE] = []
    terms = simulation.duration

    if semesters > -1
        terms = semesters
    end

    for i=1:terms
        frame[Symbol("TERM$(i)")] = []
    end

    # Populate data
    for course in simulation.curriculum.courses
        row = [course.name]
        prev = 0
        for i=1:terms
            prev += course.termpassed[i]
            row = [row round((prev/simulation.numStudents)*100, 3)]
        end
        push!(frame, row)
    end

    rates = ["GRAD RATE"]
    for i=1:terms
        rates = [rates round(simulation.termGradRates[i]*100, 3)]
    end
    push!(frame, rates)

    frame
end

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
const TIMES = 2
const TNUM = 0
const sym = ['\$', '+', '-', ':' , '*']
mutable struct Ftime
    t::Int
    num::Int
end

ftime() = Ftime(TIMES, TNUM)
@inline function fadd(t::Ftime)
	 		if t.t >= 600
				t.num += 1
			else
				t.t += TIMES
			end
		end
@inline finit(t::Ftime) = t.t > 2 && begin t.t, t.num = TIMES, TNUM end

@inline redis_test(conn::RedisConnection,test::AbstractString) = send_command(conn, echo(test) )

function redis_collect(conn::RedisConnection , test::AbstractString="test")

	res = redis_test(conn, test)
	collects(conn,test , [])

end

function collects(conn::RedisConnection , test::AbstractString, data::Vector)
	tmp = readline(conn.socket)
	syms, value = tmp[1] , tmp[2:end]
	if (syms in sym)
		val = reply(redisreply{Symbol(syms)}, value, conn.socket)
		if val != test
			push!(data, val )
		else
			return data
		end
	else
		push!(data, tmp)
	end
	collects(conn, test, data)
end

function reline(conn::RedisConnectionBase, times::Ftime, fun::Function)
    println("Failed to connect to Redis server ,Reconnect after $(times.t) seconds .")
    sleep(times.t)
    try
        reconnect(conn)
        println("""
			Success to connect to Redis server ,
			Time consuming greater than or equal to $(fun(times)) seconds .
			""")
    catch
        fadd(times)
        reline(conn, times, fun)
    end

end

macro cheak_reline(conn )

	function f(x::Ftime)::Int
		if x.num >0
			x.num * 600
		else
			(x.t/4)*x.t + (x.t/2 )
		end
	end

    con = esc(conn)
    quote
        try
            ping( $con )
        catch
            reline( $con , ftime(), $f)
        end
     end

end

macro genmacro(funname, lenfun, popfun)

	func = Expr(:call , funname, :redis,  :key ,:fun, :batch )

	nums = Expr(:(=), :num , Expr(:call, :(|>),
		    Expr(:call, lenfun, Expr(:$, :redis), Expr(:$, :key)) ,
		    Expr(:(->), :q, Expr(:block,
			Expr(:if , Expr(:call, :(>=), :q, Expr(:$, :batch)) ,
			    Expr(:$, :batch) , :q)))))

	expr = Expr(:call, :pipelines, Expr(:$, :redis),
		    Expr(:call, :repeat , Expr(:call, popfun , Expr(:$, :key)) ,:num), :num)

	body =  Expr(:block,
		    Expr(:(=), :freq, Expr(:call, :ftime)),
		    Expr(:while , true,
			Expr(:block, Expr(:macrocall, Symbol("@cheak_reline"), "", Expr(:$, :redis)) ,
			nums ,
			Expr(:if , Expr(:call, :(>=), :num ,1),
					    Expr(:try, Expr(:block,
						    Expr(:global, Expr(:(=), :data , expr )),
						    Expr(:call, Expr(:$, :fun), :data),
						    Expr(:(=), :data, 0 ),
						    Expr(:call, :finit, :freq )),
						   false ,
						Expr(:block,
						    Expr(:(&&), Expr( :call, :!, Expr(:call , :(isa) , :data , :Int)),
						    Expr(:call, Expr(:$, :fun), :data) ),
							Expr(:(=) , :tmp , Expr(:call, :redis_collect , Expr(:$, :redis))),
							Expr(:(&&), Expr(:call , :(>=) , Expr(:call, :length, :tmp) , 1),
							Expr(:call, Expr(:$, :fun), :tmp) ),
							Expr(:call, :disconnect, Expr(:$, :redis)))),
					    Expr(:block, Expr(:call, :sleep, Expr(:call, :fadd, :freq)))))))

		 esc( Expr( :macro , func,   Expr(:block, Expr(:call, :esc,Expr(:quote, body)))))

end

@genmacro spop scard spop
@genmacro lpop llen lpop
@genmacro rpop llen rpop

max_stars_repo_name: frankill/jredis

max_issues_repo_name: frankill/jredis

max_forks_repo_name: frankill/jredis


text:
# Statistical Models

@compat abstract type StatisticalModel end

coef(obj::StatisticalModel) = error("coef is not defined for $(typeof(obj)).")
coeftable(obj::StatisticalModel) = error("coeftable is not defined for $(typeof(obj)).")
confint(obj::StatisticalModel) = error("coefint is not defined for $(typeof(obj)).")

"""
    deviance(obj::StatisticalModel)

Returns the deviance of the model relative to a reference, which is usually when applicable
the saturated model. It is equal, *up to a constant*, to `-2 log L`, with `L`
the likelihood of the model.
"""
deviance(obj::StatisticalModel) = error("deviance is not defined for $(typeof(obj)).")

"""
    nulldeviance(obj::StatisticalModel)

Returns the deviance of the null model, that is the one including only the intercept.
"""
nulldeviance(obj::StatisticalModel) = error("nulldeviance is not defined for $(typeof(obj)).")
loglikelihood(obj::StatisticalModel) = error("loglikelihood is not defined for $(typeof(obj)).")
nullloglikelihood(obj::StatisticalModel) = error("nullloglikelihood is not defined for $(typeof(obj)).")

"""
    nobs(obj::StatisticalModel)

Returns the number of independent observations on which the model was fitted. Be careful
when using this information, as the definition of an independent observation may vary
depending on the model, on the format used to pass the data, on the sampling plan
(if specified), etc.
"""
nobs(obj::StatisticalModel) = error("nobs is not defined for $(typeof(obj)).")

"""
    dof(obj::StatisticalModel)

Returns the number of degrees of freedom consumed in the model, including
when applicable the intercept and the distribution's dispersion parameter.
"""
dof(obj::StatisticalModel) = error("dof is not defined for $(typeof(obj)).")
stderr(obj::StatisticalModel) = sqrt(diag(vcov(obj)))
vcov(obj::StatisticalModel) = error("vcov is not defined for $(typeof(obj)).")
fit(obj::StatisticalModel, data...) = error("fit is not defined for $(typeof(obj)).")
fit!(obj::StatisticalModel, data...) = error("fit! is not defined for $(typeof(obj)).")

"""
    aic(obj::StatisticalModel)

Akaike's Information Criterion, defined as `-2 log L + 2k`, with `L` the likelihood
of the model, and `k` its number of consumed degrees of freedom (as returned by `dof`).
"""
aic(obj::StatisticalModel) = -2loglikelihood(obj) + 2dof(obj)

"""
    aicc(obj::StatisticalModel)

Corrected Akaike's Information Criterion for small sample sizes (Hurvich and Tsai 1989),
defined as `-2 log L + 2k + 2k(k-1)/(n-k-1)`, with `L` the likelihood of the model,
`k` its number of consumed degrees of freedom (as returned by `dof`), and `n` the number
of observations (as returned by `nobs`).
"""
function aicc(obj::StatisticalModel)
    k = dof(obj)
    n = nobs(obj)
    -2loglikelihood(obj) + 2k + 2k*(k+1)/(n-k-1)
end

"""
    bic(obj::StatisticalModel)

Bayesian Information Criterion, defined as `-2 log L + k log n`, with `L`
the likelihood of the model,  `k` its number of consumed degrees of freedom
(as returned by `dof`), and `n` the number of observations (as returned by `nobs`).
"""
bic(obj::StatisticalModel) = -2loglikelihood(obj) + dof(obj)*log(nobs(obj))

"""
    r2(obj::StatisticalModel, variant::Symbol)
    r²(obj::StatisticalModel, variant::Symbol)

Coefficient of determination (R-squared).

For a linear model, the R² is defined as `ESS/TSS`, with `ESS` the explained sum of squares
and `TSS` the total sum of squares, and `variant` can be omitted.

For other models, one of several pseudo R² definitions must be chosen via `variant`.
Supported variants are:
- `:MacFadden` (a.k.a. likelihood ratio index), defined as `1 - log L/log L0`.
- `:CoxSnell`, defined as `1 - (L0/L)^(2/n)`
- `:Nagelkerke`, defined as `(1 - (L0/L)^(2/n))/(1 - L0^(2/n))`, with `n` the number
of observations (as returned by `nobs`).

In the above formulas, ``L`` is the likelihood of the model, ``L0`` that of the null model
(the model including only the intercept). These two quantities are taken to be minus half
`deviance` of the corresponding models.
"""
function r2(obj::StatisticalModel, variant::Symbol)
    ll = -deviance(obj)/2
    ll0 = -nulldeviance(obj)/2

    if variant == :McFadden
        1 - ll/ll0
    elseif variant == :CoxSnell
        1 - exp(2/nobs(obj) * (ll0 - ll))
    elseif variant == :Nagelkerke
        (1 - exp(2/nobs(obj) * (ll0 - ll)))/(1 - exp(2/nobs(obj) * ll0))
    else
        error("variant must be one of :McFadden, :CoxSnell or :Nagelkerke")
    end
end

const r² = r2

"""
    adjr2(obj::StatisticalModel, variant::Symbol)
    adjr²(obj::StatisticalModel, variant::Symbol)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as `1 - (1 - (1-R²)(n-1)/(n-p))`, with `R²`
the coefficient of determination, `n` the number of observations, and `p` the number of
coefficients (including the intercept). This definition is generally known as the Wherry Formula I.

For other models, one of the several pseudo R² definitions must be chosen via `variant`.
The only currently supported variant is `:MacFadden`, defined as `1 - (log L - k)/log L0`.
In this formula, `L` is the likelihood of the model, `L0` that of the null model
(the model including only the intercept). These two quantities are taken to be minus half
`deviance` of the corresponding models. `k` is the number of consumed degrees of freedom
of the model (as returned by `dof`).
"""
function adjr2(obj::StatisticalModel, variant::Symbol)
    ll = -deviance(obj)/2
    ll0 = -nulldeviance(obj)/2
    k = dof(obj)

    if variant == :McFadden
        1 - (ll - k)/ll0
    else
        error(":McFadden is the only currently supported variant")
    end
end

const adjr² = adjr2

@compat abstract type RegressionModel <: StatisticalModel end

fitted(obj::RegressionModel) = error("fitted is not defined for $(typeof(obj)).")
model_response(obj::RegressionModel) = error("model_response is not defined for $(typeof(obj)).")
residuals(obj::RegressionModel) = error("residuals is not defined for $(typeof(obj)).")

"""
    predict(obj::RegressionModel, [newX])

Form the predicted response of model `obj`. An object with new covariate values `newX` can be supplied,
which should have the same type and structure as that used to fit `obj`; e.g. for a GLM
it would generally be a `DataFrame` with the same variable names as the original predictors.
"""
predict(obj::RegressionModel) = error("predict is not defined for $(typeof(obj)).")
predict!(obj::RegressionModel) = error("predict! is not defined for $(typeof(obj)).")
dof_residual(obj::RegressionModel) = error("dof_residual is not defined for $(typeof(obj)).")

## coefficient tables with specialized show method

## Nms are the coefficient names, corresponding to rows in the table
type CoefTable
    cols::Vector
    colnms::Vector
    rownms::Vector
    function CoefTable(cols::Vector,colnms::Vector,rownms::Vector)
        nc = length(cols)
        nrs = map(length,cols)
        nr = nrs[1]
        length(colnms) in [0,nc] || error("colnms should have length 0 or $nc")
        length(rownms) in [0,nr] || error("rownms should have length 0 or $nr")
        all(nrs .== nr) || error("Elements of cols should have equal lengths, but got $nrs")
        new(cols,colnms,rownms)
    end

    function CoefTable(mat::Matrix,colnms::Vector,rownms::Vector,pvalcol::Int=0)
        nc = size(mat,2)
        cols = Any[mat[:, i] for i in 1:nc]
        if pvalcol != 0                         # format the p-values column
            cols[pvalcol] = [PValue(cols[pvalcol][j])
                            for j in eachindex(cols[pvalcol])]
        end
        CoefTable(cols,colnms,rownms)
    end
end

type PValue
    v::Number
    function PValue(v::Number)
        0. <= v <= 1. || isnan(v) || error("p-values must be in [0.,1.]")
        new(v)
    end
end

function show(io::IO, pv::PValue)
    v = pv.v
    if isnan(v)
        @printf(io,"%d", v)
    elseif v >= 1e-4
        @printf(io,"%.4f", v)
    else
        @printf(io,"<1e%2.2d", ceil(Integer, max(nextfloat(log10(v)), -99)))
    end
end

function show(io::IO, ct::CoefTable)
    cols = ct.cols; rownms = ct.rownms; colnms = ct.colnms;
    nc = length(cols)
    nr = length(cols[1])
    if length(rownms) == 0
        rownms = [lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
    end
    rnwidth = max(4,maximum([length(nm) for nm in rownms]) + 1)
    rownms = [rpad(nm,rnwidth) for nm in rownms]
    widths = [length(cn)::Int for cn in colnms]
    str = String[isa(cols[j][i], AbstractString) ? cols[j][i] :
        sprint(showcompact,cols[j][i]) for i in 1:nr, j in 1:nc]
    for j in 1:nc
        for i in 1:nr
            lij = length(str[i,j])
            if lij > widths[j]
                widths[j] = lij
            end
        end
    end
    widths .+= 1
    println(io," " ^ rnwidth *
            join([lpad(string(colnms[i]), widths[i]) for i = 1:nc], ""))
    for i = 1:nr
        print(io, rownms[i])
        for j in 1:nc
            print(io, lpad(str[i,j],widths[j]))
        end
        println(io)
    end
end

"""
	ConvergenceException(iters::Int)

The fitting procedure failed to converge in `iters` number of iterations.
"""
type ConvergenceException <: Exception
    iters::Int
end

Base.showerror(io::IO, ce::ConvergenceException) = print(io, "failure to converge after $(ce.iters) iterations")


text:
module BasicTest

using Talkon
using ReTest

@testset "start" begin
    @test true
end

end # module


text:
"""Reset the constants to allow for arbitrary number of fraction bits."""
function set_nfrac(n::Int)
    @assert (n > 6) && (n<max_nfrac_supported+1) "Only 6...11 fraction bits supported, $n given."

    # change pointer references based on n
    nfrac[] = n
    nint[] = 15-nfrac[]
    scale[] = 2^nfrac[]
    scale_over_logof2[] = scale[]/log(2f0)
    max_diff_resolvable[] = find_max_diff_res(scale[])

    # change Float32 lookup table for conversion LogFixPoint16 -> Float32
    f32lookup[:] = createF32LookupTable(nint[],nfrac[])

    # change addition and subtraction table
    addTable[:] = createAddLookup(scale[])
    subTable[:] = createSubLookup(scale[])

    if rounding_mode[] == :lin
        c_b[] = rounding_correction(scale[])
    else
        c_b[] = 0
    end

    @warn "LogFixPoint16 was changed to $(nint[]) integer and $(nfrac[]) fraction bits."

    return nothing
end

"""Change the rounding mode from round-to-nearest in either linear or logarithmic space."""
function set_rounding_mode(mode::Symbol=:lin)
    @assert mode in [:lin,:log] "only mode :lin or :log allowed."

    rounding_mode[] = mode

    if rounding_mode[] == :lin
        c_b[] = rounding_correction(scale[])
    else
        c_b[] = 0
    end

    # recalculate lookup tables
    max_diff_resolvable[] = find_max_diff_res(scale[])
    f32lookup[:] = createF32LookupTable(nint[],nfrac[])
    addTable[:] = createAddLookup(scale[])
    subTable[:] = createSubLookup(scale[])

    @warn "LogFixPoint16 rounding mode changed to round to nearest in $(string(mode))-space."
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# DynamicBounds.jl
# A package for compute bounds and relaxations of the solutions of
# parametric differential equations.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/DynamicBounds.jl
# Main module code.
#############################################################################

module DynamicBounds

using DocStringExtensions

include("reexport_using.jl")

@reexport using DynamicBoundsBase: HasStateBounds, HasConstantStateBounds, HasVariableStateBounds,
       HasUserJacobian, ConstantStateBounds, VariableStateBounds,
       PolyhedralConstraint, AbstractDERelaxProblem, AbstractRelaxProblemAttribute,
       Nominal, Lower, Upper, Undefined, IntegratorName, Value, Gradient,
       Subgradient, Bound, Relaxation,IsNumeric, IsSolutionSet, TerminationStatus,
       ParameterValue, ParameterBound, SupportSet, TerminationStatusCode, relax!,
       integrate!, AbstractDERelaxIntegrator, AbstractODERelaxIntegrator,
       set!, setall!,getall!, supports, ODERelaxProb

@reexport using DynamicBoundspODEsIneq: DifferentialInequality, DifferentialInequalityCond,
                                        DifferentialInequalityAffect, DifferentialInequalityAffectNeg,
                                        DifferentialInequalityf

import DynamicBoundsBase
const DBB = DynamicBoundsBase

export get
function DBB.get(integrator::DBB.AbstractDERelaxIntegrator, attr::DBB.AbstractIntegratorAttribute, idxs::Vector)
    DBB.get(integrator, attr, idxs)
end
function DBB.get(problem::DBB.AbstractDERelaxProblem, attr::DBB.AbstractRelaxProblemAttribute, idxs::Vector)
    DBB.get(problem, attr, idxs)
end
function DBB.get(m::DBB.AbstractDERelaxIntegrator, attr::DBB.AbstractIntegratorAttribute, args...)
    DBB.get(m, attr, args...)
end
function DBB.get(m::DBB.AbstractDERelaxProblem, attr::DBB.AbstractRelaxProblemAttribute, args...)
    DBB.get(m, attr, args...)
end

include("library.jl")

end


text:
"""
    Crystal{D} <: AbstractCrystal{D}

A crystal structure in `D` dimensions. Contains information about the lattice, space group, and
atoms.

Atoms can be specified in terms of the generating set, which can be used to fill a cell given
space group information, or in terms of explicitly given atomic positions.

If `sgno` is set to 0, it will be assumed that coordinates are given in Cartesian coordinates in
units of angstroms. Otherwise, it will be assumed that coordinates are given in terms of the unit
cell basis.
"""
struct Crystal{D} <: AbstractCrystal{D}
    # Lattice information (real space)
    latt::RealLattice{D}
    # Space group number - set to 0 if non-periodic/unknown
    sgno::Int
    # Space group origin (location of inversion center, if present - defaults to zeros)
    orig::SVector{D,Float64}
    # Positions of generating atoms (needed to fill the whole structure given the space group)
    gen::AtomList{D}
    # Positions of atoms explicitly generated (used to generate templates, XYZs, etc.)
    pos::AtomList{D}
    # Inner constructor
    function Crystal{D}(
        latt::AbstractLattice{D},
        sgno::Integer,
        orig::AbstractVector{<:Real},
        gen::AtomList{D},
        pos::AtomList{D}
    ) where D
        # TODO: include some validation
        # If the atom list doesn't have a basis defined (assuming Cartesian coordinates)
        # generate a new AtomList
        if basis(gen) == zeros(BasisVectors{D})
            # Do we want to use the conventional basis vectors all the time?
            # I think it's a good temporary choice, just because if a dataset is being loaded in
            # from a computation, conv(latt) should generally match the primitive cell
            gen = AtomList{D}(conv(latt), reduce_coords(conv(latt), gen, incell=true))
        end
        return new{D}(RealLattice(latt), sgno, orig, gen, pos)
    end
end

# The name `CrystalData{D}` was not used to avoid implying that this is a subtype of
# `AbstractCrystalData{D}` - which is used only for datasets
"""
    CrystalWithDatasets{D,K,V} <: AbstractCrystal{D}

A pairing of a `Crystal{D}` and a `Dict{K,V}` which allows for access to associated datasets.
"""
struct CrystalWithDatasets{D,K,V} <: AbstractCrystalData{D}
    xtal::Crystal{D}
    data::Dict{K,V}
end

# Allow for getting datasets by key; no need to reach into the Dict
# TODO: figure out how autocompletion works and how to enable it here
# UPDATE: seems like it's hard-coded into the REPL for `AbstractDict`
# So adding in that functionality doesn't seem to be feasible
function Base.getindex(xtaldata::CrystalWithDatasets{D,K,V}, key::K) where {D,K,V}
    return xtaldata.data[key]
end

# Easy way of pulling just the crystal from a `CrystalWithDatasets{D}`
function Base.convert(::Type{Crystal{D}}, xtaldata::CrystalWithDatasets{D,K,V}) where {D,K,V}
    return xtaldata.xtal
end

Crystal(xtaldata::CrystalWithDatasets{D,K,V}) where {D,K,V} = xtaldata.xtal
data(xtaldata::CrystalWithDatasets{D,K,V}) where {D,K,V} = xtaldata.data

RealLattice(xtal::Crystal) = xtal.latt
ReciprocalLattice(xtal::Crystal) = ReciprocalLattice(xtal.latt)

prim(xtal::Crystal) = prim(RealLattice(xtal))
conv(xtal::Crystal) = conv(RealLattice(xtal))

basis(xtal::Crystal; primitive::Bool=false) = primitive ? prim(xtal.latt) : conv(xtal.latt)
basis(xtaldata::CrystalWithDatasets; primitive::Bool=false) = basis(xtaldata.xtal; primitive)

volume(xtal::Crystal; primitive::Bool=false) = volume(xtal.latt; primitive)
volume(xtaldata::CrystalWithDatasets; primitive::Bool=false) = volume(xtaldata.xtal; primitive)

# TODO: fix this so that it gets the right number of atoms regardless of space group
"""
    natom_cell(xtal::Crystal)

Returns the number of atoms in a crystal's unit cell.

This function is currently not aware of space groups or settings, so if the generating set does not
contain all the atoms in the cell, it will return the wrong value.
"""
natom_cell(xtal::Crystal) = natom(xtal.gen)

"""
    natom_template(xtal::Crystal)

Returns the number of atoms in the crystal template.
"""
natom_template(xtal::Crystal) = natom(xtal.pos)

#=
"""
    all_atoms_in_cell(spgrp::Integer, l::AtomList{D}; onbounds::Bool = false) -> AtomList{D}

Generates all atoms within a unit cell given the space group information and a minimal set of
atomic positions.
"""
function all_atoms_in_cell(spgrp::Integer, l::AtomList{D}; onbounds::Bool = false) where D
    # Space group 0 effectively means no translation symmetry
    # Space group 1 should be the trivial group (just translations): do nothing
    spgrp in (0,1) && return l
    # TODO: implement this function.
    return l
end
=#


text:
using EngEconomics

# Given
i = 0.1
xInit = 12500
yInit = 8900
xAnnualBenefit = 6800
yAnnualBenefit = 2000
xSalvage = 5000
ySalvage = 8900
xLife = 2
yLife = 3

# Find: Present Worth Analysis of Both

# Determine the least common multiple of xLife and yLife to figure out how many
# times we have to repeat the PW Analysis
years = 6
xTimes = years / xLife
yTimes = years / yLife

# Determine Product X's PW
xInitPW = -xInit * (1 + presentWorthFactor(i, xLife) + presentWorthFactor(i, 2 * xLife))
xAnnualBenefitPW = (xAnnualBenefit * seriesPresentAmountFactor(i, xLife)) * (1 + presentWorthFactor(i, xLife) + presentWorthFactor(i, 2 * xLife))
xSalvagePW = xSalvage * (presentWorthFactor(i, xLife) + presentWorthFactor(i, 2 * xLife) + presentWorthFactor(i, 3 * xLife))
xPW = xInitPW + xAnnualBenefitPW + xSalvagePW

# Determine Product Y's PW
yInitPW = -yInit * (1 + presentWorthFactor(i, yLife))
yAnnualBenefitPW = (yAnnualBenefit * seriesPresentAmountFactor(i, yLife)) * (1 + presentWorthFactor(i, yLife))
ySalvagePW = ySalvage * (presentWorthFactor(i, yLife) + presentWorthFactor(i, 2 * yLife))
yPW = yInitPW + yAnnualBenefitPW + ySalvagePW

# Determine which product is better by comparing they're present value
if xPW > yPW
	println("Pick Product X")
else
	println("Pick Product Y")
end


text:
"""
    ECON 602 Problem Set 1 Question 4
"""

using Roots, Plots

function adequilibrium(π, γ, ωA, ωB, p_0)
    # solve for equilibrium price and consumption
    # p_0 is initial guess for numerical solver
    println("\n AD equilibrium:")
     @show π, γ, ωA, ωB;

    f(p, π, γ, ωA, ωB) = π*(ωA[1] + p*ωB[1]) + 1/(1 + p*(p*π/(1-π))^(-1/γ))*(ωA[2] + p*ωB[2]) - ωA[1] - ωA[2]
    cA(p, π, γ, ωA, ωB) = [π*(ωA[1] + p*ωB[1]), 1/(1 + p*(p*π/(1 - π))^(-1/γ))*(ωA[2] + p*ωB[2])]
    cB(p, π, γ, ωA, ωB) = [(1 - π)/p * (ωA[1] + p*ωB[1]), 1/(p + (p*π/(1 - π))^(1/γ))*(ωA[2] + p*ωB[2])]

    # It never hurts to plot and see what you're working with
    plot(p -> f(p, π, γ, ωA, ωB), 0, 4)
    plot!(x -> 0, 0, 4)

    @show p_star = fzero(p -> f(p, π, γ, ωA, ωB), p_0);
    #=
     The first argument, p -> f(p, π, γ, ωA, ωB), defines and 'anonymous'
     function of only p. It is a 'wrapper' for the function f that assigns
     the values to all other parameters as they are currently defined.
    =#

    @show cA = cA(p_star, π, γ, ωA, ωB);
    @show cB = cB(p_star, π, γ, ωA, ωB);
    return (p_star = p_star, cA = cA, cB = cB)
end

function adequilibrium_alt(π, γ, ωA, ωB, p_min, p_max)
    println("\n AD equilibrium (alternative function):")
    # solve for equilibrium price and consumption
    # alternative syntax and solution algorithm using bounds
    @show π, γ, ωA, ωB;

    f(p, π, γ, ωA, ωB) = π*(ωA[1] + p*ωB[1]) + 1/(1 + p*(p*π/(1-π))^(-1/γ))*(ωA[2] + p*ωB[2]) - ωA[1] - ωA[2]
    cA(p, π, γ, ωA, ωB) = [π*(ωA[1] + p*ωB[1]), 1/(1 + p*(p*π/(1 - π))^(-1/γ))*(ωA[2] + p*ωB[2])]
    cB(p, π, γ, ωA, ωB) = [(1 - π)/p * (ωA[1] + p*ωB[1]), 1/(p + (p*π/(1 - π))^(1/γ))*(ωA[2] + p*ωB[2])]

    # It never hurts to plot and see what you're working with
    plot(p -> f(p, π, γ, ωA, ωB), p_min, p_max)
    plot!(x -> 0, p_min, p_max)

    # Alternative to anonymous function
    f(p) = f(p, π, γ, ωA, ωB)
    #=
     Note that this is a distinct function from the first
     function I defined, as it has a different number of
     arguments
    =#
    @show p_star = fzero(f, p_min, p_max);
    #=
     This provides an allowable range for solutions between
     0.01 and 4, rather than a single initial guess. This
     can be useful to prevent the algorithm from searching
     for negatives prices for example.
    =#

    @show cA = cA(p_star, π, γ, ωA, ωB);
    @show cB = cB(p_star, π, γ, ωA, ωB);
    return (p_star = p_star, cA = cA, cB = cB);
end

# q4e
ωA = [2, 5];
ωB = [5, 2];
γ = 2;
π = 0.5;
q4e = adequilibrium(π, γ, ωA, ωB, 1.1);
q4e_alt = adequilibrium_alt(π, γ, ωA, ωB, 0.01, 10.);

# q4f
ωA = [2, 5];
ωB = [5, 2];
γ = 2;
# π = 0.25
π = 0.25;
q4f_π025 = adequilibrium(π, γ, ωA, ωB, 1.1);
# π = 0.75
π = 0.75;
q4f_π075 = adequilibrium(π, γ, ωA, ωB, 1.1);


text:
module StatsDiscretizations

using Reexport

using Dictionaries
using Discretizers

@reexport using Distributions
import Distributions: pdf, logpdf, cdf, ccdf

using Infinity
using Infinity.Utils

@reexport using IntervalSets

using Requires
@reexport using StatsBase

include("intervals.jl")
include("distributions_intervals.jl")
include("discretizers.jl")
include("discretized_function.jl")
include("default_discretizers.jl")

function __init__()
    @require JuMP="4076af6c-e467-56ae-b986-b466b2749572" include("discretized_function_jump.jl")
end

export RealLineDiscretizer,
       BoundedIntervalDiscretizer,
       DiscretizedFunction,
       FiniteSupportDiscretizer

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 030f5512-2a66-11ec-1a7a-21059e6facac
begin

# Import packages used

using Random
using DataFrames
using CSV
using StatFiles
using Statistics

end

# ╔═╡ af9b5960-2a65-11ec-3015-53185d92610c
# Replication of DellaVigna and Pope, 2018, "What Motivates Effort? Evidence and Expert Forecasts"

# The code in this notebook performs the estimates to replicate columns (1) (3) panel A and columns (1) (4) panel B of Table 5 using minimum distance.

# ╔═╡ 03372860-2a66-11ec-2b31-298b2c9afc9d
# Authors:
# Massimiliano Pozzi (Bocconi University, pozzi.massimiliano@studbocconi.it)
# Salvatore Nunnari (Bocconi University, https://snunnari.github.io, salvatore.nunnari@unibocconi.it)

# ╔═╡ 7c88c030-55f8-11ec-0159-99a5b81bf850
# This notebook works with both Neptune and Pluto; the "begin" and "end" tags in each cell are necessary for Pluto notebooks but do not play any role in Neptune

# ╔═╡ 7dcdea60-55f8-11ec-2fdc-1d683f29a01c
# This notebook was tested with the following packages versions:

# Pozzi: julia 1.6.1, Pluto 0.14.7, Neptune 0.14.0, DataFrames 1.2.2, CSV 0.9.5

# ╔═╡ 03008800-2a66-11ec-0c80-1767685eb4fa
begin

# 1. Data Cleaning and Data Preparation

# We import the relevant dataset containing data on the number of buttonpresses in the different treatments and for different piece rates wage that the participants received when completing the task. We then compute the means for the treatments that are used to compute the minimum distance estimates.

dt = DataFrame(load(raw"../input/mturk_clean_data_short.dta"))

emp_moments = round.(combine(groupby(dt, :treatment), :buttonpresses => mean)[2], digits=0)

end;

# ╔═╡ 944bb54e-55f9-11ec-32b3-c953374ecd8f
# We set up the Bootstrap procedure for the standard errors. The following cell creates 2000 new samples and computes the mean for the relevant treatments.

# ╔═╡ 030039e0-2a66-11ec-3766-ad9593bffbd3
begin

# Treatments are as follows:
# 1.1: benchmark specification with piece rate of 0.01
# 1.2: benchmark specification with piece rate of 0.10
# 1.3: benchmark specification with piece rate of 0.00
# 3.1: social preferences (charity) with piece rate of 0.01
# 3.2: social preferences (charity) with piece rate of 0.10
# 10 : social preferences (gift exchange) bonus of 40 cents (independently of nr buttonpresses)
# 4.1: time discounting with extra 0.01 paid two weeks later
# 4.2: time discounting with extra 0.01 paid four weeks later

# define the vector containing the treatments (used for the loop) and the vectors that will store the rounded mean buttonpresses in our 2000 new samples.
# E11 is a 1x1000 vector that will contain the average buttonpresses in each of our 1000 samples for treatment 1.1 etc.

treatment = ["1.1", "1.2", "1.3", "3.1", "3.2", "10", "4.1", "4.2"]
E11, E12, E13, E31, E32, E10, E41, E42 = Int32[], Int32[], Int32[], Int32[], Int32[], Int32[], Int32[], Int32[]

# We first get a smaller dataframe containing only observations for a specific treatment, we then resample the observations and compute the rounded mean of buttonpresses, save the result and then pass onto the next treatment. we do this 1000 times.

for i=1:1000                 # we want 1000 new samples
    for treat in treatment   # we want to compute the mean for all treatments
        db = copy(dt)        # copy the original dataset
        db = db[db[:,2] .== treat,:] # keep only observation for treatment 'treat'
        bootsample = db[vec(rand(1:size(db[!,1],1),1,size(db[!,1],1))),1] # resample the dataset by taking random rows from the first column
		if treat == "1.1" append!(E11, round(mean(bootsample))) end
        if treat == "1.2" append!(E12, round(mean(bootsample))) end
        if treat == "1.3" append!(E13, round(mean(bootsample))) end
        if treat == "3.1" append!(E31, round(mean(bootsample))) end
        if treat == "3.2" append!(E32, round(mean(bootsample))) end
        if treat == "10"  append!(E10, round(mean(bootsample))) end
        if treat == "4.1" append!(E41, round(mean(bootsample))) end
        if treat == "4.2" append!(E42, round(mean(bootsample))) end
	end
end

end

# ╔═╡ ced0482e-2a66-11ec-06ca-0b961f35a057
# For a more detailed explanation of the model we refer the reader to the python notebook

# Define the function to estimate the parameters using minimum distance:

# There are two specifications: exponential cost function and power cost function
# The parameters are found by imposing the theoretical means found from the agent's problem being equal to the empirical means found in the data.
# E11 up to E42 are the relevant empirical moments
# Specification can either be 'Exp' or 'Power'
# P is a vector containing the different piece-rates
# To simplify the computations the authors take the log of the variables to estimate gamma k and s.

function mindisest(E11,E12,E13,E31,E32,E10,E41,E42,specification)

    P = [0, 0.01, 0.1]

	# gamma, k and s are found by solving the system of N equations in N unknowns

    if specification == "Exp"
        log_k = (log(P[3]) - log(P[2])*(E12)/(E11))/(1 - (E12)/(E11))
        log_gamma = log((log(P[2]) - log_k)/(E11))
        log_s = exp(log_gamma)*(E13) + log_k
	end

    if specification == "Power"
        log_k = (log(P[3]) - log(P[2])*log(E12)/log(E11))/(1 - log(E12)/log(E11))
        log_gamma = log((log(P[2]) - log_k)/log(E11))
        log_s = exp(log_gamma)*log(E13) + log_k
	end

    k = exp(log_k)
    g = exp(log_gamma)
    s = exp(log_s)

    if specification == "Exp"
        EG31, EG32, EG10, EG41, EG42 = exp.(E31*g), exp.(E32*g), exp.(E10*g), exp.(E41*g), exp.(E42*g)

        alpha = 100/9*k*(EG32-EG31)
        a = 100*k*EG31-100*s-alpha
        s_ge = k*EG10 - s
        delta = sqrt((k*EG42-s)/(k*EG41-s))
        beta  = 100*(k*EG41-s)/(delta^2)
	end

    if specification == "Power"
        alpha = 100/9*k*(E32^g-E31^g)
        a = 100*k*E31^g-100*s-alpha
        s_ge = k*E10^g - s
        delta = sqrt((k*E42^g-s)/(k*E41^g-s))
        beta  = 100*(k*E41^g-s)/(delta^2)
	end

    return k, g, s, alpha, a, s_ge, beta, delta
end

# ╔═╡ 9ac87bc0-5608-11ec-3f39-3713d508fdfb
begin

# Table 5 minimum distance estimates: columns (1) (3) panel A and columns (1) (4) panel B

Table5Exp = mindisest(emp_moments[1],emp_moments[2],emp_moments[3],emp_moments[7],emp_moments[8],emp_moments[5],emp_moments[9],emp_moments[10],"Exp")

Table5Power = mindisest(emp_moments[1],emp_moments[2],emp_moments[3],emp_moments[7],emp_moments[8],emp_moments[5],emp_moments[9],emp_moments[10],"Power")

end;

# ╔═╡ 4b1369f0-2a70-11ec-2cf3-e514a4899e06
begin

# Store mean and standard error of estimates for the exponential cost function specification. Our function returns an array containing 1000 vectors of 8 elements each. The estimate for k is in position 1 in each 1000 vectors etc

estimatesExp = mindisest.(E11,E12,E13,E31,E32,E10,E41,E42,"Exp")
k_exp_mean = mean([item[1] for item in estimatesExp])
k_exp_sd = std([item[1] for item in estimatesExp])
g_exp_mean = mean([item[2] for item in estimatesExp])
g_exp_sd = std([item[2] for item in estimatesExp])
s_exp_mean = mean([item[3] for item in estimatesExp])
s_exp_sd = std([item[3] for item in estimatesExp])
alpha_exp_mean = mean([item[4] for item in estimatesExp])
alpha_exp_sd = std([item[4] for item in estimatesExp])
a_exp_mean = mean([item[5] for item in estimatesExp])
a_exp_sd = std([item[5] for item in estimatesExp])
s_ge_exp_mean = mean([item[6] for item in estimatesExp])
s_ge_exp_sd = std([item[6] for item in estimatesExp])
beta_exp_mean = mean([item[7] for item in estimatesExp])
beta_exp_sd = std([item[7] for item in estimatesExp])
delta_exp_mean = mean([item[8] for item in estimatesExp])
delta_exp_sd = std([item[8] for item in estimatesExp])

end;

# ╔═╡ 2e15f490-2a6f-11ec-3bec-0f631d9fbd70
begin

# Store mean and standard error of estimates for the power cost function specification

estimatesPower = mindisest.(E11,E12,E13,E31,E32,E10,E41,E42,"Power")
k_power_mean = mean([item[1] for item in estimatesPower])
k_power_sd = std([item[1] for item in estimatesPower])
g_power_mean = mean([item[2] for item in estimatesPower])
g_power_sd = std([item[2] for item in estimatesPower])
s_power_mean = mean([item[3] for item in estimatesPower])
s_power_sd = std([item[3] for item in estimatesPower])

# A couple of elements in the arrays take value "NaN". We need to remove them before computing the mean and std

alpha_power = [item[4] for item in estimatesPower]
alpha_power_mean = mean(alpha_power[(!isnan).(alpha_power)])
alpha_power_sd = std(alpha_power[(!isnan).(alpha_power)])
a_power = [item[5] for item in estimatesPower]
a_power_mean = mean(a_power[(!isnan).(a_power)])
a_power_sd = std(a_power[(!isnan).(a_power)])
s_ge_power = [item[6] for item in estimatesPower]
s_ge_power_mean = mean(s_ge_power[(!isnan).(s_ge_power)])
s_ge_power_sd = std(s_ge_power[(!isnan).(s_ge_power)])
beta_power = [item[7] for item in estimatesPower]
beta_power_mean = mean(beta_power[(!isnan).(beta_power)])
beta_power_sd = std(beta_power[(!isnan).(beta_power)])
delta_power = [item[8] for item in estimatesPower]
delta_power_mean = mean(delta_power[(!isnan).(delta_power)])
delta_power_sd = std(delta_power[(!isnan).(delta_power)])

end;

# ╔═╡ c2d506c0-2a6f-11ec-0871-fdb55e56de67
begin

# To obtain confidence intervals in panel B table 5. CI are the 2.5% and 97.5% quantiles of the distribution of our parameters vectors. We obtain CI only for alpha, a, s_ge, beta, delta as in the paper

CI_Exp, CI_Power = [], []
for ci=4:8
    a  = sort!([item[ci] for item in estimatesExp]) # sort from smallest to biggest number
    append!(CI_Exp, [quantile!(a,0.025) ,quantile!(a,0.975)])
end
vecest = [alpha_power, a_power, s_ge_power, beta_power, delta_power]
for ci in vecest
    a  = sort!(ci[(!isnan).(ci)])
    append!(CI_Power, [quantile!(a,0.025) ,quantile!(a,0.975)])
end

end

# ╔═╡ a9cedae0-2a7c-11ec-3b2f-a3577f3decf7
begin

# Print our results for Table 5 (GMM estimates)

# Standard errors are different since the seed we used for the bootstrap procedure is different from the one used by the authors since random generation across softwares/languages is not easily replicated (each software uses its own algorithm)

sd_exp = [k_exp_sd,g_exp_sd,s_exp_sd,alpha_exp_sd,a_exp_sd,s_ge_exp_sd,beta_exp_sd,delta_exp_sd]

sd_power = [k_power_sd,g_power_sd,s_power_sd,alpha_power_sd,a_power_sd,s_ge_power_sd,beta_power_sd,delta_power_sd]

params_name = ["Level k of cost of effort", "Curvature γ of cost function","Intrinsic motivation s","Social preferences α", "Warm glow coefficient a","Gift exchange Δs", "Present bias β","(Weekly) discount factor δ"]

Table5Results = DataFrame()
Table5Results.params = params_name

# Table5Power is saved as a NTuple{8,Float64} but DataFrame only accepts vectors as columns. Using Collect transforms the NTuple into a vector

Table5Results.estimatesPower = collect(Table5Power)
Table5Results.sdPower = sd_power
Table5Results.estimatesExp = collect(Table5Exp)
Table5Results.sdExp = sd_exp
Table5Results

end

# ╔═╡ 88d81640-5609-11ec-1234-23840357ca64
# Save Table5Results as a csv file

CSV.write(raw"../output/table5GMM_julia.csv", Table5Results)

# ╔═╡ Cell order:
# ╠═af9b5960-2a65-11ec-3015-53185d92610c
# ╠═03372860-2a66-11ec-2b31-298b2c9afc9d
# ╠═7c88c030-55f8-11ec-0159-99a5b81bf850
# ╠═7dcdea60-55f8-11ec-2fdc-1d683f29a01c
# ╠═030f5512-2a66-11ec-1a7a-21059e6facac
# ╠═03008800-2a66-11ec-0c80-1767685eb4fa
# ╠═944bb54e-55f9-11ec-32b3-c953374ecd8f
# ╠═030039e0-2a66-11ec-3766-ad9593bffbd3
# ╠═ced0482e-2a66-11ec-06ca-0b961f35a057
# ╠═9ac87bc0-5608-11ec-3f39-3713d508fdfb
# ╠═4b1369f0-2a70-11ec-2cf3-e514a4899e06
# ╠═2e15f490-2a6f-11ec-3bec-0f631d9fbd70
# ╠═c2d506c0-2a6f-11ec-0871-fdb55e56de67
# ╠═a9cedae0-2a7c-11ec-3b2f-a3577f3decf7
# ╠═88d81640-5609-11ec-1234-23840357ca64


text:
## Dirichlet ##

struct TuringDirichlet{T<:Real,TV<:AbstractVector,S<:Real} <: ContinuousMultivariateDistribution
    alpha::TV
    alpha0::T
    lmnB::S
end

function TuringDirichlet(alpha::AbstractVector)
    all(ai -> ai > 0, alpha) ||
        throw(ArgumentError("Dirichlet: alpha must be a positive vector."))

    alpha0 = sum(alpha)
    lmnB = sum(loggamma, alpha) - loggamma(alpha0)

    return TuringDirichlet(alpha, alpha0, lmnB)
end
TuringDirichlet(d::Integer, alpha::Real) = TuringDirichlet(Fill(alpha, d))

# TODO: remove?
TuringDirichlet(alpha::AbstractVector{<:Integer}) = TuringDirichlet(float.(alpha))
TuringDirichlet(d::Integer, alpha::Integer) = TuringDirichlet(d, float(alpha))

# TODO: remove and use `Dirichlet` only for `Tracker.TrackedVector`
Distributions.Dirichlet(alpha::AbstractVector) = TuringDirichlet(alpha)

TuringDirichlet(d::Dirichlet) = TuringDirichlet(d.alpha, d.alpha0, d.lmnB)

Base.length(d::TuringDirichlet) = length(d.alpha)

# copied from Distributions
# TODO: remove and use `Dirichlet`?
function Distributions._rand!(
    rng::Random.AbstractRNG,
    d::TuringDirichlet,
    x::AbstractVector{<:Real},
)
    @inbounds for (i, αi) in zip(eachindex(x), d.alpha)
        x[i] = rand(rng, Gamma(αi))
    end
    Distributions.multiply!(x, inv(sum(x))) # this returns x
end
function Distributions._rand!(
    rng::AbstractRNG,
    d::TuringDirichlet{<:Real,<:FillArrays.AbstractFill},
    x::AbstractVector{<:Real}
)
    rand!(rng, Gamma(FillArrays.getindex_value(d.alpha)), x)
    Distributions.multiply!(x, inv(sum(x))) # this returns x
end

function Distributions._logpdf(d::TuringDirichlet, x::AbstractVector{<:Real})
    return simplex_logpdf(d.alpha, d.lmnB, x)
end
function Distributions.logpdf(d::TuringDirichlet, x::AbstractMatrix{<:Real})
    size(x, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    return simplex_logpdf(d.alpha, d.lmnB, x)
end

ZygoteRules.@adjoint function Distributions.Dirichlet(alpha)
    return ZygoteRules.pullback(TuringDirichlet, alpha)
end
ZygoteRules.@adjoint function Distributions.Dirichlet(d, alpha)
    return ZygoteRules.pullback(TuringDirichlet, d, alpha)
end

simplex_logpdf(alpha, lmnB, x::AbstractVector) = sum(xlogy.(alpha .- 1, x)) - lmnB
function simplex_logpdf(alpha, lmnB, x::AbstractMatrix)
    return vec(sum(xlogy.(alpha .- 1, x); dims=1)) .- lmnB
end

ZygoteRules.@adjoint function simplex_logpdf(alpha, lmnB, x::AbstractVector)
    simplex_logpdf(alpha, lmnB, x), Δ -> (Δ .* log.(x), -Δ, Δ .* (alpha .- 1) ./ x)
end

ZygoteRules.@adjoint function simplex_logpdf(alpha, lmnB, x::AbstractMatrix)
    simplex_logpdf(alpha, lmnB, x), Δ -> begin
        (log.(x) * Δ, -sum(Δ), ((alpha .- 1) ./ x) * Diagonal(Δ))
    end
end

## MvNormal ##

"""
    TuringDenseMvNormal{Tm<:AbstractVector, TC<:Cholesky} <: ContinuousMultivariateDistribution

A multivariate Normal distribution whose covariance is dense. Compatible with Tracker.
"""
struct TuringDenseMvNormal{Tm<:AbstractVector, TC<:Cholesky} <: ContinuousMultivariateDistribution
    m::Tm
    C::TC
end
function TuringDenseMvNormal(m::AbstractVector, A::AbstractMatrix)
    return TuringDenseMvNormal(m, cholesky(A))
end
Base.length(d::TuringDenseMvNormal) = length(d.m)
Distributions.rand(d::TuringDenseMvNormal, n::Int...) = rand(Random.GLOBAL_RNG, d, n...)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringDenseMvNormal, n::Int...)
    return d.m .+ d.C.U' * adapt_randn(rng, d.m, length(d), n...)
end

"""
    TuringDiagMvNormal{Tm<:AbstractVector, Tσ<:AbstractVector} <: ContinuousMultivariateDistribution

A multivariate normal distribution whose covariance is diagonal. Compatible with Tracker.
"""
struct TuringDiagMvNormal{Tm<:AbstractVector, Tσ<:AbstractVector} <: ContinuousMultivariateDistribution
    m::Tm
    σ::Tσ
end

Distributions.params(d::TuringDiagMvNormal) = (d.m, d.σ)
Base.length(d::TuringDiagMvNormal) = length(d.m)
Base.size(d::TuringDiagMvNormal) = (length(d),)
Distributions.rand(d::TuringDiagMvNormal, n::Int...) = rand(Random.GLOBAL_RNG, d, n...)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringDiagMvNormal, n::Int...)
    return d.m .+ d.σ .* adapt_randn(rng, d.m, length(d), n...)
end

struct TuringScalMvNormal{Tm<:AbstractVector, Tσ<:Real} <: ContinuousMultivariateDistribution
    m::Tm
    σ::Tσ
end

Distributions.params(d::TuringScalMvNormal) = (d.m, d.σ)
Base.length(d::TuringScalMvNormal) = length(d.m)
Base.size(d::TuringScalMvNormal) = (length(d),)
Distributions.rand(d::TuringScalMvNormal, n::Int...) = rand(Random.GLOBAL_RNG, d, n...)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringScalMvNormal, n::Int...)
    return d.m .+ d.σ .* adapt_randn(rng, d.m, length(d), n...)
end

function Distributions._logpdf(d::TuringScalMvNormal, x::AbstractVector)
    σ2 = abs2(d.σ)
    return -(length(x) * log(2π * σ2) + sum(abs2.(x .- d.m)) / σ2) / 2
end
function Distributions.logpdf(d::TuringScalMvNormal, x::AbstractMatrix{<:Real})
    size(x, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    return -(size(x, 1) * log(2π * abs2(d.σ)) .+ vec(sum(abs2.((x .- d.m) ./ d.σ), dims=1))) ./ 2
end
function Distributions.loglikelihood(d::TuringScalMvNormal, x::AbstractMatrix{<:Real})
    σ2 = abs2(d.σ)
    return -(length(x) * log(2π * σ2) + sum(abs2.(x .- d.m)) / σ2) / 2
end

function Distributions._logpdf(d::TuringDiagMvNormal, x::AbstractVector)
    return -(length(x) * log(2π) + 2 * sum(log.(d.σ)) + sum(abs2.((x .- d.m) ./ d.σ))) / 2
end
function Distributions.logpdf(d::TuringDiagMvNormal, x::AbstractMatrix{<:Real})
    size(x, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    return -((size(x, 1) * log(2π) + 2 * sum(log.(d.σ))) .+ vec(sum(abs2.((x .- d.m) ./ d.σ), dims=1))) ./ 2
end
function Distributions.loglikelihood(d::TuringDiagMvNormal, x::AbstractMatrix{<:Real})
    return -(length(x) * log(2π) + 2 * size(x, 2) * sum(log.(d.σ)) + sum(abs2.((x .- d.m) ./ d.σ))) / 2
end

function Distributions._logpdf(d::TuringDenseMvNormal, x::AbstractVector)
    return -(length(x) * log(2π) + logdet(d.C) + sum(abs2.(zygote_ldiv(d.C.U', x .- d.m)))) / 2
end
function Distributions.logpdf(d::TuringDenseMvNormal, x::AbstractMatrix{<:Real})
    size(x, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    return -((size(x, 1) * log(2π) + logdet(d.C)) .+ vec(sum(abs2.(zygote_ldiv(d.C.U', x .- d.m)), dims=1))) ./ 2
end
function Distributions.loglikelihood(d::TuringDenseMvNormal, x::AbstractMatrix{<:Real})
    return -(length(x) * log(2π) + size(x, 2) * logdet(d.C) + sum(abs2.(zygote_ldiv(d.C.U', x .- d.m)))) / 2
end

function StatsBase.entropy(d::TuringScalMvNormal)
    s = log(d.σ)
    return length(d) * ((1 + oftype(s, log2π)) / 2 + s)
end
function StatsBase.entropy(d::TuringDiagMvNormal)
    s = sum(log, d.σ)
    return length(d) * (1 + oftype(s, log2π)) / 2 + s
end
function StatsBase.entropy(d::TuringDenseMvNormal)
    s = logdet(d.C)
    return (length(d) * (1 + oftype(s, log2π)) + s) / 2
end

TuringMvNormal(d::Int, σ::Real) = TuringMvNormal(zeros(d), σ)
TuringMvNormal(m::AbstractVector{<:Real}, σ::Real) = TuringScalMvNormal(m, σ)
TuringMvNormal(σ::AbstractVector) = TuringMvNormal(zeros(length(σ)), σ)
TuringMvNormal(A::AbstractMatrix) = TuringMvNormal(zeros(size(A, 1)), A)
function TuringMvNormal(m::AbstractVector{<:Real}, σ::AbstractVector{<:Real})
    return TuringDiagMvNormal(m, σ)
end
function TuringMvNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}},
) where {T <: Real}
    return TuringMvNormal(m, sqrt.(D.diag))
end
function TuringMvNormal(m::AbstractVector{<:Real}, A::AbstractMatrix{<:Real})
    return TuringDenseMvNormal(m, A)
end
function TuringMvNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvNormal(m, sqrt(A.λ))
end

# Mean and covariance
Distributions.mean(d::TuringDiagMvNormal) = d.m
Distributions.var(d::TuringDiagMvNormal) = abs2.(d.σ)
Distributions.cov(d::TuringDiagMvNormal) = Diagonal(var(d))

Distributions.mean(d::TuringScalMvNormal) = d.m
Distributions.var(d::TuringScalMvNormal) = Fill(abs2(d.σ),length(d.m))
Distributions.cov(d::TuringScalMvNormal) = Diagonal(var(d))

Distributions.mean(d::TuringDenseMvNormal) = d.m
Distributions.var(d::TuringDenseMvNormal) = diag(cov(d))
Distributions.cov(d::TuringDenseMvNormal) = Matrix(d.C) # turns cholesky to matrix

## MvLogNormal ##

struct TuringMvLogNormal{TD} <: AbstractMvLogNormal
    normal::TD
end
MvLogNormal(d::TuringDenseMvNormal) = TuringMvLogNormal(d)
MvLogNormal(d::TuringDiagMvNormal) = TuringMvLogNormal(d)
MvLogNormal(d::TuringScalMvNormal) = TuringMvLogNormal(d)
Distributions.length(d::TuringMvLogNormal) = length(d.normal)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringMvLogNormal)
    return Distributions.exp!(rand(rng, d.normal))
end
function Distributions.rand(rng::Random.AbstractRNG, d::TuringMvLogNormal, n::Int)
    return Distributions.exp!(rand(rng, d.normal, n))
end

function Distributions._logpdf(d::TuringMvLogNormal, x::AbstractVector{T}) where {T<:Real}
    if insupport(d, x)
        logx = log.(x)
        return Distributions._logpdf(d.normal, logx) - sum(logx)
    else
        return -T(Inf)
    end
end
function Distributions.logpdf(d::TuringMvLogNormal, x::AbstractMatrix{<:Real})
    size(x, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    if all(i -> DistributionsAD.insupport(d, view(x, :, i)), axes(x, 2))
        logx = log.(x)
        return Distributions.logpdf(d.normal, logx) .- vec(sum(logx; dims = 1))
    else
        return [Distributions._logpdf(d, view(x, :, i)) for i in axes(x, 2)]
    end
end
function Distributions.loglikelihood(d::TuringMvLogNormal, x::AbstractMatrix{<:Real})
    if all(i -> DistributionsAD.insupport(d, view(x, :, i)), axes(x, 2))
        logx = log.(x)
        return loglikelihood(d.normal, logx) - sum(logx)
    else
        r = Distributions._logpdf(d.normal, view(x, :, 1))
        return oftype(r, -Inf)
    end
end

function MvLogNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return MvLogNormal(MvNormal(m, D))
end

## Zygote adjoint

ZygoteRules.@adjoint function Distributions.MvNormal(
    A::Union{AbstractVector{<:Real}, AbstractMatrix{<:Real}},
)
    return ZygoteRules.pullback(TuringMvNormal, A)
end
ZygoteRules.@adjoint function Distributions.MvNormal(
    m::AbstractVector{<:Real},
    A::Union{Real, UniformScaling, AbstractVecOrMat{<:Real}},
)
    return ZygoteRules.pullback(TuringMvNormal, m, A)
end
ZygoteRules.@adjoint function Distributions.MvNormal(
    d::Int,
    A::Real,
)
    value, back = ZygoteRules.pullback(A -> TuringMvNormal(d, A), A)
    return value, x -> (nothing, back(x)[1])
end
for (f, T) in (
    (:_logpdf, :AbstractVector),
    (:logpdf, :AbstractMatrix),
    (:loglikelihood, :AbstractMatrix),
)
    @eval begin
        ZygoteRules.@adjoint function Distributions.$f(
            d::MvNormal{<:Real,<:PDMats.ScalMat},
            x::$T{<:Real},
        )
            return ZygoteRules.pullback(d, x) do d, x
                Distributions.$f(TuringScalMvNormal(d.μ, sqrt(d.Σ.value)), x)
            end
        end
        ZygoteRules.@adjoint function Distributions.$f(
            d::MvNormal{<:Real,<:PDMats.PDiagMat},
            x::$T{<:Real},
        )
            return ZygoteRules.pullback(d, x) do d, x
                Distributions.$f(TuringDiagMvNormal(d.μ, sqrt.(d.Σ.diag)), x)
            end
        end
        ZygoteRules.@adjoint function Distributions.$f(
            d::MvNormal{<:Real,<:PDMats.PDMat},
            x::$T{<:Real},
        )
            return ZygoteRules.pullback(d, x) do d, x
                Distributions.$f(TuringDenseMvNormal(d.μ, d.Σ.chol), x)
            end
        end

        ZygoteRules.@adjoint function Distributions.$f(
            d::MvLogNormal{<:Real,<:PDMats.ScalMat},
            x::$T{<:Real},
        )
            return ZygoteRules.pullback(d, x) do d, x
                dist = TuringMvLogNormal(
                    TuringScalMvNormal(d.normal.μ, sqrt(d.normal.Σ.value)),
                )
                Distributions.$f(dist, x)
            end
        end
        ZygoteRules.@adjoint function Distributions.$f(
            d::MvLogNormal{<:Real,<:PDMats.PDiagMat},
            x::$T{<:Real},
        )
            return ZygoteRules.pullback(d, x) do d, x
                dist = TuringMvLogNormal(
                    TuringDiagMvNormal(d.normal.μ, sqrt.(d.normal.Σ.diag)),
                )
                Distributions.$f(dist, x)
            end
        end
        ZygoteRules.@adjoint function Distributions.$f(
            d::MvLogNormal{<:Real,<:PDMats.PDMat},
            x::$T{<:Real},
        )
            return ZygoteRules.pullback(d, x) do d, x
                dist = TuringMvLogNormal(TuringDenseMvNormal(d.normal.μ, d.normal.Σ.chol))
                Distributions.$f(dist, x)
            end
        end
    end
end


text:
using Test
import sbp

st = sbp.Staggered
sop = sbp.OP2019

nx = 10
ny = 20

desc = sop.load_description()
xp, xm, h, Dxp, Dxm, Hxp, Hxm, Pxp, Pxm = sop.build_operators(desc, nx)
yp, ym, h, Dyp, Dym, Hyp, Hym, Pyp, Pym = sop.build_operators(desc, ny)

@testset "Boundary matrix"  begin
        Bp = st.boundary_matrix_p(nx)
        @test size(Bp, 1) == nx
        @test size(Bp, 2) == nx + 1
        @test isapprox(sum(Bp), 0)

        Bp = st.boundary_matrix_p(nx, false)
        @test size(Bp, 1) == nx
        @test size(Bp, 2) == nx + 1
        @test isapprox(sum(Bp), 0)

        Bm = st.boundary_matrix_m(nx+1)
        @test size(Bm, 1) == nx + 1
        @test size(Bm, 2) == nx
        @test isapprox(sum(Bm), 0)
end

@testset "2D operators" begin
        Bxp = st.boundary_matrix_p(nx)
        Bxm = st.boundary_matrix_m(nx+1)
        Byp = st.boundary_matrix_p(ny)
        Bym = st.boundary_matrix_m(ny+1)
        pp = st.build_operators_2d(Dxp, Dyp, Pxp, Pyp, Hxp, Hyp, Bxp, Byp)
        @test pp.nx == nx
        @test pp.ny == ny
        @test size(pp.Dx, 1) == nx * ny
        @test size(pp.Dx, 2) == (nx + 1) * ny
        @test size(pp.Dy, 1) == nx * ny
        @test size(pp.Dy, 2) == nx * (ny + 1)

        builder = n -> sop.build_operators(desc, n)
        pp, mm, pm, mp = st.build_all_operators_2d(builder, nx, ny)
        @test pp.nx == nx
        @test pp.ny == ny

        @test mm.nx == nx + 1
        @test mm.ny == ny + 1

        @test pm.nx == nx
        @test pm.ny == ny + 1

        @test mp.nx == nx + 1
        @test mp.ny == ny

        @test size(pp.Dx) == size(pp.Px)
        @test size(pp.Dy) == size(pp.Py)
        @test size(pp.Bx) == size(pp.Bx)
        @test size(pp.By) == size(pp.By)

        @test size(mm.Dx) == size(mm.Px)
        @test size(mm.Dy) == size(mm.Py)
        @test size(mm.Bx) == size(mm.Bx)
        @test size(mm.By) == size(mm.By)

        @test size(pm.Dx) == size(pm.Px)
        @test size(pm.Dy) == size(pm.Py)
        @test size(pm.Bx) == size(pm.Bx)
        @test size(pm.By) == size(pm.By)

        @test size(mp.Dx) == size(mp.Px)
        @test size(mp.Dy) == size(mp.Py)
        @test size(mp.Bx) == size(mp.Bx)
        @test size(mp.By) == size(mp.By)

end

@testset "Covariant basis (cell-centers)" begin

        builder = n -> sop.build_operators(desc, n)
        pp, mm, pm, mp = st.build_all_operators_2d(builder, nx, ny)
        fx = ones(pp.nx*pp.ny)
        fy = ones(pp.nx*pp.ny)
        a = st.build_covariant_basis_mm(fx, fy, mm, pm)

        x_r1 = zeros(mm.nx * mm.ny)
        x_r2 = zeros(mm.nx * mm.ny)
        y_r1 = zeros(mm.nx * mm.ny)
        y_r2 = zeros(mm.nx * mm.ny)

        @test size(a.x_r1) == size(a.x_r2)
        @test size(a.x_r1) == size(a.y_r1)
        @test size(a.x_r1) == size(a.y_r2)

        @test isapprox(a.x_r1, x_r1, atol=1.0)
        @test isapprox(a.x_r2, x_r2, atol=1.0)
        @test isapprox(a.y_r1, y_r1, atol=1.0)
        @test isapprox(a.y_r2, y_r2, atol=1.0)
end

@testset "Covariant basis (v1)" begin

        builder = n -> sop.build_operators(desc, n)
        pp, mm, pm, mp = st.build_all_operators_2d(builder, nx, ny)
        fx = ones(pp.nx * pp.ny)
        fy = ones(pp.nx * pp.ny)
        a = st.build_covariant_basis_pm(fx, fy, mm, pm, mp)

        x_r1 = zeros(pm.nx * pm.ny)
        x_r2 = zeros(pm.nx * pm.ny)
        y_r1 = zeros(pm.nx * pm.ny)
        y_r2 = zeros(pm.nx * pm.ny)

        @test size(a.x_r1) == size(a.x_r2)
        @test size(a.x_r1) == size(a.y_r1)
        @test size(a.x_r1) == size(a.y_r2)

        @test isapprox(a.x_r1, x_r1, atol=1.0)
        @test isapprox(a.x_r2, x_r2, atol=1.0)
        @test isapprox(a.y_r1, y_r1, atol=1.0)
        @test isapprox(a.y_r2, y_r2, atol=1.0)
end

@testset "Covariant basis (v2)" begin

        builder = n -> sop.build_operators(desc, n)
        pp, mm, pm, mp = st.build_all_operators_2d(builder, nx, ny)
        fx = ones(pp.nx * pp.ny)
        fy = ones(pp.nx * pp.ny)
        a = st.build_covariant_basis_mp(fx, fy, mm, pm, mp)

        x_r1 = zeros(mp.nx * mp.ny)
        x_r2 = zeros(mp.nx * mp.ny)
        y_r1 = zeros(mp.nx * mp.ny)
        y_r2 = zeros(mp.nx * mp.ny)

        @test size(a.x_r1) == size(a.x_r2)
        @test size(a.x_r1) == size(a.y_r1)
        @test size(a.x_r1) == size(a.y_r2)

        @test isapprox(a.x_r1, x_r1, atol=1.0)
        @test isapprox(a.x_r2, x_r2, atol=1.0)
        @test isapprox(a.y_r1, y_r1, atol=1.0)
        @test isapprox(a.y_r2, y_r2, atol=1.0)
end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
include("InertiaReserve.jl")
include("CleanEnergyReserve.jl")
include("ThermalCleanEnergy.jl")

export get_active_power
export get_active_power_limits
export get_available
export get_base_power
export get_bus
export get_dynamic_injector
export get_ext
export get_fuel
export get_name
export get_operation_cost
export get_prime_mover
export get_ramp_limits
export get_rating
export get_reactive_power
export get_reactive_power_limits
export get_requirement
export get_services
export get_status
export get_time_at_status
export get_time_frame
export get_time_limits
export get_time_series_container
export set_active_power!
export set_active_power_limits!
export set_available!
export set_base_power!
export set_bus!
export set_dynamic_injector!
export set_ext!
export set_fuel!
export set_name!
export set_operation_cost!
export set_prime_mover!
export set_ramp_limits!
export set_rating!
export set_reactive_power!
export set_reactive_power_limits!
export set_requirement!
export set_services!
export set_status!
export set_time_at_status!
export set_time_frame!
export set_time_limits!
export set_time_series_container!

max_stars_repo_name: sourabhdalvi/PowerSystemExtensions

max_issues_repo_name: sourabhdalvi/PowerSystemExtensions

max_forks_repo_name: sourabhdalvi/PowerSystemExtensions


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d1e8bbf9-c687-42da-88e4-95d731fc85ad
md"# 让我们开始吧"

# ╔═╡ 48e4cc4c-eec8-4768-a7b4-339607b323bc

# ╔═╡ dfe13ef0-06de-11ec-36bb-e97e67535626
@bind x html"<input type=range>"

# ╔═╡ 014a19b3-c9e2-4e0d-a100-367c04a72bc6
x

# ╔═╡ Cell order:
# ╟─d1e8bbf9-c687-42da-88e4-95d731fc85ad
# ╠═48e4cc4c-eec8-4768-a7b4-339607b323bc
# ╟─dfe13ef0-06de-11ec-36bb-e97e67535626
# ╠═014a19b3-c9e2-4e0d-a100-367c04a72bc6


text:
using Documenter, PerronFrobenius
ENV["GKSwstype"] = "100"
push!(LOAD_PATH,"../src/")

PAGES = [
    "Overview" => "index.md",
    "Estimators" =>
        ["Overview" => "estimators/estimators.md",
        "Grid estimator" => "estimators/direct.md"
        ]
]

makedocs(
    modules = [PerronFrobenius],
    format = :markdown,
    sitename = "PerronFrobenius.jl",
    authors = "Kristian Agasøster Haaga",
    pages = PAGES
    # Use clean URLs, unless built as a "local" build
    #html_prettyurls = !("local" in ARGS),
    #html_canonical = "https://kahaaga.github.io/TimeseriesSurrogates.jl/latest/"
)

# deploydocs(
#     repo   = "github.com/kahaaga/PerronFrobenius.jl.git",
#     julia  = "0.6",
#     target = "build",
#     deps = nothing,
#     make = nothing,
#     osname = "linux"
# )


text:
using Geodesy
using Base.Test

# Construction

x, y = (rand(2) - .5) * 10_000

@test ENU(x, y) == ENU(x, y, 0.0)

@test LLA(x, y) == LLA(x, y, 0.0)

ECEF(x, y, 0.0)

# get* methods

ll = LL(x, y)
lla = LLA(x, y, rand())
@test LLA(getX(lla), getY(lla), getZ(lla)) == lla
@test getY(ll) == y
@test getX(ll) == x

enu = ENU(x, y, rand())
@test ENU(getX(enu), getY(enu), getZ(enu)) == enu


text:
export replace_blurbs

"""
    replace_blurbs(path::AbstractString)
    replace_blurbs(file::File)
    replace_blurbs(dir::Dir)

    Replace sections looking like this:

        > **INFO Foobar**
        >
        > main text

    With AsciiDoc Sections looking like this:

        [NOTE]
        ====
        main text
        ====
"""
replace_blurbs(path::AbstractString) = replace_blurbs(Path(path))

@enum State FindStart FindEnd

function replace_blurbs(file::File)
    state = FindStart
    lines = readlines(file.path)

    open(file.path, "w") do out
         i = 1
         while i <= length(lines)
            line = lines[i]

            if state == FindStart && !startswith(line, "> **")
                println(out, line)
            elseif state == FindStart
                state = FindEnd
                tag = strip(in("> *"), line)
                println(out, "[$tag]")
                println(out, "====")

                # skip next line if it is empty
                # Pandoc blurb needs empty line to look like header
                if isempty(lstrip(in("> "), lines[i+1]))
                    i += 1
                end
            elseif state == FindEnd && !startswith(line, '>')
                state = FindStart
                println(out, "====")
                println(out, line)
            elseif state == FindEnd
                println(out, lstrip(in("> "), line))
            else
                error("$state is not handled")
            end

            i += 1
        end
    end
end

replace_blurbs(dir::Dir) = visitdir(replace_blurbs, dir.path)

max_stars_repo_name: ordovician/AsciiDocTools

max_issues_repo_name: ordovician/AsciiDocTools

max_forks_repo_name: ordovician/AsciiDocTools


text:
# ---
# title: 239. Sliding Window Maximum
# id: problem239
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap, Sliding Window, Dequeue
# link: <https://leetcode.com/problems/sliding-window-maximum/description/>
# hidden: true
# ---
#
# You are given an array of integers `nums`, there is a sliding window of size
# `k` which is moving from the very left of the array to the very right. You can
# only see the `k` numbers in the window. Each time the sliding window moves
# right by one position.
#
# Return _the max sliding window_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,-1,-3,5,3,6,7], k = 3
#     Output: [3,3,5,5,6,7]
#     Explanation:
#     Window position                Max
#     ---------------               -----
#     [1  3  -1] -3  5  3  6  7       **3**
#      1 [3  -1  -3] 5  3  6  7       **3**
#      1  3 [-1  -3  5] 3  6  7      **5**
#      1  3  -1 [-3  5  3] 6  7       **5**
#      1  3  -1  -3 [5  3  6] 7       **6**
#      1  3  -1  -3  5 [3  6  7]      **7**
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1], k = 1
#     Output: [1]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,-1], k = 1
#     Output: [1,-1]
#
#
# **Example 4:**
#
#
#
#     Input: nums = [9,11], k = 2
#     Output: [11]
#
#
# **Example 5:**
#
#
#
#     Input: nums = [4,-2], k = 2
#     Output: [4]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `-104 <= nums[i] <= 104`
#   * `1 <= k <= nums.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using CartesianGeneticProgramming

include("functions.jl")
include("individual.jl")
include("mutation.jl")
include("crossover.jl")
include("evolution.jl")


text:
module RGFPrecompTest2
    using RuntimeGeneratedFunctions
    RuntimeGeneratedFunctions.init(@__MODULE__)

    y_in_RGFPrecompTest2 = 2

    # Simulates a helper function which generates an RGF, but caches it in a
    # different module.
    function generate_rgf(cache_module)
        context_module = @__MODULE__
        RuntimeGeneratedFunction(cache_module, @__MODULE__, :((x)->y_in_RGFPrecompTest2+x))
    end
end


text:
using Documenter, Example

makedocs(modules = [Example], sitename = "Example.jl")

deploydocs(repo = "github.com/mortenpi/Example.jl.git")


text:
module RowMajorArrays

export RowMajorArray

using LinearAlgebra

"""
Wrapper of a column major array (e.g. `Array` or `CuArray`) to make it a row-major array.

The default constructor takes a column major array as input, it is interpreted as row-major array.
This causes an implicit transpose of the input array, e.g. a `(2, 3)` column major array is converted to a
`(3, 2)` row major array.

For more than 2 dimensions, the order of the dimensions is reversed, e.g. a `(2, 3, 4, 5)` element column major array
is converted to a `(5, 4, 3, 2)` RowMajorArray.

Example:

    a = [1 2 3; 4 5 6] # 2x3 (column major) array
    a_r = RowMajorArray(a) # 3x2 row-major array

"""
struct RowMajorArray{T, N, A <: AbstractArray{T, N}} <: AbstractArray{T, N}
    data:: A
end

"""
Constructor which allows modification of array eltype.

Example:

    a = [1 2 3; 4 5 6] # eltype: Int64
    a_r = RowMajorArray{Float64}(a) # eltype Float64

"""
function RowMajorArray{T}(a:: AbstractArray{X, N}) where {X, T, N}
    a_converted = convert.(T, a)
    RowMajorArray{T, N, typeof(a_converted)}(a_converted)
end

# initializing with undefined fields
RowMajorArray{T, N, A}(::UndefInitializer, I::Vararg{<: Int, N}) where {T, N, A <: AbstractArray} = RowMajorArray(A(undef, reverse(I)...))

import Base

# getindex and setindex! methods
Base.getindex(a:: RowMajorArray, I::Vararg{<: Int, N}) where {N} = getindex(a.data, reverse(I)...)
Base.getindex(a:: RowMajorArray, i:: Int) = getindex(a.data, i) # linear indexing is done row-major here to have optimal memory alignment

Base.setindex!(a:: RowMajorArray, v, i:: Int) = setindex!(a.data, v, i)
Base.setindex!(a:: RowMajorArray, v, I::Vararg{<: Int, N}) where {N} = setindex!(a.data, v, reverse(I)...)

# standard function definitions
Base.size(a:: RowMajorArray) = reverse(size(a.data))
Base.size(a:: RowMajorArray, dim:: Int) = reverse(size(a.data))[dim]

Base.:(==)(a1:: RowMajorArray, a2:: RowMajorArray) = a1.data == a2.data

# for display for 2d RowMajorArrays, transpose data for printing so that the dimensions are printed in the right order
Base.show_nd(io::IO, a::RowMajorArray{T, 2, A}, print_matrix::Function, label_slices::Bool) where {T, A} = Base.show_nd(io, collect(transpose(a)), print_matrix, label_slices)

forward_methods = (:length, :eltype)
for m in forward_methods
    @eval Base.$m(a:: RowMajorArray) = Base.$m(a.data)
end

Base.similar(a:: RowMajorArray{T, N, A}) where {T, N, A} = similar(a, T, size(a))
Base.similar(a:: RowMajorArray{T, N, A}, dims:: Dims) where {T, N, A} =  RowMajorArray{T, N, A}(A(undef, reverse(dims)...))

function Base.similar(a:: RowMajorArray{X, N, A}, ::Type{T}, dims:: Dims) where {X, T, N, A}
    data = similar(a.data, T, reverse(dims))
    RowMajorArray(data)
end

# defining broadcasting
Base.BroadcastStyle(::Type{<: RowMajorArray}) = Broadcast.ArrayStyle{RowMajorArray}()

# see https://docs.julialang.org/en/v1/manual/interfaces/#Selecting-an-appropriate-output-array
function Base.similar(bc:: Broadcast.Broadcasted{Broadcast.ArrayStyle{RowMajorArray}}, ::Type{ElType}) where ElType
    # Scan the inputs for the ArrayAndChar:
    A = find_row_major_array(bc)
    # Use the char field of A to create the output
    similar(A, ElType, size(A))
end

# utility function for similar
# `A = find_row_major_array(As)` returns the first RowMajorArray among the arguments.
find_row_major_array(bc::Base.Broadcast.Broadcasted) = find_row_major_array(bc.args)
find_row_major_array(args::Tuple) = find_row_major_array(find_row_major_array(args[1]), Base.tail(args))
find_row_major_array(a::RowMajorArray) = a
find_row_major_array(::Tuple{}) = nothing
find_row_major_array(a::RowMajorArray, rest) = a
find_row_major_array(::Any, rest) = find_row_major_array(rest)

# constructor methods
# separate from `Base.zeros`, etc. - execute with qualified module name, e.g. `RowMajorArrays.zeros`, etc.
forward_methods = (:zeros, :ones, :rand, :randn)
for m in forward_methods
    @eval $m(T:: Type, I::Vararg{<: Int, N}) where {N} = RowMajorArray(Base.$m(T, reverse(I)...))
    @eval $m(I::Vararg{<: Int, N}) where {N} = RowMajorArray(Base.$m(reverse(I)...))
end

fill(x, I::Vararg{<: Int, N}) where {N} = RowMajorArray(Base.fill(x, reverse(I)...))
fill(x, I::Tuple) = RowMajorArray(Base.fill(x, reverse(I)))

# operations
forward_methods = (:(+), :(-))
for m in forward_methods
    @eval Base.$m(a1:: RowMajorArray, a2:: RowMajorArray) = RowMajorArray(Base.$m(a1.data, a2.data))
end

# Matrix multiplication - very crude implementation by converting to standard arrays.
# This should not be too slow for large matrices because converting is O(n²) whereas matrix multiplication is O(n³).
# Tdeally, this should be replaced with calls to OpenBlas Row-Major functions
Base.:(*)(a1:: RowMajorArray, a2:: RowMajorArray) = RowMajorArray(collect(transpose(a1.data)) * collect(transpose(a2.data)))

end


text:
struct SSARename
    vars::Dict{Symbol,Symbol}
end