text:
using Phylo

function _newlabel(ids::Vector{Label}) where Label <: Integer
    return isempty(ids) ? 1 : maximum(ids) + 1
end

function _newlabel(names::Vector{String}, prefix)
    names = collect(Iterators.filter(n -> length(n) > length(prefix) &&
                                     n[1:length(prefix)]==prefix, names))
    start = length(names) + 1
    name = prefix * "$start"
    while (name ∈ names)
        start += 1
        name = prefix * "$start"
    end
    return name
end

function _ntrees(::AbstractTree)
    return 1
end

function _extractnode(::T, node::N) where {T <: AbstractTree, N <: AbstractNode}
    return node
end

function _extractnode(tree::T, nodename::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    return getnode(tree, nodename)
end

function _extractnode(::T, pair::Pair{NL, N}) where {NL, BL, T <: AbstractTree{NL, BL},
                                                     N <: AbstractNode}
    return pair[2]
end

function _extractnodename(::T, pair::Pair{NL, N}) where {NL, BL, T <: AbstractTree{NL, BL},
                                                         N <: AbstractNode}
    return pair[1]
end

function _extractbranch(::T, branch::B) where {T <: AbstractTree, B <: AbstractBranch}
    return branch
end

function _extractbranch(tree::T, branchname::BL) where {NL, BL, T <: AbstractTree{NL, BL}}
    return getbranch(tree, branchname)
end

function _extractbranch(::T, pair::Pair{BL, B}) where {NL, BL, T <: AbstractTree{NL, BL},
                                                       B <: AbstractBranch}
    return pair[2]
end

function _extractbranchname(::T, pair::Pair{BL, B}) where {NL, BL, T <: AbstractTree{NL, BL},
                                                           B <: AbstractBranch}
    return pair[1]
end

# AbstractTree methods
"""
    _nodetype(::AbstractTree)

Returns type of nodes in a tree.
"""
function _nodetype end

"""
    _branchtype(::AbstractTree)

Returns type of branches in a tree.
"""
function _branchtype end

"""
    _newbranchlabel(tree::AbstractTree)

"""
function _newbranchlabel end

function _newbranchlabel(tree::AbstractTree{NL, String}) where NL
    return _newlabel(_getbranchnames(tree), "Branch ")
end

function _newbranchlabel(tree::AbstractTree{NL, I}) where {NL, I <: Integer}
    return _newlabel(_getbranchnames(tree))
end

"""
    _addbranch!(tree::AbstractTree, source, destination;
                length::Float64 = NaN,
                branchname = _newbranchlabel(tree))

Must be implemented for any AbstractTree subtype.
"""
function _addbranch! end

"""
    _deletebranch!(tree::AbstractTree, branchname)

Must be implemented for any AbstractTree subtype.
"""
function _deletebranch! end

"""
    _branch!(tree::AbstractTree, source, length::Float64, destination, branchname)

"""
function _branch!(tree::AbstractTree, source, length::Float64, destination, branchname)
    _addbranch!(tree, source, _addnode!(tree, destination), length, branchname)
    return destination
end

"""
    _newnodelabel(tree::AbstractTree)

"""
function _newnodelabel end

function _newnodelabel(tree::T) where {BL, T <: AbstractTree{String, BL}}
    return _newlabel(_getnodenames(tree), "Node ")
end

function _newnodelabel(tree::T) where {I <: Integer, BL, T <: AbstractTree{I, BL}}
    return _newlabel(_getnodenames(tree))
end

"""
    _getnodenames(tree::AbstractTree)

Must be implemented for any AbstractTree subtype.
"""
function _getnodenames end

"""
    _getnodes(tree::AbstractTree)

Must be implemented for any AbstractTree subtype.
"""
function _getnodes end

"""
    _addnode!(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _addnode! end

"""
    _addnodes!(tree::AbstractTree, nodenames::AbstractVector)
    _addnodes!(tree:AbstractTree, count::Integer)

"""
function _addnodes! end

function _addnodes!(tree::AbstractTree, nodenames::AbstractVector)
    return map(name -> _addnode!(tree, name), nodenames)
end

function _addnodes!(tree::AbstractTree, count::Integer)
    return map(name -> addnode!(tree), 1:count)
end

"""
    _deletenode!(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _deletenode! end

"""
    _hasnode(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _hasnode end

"""
    _getnode(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _getnode end

"""
    _getbranchnames(tree::AbstractTree)

Must be implemented for any AbstractTree subtype.
"""
function _getbranchnames end

"""
    _getbranches(tree::AbstractTree)

Must be implemented for any AbstractTree subtype.
"""
function _getbranches end

"""
    _hasbranch(tree::AbstractTree, branchname)

Must be implemented for any AbstractTree subtype.
"""
function _hasbranch end

"""
    _getbranch(tree::AbstractTree, branchname)

Must be implemented for any AbstractTree subtype.
"""
function _getbranch end

"""
    _hasrootheight(::AbstractTree)

"""
function _hasrootheight(::AbstractTree)
    return false
end

"""
    _getrootheight(::AbstractTree)

"""
function _getrootheight(::AbstractTree)
    throw(NullException())
    return NaN
end

"""
    _setrootheight!(::AbstractTree, value)

"""
function _setrootheight!(::AbstractTree, value)
    throw(NullException())
    return value
end
"""
    _validate(::AbstractTree)

"""
function _validate(::AbstractTree)
    return true
end

# AbstractNode methods
"""
    _isleaf(node::AbstractNode)

"""
_isleaf(node::AbstractNode) = _outdegree(node) == 0

"""
    _isroot(node::AbstractNode)

"""
_isroot(node::AbstractNode) = !_hasinbound(node)

"""
    _isinternal(node::AbstractNode)

"""
_isinternal(node::AbstractNode) = _outdegree(node) > 0 && _hasinbound(node)

"""
    _isunattached(node::AbstractNode)

"""
_isunattached(node::AbstractNode) = _outdegree(node) == 0 && !_hasinbound(node)

"""
    _indegree(node::AbstractNode)

"""
_indegree(node::AbstractNode) = _hasinbound(node) ? 1 : 0

"""
    _hasinboundspace(node::AbstractNode)

"""
_hasinboundspace(node::AbstractNode) = !_hasinbound(node)

"""
    _outdegree(node::AbstractNode)

Must be implemented for any AbstractNode subtype.
"""
function _outdegree end

"""
    _hasoutboundspace(node::AbstractNode)

"""
function _hasoutboundspace end

"""
    _hasinbound(node::AbstractNode)

Must be implemented for any AbstractNode subtype.
"""
function _hasinbound end

"""
    _getinbound(node::AbstractNode)

Must be implemented for any AbstractNode subtype.
"""
function _getinbound end

"""
    _setinbound!(node::AbstractNode, inbound)

Must be implemented for any AbstractNode subtype.
"""
function _setinbound! end

"""
    _deleteinbound!(node::AbstractNode, inbound)

Must be implemented for any AbstractNode subtype.
"""
function _deleteinbound! end

"""
    _getoutbounds(node::AbstractNode)

Must be implemented for any AbstractNode subtype.
"""
function _getoutbounds end

"""
    _addinbound!(node::AbstractNode, outbound)

Must be implemented for any AbstractNode subtype.
"""
function _addoutbound! end

"""
    _deleteoutbound!(node::AbstractNode, outbound)

Must be implemented for any AbstractNode subtype.
"""
function _deleteoutbound! end

"""
    _hasheight(tree::AbstractTree, nodename)

"""
function _hasheight(::AbstractTree, _)
    return false
end

"""
    _getheight(tree::AbstractTree, nodename)

"""
function _getheight(::AbstractTree, _)
    throw(NullException())
    return NaN
end

"""
    _setheight!(::AbstractTree, nodename, value)

"""
function _setheight!(::AbstractTree, _, value)
    throw(NullException())
    return value
end

# AbstractBranch methods
"""
    _src

Return source node for a branch. Must be implemented for any
AbstractBranch subtype.
"""
function _src end

"""
    _dst

Return destination node for a branch. Must be implemented for any
AbstractBranch subtype.
"""
function _dst end

"""
    _getlength

Return length of a branch. May be implemented for any AbstractBranch
subtype.
"""
function _getlength end
function _getlength(::AbstractBranch)
    return NaN
end

"""
    _setsrc!

Set source node for a branch. Must be implemented for any
AbstractBranch subtype.
"""
function _setsrc! end

"""
    _setdst!

Set destination node for a graph. Must be implemented for any
AbstractBranch subtype.
"""
function _setdst! end

#  - _getleafnames()
function _getleafnames(tree::AbstractTree)
    return collect(nodenamefilter(_isleaf, tree))
end

function _getleafinfo end
function _setleafinfo! end
function _getnoderecord end
function _setnoderecord! end
function _resetleaves! end
function _clearrootheight! end
function _setnode! end
function _setbranch! end
function _leafinfotype end
function _nleaves end


text:
function ueberpruefeeindeutigeLoesung2x2(A,b)

if size(A,1) == 2 && size(A,2) == 2 && size(b,1) == 2 && size(b,2) == 1

    homogen = nullspace(A);
    # inhomogene Lösung berechnen
    inhomogen = pinv(A)*b;

    if size(homogen,2) > 0
        println("Die Größe eures Gleichungssystems stimmt zwar, allerdings lässt es sich nicht eindeutig lösen.")
    else
        println("Super, ihr habt ein lineares Gleichungssystem aufgestellt, das sich eindeutig lösen lässt!");
        y1 = inv(transpose(A)*A)*transpose(A)*b;
        println("Eure Lösung lautet:");
        println("x1 = ",round(y1[1];digits=2));
        println("x2 = ",round(y1[2];digits=2));
    end
else
    println("Ihr habt noch einen Fehler bei der Anzahl der Einträge in der Matrix oder dem Vektor. Stellt noch einmal sicher, dass ihr ein lineares Gleichungssystem mit einer 2x2 Matrix augestellt habt.")
end

end

ext: jl


text:
closesttozero(ev) = ev[sortperm(ev, by = abs∘real)]
getinterval(a,b) = (min(a,b), max(a,b))
####################################################################################################
function displayIteration(i, residual, itlinear = 0, lastRow = false)
	if lastRow
		lastRow && println("└─────────────┴──────────────────────┴────────────────┘")
	else
		if i == 0
			println("\n┌─────────────────────────────────────────────────────┐")
			  println("│ Newton Iterations      f(x)      Linear Iterations  │")
			  println("├─────────────┬──────────────────────┬────────────────┤")
		end
		if length(itlinear) == 1
			@printf("│%8d     │ %16.4e     │ %8d       │\n", i, residual, itlinear);
		else
			if itlinear isa Tuple{Int64, Int64}
				@printf("|%8d     │ %16.4e     │ (%4d, %4d)   |\n", i, residual, itlinear[1], itlinear[2]);
			else
				# used for nested linear solves
				@printf("%11d %19.4e  ", i, residual); println(itlinear);
			end
		end
	end
end
####################################################################################################
function computeEigenvalues(it::ContIterable, u0, par, nev = it.contParams.nev; kwargs...)
	return it.contParams.newtonOptions.eigsolver(it.J(u0, par), nev; kwargs...)
end

function computeEigenvalues(iter::ContIterable, state::ContState; kwargs...)
	# we compute the eigen-elements
	n = state.n_unstable[2]
	nev_ = max(n + 5, iter.contParams.nev)
	eiginfo = computeEigenvalues(iter, getx(state), setParam(iter, getp(state)), nev_; kwargs...)
	_isstable, n_unstable, n_imag = isStable(iter.contParams, eiginfo[1])
	return eiginfo, _isstable, n_unstable, n_imag
end

# same as previous but we save the eigen-elements in state
function computeEigenvalues!(iter::ContIterable, state::ContState; saveEigenVec = true, kwargs...)
	eiginfo, _isstable, n_unstable, n_imag = computeEigenvalues(iter, state; kwargs...)
	# we update the state
	updateStability!(state, n_unstable, n_imag)
	if isnothing(state.eigvals) == false
		state.eigvals = eiginfo[1]
	end
	if saveEigenVec && iter.contParams.saveEigenvectors
		state.eigvecs = eiginfo[2]
	end
	# iteration number in eigen-solver
	it_number = eiginfo[end]
	return it_number
end

####################################################################################################
"""
	finiteDifferences(F, x::AbstractVector; δ = 1e-9)

Compute a Jacobian by Finite Differences
"""
function finiteDifferences(F, x::AbstractVector; δ = 1e-9)
	N = length(x)
	Nf = length(F(x))
	J = zeros(eltype(x), Nf, N)
	finiteDifferences!(F, J, x; δ = δ)
	return J
end

"""
	finiteDifferences!(F, J, x::AbstractVector; δ = 1e-9)

Compute a Jacobian by Finite Differences, update J. Use the centered formula (f(x+δ)-f(x-δ))/2δ
"""
function finiteDifferences!(F, J, x::AbstractVector; δ = 1e-9)
	f = F(x)
	x1 = copy(x)
	@inbounds for i ∈ eachindex(x)
		x1[i] += δ
		J[:, i] .= F(x1)
		x1[i] -= 2δ
		J[:, i] .-= F(x1)
		J[:, i] ./= 2δ
		x1[i] += δ
	end
	return J
end
####################################################################################################
using BlockArrays, SparseArrays
"""
Function waiting to be accepted to BlockArrays.jl
"""
function blockToSparse(J::AbstractBlockArray)
	nl, nc = size(J.blocks)
	# form the first line of blocks
	res = J[Block(1,1)]
	@inbounds for j in 2:nc
		res = hcat(res, J[Block(1,j)])
	end
	# continue with the other lines
	@inbounds for i in 2:nl
		line = J[Block(i,1)]
		for j in 2:nc
			line = hcat(line, J[Block(i,j)])
		end
		res = vcat(res,line)
	end
	return res
end
####################################################################################################
"""
$(SIGNATURES)

This function extracts the indices of the blocks composing the matrix A which is a M x M Block matrix where each block N x N has the same sparsity.
"""
function getBlocks(A::SparseMatrixCSC, N, M)
	I,J,K = findnz(A)
	out = [Vector{Int}() for i in 1:M+1, j in 1:M+1];
	for k in eachindex(I)
		m, l = div(I[k]-1, N), div(J[k]-1, N)
		push!(out[1+m,1+l], k)
	end
	res = [length(m) for m in out]
	out
end
####################################################################################################
"""
$(SIGNATURES)

This function implements a counter. If `everyN == 0`, it returns false. Otherwise, it returns `true` when `step` is a multiple of `everyN`
"""
function modCounter(step, everyN)
	if step == 0; return false; end
	if everyN == 0; return false; end
	if everyN == 1; return true; end
	return mod(step, everyN) == 0
end
####################################################################################################
# this trick is extracted from KrylovKit. It allows for the Jacobian to be specified as a matrix (sparse / dense) or as a function.
apply(A::AbstractMatrix, x::AbstractVector) = A * x
apply(f, x) = f(x)

# apply!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
# apply!(y, f, x) = copyto!(y, f(x))

# empty eigenvectors to save memory
# _empty(a::AbstractVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()
# _empty(a::AbstractMatrix{T}, ::Type{U}=T) where {T,U} = similar(a, (0,0))
####################################################################################################
# the following structs are a machinery to extend multilinear mapping from Real valued to Complex valued Arrays
# this is done so as to use AD (ForwardDiff.jl,...) to provide the differentials which only works on reals (usually).

"""
$(TYPEDEF)

This structure wraps a bilinear map to allow evaluation on Complex arguments. This is especially useful when these maps are produced by ForwardDiff.jl.
"""
struct BilinearMap{Tm}
	bl::Tm
end

function (R2::BilinearMap)(dx1, dx2)
	dx1r = real.(dx1); dx2r = real.(dx2)
	dx1i = imag.(dx1); dx2i = imag.(dx2)
	return R2(dx1r, dx2r) .- R2(dx1i, dx2i) .+ im .* (R2(dx1r, dx2i) .+ R2(dx1i, dx2r))
end
(b::BilinearMap)(dx1::T, dx2::T) where {T <: AbstractArray{<: Real}} = b.bl(dx1, dx2)

"""
$(TYPEDEF)

This structure wraps a trilinear map to allow evaluation on Complex arguments. This is especially useful when these maps are produced by ForwardDiff.jl.
"""
struct TrilinearMap{Tm}
	tl::Tm
end

function (R3::TrilinearMap)(dx1, dx2, dx3)
	dx1r = real.(dx1); dx2r = real.(dx2); dx3r = real.(dx3)
	dx1i = imag.(dx1); dx2i = imag.(dx2); dx3i = imag.(dx3)
	outr =  R3(dx1r, dx2r, dx3r) .- R3(dx1r, dx2i, dx3i) .-
			R3(dx1i, dx2r, dx3i) .- R3(dx1i, dx2i, dx3r)
	outi =  R3(dx1r, dx2r, dx3i) .+ R3(dx1r, dx2i, dx3r) .+
			R3(dx1i, dx2r, dx3r) .- R3(dx1i, dx2i, dx3i)
	return Complex.(outr, outi)
end

(b::TrilinearMap)(dx1::T, dx2::T, dx3::T) where {T <: AbstractArray{<: Real}} = b.tl(dx1, dx2, dx3)
####################################################################################################
"""
$(SIGNATURES)

Return the 3-Jet of f, namely (f, df, d2f, d3f). If J is passed, the function returns (f, J, d2f, d3f). These differentials are computed using ForwardDiff.jl
"""
function get3Jet(f, J = nothing)
	d1f(x,p,dx1) = ForwardDiff.derivative(t -> f(x .+ t .* dx1, p), 0.)
	d2f(x,p,dx1,dx2) = ForwardDiff.derivative(t -> d1f(x .+ t .* dx2, p, dx1), 0.)
	d3f(x,p,dx1,dx2,dx3) = ForwardDiff.derivative(t -> d2f(x .+ t .* dx3, p, dx1, dx2), 0.)
	if isnothing(J)
	 	return (f, d1f, d2f, d3f)
	else
		return (f, J, d2f, d3f)
	end
end


text:
module SumOfSquaresBench

using DynamicPolynomials

import JuMP

import SumOfSquares
include(joinpath(dirname(dirname(pathof(SumOfSquares))), "test", "Tests", "Tests.jl"))

function dummy_problem(opt; T = Float64)
    T == Float64 || return
    @polyvar x y
    f = x^2 + 2x*y + y^2
    model = JuMP.Model(opt)
    # All solvers should be able to solve this, even LP solvers.
    JuMP.@constraint(model, f in SumOfSquares.DSOSCone())
    JuMP.optimize!(model)
end

const PROBLEMS = Dict(
    "lp" => Tests.linear_tests,
    "socp" => Tests.soc_tests,
    "sdp" => Tests.sd_tests
)

handle_problem_function(opt) = opt

function _test(test_func::Function, opt, ::Val{true}, atol, rtol, T)
    T == Float64 || return
    config = JuMP.MOI.Test.TestConfig(atol=atol, rtol=rtol, query=false)
    test_func(opt, config)
end

function foreach_problem(apply::Function, class, problems; exclude = exclude)
    for (name, test_func) in PROBLEMS[class]
        apply(name, (args...) -> _test(test_func, args...))
    end
end

end

#function do_tests(opt, problems=nothing; exclude = Regex[], T = Float64,
#                  atol=1e-3, rtol=1e-3, to)
#    for (class, tests) in [
#        ("lp", Tests.linear_tests),
#        ("socp", Tests.soc_tests),
#        ("sdp", Tests.sd_tests)
#    ]
#        any(occursin.(exclude, Ref(class))) && continue
#        @timeit to "$class" begin
#            @testset "$class" begin
#                for (name, problem_func) in tests
#                    @timeit to "$name" begin
#                        @testset "$name" begin
#                            problem_func(handle_problem!, Val(true), atol, rtol, T)
#                        end
#                    end
#                end
#            end
#        end
#    end
#end


text:
@testset "PointwiseLearn" begin
  # synthetic data
  Random.seed!(1234)
  f(x,y) = sin(4*(abs(x)+abs(y))) < 0 ? 1 : 0
  X = [sin(i/10) for i in 1:100, j in 1:100]
  Y = [sin(j/10) for i in 1:100, j in 1:100]
  Z = categorical(f.(X,Y))
  ϵ₁ = 0.1randn(Float64, size(X))
  ϵ₂ = 0.1randn(Float64, size(Y))

  # source and target data
  S = georef((X=X,Y=Y,Z=Z))
  T = georef((X=X+ϵ₁,Y=Y+ϵ₂))

  # view versions
  inds = shuffle(1:nelements(S))
  Sv = view(S, inds)
  Tv = view(T, inds)

  # classification task
  𝓉 = ClassificationTask((:X,:Y), :Z)

  # learning problems
  𝒫₁ = LearningProblem(S, T, 𝓉)
  𝒫₂ = LearningProblem(Sv, Tv, 𝓉)

  # pointwise solver
  ℒ = PointwiseLearn(dtree())

  R₁ = solve(𝒫₁, ℒ)
  R₂ = solve(𝒫₂, ℒ)

  # error is small
  @test mean(S[:Z] .!= R₁[:Z]) < 0.15
  @test mean(Sv[:Z] .!= R₂[:Z]) < 0.15

  if visualtests
    for (i,s) in enumerate([(S,R₁), (Sv,R₂)])
      p1 = plot(s[1], (:Z,))
      p2 = plot(s[2])
      plt = plot(p1, p2, size=(800,400))
      @test_reference "data/pointwise$i.png" plt
    end
  end
end


text:
### Calculate the number of passengers on each flight
NMR.clean_split_dir()
jobid = 3
NMR.runjob(NMR.Nmr(jobid, ["AComp_Passenger_data.csv"], NMR.mapper_parserecordacomp, "acomp.csv"))
NMR.runjob(NMR.Nmr(jobid, ["acomp.csv"], NMR.reducer_numberofpassengers, "acomp_reduced.json"))
output = NMR.runcombiner("acomp_reduced.json", NMR.combiner_parsejson)
output = map(x -> string(JSON.parse(x)[1]) * "," * string(JSON.parse(x)[2]), output)

# store output
header = "flightid,passengers"
NMR.save_job_output(jobid, header, output)


text:
module NRELssc

hello() = print("Hello World adfa!")

ENV["SSC_DIR"]="e:\\develop\\YangSimOpt.jl\\external\\nrelssc\\"

const ssc_location = ENV["SSC_DIR"]

const sscdll = joinpath(ssc_location, "ssc.dll")

const INVALID = 0
const STRING = 1
const NUMBER = 2
const ARRAY = 3
const MATRIX = 4
const INPUT = 1
const OUTPUT = 2
const INOUT = 3

function ssc_version()
    ccall((:ssc_version, sscdll), Cint, ())
end

function ssc_data_create()
    ccall((:ssc_data_create, sscdll), Ptr{Cvoid}, ())
end

function ssc_build_info()
    build_info = ccall((:ssc_build_info, sscdll), Cstring, ())
    if build_info == C_NULL
        error("call ssc_build_info error.")
    end
    return unsafe_string(build_info)
end

end # module


text:
using BinDeps
using Compat

if isdefined((@static VERSION < v"0.7.0-DEV.484" ? current_module() : @__MODULE__), :Compat)
    import Compat.Libdl
elseif VERSION >= v"0.7.0-DEV.3382"
    import Libdl
end

@BinDeps.setup

libopus = library_dependency("libopus", aliases = ["libopus"])

@static if Sys.isapple()
  using Homebrew
  provides( Homebrew.HB, "opus", libopus, os = :Darwin )
end

provides( AptGet, "libopus0", libopus )

# Source build
provides( Sources,
          URI("http://downloads.xiph.org/releases/opus/opus-1.1.2.tar.gz"),
          SHA="0e290078e31211baa7b5886bcc8ab6bc048b9fc83882532da4a1a45e58e907fd",
          libopus )
provides( BuildProcess,
          Autotools(libtarget=".libs/libopus."*Libdl.dlext, configure_options = ["--libdir=$(BinDeps.libdir(libopus))"]),
          libopus )

@compat @BinDeps.install Dict(:libopus => :libopus)


text:
module Symbolics

using DocStringExtensions

using LinearAlgebra

using Reexport

@reexport using SymbolicUtils

import SymbolicUtils: Term, Add, Mul, Pow, Sym, symtype,
                      FnType, @rule, Rewriters, substitute, similarterm,
                      promote_symtype, istree, operation, arguments

import SymbolicUtils.Rewriters: Chain, Prewalk, Postwalk, Fixpoint

import SymbolicUtils.Code: toexpr

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

# re-export

export simplify, substitute

using SciMLBase, IfElse
export Num
include("num.jl")

export Equation, ConstrainedEquation
include("equations.jl")

include("utils.jl")
export degree

using MacroTools
import MacroTools: splitdef, combinedef, postwalk, striplines
export @register
include("register.jl")

using TreeViews
export @variables, Variable
include("variable.jl")

include("linearity.jl")

using DiffRules, SpecialFunctions, NaNMath

using SparseArrays

export Differential, expand_derivatives

include("diff.jl")

include("linear_algebra.jl")

import Libdl
include("build_function.jl")
export build_function

using Distributions
include("extra_functions.jl")

using Latexify
include("latexify_recipes.jl")

using RecipesBase
include("plot_recipes.jl")

include("domains.jl")

using Requires

export symbolics_to_sympy
include("init.jl")

end # module


text:
"""
    test_init_field(N, L, ftf)

Test that the field initialized by the field type function `ftf` on the grid g
has the correct size.
"""
function correct_field_size(arch::Architecture, g::Grid, field_type)
    f = field_type(arch, g)
    size(f) == size(g)
end

"""
    test_set_field(N, L, ftf, val)

Test that the field initialized by the field type function `ftf` on the grid g
can be correctly filled with the value `val` using the `set!(f::Field, v)`
function.
"""
function correct_field_value_was_set(arch::Architecture, g::Grid, field_type, val::Number)
    f = field_type(arch, g)
    set!(f, val)
    f.data ≈ val * ones(size(f))
end

function correct_field_addition(arch::Architecture, g::Grid, field_type, val1::Number, val2::Number)
    f1 = field_type(arch, g)
    f2 = field_type(arch, g)

    set!(f1, val1)
    set!(f2, val2)
    f3 = f1 + f2

    val3 = convert(mm.float_type, val1) + convert(mm.float_type, val2)
    f_ans = val3 * ones(size(f1))
    f3.data ≈ f_ans
end


text:
"""
    kerbal(c::KRPCConnection, call::Request{S, P, R}) where {S, P, R}

Send a single message `call` to the server with connection `c`. See the
call stubs generated in the KRPC.Interface.[service] module for valid requests.

# Examples

```
active_vessel = kerbal(conn, KRPC.Interface.SpaceCenter.get_ActiveVessel())
```
"""
function kerbal(c::KRPCConnection, call::Request{S, P, R}) where {S, P, R}
	return kerbal(c, (call, ))[1]
end

@generated function make_kerbal_call(call::RT) where {S, P, R, RT<:Request{S, P, R}}
	argument_asts = Any[]
	for (i, arg) in enumerate(fieldnames(RT))
		push!(argument_asts, :(krpc.schema.Argument(position=$(i-1), value=getWireValue(getfield(call, $(QuoteNode(arg)))))))
	end
	return :(krpc.schema.ProcedureCall(service=string(S), procedure=string(P), arguments=[$(argument_asts...)]))
end

process_return(conn, value, call::Request{S, P, Nothing}) where {S, P} = Nothing()
process_return(conn, value, call::Request{S, P, R}) where {S, P, R} = getJuliaValue(conn, value, R)

"""
    kerbal(c::KRPCConnection, calls::T) where {K, T<:Tuple{Vararg{RT where {S, P, R, RT<:Request{S, P, R}}, K}}}

Send multiple messages `calls` to the server with connection `c`. See the
call stubs generated in the KRPC.Interface.[service] module for valid requests.

# Examples

```
active_vessel, gamemode = kerbal(conn, (KRPC.Interface.SpaceCenter.get_ActiveVessel(), KRPC.Interface.SpaceCenter.get_GameMode()))
```
"""
@generated function kerbal(c::KRPCConnection, calls::T) where {K, T<:Tuple{Vararg{RT where {S, P, R, RT<:Request{S, P, R}}, K}}}
	call_asts = [:(make_kerbal_call(calls[$i])) for i=1:K]
	return_asts = [:(process_return(c,res.results[$i].value, calls[$i])) for i=1:K]
    return quote
    	res = SendBiMessage(c, krpc.schema.Request(calls=[$(call_asts...)]))
    	handle_potential_error(res)
    	return ($(return_asts...), )
    end
end


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use Knitro to solve the following
# simple nonlinear optimization problem.  This model is test problem
# HS40 from the Hock & Schittkowski collection.
#
# max   x0*x1*x2*x3        (obj)
# s.t.  x0^3 + x1^2 = 1    (c0)
#       x0^2*x3 - x2 = 0   (c1)
#       x3^2 - x1 = 0      (c2)
#
# This is the same as exampleNLP2.c, but it demonstrates using
# multiple callbacks for the nonlinear evaluations and computing
# some gradients using finite-differencs, while others are provided
# in callback routines.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO, Test

#*------------------------------------------------------------------*
#*     FUNCTION EVALUATION CALLBACKS                                *
#*------------------------------------------------------------------*

# The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
# Only "obj" is set in the KNITRO.KN_eval_result structure.
function callbackEvalObj(kc, cb, evalRequest, evalResult, userParams)
    xind = userParams[:data]

    if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALFC
        println("*** callbackEvalObj incorrectly called with eval type ", evalRequest.evalRequestCode)
        return -1
    end
    x = evalRequest.x

    # Evaluate nonlinear term in objective
    evalResult.obj[1] = x[xind[1]] * x[xind[2]] * x[xind[3]] * x[xind[4]]

    return 0
end

# The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
# Only "c0" is set in the KNITRO.KN_eval_result structure.
function callbackEvalC0(kc, cb, evalRequest, evalResult, userParams)
    xind = userParams[:data]

    if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALFC
        println("*** callbackEvalC0 incorrectly called with eval type ", evalRequest.evalRequestCode)
        return -1
    end
    x = evalRequest.x

    # Evaluate nonlinear terms in constraint, c0
    evalResult.c[1] = x[xind[1]] * x[xind[1]] * x[xind[1]]

    return 0
end

# The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
# Only "c1" is set in the KNITRO.KN_eval_result structure.
function callbackEvalC1(kc, cb, evalRequest, evalResult, userParams)
    xind = userParams[:data]

    if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALFC
        println("*** callbackEvalC1 incorrectly called with eval type %d" % evalRequest.evalRequestCode)
        return -1
    end
    x = evalRequest.x

    # Evaluate nonlinear terms in constraint, c1
    evalResult.c[1] = x[xind[1]] * x[xind[1]] * x[xind[4]]

    return 0
end

#*------------------------------------------------------------------*
#*     GRADIENT EVALUATION CALLBACKS                                *
#*------------------------------------------------------------------*

# The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
# Only "objGrad" is set in the KNITRO.KN_eval_result structure.
function callbackEvalObjGrad(kc, cb, evalRequest, evalResult, userParams)
    xind = userParams[:data]

    if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALGA
        println("*** callbackEvalObjGrad incorrectly called with eval type %d" % evalRequest.evalRequestCode)
        return -1
    end
    x = evalRequest.x

    # Evaluate nonlinear terms in objective gradient
    evalResult.objGrad[xind[1]] = x[xind[2]] * x[xind[3]] * x[xind[4]]
    evalResult.objGrad[xind[2]] = x[xind[1]] * x[xind[3]] * x[xind[4]]
    evalResult.objGrad[xind[3]] = x[xind[1]] * x[xind[2]] * x[xind[4]]
    evalResult.objGrad[xind[4]] = x[xind[1]] * x[xind[2]] * x[xind[3]]

    return 0
end

# The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
# Only gradient of c0 is set in the KNITRO.KN_eval_result structure.
function callbackEvalC0Grad(kc, cb, evalRequest, evalResult, userParams)
    xind = userParams[:data]

    if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALGA
        println("*** callbackEvalC0Grad incorrectly called with eval type ", evalRequest.evalRequestCode)
        return -1
    end
    x = evalRequest.x

    # Evaluate nonlinear terms in c0 constraint gradients
    evalResult.jac[1] = 3.0 * x[xind[1]] * x[xind[1]] # *  derivative of x0^3 term  wrt x0

    return 0
end

#*------------------------------------------------------------------*
#*     main                                                         *
#*------------------------------------------------------------------*

# Create a new Knitro solver instance.
kc = KNITRO.KN_new()

# Initialize Knitro with the problem definition.

# Add the variables and specify initial values for them.
# Note: any unset lower bounds are assumed to be
# unbounded below and any unset upper bounds are
# assumed to be unbounded above.
xIndices = KNITRO.KN_add_vars(kc, 4)
for x in xIndices
    KNITRO.KN_set_var_primal_init_values(kc, x, 0.8)
end

# Add the constraints and set the rhs and coefficients
cIndices = KNITRO.KN_add_cons(kc, 3)
KNITRO.KN_set_con_eqbnds(kc, cIndices[1], 1.0)
KNITRO.KN_set_con_eqbnds(kc, cIndices[2], 0.0)
KNITRO.KN_set_con_eqbnds(kc, cIndices[3], 0.0)

# Coefficients for 2 linear terms
lconIndexCons = Int32[1, 2]
lconIndexVars = Int32[2, 1]
lconCoefs = [-1.0, -1.0]
KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

# Coefficients for 2 quadratic terms

#* x1^2 term in c0
#* x3^2 term in c2
qconIndexCons = Int32[0, 2]
qconIndexVars1 = Int32[1, 3]
qconIndexVars2 = Int32[1, 3]
qconCoefs = [1.0, 1.0]

KNITRO.KN_add_con_quadratic_struct(kc, qconIndexCons, qconIndexVars1, qconIndexVars2, qconCoefs)

# Add separate callbacks.

# Set callback data for nonlinear objective term.
cbObj = KNITRO.KN_add_eval_callback(kc, callbackEvalObj)
KNITRO.KN_set_cb_grad(kc, cbObj, callbackEvalObjGrad,
                      objGradIndexVars=xIndices, nV=length(xIndices))

# Set callback data for nonlinear constraint 0 term.
cbC0 = KNITRO.KN_add_eval_callback(kc, false, [cIndices[1]], callbackEvalC0)
indexCons = cIndices[1]  # constraint c0
indexVars = xIndices[1]  # variable x0
KNITRO.KN_set_cb_grad(kc, cbC0, callbackEvalC0Grad, jacIndexCons=[indexCons],
                      jacIndexVars=[indexVars])

# Set callback data for nonlinear constraint 1 term
cbC1 = KNITRO.KN_add_eval_callback(kc, false, [cIndices[2]], callbackEvalC1)
indexCons = [cIndices[2], cIndices[2]]  # constraint c1
indexVars = [xIndices[1], xIndices[1]]  # variables x0 and x3
KNITRO.KN_set_cb_grad(kc, cbC1, nothing, jacIndexCons=indexCons, jacIndexVars=indexVars)
# This one will be approximated via forward finite differences.
KNITRO.KN_set_cb_gradopt(kc, cbC1, KNITRO.KN_GRADOPT_FORWARD)

# Demonstrate passing a userParams structure to the evaluation
# callbacks.  Here we pass back the variable indices set from
# KNITRO.KN_add_vars() for use in the callbacks. Here we pass the same
# userParams structure to each callback but we could define different
# userParams for different callbacks.  This could be useful, for
# instance, if different callbacks operate on different sets of
# variables.
# We increment the indices as Julia is 1-indexed
userEval = xIndices .+ 1
KNITRO.KN_set_cb_user_params(kc, cbObj, userEval)
KNITRO.KN_set_cb_user_params(kc, cbC0, userEval)
KNITRO.KN_set_cb_user_params(kc, cbC1, userEval)

# Set minimize or maximize(if not set, assumed minimize)
KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)

# Approximate hessian using BFGS
KNITRO.KN_set_param(kc, "hessopt", KNITRO.KN_HESSOPT_BFGS)

# Solve the problem.
#
# Return status codes are defined in "knitro.h" and described
# in the Knitro manual.
nStatus = KNITRO.KN_solve(kc)

println()
println("Knitro converged with final status = ", nStatus)

# An example of obtaining solution information.
nStatus, objSol, x, lambda_ =  KNITRO.KN_get_solution(kc)
println("  optimal objective value  = ", objSol)
println("  optimal primal values x  = ", x)
println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))

# Delete the Knitro solver instance.
KNITRO.KN_free(kc)

@testset "Exemple multipleCB" begin
    @test nStatus == 0
    @test_broken objSol ≈ 0.25
    @test_broken x ≈ [0.793701, 0.707107, 0.529732, 0.840896] atol=1e-5
end


text:
module HCIToolbox

export derotate,
       derotate!,
       expand,
       flatten,
       collapse,
       collapse!,
       crop,
       cropview,
       mask_circle!,
       mask_circle,
       mask_annulus!,
       mask_annulus,
       normalize_par_angles,
       normalize_par_angles!,
       inject,
       inject!,
       scale,
       invscale,
       shift_frame,
       shift_frame!,
       scale_and_stack,
       invscale_and_collapse,
       scale_list,
       # geometry
       AnnulusView,
       MultiAnnulusView,
       eachannulus,
       inverse!,
       inverse,
       # profiles
       radial_profile

# Utilities for dealing with cubes like derotating and median combining
include("morphology.jl")
include("geometry/geometry.jl")
include("inject.jl")

# # Utilities for mask design
include("masking.jl")
include("angles.jl")

include("scaling.jl")
include("profiles.jl")

end


text:
function getAhornModDirs()
    if !get(config, "load_plugins_ahorn", true)
        return String[]
    end

    modsPath = joinpath(storageDirectory, "plugins")
    targetFolders = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if isdir(joinpath(modsPath, fn))
                push!(targetFolders, joinpath(modsPath, fn))
            end
        end
    end

    return targetFolders
end

function getAhornModZips()
    if !get(config, "load_plugins_ahorn_zip", true)
        return String[]
    end

    modsPath = joinpath(storageDirectory, "plugins")
    targetFolders = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if isfile(joinpath(modsPath, fn)) && hasExt(fn, ".zip")
                push!(targetZips, joinpath(modsPath, fn))
            end
        end
    end

    return targetFolders
end

function getEverestBlacklist()
    entries = Set()

    celesteDir = get(config, "celeste_dir", "")
    blacklistPath = joinpath(celesteDir, "Mods", "blacklist.txt")

    if isfile(blacklistPath)
        for line in readlines(blacklistPath)
            # Lines starting by # are comments
            if !startswith(line, "#")
                push!(entries, strip(line))
            end
        end
    end

    return entries
end

# Good enough for now, can be expanded later
function getModLoadingFilter()
    useEverestBlacklist = get(config, "use_everest_blacklist", true)

    if useEverestBlacklist
        return getEverestBlacklist()
    end

    return Set()
end

function getCelesteModDirs(ignoreFilters::Bool=false)
    if !get(config, "load_plugins_celeste", true)
        return String[]
    end

    celesteDir = get(config, "celeste_dir", "")
    modsPath = joinpath(celesteDir, "Mods")
    modFilters = getModLoadingFilter()

    targetFolders = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if !ignoreFilters && fn in modFilters
                continue
            end

            if isdir(joinpath(modsPath, fn))
                push!(targetFolders, joinpath(modsPath, fn))
            end
        end
    end

    return targetFolders
end

function samePath(path1::String, path2::String, ignoreCase::Bool=Sys.iswindows())
    if ignoreCase
        return lowercase(path1) == lowercase(path2)

    else
        return path1 == path2
    end
end

function getModRoot(fn::String)
    celesteDir = get(config, "celeste_dir", "")
    modsPath = normpath(joinpath(celesteDir, "Mods"))

    path = normpath(fn)

    while true
        parent = dirname(path)

        if samePath(parent, modsPath)
            if isdir(path)
                return true, path

            else
                return true, parent
            end
        end

        if samePath(path, parent)
            return false, ""
        end

        path = parent
    end
end

function getCelesteModZips(ignoreFilters::Bool=false)
    if !get(config, "load_plugins_celeste_zip", true)
        return String[]
    end

    celesteDir = get(config, "celeste_dir", "")
    modsPath = joinpath(celesteDir, "Mods")
    modFilters = getModLoadingFilter()

    targetZips = String[]

    if isdir(modsPath)
        for fn in readdir(modsPath)
            if !ignoreFilters && fn in modFilters
                continue
            end

            if isfile(joinpath(modsPath, fn)) && hasExt(fn, ".zip")
                push!(targetZips, joinpath(modsPath, fn))
            end
        end
    end

    return targetZips
end

function findExternalSprite(resource::String, atlas::String="Gameplay")
    targetFolders = getCelesteModDirs()
    targetZips = getCelesteModZips()
    atlasPath = joinpath("Graphics", "Atlases", atlas)

    for target in targetFolders
        fn = joinpath(target, atlasPath, resource * ".png")
        if isfile(fn)
            return fn
        end
    end

    # ZipFile uses unix paths
    gameplayUnixPath = replace(atlasPath, "\\" => "/")

    for target in targetZips
        fn = gameplayUnixPath * "/" * resource * ".png"
        zipfh = ZipFile.Reader(target)

        for file in zipfh.files
            if file.name == fn
                return target
            end
        end

        close(zipfh)
    end
end

const warnedFilenames = Set{String}()

function warnBadExtensions(fixable::Array{String, 1}, nonFixable::Array{String, 1}, correctExt::String)
    res = false

    fixable = String[fn for fn in fixable if !(fn in warnedFilenames)]
    nonFixable = String[fn for fn in nonFixable if !(fn in warnedFilenames)]

    if !isempty(fixable)
        fixableCount = length(fixable)
        plural = fixableCount == 1 ? "" : "s"
        dialogText = "You have $fixableCount image$plural with bad extension$plural.\n" *
            "These are not loadable by Celeste (Everest).\n" *
            "Do you want to rename them?\n\n" *
            join(fixable, "\n")
        confirmed = ask_dialog(dialogText, Ahorn.window)
        res = confirmed

        if confirmed
            for fn in fixable
                path, ext = splitext(fn)
                out = path * correctExt
                tmp = out * ".tmp"

                # Workaround for Windows...
                mv(fn, tmp)
                mv(tmp, out)
            end

        else
            push!(warnedFilenames, fixable...)
        end
    end

    if !isempty(nonFixable)
        nonFixableCount = length(nonFixable)
        plural = nonFixableCount == 1 ? "" : "s"
        dialogText = "You have $nonFixableCount file$plural with bad extension$plural that cannot be automatically fixed.\n" *
            "These are not loadable by Celeste (Everest) or Ahorn.\n" *
            "Please rename the file$plural manually with the proper extension '$correctExt' or contact the mod maker.\n\n" *
            join(nonFixable, "\n")
        info_dialog(dialogText, Ahorn.window)

        push!(warnedFilenames, nonFixable...)
    end

    return res
end

const resourceZipCache = Dict{String, Any}()

function getSpriteAtlasAndResource(filename::String)
    parts = split(filename, "/")
    atlas = parts[3]
    resource = join(parts[3:end], "/")

    return atlas, resource
end

# Resource prefix filters uneeded files, uses unix paths
function cacheZipContent(filename::String, resourcePrefix::String="Graphics/Atlases/Gameplay/", ahornPrefix::String="Ahorn/", languagePrefix::String="Ahorn/lang/")
    loadMetaYaml = get(config, "load_image_meta_yaml", false)

    debug.log("Caching zip content for $filename", "ZIP_CACHING_VERBOSE")

    resourceZipCache[filename] = Dict{String, Any}()

    # Always uses Unix path
    try
        zipfh = ZipFile.Reader(filename)

        for file in zipfh.files
            name = file.name

            if startswith(name, resourcePrefix)
                if hasExt(name, ".png")
                    atlas, resource = getSpriteAtlasAndResource(name)
                    resourceZipCache[filename][resource] = Cairo.read_from_png(fastZipFileRead(file, IOBuffer))

                elseif loadMetaYaml && hasExt(name, ".yaml")
                    try
                        atlas, resource = getSpriteAtlasAndResource(name)
                        resourceZipCache[filename][resource * ".meta.yaml"] = YAML.load(fastZipFileRead(file, String))

                    catch
                        # Bad yaml file
                    end
                end

            elseif startswith(name, ahornPrefix) && hasExt(name, ".jl")
                resourceZipCache[filename][name] = fastZipFileRead(file, String)

            elseif startswith(name, languagePrefix) && hasExt(name, ".lang")
                resourceZipCache[filename][name] = fastZipFileRead(file, String)
            end
        end

        close(zipfh)

    catch e
        println(Base.stderr, "Failed to get content of zip file: $filename")
        for (exc, bt) in Base.catch_stack()
            showerror(Base.stderr, exc, bt)
            println(Base.stderr, "")
        end
    end
end

# Cleanup and destroy surfaces
function uncacheZipContent(filename::String)
    debug.log("Uncaching zip content for $filename", "ZIP_CACHING_VERBOSE")

    fileCache = get(resourceZipCache, filename, nothing)

    if fileCache !== nothing
        for (resource, data) in fileCache
            if hasExt(resource, ".png")
                deleteSurface(data)
            end
        end
    end
end

function findExternalSpritesInZip(filename::String, atlasesPath::String, atlasesUnixPath::String, badExt::Array{String, 1}, allowRetry::Bool=true)
    res = Tuple{String, String, String}[]
    fileCache = get(resourceZipCache, filename, nothing)

    if fileCache !== nothing
        for (resource, data) in fileCache
            if hasExt(resource, ".png")
                if !hasExt(resource, ".png", false)
                    # Case sensitive check for warnings
                    push!(badExt, resource)
                end

                parts = split(resource, "/")
                atlas = parts[1]
                resourceNoAtlas = join(parts[2:end], "/")

                push!(res, (resourceNoAtlas, filename, atlas))
            end
        end

    else
        if allowRetry
            cacheZipContent(filename)

            return findExternalSpritesInZip(filename, atlasesPath, atlasesUnixPath, badExt, false)
        end
    end

    return res
end

function findExternalSpritesInFolder(filename::String, atlasesPath::String, badExt::Array{String, 1}, targetAtlases::Array{String, 1}=["Gameplay"])
    res = Tuple{String, String, String}[]
    atlasPath = joinpath(filename, atlasesPath)

    if !isdir(atlasPath)
        return res
    end

    for atlas in readdir(atlasPath)
        if !(atlas in targetAtlases)
            continue
        end

        path = joinpath(atlasPath, atlas)

        if isdir(path)
            for (root, dir, files) in walkdir(path)
                for file in files
                    if hasExt(file, ".png")
                        if !hasExt(file, ".png", false)
                            # Case sensitive check for warnings
                            push!(badExt, joinpath(root, file))
                        end

                        rawpath = joinpath(root, file)
                        push!(res, (relpath(rawpath, path), rawpath, atlas))
                    end
                end
            end
        end
    end

    return res
end

function findAllExternalSprites()
    warnOnBadExts = get(config, "prompt_bad_resource_exts", true)
    res = Tuple{String, String, String}[]

    folderFileBadExt = String[]
    zipFileBadExt = String[]

    targetFolders = getCelesteModDirs()
    targetZips = getCelesteModZips()
    atlasesPath = joinpath("Graphics", "Atlases")

    # ZipFile uses unix paths
    atlasesUnixPath = replace(atlasesPath, "\\" => "/")

    for target in targetFolders
        append!(res, findExternalSpritesInFolder(target, atlasesPath, folderFileBadExt))
    end

    for target in targetZips
        append!(res, findExternalSpritesInZip(target, atlasesPath, atlasesUnixPath, zipFileBadExt))
    end

    if warnOnBadExts
        confirmed = warnBadExtensions(folderFileBadExt, zipFileBadExt, ".png")

        if confirmed
            return findAllExternalSprites()
        end
    end

    return res
end

function findExternalModulesInZip(filename::String, pluginPath::String, allowRetry::Bool=true)
    res = Tuple{String, String, String}[]
    fileCache = get(resourceZipCache, filename, nothing)

    if fileCache !== nothing
        for (resource, data) in fileCache
            if startswith(resource, pluginPath) && hasExt(resource, ".jl")
                push!(res, (filename, resource, data))
            end
        end

    else
        if allowRetry
            cacheZipContent(filename)

            return findExternalModulesInZip(filename, pluginPath, false)
        end
    end

    return res
end

function findExternalModulesInFolder(filename::String, pluginPath::String)
    res = Tuple{String, String, String}[]

    pluginDir = joinpath(filename, pluginPath)

    if !isdir(pluginDir)
        return res
    end

    for (root, dirs, files) in walkdir(pluginDir)
        for file in files
            if hasExt(file, ".jl")
                rawpath = joinpath(root, file)

                push!(res, (rawpath, rawpath, open(f -> read(f, String), rawpath)))
            end
        end
    end

    return res
end

function findAllExternalModules(pluginPath::String)
    res = Tuple{String, String, String}[]

    targetFolders = getCelesteModDirs()
    targetZips = getCelesteModZips()

    for target in targetFolders
        append!(res, findExternalModulesInFolder(target, pluginPath))
    end

    for target in targetZips
        append!(res, findExternalModulesInZip(target, pluginPath))
    end

    return res
end

function findChangedExternalSprites()
    warnOnBadExts = get(config, "prompt_bad_resource_exts", true)
    res = Tuple{String, String, String}[]

    changedFilenames = FileWatcher.processWatchEvents(FileWatcher.basePath)

    folderFileBadExt = String[]
    zipFileBadExt = String[]

    # ZipFile uses unix paths
    atlasesPath = joinpath("Graphics", "Atlases")
    atlasesUnixPath = replace(atlasesPath, "\\" => "/")

    for filename in changedFilenames
        if hasExt(filename, ".zip")
            append!(res, findExternalSpritesInZip(filename, atlasesPath, atlasesUnixPath, zipFileBadExt))

        elseif hasExt(filename, ".png")
            hasModRoot, modRoot = getModRoot(filename)

            if hasModRoot
                relative = relpath(filename, modRoot)
                parts = splitpath(relative)

                if parts[1] == "Graphics" && parts[2] == "Atlases"
                    push!(res, (joinpath(parts[4:end]...), filename, parts[3]))
                end
            end
        end
    end

    if warnOnBadExts
        confirmed = warnBadExtensions(folderFileBadExt, zipFileBadExt, ".png")

        if confirmed
            return findChangedExternalSprites()
        end
    end

    return res
end

function findLanguageFileInZip(filename::String, languageFilename::String, allowRetry::Bool=true)
    fileCache = get(resourceZipCache, filename, nothing)

    if fileCache !== nothing
        for (resource, data) in fileCache
            if resource == languageFilename
                return data
            end
        end

    else
        if allowRetry
            cacheZipContent(filename)

            return findLanguageFileInZip(filename, languagePath, false)
        end
    end

    return nothing
end

function loadLangdata()
    targetFolders = joinpath.(getCelesteModDirs(), "Ahorn")
    targetZips = getCelesteModZips()

    lang = get(config, "language", "en_gb")
    internalFilename = abspath("../lang/$lang.lang")

    res = parseLangfile(read(internalFilename, String))

    for folder in targetFolders
        fn = joinpath(folder, "lang", "$lang.lang")
        if ispath(fn)
            content = read(fn, String)
            parseLangfile(content, init=res)
        end
    end

    for fn in targetZips
        content = findLanguageFileInZip(fn, "Ahorn/lang/$(lang).lang")

        if content !== nothing
            parseLangfile(content, init=res)
        end
    end

    return res
end

function findExternalModules(args::String...)
    res = String[]

    targetFolders = vcat(
        joinpath.(getCelesteModDirs(), "Ahorn"),
        getAhornModDirs()
    )

    for folder in targetFolders
        path = joinpath(folder, args...)
        if isdir(path)
            for (root, dirs, files) in walkdir(path)
                for file in files
                    if hasExt(file, ".jl")
                        push!(res, joinpath(root, file))
                    end
                end
            end
        end
    end

    return res
end

function loadExternalZipModules!(loadedModules::Dict{String, Module}, loadedNames::Array{String, 1}, args::String...)
    # ZipFile uses unix paths
    path = joinpath("Ahorn", args...)
    path = replace(path, "\\" => "/")

    targets = vcat(
        getCelesteModZips(),
        getAhornModZips()
    )

    for target in targets
        targetMtime = stat(target).mtime
        moduleData = findExternalModulesInZip(target, path)

        for (rawpath, relative, content) in moduleData
            fakeFn = rawpath * ":" * relative * ".from_zip"

            try
                if targetMtime > get(loadedModulesTimes, fakeFn, 0)
                    loadedModules[fakeFn] = Base.eval(Ahorn, Meta.parse(strip(content)))
                    loadedModulesTimes[fakeFn] = targetMtime
                end

                if !(fakeFn in loadedNames)
                    push!(loadedNames, fakeFn)
                end

            catch e
                println(Base.stderr, "! Failed to load \"$relative\" from \"$target\"")
                println(Base.stderr, e)
            end
        end
    end
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# Use Chemfiles_jll & the artifact system on julia >= 1.3
if VERSION >= v"1.3.0"
    Sys.exit()
end

include("build_Chemfiles.v0.10.2.jl")


text:
import Base.:+
import Base.:*

export Biquaternion
export getrotation
export gettranslation

"""
    Represents a bi-quaternion for describing rotation and translation in a unified way.

fields: real and dual.
"""
struct Biquaternion
    real::Quaternion
    dual::Quaternion
    Biquaternion(qr::Quaternion, qd::Quaternion) = new(qr, qd)
    Biquaternion(q::Quaternion, t::ℝ³) = begin
        qr = normalize(q)
        qd = 0.5 * Quaternion([0; vec(t)]) * qr
        new(qr, qd)
    end
end

"""
    show(b)

Print a string representation of the given Biquaternion `b`.
"""
Base.show(io::IO, q::Biquaternion) = print(io, "$(q.real) + $(q.dual)𝜺")

"""
    Biquaternion(q)

Construct a Biquaternion with the given Biquaternion `q` and also normalize.
"""
Biquaternion(q::Biquaternion) = begin
    q̂ = normalize(q)
    Biquaternion(q̂.real, q̂.dual)
end

"""
    Biquaternion(rotation)

Construct a Biquaternion with the given Quaternion `rotation` and also normalize.
"""
Biquaternion(rotation::Quaternion) = Biquaternion(normalize(rotation), ℝ³(0, 0, 0))

"""
    Biquaternion(translation)

Construct a Biquaternion with the given ℝ³ `translation`.
"""
Biquaternion(translation::ℝ³) = Biquaternion(Quaternion(1, 0, 0, 0), translation)

vec(q::Biquaternion) = [vec(q.real); vec(q.dual)]

+(q₁::Biquaternion, q₂::Biquaternion) = Biquaternion(q₁.real + q₂.real, q₁.dual + q₂.dual)
-(q₁::Biquaternion, q₂::Biquaternion) = Biquaternion(q₁.real - q₂.real, q₁.dual - q₂.dual)
*(q₁::Biquaternion, q₂::Biquaternion) = Biquaternion(q₁.real * q₂.real,
                                                     q₁.real * q₂.dual + q₁.dual * q₂.real)
*(q::Biquaternion, λ::Real) = Biquaternion(λ * q.real, λ * q.dual)
*(λ::Real, q::Biquaternion) = q * λ
Base.conj(q::Biquaternion) = Biquaternion(conj(q.real), conj(q.dual))
norm(q::Biquaternion) = begin
    qrnorm = norm(q.real)
    Q = (conj(q.real) * q.dual + q.real * conj(q.dual)) * (1 / 2qrnorm)
    scalar = vec(Q)[1]
    sqrt(qrnorm^2 + scalar^2)
end
normalize(q::Biquaternion) = begin
    magnitude = norm(q)
    @assert(magnitude > 1e-5, "The magnitude is almost equal to zero, too small.")
    q * (1 / magnitude)
end
getrotation(q::Biquaternion) = q.real
gettranslation(q::Biquaternion) = ℝ³(vec(2q.dual * conj(q.real))[2:4])

Base.isapprox(q1::Biquaternion, q2::Biquaternion) = isapprox(q1.real, q2.real) &&
                                                    isapprox(q1.dual, q2.dual)

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
@safetestset "isean" begin
    using Validator

    @test Validator.isean("9421023610112") === true
    @test Validator.isean("1234567890128") === true
    @test Validator.isean("4012345678901") === true
    @test Validator.isean("9771234567003") === true
    @test Validator.isean("9783161484100") === true
    @test Validator.isean("73513537") === true
    @test Validator.isean("00012345600012") === true
    @test Validator.isean("10012345678902") === true
    @test Validator.isean("20012345678909") === true

    @test Validator.isean("5901234123451") === false
    @test Validator.isean("079777681629") === false
    @test Validator.isean("0705632085948") === false
end


text:
#
# This function fits penalized variance component model with responses y, design matrix X
# and variance components in V by a multiplicative MM algorithm.
#
using Polynomial
function VClasso(y::Vector{Float64}, X::Matrix{Float64},
                 V::Array{Matrix{Float64}, 1}, # Z::Matrix{Float64}, # add Z for Longitudinal design
                 λfactor::Vector{Int64}, # 0 is not penalized, 1 is penalized to 1.0
                       print::Bool=false, λgrid::Vector{Float64}=[0.0])
  n = length(y)   # no. observations
  n2 = n * n
  p = size(X, 2)  # no. predictors
  m = length(V)   # no. variance components
  ngrid = length(λgrid);
  σ2grid = zeros(m, ngrid);
  Itergrid = zeros(ngrid);
  LogLgrid = zeros(ngrid);
  Objgrid = zeros(ngrid);
  # initialize algorithm
  for k = 1:ngrid
   β = X \ y ;
   σ2 = 1e-2 * ones(m);
   σ = sqrt(σ2);
   Ω = zeros(n, n);
   for j = 1:m
     BLAS.axpy!(n2, σ2[j], V[j], 1, Ω, 1)
   end
   Ωchol = cholfact(Ω);
   Ωinv = inv(Ωchol);
   res = y - X * β;
   v = Ωchol \ res ;
   ##
   loglConst = - 0.5 * n * log(2.0 * pi);
   logl = loglConst - 0.5 * logdet(Ωchol) - 0.5 * dot(res, v)
   Iter = 1 :: Int
   logL = 0.0
   Obj = 0.0

   # MM loop
   maxiter = 10000::Int
   tol = 1e-4
   λ = λgrid[k];
   logllasso = -logl + λ * sumabs(λfactor[1:end].*σ[1:end-1]);
   for i = 1:maxiter
    if print
      println(i, " ",  "Obj ", logllasso, " ", "logL", logl)
    end

    # update variance components
    fill!(Ω, 0.0);
    for l = 1: m - 1
      quartic = BLAS.dot(n2, Ωinv, 1, V[l], 1)
      cubic = λ * λfactor[l]
      constant = -σ[l]^4 * ((v'* V[l])*v)[]
      @inbounds σ[l] = min_root(quartic, cubic, constant, l, V, σ, res)
      @inbounds σ2[l] = σ[l] .^2
      @inbounds BLAS.axpy!(n2, σ2[l], V[l], 1, Ω, 1)
    end
    σ2[m] *= BLAS.nrm2(v) / sqrt(trace(Ωinv))
    σ[m] = sqrt(σ2[m])
    BLAS.axpy!(n2, σ2[m], V[m], 1, Ω, 1)
    Ωchol = cholfact(Ω, :L);
    Ωinv = inv(Ωchol);
    Xnew = Ωchol[:L] \ X;
    ynew = Ωchol[:L] \ y;
    # update fixed effects
    β = Xnew \ ynew;
    res = y - X * β;
    v = Ωchol \ res;
    # check convergence
    logllassoold = logllasso
    logl = loglConst - 0.5 * logdet(Ωchol) - 0.5 * dot(res, v)
    logllasso = -logl + λ * sumabs(λfactor[1:end].*σ[1:end-1]);
    if abs(logllasso - logllassoold) < tol * (abs(logllassoold) + 1.0)
      #println(i," ", logl)
      Obj = logllasso
      logL = logl;
      Iter = i
      break
    end
   end # end of loop over maxiter
   σ2grid[:, k] = σ2;
   Itergrid[k] = Iter;
   LogLgrid[k] = logL;
   Objgrid[k] = Obj;
  end # end of loop over λgrid
    # output
  return σ2grid, Itergrid, LogLgrid, Objgrid
end

# this function returns the \sigma that minimize the objective function: negative loglikelihood plus penalty term
function min_root(quartic, cubic, constant, index, V, σ, res)
  # find the real positive root
  complex_root = roots(Poly([quartic, cubic, 0, 0, constant]))
  # keep real roots
  real_root = real(complex_root[imag(complex_root) .== 0])
  real_pos_root = real_root[real_root .> 0]
  # find the real negative root
  complex_root = roots(Poly([quartic, -cubic, 0, 0, constant]))
  real_root = real(complex_root[imag(complex_root) .== 0])
  real_neg_root = real_root[real_root .< 0]
  root = [real_pos_root, 0.0, real_neg_root]
  obj = copy(root);
  m = length(σ)
  n = size(V[1])[1]
  n2 = n^2
  for i = 1:length(root)
    σ[index] = root[i]
    σ2 = σ .^2
    Ω = zeros(n, n);
    for j = 1:m
      BLAS.axpy!(n2, σ2[j], V[j], 1, Ω, 1)
    end
    Ωchol = cholfact(Ω);
    Ωinv = inv(Ωchol);
    v = Ωchol \ res ;
    ##
    obj[i] =  0.5 * logdet(Ωchol) + 0.5 * dot(res, v) + cubic * sumabs(σ[1:end-1])
  end
  return root[indmin(obj)]

end


text:
function get_comp(c::Char)
    if c=='A' return 'T'
    elseif c=='T' return 'A'
    elseif c=='C' return 'G'
    elseif c=='G' return 'C'
    else
        println("$c is not A, T, C, or G.")
        return 'N'
    end
end

get_revcomp(s::String) = reverse(join([get_comp(c) for c in s]))


text:
@testset "Constructors" begin
    h1 = SwissDict()
    @test length(h1) == 0
    @test isempty(h1) == true
    @test h1.idxfloor == 1
    @test length(h1.keys) == 16
    @test length(h1.vals) == 16
    @test length(h1.slots) == 1
    @test eltype(h1) == Pair{Any, Any}
    @test keytype(h1) == Any
    @test valtype(h1) == Any
end

@testset "SwissDict" begin
    h = SwissDict()
    for i=1:10000
        h[i] = i+1
    end
    for i=1:10000
        @test (h[i] == i+1)
    end
    for i=1:2:10000
        delete!(h, i)
    end
    for i=1:2:10000
        h[i] = i+1
    end
    for i=1:10000
        @test (h[i] == i+1)
    end
    for i=1:10000
        delete!(h, i)
    end
    @test isempty(h)
    h[77] = 100
    @test h[77] == 100
    for i=1:10000
        h[i] = i+1
    end
    for i=1:2:10000
        delete!(h, i)
    end
    for i=10001:20000
        h[i] = i+1
    end
    for i=2:2:10000
        @test h[i] == i+1
    end
    for i=10000:20000
        @test h[i] == i+1
    end
    h = SwissDict{Any,Any}("a" => 3)
    @test h["a"] == 3
    h["a","b"] = 4
    @test h["a","b"] == h[("a","b")] == 4
    h["a","b","c"] = 4
    @test h["a","b","c"] == h[("a","b","c")] == 4

    @testset "eltype, keytype and valtype" begin
        @test eltype(h) == Pair{Any,Any}
        @test keytype(h) == Any
        @test valtype(h) == Any

        td = SwissDict{AbstractString,Float64}()
        @test eltype(td) == Pair{AbstractString,Float64}
        @test keytype(td) == AbstractString
        @test valtype(td) == Float64
        @test keytype(Dict{AbstractString,Float64}) === AbstractString
        @test valtype(Dict{AbstractString,Float64}) === Float64
    end
    # test rethrow of error in ctor
    @test_throws DomainError SwissDict((sqrt(p[1]), sqrt(p[2])) for p in zip(-1:2, -1:2))
end

@testset "SwissDict on pairs" begin
	let x = SwissDict(3=>3, 5=>5, 8=>8, 6=>6)
	    pop!(x, 5)
	    for k in keys(x)
	        SwissDict{Int,Int}(x)
	        @test k in [3, 8, 6]
	    end
	end

    let y = SwissDict{Any, Int}(3=>3, 5=>5, "8"=>8, 6=>6)
        pop!(y, "8")
        for k in keys(y)
            SwissDict{Int,Int}(y)
            @test k in [3, 5, 6]
        end
    end

    d = @inferred SwissDict(Pair(1,1), Pair(2,2), Pair(3,3))
    @test isa(d, SwissDict)
    @test d == SwissDict(1=>1, 2=>2, 3=>3)
    @test eltype(d) == Pair{Int,Int}
end

@testset "KeyError" begin
    let z = SwissDict()
        get_KeyError = false
        try
            z["a"]
        catch _e123_
            get_KeyError = isa(_e123_,KeyError)
        end
        @test get_KeyError
    end
end

@testset "Filter function" begin
    _d = SwissDict("a"=>0)
    @test isa([k for k in filter(x->length(x)==1, collect(keys(_d)))], Vector{String})

    h = SwissDict{Int, Int}()
    for i in 1:100
        h[i] = i
    end
    filter!(x->isodd(x.first), h)
    for i in 1:2:100
        @test h[i] == i
    end
end

@testset "typeof" begin
    d = SwissDict(((1, 2), (3, 4)))
    @test d[1] === 2
    @test d[3] === 4
    d2 = SwissDict(1 => 2, 3 => 4)
    d3 = SwissDict((1 => 2, 3 => 4))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == SwissDict{Int,Int}

    d = SwissDict(((1, 2), (3, "b")))
    @test d[1] === 2
    @test d[3] == "b"
    d2 = SwissDict(1 => 2, 3 => "b")
    d3 = SwissDict((1 => 2, 3 => "b"))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == SwissDict{Int,Any}

    d = SwissDict(((1, 2), ("a", 4)))
    @test d[1] === 2
    @test d["a"] === 4
    d2 = SwissDict(1 => 2, "a" => 4)
    d3 = SwissDict((1 => 2, "a" => 4))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == SwissDict{Any,Int}

    d = SwissDict(((1, 2), ("a", "b")))
    @test d[1] === 2
    @test d["a"] == "b"
    d2 = SwissDict(1 => 2, "a" => "b")
    d3 = SwissDict((1 => 2, "a" => "b"))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == SwissDict{Any,Any}
end

@testset "type of SwissDict constructed from varargs of Pairs" begin
    @test SwissDict(1=>1, 2=>2.0) isa SwissDict{Int,Real}
    @test SwissDict(1=>1, 2.0=>2) isa SwissDict{Real,Int}
    @test SwissDict(1=>1.0, 2.0=>2) isa SwissDict{Real,Real}

    for T in (Nothing, Missing)
        @test SwissDict(1=>1, 2=>T()) isa SwissDict{Int,Union{Int,T}}
        @test SwissDict(1=>T(), 2=>2) isa SwissDict{Int,Union{Int,T}}
        @test SwissDict(1=>1, T()=>2) isa SwissDict{Union{Int,T},Int}
        @test SwissDict(T()=>1, 2=>2) isa SwissDict{Union{Int,T},Int}
    end
end

@testset "equality" for eq in (isequal, ==)
    @test  eq(SwissDict(), SwissDict())
    @test  eq(SwissDict(1 => 1), SwissDict(1 => 1))
    @test !eq(SwissDict(1 => 1), SwissDict())
    @test !eq(SwissDict(1 => 1), SwissDict(1 => 2))
    @test !eq(SwissDict(1 => 1), SwissDict(2 => 1))

    # Generate some data to populate dicts to be compared
    data_in = [ (rand(1:1000), randstring(2)) for _ in 1:1001 ]

    # Populate the first dict
    d1 = SwissDict{Int, AbstractString}()
    for (k, v) in data_in
        d1[k] = v
    end
    data_in = collect(d1)
    # shuffle the data
    for i in 1:length(data_in)
        j = rand(1:length(data_in))
        data_in[i], data_in[j] = data_in[j], data_in[i]
    end
    # Inserting data in different (shuffled) order should result in
    # equivalent dict.
    d2 = SwissDict{Int, AbstractString}()
    for (k, v) in data_in
        d2[k] = v
    end

    @test eq(d1, d2)
    d3 = copy(d2)
    d4 = copy(d2)
    # Removing an item gives different dict
    delete!(d1, data_in[rand(1:length(data_in))][1])
    @test !eq(d1, d2)
    # Changing a value gives different dict
    d3[data_in[rand(1:length(data_in))][1]] = randstring(3)
    !eq(d1, d3)
    # Adding a pair gives different dict
    d4[1001] = randstring(3)
    @test !eq(d1, d4)

    @test eq(SwissDict(), sizehint!(SwissDict(),96))

    # Dictionaries of different types
    @test !eq(SwissDict(1 => 2), SwissDict("dog" => "bone"))
    @test eq(SwissDict{Int,Int}(), SwissDict{AbstractString,AbstractString}())
end

@testset "equality special cases" begin
    @test SwissDict(1=>0.0) == SwissDict(1=>-0.0)
    @test !isequal(SwissDict(1=>0.0), SwissDict(1=>-0.0))

    @test SwissDict(0.0=>1) != SwissDict(-0.0=>1)
    @test !isequal(SwissDict(0.0=>1), SwissDict(-0.0=>1))

    @test SwissDict(1=>NaN) != SwissDict(1=>NaN)
    @test isequal(SwissDict(1=>NaN), SwissDict(1=>NaN))

    @test SwissDict(NaN=>1) == SwissDict(NaN=>1)
    @test isequal(SwissDict(NaN=>1), SwissDict(NaN=>1))

    @test ismissing(SwissDict(1=>missing) == SwissDict(1=>missing))
    @test isequal(SwissDict(1=>missing), SwissDict(1=>missing))

    @test SwissDict(missing=>1) == SwissDict(missing=>1)
    @test isequal(SwissDict(missing=>1), SwissDict(missing=>1))
end

@testset "get!" begin
    f(x) = x^2
    d = SwissDict(8=>19)
    @test get!(d, 8, 5) == 19
    @test get!(d, 19, 2) == 2

    @test get!(d, 42) do  # d is updated with f(2)
        f(2)
    end == 4

    @test get!(d, 42) do  # d is not updated
        f(200)
    end == 4

    @test get(d, 13) do   # d is not updated
        f(4)
    end == 16

    @test d == SwissDict(8=>19, 19=>2, 42=>4)

    d1 = SwissDict{Int8, Char}()
    d1[1%Int8] = 'a'
    @test get!(d1, 1, 'b') == 'a'
    @test get!(d1, 2, 'c') == 'c'
end

@testset "push!" begin
    d = SwissDict()
    @test push!(d, 'a' => 1) === d
    @test d['a'] == 1
    @test push!(d, 'b' => 2, 'c' => 3) === d
    @test d['b'] == 2
    @test d['c'] == 3
    @test push!(d, 'd' => 4, 'e' => 5, 'f' => 6) === d
    @test d['d'] == 4
    @test d['e'] == 5
    @test d['f'] == 6
    @test length(d) == 6
end

@testset "pop!" begin
    d = SwissDict(1=>2, 3=>4)
    @test pop!(d, 1) == 2
    @test_throws KeyError pop!(d, 1)
    @test pop!(d, 1, 0) == 0
    @test pop!(d) == (3=>4)
    @test_throws ArgumentError pop!(d)
end

@testset "keys as a set" begin
    d = SwissDict(1=>2, 3=>4)
    @test keys(d) isa AbstractSet
    @test empty(keys(d)) isa AbstractSet
    let i = keys(d) ∩ Set([1,2])
        @test i isa AbstractSet
        @test i == Set([1])
    end
    @test Set(string(k) for k in keys(d)) == Set(["1","3"])
end

@testset "find" begin
    @test findall(isequal(1), SwissDict(:a=>1, :b=>2)) == [:a]
    @test sort(findall(isequal(1), SwissDict(:a=>1, :b=>1))) == [:a, :b]
    @test isempty(findall(isequal(1), SwissDict()))
    @test isempty(findall(isequal(1), SwissDict(:a=>2, :b=>3)))

    @test findfirst(isequal(1), SwissDict(:a=>1, :b=>2)) == :a
    @test findfirst(isequal(1), SwissDict(:a=>1, :b=>1, :c=>3)) in (:a, :b)
    @test findfirst(isequal(1), SwissDict()) === nothing
    @test findfirst(isequal(1), SwissDict(:a=>2, :b=>3)) === nothing
end

@testset "haskey" begin
    h = SwissDict(1=>2, 2=>3)
    @test haskey(h, 1) == true
    @test haskey(h, 2) == true
    @test haskey(h, 3) == false
    @test !haskey(h, "1")
end

@testset "getkey" begin
    h = SwissDict(1=>2, 3 => 6, 5=>10)
    @test getkey(h, 1, 7) == 1
    @test getkey(h, 4, 6) == 6
    @test getkey(h, "1", 8) == 8
end

@testset "ArgumentError" begin
    @test_throws ArgumentError SwissDict(0)
    @test_throws ArgumentError SwissDict([1])
    @test_throws ArgumentError SwissDict([(1,2),0])
end

@testset "empty tuple" begin
    h = SwissDict(())
    @test length(h) == 0
end

@testset "empty" begin
    h = SwissDict()
    for i=1:10000
        h[i] = i+1
    end
    prev_sz = length(h.keys)
    @test length(h) != 0
    empty!(h)
    @test length(h) == 0
    @test length(h.keys) == length(h.vals) == prev_sz
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Start by checking the byte-order.
ordered = 0x0123456789ABCDEF
@test 3512700564088504e-318 == reinterpret(Float64,ordered)

min_double64 = 0x0000000000000001
@test 5e-324 == reinterpret(Float64,min_double64)

max_double64 = 0x7fefffffffffffff
@test 1.7976931348623157e308 == reinterpret(Float64,max_double64)

# Start by checking the byte-order.
ordered = 0x01234567
@test float32(2.9988165487136453e-38) == reinterpret(Float32,ordered)

min_float32 = 0x00000001
@test float32(1.4e-45) == reinterpret(Float32,min_float32)

max_float32 = 0x7f7fffff
@test float32(3.4028234e38) == reinterpret(Float32,max_float32)

ordered = 0x0123456789ABCDEF
diy_fp  = Grisu.Float(reinterpret(Float64,ordered))
@test 0x12 - 0x3FF - 52 == uint64(diy_fp.e)
# The 52 mantissa bits, plus the implicit 1 in bit 52 as a UINT64.
@test 0x0013456789ABCDEF== diy_fp.s

min_double64 = 0x0000000000000001
diy_fp  = Grisu.Float(reinterpret(Float64,min_double64))
@test -0x3FF - 52 + 1 == uint64(diy_fp.e)
# This is a denormal so no hidden bit.
@test 1 == diy_fp.s

max_double64 = 0x7fefffffffffffff
diy_fp  = Grisu.Float(reinterpret(Float64,max_double64))
@test 0x7FE - 0x3FF - 52 == uint64(diy_fp.e)
@test 0x001fffffffffffff== diy_fp.s

ordered = 0x01234567
diy_fp  = Grisu.Float(reinterpret(Float32,ordered))
@test 0x2 - 0x7F - 23 == uint64(diy_fp.e)
# The 23 mantissa bits, plus the implicit 1 in bit 24 as a uint32_t.
@test 0xA34567 == uint64(diy_fp.s)

min_float32 = 0x00000001
diy_fp  = Grisu.Float(reinterpret(Float32,min_float32))
@test -0x7F - 23 + 1 == uint64(diy_fp.e)
# This is a denormal so no hidden bit.
@test 1 == uint64(diy_fp.s)

max_float32 = 0x7f7fffff
diy_fp  = Grisu.Float(reinterpret(Float32,max_float32))
@test 0xFE - 0x7F - 23 == uint64(diy_fp.e)
@test 0x00ffffff == uint64(diy_fp.s)

ordered = 0x0123456789ABCDEF
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float64,ordered)))
@test 0x12 - 0x3FF - 52 - 11 == uint64(diy_fp.e)
@test 0x0013456789ABCDEF<< 11 == diy_fp.s

min_double64 = 0x0000000000000001
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float64,min_double64)))
@test -0x3FF - 52 + 1 - 63 == uint64(diy_fp.e)
# This is a denormal so no hidden bit.
@test 0x8000000000000000== diy_fp.s

max_double64 = 0x7fefffffffffffff
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float64,max_double64)))
@test 0x7FE - 0x3FF - 52 - 11 == uint64(diy_fp.e)
@test (0x001fffffffffffff<< 11) == diy_fp.s

min_double64 = 0x0000000000000001
@test Grisu.isdenormal(reinterpret(Float64,min_double64))
bits = 0x000FFFFFFFFFFFFF
@test Grisu.isdenormal(reinterpret(Float64,bits))
bits = 0x0010000000000000
@test !Grisu.isdenormal(reinterpret(Float64,bits))

min_float32 = 0x00000001
@test Grisu.isdenormal(reinterpret(Float32,min_float32))
bits = 0x007FFFFF
@test Grisu.isdenormal(reinterpret(Float32,bits))
bits = 0x00800000
@test !Grisu.isdenormal(reinterpret(Float32,bits))

diy_fp = Grisu.normalize(Grisu.Float(1.5))
boundary_minus, boundary_plus = Grisu.normalizedbound(1.5)
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# 1.5 does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (1 << 10) == diy_fp.s - boundary_minus.s

diy_fp = Grisu.normalize(Grisu.Float(1.0))
boundary_minus, boundary_plus = Grisu.normalizedbound(1.0)
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# 1.0 does have a significand of the form 2^p (for some p).
# Therefore its lower boundary is twice as close as the upper boundary.
@test boundary_plus.s - diy_fp.s > diy_fp.s - boundary_minus.s
@test (1 << 9) == diy_fp.s - boundary_minus.s
@test (1 << 10) == boundary_plus.s - diy_fp.s

min_double64 = 0x0000000000000001
diy_fp = Grisu.normalize(Grisu.Float(reinterpret(Float64,min_double64)))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float64,min_double64))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# min-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
# Denormals have their boundaries much closer.
@test (uint64(1) << 62) == diy_fp.s - boundary_minus.s

smallest_normal64 = 0x0010000000000000
diy_fp = Grisu.normalize(reinterpret(Float64,smallest_normal64))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float64,smallest_normal64))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# Even though the significand is of the form 2^p (for some p), its boundaries
# are at the same distance. (This is the only exception).
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (1 << 10) == diy_fp.s - boundary_minus.s

largest_denormal64 = 0x000FFFFFFFFFFFFF
diy_fp = Grisu.normalize(reinterpret(Float64,largest_denormal64))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float64,largest_denormal64))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (1 << 11) == diy_fp.s - boundary_minus.s

max_double64 = 0x7fefffffffffffff
diy_fp = Grisu.normalize(reinterpret(Float64,max_double64))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float64,max_double64))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# max-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (1 << 10) == diy_fp.s - boundary_minus.s

kOne64 = uint64(1)
diy_fp  = Grisu.normalize(Grisu.Float(float32(1.5)))
boundary_minus, boundary_plus = Grisu.normalizedbound(float32(1.5))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# 1.5 does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
# Normalization shifts the significand by 8 bits. Add 32 bits for the bigger
# data-type, and remove 1 because boundaries are at half a ULP.
@test (kOne64 << 39) == diy_fp.s - boundary_minus.s

diy_fp  = Grisu.normalize(Grisu.Float(float32(1.0)))
boundary_minus, boundary_plus = Grisu.normalizedbound(float32(1.0))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# 1.0 does have a significand of the form 2^p (for some p).
# Therefore its lower boundary is twice as close as the upper boundary.
@test boundary_plus.s - diy_fp.s > diy_fp.s - boundary_minus.s
@test (kOne64 << 38) == diy_fp.s - boundary_minus.s
@test (kOne64 << 39) == boundary_plus.s - diy_fp.s

min_float32 = 0x00000001
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float32,min_float32)))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float32,min_float32))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# min-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
# Denormals have their boundaries much closer.
@test (kOne64 << 62) == diy_fp.s - boundary_minus.s

smallest_normal32 = 0x00800000
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float32,smallest_normal32)))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float32,smallest_normal32))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# Even though the significand is of the form 2^p (for some p), its boundaries
# are at the same distance. (This is the only exception).
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (kOne64 << 39) == diy_fp.s - boundary_minus.s

largest_denormal32 = 0x007FFFFF
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float32,largest_denormal32)))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float32,largest_denormal32))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (kOne64 << 40) == diy_fp.s - boundary_minus.s

max_float32 = 0x7f7fffff
diy_fp  = Grisu.normalize(Grisu.Float(reinterpret(Float32,max_float32)))
boundary_minus, boundary_plus = Grisu.normalizedbound(reinterpret(Float32,max_float32))
@test diy_fp.e == boundary_minus.e
@test diy_fp.e == boundary_plus.e
# max-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
@test diy_fp.s - boundary_minus.s == boundary_plus.s - diy_fp.s
@test (kOne64 << 39) == diy_fp.s - boundary_minus.s


text:
include("graphics.jl")

function p1() #Player 1 gets to enter the word they want
    word = string(input("Player 1: Please enter the word to be guessed by Player 2\n"))
    for char in word
        if isnumeric(char)
            println("Please enter a word with no numbers!")
            p1()
        end
    end
    return lowercase(word)
end

function p2guess(already_guessed) #Function where P2 guesses a letter
    guess = string(input("Player 2: Please guess a letter\n"))
    if length(guess) != 1 #check whether its length is more or less than 1
        println("Please enter a letter, not word!")
        p2guess(already_guessed)
    end
    for char in guess #check whether it has numbers
        if isnumeric(char)
            println("Please enter a letter, not number!")
            p2guess(already_guessed)
        end
    end
    if guess in already_guessed #check whether the letter is already guessed
        println("Please enter a letter that's not guessed yet!")
        p2guess(already_guessed)
    end
    return lowercase(guess)
end

function display_hangman(lives, correctletters, incorrectletters)
    println(hangman[length(incorrectletters)+1]) #prints the correct hangman picture from graphics.jl
    println("\nLives left: $(lives)") #prints the lives left
    blanks = ["_" for i in 1:length(word)]
    for num in 1:length(word)
        letter = string(word[num])
        for i in correctletters
            if i == letter
                tempblanks = blanks[(num+1):end]
                blanks = append!(blanks[1:num-1], [letter])
                blanks = append!(blanks, tempblanks) #creating the new blanks array with correctletters filled in
            end
        end

    end
    print("Word: ")
    for j in blanks
        print(j, " ")
    end
    println()
    print("Missed Letters: ")
    for h in incorrectletters
        print(h, ", ")
    end
    println()
end

function playagain() #play again functionality
    playagain = lowercase(string(input("Do you want to play again? (Y/N)\n")))
    if playagain == "y"
        return true
    else
        return false
    end
end

function main()
    global word = p1()
    correctletters = []
    incorrectletters = []
    lettersinword = unique(split(word, ""))
    lives = 6
    gamedone = false
    run(`clear`)
    while true
        display_hangman(lives, correctletters, incorrectletters)
        guess = p2guess(append!(correctletters, incorrectletters))
        if (guess in lettersinword)
            push!(correctletters, guess)
            #checks whether all the letters have been found or not
            allfound = true #assuming that it is true until an exception is found (ie. whether player2 has won or not)
            for i in 1:length(word)
                if (string(word[i]) in correctletters) == false
                    allfound = false #exception is found here
                    break
                end
            end
            if allfound == true
                println("All letters have been found. Well done!")
                println("Player 2 has won! The word was $(word)!")
                gamedone = true
            end
        else
            push!(incorrectletters, guess)
            lives = lives - 1
            #checks whether they've ran out of lives
            if lives == 0
                display_hangman(lives, correctletters, incorrectletters)
                println("The word is $(word)")
                println("Player 1 has won! Player 2 didn't manage to guess the word!")
                gamedone = true
            end
        end
        #check whether they want to play again
        if gamedone == true
            if playagain()
                global word = p1()
                correctletters = []
                incorrectletters = []
                lettersinword = unique(split(word, ""))
                lives = 6
                gamedone = false
                run(`clear`)
            else
                break
            end
        end
    end
end

main()

max_stars_repo_name: fungktt/JuliaGames

max_issues_repo_name: fungktt/JuliaGames

max_forks_repo_name: fungktt/JuliaGames


text:
page_years = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
page_year_0 = 2015

function getpageindexfromyear(year)
    i = findfirst(isequal(year), page_years)
    if i == 0
        error("Invalid PAGE year: $year.")
    end
    return i
end

function getperiodlength(year)      # same calculations made for yagg_periodspan in the model
    i = getpageindexfromyear(year)

    if year == page_years[1]
        start_year = page_year_0
    else
        start_year = page_years[i - 1]
    end

    if year == page_years[end]
        last_year = page_years[end]
    else
        last_year = page_years[i + 1]
    end

    return (last_year - start_year) / 2
end

"""
Applies undiscounting factor to get the SCC, discounted to the emissions year instead of the base year.
"""
function undiscount_scc(m::Model, year::Int)
    df = m[:EquityWeighting, :df_utilitydiscountfactor]
    consfocus0 = m[:GDP, :cons_percap_consumption_0][1]
    consfocus = m[:GDP, :cons_percap_consumption][:, 1]
    emuc = m[:EquityWeighting, :emuc_utilityconvexity]
    sccii = getpageindexfromyear(year)

    return df[sccii] * ((consfocus[sccii] / consfocus0)^-emuc)
end

@defcomp ExtraEmissions begin
    e_globalCO2emissions = Parameter(index=[time], unit="Mtonne/year")
    pulse_size = Parameter(unit="Mtonne CO2")
    pulse_year = Parameter()
    e_globalCO2emissions_adjusted = Variable(index=[time], unit="Mtonne/year")

    function run_timestep(p, v, d, t)
        if gettime(t) == p.pulse_year
            # pulse is applied to the years around this year, as a triangular distribution
            v.e_globalCO2emissions_adjusted[t] = p.e_globalCO2emissions[t] + p.pulse_size / getperiodlength(p.pulse_year)
        else
            v.e_globalCO2emissions_adjusted[t] = p.e_globalCO2emissions[t]
        end
    end
end

"""
    compute_scc(
        m::Model = get_model();
        year::Union{Int, Nothing} = nothing,
        eta::Union{Float64, Nothing} = nothing,
        prtp::Union{Float64, Nothing} = nothing,
        equity_weighting::Bool = true,
        pulse_size = 100_000.,
        n::Union{Int,Nothing}=nothing,
        trials_output_filename::Union{String, Nothing} = nothing,
        seed::Union{Int, Nothing} = nothing)

Computes the social cost of CO2 for an emissions pulse in `year` for the provided the model `m`.
If no model is provided, the default model from get_model() is used.
Units of the returned value are \$ per metric tonne of CO2.

The discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity
and ptp_timepreference in the model. If no values are provided, the discount factors will be computed using the default
PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.

The size of the marginal emission defaults to 100_000 metric megatonnes of CO2 (Mtonne CO2), and this
pulse can be modified with the `pulse_size` keyword argument, in metric megatonnes of CO2 (Mtonne CO2)
(this does not change the units of the returned value, which is always normalized by the
`pulse_size` used). A pulse in `year` is actually a gradual increase throughout the timestep preceeding `year`,
followed by a gradual decrease in emissions in the timestep superseeding `year`. Emissions are linearly interpolated
between the points given by the years.

By default, `n = nothing`, and a single value for the "best guess" social cost of CO2 is returned. If a positive
value for keyword `n` is specified, then a Monte Carlo simulation with sample size `n` will run, sampling from
all of PAGE's random variables, and a vector of `n` social cost values will be returned.
Optionally providing a CSV file path to `trials_output_filename` will save all of the sampled trial data as a CSV file.
Optionally providing a `seed` value will set the random seed before running the simulation, allowing the
results to be replicated.
"""
function compute_scc(
        m::Model=get_model();
        year::Union{Int,Nothing}=nothing,
        eta::Union{Float64,Nothing}=nothing,
        prtp::Union{Float64,Nothing}=nothing,
        equity_weighting::Bool=true,
        pulse_size=75000.,
        n::Union{Int,Nothing}=nothing,
        trials_output_filename::Union{String,Nothing}=nothing,
        seed::Union{Int,Nothing}=nothing
        )

    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    if eta !== nothing
        try
            set_param!(m, :emuc_utilityconvexity, eta)      # since eta is a default parameter in PAGE, we need to use `set_param!` if it hasn't been set yet
        catch e
            update_param!(m, :emuc_utilityconvexity, eta)   # or update_param! if it has been set
        end
    end

    if prtp !== nothing
        try
            set_param!(m, :ptp_timepreference, prtp * 100)      # since prtp is a default parameter in PAGE, we need to use `set_param!` if it hasn't been set yet
        catch e
            update_param!(m, :ptp_timepreference, prtp * 100)   # or update_param! if it has been set
        end
    end

    if !equity_weighting
        try
            set_param!(m, :equity_proportion, 0)
        catch e
            update_param!(m, :equity_proportion, 0)
        end
    end

    # note here that we use `pulse_size` as the `delta` keyword argument for
    # the marginal model so we can normalize to $ per ton
    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run

    if n === nothing
        # Run the "best guess" social cost calculation
        run(mm)
        scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
    elseif n < 1
        error("Invalid `n` value, only values >=1 allowed.")
    else
        # Run a Monte Carlo simulation

        simdef = getsim()   # get the default simulation, need to remove :emuc_utilityconvexity and :ptp_timepreference RVs if user specified values for these
        eta !== nothing ? _remove_RV!(simdef, "EquityWeighting.emuc_utilityconvexity") : nothing
        prtp !== nothing ? _remove_RV!(simdef, "EquityWeighting.ptp_timepreference") : nothing

        seed !== nothing ? Random.seed!(seed) : nothing
        Mimi.set_payload!(simdef, (Vector{Float64}(undef, n), year))  # pass the year and a vector for storing SCC values to the `run` function
        si = run(simdef, mm, n, trials_output_filename=trials_output_filename, post_trial_func=_scc_post_trial_func)
        scc = Mimi.payload(si)[1]   # collect the values computed during the post-trial function
    end

    return scc
end

# Post trial function for calculating the scc in Monte Carlo mode
function _scc_post_trial_func(si::SimulationInstance, trial::Int, ntimesteps::Int, tup::Nothing)
    (scc_array, year) = Mimi.payload(si)
    mm = si.models[1]
    scc_array[trial] = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
end

# Helper function for removing a random variable by its parameter name in the model.
#   Their random variable names have been appended with a unique number by Mimi,
#   so need to look them up this way to find the RV name in the simulation definition
#   in order to use the Mimi.delete_RV! function.
#   If it is an array variable, it will be represented by multiple RVs; this deletes them all.
function _remove_RV!(simdef, _name)
    all_rv_names = collect(keys(simdef.rvdict))
    rv_names = all_rv_names[findall(startswith(String(_name)), map(String, all_rv_names))]
    [Mimi.delete_RV!(simdef, rv_name) for rv_name in rv_names]
end

"""
    compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing, pulse_size = 100000.)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided Mimi-PAGE model.

This pulse defaults to 100_000 metric megatonnes of CO2 (Mtonne CO2), and is spread over all years within the
period following `year`. A pulse in `year` is actually a gradual increase throughout the timestep preceeding `year`,
followed by a gradual decrease in emissions in the timestep superseeding `year`. Emissions are linearly interpolated
between the points given by the years.
The SCC will be returned in dollars per ton regardless of pulse size, because it is normalized over the pulse size.
If no model is provided, the default model from MimiPAGE2020.get_model() is used. Discounting scheme can be
specified by the `eta` and `prtp` parameters, which will  update the values of emuc_utilitiyconvexity and
ptp_timepreference in the model.  If no values are provided, the discount factors will be computed using the
default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
"""
function compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, eta::Union{Float64,Nothing}=nothing, prtp::Union{Float64,Nothing}=nothing, pulse_size=75000.)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    if eta != nothing
        try
            set_param!(m, :emuc_utilityconvexity, eta)      # since eta is a default parameter in PAGE, we need to use `set_param!` if it hasn't been set yet
        catch e
            update_param!(m, :emuc_utilityconvexity, eta)   # or update_param! if it has been set
        end
    end

    if prtp != nothing
        try
            set_param!(m, :ptp_timepreference, prtp * 100)      # since prtp is a default parameter in PAGE, we need to use `set_param!` if it hasn't been set yet
        catch e
            update_param!(m, :ptp_timepreference, prtp * 100)   # or update_param! if it has been set
        end
    end

    # note here that we use `pulse_size` as the `delta` keyword argument for
    # the marginal model so we can normalize to $ per ton
    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
    scc_disaggregated = mm[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated] / undiscount_scc(mm.base, year)

    return (scc = scc, scc_disaggregated = scc_disaggregated, mm = mm)
end

"""
    get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing, pulse_size = 100000.)

Returns a Mimi MarginalModel where the provided m is the base model, and the
marginal model has additional emissions of CO2 in year `year`. If no Model m is
provided, the default model from MimiPAGE2020.get_model() is used as the base model.
Note that the returned MarginalModel has already been run. The `pulse_size` defaults
to 100_000 metric megatonnes of CO2 (Mtonne CO2), and is spread over all years within the
period following `year`.
"""
function get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing, pulse_size=75000.)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $page_years.") : nothing

    # note here that we use `pulse_size` as the `delta` keyword argument for
    # the marginal model so we can normalize to $ per ton
    mm = create_marginal_model(m, pulse_size)

    add_comp!(mm.modified, ExtraEmissions, :extra_emissions; after=:co2emissions)
    connect_param!(mm.modified, :extra_emissions => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
    set_param!(mm.modified, :extra_emissions, :pulse_size, pulse_size)
    set_param!(mm.modified, :extra_emissions, :pulse_year, year)

    connect_param!(mm.modified, :CO2Cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted)

    run(mm)
    return mm
end

function compute_scc_mcs(m::Model, samplesize::Int; year::Union{Int,Nothing}=nothing, eta::Union{Float64,Nothing}=nothing, prtp::Union{Float64,Nothing}=nothing, pulse_size=100000.)
    # Setup of location of final results
    scc_results = zeros(samplesize)

    function mc_scc_calculation(sim_inst::SimulationInstance, trialnum::Int, ntimesteps::Int, ignore::Nothing)
        marginal = sim_inst.models[1]
        marg_damages = marginal[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
        scc_results[trialnum] = marg_damages
    end

    # get simulation
    mcs = getsim()

    # Setup models
    eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run

    # Run
    res = run(mcs, mm, samplesize; post_trial_func=mc_scc_calculation)

    scc_results
end