text:
# make QuasiNewton{SR1, approx} callable and then pass that as the
# "update" parameter
# the minimize can take a single SR1 input as -> QuasiNewton{SR1, Approx}
struct SR1{T1} <: QuasiNewton{T1}
   approx::T1
end
SR1() = SR1(Direct())
hasprecon(::SR1) = NoPrecon()

summary(::SR1) = "SR1"
function update(scheme::SR1{<:Inverse}, H, s, y)
   T = real(eltype(s))
   w = s - H*y
   θ = dot(w, y) # angle between residual and change in gradient
   ρy = norm(y)
   if abs(θ) ≥ T(1e-12)*norm(w)*ρy && !iszero(ρy)
      H = H + (w*w')/θ
   end
   H
end
function update(scheme::SR1{<:Direct}, B, s, y)
   T = real(eltype(s))
   res = y - B*s # resesidual in secant equation
   θ = dot(res, s) # angle between residual and change in state
   if abs(inv(θ)) ≥ T(1e-12)*norm(res)*norm(s)
      if true #abs(θ) ≥ 1e-12
         B = B + (res*res')/θ
      end
   end
   B
end
function update!(scheme::SR1{<:Inverse}, H, s, y)
   T = real(eltype(s))
   w = s - H*y
   θ = dot(w, y) # angle between residual and change in gradient
   ρy = norm(y)
   if abs(θ) ≥ T(1e-12)*norm(w)*ρy && !iszero(ρy)
      H .= H .+ (w*w')/θ
   end
   H
end
function update!(scheme::SR1{<:Direct}, B, s, y)
   T = real(eltype(s))
   res = y - B*s
   θ = dot(res, s) # angle between residual and change in state
   if abs(inv(θ)) ≥ max(T(1e-12)*norm(res, 2)*norm(s, 2), sqrt(eps(T)))
      if true #abs(θ) ≥ 1e-12
         B .= B .+ (res*res')/θ
      end
   end
   B
end
function update!(scheme::SR1{<:Inverse}, A::UniformScaling, s, y)
   update(scheme, A, s, y)
end
function update!(scheme::SR1{<:Direct}, A::UniformScaling, s, y)
   update(scheme, A, s, y)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ResourceRequirements describes the compute resource requirements.

    IoK8sApiCoreV1ResourceRequirements(;
        limits=nothing,
        requests=nothing,
    )

    - limits::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Limits describes the maximum amount of compute resources allowed. More info: https://kubernetes.io/docs/concepts/configuration/manage-compute-resources-container/
    - requests::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Requests describes the minimum amount of compute resources required. If Requests is omitted for a container, it defaults to Limits if that is explicitly specified, otherwise to an implementation-defined value. More info: https://kubernetes.io/docs/concepts/configuration/manage-compute-resources-container/
"""
mutable struct IoK8sApiCoreV1ResourceRequirements <: SwaggerModel
    limits::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: limits
    requests::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: requests

    function IoK8sApiCoreV1ResourceRequirements(;limits=nothing, requests=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ResourceRequirements, Symbol("limits"), limits)
        setfield!(o, Symbol("limits"), limits)
        validate_property(IoK8sApiCoreV1ResourceRequirements, Symbol("requests"), requests)
        setfield!(o, Symbol("requests"), requests)
        o
    end
end # type IoK8sApiCoreV1ResourceRequirements

const _property_map_IoK8sApiCoreV1ResourceRequirements = Dict{Symbol,Symbol}(Symbol("limits")=>Symbol("limits"), Symbol("requests")=>Symbol("requests"))
const _property_types_IoK8sApiCoreV1ResourceRequirements = Dict{Symbol,String}(Symbol("limits")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("requests")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}")
Base.propertynames(::Type{ IoK8sApiCoreV1ResourceRequirements }) = collect(keys(_property_map_IoK8sApiCoreV1ResourceRequirements))
Swagger.property_type(::Type{ IoK8sApiCoreV1ResourceRequirements }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ResourceRequirements[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ResourceRequirements }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ResourceRequirements[property_name]

function check_required(o::IoK8sApiCoreV1ResourceRequirements)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ResourceRequirements }, name::Symbol, val)
end

ext: jl


text:
#Parses header files for the GLDefinitions

using Match, DictUtils

include("parseGLFile.jl")

function writeDictFunctionsGLU(fs, dict)
	write(fs, "#function bodies\n")
	for elem in dict
		write(fs, "@getCFun \"libGLU\" "*elem[1]* " " *elem[1]*elem[2]*"\n")
	end
end
function writeDictFunctions(fs, dict)
	write(fs, "#function bodies\n")
	for elem in dict
		write(fs, "@createOpenGLFun " * elem[1]*elem[2]*"\n")
	end
end
function writeDictConstants(fs, dict)
	write(fs, "#constants\n")
	for elem in dict
		spaces = 64 - length(elem[1])
		write(fs, "const " * elem[1] * " "^spaces *"= " *elem[2]*"\n")
	end
end
function writeDictTypes(fs, dict)
	write(fs, "#types\n")
	for elem in dict
		spaces = 60 - length(elem[1])
		write(fs,"typealias " * elem[1] * " "^spaces * elem[2]*"\n")
	end
end

sourceLocation = "deprecated/"
openGLFolder = "../"

GL_FILES =  ["gl10","gl21", "gl33", "gl42", "gl43", "glu", "wgl", "nv", "glx", "amd", "arb", "ati", "ext"]

types 		= Dict{ASCIIString, Dict{ASCIIString, ASCIIString}}()
functions 	= Dict{ASCIIString, Dict{ASCIIString, ASCIIString}}()
constants 	= Dict{ASCIIString, Dict{ASCIIString, ASCIIString}}()

for file in GL_FILES
	t, f, c = parseJLGLFile(sourceLocation*file*"/"*file*".jl")

	types[file] 	= t
	functions[file] = f
	constants[file] = c
	println("parsed: $(file)")
end

#######################################################################################
#gl10 and glu types are not translated to the platform independent julia C types
gl10_types 	= get(types, "gl10", error)
glu_types 	= get(types, "glu", error)
gl43_types 	= get(types, "gl43", error)

for elem in gl10_types
	gl10_types[elem[1]] = get(gl43_types, elem[1], elem[2])
end
for elem in glu_types
	glu_types[elem[1]] = get(gl43_types, elem[1], elem[2])
end

#We just create one file with all type definitions...
#Types need to be defined in the right order
mergedTypes = Dict{ASCIIString, ASCIIString}()
for elem in types
	merge!(mergedTypes, elem[2])
end
typeValues 	= collect(values(mergedTypes))
allKeys 	= collect(keys(mergedTypes))
orderedKeys = ASCIIString[]
index = 0
# find out if value is defined, if yes it can be pushed into the ordered key Array
while length(orderedKeys) < length(typeValues)
	for elem in mergedTypes
		a = parse(elem[2])
		if isa(a, Symbol)
			if  isdefined(a) && !in(elem[1], orderedKeys)
				eval(parse("typealias $(elem[1]) $(elem[2])"))
				push!(orderedKeys, elem[1])
			end
		elseif isa(a, Expr)
			if  isdefined(a.args[2]) && !in(elem[1], orderedKeys)
				eval(parse("typealias $(elem[1]) $(elem[2])"))
				push!(orderedKeys, elem[1])
			end
		end
	end
end

# now we can write the files
fs = open(openGLFolder*"types.jl", "w")
write(fs, "#all OpenGL Types\n")
for elem in orderedKeys
	spaces = 60 - length(elem)
	write(fs,"typealias " * elem * " "^spaces * get(mergedTypes, elem, error)*"\n")
end
exports		= join(orderedKeys, ", ")
write(fs, "export "*exports*"\n")
close(fs)

#Writing all the other files
for file in GL_FILES

	_types 		= get(types, 		file, error)
	_functions 	= get(functions, 	file, error)
	_constants 	= get(constants, 	file, error)
	exports		= join([collect(keys(_functions)), collect(keys(_constants))], ", ")

	fs = open(openGLFolder*"$(file).jl", "w")
	if file == "glu"
		write(fs, "import GetC.@getCFun\n")
		writeDictFunctionsGLU(fs, 	_functions)
	else
		writeDictFunctions(fs, 	_functions)
	end
	writeDictConstants(fs, 	_constants)
	write(fs, "# Export everything!\n")
	write(fs, "export "*exports*"\n")
	close(fs)
	println("written: $(openGLFolder)$(file).jl")
end

# Functions to determine what the differences between the OpenGL versions are...
# Not used yet, But i might want to use them to create an info File
#
function sameArguments(arguments1, arguments2)
	matches = true
	i = 1
	while matches && !done(arguments1, i)
		matches = arguments1[i] == arguments2[i]
		i += 1
	end
	matches
end
function reduceToCommonEntries(dicts)
	allTypesResult = deepcopy(dicts)
	map!(dict -> begin mapValues!(x -> "0", dict) ; dict end, allTypesResult)
	index = 1
	for dict in dicts
		allKeys = keys(dict)
		for _key in allKeys
			for i = index:length(dicts)
				dict2 = allTypesResult[i]
				if haskey(dict2, _key)

					if get(dict, _key, error) == get(dicts[i], _key, "error2")
						dict2[_key] = string(int(get(dict2, _key, "0")) + 1)
					elseif isa(parse(get(dict, _key, error)), Expr)
						expr1 = get(dict, _key, "")
						expr2 = get(dicts[i], _key, "")
						@assert ~isempty(expr1) && ~isempty(expr2)
						expr1 = parse(expr1)
						expr2 = parse(expr2)
						if expr1.head == :(::) && expr2.head == :(::) &&
							length(expr1.args[1].args) == length(expr2.args[1].args)
							arguments1 = expr1.args[1].args
							arguments2 = expr2.args[1].args
							while matches && !done(arguments1, i)
								matches = arguments1[i] == arguments2[i]
								i += 1
							end
							if sameArguments(expr1.args[1].args, expr2.args[1].args)
								dict2[_key] = string(int(get(dict2, _key, "0")) + 1)
							end
						end

					end
				end
			end
		end
		index += 1
	end
	commonEntries = Dict{ASCIIString, ASCIIString}()
	for dict in allTypesResult
		allKeys = keys(dict)
		for _key in allKeys
			if int(get(dict, _key, "0")) == length(dicts)
				delete!(dict, _key)
				commonEntries[_key] = get(dicts[1], _key, error)
			end
		end
	end
	commonEntries
end
function dictDiff(subtractDict, dicts)
	subtractSet = Set(keys(subtractDict))
	setArray = map(x-> Set(keys(x)), dicts)
	map(s -> setdiff(s, subtractSet), setArray)
end
function dictDiff(subtractDict::Dict{ASCIIString, ASCIIString}, dict::Dict{ASCIIString, ASCIIString})
	subtractSet = Set(collect(keys(subtractDict))...)
	dictSet 	= Set(collect(keys(dict))...)

	tmp = setdiff(dictSet, subtractSet)
	returnDict = Dict{ASCIIString, ASCIIString}()

	for elem in tmp
		returnDict[elem] = get(dict, elem, error)
	end
	returnDict
end


text:
# productmap.jl

"""
A product map is diagonal and acts on each of the components of x separately:
`y = f(x)` becomes `y_i = f_i(x_i)`.
"""
struct ProductMap{MAPS,S <: Tuple,T <: Tuple} <: AbstractMap{S,T}
    # maps has an indexable and iterable type, for example a tuple of maps
    maps    ::  MAPS
end

function ProductMap(maps...)
    MAPS = typeof(maps)
    S = typeof(map(x->zero(domaintype(x)), maps))
    T = typeof(map(x->zero(codomaintype(x)), maps))

    ProductMap{MAPS,S,T}(maps)
end

elements(pm::ProductMap) = pm.maps

(m::ProductMap)(x) = applymap(m, x)

applymap(pm::ProductMap, x) = map(applymap, elements(pm), x)

tensorproduct(map1::AbstractMap, map2::AbstractMap) = ProductMap(map1, map2)
tensorproduct(map1::ProductMap, map2::AbstractMap) = ProductMap(elements(map1)..., map2)
tensorproduct(map1::AbstractMap, map2::ProductMap) = ProductMap(map1, elements(map2)...)
tensorproduct(map1::ProductMap, map2::ProductMap) = ProductMap(elements(map1)..., elements(map2)...)

for op in (:inv, :left_inverse, :right_inverse, :jacobian)
    @eval $op(pm::ProductMap) = ProductMap(map($op, elements(pm))...)
end


text:
keys = open("public_keys.txt") do f
    parse.(Int, readlines(f))
end

function transform(subject_number::Int64, loop_size::Int64)::Int64
    val = 1
    for i in 1:loop_size
        val = (val * subject_number) % 20201227
    end

    return val
end

function crack(public_keys::Array{Int64, 1}, sn::Int64, ls_range::UnitRange{Int64})::Tuple{Int64, Int64}
    val = 1
    for ls in ls_range
        val = (val * sn) % 20201227
        if val == public_keys[1]
            # println("sn $(sn), ls $(ls) match first")
            val2 = 1
            for ls2 in ls_range
                val2 = (val2 * sn) % 20201227
                if val2 == public_keys[2] && ls != ls2
                    # println("sn $(sn), ls $(ls2) match second")

                    return ls, ls2
                end
            end
            # println("no matching ls2 found")
        end
    end

    return nothing, nothing
end

ls1, ls2 = crack(keys, 7, 1:100_000_000)

println(transform(keys[1], ls2))
# println(transform(keys[2], ls1))


text:
macro disjunction(args...)
    pos_args, kw_args, _ = Containers._extract_kw_args(args)

    #get kw_args
    reformulation = filter(i -> i.args[1] == :reformulation, kw_args)
    if !isempty(reformulation)
        reformulation = reformulation[1].args[2]
    else
        throw(UndefKeywordError(:reformulation))
    end
    M = filter(i -> i.args[1] == :M, kw_args)
    M = !isempty(M) ? esc(M[1].args[2]) : :(missing)
    eps = filter(i -> i.args[1] == :eps, kw_args)
    eps = !isempty(eps) ? esc(eps[1].args[2]) : :(1e-6)
    name = filter(i -> i.args[1] == :name, kw_args)
    name = !isempty(name) ? esc(name[1].args[2]) : :(missing)

    #get args
    m = esc(pos_args[1])
    disj = [esc(a) for a in pos_args[2:end]]

    #build disjunction
    :(add_disjunction($m, $(disj...), reformulation = $reformulation, M = $M, eps = $eps, name = $name))
end


text:
using ConstrainedDynamics

path = "examples/examples_files/doublependulum.urdf"
mech, shapes = Mechanism(path)


text:
@testset "Residue" begin

    @testset "Comparisons" begin

        @test GAP == GAP
        @test XAA == XAA
        @test GAP != XAA

        for i in 1:22, j in 1:22
            if j != i
                @test Residue(i) != Residue(j)
            else
                @test Residue(i) == Residue(j)
            end
        end
    end

    @testset "Convert" begin

        alphabet = "ARNDCQEGHILKMFPSTWYV"

        @test Int[ Residue(char) for char in alphabet] == Int[ i for i in 1:20 ]
        @test Int[ Residue(char) for char in lowercase(alphabet)] == Int[ 21 for i in 1:20 ]

        @testset "GAP" begin

            @test Int(GAP) == 21

            @test Char(GAP) == '-'

            @test GAP == Residue('-')
            @test GAP == Residue('*')
            @test GAP == Residue('.')

            for X in "UOBZJX"
                @test GAP != Residue(X)
                @test GAP == Residue(lowercase(X))
            end

            @test GAP != Residue(42)
            @test GAP != Residue(-1)
        end

        @testset "XAA" begin

            @test Int(XAA) == 22

            @test Char(XAA) == 'X'

            for X in "UOBZJX"
                @test XAA == Residue(X)
                @test XAA != Residue(lowercase(X))
            end

            @test XAA == Residue(42)
            @test XAA == Residue(-1)
        end
    end

    @testset "Valid Residues" begin

        for res in res"ARNDCQEGHILKMFPSTWYV-X"
            @test isvalid(res)
        end

        @test isvalid(Residue('ñ'))
        @test isvalid(Residue(42))
        @test !isvalid(reinterpret(Residue, 42))
    end

    @testset "Strings" begin

        alphabet = "ARNDCQEGHILKMFPSTWYV-X"

        @test res"ARNDCQEGHILKMFPSTWYV-X" == Residue[ char for char in alphabet ]
        @test String(res"ARNDCQEGHILKMFPSTWYV-X") == alphabet
    end

    @testset "String vector" begin

        msa =  ["DAWAEF",
                "DAWAED",
                "DAYCMD" ]

        badmsa=["DAWAEF",
                "DAWAED",
                "DAM"    ]

        @test convert(Matrix{Residue}, msa) == Residue['D' 'A' 'W' 'A' 'E' 'F'
                                                       'D' 'A' 'W' 'A' 'E' 'D'
                                                       'D' 'A' 'Y' 'C' 'M' 'D']

        @test_throws ErrorException convert(Matrix{Residue}, String[])
        @test_throws ErrorException convert(Matrix{Residue}, badmsa)
    end

    @testset "Random" begin

        for i in 1:5000
            res = rand(Residue)
            @test res != GAP
            @test res != XAA
            @test isvalid(res)
            @test typeof(res) == Residue
        end
    end

    @testset "Initialized arrays" begin

        @test size(rand(Residue, 20)) == (20,)
        @test typeof(rand(Residue, 20)) == Array{Residue, 1}
        @test size(rand(Residue, 20,30)) == (20,30)
        @test typeof(rand(Residue, 20, 30)) == Array{Residue, 2}

        @test zero(Residue) == GAP
        @test one( Residue) == XAA

        @test zeros(Residue,4) == [GAP, GAP, GAP, GAP]
        @test ones(Residue,4)  == [XAA, XAA, XAA, XAA]

        @test zeros(Residue,2,2) == [   GAP GAP
                                        GAP GAP ]

        @test ones(Residue,2,2)  == [   XAA XAA
                                        XAA XAA ]
    end

    @testset "Other base methods" begin

        for i in 1:22
            @test bitstring(i) == bitstring(Residue(i))
        end

        @test typemin(Residue) == Residue(1)
        @test typemax(Residue) == XAA
    end

    @testset "Show" begin

        io = IOBuffer()
        alphabet = "ARNDCQEGHILKMFPSTWYV-X"
        for char in alphabet
            show(io, Residue(char))
            @test String(take!(io)) == String([char])
        end

        show(io, reinterpret(Residue, -30))
        @test String(take!(io)) == "�"
    end

    @testset "Print" begin

        io = IOBuffer()
        alphabet = "ARNDCQEGHILKMFPSTWYV-X"
        for char in alphabet
            print(io, Residue(char))
            @test String(take!(io)) == String([char])
        end

        print(io, reinterpret(Residue, -30))
        @test String(take!(io)) == "X"
    end
end


text:
# Functions for reading, using, analysing gensel outputs (alphas)
# alphas come in binary format (m=Int32, m*{Int32, Float32})
function read_bin_alphas(alphas_filename, m)
    stream = open(alphas_filename, "r")
    alphas = Array{Float64}(1:m)
    while !eof(stream)
        mx = read(stream, Int32)
        alphas_col = zeros(Float64, m)
        for i in 1:mx
            index = read(stream, Int32) + 1
            alpha = read(stream, Float32)
            alphas_col[index] = alpha
        end
        alphas = hcat(alphas, alphas_col)
    end
    return alphas
end
# mrkRes files
function read_alphanames(alphanames_filename)
    alpharesults = readdlm(alphanames_filename, ' ', String, '\n', header = true)
    return alpharesults[1][:,1]
end
# keep every 'freq'th sample between 'first' and 'last' index.
# Note every 'freq'th means if modulo freq == 0. NOT from first.
function filter_alphas(alphas, freq, first, final)
    filterindex = Int64[]
    for i in first:final
        if (i % freq) == 0
            push!(filterindex, i)
        end
    end
    alphas.alphas = alphas.alphas[:,filterindex]
    alphas.k = size(alphas.alphas, 2)
end

function print_alphas(alphas, output_filename)
    output_alphas = hcat(alphas.markerNames, alphas.alphas)
    writedlm(output_filename, output_alphas)
end
# Make Alphas type
function make_alphas(alphas, alphanames_filename, alphas_filename)
    alphas.markerNames = read_alphanames(alphanames_filename)
    alphas.m = size(alphas.markerNames, 1)
    alphas.alphas = read_bin_alphas(alphas_filename, alphas.m)
    alphas.k = size(alphas.alphas, 2)
end
# geweke diagnostics of MCMC chains
function get_mcmc_diagnostics(alphas, outputprefix)
    geweke_pvalues = zeros(Float64, alphas.m)
    try
        geweke_pvalues = gewekediag(Chains(alphas.alphas[:,2:end]'))[1][:,3]
    catch
        for i in 1:alphas.m
            geweke_pvalues[i] = gewekediag(alphas.alphas[i,2:end])[2]
        end
    end
    geweke_by_marker = hcat(alphas.markerNames, geweke_pvalues)
    writedlm(outputprefix .* ".geweke.stats", geweke_by_marker)

    geweke_pass = sum(geweke_pvalues .> 0.05)
    writedlm(outputprefix .* "geweke.summary", hcat(alphas.m, geweke_pass, geweke_pass/alphas.m))
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function efc(d)
    if "EFC" in keys(d["ElectricStorage"])
        return d["ElectricStorage"]["EFC"]
    end
    # TODO use REoptInputs for initial SOC
    power = diff(append!([0.5], d["ElectricStorage"]["year_one_soc_series_pct"]))
    abs_power = abs.(power)
    efc = zeros(365)
    # TODO handle time resolution other than hourly
    for d in 1:365
        ts0 = Int(round((24 * (d - 1) + 1)))
        tsF = Int(round(24 * d))
        efc[d] = sum(abs_power[ts0:tsF])/2
    end
    return efc
end

function plot_violin_compare_soc_efc(d1::Dict, d2::Dict;
    title="Violin distribution of SOC and EFC",
    name_soc_1="SOC no degr.",
    name_soc_2="SOC with degr.",
    name_efc_1="EFC no degr.",
    name_efc_2="EFC with degr.",
    )

    soc1 = d1["ElectricStorage"]["year_one_soc_series_pct"]
    soc2 = d2["ElectricStorage"]["year_one_soc_series_pct"]

    efc1 = efc(d1)
    efc2 = efc(d2)

    traces = [
        PlotlyJS.violin(
            y=soc1,
            side="negative",
            x=repeat(["SOC"], length(soc1)),
            spanmode="hard",
            points=false,
            name=name_soc_1,
        ),
        PlotlyJS.violin(
            y=soc2,
            side="positive",
            x=repeat(["SOC"], length(soc2)),
            spanmode="hard",
            points=false,
            name=name_soc_2,
        ),
        PlotlyJS.violin(
            y=efc1,
            side="negative",
            x=repeat(["EFC"], length(efc1)),
            spanmode="hard",
            points=false,
            name=name_efc_1,
        ),
        PlotlyJS.violin(
            y=efc2,
            side="positive",
            x=repeat(["EFC"], length(efc2)),
            spanmode="hard",
            points=false,
            name=name_efc_2,
        ),
    ]

    layout = PlotlyJS.Layout(
        title_text = title,
        font_size=20,
        # tickfont_size=20
    )

    PlotlyJS.plot(traces, layout)

end


text:
# This file contains a script to run posterior checks, deprecated.

using Plots
using PyPlot
using StatPlots
using KernelDensity
using Distributions
using DataFrames
using StatsBase
# set correct folder

try
  cd("DWPSDE model")
catch
  warn("Already in the DWPSDE model folder")
end

# load all functions
cd("..")
include(pwd()*"\\select case\\selectcase.jl")
cd("DWPSDE model")
include("set_up.jl")

################################################################################
#   Set up models
################################################################################

# nbr parameters
set_nbr_params = 7

# nbr cores
nbr_of_cores = 4 # was 10

# length burn-in
burn_in = 1

# nbr iterations
nbr_iterations = 5000 # should be 20000

# length training data
length_training_data = 5000 # thid should ne 5000

# log-scale priors
log_scale_prior = false

# algorithm
mcmc_alg = "MCWM"  # set MCWM or PMCMC

# prob run MH update
beta_MH = 1 # should be 0.1

# load training data
load_tranining_data = true

# type of job
job = "new_data" # set work to simdata or new_data

# set jod dep. parameters
if job == "simdata"

	# nbr particels
	nbr_particels = 200

	# use simulated data
	sim_data = true # set to true to use sim data

	# data set
	data_set = "old" # was "old"

	# dt
	dt = 0.035 # new = 0.35 old = 0.035

	# dt_U
	dt_U = 1. # new = 1 old = 1

elseif job == "new_data"

	# nbr particels
	nbr_particels = 250

	# use simulated data
	sim_data = false # set to true to use sim data

	# data set
	data_set = "new" # was "old"

	# dt
	dt = 0.35 # new = 0.35 old = 0.035

	# dt_U
	dt_U = 1. # new = 1 old = 1

end

# load stored data
load_tranining_data = true

################################################################################
##                         training data                                      ##
################################################################################

# set parameters
burn_in_tranining = 10000 # this should be 2000 when estimating 2 parameters

if set_nbr_params == 2
	burn_in_tranining = 1000 # this should be 2000
	length_training_data = 2000
end

nbr_iterations_tranining = burn_in_tranining +length_training_data
nbr_particels_tranining = 25 # should be 200
nbr_of_cores_tranining = 2 # should be 8

################################################################################
##      create (A)DA problem                                     					    ##
################################################################################

problem = set_up_gp_problem(nbr_of_unknown_parameters=set_nbr_params, use_sim_data = sim_data, data_set = data_set)
problem.alg_param.R = nbr_iterations
problem.alg_param.N = nbr_particels
problem.alg_param.burn_in = 1
problem.alg_param.nbr_of_cores = nbr_of_cores
problem.alg_param.length_training_data = length_training_data

problem.adaptive_update =  AMUpdate_gen(eye(set_nbr_params), 1/sqrt(set_nbr_params), 0.2, 1, 0.8, 25) # was 0.3

problem.alg_param.dt = dt
problem.alg_param.dt_U = dt_U
problem.alg_param.alg = mcmc_alg
problem.alg_param.compare_GP_and_PF = false
problem.alg_param.noisy_est = false
problem.alg_param.pred_method = "sample"
problem.alg_param.length_training_data = length_training_data
problem.alg_param.beta_MH = beta_MH

################################################################################
##                         set model parameters  for training                               ##
################################################################################

problem_training = set_up_problem(nbr_of_unknown_parameters=set_nbr_params, use_sim_data = sim_data, data_set = data_set)
problem_training.alg_param.alg = mcmc_alg
problem_training.alg_param.R = nbr_iterations_tranining
problem_training.alg_param.N = nbr_particels_tranining
problem_training.alg_param.burn_in = burn_in_tranining
problem_training.alg_param.nbr_of_cores = nbr_of_cores_tranining
problem_training.alg_param.dt = dt
problem_training.alg_param.dt_U = dt_U
problem_training.alg_param.alg = mcmc_alg
problem_training.adaptive_update =  AMUpdate_gen(eye(set_nbr_params), 1/sqrt(set_nbr_params), 0.15, 1, 0.8, 25)

################################################################################
#    load results
################################################################################

load_data_from_files = true # load data from files or form some  workspace
#dagp = "_dagp" #  set to _dagp to load ER-GP file  o.w. use ""
dagp = true

#=

jobnames

simdata:
jobname_mcwm = "gp_training_7_par_lunarc_simdata_4_coressimdata" # jobname for mcwm
jobname_da = "_dagpest7simdatada_gp_mcmc"
jobname_ada = "_dagpest7simdataada_gp_mcmc_dt"

new_data:
jobname_mcwm = "gp_training_7_par_lunarc_new_data_4_coresnew_data" # jobname for mcwm
jobname_da = "_dagpest7new_datada_gp_mcmc"
jobname_ada = "_dagpest7new_dataada_gp_mcmc_dt"
=#

jobname = "_dagpest7new_dataada_gp_mcmc_dt" # set to jobname string

if load_data_from_files

    data_res = convert(Array,readtable("Results/output_res"*jobname*".csv"))

    M, N = size(data_res)

    data_param = convert(Array,readtable("Results/output_param"*jobname*".csv"))

    theta_true = data_param[1:N-2]
    burn_in = Int64(data_param[N-2+1])

    data_prior_dist = convert(Array,readtable("Results/output_prior_dist"*jobname*".csv"))

    data_prior_dist_type = convert(Array,readtable("Results/output_prior_dist_type"*jobname*".csv"))
    data_prior_dist_type = data_prior_dist_type[2]

    Z = convert(Array,readtable("Results/data_used"*jobname*".csv"))
    Z = Z[:,1]

else

    # this option should be used to load from stored .jld files

end

if dagp
  burn_in = 1
end

Theta = data_res[:,1:N-2]' # stor data in column-major order
Theta = Theta[:, burn_in:end]

################################################################################
##  Posterior predictiv distribtuion
################################################################################

# generate samples from the posterior pred distribtuion
include("run_pf_paralell.jl")

N = nbr_particels

N_sample_from_posterior = 100

posterior_dist = Categorical(1/size(Theta,2)*ones(size(Theta,2)))

posterior_pred_samples = zeros(length(Z), N_sample_from_posterior)

for i = 1:N_sample_from_posterior

  # sample from posterior
  idx = rand(posterior_dist)

  # set parameters for pf
  theta_star = Theta[:,idx]
  (Κ, Γ, A, A_sign, B,c,d,g,f,power1,power2,sigma) = set_parameters(theta_star, problem.model_param.theta_known,length(theta_star))
  A = A*A_sign
  # set value for constant b function
  b_const = sqrt(2.*sigma^2 / 2.)
  (subsample_interval_calc, nbr_x0_calc, nbr_x_calc,N_calc) = map(Float64, (problem_training.alg_param.subsample_int, problem_training.alg_param.nbr_x0, problem_training.alg_param.nbr_x,N))

  # generate data from the boostrap filter
  Z_star = gen_trejectory_using_pf(Z,theta_star,problem.model_param.theta_known, N, N_calc, problem_training.alg_param.dt, problem_training.alg_param.dt_U, problem_training.alg_param.nbr_x0, nbr_x0_calc, problem_training.alg_param.nbr_x, nbr_x_calc, problem_training.alg_param.subsample_int, subsample_interval_calc, false, true, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)

  # store generated data set
  posterior_pred_samples[:,i] = Z_star

end

# analysis

dist_post_pred = Categorical(1/N_sample_from_posterior*ones(N_sample_from_posterior))
samples_dist_post_pred = rand(dist_post_pred,8)

text_size = 25
label_size = 20

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.subplot(331)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[1]])
PyPlot.subplot(332)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[2]])
PyPlot.subplot(333)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[3]])
PyPlot.subplot(334)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[4]])
PyPlot.subplot(335)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[5]])
PyPlot.subplot(336)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[6]])
PyPlot.subplot(337)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[7]])
#PyPlot.xlabel("Index",fontsize=text_size)
PyPlot.subplot(338)
PyPlot.plot(posterior_pred_samples[:,samples_dist_post_pred[8]])
#PyPlot.xlabel("Index",fontsize=text_size)
PyPlot.subplot(339)
PyPlot.plot(Z, "k")
#PyPlot.xlabel("Index",fontsize=text_size)
ax[:tick_params]("both",labelsize = label_size)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.subplot(331)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[1]],50)
#PyPlot.ylabel("Freq",fontsize=text_size)
PyPlot.subplot(332)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[2]],50)
PyPlot.subplot(333)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[3]],50)
PyPlot.subplot(334)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[4]],50)
#PyPlot.ylabel("Freq",fontsize=text_size)
PyPlot.subplot(335)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[5]],50)
PyPlot.subplot(336)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[6]],50)
PyPlot.subplot(337)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[7]],50)
#PyPlot.ylabel("Freq",fontsize=text_size)
PyPlot.subplot(338)
PyPlot.plt[:hist](posterior_pred_samples[:,samples_dist_post_pred[8]],50)
PyPlot.subplot(339)
PyPlot.plt[:hist](Z,50, color = "k")
ax[:tick_params]("both",labelsize = label_size)

PyPlot.figure()
for i = 1:8
  h = kde(posterior_pred_samples[:,samples_dist_post_pred[8]])
  PyPlot.plot(h.x,h.density)
end
h = kde(Z)
PyPlot.plot(h.x,h.density)

PyPlot.figure()
h = PyPlot.plt[:hist](mean(posterior_pred_samples,1)[:],20)
PyPlot.plot((mean(Z), mean(Z)), (0, maximum(h[1])), "k")

PyPlot.figure()
h = PyPlot.plt[:hist](std(posterior_pred_samples,1)[:],20)
PyPlot.plot((std(Z), std(Z)), (0, maximum(h[1])), "k")

posterior_pred_auto_cov = autocov(posterior_pred_samples, 1:5)
data_auto_cov = autocov(Z, 1:5)

PyPlot.figure()
h = PyPlot.plt[:hist](posterior_pred_auto_cov[1,:],20)
PyPlot.plot((data_auto_cov[1], data_auto_cov[1]), (0, maximum(h[1])), "k")

PyPlot.figure()
h = PyPlot.plt[:hist](posterior_pred_auto_cov[2,:],20)
PyPlot.plot((data_auto_cov[2], data_auto_cov[2]), (0, maximum(h[1])), "k")

PyPlot.figure()
h = PyPlot.plt[:hist](posterior_pred_auto_cov[3,:],20)
PyPlot.plot((data_auto_cov[3], data_auto_cov[3]), (0, maximum(h[1])), "k")

PyPlot.figure()
h = PyPlot.plt[:hist](posterior_pred_auto_cov[4,:],20)
PyPlot.plot((data_auto_cov[4], data_auto_cov[4]), (0, maximum(h[1])), "k")

PyPlot.figure()
h = PyPlot.plt[:hist](posterior_pred_auto_cov[5,:],20)
PyPlot.plot((data_auto_cov[5], data_auto_cov[5]), (0, maximum(h[1])), "k")

# calc hist for nbr break points

# print path for real data
print("Results/data_used"*dagp*jobname*".csv")

# write data from posterior pred dist
writetable("./Results/post_pred_data.csv", convert(DataFrame, posterior_pred_samples))

nbr_cp_real_data = length(readdlm("Results/cp_data_real_data.txt")[2:end])

nbr_cp_post_pred = readdlm("Results/cp_data_post_pred.txt")[2:end]

# all 15 for mcwm

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
########################### Chisq Diagnostic ###########################

function update_hangartner_temp_vars!(u::Matrix{Int64}, X::Matrix{Int64}, t::Int64)
    d = size(X, 2)

    for j = 1:d
        u[X[t, j], j] += 1
    end
end

function hangartner_inner(Y::AbstractMatrix, m::Int)
    ## setup temp vars
    n, d = size(Y)

    # Count for each category in each chain
    u = zeros(Int64, m, d)
    v = zeros(Int64, m, d)

    for t = 1:n
        # fill out temp vars
        update_hangartner_temp_vars!(u, Y, t)
    end
    phia, chi_stat, m_tot = weiss_sub(u, v, n)

    return (n * sum(chi_stat), m_tot)
end

function weiss(X::AbstractMatrix{U}) where {U<:Any}
    ## number of iterations, number of chains
    n, d = size(X)

    ## mapping of values to integers
    v_dict = Dict{U,Int}()

    ## max unique categories
    mc = map(c -> length(unique(X[:, c])), 1:d)
    m = length(unique(X))

    ## counter of number of unique values in each chain
    r0 = 0

    ## Count for each category in each chain
    u = zeros(Int, m, d)

    ## Number of times a category did not transition in each chain
    v = zeros(Int, m, d)

    for t = 1:n
        for c = 1:d
            if !(X[t, c] in keys(v_dict))
                r0 += 1
                v_dict[X[t, c]] = r0
            end
            idx1 = v_dict[X[t, c]]
            u[idx1, c] += 1

            if t > 1
                if X[t-1, c] == X[t, c]
                    v[idx1, c] += 1
                end
            end
        end
    end
    phia, chi_stat, m_tot = weiss_sub(u, v, n)
    ca = (1 + phia) / (1 - phia)
    stat = (n / ca) * sum(chi_stat)
    pval = NaN
    if ((m_tot - 1) * (d - 1)) >= 1
        pval = 1 - cdf(Chisq((m_tot - 1) * (d - 1)), stat)
    end

    return (stat, m_tot, pval, ca)
end

function weiss_sub(u::Matrix{Int64}, v::Matrix{Int64}, t::Int)
    m, d = size(u)
    nt = 0.0
    dt = 0.0
    m_tot = 0

    mp = zeros(Float64, m, d)
    ma = zeros(Float64, m)
    phia = 0.0
    ca = 0.0
    df = 0.0

    chi_stat = zeros(Float64, d)

    for j = 1:m
        p1 = 0.0
        p2 = 0.0
        for l = 1:d
            #aggregate
            p1 += v[j, l] / (d * (t - 1))
            p2 += u[j, l] / (d * t)

            #per chain
            mp[j, l] = u[j, l] / t
            ma[j] += u[j, l] / (d * t)
        end
        nt += p1
        dt += p2^2

        if ma[j] > 0
            m_tot += 1
            for l = 1:d
                chi_stat[l] += (mp[j, l] - ma[j])^2 / ma[j]
            end
        end
    end
    phia = 1.0 + (1.0 / t) - ((1 - nt) / (1 - dt))
    phia = min(max(phia, 0.0), 1.0 - eps())
    return (phia, chi_stat, m_tot)
end

function update_billingsley_temp_vars!(f::Array{Int64,3}, X::Matrix{Int64}, t::Int64)
    d = size(X, 2)
    for j = 1:d
        if t > 1
            f[X[t-1, j], X[t, j], j] += 1
        end
    end
end

function billingsley_sub(f::Array{Int64,3})
    df = 0.0
    stat = 0.0

    m, d = size(f, 2), size(f, 3)

    # marginal transitions, i.e.
    # number of transitions from each category
    mf = mapslices(sum, f, dims = 2)

    # For each category, number of chains for which
    # that category was present
    A = vec(mapslices((x) -> sum(x .> 0), mf, dims = 3))

    # For each category, number of categories it
    # transitioned to
    B = vec(mapslices((x) -> sum(x .> 0), mapslices(sum, f, dims = 3), dims = 2))

    # transition probabilities in each chain
    P = f ./ mf

    # transition probabilities
    mP = reshape((mapslices(sum, f, dims = 3) ./ mapslices(sum, mf, dims = 3)), Val(2))

    idx = findall(x -> x != 0, A .* B)
    for j in idx
        #df for billingsley
        df += (A[j] - 1) * (B[j] - 1)

        #billingsley
        for k in idx
            if (mP[j, k] > 0.0)
                for l = 1:d
                    if mf[j, 1, l] > 0 && isfinite(P[j, k, l])
                        stat += mf[j, 1, l] * (P[j, k, l] - mP[j, k])^2 / mP[j, k]
                    end
                end
            end
        end
    end
    return (stat, df, mP)
end

function bd_inner(Y::AbstractMatrix, m::Int)
    num_iters, num_chains = size(Y)
    # Transition matrix for each chain
    f = zeros(Int64, m, m, num_chains)

    for t = 1:num_iters
        # fill out temp vars
        update_billingsley_temp_vars!(f, Y, t)
    end
    billingsley_sub(f)
end

function simulate_NDARMA(
    N::Int,
    p::Int,
    q::Int,
    prob::Vector{Float64},
    phi::Vector{Float64},
)
    X = zeros(Int64, N)
    X[1:p] = rand(Categorical(prob), p)
    d1 = Multinomial(1, phi)
    d2 = Categorical(prob)
    for t = (p+1):N
        alphabeta = rand(d1)
        eps = rand(d2, q + 1)
        X[t] = sum([X[(t-p):(t-1)]; eps] .* alphabeta)
    end
    return X
end

function simulate_MC(N::Int, P::Matrix{Float64})
    X = zeros(Int64, N)
    n, m = size(P)
    X[1] = sample(1:n)
    for i = 2:N
        X[i] = wsample(1:n, vec(P[X[i-1], :]))
    end
    return X
end

function diag_all(
    X::AbstractMatrix{U},
    method::Symbol,
    nsim::Int,
    start_iter::Int,
    step_size::Int,
) where {U<:Any}

    ## number of iterations, number of chains
    n, d = size(X)

    ## mapping of values to integers
    v_dict = Dict{U,Int}()

    ## max unique categories
    mc = map(c -> length(unique(X[:, c])), 1:d)
    m = length(unique(X))

    ## counter of number of unique values in each chain
    r0 = 0

    ## Count for each category in each chain
    u = zeros(Int, m, d)

    ## Number of times a category did not transition in each chain
    v = zeros(Int, m, d)

    ## transition matrix for each chain
    f = zeros(Int, m, m, d)

    result = zeros(Float64, 3, length(start_iter:step_size:n))
    result_iter = 1
    for t = 1:n
        for c = 1:d
            if !(X[t, c] in keys(v_dict))
                r0 += 1
                v_dict[X[t, c]] = r0
            end
            idx1 = v_dict[X[t, c]]
            u[idx1, c] += 1

            if t > 1
                idx2 = v_dict[X[t-1, c]]
                f[idx1, idx2, c] += 1

                if X[t-1, c] == X[t, c]
                    v[idx1, c] += 1
                end
            end
        end

        if ((t >= start_iter) && (((t - start_iter) % step_size) == 0))
            phia, chi_stat, m_tot = weiss_sub(u, v, t)
            hot_stat, df, mP = billingsley_sub(f)

            phat = mapslices(sum, u, dims = 2)[:, 1] / sum(mapslices(sum, u, dims = 2))
            ca = (1 + phia) / (1 - phia)
            stat = NaN
            pval = NaN
            df0 = NaN

            if method == :hangartner
                stat = t * sum(chi_stat)
                df0 = (m - 1) * (d - 1)
                if m > 1 && !isnan(stat)
                    pval = 1 - cdf(Chisq((m - 1) * (d - 1)), stat)
                end
            elseif method == :weiss
                stat = (t / ca) * sum(chi_stat)
                df0 = (m - 1) * (d - 1)
                pval = NaN
                if m > 1 && !isnan(stat)
                    pval = 1 - cdf(Chisq((m - 1) * (d - 1)), stat)
                end
            elseif method == :DARBOOT
                stat = t * sum(chi_stat)
                bstats = zeros(Float64, nsim)
                for b = 1:nsim
                    Y = hcat(
                        [simulate_NDARMA(t, 1, 0, phat, [phia, 1 - phia]) for j = 1:d]...,
                    )
                    s = hangartner_inner(Y, m)[1]
                    bstats[b] = s
                end
                idx = findall(!isnan(bstats))
                df0 = mean(bstats[idx])
                pval = sum(stat .<= bstats[idx]) / length(idx)
            elseif method == :MCBOOT
                bstats = zeros(Float64, nsim)
                for b = 1:nsim
                    Y = hcat([simulate_MC(t, mP) for j = 1:d]...)
                    s = hangartner_inner(Y, m)[1]
                    bstats[b] = s
                end
                idx = findall(!isnan(bstats))
                df0 = mean(bstats[idx])
                pval = sum(stat .<= bstats[idx]) / length(idx)
            elseif method == :billingsley
                stat = hot_stat
                df0 = df
                if df > 0 && !isnan(hot_stat)
                    pval = 1 - cdf(Chisq(df), hot_stat)
                end
            elseif method == :billingsleyBOOT
                stat = hot_stat
                bstats = zeros(Float64, nsim)
                for b = 1:nsim
                    Y = hcat([simulate_MC(t, mP) for j = 1:d]...)
                    (s, sd) = bd_inner(Y, m)[1:2]
                    bstats[b] = s / sd
                end
                idx = findall(!isnan(bstats))
                df0 = mean(bstats[idx])
                pval = sum(stat / df .<= bstats[idx]) / length(idx)
            else
                error("Unexpected")
            end
            result[:, result_iter] = [stat, df0, pval]
            result_iter += 1
        end
    end
    return result
end

function discretediag_sub(
    c::AbstractChains,
    frac::Real,
    method::Symbol,
    nsim::Int,
    start_iter::Int,
    step_size::Int,
)

    num_iters, num_vars, num_chains = size(c.value)

    vals = zeros(Float64, 3 * (num_chains + 1), num_vars)
    plot_vals_stat = zeros(length(start_iter:step_size:num_iters), num_vars)
    plot_vals_pval = zeros(length(start_iter:step_size:num_iters), num_vars)

    ## Between-chain diagnostic
    X = zeros(Int64, num_iters, num_chains)
    for j = 1:length(num_vars)
        X = convert(Array{Int64,2}, c.value[:, j, :])
        result = diag_all(X, method, nsim, start_iter, step_size)
        plot_vals_stat[:, j] = result[1, :] ./ result[2, :]
        plot_vals_pval[:, j] = result[3, :]
        vals[1:3, j] = result[:, end]
    end

    ## Within-chain diagnostic
    x = zeros(Int64, num_iters)
    Y = zeros(Int64, num_iters, 2)
    for j = 1:num_vars
        for k = 1:num_chains
            x = convert(Array{Int64,1}, c.value[:, j, k])

            idx1 = 1:round(Int64, frac * num_iters)
            idx2 = round(Int64, num_iters - frac * num_iters + 1):num_iters
            x1 = x[idx1]
            x2 = x[idx2]
            n_min = min(length(x1), length(x2))
            Y = [x1[1:n_min] x2[(end-n_min+1):end]]

            vals[(3+3*(k-1)+1):(3+3*(k-1)+3), j] =
                diag_all(Y, method, nsim, n_min, step_size)[:, end]
        end
    end
    return (collect(1:num_vars), vals, plot_vals_stat, plot_vals_pval)

end

function discretediagplot(
    c::AbstractChains;
    frac::Real = 0.3,
    method::Symbol = :weiss,
    nsim::Int = 1000,
    start_iter::Int = 100,
    step_size::Int = 10000,
)

    num_iters, num_vars, num_chains = size(c.value)

    valid_methods = [:hangartner, :weiss, :DARBOOT, :MCBOOT, :billingsley, :billingsleyBOOT]
    if !(method in valid_methods)
        methods_str = join([":$f" for f in valid_methods], ", ")
        throw(ArgumentError("method must be one of ", methods_str))
    end

    if !(0.0 < frac < 1.0)
        throw(ArgumentError("frac must be in (0,1)"))
    end

    if (start_iter > num_iters) || (step_size > num_iters)
        throw(ArgumentError("start_iter, step_size must be less than $num_iters"))
    end

    V, vals, plot_vals_stat, plot_vals_pval =
        discretediag_sub(c, frac, method, nsim, start_iter, step_size)

    p1 = Plots.plot(
        repeat(
            collect(c.range[start_iter:step_size:num_iters]) / 1000,
            outer = [length(V)],
        ),
        vcat([plot_vals_stat[:, j] for j = 1:length(V)]...),
        seriestype = :line,
        group = repeat(c.names[V], inner = [length(start_iter:step_size:num_iters)]),
        xlabel = "Iteration (thousands)",
        ylabel = "stat/df",
        legendtitle = "Variable",
    )

    p2 = Plots.plot(
        repeat(
            collect(c.range[start_iter:step_size:num_iters]) / 1000,
            outer = [length(V)],
        ),
        vcat([plot_vals_pval[:, j] for j = 1:length(V)]...),
        seriestype = :line,
        group = repeat(c.names[V], inner = [length(start_iter:step_size:num_iters)]),
        xlabel = "Iteration (thousands)",
        ylabel = "pval",
        legendtitle = "Variable",
    )

    return [p1, p2]
end
"""
    discretediag(c::AbstractChains; frac::Real=0.3,
                  method::Symbol=:weiss, nsim::Int=1000)

Compute the convergence diagnostic for a discrete variable. Several options are available by choosing method to be one of :hangartner, :weiss, :DARBOOT, MCBOOT, :billinsgley, :billingsleyBOOT. The first four are based off of Pearson’s chi-square test of homogeneity. The diagnostic tests whether the proportion of the categories of the discrete variable are similar in each chain. The last two methods test whether the transition probabilities between each category are similar between each chain. Along with a between chain assessment of convergence, a within-chain assessment is carried out by comparing a specified fraction (frac), or window, of the beginning of a chain to the specified fraction of the end of the chain. For within-chain assessment, users should ensure that there is sufficient separation between the windows to assume that their samples are independent. A non-significant test p-value indicates convergence. Significant p-values indicate non-convergence and the possible need to discard initial samples as a burn-in sequence or to simulate additional samples.

A `ChainSummary` type object with parameters contained in the rows of the `value` field. The first three columns correspond to the test statistic, degrees of freedom, and p-value of the between-chain assessment. The next columns are the test statistic, degrees of freedom, and p-value for each chain of the within-chain assessment.

* `c` : sampler output on which to perform calculations.

* `frac` : proportion of iterations to include in the first window.

* `method` : Specify which method to use. One of :hangartner, :weiss, :DARBOOT, MCBOOT, :billinsgley, :billingsleyBOOT`.

* `nsim` : For the bootstrap methods (:DARBOOT, :MCBOOT, and :billingsleyBOOT) the number of bootstrap simulations.
"""
function discretediag(
    c::AbstractChains;
    frac::Real = 0.3,
    method::Symbol = :weiss,
    nsim::Int = 1000,
)

    num_iters, num_vars, num_chains = size(c.value)

    valid_methods = [:hangartner, :weiss, :DARBOOT, :MCBOOT, :billingsley, :billingsleyBOOT]
    if !(method in valid_methods)
        methods_str = join([":$f" for f in valid_methods], ", ")
        throw(ArgumentError("method must be one of ", methods_str))
    end

    if !(0.0 < frac < 1.0)
        throw(ArgumentError("frac must be in (0,1)"))
    end

    V, vals =
        discretediag_sub(c, frac, method, nsim, size(c.value, 1), size(c.value, 1))[1:2]

    hdr = header(c) * "\nChisq Diagnostic:\nEnd Fractions = $frac\n" * "method = $method\n"
    ChainSummary(
        copy(round.(vals, digits = 3)'),
        c.names[V],
        convert(
            Array{AbstractString,1},
            vcat([["stat", "df", "p-value"] for k = 1:(num_chains+1)]...),
        ),
        hdr,
    )
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
using FileIO, ImageIO, Images

root_dir = dirname(Base.current_project())

epochs = load(joinpath(root_dir, "plots/training/epoch1.png"))
base_shape = (size(epochs,1),size(epochs,2),1)
let epochs = reshape(epochs, base_shape)
for i in 2:42
    epochs = cat(epochs, reshape(load(joinpath(root_dir, "plots/training/epoch$i.png")),base_shape), dims=3)
end

# Save animation as GIF
save(joinpath(root_dir, "plots/training/UDE_training.gif"), epochs, fps=5)
end # let


text:
begin
pushfirst!(LOAD_PATH, "@stdlib")
import Pkg
popfirst!(LOAD_PATH)

# We need to Pkg.instantiate the package environment that this notebook worker process will launch in
local my_dir = @__DIR__
local pluto_dir = joinpath(my_dir, "..", "..")

local runner_env_dir = mkpath(joinpath(Pkg.envdir(Pkg.depots()[1]), "__pluto_boot_" * string(VERSION)))
local new_ptoml_path = joinpath(runner_env_dir, "Project.toml")

local ptoml_contents = read(joinpath(pluto_dir, "Project.toml"), String)
write(new_ptoml_path, ptoml_contents[findfirst("[deps]", ptoml_contents)[1]:end])

local pkg_ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(new_ptoml_path))

try
    Pkg.resolve(pkg_ctx; io=devnull) # supress IO
catch
    # if it failed, do it again without suppressing io
    try
        Pkg.resolve(pkg_ctx)
    catch e
        @error "Failed to resolve notebook boot environment" exception=(e, catch_backtrace())
    end
end
try
    # we don't suppress IO for this one because it can take very long, and that would be a frustrating experience without IO
    if VERSION >= v"1.6.0-a"
        # precompilation switched off because of https://github.com/fonsp/Pluto.jl/issues/875
        Pkg.instantiate(pkg_ctx; update_registry=false, allow_autoprecomp=false)
    elseif VERSION >= v"1.3.0"
        # registry update is not required here and may save some time for startup
        Pkg.instantiate(pkg_ctx; update_registry=false)
    else
        Pkg.instantiate(pkg_ctx)
    end
catch e
    @error "Failed to instantiate notebook boot environment" exception=(e, catch_backtrace())
end

pushfirst!(LOAD_PATH, runner_env_dir)

#
include(joinpath(my_dir, "PlutoRunner.jl"))
#

popfirst!(LOAD_PATH)
end


text:
##COMPARE TIME TAKEN BY ADAPTIVE AND NON-ADAPTIVE RE-ABC IN NORMAL EXAMPLE
using RareEventABC
using Distributions
using DataFrames
using ProgressMeter
using Plots; pyplot();
using StatPlots;

"
A pseudo-marginal Metropolis-Hastings algorithm which uses RE-ABC likelihood estimates.

## Arguments

- `a` Details of the application.
- `schedule` Sequence of epsilon values. Use this option for non-adaptive RE-ABC.
- `targetε` Target ε value. Use this option for adaptive RE-ABC.
- `N` Number of iterations to do.
- `nparticles` Number of particles to use.
- `prior` Prior distribution.
- `θinit` Initial state of chain.
- `inc` Distribution of proposal increments
- `k` Number of particles to accept at each iteration of adaptive RE-ABC.
"

function timedPMMH(a::RareEventABC.ApplicationDetails, schedule::Array{Float64, 1}, N::Int, nparticles::Int, prior::UnivariateDistribution, θinit::Float64, inc::UnivariateDistribution)
    states = zeros(Float64, N)
    loglikelihoods = zeros(Float64, N)
    states[1] = θ = θinit
    times = DataFrame(time=Float64[], early=Bool[])
    tic()
    loglikelihoods[1] = ll = doSMC(a, [θ], schedule, nparticles=nparticles, quiet=true).llest
    t = toq()
    push!(times, (t, false))
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            tic()
            p = doSMC(a, [θprop], schedule, nparticles=nparticles, llBound=llBound, quiet=true)
            t = toq()
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
                push!(times, (t, false))
            else
                push!(times, (t, true))
            end
        end
        states[i], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods, times
end

function timedPMMH(a::RareEventABC.ApplicationDetails, targetε::Float64, N::Int, nparticles::Int, prior::UnivariateDistribution, θinit::Float64, inc::UnivariateDistribution; k=ceil(Int, nparticles/2))
    states = zeros(Float64, N)
    loglikelihoods = zeros(Float64, N)
    states[1] = θ = θinit
    times = DataFrame(time=Float64[], early=Bool[])
    tic()
    loglikelihoods[1] = ll = doSMC(a, [θ], targetε, nparticles=nparticles, quiet=true, k=k).llest
    t = toq()
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            tic()
            p = doSMC(a, [θprop], targetε, nparticles=nparticles, llBound=llBound, quiet=true, k=k)
            t = toq()
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
                push!(times, (t, false))
            else
                push!(times, (t, true))
            end
        end
        states[i], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods, times
end

##Generate observations
srand(1);
d = 25;
yobs = randn(d) .* 3.0;
a = RareEventABC.IIDNormalDetails(yobs, d);

rmse(x::Vector, y::Float64) = (x .- y) .^ 2 |> mean |> sqrt;
rmse(x::Vector, y::Float64, w::Vector) = (w .* ((x .- y) .^ 2)) |> sum |> a->a/sum(w) |> sqrt;

##NON-ADAPTIVE RE-ABC
##Tuning
schedule = doSMC(a, [3.0], 5.0, nparticles=300, quiet=false).schedule;
thetas, ll, times = timedPMMH(a, schedule, 200, 300, Uniform(0.0,10.0), 3.0, Normal(0.0,0.5));
sd_hd = sqrt(2.562) * std(thetas)
##Main run
thetas, ll, times = timedPMMH(a, schedule, 2000, 300, Uniform(0.0,10.0), 3.0, Normal(0.0,sd_hd));

##ADAPTIVE RE-ABC
thetas2, ll2, times2 = timedPMMH(a, 5.0, 2000, 300, Uniform(0.0,10.0), 3.0, Normal(0.0,sd_hd));

alltimes = vcat(hcat(times, fill("non-adaptive", size(times)[1])), hcat(times2, fill("adaptive", size(times2[1]))));
rename!(alltimes, :x1, :method);

writetable("pmmh_normal_timed.csv", alltimes) ##To plot (using R)

isearly = alltimes[:early].data;
isadaptive = alltimes[:method].data == "adaptive";

default(xlim=(0., 0.26), ylim=(0, 250), legend=false)
p1 = histogram(alltimes[isearly  & isadaptive, :time], xlim=(0.,0.26), title="Adaptive, stopped early");
p2 = histogram(alltimes[!isearly & isadaptive, :time], xlim=(0.,0.26), title="Adaptive, completed");
p3 = histogram(alltimes[isearly  & !isadaptive, :time], xlim=(0.,0.26), title="Non-adaptive, stopped early");
p4 = histogram(alltimes[!isearly & !isadaptive, :time], xlim=(0.,0.26), title="Non-adaptive, completed");

plot(p1,p2,p3,p4,layout=4)


text:
module MasslessChargedParticleSPARK

    const Δt = 0.1
    const nt = 100000

    const PLOT_DIR = "figures"
    const SYMP_DIR = "symplecticity"

    using GeometricIntegrators
    using GeometricProblems.MasslessChargedParticle

    set_config(:nls_atol, 1E-14)
    set_config(:nls_rtol, 1E-14)
    set_config(:nls_atol_break, 1E6)
    set_config(:nls_rtol_break, 1E6)
    set_config(:nls_stol_break, 1E6)
    set_config(:nls_nmax, 100);

    include("common.jl")
    include("tableau_lists.jl")

    function plot(sol, stages, equ, dir, file, fig_suff)
        if !isdir(dir)
            mkdir(dir)
        end

        plot_massless_charged_particle(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * fig_suff)

        plot_massless_charged_particle_solution(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_solution" * fig_suff)

        plot_massless_charged_particle_traces(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_traces" * fig_suff)

        _plot(sol, stages, equ, dir, file, fig_suff)
    end

end


text:
#WARNING- look, think about market caps and scaling v log
#treat thsi table as a one-off since it has no common elements w/ other tables
function maketable9(cumex::AbstractDataFrame;
    panela::Bool=true,
    panelb::Bool=true,
    panelc::Bool=true)

  scumex::SubDataFrame = view(cumex, cumex.primary, :)

  ###################PANEL A############

  T9_PANEL_A_NAME::String = "t9panela-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"

  T9_PANEL_A_ROWKEYS::Vector{Symbol} = [
    :all,
    :lowyield,
    :midyield,
    :highyield]

  T9_PANEL_A_ROWLABEL::Dict = Dict(
    :all => "All",
    :lowyield=> "\\midrule Low Div Yield \$(\\delta < 0.75\\%)\$",
    :midyield=> "Mid Div Yield \$(0.75\\% < \\delta < 1.5\\%)\$",
    :highyield=> "High Div Yield \$(1.5\\%) < \\delta\$",)

  T9_PANEL_A_COLKEYS::Vector{Symbol} = [
    :yield,
    :riskcum,
    :riskex,
    :riskdiff,
    :rewardcum,
    :rewardex,
    :rewarddiff,
    :N]

  T9_PANEL_A_COLLABEL::Dict = Dict(
    :yield=>"DivYield",
    :riskcum=>"Cum",
    :riskex=>"Ex",
    :riskdiff=>"Diff",
    :rewardcum=>"Cum",
    :rewardex=>"Ex",
    :rewarddiff=>"Diff",
    :N=>"N")

  T9_PANEL_A_COLNAMES::Vector{Vector{String}} =[
    ["", "``Risk'': \$\\sum{|r|}\$", "``Reward'': \$\\sum{r}\$", ""],
    (s::Symbol->T9_PANEL_A_COLLABEL[s]).(T9_PANEL_A_COLKEYS)]

  T9_PANEL_A_WIDTHCOLNAMES::Vector{Vector{Int}} = [[1,3,3,1], ones(Int,length(T9_PANEL_A_COLKEYS))]

  panela && table9panela(scumex,
      tablename = T9_PANEL_A_NAME,
      rowkeys = T9_PANEL_A_ROWKEYS,
      rowlabel = T9_PANEL_A_ROWLABEL,
      colkeys = T9_PANEL_A_COLKEYS,
      colnames = T9_PANEL_A_COLNAMES,
      widthcolnames = T9_PANEL_A_WIDTHCOLNAMES
      )

  ###################PANEL B/C############

  T9_PANEL_BC_ROWKEYS::Vector{Symbol} = [
    :interceptδ,
    :yieldδ,
    :interceptδ1δ,
    :yieldδ1δ,
    :exdayret]

  T9_PANEL_BC_ROWLABEL::Dict = Dict(
    :interceptδ => "Intercept",
    :yieldδ => "Div Yield \$\\delta\$",
    :interceptδ1δ => "\\midrule Intercept",
    :yieldδ1δ => "\$\\delta/(1-\\delta)\$",
    :exdayret => "Ex-Day Return"
    )

  T9_PANEL_BC_COLKEYS::Vector{Symbol} = [
    :coef,
    :se,
    :tstat,
    :stdcoef,
    :mwse,
    :mwtstat]

  T9_PANEL_BC_COLLABEL::Dict = Dict(
    :coef=>"Coef",
    :se=>"StdErr",
    :tstat=>"T-stat",
    :stdcoef=>"StdCoef",
    :mwse=>"StdErr",
    :mwtstat=>"T-stat")

  T9_PANEL_BC_COLNAMES::Vector{Vector{String}} =[
    ["", "Heteroskedastic"],
    (s::Symbol->T9_PANEL_BC_COLLABEL[s]).(T9_PANEL_BC_COLKEYS)]

  T9_PANEL_BC_WIDTHCOLNAMES::Vector{Vector{Int}} = [[4,2], ones(Int,length(T9_PANEL_BC_COLKEYS))]

  T9_PANEL_BC_DECIMALS = 3

  #Now do panel b
  T9_PANEL_B_YSYM = :Dabs
  T9_PANEL_B_YSYM_UNIT = :UDabs
  T9_PANEL_B_NAME::String = "t9panelb-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"

  panelb && table9panelbc(scumex,
    Ysym = T9_PANEL_B_YSYM,
    Ysymunit = T9_PANEL_B_YSYM_UNIT,
    tablename = T9_PANEL_B_NAME,
    rowkeys = T9_PANEL_BC_ROWKEYS,
    rowlabel = T9_PANEL_BC_ROWLABEL,
    colkeys = T9_PANEL_BC_COLKEYS,
    colnames = T9_PANEL_BC_COLNAMES,
    widthcolnames = T9_PANEL_BC_WIDTHCOLNAMES,
    decimals = T9_PANEL_BC_DECIMALS
    )

  #Now do panel c
  T9_PANEL_C_YSYM = :Dret
  T9_PANEL_C_YSYM_UNIT = :UDret
  T9_PANEL_C_NAME::String = "t9panelc-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"
  panelc && table9panelbc(scumex,
    Ysym = T9_PANEL_C_YSYM,
    Ysymunit = T9_PANEL_C_YSYM_UNIT,
    tablename = T9_PANEL_C_NAME,
    rowkeys = T9_PANEL_BC_ROWKEYS,
    rowlabel = T9_PANEL_BC_ROWLABEL,
    colkeys = T9_PANEL_BC_COLKEYS,
    colnames = T9_PANEL_BC_COLNAMES,
    widthcolnames = T9_PANEL_BC_WIDTHCOLNAMES,
    decimals = T9_PANEL_BC_DECIMALS
    )
end

function table9panelbc(cumex::AbstractDataFrame;
  Ysym::Symbol = error("Ysym is required"),
  Ysymunit::Symbol = error("Ysymunit is required"),
  tablename::String = error("tablename is required"),
  rowkeys::Vector{Symbol} = error("rowkeys is required"),
  rowlabel::Dict = error("rowlabel is required"),
  colkeys::Vector{Symbol} = error("colkeys is required"),
  colnames::Vector{Vector{String}} = error("colnames is required"),
  widthcolnames::Vector{Vector{Int}} = error("widthcolnames is required"),
  tablepath::String = TABLE_PATH,
  decimals::Int = DEFAULT_DECIMALS,
  )

  #pre-declare some regresion variables taht will be assigned multiple times
  local Xexpr::FMExpr
  local Xnames::Vector{Symbol}
  local reg::FMLM
  local se::Matrix{Float64}
  local semw::Matrix{Float64}
  local scumex::SubDataFrame

  #housekeeping and pre-allocation
  colidx::Dict = Dict(colkey=>i for (i,colkey) ∈ enumerate(colkeys))
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, sum(widthcolnames[1]))).(1:length(rowkeys))
  rowidx::Dict = Dict(rowkey=>desccontent[i] for (i,rowkey) ∈ enumerate(rowkeys))
  descrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys)

  scumex = view(cumex, :, [
    :yield, :yieldd1d, :exdayret, :Uyield, :Uyieldd1d, :Uexdayret, Ysym, Ysymunit])
  scumex = view(scumex, completecases(scumex), :)

  #constructs a row for the table from the regression object and the name of the coefficient
  function regrow!(reg::FMLM,  Fcoef::Symbol, descrow::Vector{String},
    Σ::Matrix{Float64}, Σmw::Matrix{Float64})

    ind::Int = findfirst(isequal(Fcoef), reg.Xnames)
    @assert Fcoef == reg.Xnames[ind]

    #display(Σ)

    coef::Float64 = reg.β[ind]
    sigma::Float64 = Σ[ind,ind] ^0.5
    tstat::Float64 = coef/sigma
    mwsigma::Float64 = Σmw[ind,ind]^0.5
    mwtstat::Float64 = coef/mwsigma

    descrow[colidx[:coef]] = brnum(coef, decimals, scalefactor=100.)
    descrow[colidx[:se]] = num2str(sigma, decimals, scalefactor=100.)
    descrow[colidx[:tstat]] = brnum(tstat, decimals, scalefactor=1.)
    descrow[colidx[:mwse]] = num2str(mwsigma, decimals, scalefactor=100.)
    descrow[colidx[:mwtstat]] = brnum(mwtstat, decimals, scalefactor=1.)
  end

  #do the bivariate regressions
  Xexpr = :yield
  Xnames= [:intercept, :yield]
  reg = FMLM(scumex,  Xexpr, Ysym, Xnames=Xnames, Yname=Ysym, containsmissings=false)
  se = homoskedasticΣ!(reg)
  semw = modifiedwhiteΣ!(reg)
  regrow!(reg, :intercept, rowidx[:interceptδ], se, semw)
  regrow!(reg, :yield, rowidx[:yieldδ], se, semw)

  #unit-standardized bivariate regression
  Xexpr = :Uyield
  reg = FMLM(scumex,  Xexpr, Ysymunit, Xnames=Xnames, Yname=Ysym, containsmissings=false)
  @assert (reg.Xnames[2] == :yield) #name check
  rowidx[:yieldδ][colidx[:stdcoef]] = brnum(reg.β[2], decimals, scalefactor=100.)
  rowidx[:interceptδ][colidx[:stdcoef]] = ""

  #multi-variate regression
  Xexpr = Meta.parse("yieldd1d + exdayret")
  Xnames= [:intercept, :yieldd1d, :exdayret]
  reg = FMLM(scumex,  Xexpr, Ysym, Xnames=Xnames, Yname=Ysym, containsmissings=false)
  se = homoskedasticΣ!(reg)
  semw = modifiedwhiteΣ!(reg)
  regrow!(reg, :intercept, rowidx[:interceptδ1δ], se, semw)
  regrow!(reg, :yieldd1d, rowidx[:yieldδ1δ], se, semw)
  regrow!(reg, :exdayret, rowidx[:exdayret], se, semw)

  #unit-standardized multivariate regression
  Xexpr = Meta.parse("Uyieldd1d + Uexdayret")
  reg = FMLM(scumex,  Xexpr, Ysymunit, Xnames=Xnames, Yname=Ysym, containsmissings=false)
  @assert (reg.Xnames[2] == :yieldd1d) && (reg.Xnames[3] == :exdayret)
  rowidx[:yieldδ1δ][colidx[:stdcoef]] = brnum(reg.β[2], decimals, scalefactor=100.)
  rowidx[:exdayret][colidx[:stdcoef]] = brnum(reg.β[3], decimals, scalefactor=100.)
  rowidx[:interceptδ1δ][colidx[:stdcoef]] = ""

  tt::String = textable(colnames=colnames,
    descrownames=descrownames,
    desccontent=desccontent,
    #alignmentstring=alignmentstring,
    widthcolnames=widthcolnames)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

function table9panela(cumex::AbstractDataFrame;
  tablename::String = error("tablename is required"),
  rowkeys::Vector{Symbol} = error("rowkeys is required"),
  rowlabel::Dict = error("rowlabel is required"),
  colkeys::Vector{Symbol} = error("colkeys is required"),
  colnames::Vector{Vector{String}} = error("colnames is required"),
  widthcolnames::Vector{Vector{Int}} = error("widthcolnames is required"),
  tablepath::String = TABLE_PATH,
  decimals::Int = DEFAULT_DECIMALS,
  )

  colidx::Dict = Dict(colkey=>i for (i,colkey) ∈ enumerate(colkeys))
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, sum(widthcolnames[1]))).(1:length(rowkeys))
  rowidx::Dict = Dict(rowkey=>desccontent[i] for (i,rowkey) ∈ enumerate(rowkeys))
  descrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys)

  @inline function tablerow!(df::AbstractDataFrame, descrow::Vector{String})
    #first get our datasets

    #get the values
    yield::Float64 = mean(df.yield)
    cumrisk::Float64 = mean(df.cumabs)
    exrisk::Float64 = mean(df.exabs)
    Δrisk::Float64 = cumrisk - exrisk

    cumret::Float64 = mean(df.cumret)
    exret::Float64 = mean(df.exret)
    Δret::Float64 = cumret - exret
    N::Int = size(df,1)

    #write the row
    descrow[colidx[:yield]] = num2str(yield, decimals, scalefactor=100.0)
    descrow[colidx[:riskcum]] = brnum(cumrisk, decimals, scalefactor=100.0)
    descrow[colidx[:riskex]] = brnum(exrisk, decimals, scalefactor=100.0)
    descrow[colidx[:riskdiff]] = brnum(Δrisk, decimals, scalefactor=100.0)
    descrow[colidx[:rewardcum]] = brnum(cumret, decimals, scalefactor=100.0)
    descrow[colidx[:rewardex]] = brnum(exret, decimals, scalefactor=100.0)
    descrow[colidx[:rewarddiff]] = brnum(Δret, decimals, scalefactor=100.0)
    descrow[colidx[:N]] = num2str(N, decimals, Ints=true)
  end

  #write the row using the appropriate data
  tablerow!(cumex, rowidx[:all])

  scumex = view(cumex,cumex.Gyield .== :lowyield, :)
  tablerow!(scumex, rowidx[:lowyield])

  scumex = view(cumex,cumex.Gyield .== :medyield, :)
  tablerow!(scumex, rowidx[:midyield])

  scumex = view(cumex,cumex.Gyield .== :highyield, :)
  tablerow!(scumex, rowidx[:highyield])

  tt::String = textable(colnames=colnames,
    descrownames=descrownames,
    desccontent=desccontent,
    #alignmentstring=alignmentstring,
    widthcolnames=widthcolnames)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

function aggregateoverdays(F::Function, dist::AbstractDataFrame, dayrange::AbstractRange;
  Fgroup::Symbol = :eid, selectcols::Vector{Symbol} = names(dist), suffix::String = "",
  colnames::Vector{Symbol} = names(dist), Fcriteria::Symbol = :day)

  sdist::SubDataFrame = view(dist, (g::Union{MInt, Symbol}->
    (!ismissing(g)) && g ∈ dayrange).(dist[!,Fcriteria]), selectcols)

  agg::DataFrame = aggregate(groupby(sdist, Fgroup), F)
  names!(agg, colnames)

  return agg
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Oceananigans
using OceanTurbulenceParameterEstimation
using Distributions
using DataDeps
using GLMakie

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

θ_variable_Pr = (
                 CᵂwΔ  = 8.20,
                 Cᵂu★  = 6.77,
                 Cᴰ    = 1.69,
                 Cᴸᵇ   = 1.15,
                 Cᴷu⁻  = 0.0924,
                 Cᴷc⁻  = 0.288,
                 Cᴷe⁻  = 2.46,
                 Cᴷuʳ  = 1.32,
                 Cᴷcʳ  = 0.0404,
                 Cᴷeʳ  = 4.27,
                 CᴷRiʷ = 3.15,
                 CᴷRiᶜ = 2.08,
                 Cᴬu   = 0.0,
                 Cᴬc   = 0.0,
                 Cᴬe   = 0.0,
                 Cᵟu   = 0.5,
                 Cᵟc   = 0.5,
                 Cᵟe   = 0.5,
                )

function perturbation_prior(θ★, ϵ=0.5)
    L = (1 - ϵ) * θ★
    U = (1 + ϵ) * θ★
    return ScaledLogitNormal(bounds=(L, U))
end

θ★ = θ_constant_Ri
priors = NamedTuple(name => perturbation_prior(θ★[name]) for name in keys(θ★))
free_parameters = FreeParameters(priors)

case_path(case) = @datadep_str("six_day_suite_1m/$(case)_instantaneous_statistics.jld2")

cases = ["free_convection",
         "weak_wind_strong_cooling",
         "strong_wind_weak_cooling",
         "strong_wind",
         "strong_wind_no_rotation"]

field_names = (:u, :v, :T, :e)
regrid_size = nothing
observations = [SyntheticObservations(case_path(case); field_names, regrid_size)
                for case in cases]

observation = first(observations)
α = observation.metadata.parameters.thermal_expansion_coefficient
equation_of_state = LinearEquationOfState(; α)

mixing_length = MixingLength(Cᴬu   = 0.0,
                             Cᴬc   = 0.0,
                             Cᴬe   = 0.0,
                             CᴷRiᶜ = 2.0,
                             Cᴷuʳ  = 0.0,
                             Cᴷcʳ  = 0.0,
                             Cᴷeʳ  = 0.0)

catke = CATKEVerticalDiffusivity(; mixing_length)

Nensemble = 20

simulation = ensemble_column_model_simulation(observations;
                                              Nensemble,
                                              architecture = CPU(),
                                              buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=true),
                                              tracers = (:T, :e),
                                              closure = catke)

simulation.Δt = 30.0

progress(sim) = @info "Iter: $(iteration(sim)), time: $(prettytime(sim))"
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
Qᵀ = simulation.model.tracers.T.boundary_conditions.top.condition
dTdz = simulation.model.tracers.T.boundary_conditions.bottom.condition

for (case, obs) in enumerate(observations)
    @show case cases[case]
    @show obs.metadata.parameters.momentum_flux
    @show obs.metadata.parameters.temperature_flux
    @show obs.metadata.parameters.buoyancy_flux
    @show f = obs.metadata.parameters.coriolis_parameter

    view(Qᵘ, :, case) .= obs.metadata.parameters.momentum_flux
    view(Qᵀ, :, case) .= obs.metadata.parameters.temperature_flux
    view(dTdz, :, case) .= obs.metadata.parameters.dθdz_deep
    view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
end

calibration = InverseProblem(observations, simulation, free_parameters)

model = simulation.model

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, model.diffusivity_fields),
                     prefix = "catke_simulation",
                     schedule = SpecifiedTimes(round.(observation_times(obs))...),
                     force = true)

θ_perturbed = [NamedTuple(name => rand(priors[name]) for name in keys(priors)) for k=1:Nensemble-1]
θ = [θ★, θ_perturbed...]

forward_run!(calibration, θ; suppress=false)

output_path = simulation.output_writers[:fields].filepath
u = FieldTimeSeries(output_path, "u")
v = FieldTimeSeries(output_path, "v")
T = FieldTimeSeries(output_path, "T")
e = FieldTimeSeries(output_path, "e")

z = znodes(u)

fig = Figure(resolution=(1200, 1200))

ax_T = [Axis(fig[1, i], xlabel = "Temperature (ᵒC)", ylabel = "z (m)") for i = 1:2]
ax_u = [Axis(fig[2, i], xlabel = "Velocity (m s⁻¹)", ylabel = "z (m)") for i = 1:2]
ax_e = [Axis(fig[3, i], xlabel = "Turbulent Kinetic Energy (m² s⁻²)", ylabel = "z (m)") for i = 1:2]

Nt = minimum(length(obs.times) for obs in observations) - 5
slider = Slider(fig[0, :], range=1:Nt, horizontal=true, startvalue=1)
n = slider.value

profile_max(u, i, n) = Point2f.(maximum(interior(u[n])[:, i, :], dims=1)[:], z)
profile_min(u, i, n) = Point2f.(minimum(interior(u[n])[:, i, :], dims=1)[:], z)

for i = 1:length(cases)
    u★ = @lift interior(u[$n])[1, i, :]
    v★ = @lift interior(v[$n])[1, i, :]
    T★ = @lift interior(T[$n])[1, i, :]
    e★ = @lift interior(e[$n])[1, i, :]

    u_max = @lift profile_max(u, i, $n)
    v_max = @lift profile_max(v, i, $n)
    T_max = @lift profile_max(T, i, $n)
    e_max = @lift profile_max(e, i, $n)

    u_min = @lift profile_min(u, i, $n)
    v_min = @lift profile_min(v, i, $n)
    T_min = @lift profile_min(T, i, $n)
    e_min = @lift profile_min(e, i, $n)

    obs = observations[i]
    u_truth = @lift interior(obs.field_time_serieses.u[$n])[1, 1, :]
    v_truth = @lift interior(obs.field_time_serieses.v[$n])[1, 1, :]
    T_truth = @lift interior(obs.field_time_serieses.T[$n])[1, 1, :]
    e_truth = @lift interior(obs.field_time_serieses.e[$n])[1, 1, :]

              lines!(ax_T[i], T_truth, z, color=(:gray23, 0.6), linewidth=5,   label="LES")
              lines!(ax_u[i], u_truth, z, color=(:gray23, 0.6), linewidth=5,   label="u, LES")
    i == 1 && lines!(ax_u[i], v_truth, z, color=(:gray23, 0.6), linewidth=1.5, label="v, LES")
              lines!(ax_e[i], e_truth, z, color=(:gray23, 0.6), linewidth=5,   label="LES")

    color = (:pink, 0.4)
    band!(ax_T[i], T_max, T_min; color)
    band!(ax_e[i], e_max, e_min; color)
    band!(ax_u[i], u_max, u_min; color)
    i == 1 && band!(ax_u[i], v_max, v_min; color)

              lines!(ax_T[i], T★, z, color=:purple1, linewidth=3,   label="CATKE")
              lines!(ax_u[i], u★, z, color=:purple1, linewidth=3,   label="u, CATKE")
    i == 1 && lines!(ax_u[i], v★, z, color=:purple1, linewidth=1.5, label="v, CATKE")
              lines!(ax_e[i], e★, z, color=:purple1, linewidth=3,   label="CATKE")
end

for i = 1:length(cases)
    xlims!(ax_T[i], 19.6, 20.05)
    xlims!(ax_u[i], -0.6, 0.6)
    xlims!(ax_e[i], -0.0005, 0.005)
    axislegend(ax_T[i], position=:rb)
    axislegend(ax_u[i], position=:rb)
    axislegend(ax_e[i], position=:rb)
end

display(fig)


text:
# schroedinger.jl
##########################
using FFTW, DSP, LinearAlgebra, Plots

# initial data
N=256;
n=128;
M=1000;
u=zeros(ComplexF64, N,M);
for p=65:192
    u[p,1]=exp(2*pi*im*p/128)+exp(2*pi*im*p^2/12800);
end

# setting
E=zeros(ComplexF64, N, N);
for k=1:N
    E[k,k]=1
end

t=1/100;

# sqrt(M) times Fourier matrix
omega=exp(im*2*pi/N);
F=zeros(ComplexF64, N, N);
for k=1:N
    for l=1:N
    F[k,l]=omega^(-(k-1)*(l-1));
    end
end

D=zeros(ComplexF64, N, N)
for k=1:n
    D[k,k]=(k-1)^2;
    D[k+n,k+n]=(n-k)^2;
end

G=-im*(4*pi^2*t*N^(-3))*F'*D*F
A=E+G+G^2/2+G^3/6;

for l=2:M
    u[:,l]=A*u[:,l-1];
end

v=zeros(ComplexF64, N,M)
for l=1:M
    v[:,l]=u[:,M-l+1];
end

p1=plot(real.(u'),
   title="The real part of Schroedinger flow",
    st=:heatmap,
    color=:jet,
    yaxis="time t",
    xaxis="position x")
plot!(yticks = ([0 1000;], [0 10]))
plot!(xticks = ([0 256;], [0 256]))

p2=plot(imag.(-u'),
   title="The imaginary part of Schredinger flow",
    st=:heatmap,
    color=:jet,
    yaxis="time t",
    xaxis="position x")
plot!(yticks = ([0 1000;], [0 10]))
plot!(xticks = ([0 256;], [0 256]))

q1=surface(real.(v),
        title="The real part of Schroedinger flow",
        camera=(80,60),
        color=:jet,
        xaxis="time t",
        yaxis="position x")
plot!(xticks = ([0 1000;], [10 0]))
plot!(yticks = ([0 256;], [0 256]))

q2=surface(imag.(v),
        title="The imaginary part of Schroedinger flow",
        camera=(80,60),
        color=:jet,
        xaxis="time t",
        yaxis="position x")
plot!(xticks = ([0 1000;], [10 0]))
plot!(yticks = ([0 256;], [0 256]))

plot(p1, q1, p2, q2,
     layout=(2,2),
     size=(1600,1300),
     margin=Plots.Measures.Length(:mm, 10.0))
savefig("schroedinger.png")


text:
using DelimitedFiles

function output_ttvfast_times()
    data = readdlm("ttvfast/Times", ' ')
    data = sortslices(data, by=x->x[1], dims=1)
    masks = [data[:,1] .== x for x in collect(0.0:1.0:6.0)]
    times = [data[mask,3] for mask in masks]
    tts = zeros(7, length(times[1]))
    for i in 1:7
        tts[i,1:length(times[i])] = times[i][1:length(times[i])]
    end
    open("ttvfast_tts.txt", "w") do io
        writedlm(io, eachrow(tts), ',')
    end
end

max_forks_repo_name: langfzac/nbg-papers


text:
function _dft_mtime_old_state(e::CustomEvent, file::String)
    get!(e, file, 0.0)
end

function _dft_mtime_new_state(::CustomEvent, file::String)
    mtime(file)
end

function _dft_mtime_update!(e::CustomEvent, file::String)
    new_mtime = new_state(e, file)
    setindex!(e, new_mtime, file)
    return e
end

function _dft_mtime_trigger(::CustomEvent, old_mtime, new_mtime)
    return old_mtime < new_mtime
end

function FileMTimeEvent()
    return CustomEvent{String, Float64}(;
        old_state = _dft_mtime_old_state,
        new_state = _dft_mtime_new_state,
        update! = _dft_mtime_update!,
        trigger = _dft_mtime_trigger,
    )
end


text:
using Compat, Compat.Test, PyCall

py"""
def mklist(*args, **kwargs):
    l = list(args)
    l.extend(kwargs.items())
    return l
"""
@testset "pycall!" begin
    pymklist = py"mklist"
    ret = PyNULL()

    function pycall_checks(res, pyf, RetType, args...; kwargs...)
        pycall_res = pycall(pyf, RetType, args...; kwargs...)
        res_converted = pycall!(res, pyf, RetType, args...; kwargs...)
        @test res_converted == pycall_res
        @test convert(RetType, res) == res_converted
        RetType != PyAny && @test res_converted isa RetType
    end

    @testset "basics" begin
        args = ("a", 2, 4.5)
        for RetType in (PyObject, PyAny, Tuple)
            pycall_checks(ret, pymklist, RetType, args...)
            GC.gc()
        end
    end

    @testset "kwargs" begin
        args = ("a", 2, 4.5)
        for RetType in (PyObject, PyAny, Tuple)
            pycall_checks(ret, pymklist, RetType, args...; b=19610, c="5")
            GC.gc()
        end
    end

end


text:
export
    LoopyBP,
    loopy_BP

"""
    LoopyBP

An instant algorithm that runs loopy belief propagation.

# Arguments
- default_range_size: The size to use as default when calling `support` on a node.
- epsilon: The allowable difference between beliefs on successive iterations
for termination.
- maxiterations: The maximum number of iterations to run. `infer` will terminate if
this number of iterations is reached, even if it has not converged.
"""
struct LoopyBP <: BP
    default_range_size::Int
    epsilon::Float64
    maxiterations::Int
    LoopyBP(drs = 10, epsilon = 0.0001, maxiterations = 10) = new(drs, epsilon, maxiterations)
end

function loopy_BP(runtime::Runtime; default_range_size = 10, epsilon = 0.0001, maxiterations = 10)
    network = get_network(runtime)
    for node in topsort(get_initial_graph(network))
        remove_messages!(runtime, node, :pi_message)
        remove_messages!(runtime, node, :lambda_message)
    end
    run_bp(LoopyBP(default_range_size, epsilon, maxiterations), runtime)
end

function run_bp(algorithm::LoopyBP, runtime)
    network = get_network(runtime)
    ranges = Dict()
    nodes = get_nodes(network)

    for node in nodes
        inst = current_instance(runtime, node)
        rng = get_range(runtime, inst)
        ranges[node.name] = rng
    end

    for node in nodes
        for par in get_parents(network, node)
            if par isa Variable
                set_message!(runtime, par, node, :pi_message,
                    ones(length(ranges[par.name])))
            end
        end
    end

    newbeliefs = initial_pass(runtime, network, ranges)
    conv = false
    iteration = 0

    while !conv && iteration < algorithm.maxiterations
        oldbeliefs = copy(newbeliefs)
        backward_pass(runtime, network, ranges)
        forward_pass(runtime, network, ranges, newbeliefs)
        conv = converged_loopy(get_variables(network), ranges, newbeliefs, oldbeliefs, algorithm.epsilon)
        iteration += 1
    end
end

function initialize(runtime, network, ranges)
end

function initial_pass(runtime, network, ranges)
    newbeliefs = Dict{Variable, Dist}()
    variables = [v for v in topsort(get_initial_graph(network)) if v isa Variable]
    for var in variables
        inst = current_instance(runtime, var)
        sf = get_sfunc(inst)
        pars = get_parents(network, var)
        incoming_pis :: Vector{Dist} =
            collect_messages(runtime, pars, var, :pi_message)
        range = ranges[var.name]
        parranges = [ranges[p.name] for p in pars]
        pi = compute_pi(sf, range, tuple(parranges...), tuple(incoming_pis...))
        # on the first pass, we interpret evidence as lambda message
        # coming from elsewhere, so include it in the pi
        if has_evidence(runtime, inst)
            evidence = get_evidence(runtime, inst)
            pi = Cat(range, [get_score(evidence, x) for x in range])
        end
        if has_intervention(runtime, inst)
            intervention = get_intervention(runtime, inst)
            pi = intervention
        end
        set_value!(runtime, inst, :pi, pi)
        newbeliefs[var] = pi
        for ch in get_children(network, var)
            set_message!(runtime, var, ch, :pi_message, pi)
        end
    end
    return newbeliefs
end

function backward_pass(runtime, network, ranges)
    variables = [v for v in topsort(get_initial_graph(network)) if v isa Variable]
    for var in reverse(variables)
        _backstep(runtime, var, ranges)
    end
end

function forward_pass(runtime, network, ranges, newbeliefs)
    variables = [v for v in topsort(get_initial_graph(network)) if v isa Variable]
    for var in variables
        _forwardstep(runtime, var, ranges, false)
        inst = current_instance(runtime, var)
        newbeliefs[var] = get_belief(runtime, inst)
    end
end

function converged_loopy(variables, ranges, new_beliefs, old_beliefs, epsilon::Float64)
    total_diff = 0.0
    total_len = 0
    for var in variables
        range = ranges[var.name]
        diffs = [abs(cpdf(new_beliefs[var], (), x) - cpdf(old_beliefs[var], (), x)) for x in range]
        total_diff += sum(diffs)
        total_len += length(range)
    end
    return total_diff / total_len < epsilon
end


text:
function Base.show(io::IO, K::KummerExt)
  if isdefined(K.AutG, :snf)
    print(io, "KummerExt with structure $(K.AutG.snf)")
  else
    print(io, "KummerExt with structure $([K.AutG.rels[i, i] for i=1:ngens(K.AutG)])")
  end
end

@doc Markdown.doc"""
    kummer_extension(n::Int, gens::Vector{FacElem{nf_elem, AnticNumberField}}) -> KummerExt

Creates the Kummer extension of exponent $n$ generated by the elements in `gens`.
"""
function kummer_extension(n::Int, gen::Vector{FacElem{nf_elem, AnticNumberField}})
  K = KummerExt()
  k = base_ring(gen[1])
  zeta, o = torsion_units_gen_order(k)
  @assert o % n == 0
  K.zeta = zeta^div(o, n)
  K.n = n
  K.gen = gen
  K.AutG = GrpAbFinGen(fmpz[n for i=gen])
  K.frob_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
  return K
end

function kummer_extension(exps::Array{Int, 1}, gens::Vector{FacElem{nf_elem, AnticNumberField}})
  K = KummerExt()
  k = base_ring(gens[1])
  zeta, o = torsion_units_gen_order(k)
  n = lcm(exps)
  @assert o % n == 0

  K.zeta = zeta^div(o, n)
  K.n = n
  K.gen = gens
  K.AutG = abelian_group(exps)
  K.frob_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
  return K
end

function kummer_extension(n::Int, gen::Array{nf_elem, 1})
  g = FacElem{nf_elem, AnticNumberField}[FacElem(x) for x in gen]
  return kummer_extension(n, g)
end

###############################################################################
#
#  Base Field
#
###############################################################################

function base_field(K::KummerExt)
  return base_ring(K.gen[1])
end

###############################################################################
#
#  Exponent of a Kummer extension
#
###############################################################################

function exponent(K::KummerExt)
  return Int(exponent(K.AutG))
end

###############################################################################
#
#  Degree
#
###############################################################################

function degree(K::KummerExt)
  return Int(order(K.AutG))
end

###############################################################################
#
#  IsCyclic
#
###############################################################################

function iscyclic(K::KummerExt)
  return isone(length(K.gen)) || iscyclic(K.AutG)
end

###############################################################################
#
#  From Kummer Extension to Number Field
#
###############################################################################

function number_field(K::KummerExt)
  k = base_field(K)
  kt = PolynomialRing(k, "t", cached = false)[1]
  pols = Array{elem_type(kt), 1}(undef, length(K.gen))
  for i = 1:length(pols)
    p = Vector{nf_elem}(undef, Int(order(K.AutG[i]))+1)
    p[1] = -evaluate(K.gen[i])
    for i = 2:Int(order(K.AutG[i]))
      p[i] = zero(k)
    end
    p[end] = one(k)
    pols[i] = kt(p)
  end
  return number_field(pols, check = false, cached = false)
end

###############################################################################
#
#  Computation of Frobenius automorphisms
#
###############################################################################

function assure_gens_mod_nth_powers(K::KummerExt)
  if isdefined(K, :gen_mod_nth_power)
    return nothing
  end
  gens = Vector{FacElem{nf_elem, AnticNumberField}}(undef, length(K.gen))
  for i = 1:length(gens)
    gens[i] = RelSaturate._mod_exponents(K.gen[i], K.n)
  end
  K.gen_mod_nth_power = gens
  return nothing
end

@doc Markdown.doc"""
    canonical_frobenius(p::NfOrdIdl, K::KummerExt) -> GrpAbFinGenElem

Computes the element of the automorphism group of $K$ corresponding to the
Frobenius automorphism induced by the prime ideal $p$ of the base field of $K$.
It fails if the prime is an index divisor or if $p$ divides the given generators
of $K$
"""
function canonical_frobenius(p::NfOrdIdl, K::KummerExt)
  @assert norm(p, copy = false) % K.n == 1
  if haskey(K.frob_cache, p)
    return K.frob_cache[p]
  end
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end
  if !fits(Int, minimum(p, copy = false))
    return canonical_frobenius_fmpz(p, K)
  end
  assure_gens_mod_nth_powers(K)
  if degree(p) != 1
    F, mF = ResidueFieldSmall(Zk, p)
    mF1 = NfToFqNmodMor_easy(mF, number_field(Zk))
    aut = _compute_frob(K, mF1, p)
  else
    F2, mF2 = ResidueFieldSmallDegree1(Zk, p)
    mF3 = NfToGFMor_easy(mF2, number_field(Zk))
    aut = _compute_frob(K, mF3, p)
  end
  z = K.AutG(aut)
  K.frob_cache[p] = z
  return z
end

function _compute_frob(K, mF, p)
  z_p = image(mF, K.zeta)^(K.n-1)

  # K = k(sqrt[n_i](gen[i]) for i=1:length(gen)), an automorphism will be
  # K[i] -> zeta^divexact(n, n_i) * ? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  aut = Array{fmpz, 1}(undef, length(K.gen))
  for j = 1:length(K.gen)
    ord_genj = Int(order(K.AutG[j]))
    ex = div(norm(p, copy = false)-1, ord_genj)
    mu = image(mF, K.gen_mod_nth_power[j])^ex
    i = 0
    z_pj = z_p^divexact(K.n, ord_genj)
    while !isone(mu)
      i += 1
      @assert i <= K.n
      mu = mul!(mu, mu, z_pj)
    end
    aut[j] = fmpz(i)
  end
  return aut
end

function canonical_frobenius_fmpz(p::NfOrdIdl, K::KummerExt)
  @assert norm(p, copy = false) % K.n == 1
  if haskey(K.frob_cache, p)
    return K.frob_cache[p]
  end
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end

  F, mF = ResidueField(Zk, p)
  #_mF = extend_easy(mF, number_field(Zk))
  mF1 = NfToFqMor_easy(mF, number_field(Zk))
  z_p = image(mF1, K.zeta)^(K.n-1)

  # K = k(sqrt[n_i](gen[i]) for i=1:length(gen)), an automorphism will be
  # K[i] -> zeta^divexact(n, n_i) * ? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  aut = Array{fmpz, 1}(undef, length(K.gen))
  for j = 1:length(K.gen)
    ord_genj = Int(order(K.AutG[j]))
    ex = div(norm(p, copy = false)-1, ord_genj)
    mu = image(mF1, K.gen[j], K.n)^ex  # can throw bad prime!
    i = 0
    z_pj = z_p^divexact(K.n, ord_genj)
    while !isone(mu)
      i += 1
      @assert i <= K.n
      mul!(mu, mu, z_pj)
    end
    aut[j] = fmpz(i)
  end
  z = K.AutG(aut)
  K.frob_cache[p] = z
  return z
end

#In this function, we are computing the image of $sqrt[n](g) under the Frobenius automorphism of p
function canonical_frobenius(p::NfOrdIdl, K::KummerExt, g::FacElem{nf_elem})
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    throw(BadPrime(p))
  end

  if !fits(Int, minimum(p, copy = false))
    error("Oops")
  end

  @assert norm(p, copy = false) % K.n == 1
  ex = div(norm(p, copy = false)-1, K.n)

  #K = sqrt[n](gen), an automorphism will be
  # K[i] -> zeta^? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  if degree(p) != 1
    F, mF = ResidueFieldSmall(Zk, p)
    mF1 = extend_easy(mF, nf(Zk))
    z_p = inv(mF1(K.zeta))
    mu = image(mF1, g, K.n)^ex  # can throw bad prime!
    i = 0
    while true
      if isone(mu)
        break
      end
      i += 1
      @assert i <= K.n
      mu = mul!(mu, mu, z_p)
    end
    return i
  else
    F2, mF2 = ResidueFieldSmallDegree1(Zk, p)
    mF3 = extend_easy(mF2, nf(Zk))
    z_p1 = inv(mF3(K.zeta))
    mu1 = image(mF3, g, K.n)^ex  # can throw bad prime!
    i = 0
    while true
      if isone(mu1)
        break
      end
      i += 1
      @assert i <= K.n
      mu1 = mul!(mu1, mu1, z_p1)
    end
    return i
  end
end

################################################################################
#
#  Frobenius for cft
#
################################################################################

# In this context, we are computing the Frobenius for conjugate prime ideals
# We save the projection of the factor base, we can reuse them
#Computes a set of prime ideals of the base field of K such that the corresponding Frobenius
#automorphisms generate the automorphism group
function find_gens(K::KummerExt, S::PrimesSet, cp::fmpz=fmpz(1))
  if isdefined(K, :frob_gens)
    return K.frob_gens[1], K.frob_gens[2]
  end
  k = base_field(K)
  ZK = maximal_order(k)
  R = K.AutG
  sR = Vector{GrpAbFinGenElem}(undef, length(K.gen))
  lp = Vector{NfOrdIdl}(undef, length(K.gen))

  indZK = index(ZK)
  q, mq = quo(R, GrpAbFinGenElem[], false)
  s, ms = snf(q)
  ind = 1
  threshold = max(div(degree(k), 5), 5)

  for p in S
    if cp % p == 0 || indZK % p == 0
      continue
    end
    @vprint :ClassField 2 "Computing Frobenius over $p\n"
    lP = prime_decomposition(ZK, p)
    LP = NfOrdIdl[P for (P, e) in lP if degree(P) < threshold]
    if isempty(LP)
      continue
    end
    #Compute the projections of the gens as gfp_poly.
    #I can use these projections for all the prime ideals, saving some time.
    f = R[1]
    D = Vector{Vector{gfp_poly}}(undef, length(K.gen))
    for i = 1:length(D)
      D[i] = Vector{gfp_poly}(undef, length(K.gen[i].fac))
    end

    first = false
    for P in LP
      try
        f = _canonical_frobenius_with_cache(P, K, first, D)
        @hassert :ClassField 1 f == canonical_frobenius(P, K)
        first = true
      catch e
        if !isa(e, BadPrime)
          rethrow(e)
        end
        continue
      end
      if iszero(mq(f))
        continue
      end
      #At least one of the coefficient of the element
      #must be invertible in the snf form.
      el = ms\f
      to_be = false
      for w = 1:ngens(s)
        if gcd(s.snf[w], el.coeff[w]) == 1
          to_be = true
          break
        end
      end
      if !to_be
        continue
      end
      sR[ind] = f
      lp[ind] = P
      ind += 1
      q, mq = quo(R, sR[1:ind-1], false)
      s, ms = snf(q)
    end
    if order(s) == 1
      break
    end
    @vprint :ClassField 3 "Index: $(exponent(s))^($(valuation(order(s), exponent(s))))\n"
  end
  K.frob_gens = (lp, sR)
  return lp, sR
end

function _canonical_frobenius_with_cache(p::NfOrdIdl, K::KummerExt, cached::Bool, D::Vector{Vector{gfp_poly}})
  @assert norm(p, copy = false) % K.n == 1
  if haskey(K.frob_cache, p)
    return K.frob_cache[p]
  end
  Zk = order(p)

  assure_gens_mod_nth_powers(K)

  if degree(p) != 1
    F, mF = ResidueFieldSmall(Zk, p)
    mF1 = NfToFqNmodMor_easy(mF, number_field(Zk))
    aut = _compute_frob(K, mF1, p, cached, D)
  else
    F2, mF2 = ResidueFieldSmallDegree1(Zk, p)
    mF3 = NfToGFMor_easy(mF2, number_field(Zk))
    aut = _compute_frob(K, mF3, p, cached, D)
  end
  z = K.AutG(aut)
  K.frob_cache[p] = z
  return z
end

function _compute_frob(K, mF, p, cached, D)
  z_p = image(mF, K.zeta)^(K.n-1)

  # K = k(sqrt[n_i](gen[i]) for i=1:length(gen)), an automorphism will be
  # K[i] -> zeta^divexact(n, n_i) * ? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p
  aut = Array{fmpz, 1}(undef, length(K.gen))
  for j = 1:length(K.gen)
    ord_genj = Int(order(K.AutG[j]))
    ex = div(norm(p, copy = false)-1, ord_genj)
    mu = image(mF, K.gen_mod_nth_power[j], D[j], cached, K.n)^ex  # can throw bad prime!
    i = 0
    z_pj = z_p^divexact(K.n, ord_genj)
    while !isone(mu)
      i += 1
      @assert i <= K.n
      mu = mul!(mu, mu, z_pj)
    end
    aut[j] = fmpz(i)
  end
  return aut
end

################################################################################
#
#  IsSubfield
#
################################################################################

@doc Markdown.doc"""
    issubfield(K::KummerExt, L::KummerExt) -> Bool, Vector{Tuple{nf_elem, Vector{Int}}}

Given two kummer extensions of a base field $k$, returns true and the data
to define an injection from $K$ to $L$ if $K$ is a subfield of $L$. Otherwise
the function returns false and some meaningless data.
"""
function issubfield(K::KummerExt, L::KummerExt)
  @assert base_field(K) == base_field(L)
  @assert divisible(exponent(L), exponent(K))
  #First, find prime number that might be ramified.
  norms = Vector{fmpz}(undef, length(K.gen)+length(L.gen)+1)
  for i = 1:length(K.gen)
    norms[i] = numerator(norm(K.gen[i]))
  end
  for i = 1:length(L.gen)
    norms[i+length(K.gen)] = numerator(norm(L.gen[i]))
  end
  norms[end] = fmpz(exponent(L))
  norms = coprime_base(norms)
  coprime_to = lcm(norms)
  res = Vector{Tuple{FacElem{nf_elem, AnticNumberField}, Vector{Int}}}(undef, length(K.gen))
  lP = find_gens(L, Vector{FacElem{nf_elem, AnticNumberField}}[K.gen], coprime_to)
  for i = 1:length(K.gen)
    fl, coord, rt = _find_embedding(L, K.gen[i], Int(order(K.AutG[i])), lP)
    if !fl
      return fl, res
    end
    res[i] = (rt, Int[Int(coord[j]) for j = 1:length(L.gen)])
  end
  return true, res
end

################################################################################
#
#  Kummer failure
#
################################################################################

@doc Markdown.doc"""
    kummer_failure(x::nf_elem, M::Int, N::Int) -> Int

Computes the quotient of $N$ and $[K(\zeta_M, \sqrt[N](x))\colon K(\zeta_M)]$,
where $K$ is the field containing $x$ and $N$ divides $M$.
"""
function kummer_failure(x::nf_elem, M::Int, N::Int)
  @assert divisible(M, N)
  K = parent(x)
  CE = cyclotomic_extension(K, M)
  el = CE.mp[2](x)
  lp = factor(N)
  deg = 1
  for (p, v) in lp
    e = 1
    y = x
    for i = v:-1:1
      fl, y = ispower(y, Int(p), with_roots_unity = true)
      if !fl
        e = v
        break
      end
    end
    deg *= Int(p)^e
  end
  @assert divisible(N, deg)
  return divexact(N, deg)
end

################################################################################
#
#  Reduction of Kummer generator
#
################################################################################

@doc Markdown.doc"""
    reduce_mod_powers(a::nf_elem, n::Int) -> nf_elem
    reduce_mod_powers(a::nf_elem, n::Int, primes::Array{NfOrdIdl, 1}) -> nf_elem

Given some non-zero algebraic integeri $\alpha$, try to find  $\beta$ s.th.
$\beta$ is "small" and $\alpha/\beta$ is an $n$-th power.
If the factorisation of $a$ into prime ideals is known, the ideals
should be passed in.
"""
function reduce_mod_powers(a::nf_elem, n::Int, primes::Array{NfOrdIdl, 1})
  @vprint :ClassField 2 "reducing modulo $(n)-th powers\n"
  @vprint :ClassField 3 "starting with $a\n"
  return reduce_mod_powers(FacElem(a), n, primes)
end

function reduce_mod_powers(a::nf_elem, n::Int)
  return reduce_mod_powers(FacElem(a), n)
end

function reduce_mod_powers(a::FacElem{nf_elem, AnticNumberField}, n::Int, decom::Dict{NfOrdIdl, fmpz})
  a1 = RelSaturate._mod_exponents(a, n)
  if nbits(maximum(values(a.fac))) > 30000
    b = compact_presentation(a1, n)
  else
    c = conjugates_arb_log(a, 64)
    c1 = conjugates_arb_log(a1, 64)
    bn = maximum(fmpz[upper_bound(abs(x), fmpz) for x in c])
    bn1 = maximum(fmpz[upper_bound(abs(x), fmpz) for x in c1])
    if bn1 < root(bn, 2)
      b = compact_presentation(a1, n)
    else
      b = compact_presentation(a, n, decom = decom)
    end
  end
  if any(!iszero(v % n) for (k, v) = b.fac)
    b1 = prod(nf_elem[k^(v % n) for (k, v) = b.fac if !iszero(v % n)])
  else
    b1 = one(base_ring(a))
  end
  d = denominator(b1, maximal_order(parent(b1)))
  k, d1 = ispower(d)
  if k > 1
    d = d1^(div(k, n) + 1)
  end
  b1 *= d^n  #non-optimal, but integral...
  return FacElem(b1)
end

function reduce_mod_powers(a::FacElem{nf_elem, AnticNumberField}, n::Int, primes::Array{NfOrdIdl, 1})
  vals = fmpz[valuation(a, p) for p in primes]
  lp = Dict{NfOrdIdl, fmpz}(primes[i] => vals[i] for i = 1:length(primes) if !iszero(vals[i]))
  return reduce_mod_powers(a, n, lp)
end

function reduce_mod_powers(a::FacElem{nf_elem, AnticNumberField}, n::Int)
  Zk = maximal_order(base_ring(a))
  lp = factor_coprime(a, IdealSet(Zk))
  lp1 = Dict{NfOrdIdl, fmpz}((x, fmpz(y)) for (x, y) in lp)
  return reduce_mod_powers(a, n, lp1)
end