text:
"Run given cells and all the cells that depend on them."
function run_reactive!(notebook::Notebook, cells::Array{Cell, 1})
	# make sure that we're the only run_reactive! being executed - like a semaphor
	token = take!(notebook.executetoken)

	# save the old topology - we'll delete variables assigned from it and re-evalutate its cells
	old_topology = dependent_cells(notebook, cells)

	old_runnable = old_topology.runnable
	to_delete_vars = union(Set{Symbol}(), (runnable.symstate.assignments for runnable in old_runnable)...)
	to_delete_funcs = union(Set{Symbol}(), (Set(keys(runnable.symstate.funcdefs)) for runnable in old_runnable)...)

	# update the cache using the new code and compute the new topology
	for cell in cells
		start_cache!(notebook, cell)
	end
	update_funcdefs!(notebook)
	for cell in cells
		finish_cache!(notebook, cell)
	end

	new_topology = dependent_cells(notebook, union(cells, keys(old_topology.errable)))
	to_run = setdiff(union(new_topology.runnable, old_topology.runnable), keys(new_topology.errable)) # TODO: think if old error cell order matters

	# change the bar on the sides of cells to "running"
	for cell in to_run
		putnotebookupdates!(notebook, clientupdate_cell_running(notebook, cell))
	end
	for (cell, error) in new_topology.errable
		relay_reactivity_error!(cell, error)
		putnotebookupdates!(notebook, clientupdate_cell_output(notebook, cell))
	end

	# delete new variables in case a cell errors (then the later cells show an UndefVarError)
	new_runnable = new_topology.runnable
    to_delete_vars = union(to_delete_vars, (runnable.symstate.assignments for runnable in new_runnable)...)
	to_delete_funcs = union(to_delete_funcs, (Set(keys(runnable.symstate.funcdefs)) for runnable in new_runnable)...)

	new_errable = keys(new_topology.errable)
	to_delete_vars = union(to_delete_vars, (errable.symstate.assignments for errable in new_errable)...)
	to_delete_funcs = union(to_delete_funcs, (Set(keys(errable.symstate.funcdefs)) for errable in new_errable)...)

	workspace = WorkspaceManager.get_workspace(notebook)
	WorkspaceManager.delete_vars(workspace, to_delete_vars)

	local any_interrupted = false
	for cell in to_run
		if any_interrupted
			relay_reactivity_error!(cell, InterruptException())
		else
			any_interrupted |= run_single!(notebook, cell)
		end
		putnotebookupdates!(notebook, clientupdate_cell_output(notebook, cell))
	end

	# allow other run_reactive! calls to be executed
	put!(notebook.executetoken, token)
	return new_topology
end

"See `run_reactive`."
function run_reactive_async!(notebook::Notebook, cells::Array{Cell, 1})::Task
	@async begin
		# because this is being run async, we need to catch exceptions manually
		try
			run_reactive!(notebook, cells)
		catch ex
			bt = stacktrace(catch_backtrace())
			showerror(stderr, ex, bt)
		end
	end
end

run_reactive!(notebook::Notebook, cell::Cell) = run_reactive!(notebook, [cell])
run_reactive_async!(notebook::Notebook, cell::Cell) = run_reactive_async!(notebook, [cell])

"Run a single cell non-reactively, return whether the run was Interrupted."
function run_single!(notebook::Notebook, cell::Cell)::Bool
	starttime = time_ns()
	run = WorkspaceManager.eval_fetch_in_workspace(notebook, cell.parsedcode)
	cell.runtime = time_ns() - starttime

	if run.errored
		cell.output_repr = nothing
		cell.error_repr = run.output_formatted[1]
		cell.repr_mime = run.output_formatted[2]
	else
		cell.output_repr = run.output_formatted[1]
		cell.error_repr = nothing
		cell.repr_mime = run.output_formatted[2]
	end

	return run.interrupted
	# TODO: capture stdout and display it somehwere, but let's keep using the actual terminal for now
end

"Update a single cell's cache - parsed code etc"
function start_cache!(notebook::Notebook, cell::Cell)
	cell.parsedcode = Meta.parse(cell.code, raise=false)
	cell.module_usings = ExploreExpression.compute_usings(cell.parsedcode)
	cell.symstate = try
		ExploreExpression.compute_symbolreferences(cell.parsedcode)
	catch ex
		@error "Expression explorer failed on: " cell.code
		showerror(stderr, ex, stacktrace(backtrace()))
		ExploreExpression.SymbolsState()
	end
end

"Account for globals referenced in function calls by including `SymbolsState`s from called functions in the cell itself."
function finish_cache!(notebook::Notebook, cell::Cell)
	calls = all_recursed_calls!(notebook, cell.symstate)
	calls = union(calls, keys(cell.symstate.funcdefs)) # _assume_ that all defined functions are called inside the cell.
	filter!(calls) do call
		call in keys(notebook.combined_funcdefs)
	end
	cell.symstate.references = union(cell.symstate.references, (notebook.combined_funcdefs[func].references for func in calls)...)
	cell.symstate.assignments = union(cell.symstate.assignments, (notebook.combined_funcdefs[func].assignments for func in calls)...)
end


text:
@testset "High-level tools" begin

    W = addprocs(2)
    @everywhere using DistributedData

    Random.seed!(1)
    dd = rand(11111, 5)
    buckets = rand([1, 2, 3], 11111)

    di1 = scatter_array(:test1, dd, W[1:1])
    di2 = scatter_array(:test2, dd, W)
    buckets1 = scatter_array(:buckets1, buckets, W[1:1])
    buckets2 = scatter_array(:buckets2, buckets, W)

    @testset "Distribution works as expected" begin
        @test gather_array(di1) == dd
        @test gather_array(di2) == dd
    end

    subsel = [1, 3, 4]
    dd = dd[:, subsel]
    dselect(di1, subsel)
    dselect(di2, subsel)

    @testset "dapply_cols, dapply_rows" begin
        diV = scatter_array(:testV, dd, W[1:1])
        diH = scatter_array(:testH, Matrix(dd'), W[1:1])
        dapply_cols(diV, (x, _) -> x ./ sum(x), Vector(1:3))
        dapply_rows(diH, x -> x ./ sum(x))
        @test isapprox(gather_array(diV), Matrix(gather_array(diH)'))
        unscatter(diV)
        unscatter(diH)
    end

    @testset "dselect" begin
        @test gather_array(di1) == dd
        @test gather_array(di2) == dd
    end

    di3 = dcopy(di2, :test3)

    @testset "dcopy" begin
        @test di2.workers == di3.workers
        @test dd == gather_array(di3)
    end

    unscatter(di3)

    (means1, sds1) = dstat(di1, [1, 3])
    (means2, sds2) = dstat(di2, [1, 3])

    @testset "dstat" begin
        @test length(means1) == 2
        @test length(sds1) == 2
        @test isapprox(means1, means2, atol = 1e-8)
        @test isapprox(sds1, sds2, atol = 1e-8)
        @test all([isapprox(0.5, m, atol = 0.1) for m in means1])
        @test all([isapprox(0.29, s, atol = 0.1) for s in sds1])
    end

    medians1 = dmedian(di1, [1, 3])
    medians2 = dmedian(di2, [1, 3])

    @testset "dmedian" begin
        @test length(medians1) == 2
        @test isapprox(medians1, medians2, atol = 1e-4)
        @test all([isapprox(0.5, m, atol = 0.1) for m in medians1])
    end

    (bmeans1, bsds1) = dstat_buckets(di1, 3, buckets1, [1, 3])
    (bmeans2, bsds2) = dstat_buckets(di2, 3, buckets2, [1, 3])

    @testset "bucketed dstat" begin
        @test size(bmeans1) == (3, 2)
        @test size(bsds1) == (3, 2)
        @test isapprox(bmeans1, bmeans2, atol = 1e-8)
        @test isapprox(bsds1, bsds2, atol = 1e-8)
        @test all([isapprox(0.5, m, atol = 0.1) for m in bmeans1])
        @test all([isapprox(0.29, s, atol = 0.1) for s in bsds1])
    end

    medians1 = dmedian_buckets(di1, 3, buckets1, [1, 3])
    medians2 = dmedian_buckets(di2, 3, buckets2, [1, 3])

    @testset "dmedian" begin
        @test size(medians1) == (3, 2)
        @test isapprox(medians1, medians2, atol = 1e-4)
        @test all([isapprox(0.5, m, atol = 0.1) for m in medians1])
    end

    dscale(di1, [2, 3])
    dscale(di2, [2, 3])
    (means1, sds1) = dstat(di1, [1, 2, 3])
    (means2, sds2) = dstat(di2, [1, 2, 3])

    @testset "dscale" begin
        @test isapprox(means1, means2, atol = 1e-8)
        @test isapprox(sds1, sds2, atol = 1e-8)
        @test isapprox(0.5, means1[1], atol = 0.1)
        @test all([isapprox(0.0, m, atol = 1e-3) for m in means1[2:3]])
        @test all([isapprox(1.0, s, atol = 1e-3) for s in sds1[2:3]])
    end

    unscatter(di1)
    unscatter(di2)
    rmprocs(W)

    d = ones(Float64, 2, 2)
    di = scatter_array(:test, d, [1])
    dscale(di, [1, 2])
    @testset "dscale does not produce NaNs on sdev==0" begin
        @test gather_array(di) == zeros(Float64, 2, 2)
    end
    unscatter(di)

end

ext: jl


text:
function dryair_check(p,t)
    if !(173.1 < t < 473.2)
        return DomainError("temperature should be between 173.15 K and 473.15 K")
    elseif !(0.0 <= p < 5e6)
        return DomainError("Pressure should be below 5×10^5 Pa")
    else
        return nothing
    end
end
function mass_volume_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return volumeair(t, p)
end

function mol_volume_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return molarvolumeair(t, p)
end

function mass_density_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    return one(p)/mass_volume_impl(mt,model,p,t)
end

function mol_density_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    return one(p)/mol_volume_impl(mt,model,p,t)
end

function mass_enthalpy_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return enthalpyair(t, p)
end

function mol_enthalpy_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return molarenthalpyair(t, p)
end

function mass_entropy_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return entropyair(t, p)
end

function mol_entropy_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return molarentropyair(t, p)
end

function compressibility_factor_impl(mt::SinglePT,model::ASHRAEDryAir,p, t)
    dryair_check(p,t)
    return p*molarvolumeair(t, p) / (R*t)
end

function water_vapor_check(t)
    if !(173.1 < t < 473.2)
        return DomainError("temperature should be between 173.15 K and 473.15 K")
    else
        return nothing
    end
end
function mass_volume_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return volumevapor(t)
end

function mol_volume_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return molarvolvapor(t)
end

function mass_density_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return one(t)/volumevapor(t)
end

function mol_density_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return one(t)/molarvolvapor(t)
end

function mass_enthalpy_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    enthalpyvapor(t)
end

function mol_enthalpy_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return enthalpyvapor(t)*Mv
end

function mol_entropy_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return molarentropyvapor(t)
end

function mass_entropy_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return entropyvapor(t)
end

function compressibility_factor_impl(mt::SingleSatT, model::ASHRAEWaterVapor,t)
    water_vapor_check(t)
    return Zvapor(t)
end

#=
,WebBulbPT  => :hum_wetbulb
,HumidityRatioPT => :hum_ratio
,MolHumidityPT => :hum_molfrac
,MassHumidityPT => :hum_massfrac
,RelativeHumidityPT => :rel_hum
,DewPointPT => :hum_dewpoint

=#

const HumPT = Tuple{Pressure,Temperature,HumiditySpec}
const MolHumidityPT = Tuple{Pressure,Temperature,HumiditySpec{MolarHumidity}}
const MassHumidityPT = Tuple{Pressure,Temperature,HumiditySpec{MassHumidity}}
const HumidityRatioPT = Tuple{Pressure,Temperature,HumiditySpec{HumidityRatio}}
const RelativeHumidityPT = Tuple{Pressure,Temperature,HumiditySpec{RelativeHumidity}}
const DewPointPT = Tuple{Pressure,Temperature,HumiditySpec{HumidityDewPoint}}
const WebBulbPT = Tuple{Pressure,Temperature,HumiditySpec{WetBulbTemperature}}

"mixture molecular weight of moist air"
moist_mw(xv) = xv*Mv + (1-xv)*Ma
#for interaction with ThermoState
function hum_molfrac_impl(mt::MultiPT,model::ASHRAEMoistAir,p,t,x)
    return last(x)
end

function hum_molfrac_impl(mt::MolHumidityPT,model::ASHRAEMoistAir,p,t,xv)
    return xv
end

function hum_molfrac_impl(mt::HumidityRatioPT,model::ASHRAEMoistAir, p,t,hum)
    return hum  / (Mv/Ma + hum)
end

function hum_molfrac_impl(mt::MassHumidityPT,model::ASHRAEMoistAir,p,t,hum)
    return hum * Ma / (Mv + hum*(Ma - Mv))
end

function hum_molfrac_impl(mt::RelativeHumidityPT,model::ASHRAEMoistAir,p,t,hum)
    return hum * efactor(t, p) * Pws(t) / p
end

function hum_molfrac_impl(mt::DewPointPT,model::ASHRAEMoistAir,p,t,hum)
    return efactor(hum,p) * Pws(hum) / p
end

function hum_molfrac_impl(mt::WebBulbPT,model::ASHRAEMoistAir,p,t,B)
    w = calc_W_from_B(t, B, p)
    return w / (Mv/Ma + w)
end

function dry_volume_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return volumemoist(t, p, xv)
end

function mass_volume_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    mw = moist_mw(xv)
    vm = molarvolumemoist(t, p, xv)
    return vm/mw
end

function mol_volume_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return molarvolumemoist(t, p, xv)
end

function mol_density_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    vm =  molarvolumemoist(t, p, xv)
    return one(vm)/vm
end

function mass_density_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    mw = moist_mw(xv)
    vm = molarvolumemoist(t, p, xv)
    return mw/vm
end

function dry_enthalpy_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return enthalpymoist(t, p, xv)
end

function mol_enthalpy_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return molarenthalpymoist(t, p, xv)
end

function mass_enthalpy_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    mw = moist_mw(xv)
    hn = molarenthalpymoist(t, p, xv)
    return hn/mw
end

function dry_entropy_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return entropymoist(t, p, xv)
end

function mass_entropy_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    mw = moist_mw(xv)
    sn = molarentropymoist(t, p, xv)
    return sn/mw
end

function mol_entropy_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return molarentropymoist(t, p, xv)
end

function compressibility_factor_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return Zmoist(t, p, xv)
end

function calcdewpoint(Tk, P, xv, EPS=1e-9, MAXITER=100)
    # Use Ideal Gas to
    D = Tws(xv*P)
    Dnew = D
    i = 0
    err = 0.0
    for i = 1:MAXITER
        f = efactor(D, P)
        Dnew = Tws(xv*P/f)
        err = abs(D-Dnew)
        if err < EPS
            return Dnew
        end
        D = Dnew
    end
    return Dnew
end

function hum_dewpoint_impl(mt::DewPointPT, model::ASHRAEMoistAir,p,t,hum)
    return hum
end

function hum_dewpoint_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    D = calcdewpoint(t, p, xv)
    return D
end

function rel_hum_impl(mt::RelativeHumidityPT, model::ASHRAEMoistAir,p,t,hum)
    return hum
end

function rel_hum_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    xv * p / (efactor(t, p) * Pws(t))
end

function hum_wetbulb_impl(mt::WebBulbPT, model::ASHRAEMoistAir,p,t,hum)
    return hum
end

function hum_wetbulb_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    B = calcwetbulb(t, p, xv)
    return B
end

_humrat(xv) = xv / (1-xv) * (Mv/Ma)

function hum_ratio_impl(mt::HumidityRatioPT, model::ASHRAEMoistAir,p,t,hum)
    return hum
end
function hum_ratio_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return _humrat(xv)
end

function hum_massfrac_impl(mt::MassHumidityPT, model::ASHRAEMoistAir,p,t,hum)
    return hum
end

function hum_massfrac_impl(mt::HumPT, model::ASHRAEMoistAir,p,t,hum)
    dryair_check(p,t)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return xv*Mv / ( (1-xv)*Ma + xv*Mv )
end

#implementations done, now with thermostate conections:

#extracts the value of the humidity spec
function humidity_spec(::FromState,st::ThermodynamicState)
    return value(get_spec(HumiditySpec,st))
end

function mol_fraction_impl(mt::HumPT,model::ASHRAEMoistAir,p,t,hum)
    xv = hum_molfrac_impl(mt,model,p,t,hum)
    return SVector(one(xv)-xv,xv)
end

function mass_fraction_impl(mt::HumPT,model::ASHRAEMoistAir,p,t,hum)
    xv = hum_massfrac_impl(mt,model,p,t,hum)
    return SVector(one(xv)-xv,xv)
end

#hum functions
for (hum_op,humunit) in zip(
    [:hum_dewpoint,:hum_massfrac,
    :hum_molfrac,:rel_hum,
    :hum_wetbulb,:hum_ratio,
    :dry_volume,:dry_enthalpy,:dry_entropy],
    [u"K",Unitful.NoUnits,
    Unitful.NoUnits,Unitful.NoUnits,
    u"K",Unitful.NoUnits,
    u"m^3/kg",u"J/kg",u"J/(kg*K)"])
    hum_op_impl = Symbol(hum_op,:_impl)

    @eval begin

        ThermoState.default_units(x::typeof($hum_op)) = $humunit

        function $hum_op(model::ASHRAEMoistAir,st::ThermodynamicState,unit=$humunit)
            return $hum_op(state_type(st),model,st,unit)
        end

        function $hum_op(mt::HumPT,model::ASHRAEMoistAir,st::ThermodynamicState,unit)
            hum = humidity_spec(FromState(),st)
            p = pressure(FromState(),st)
            t = temperature(FromState(),st)
            res  = $hum_op_impl(mt,model,p,t,hum)
            return convert_unit($humunit,unit,res)
        end

        function $hum_op(mt::MultiPT,model::ASHRAEMoistAir,st::ThermodynamicState,unit)
            x = mol_fraction(FromState(),st,nothing,molecular_weight(model))
            hum = last(x)
            _mt = MolHumidityPT() #MolHumidityPT is a complete state type, so it can be called directly instead of using QuickTypes
            p = pressure(FromState(),st)
            t = temperature(FromState(),st)
            res  = $hum_op_impl(_mt,model,p,t,hum)
            return convert_unit($humunit,unit,res)
        end
    end
end

for (op,_unit) in zip([:mol_volume,:mass_volume,
        :mol_density,:mass_density,
        :mol_entropy,:mass_entropy,
        :mol_enthalpy,:mass_enthalpy,
        :compressibility_factor],
        [u"m^3/mol",u"m^3/kg",
        u"mol/(m^3)",u"kg/(m^3)",
        u"J/(mol*K)",u"J/(kg*K)",
        u"J/mol",u"J/kg",
        Unitful.NoUnits])

    op_impl = Symbol(op,:_impl)
    @eval begin
        function $op(model::ASHRAEMoistAir,st::ThermodynamicState,unit=$_unit)
            return $op(state_type(st),model,st,unit)
        end

        function $op(model::ASHRAEWaterVapor,st::ThermodynamicState,unit=$_unit)
            return $op(state_type(st),model,st,unit)
        end

        function $op(model::ASHRAEDryAir,st::ThermodynamicState,unit=$_unit)
            return $op(state_type(st),model,st,unit)
        end

        function $op(mt::SinglePT,model::ASHRAEWaterVapor,st::ThermodynamicState,unit)
            p = pressure(FromState(),st)
            t = temperature(FromState(),st)
            res  = $op_impl(mt,model,p,t)
            return convert_unit($_unit,unit,res)
        end

        function $op(mt::SinglePT,model::ASHRAEDryAir,st::ThermodynamicState,unit)
            p = pressure(FromState(),st)
            t = temperature(FromState(),st)
            res  = $op_impl(mt,model,p,t)
            return convert_unit($_unit,unit,res)
        end

        function $op(mt::HumPT,model::ASHRAEMoistAir,st::ThermodynamicState,unit)
            hum = humidity_spec(FromState(),st)
            p = pressure(FromState(),st)
            t = temperature(FromState(),st)
            res  = $op_impl(mt,model,p,t,hum)
            return convert_unit($_unit,unit,res)
        end

        function $op(mt::MultiPT,model::ASHRAEMoistAir,st::ThermodynamicState,unit)
            x = mol_fraction(FromState(),st,nothing,molecular_weight(model))
            hum = last(x)
            _mt = MolHumidityPT() #MolHumidityPT is a complete state type, so it can be called directly instead of using QuickTypes
            p = pressure(FromState(),st)
            t = temperature(FromState(),st)
            res  = $op_impl(mt,model,p,t,hum)
            return convert_unit($_unit,unit,res)
        end
    end
end

#total functions
function total_volume(model::ASHRAEMoistAir,st::ThermodynamicState,unit=u"m^3")
    val = mass_volume(model,st)
    m = mass(FromState(),st,u"kg",molecular_weight(model))
    res = val*m
    return convert_unit(u"m^3",unit,res)
end

function total_entropy(model::ASHRAEMoistAir,st::ThermodynamicState,unit=u"J/K")
    val = mass_entropy(model,st)
    m = mass(FromState(),st,u"kg",molecular_weight(model))
    res = val*m
    return convert_unit(u"J/K",unit,res)
end

function total_enthalpy(model:: AbstractMoistAirModel,st::ThermodynamicState,unit=u"J")
    val = mass_enthalpy(model,st)
    m = mass(FromState(),st,u"kg",molecular_weight(model))
    res = val*m
    return convert_unit(u"J",unit,res)
end

function mol_gibbs(model::AbstractMoistAirModel,st::ThermodynamicState,unit=u"J/mol")
    h = mol_enthalpy(model,st)
    t = temperature(FromState(),st)
    s = mol_entropy(model,st)
    res = h-ts
    return convert_unit("J/mol",unit,res)
end

function mass_gibbs(model::AbstractMoistAirModel,st::ThermodynamicState,unit=u"J/kg")
    h = mass_enthalpy(model,st)
    t = temperature(FromState(),st)
    s = mass_entropy(model,st)
    res = h-ts
    return convert_unit("J/kg",unit,res)
end

function total_gibbs(model::AbstractMoistAirModel,st::ThermodynamicState,unit=u"J")
    h = mol_enthalpy(model,st)
    t = temperature(FromState(),st)
    s = mol_entropy(model,st)
    n = moles(FromState,st,u"mol",molecular_weight(model))
    res = n*(h-ts)
    return convert_unit("J",unit,res)
end

function mol_fraction(model::ASHRAEMoistAir,st::ThermodynamicState)
    xv = hum_molfrac(model,st)
    return SVector(1-xv,xv)
end

function mass_fraction(model::ASHRAEMoistAir,st::ThermodynamicState)
    xv = hum_massfrac(model,st)
    return SVector(1-xv,xv)
end

function mol_number(model::ASHRAEMoistAir,st::ThermodynamicState)
    xv = hum_molfrac(model,st)
    n = moles(FromState(),st,u"mol",molecular_weight(model))
    return SVector(n*(1-xv),n*xv)
end

function mass_number(model::ASHRAEMoistAir,st::ThermodynamicState)
    xv = hum_massfrac(model,st)
    n = moles(FromState(),st,u"mol",molecular_weight(model))
    mw = moist_mw(xv)
    m = n*mw
    return SVector(m*(1-xv),m*xv)
end


text:
### MATERIAL STATES ###
"""
`MaterialState` is the abstract supertype of objets storing the stresses and strains states of the material at the location of one quadrature point.
The states of the whole material are then stored in a `Vector{Vector{<:MaterialState}} (`MaterialState`s contained in a cell, and the cell in a grid).
Stores stress and strain informations as tensorial quantities.
"""
abstract type MaterialState{S <: Real} end

function update_material_state!(model)
    for cell_states in model.material_state
        foreach(update_state!, cell_states)
    end
end

function copy_temp_state!(model,mp)
    for (cell_id, cell_states) in enumerate(model.material_state)
        for (qp_id,state) in enumerate(cell_states)
            for field in propertynames(state)
                str_field = String(field)
                if length(str_field) > 4
                    (str_field[1:4] == "temp") && setproperty!(state,field,getproperty(mp[cell_id][qp_id],field))
                end
            end
        end
    end
end
##### Elastic or viscous #####
"""
`BasicMaterialState` is used for visco-elastic materials. It stores the elastic strains and stresses.

# Fields
ϵ : strain tensor
σ : stress tensor
temp_ϵ : temporary strain tensor
temp_σ : temporary stress tensor
"""
mutable struct BasicMaterialState{T} <: MaterialState{T}
    # Store converged values
    ϵ::SymmetricTensor{2, 3, T} # total strain
    σ::SymmetricTensor{2, 3, T} # stress

    # Store temporary values used during equilibrium iterations
    temp_ϵ::SymmetricTensor{2, 3, T}
    temp_σ::SymmetricTensor{2, 3, T}
end

function BasicMaterialState()
    return BasicMaterialState(
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}))
end

function update_state!(state::BasicMaterialState)
    state.ϵ = state.temp_ϵ
    state.σ = state.temp_σ
end

##### Plastic #####
"""
`PlasticMaterialState` is used for visco-elasto-plastic materials. It stores the elastic strains and stresses as well as the plastic strains and the accumulated plastic strain.

# Fields
ϵᵖ : plastic strain tensor
ϵ̅ᵖ : accumulated plastic strain
ϵ  : strain tensor
σ  : stress tensor
temp_ϵᵖ : temporary plastic strain tensor
temp_ϵ̅ᵖ : temporary accumulated plastic strain
temp_ϵ : temporary strain tensor
temp_σ : temporary stress tensor
"""
mutable struct PlasticMaterialState{T} <: MaterialState{T}
    # Store converged values
    ϵᵖ::SymmetricTensor{2, 3, T} # plastic strain
    ϵ̅ᵖ::T # accumulated plastic strain
    ϵ::SymmetricTensor{2, 3, T} # total strain
    σ::SymmetricTensor{2, 3, T} # stress

    # Store temporary values used during equilibrium iterations
    temp_ϵᵖ::SymmetricTensor{2, 3, T}
    temp_ϵ̅ᵖ::T
    temp_ϵ::SymmetricTensor{2, 3, T}
    temp_σ::SymmetricTensor{2, 3, T}
end
function PlasticMaterialState{T}() where {T<:Real}
    return PlasticMaterialState(
                zero(SymmetricTensor{2, 3, T}),
                zero(T),
                zero(SymmetricTensor{2, 3, T}),
                zero(SymmetricTensor{2, 3, T}),
                zero(SymmetricTensor{2, 3, T}),
                zero(T),
                zero(SymmetricTensor{2, 3, T}),
                zero(SymmetricTensor{2, 3, T}))
end
function PlasticMaterialState()
    return PlasticMaterialState(
                zero(SymmetricTensor{2, 3}),
                0.0,
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                0.0,
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}))
end

function update_state!(state::PlasticMaterialState)
    state.ϵᵖ = state.temp_ϵᵖ
    state.ϵ̅ᵖ = state.temp_ϵ̅ᵖ
    state.ϵ = state.temp_ϵ
    state.σ = state.temp_σ
end
function add_state!(stateout::PlasticMaterialState,statein::PlasticMaterialState)
    stateout.ϵᵖ += statein.ϵᵖ
    stateout.ϵ̅ᵖ += statein.ϵ̅ᵖ
    stateout.ϵ += statein.ϵ
    stateout.σ += statein.σ
end

function set_temp_state!(r::Rheology{TT,TD,Nothing,TE,TP},clock,state,σ,ϵ) where {TT,TD,TE,TP}
    Dᵉ = r.elasticity.Dᵉ
    Δϵ = ϵ - state.ϵ
    invDᵉ = get_elastic_compliance_tensor(r)
    Δϵᵖ = Δϵ - invDᵉ⊡(σ - state.σ)
    ϵ̇ᵖ = Δϵᵖ/clock.Δt
    Δϵ̅ᵖ = sqrt((2/3) * ϵ̇ᵖ ⊡ ϵ̇ᵖ) * clock.Δt

    state.temp_ϵᵖ = state.ϵᵖ + Δϵᵖ
    state.temp_ϵ̅ᵖ = state.ϵ̅ᵖ + Δϵ̅ᵖ
    state.temp_σ = σ
    state.temp_ϵ = ϵ
end
##### Damaged plastic #####
"""
`DamagedPlasticMaterialState` is used for damaged visco-elasto-plastic materials.
It stores the elastic strains and stresses as well as the plastic strains and the accumulated plastic strain and a damage state variable.

# Fields
D  : damage state variable
ϵᵖ : plastic strain tensor
ϵ̅ᵖ : accumulated plastic strain
ϵ  : strain tensor
σ  : stress tensor
temp_D  : temporary damage state variable
temp_ϵᵖ : temporary plastic strain tensor
temp_ϵ̅ᵖ : temporary accumulated plastic strain
temp_ϵ  : temporary strain tensor
temp_σ  : temporary stress tensor
"""
mutable struct DamagedPlasticMaterialState{T} <: MaterialState{T}
    # Store converged values
    D::T
    ϵᵖ::SymmetricTensor{2, 3, T} # plastic strain
    ϵ̅ᵖ::T # cumulated plastic strain
    ϵ::SymmetricTensor{2, 3, T} # total strain
    σ::SymmetricTensor{2, 3, T} # stress

    # Store temporary values used during equilibrium iterations
    temp_D::T
    temp_ϵᵖ::SymmetricTensor{2, 3, T}
    temp_ϵ̅ᵖ::T
    temp_ϵ::SymmetricTensor{2, 3, T}
    temp_σ::SymmetricTensor{2, 3, T}
end

function DamagedPlasticMaterialState(r::Rheology{T}) where {T}
    return DamagedPlasticMaterialState(
                r.damage.D₀,
                zero(SymmetricTensor{2, 3, T}),
                zero(T),
                zero(SymmetricTensor{2, 3, T}),
                zero(SymmetricTensor{2, 3, T}),
                r.damage.D₀,
                zero(SymmetricTensor{2, 3, T}),
                zero(T),
                zero(SymmetricTensor{2, 3, T}),
                zero(SymmetricTensor{2, 3, T}))
end

function update_state!(state::DamagedPlasticMaterialState)
    state.D = state.temp_D
    state.ϵᵖ = state.temp_ϵᵖ
    state.ϵ̅ᵖ = state.temp_ϵ̅ᵖ
    state.ϵ = state.temp_ϵ
    state.σ = state.temp_σ
end


text:
sourcedir = joinpath(dirname(@__FILE__), "src", "word2vec-c")
cd(sourcedir)
run(`make`)


text:
using DataStructures
using Random: rand, AbstractRNG
using Distributions: Uniform, Exponential, params

"""
Classic Direct method for exponential transitions.
This doesn't do any caching of rates.
"""
struct MarkovDirect
end

function next(rm::MarkovDirect, process, when, rng::AbstractRNG)
    total = 0.0
    cumulative = zeros(Float64, 0)
    keys = Array{Any,1}()
    hazards(process, rng) do clock, distribution::Exponential, enabled::Bool
        if enabled
            total += params(distribution)[1]
            push!(cumulative, total)
            push!(keys, clock)
        end
    end

    if total > eps(Float64)
        chosen = searchsortedfirst(cumulative, rand(rng, Uniform(0, total)))
        @assert chosen < length(cumulative) + 1
        return (when - log(rand(rng)) / total, keys[chosen])
    else
        return (Inf, nothing)
    end
end

Observer(fr::MarkovDirect) = (hazard, time, updated, rng) -> nothing

struct DirectCall{T}
    key::Dict{T, Int64}
    propensity::Vector{Float64}
end

function DirectCall(::Type{T}) where {T}
    DirectCall{T}(Dict{T, Int64}(), zeros(Float64, 0))
end

function set_clock!(dc::DirectCall{T}, clock::T, distribution::Exponential,
        enabled::Symbol, rng::AbstractRNG) where {T}
    if enabled == :Enabled
        hazard = params(distribution)[1]
        if !haskey(dc.key, clock)
            dc.key[clock] = length(push!(dc.propensity, hazard))
        else
            dc.propensity[dc.key[clock]] = hazard
        end
    elseif enabled == :Changed  # Why is == not found otherwise?
        dc.propensity[dc.key[clock]] = params(distribution)[1]
    else  # else it's disabled.
        dc.propensity[dc.key[clock]] = 0.0
    end
end

function next(dc::DirectCall, when::Float64, rng::AbstractRNG)
    total = sum(dc.propensity)
    if total > eps(Float64)
        chosen = searchsortedfirst(cumsum(dc.propensity), rand(rng, Uniform(0, total)))
        @assert chosen < length(dc.propensity) + 1
        key_name = [x for (x, y) in pairs(dc.key) if y == chosen][1]
        return (when - log(rand(rng)) / total, key_name)
    else
        return (Inf, nothing)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Gnome_themes_extra"
version = v"3.28"

# Collection of sources required to build gnome-themes-extra
sources = [
    "http://ftp.gnome.org/pub/gnome/sources/gnome-themes-extra/$(version.major).$(version.minor)/gnome-themes-extra-$(version.major).$(version.minor).tar.xz" =>
    "7c4ba0bff001f06d8983cfc105adaac42df1d1267a2591798a780bac557a5819",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gnome-themes-extra-*/
apk add intltool

# Clear out `.la` files since they're often wrong and screw us up
rm -f ${prefix}/lib/*.la

FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    FLAGS+=(ac_cv_path_GTK_UPDATE_ICON_CACHE=gtk-update-icon-cache.exe)
fi

if [[ "${target}" == powerpc64le-* ]]; then
    export CFLAGS="-Wl,-rpath-link,${prefix}/lib64"
fi

./configure --prefix=${prefix} --host=${target} \
    --disable-gtk2-engine \
    "${FLAGS[@]}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
#platforms = supported_platforms()

# Limit to the same platforms as Gtk for now
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    FileProduct("share/themes/Adwaita/index.theme", :adwaita_index),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "GTK3_jll",
    "Librsvg_jll",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using InformativePrior
using Test

@testset "InformativePrior.jl" begin
    # Write your tests here.
end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
import GeometricFlux: message, update, propagate

in_channel = 10
out_channel = 5
N = 6
adj = [0. 1. 0. 0. 0. 0.;
       1. 0. 0. 1. 1. 1.;
       0. 0. 0. 0. 0. 1.;
       0. 1. 0. 0. 1. 0.;
       0. 1. 0. 1. 0. 1.;
       0. 1. 1. 0. 1. 0.]
ne = [[2], [1,4,5,6], [6], [2,5], [2,4,6], [2,3,5]]

struct NewLayer <: MessagePassing
    adjlist::AbstractVector{<:AbstractVector}
    weight
end

NewLayer(adjm::AbstractMatrix, m, n) = NewLayer(neighbors(adjm), randn(m,n))

(l::NewLayer)(X) = propagate(l, X=X, aggr=:add)
message(n::NewLayer; x_i=zeros(0), x_j=zeros(0)) = x_j
update(::NewLayer; X=zeros(0), M=zeros(0)) = M

X = Array(reshape(1:N*in_channel, in_channel, N))
l = NewLayer(adj, in_channel, out_channel)

@testset "Test neighboring" begin
    Y, clst = neighboring(l; X=X)
    @test adjlist(l) == ne
    @test size(Y) == (in_channel, 14)
    @test Y[1, :] == [11, 1, 31, 41, 51, 51, 11, 41, 11, 31, 51, 11, 21, 41]
    @test clst == [1, 2, 2, 2, 2, 3, 4, 4, 5, 5, 5, 6 ,6, 6]
end

message(n::NewLayer; x_i=zeros(0), x_j=zeros(0)) = n.weight' * x_j

@testset "Test MessagePassing layer" begin
    Y = l(X)
    @test size(Y) == (out_channel, N)
end

in_channel = 100
X = Array(reshape(1:N*in_channel, in_channel, N))
l = NewLayer(adj, in_channel, out_channel)

@testset "Test multi-thread MessagePassing layer" begin
    Y = l(X)
    @test size(Y) == (out_channel, N)
end


text:
######################################################################
# corner_initializer.jl
######################################################################

"""
`ci = CornerInitializer()`

Takes a trajectory to the farthest corner.
"""
mutable struct CornerInitializer <: Initializer end

# Assumes we our domain is 2d
# corner 1 = 0,0
# corner 2 = 1,0
# corner 3 = 0,1
# corner 4 = 1,1
function initialize(ci::CornerInitializer, em::ErgodicManager, tm::TrajectoryManager)

	# dimensionality of state space
	x0 = tm.x0
	n = length(x0)

	# first corner: 0,0
	dx1 = 0.0 - x0[1]
	dy1 = 0.0 - x0[2]
	dc1 = sqrt(dx1*dx1 + dy1*dy1)

	# second corner: 1,0
	dx2 = 1.0 - x0[1]
	dy2 = 0.0 - x0[2]
	dc2 = sqrt(dx2*dx2 + dy2*dy2)

	# third corner: 0,1
	dx3 = 0.0 - x0[1]
	dy3 = 1.0 - x0[2]
	dc3 = sqrt(dx3*dx3 + dy3*dy3)

	# fourth corner: 1,1
	dx4 = 1.0 - x0[1]
	dy4 = 1.0 - x0[2]
	dc4 = sqrt(dx4*dx4 + dy4*dy4)

	bi = indmax([dc1,dc2,dc3,dc4])

	bc = [0.,0.]
	if bi == 2
		bc = [1.,0.]
	elseif bi == 3
		bc = [0.,1.]
	elseif bi == 4
		bc = [1.,1.]
	end

	# could now do point initializer to the selected corner
	return initialize(PointInitializer(bc), em, tm)
end


text:
using Test

m = MimiPAGE2009.get_model()
run(m)

m = MimiPAGE2009.get_model("policy-b")
run(m)

while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
    run(m)
end

# Climate component
temp = m[:ClimateTemperature, :rt_g_globaltemperature]
temp_compare = readpagedata(m,"test/validationdata/policy-b/rt_g_globaltemperature.csv")
@test temp ≈ temp_compare rtol=1e-4

# Abatement Costs
abatement = m[:TotalAbatementCosts, :tct_totalcosts]
abatement_compare = readpagedata(m,"test/validationdata/policy-b/tct_totalcosts.csv")
@test abatement ≈ abatement_compare rtol=1e-2

te = m[:EquityWeighting, :te_totaleffect]
te_compare = 224812034.49558273
@test te ≈ te_compare rtol=1e4


text:
using PiecewiseDeterministicMarkovProcesses, Test, LinearAlgebra, Random, DifferentialEquations

macro testS(label, args...)
	:(@testset $label begin @test $(args...); end)
end

include("simpleTests.jl")

@testset "Example TCP" begin
	include("../examples/tcp.jl")
	@test norm(errors[6:end], Inf64) < 1e-4
end

@testset "Example explosion of jump times" begin
	include("../examples/pdmpExplosion.jl")
	@test norm(errors[1:6], Inf64) < 1e-4
end

@testset "Example with stiff ODE part" begin
	include("pdmpStiff.jl")
	@test norm(errors, Inf64) < 1e-3
	@testS "Call many times the same problem" restime1 == res12.time
end

@testset "Controlling allocations" begin
	# @test alloc1 == alloc2
end

@testset "Test Rate structures 1/2" begin
	include("testRatesCst.jl")
end

@testset "Test Rate structures 2/2" begin
	include("testRatesComposite.jl")
end

@testset "Example with 2d example, for autodiff" begin
	include("../examples/tcp2d.jl")
end

@testset "Rejection method" begin
	include("../examples/tcp_rejection.jl")
	@test norm(errors, Inf64) < 1e-5
end

@testset "Test number of fictitious jumps" begin
	@test res1.njumps == res2.njumps
	@test res1.nrejected == res2.nrejected
end

@testset "Neuron model" begin
	include("../examples/pdmp_example_eva.jl")
	@test result1.time[end] == 100.
	@test result1.xd[2,end] == 113
end

@testset "Example SIR" begin
	include("../examples/sir.jl")
	@test result.xd[1,end] == 0
	@test result.xd[2,end] == 36
	@test result.xd[3,end] == 73
end

@testset "Example SIR(rejection)" begin
	include("../examples/sir-rejection.jl")
	@test result.xd[1,end] == 0
	@test result.xd[2,end] == 33
	@test result.xd[3,end] == 76
end

@testset "Neural network" begin
	include("../examples/neuron_rejection_exact.jl")
	@test result.xd[1,end] == 97
	@test size(result.xd)[1] == 100
end

@testset "DiffEqJump Wrap" begin
	include("../examples/examplediffeqjumpwrapper.jl")
end


text:
"""
    op = JopReal(dom)

Extract the real part of a complex input array where `dom::JetSpace{<:Complex}`.
"""
JopReal(dom::JetSpace{Complex{T}}) where {T} = JopLn(dom = dom, rng = JetSpace(T,size(dom)), df! = JopReal_df!, df′! = JopReal_df′!)
export JopReal

JopReal_df!(d::AbstractArray{T}, m::AbstractArray{Complex{T}}; kwargs...) where {T} = d .= real.(m)
JopReal_df′!(m::AbstractArray{Complex{T}}, d::AbstractArray{T}; kwargs...) where {T} = m .= d .+ 0 * im


text:
# this file contains some basic extensions to the FastChebInterp package
# https://github.com/stevengj/FastChebInterp.jl (MIT license). I am keeping them
# here for the moment, but at some point a PR may make sense. Because the IFGF
# package needs to evaluate Cheb polynomials of relatively low order, and
# because the order of the chebpoly the same for all nodes in the tree, it makes
# sense to:
# (a) pass the size statically so that evaluate can skip some conditionals and
# unroll loops. This yields significantly speed-ups in the evaluation of a
# chebpoly for given point
# (b) precompute a plan for the cheb transforms once and pass it to chebinterp.
# In fact the orders are small enough that the plan should just be a "direct
# summation" algorithm.
# (c) implement an in-place `chebinterp!` and `chebcoeffs!` which will use the
# array of values to store the chebyshev coefficiets
# (d) statically compute the reference tensor-product chebnodes given the order

"""
    cheb2nodes(n,[domain])
    cheb2nodes(p::NTuple,[domain])
    cheb2nodes(::Val{P},[domain])

Return the `n` Chebyshev nodes of the second kind on the interval `[-1,1]`.

If passed a tuple `p`, return the tensor product of one-dimensinal nodes of size
`p`.

If passed a `Val{P}`, statically compute the tensor product nodes.

Finally, passing a `domain::HyperRectangle` shifts the interpolation nodes to
it.

```julia
lb = (0.,0.)
ub = (3.,2.)
domain = HyperRectangle(lb,ub)
cheb2nodes((10,5),domain) # generate `10×5` Chebyshev nodes on the `[0,3]×[0,2]` rectangle.
```
"""
function cheb2nodes(n)
    [cos((i-1)*π /(n-1)) for i in 1:n]
end

function cheb2nodes(p::NTuple)
    nodes1d = map(i->cheb2nodes(i),p)
    iter    = Iterators.product(nodes1d...)
    return iter
end

@generated function cheb2nodes(::Val{P}) where {P}
    nodes1d = map(i->cheb2nodes(i),P)
    iter    = Iterators.product(nodes1d...)
    nodes   = SVector.(collect(iter))
    if prod(P) < 1e3
        nodes  = SArray{Tuple{P...}}(nodes)
    end
    return :($nodes)
end

function cheb2nodes(p::Val{P},rec::HyperRectangle) where {P}
    refnodes = cheb2nodes(p)
    lc,hc    = low_corner(rec), high_corner(rec)
    c0       = (lc+hc)/2
    c1       = (hc-lc)/2
    map(x-> c0 + c1.*x,refnodes)
end

function cheb2nodes(p::NTuple,rec::HyperRectangle)
    refnodes = cheb2nodes(p)
    lc,hc    = low_corner(rec), high_corner(rec)
    c0       = (lc+hc)/2
    c1       = (hc-lc)/2
    map(x-> c0 + c1.*x,refnodes)
end

function chebcoefs!(vals::AbstractArray)
    plan = FFTW.plan_r2r!(copy(vals),FFTW.REDFT00)
    chebcoefs!(vals,plan)
end

function chebcoefs!(vals::AbstractArray{<:Number,N},plan::FFTW.FFTWPlan) where {N}
    FFTW.mul!(vals,plan,vals) # type-I DCT
    # renormalize the result to obtain the conventional
    # Chebyshev-polnomial coefficients
    s = size(vals)
    # coefs ./= prod(n -> 2(n-1), s)
    scale = 2^N/prod(n -> 2(n-1), s)
    @. vals *= scale
    for dim = 1:N
        I = CartesianIndices(ntuple(i -> i == dim ? (1:1) : (1:s[i]),N))
        @views @. vals[I] /= 2
        I = CartesianIndices(ntuple(i -> i == dim ? (s[i]:s[i]) : 1:s[i],N))
        @views @. vals[I] /= 2
    end
    return vals
end

function chebcoefs!(vals::AbstractArray{<:SVector{K}},plan) where {K}
    coefs = ntuple(i -> chebcoefs!([v[i] for v in vals],plan), Val{K}())
    copyto!(vals,SVector{K}.(coefs...))
end

@fastmath function chebeval(coefs,x::SVector{N,<:Real},rec::HyperRectangle,sz::Val{SZ}) where {N,SZ}
    x0 = @. (x - rec.low_corner) * 2 / (rec.high_corner - rec.low_corner) - 1
    return evaluate(x0, coefs, Val{N}(), 1, length(coefs),sz)
end

@fastmath function evaluate(x::SVector{N}, c, ::Val{dim}, i1, len, sz::Val{SZ}) where {N,dim,SZ}
    @inbounds n = SZ[dim]
    @inbounds xd = x[dim]
    if dim == 1
        @inbounds c₁ = c[i1]
        if n ≤ 2
            n == 1 && return c₁ + one(xd) * zero(c₁)
            return muladd(xd, c[i1+1], c₁)
        end
        @inbounds bₖ = muladd(2xd, c[i1+(n-1)], c[i1+(n-2)])
        @inbounds bₖ₊₁ = oftype(bₖ, c[i1+(n-1)])
        for j = n-3:-1:1
            @inbounds bⱼ = muladd(2xd, bₖ, c[i1+j]) - bₖ₊₁
            bₖ, bₖ₊₁ = bⱼ, bₖ
        end
        return muladd(xd, bₖ, c₁) - bₖ₊₁
    else
        Δi = len ÷ n # column-major stride of current dimension

        # we recurse downward on dim for cache locality,
        # since earlier dimensions are contiguous
        dim′ = Val{dim-1}()

        c₁ = evaluate(x, c, dim′, i1, Δi,sz)
        if n ≤ 2
            n == 1 && return c₁ + one(xd) * zero(c₁)
            c₂ = evaluate(x, c, dim′, i1+Δi, Δi, sz)
            return c₁ + xd*c₂
        end
        cₙ₋₁ = evaluate(x, c, dim′, i1+(n-2)*Δi, Δi,sz)
        cₙ = evaluate(x, c, dim′, i1+(n-1)*Δi, Δi,sz)
        bₖ = muladd(2xd, cₙ, cₙ₋₁)
        bₖ₊₁ = oftype(bₖ, cₙ)
        for j = n-3:-1:1
            cⱼ = evaluate(x, c, dim′, i1+j*Δi, Δi,sz)
            bⱼ = muladd(2xd, bₖ, cⱼ) - bₖ₊₁
            bₖ, bₖ₊₁ = bⱼ, bₖ
        end
        return muladd(xd, bₖ, c₁) - bₖ₊₁
    end
end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
data = readline()

function move(mov)
    if mov == '^'
        [-1, 0]
    elseif mov == '>'
        [0, 1]
    elseif mov == 'v'
        [1, 0]
    else
        [0, -1]
    end
end

function solve_from(data)
    curr_position = [0, 0]
    positions = Set()
    push!(positions, curr_position)
    for mov in data
        curr_position += move(mov)
        push!(positions, curr_position)
    end
    positions
end

function solve_1()
    solve_from(data) |> length
end

function solve_2()
    union(solve_from(data[1:2:end]), solve_from(data[2:2:end])) |> length
end

solve_1() |> println
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
function J_hypot_0(x0, x1)
    x2 = DIO_Undef
    x3 = DIO_Undef
    x4 = DIO_Undef

@label _1
    DIO_DecRef(x2)
    x2 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dff4f3360) #= <built-in function isinstance> =#, x0, @DIO_Obj(0x00007ffd6a132140) #= <class 'str'> =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x00007ffd6a1239d0) #= <class 'bool'> =#, x2))
    if 0x00007ffd6a123868 === reinterpret(UInt64, x3)
        @goto _2
    else
        @goto _7
    end

@label _10
    DIO_DecRef(x2)
    x2 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x0, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x1, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x2)
    x2 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffbf7c20) #= <built-in function add> =#, x2, x3))
    DIO_DecRef(x4)
    x4 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011d81cdb310) #= <built-in function sqrt> =#, x2))
    @DIO_Return DIO_IncRef(x4)

@label _2
    DIO_DecRef(x2)
    x2 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x00007ffd6a12efd0) #= <class 'int'> =#, x0))
    @goto _3

@label _3
    DIO_DecRef(x0)
    x0 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dff4f3360) #= <built-in function isinstance> =#, x1, @DIO_Obj(0x00007ffd6a132140) #= <class 'str'> =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x00007ffd6a1239d0) #= <class 'bool'> =#, x0))
    if 0x00007ffd6a123868 === reinterpret(UInt64, x3)
        @goto _4
    else
        @goto _6
    end

@label _4
    DIO_DecRef(x0)
    x0 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x00007ffd6a12efd0) #= <class 'int'> =#, x1))
    @goto _5

@label _5
    DIO_DecRef(x1)
    x1 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x2, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x0, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x1)
    x1 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffbf7c20) #= <built-in function add> =#, x1, x3))
    DIO_DecRef(x4)
    x4 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011d81cdb310) #= <built-in function sqrt> =#, x1))
    @DIO_Return DIO_IncRef(x4)

@label _6
    DIO_DecRef(x0)
    x0 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x2, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x1, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x0)
    x0 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffbf7c20) #= <built-in function add> =#, x0, x3))
    DIO_DecRef(x4)
    x4 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011d81cdb310) #= <built-in function sqrt> =#, x0))
    @DIO_Return DIO_IncRef(x4)

@label _7
    DIO_DecRef(x2)
    x2 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dff4f3360) #= <built-in function isinstance> =#, x1, @DIO_Obj(0x00007ffd6a132140) #= <class 'str'> =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x00007ffd6a1239d0) #= <class 'bool'> =#, x2))
    if 0x00007ffd6a123868 === reinterpret(UInt64, x3)
        @goto _8
    else
        @goto _10
    end

@label _8
    DIO_DecRef(x2)
    x2 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x00007ffd6a12efd0) #= <class 'int'> =#, x1))
    @goto _9

@label _9
    DIO_DecRef(x1)
    x1 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x0, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x3)
    x3 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffc01770) #= <built-in function pow> =#, x2, @DIO_Obj(0x0000011dff4a6950) #= 2 =#))
    DIO_DecRef(x1)
    x1 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011dffbf7c20) #= <built-in function add> =#, x1, x3))
    DIO_DecRef(x4)
    x4 = @DIO_ChkExcAndDecRefSubCall(Py_CallFunction(@DIO_Obj(0x0000011d81cdb310) #= <built-in function sqrt> =#, x1))
    @DIO_Return DIO_IncRef(x4)

@label except
    DIO_Return = Py_NULL

@label ret
    DIO_DecRef(x2)
    DIO_DecRef(x3)
    DIO_DecRef(x4)
    return DIO_Return
end
DIO.DIO_ExceptCode(::typeof(J_hypot_0)) = Py_NULL
const DOC_J_hypot_0 = Base.unsafe_convert(Cstring, "Top J_hypot_0( D0 : int, D1 : int ) {\nlabel _1:\n  D2 = Py_CallFunction(builtins.isinstance, D0 : int, str)\n  D3 : bool = Py_CallFunction(bool, D2)\n  if D3 : bool\n  then goto _2\n  else goto _7\nlabel _10:\n  D2 = Py_CallFunction(_operator.pow, D0 : int, 2)\n  D3 = Py_CallFunction(_operator.pow, D1 : int, 2)\n  D2 = Py_CallFunction(_operator.add, D2, D3)\n  D4 = Py_CallFunction(math.sqrt, D2)\n  return D4\nlabel _2:\n  D2 = Py_CallFunction(int, D0 : int)\n  goto _3\nlabel _3:\n  D0 = Py_CallFunction(builtins.isinstance, D1 : int, str)\n  D3 : bool = Py_CallFunction(bool, D0)\n  if D3 : bool\n  then goto _4\n  else goto _6\nlabel _4:\n  D0 = Py_CallFunction(int, D1 : int)\n  goto _5\nlabel _5:\n  D1 = Py_CallFunction(_operator.pow, D2, 2)\n  D3 = Py_CallFunction(_operator.pow, D0, 2)\n  D1 = Py_CallFunction(_operator.add, D1, D3)\n  D4 = Py_CallFunction(math.sqrt, D1)\n  return D4\nlabel _6:\n  D0 = Py_CallFunction(_operator.pow, D2, 2)\n  D3 = Py_CallFunction(_operator.pow, D1 : int, 2)\n  D0 = Py_CallFunction(_operator.add, D0, D3)\n  D4 = Py_CallFunction(math.sqrt, D0)\n  return D4\nlabel _7:\n  D2 = Py_CallFunction(builtins.isinstance, D1 : int, str)\n  D3 : bool = Py_CallFunction(bool, D2)\n  if D3 : bool\n  then goto _8\n  else goto _10\nlabel _8:\n  D2 = Py_CallFunction(int, D1 : int)\n  goto _9\nlabel _9:\n  D1 = Py_CallFunction(_operator.pow, D0 : int, 2)\n  D3 = Py_CallFunction(_operator.pow, D2, 2)\n  D1 = Py_CallFunction(_operator.add, D1, D3)\n  D4 = Py_CallFunction(math.sqrt, D1)\n  return D4\n}\n")

CFunc_J_hypot_0(self :: PyPtr, args::Ptr{PyPtr}, n::Py_ssize_t) = @DIO_MakePyFastCFunc(J_hypot_0, args, n, 2)
const CFuncPtr_J_hypot_0 = @cfunction(CFunc_J_hypot_0, PyPtr, (PyPtr, Ptr{PyPtr}, Py_ssize_t))
const PyMeth_J_hypot_0 = PyMethodDef(
    Base.unsafe_convert(Cstring, :hypot),
    CFuncPtr_J_hypot_0,
    METH_FASTCALL,
    DOC_J_hypot_0
)
const PyFunc_J_hypot_0 = PyCFunction_New(pointer_from_objref(PyMeth_J_hypot_0), Py_NULL)


text:
"""
    readlog(file)

Read standard LWPC `.log` files and return `distance`, `amplitude`, and `phase` vectors.

The amplitude and phase at the receiver is not returned by this function.
"""
function readlog(file)
    # Search for first and last line of data
    lines = readlines(file)
    firstdataline = findfirst(startswith.(lines, "  dist   amplitude  phase")) + 1
    lastdataline = findfirst(startswith.(lines, "nc nrpt bearng")) - 1
    numdatalines = lastdataline - firstdataline + 1

    # Read log file
    raw = CSV.File(file; datarow=firstdataline, limit=numdatalines,
                   delim=' ', ignorerepeated=true, header=false,
                   silencewarnings=true, threaded=false)

    # Stack the columns together
    dist = vcat(raw.Column1, raw.Column4, raw.Column7)
    amplitude = vcat(raw.Column2, raw.Column5, raw.Column8)
    phase = vcat(raw.Column3, raw.Column6, raw.Column9)

    # Clean up end of the last column
    delidxs = ismissing.(dist)  # assuming these are at the end
    delidxs[findfirst(delidxs)-1] = true  # last valid entry is at receiver distance
    deleteat!(dist, delidxs)
    deleteat!(amplitude, delidxs)
    deleteat!(phase, delidxs)

    # If phase gets above 9999 deg in log file, there is no space between amplitude and phase
    if count(ismissing.(amplitude)) != count(ismissing.(phase))
        for i in eachindex(phase)
            if ismissing(phase[i])
                phase[i] = parse(Float64, amplitude[i][end-9:end])  # phase is 10000.0000
                amplitude[i] = parse(Float64, amplitude[i][1:end-10])
            end
        end
        # Other elements in the same column will also be string type
        for i in eachindex(amplitude)
            if amplitude[i] isa String
                amplitude[i] = parse(Float64, amplitude[i])
            end
        end
    end

    return dist, amplitude, phase
end

"""
    isroot(x; atol=1e-2)

Return `true` if `x` is approximately equal to 0 with the absolute tolerance `atol`.
"""
isroot(x; atol=1e-2, kws...) = isapprox(x, 0; atol=atol, kws...)


text:
@testset "sums" begin

    A = rand(10, 10)
    B = rand(Bool, 5, 5)
    C = randn(10, 10)
    D = rand(1:10, 5, 5)

    KA = A ⊗ B
    KB = C ⊗ D

    K = KA + KB
    Kdense = kron(A, B) .+ kron(C, D)

    v = randn(50)
    V = randn(50, 3)

    @test K ≈ Kdense
    @test size(K) == size(Kdense)
    @test collect(K) ≈ Kdense
    @test sum(K) ≈ sum(Kdense)
    @test sum(K, dims=1) ≈ sum(Kdense, dims=1)
    @test K' ≈ Kdense'
    @test transpose(K) ≈ transpose(Kdense)
    @test conj(K) ≈ conj(Kdense)

    @test K * v ≈ Kdense * v
    @test K * V ≈ Kdense * V

    @test 5 * K ≈ 5 * Kdense
    @test K * 2.1 ≈ Kdense * 2.1

    # recursive, kronecker products of Kronecker sums

    Kl = (KA + KB) ⊗ C
    Kr = C ⊗ (KA + KB)
    Kb = ((B ⊗ C) + (B ⊗ C)) ⊗ ((B ⊗ C) + (B ⊗ C))

    #=
    Kl = kronecker(KA + KB, C)
    Kr = kronecker(C, KA + KB)
    Kb = kronecker(kronecker(B, C) + kronecker(B, C), kronecker(B, C) + kronecker(B, C))
    =#

    n = size(Kl, 2)
    m = size(Kb, 2)

    v = randn(n)
    u = randn(m)

    @test Kl * v ≈ collect(Kl) * v
    @test Kr * v ≈ collect(Kr) * v
    @test Kb * u ≈ collect(Kb) * u

end


text:
using ImageFiltering
using Images

export rotation_matrix, ellipse_to_points, kernelize

function rotation_matrix(angle::Real; ccw=true)
    #= Helper for rotating points in 2d space

    Builds the rotation matrix specified by the angle and the direction of rotation

    Args :
        angle : angle in radians of rotation
        ccw : whether the angle is measured counter clockwise wrt the positive x axis

    Returns :
        A 2x2 rotation matrix
    =#

    rotate_mat = zeros(2,2);
    if ccw
        rotate_mat = [cos(angle) -sin(angle);
                      sin(angle) cos(angle)];
    else
        rotate_mat = [cos(angle) sin(angle);
                      -sin(angle) cos(angle)];
    end

    return rotate_mat;
end

# might add method for ConicFormEllipse or unite both Conic and Parametric in another structure
function ellipse_to_points(ellipse::ParametricFormEllipse; n=1000::Int)

    theta_plot_vals = range(0, 2*pi, length=n)

    unit_circle = [cos.(theta_plot_vals) sin.(theta_plot_vals)]

    onaxis_ellipse = reshape(ellipse.axis_lengths ./ 2, (1,2)) .* unit_circle

    U = rotation_matrix(deg2rad(ellipse.ccw_angle))

    rotated_ellipse = (U * onaxis_ellipse')'

    shifted_ellipse = vec(ellipse.center)' .+ rotated_ellipse

    return shifted_ellipse
end

function kernelize(ellipse::ParametricFormEllipse) #; gaussian = true)
    X = round.(Int, ellipse_to_points(ellipse, n = 1000))

# TODO change translateImage so that the new image is centered at x(x, y)
    ell = centered(zeros(Float64, (301, 301)))
    for i ∈ 1:size(X, 1)
        ell[X[i, 2], -X[i, 1]] = 1.0
    end
    ell = parent(ell)

    # if gaussian
    gauss3x3 = centered(1/16 .* [1 2 1; 2 4 2; 1 2 1])
    ell = imfilter(ell, gauss3x3)
    # end

    slice1 = mapslices(col -> any(col .!= 0), ell, dims = 1)
    slice2 = mapslices(col -> any(col .!= 0), ell, dims = 2)

    return centered(ell[:, vec(slice1)][vec(slice2), :])
end

# i'll probably never use this
function kernelize(ellipse_image::OffsetArray{Float64,2,Array{Float64,2}})
    ell = parent(ellipse_image)
    slice1 = mapslices(col -> any(col .!= 0), ell, dims = 1)
    slice2 = mapslices(col -> any(col .!= 0), ell, dims = 2)

    return centered(ell[:, vec(slice1)][vec(slice2), :])
end


text:
"""
    MaskMulti(a, [classes])

An `N`-dimensional multilabel mask with labels `classes`.

## Examples

{cell=MaskMulti}
```julia
using DataAugmentation

mask = MaskMulti(rand(1:3, 100, 100))
```
{cell=MaskMulti}
```julia
showitems(mask)
```
"""
struct MaskMulti{N, T<:Integer, U} <: AbstractArrayItem{N, T}
    data::AbstractArray{T, N}
    classes::AbstractVector{U}
    bounds::Bounds{N}
end

function MaskMulti(a::AbstractArray, classes = unique(a))
    bounds = Bounds(size(a))
    minimum(a) >= 1 || error("Class values must start at 1")
    return MaskMulti(a, classes, bounds)
end

MaskMulti(a::AbstractArray{<:Gray{T}}, args...) where T = MaskMulti(reinterpret(T, a), args...)
MaskMulti(a::AbstractArray{<:Normed{T}}, args...) where T = MaskMulti(reinterpret(T, a), args...)
MaskMulti(a::IndirectArray, classes = a.values, bounds = Bounds(size(a))) =
    MaskMulti(a.index, classes, bounds)

Base.show(io::IO, mask::MaskMulti{N, T}) where {N, T} =
    print(io, "MaskMulti{$N, $T}() with size $(size(itemdata(mask))) and $(length(mask.classes)) classes")

getbounds(mask::MaskMulti) = mask.bounds

function project(P, mask::MaskMulti, bounds::Bounds)
    a = itemdata(mask)
    etp = mask_extrapolation(a)
    res = warp(etp, inv(P), bounds.rs)
    return MaskMulti(
        res,
        mask.classes,
        bounds
    )
end

function project!(bufmask::MaskMulti, P, mask::MaskMulti, bounds)
    a = OffsetArray(parent(itemdata(bufmask)), bounds.rs)
    warp!(
        a,
        mask_extrapolation(itemdata(mask)),
        inv(P),
    )
    return MaskMulti(
        a,
        mask.classes,
        bounds
    )
end

function showitem!(img, mask::MaskMulti)
    colors = distinguishable_colors(length(mask.classes))
    maskimg = map(itemdata(mask)) do val
        colors[findfirst(==(val), mask.classes)]
    end
    showimage!(img, maskimg)
end

# ## Binary masks

"""
    MaskBinary(a)

An `N`-dimensional binary mask.

## Examples

{cell=MaskMulti}
```julia
using DataAugmentation

mask = MaskBinary(rand(Bool, 100, 100))
```
{cell=MaskMulti}
```julia
showitems(mask)
```
"""
struct MaskBinary{N} <: AbstractArrayItem{N, Bool}
    data::AbstractArray{Bool, N}
    bounds::Bounds{N}
end

function MaskBinary(a::AbstractArray{Bool, N}, bounds = Bounds(size(a))) where N
    return MaskBinary(a, bounds)
end

Base.show(io::IO, mask::MaskBinary{N}) where {N} =
    print(io, "MaskBinary{$N}() with size $(size(itemdata(mask)))")

getbounds(mask::MaskBinary) = mask.bounds

function project(P, mask::MaskBinary, bounds::Bounds)
    etp = mask_extrapolation(itemdata(mask))
    return MaskBinary(
        warp(etp, inv(P), bounds.rs),
        bounds,
    )
end

function project!(bufmask::MaskBinary, P, mask::MaskBinary, bounds)
    a = OffsetArray(parent(itemdata(bufmask)), bounds.rs)
    res = warp!(
        a,
        mask_extrapolation(itemdata(mask)),
        inv(P),
    )
    return MaskBinary(
        a,
        bounds
    )
end

function showitem!(img, mask::MaskBinary)
    showimage!(img, colorview(Gray, itemdata(mask)))
end
# ## Helpers

function mask_extrapolation(
        mask::AbstractArray{T};
        t = T,
        degree = Constant(),
        boundary = Flat()) where T
    itp = interpolate(t, t, mask, BSpline(degree))
    etp = extrapolate(itp, Flat())
    return etp
end


text:
using ACEtb.SlaterKoster, Test, LinearAlgebra, Calculus
import ACEtb.SlaterKoster.CodeGeneration

SK = ACEtb.SlaterKoster
CG = ACEtb.SlaterKoster.CodeGeneration

CG.sk_table(2)

@show CG.Gsym(1, 0, 1, 0, 0)
CG.Gsym(0,0,0,0,0)
CG.Gsym(0,1,0,-1,0)

@show CG.Gsym(1,1,1,1,0)

# function expr(l1, l2, m1, m2)
#    ex = :(0.0)
#    for M = 0:min(l1,l2)
#       ex1 = Meta.parse(CG.Gsym(l1, l2, m1, m2, M))
#       # Vsym = Meta.parse("V$(SK.LSym[l1])$(SK.LSym[l2])$(SK.MSym[M])")
#       Vsym = Meta.parse("V[$(SK.bondintegral_index(l1, l2, M))]")
#       ex = Calculus.simplify(:($(ex)  + $(ex1) * $(Vsym)))
#    end
#    I1 = SK.orbital_index(l1, m1)
#    I2 = SK.orbital_index(l2, m2)
#    println("mat[$I1, $I2]")
#    println("   ", ex)
# end
#
# for l1 = 0:1, l2 = l1:1 # s,p,d
#    for m1=-l1:l1, m2=-l2:l2
#       ex = :(0.0)
#       for M = 0:min(l1,l2)
#          ex1 = Meta.parse(CG.Gsym(l1, l2, m1, m2, M))
#          # Vsym = Meta.parse("V$(SK.LSym[l1])$(SK.LSym[l2])$(SK.MSym[M])")
#          Vsym = Meta.parse("V[$(SK.bondintegral_index(l1, l2, M))]")
#          ex = Calculus.simplify(:($(ex)  + $(ex1) * $(Vsym)))
#       end
#       I1 = SK.orbital_index(l1, m1)
#       I2 = SK.orbital_index(l2, m2)
#       println("mat[$I1, $I2]")
#       println("   ", ex)
#    end
# end

# joinpath(@__DIR__(), "..", "src")
# using PyCall
# pysys = pyimport("sys")
# push!(pysys."path", joinpath(@__DIR__(), "..", "src"))
# codegen = pyimport("codegen")
# @show codegen.Gsym(1,1,0,0,0)
# ex = Meta.parse(replace(codegen.Gsym(1,1,0,0,0), "**" => "^"))
# ex1 = :( $ex + $ex )
# f = eval(:(beta -> $ex1))
# f(0.0)

CG.Gsym(2, 0, 1, 0, 0)
# 2 * sin(th) * cos(th) * cos(phi)
# = 2 * z * x
CG.Gsym(0, 2, 0, 2, 0)
# sin(th)^2*cos(2phi)
# = sin(th)^2*[cos^2(phi)-sin^2(phi)]
# = x^2 - y^2

## matrix block assembly

V = rand(10) #.- 0.5
U = rand(3) .- 0.5
U = U / norm(U)
Hold = zeros(9,9)
Hnew = zeros(9,9)
SK.OldSK.sk9!(U, V, Hold)
SK._sk!(Hnew, Val(2), U, V)
# perm = [1,2,4,3]
# @test Hold[perm, perm] ≈ Hnew

CG.sksign(0, 2)

perm = [1,3,4,2,5,6,9,7,8]
@test Hold[perm,perm] ≈ Hnew

set1 = JSON.parsefile(@__DIR__() * "/data/sp_o3_offsite_data.json")
@show set1["1"]["V"]
display(hcat(set1["1"]["E"]...))

set2 = JSON.parsefile(@__DIR__() * "/data/spdf_au2_offsite_data.json")
@show set2["1"]["V"]


text:
using Random
using Wavelets
using ContinuousWavelets
Random.seed!(135)
using Test
@testset "st using cwt tests" begin
    sharpExample = zeros(100, 2)
    sharpExample[26:75, :] .= 1
    layers = stParallel(2, size(sharpExample, 1))
    resSt = st(sharpExample, layers, abs)
    @test eltype(resSt) <: eltype(sharpExample) # check the type
    otherResSt = layers(sharpExample)
    compare = roll(resSt, layers, sharpExample)
    @test compare == otherResSt
    @test compare[2][:, end, end, 1] == resSt[end-12:end, 1]
    multiDExample = 1000 * randn(50, 10, 43)
    layers = stParallel(2, size(multiDExample, 1))
    resSt = st(multiDExample, layers, abs)
    @test eltype(resSt) <: eltype(multiDExample)
    @test 0 < minimum(abs.(resSt[:, 1, 1]))  # input is always non-zero so nothing should be *Exactly* zero, except by chance (which is why we set the seed above)
    @test 0 < minimum(abs.(resSt[:, end, end]))
    resSt = st(multiDExample, layers, abs, outputSubsample = (3, -1))
    @test eltype(resSt) <: eltype(multiDExample)
    @test 0 < minimum(abs.(resSt[:, 1, 1]))  # input is always non-zero so nothing should be *Exactly* zero, except by chance (which is why we set the seed above)
    @test 0 < minimum(abs.(resSt[:, end, end]))
end

f = randn(102)
lay = stParallel(1, size(f)[1])
ScatteredFull(lay, f)
# some basic tests for the Scattered type
# m=2, 3extra, 2transformed
fixDims = [3, 5, 2]
m = 2;
k = 2;
n = [100 200; 50 100; 25 50];
q = [5, 5, 5];
ex1 = ScatteredFull(m, k, fixDims, n, q, Float64)
@test ex1.m == m
@test ex1.k == k
@test size(ex1.data, 1) == m + 1
@test size(ex1.output, 1) == m + 1
@testset "correct size for output" begin
    for i = 1:size(ex1.output, 1)
        @test size(ex1.output[i])[4:end] == tuple(fixDims...)
        @test size(ex1.data[i])[4:end] == tuple(fixDims...)
        @test size(ex1.output[i])[1:2] == tuple(n[i, :]...)
        @test size(ex1.data[i])[1:2] == tuple(n[i, :]...)
    end
end


text:
module Beep

using WAV

include("sounds.jl")
include("functions.jl")
include("macros.jl")

export beep_sounds
export beep
export @beep
export load_sound

end


text:
const CENSUSPATH = joinpath(dirname(@__FILE__), "..", "data", "Census_data.dat")
const GENCOORDPATH = joinpath(dirname(@__FILE__), "..", "data", "Generator_coord.dat")
const GENDATAPATH = joinpath(dirname(@__FILE__), "..", "data", "Generator_data.dat")
const GENJSONPATH = joinpath(dirname(@__FILE__), "..", "data", "GenData.json")

function zipcode_builder(datapath = CENSUSPATH)
    df = CSV.read(datapath, delim='\t')
    zipcodes = Vector(size(df, 1))
    for i in 1:size(df, 1)
        zip = [df[i, 1], df[i, 2], df[i, 8], df[i, 9]]
        zipcodes[i] = zip
    end
    return zipcodes
end

function linear_load(bus::LoadBus)
    return 0.002 * bus.population # Overbye's first approxmation
end

function rand_volt(bus::Bus)
    volts = [100, 200, 350, 500]
    return volts[rand(1:length(volts))]
end

"""
    place_loads_from_zips!(
        grid::Grid;
        latlim = (-Inf,Inf),
        longlim = (-Inf,Inf),
        datapath = CENSUSPATH,
        loadfunc = linearload,
        voltfunc = randvolt
    )

Create load buses based on zipcode information.

# Arguments
* grid: The synthetic grid instance that will be altered.
* latlim: Tuple of the form (minimum latitude, maximum latitude).
* longlim: Tuple of the form (minimum longitude, maximum longitude).
* datapath: Path to the text file containing census information.
* loadfunc: Function for computing bus load based on some rule.
* voltfunc: Function for computing bus voltage based on some rule.

REFERENCE: https://www.census.gov/geo/maps-data/data/gazetteer2010.html
"""
function place_loads_from_zips!(
    grid::Grid;
    latlim = (-Inf, Inf),
    longlim = (-Inf, Inf),
    datapath = CENSUSPATH,
    loadfunc = linear_load,
    voltfunc = rand_volt,
)
    allzips = zipcode_builder(datapath)
    fzips = [
        z for z in allzips if z[3]>latlim[1] &&
        z[3]<latlim[2] &&
        z[4]>longlim[1] &&
        z[4]<longlim[2]
    ]

    srand(grid.seed + 23) # Just so we don't keep returning to the same point.
    for zip in fzips
        dummy = LoadBus(-1, (zip[3],zip[4]), -1, -1, zip[2])
        push!(
            buses(grid),
            LoadBus(
                length(buses(grid)) + 1,
                (zip[3],zip[4]),
                loadfunc(dummy), # This is pretty convoluted for rigid values, but makes it
                voltfunc(dummy), # extremely easy to extend in the future for more
                zip[2]           # ellaborate functions.
            )
        )
    end
end

"""
    place_loads_from_zips!(
        grid::Grid,
        geolim::Function;
        datapath = CENSUSPATH,
        loadfunc = linearload,
        voltfunc = randvolt
    )

Create load buses based on zipcode information.

# Arguments
* grid: The synthetic grid instance that will be altered.
* geolim: Function that receives a latitude-longitude pair and returns `true` or `false`.
Specifies the region of interest.
* datapath: Path to the text file containing census information.
* loadfunc: Function for computing bus load based on some rule.
* voltfunc: Function for computing bus voltage based on some rule.

REFERENCE: https://www.census.gov/geo/maps-data/data/gazetteer2010.html
"""
function place_loads_from_zips!(
    grid::Grid,
    geolim::Function;
    datapath = CENSUSPATH,
    loadfunc = linear_load,
    voltfunc = rand_volt,
)
    allzips = zipcode_builder(datapath)
    fzips = [z for z in allzips if geolim((z[3], z[4]))]

    srand(grid.seed + 23) # Just so we don't keep returning to the same point.
    for zip in fzips
        dummy = LoadBus(-1, (zip[3],zip[4]), -1, -1, zip[2])
        push!(
            buses(grid),
            LoadBus(
                length(buses(grid)) + 1,
                (zip[3],zip[4]),
                loadfunc(dummy), # This is pretty convoluted for rigid values, but makes it
                voltfunc(dummy), # extremely easy to extend in the future for more
                zip[2]           # ellaborate functions.
            )
        )
    end
end

function get_plant_data(coordpath = GENCOORDPATH)
    df = CSV.read(coordpath, rows_for_type_detect=typemax(Int))
    pcodes = sizehint!(Int[], size(df, 1))
    pcoords = sizehint!([], size(df, 1))
    pvolts = sizehint!([], size(df, 1))
    for i in 1:size(df, 1)
        push!(pcodes, df[i, Symbol("Plant Code")])
        push!(
            pcoords,
            (df[i, Symbol("Latitude")], df[i, Symbol("Longitude")])
        )
        vs = Real[]
        for s in [
            Symbol("Grid Voltage (kV)"),
            Symbol("Grid Voltage 2 (kV)"),
            Symbol("Grid Voltage 3 (kV)")
        ]
            if !ismissing(df[i, s])
                push!(vs, df[i, s])
            end
        end
        push!(pvolts, vs)
    end
    return pcodes, pcoords, pvolts
end

function get_gen_data(
    datapath = GENDATAPATH,
    coordpath = GENCOORDPATH,
)
    # Generator data is separated in two spreadsheets.
    # We have to parse both and then combine.

    #Grabbing plant information
    pcodes, pcoords, pvolts = get_plant_data(coordpath)
    #Grabbing generator information
    keys2keep = [
        "Plant Code",
        "Technology",
        "Nameplate Capacity (MW)",
        "Nameplate Power Factor",
        "Minimum Load (MW)",
        "Summer Capacity (MW)",
        "Winter Capacity (MW)",
        "Time from Cold Shutdown to Full Load",
        "Status"
    ]
    df = CSV.read(datapath, rows_for_type_detect=typemax(Int))
    tempdata = sizehint!([], size(df, 1))
    for i in 1:size(df, 1)
        gdata = []
        for k in keys2keep[2:end]
            push!(gdata, df[i, Symbol(k)])
        end
        for j in 1:length(pcodes)
            if pcodes[j] == df[i, Symbol(keys2keep[1])]
                gdata = vcat([pcoords[j], pvolts[j]], gdata)
                break
            end
        end
        push!(tempdata, gdata)
    end
    plants = []
    for c in pcoords
        p = [g for g in tempdata if g[1] == c]
        length(p) > 0 ? push!(plants, p) : nothing #lots of plants have no
        # corresponding generators. This means our survey data is not so good.
    end
    return plants
end

function format_gen!(gen::Dict)
    gen["cap"] = gen["cap"]
    gen["pfactor"] = !ismissing(gen["pfactor"]) ? gen["pfactor"] : 1
    gen["minload"] = !ismissing(gen["minload"]) ? gen["minload"] : 0
    gen["scap"] = !ismissing(gen["scap"]) ? gen["scap"] : gen["cap"]
    gen["wcap"] = !ismissing(gen["wcap"]) ? gen["wcap"] : gen["cap"]
    gen["shut2loadtime"] = !ismissing(gen["shut2loadtime"]) ? gen["shut2loadtime"] : " "
    gen["status"] = !ismissing(gen["status"]) ? gen["status"] : " "
end

"""
    prepare_gen_data(indata = "Generator_data.dat",
                    incoord = "Generator_coord.dat",
                    outfile = "GenData.json")

Prepare json (outfile) with generator data gathered from the survey spreadsheets
(indata and incoord). Output is organized as an array of plants. Each plant is
an array of generators. Each generator is a dictionary.
"""
function prepare_gen_data(
    indata = GENDATAPATH,
    incoord = GENCOORDPATH,
    outfile = GENJSONPATH,
)
    plants = get_gen_data(indata, incoord)
    fplants = sizehint!([], length(plants))
    keys = ["coords", "volt", "tech", "cap", "pfactor", "minload", "scap",
            "wcap", "shut2loadtime", "status"]
    for p in plants
        push!(fplants, [Dict(k => v for (k,v) in zip(keys,g)) for g in p])
        for g in fplants[end]
            format_gen!(g)
        end
    end

    open(outfile, "w") do f
        JSON.print(f, fplants)
    end
end

function format_gen!(gen::Generator)
    gen.cap = gen.cap
    gen.pfactor = !ismissing(gen.pfactor) ? gen.pfactor : 1
    gen.minload = !ismissing(gen.minload) ? gen.minload : 0
    gen.scap = !ismissing(gen.scap) ? gen.scap : gen.cap
    gen.wcap = !ismissing(gen.wcap) ? gen.wcap : gen.cap
    gen.shut2loadtime = !ismissing(gen.shut2loadtime) ? gen.shut2loadtime : " "
    gen.status = !ismissing(gen.status) ? gen.status : " "
end

function pfacwavg(plant::AbstractArray{Generator})
    on = filter(g -> g.status in ["OP", "SP"], plant)
    if length(on) > 0
        caps = [g.cap for g in on]
        pfacts = [g.pfactor for g in on]
        pfactor = sum(pfacts .* caps)
        return pfactor / sum(caps)
    else
        return 0
    end
end

"""
    place_gens_from_data!(
        grid::Grid,
        plants::Array,
        pfactorfunc::Function = pfacwavg
    )

Create generation buses based on power plant data.

# Arguments:
* grid: Grid instance which will be populated.
* plants: Power plant information already parsed.
* pfactorfunc: Function for computing bus power factor based on some rule.
"""
function place_gens_from_data!(
    grid::Grid,
    plants::Array,
    pfactorfunc::Function = pfacwavg,
)
    # Creating generators
    genbs = sizehint!(Vector{Vector{Generator}}(), length(plants))
    for p in plants
        push!(
            genbs,
            [
                Generator(
                    (g["coords"][1],g["coords"][2]),
                    g["volt"],
                    g["tech"],
                    g["cap"],
                    g["pfactor"],
                    g["minload"],
                    g["scap"],
                    g["wcap"],
                    g["shut2loadtime"],
                    g["status"]
                ) for g in p
            ]
        )
    end
    # Creating generation buses
    for p in genbs
        op = filter(g -> g.status in ["OP", "SP"], p)
        push!(buses(grid), GenBus(length(buses(grid)) + 1, p, op, pfactorfunc))
    end
end

"""
    place_gens_from_data!(
        grid::Grid;
        latlim::Tuple = (-Inf, Inf),
        longlim::Tuple = (-Inf, Inf),
        jsonpath = GENJSONPATH,
        pfactorfunc::Function = pfacwavg,
	)

Create generation buses based on power plant data.

# Arguments:
* grid: Grid instance which will be populated.
* latlim: Tuple of the form (minimum latitude, maximum latitude).
* longlim: Tuple of the form (minimum longitude, maximum longitude).
* jsonpath: Path to the JSON file with the generator data.
* pfactorfunc: Function for computing bus power factor based on some rule.

REFERENCE: https://www.eia.gov/electricity/data/eia860/index.html
"""
function place_gens_from_data!(
    grid::Grid;
    latlim::Tuple = (-Inf, Inf),
    longlim::Tuple = (-Inf, Inf),
    jsonpath = GENJSONPATH,
    pfactorfunc::Function = pfacwavg,
)
    #Loading data from json
    plants = JSON.parsefile(jsonpath)
    #Filtering per latitude and longitude
    plants = [
        p for p in plants if p[1]["coords"][1] > latlim[1] &&
        p[1]["coords"][1] < latlim[2] &&
        p[1]["coords"][2] > longlim[1] &&
        p[1]["coords"][2] < longlim[2]
    ]
    place_gens_from_data!(grid, plants, pfactorfunc)
end

"""
    place_gens_from_data!(
        grid::Grid,
        geolim::Function;
        jsonpath = GENJSONPATH,
        pfactorfunc::Function = pfacwavg,
	)

Create generation buses based on power plant data.

# Arguments:
* grid: Grid instance which will be populated.
* geolim: Function that receives a latitude-longitude pair and returns `true` or `false`.
Specifies the region of interest.
* jsonpath: Path to the JSON file with the generator data.
* pfactorfunc: Function for computing bus power factor based on some rule.

REFERENCE: https://www.eia.gov/electricity/data/eia860/index.html
"""
function place_gens_from_data!(
    grid::Grid,
    geolim::Function;
    jsonpath = GENJSONPATH,
    pfactorfunc::Function = pfacwavg,
)
    #Loading data from json
    plants = JSON.parsefile(jsonpath)
    #Filtering per latitude and longitude
    plants = [p for p in plants if geolim((p[1]["coords"][1], p[1]["coords"][2]))]
    place_gens_from_data!(grid, plants, pfactorfunc)
end

# """
#     place_gens_from_data_old!(grid::Grid, latlim=(-Inf,Inf),
#                       longlim=(-Inf,Inf), datapath="Generator_data.dat",
#                       coordpath = "Generator_coord.dat", pfactorfunc = pfacwavg)
#
# OBSOLETE. Create generation buses based on power plant data.
# Obsolete now that we have a standalone function for creating a json with the
# data. Reading in the json file via place_gens_from_data!() is considerably faster.
#
# # Arguments:
# * grid: Grid instance which will be populated.
# * latlim: Tuple of the form (minimum latitude, maximum latitude).
# * longlim: Tuple of the form (minimum longitude, maximum longitude).
# * datapath: Path to the text file containing detailed generator information.
# * coordpath: Path to the text file containing power plant information.
# * pfactorfunc: Function for computing bus power factor based on some rule.
#
# REFERENCE: https://www.eia.gov/electricity/data/eia860/index.html
# """
# function place_gens_from_data_old!(
#     grid::Grid;
#     latlim = (-Inf, Inf),
#     longlim = (-Inf, Inf),
#     datapath = GENDATAPATH,
#     coordpath = GENCOORDPATH,
#     pfactorfunc = pfacwavg,
# )
#     plants = get_gen_data(datapath, coordpath)
#     plants = [
#                 p for p in plants if p[1][1][1] > latlim[1] &&
#                 p[1][1][1] < latlim[2] && p[1][1][2] > longlim[1] &&
#                 p[1][1][2] < longlim[2]
#             ]
#     # Creating generators
#     genbs = Vector(plants)
#     for p in plants
#         push!(
#             genbs,
#             [
#                 Generator(
#                     g[1],
#                     g[2],
#                     g[3],
#                     g[4],
#                     g[5],
#                     g[6],
#                     g[7],
#                     g[8],
#                     g[9],
#                     g[10]
#                 ) for g in p
#             ]
#         )
#     end
#     # Formatting generators
#     for p in genbs
#         for g in p
#             format_gen!(g)
#         end
#     end
#     # Creating generation buses
#     for p in genbs
#         op = filter(g -> g.status in ["OP", "SP"], p)
#         push!(buses(grid), GenBus(length(buses(grid)) + 1, p, op, pfactorfunc))
#     end
# end

# Here we will adopt transmission lines of type 'Line' with predefined values. We need
# their predefined values because we do not have all properties required to create a
# 'Line'. The use of 'DC Line' is precluded by the fact that, in pandapower, it is only
# capable of unidirectional flow.
const SAFE_LOAD_PERCENT = 100 # Maximum load allowed at lines and transformers.
function add_line(pgrid::PyObject, b1::LoadBus, b2::LoadBus)
    pp[:create_line](
        pgrid,
        from_bus = (b1.id - 1), # Adopting zero indexing because this is for Python.
        to_bus = (b2.id - 1), # Adopting zero indexing because this is for Python.
        length_km = haversine(b1, b2),
        std_type = "305-AL1/39-ST1A 110.0",
        max_loading_percent = SAFE_LOAD_PERCENT,
    )
end
function add_line(pgrid::PyObject, b1::LoadBus, b2::GenBus)
    pp[:create_transformer](
        pgrid,
        hv_bus = (b2.id - 1), # Adopting zero indexing because this is for Python.
        lv_bus = (b1.id - 1), # Adopting zero indexing because this is for Python.
        std_type = "160 MVA 380/110 kV",
        max_loading_percent = SAFE_LOAD_PERCENT,
    )
end
function add_line(pgrid::PyObject, b1::GenBus, b2::LoadBus)
    pp[:create_transformer](
        pgrid,
        hv_bus = (b1.id - 1), # Adopting zero indexing because this is for Python.
        lv_bus = (b2.id - 1), # Adopting zero indexing because this is for Python.
        std_type = "160 MVA 380/110 kV",
        max_loading_percent = SAFE_LOAD_PERCENT,
    )
end
function add_line(pgrid::PyObject, b1::GenBus, b2::GenBus)
    pp[:create_line](
        pgrid,
        from_bus = (b1.id - 1), # Adopting zero indexing because this is for Python.
        to_bus = (b2.id - 1), # Adopting zero indexing because this is for Python.
        length_km = haversine(b1, b2),
        std_type = "490-AL1/64-ST1A 380.0",
        max_loading_percent = SAFE_LOAD_PERCENT,
    )
end

"""
    to_pandapower(grid::Grid)

Export a 'grid' to pandapower format.

Important:
Currently, grid voltages and line properties are ignored. This function places all loads at
100kV and all generation at 380kV. Transmission lines use standard values for 110kV and
380kV lines. Transformers are automatically placed between buses that operate at different
voltage values.
"""
function to_pandapower(grid::Grid,)
    const LOAD_VOLT = 110 # kV # This is a crude approximation
    const GEN_VOLT = 380 # kV # This is a crude approximation
    voltage(bus::LoadBus) = LOAD_VOLT
    voltage(bus::GenBus) = GEN_VOLT
    # Create empty pandapower network
    pgrid = pycall(pp[:create_empty_network], PyObject)
    # Create buses
    for i in 1:length(buses(grid))
        pp[:create_bus](
            pgrid,
            voltage(buses(grid)[i]),
            index = (i - 1), # Adopting zero indexing because this is for Python.
            geodata = (buses(grid)[i].coords.lat, buses(grid)[i].coords.lon)
        ) # Despite what the documentation says, 'create_bus' requires voltage.
    end
    # Create loads and generators
    for i in 1:length(buses(grid))
        if isa(buses(grid)[i], LoadBus)
            pp[:create_load](
                pgrid,
                bus = (i-1), # Adopting zero indexing because this is for Python.
                p_kw = 1000*buses(grid)[i].load, # MW to kW
                controllable = false # for OPF
            )
        else
        # The type of the generator is being stored in 'name' instead
        # of in 'type' because Julia seems to have problems with a parameter called 'type'.
            for gen in buses(grid)[i].gens
                pp[:create_gen](
                    pgrid,
                    bus = (i-1), # Adopting zero indexing because this is for Python.
                    p_kw = -1000 * gen.cap, # negative by pandapower convention
                    name = gen.tech,
                    max_p_kw = -1000 * gen.minload, # min inverts with max because negative
                    min_p_kw = -1000 * gen.cap,
                    controllable = true, # for OPF
                )
            end
        end
    end
    # Create lines and transformers
    line_types = ( # All pre-defined overhead line types. Seem to be for lower voltages.
        "149-AL1/24-ST1A 10.0",
        "149-AL1/24-ST1A 110.0",
        "149-AL1/24-ST1A 20.0",
        "15-AL1/3-ST1A 0.4",
        "184-AL1/30-ST1A 110.0",
        "184-AL1/30-ST1A 20.0",
        "24-AL1/4-ST1A 0.4",
        "243-AL1/39-ST1A 110.0",
        "243-AL1/39-ST1A 20.0",
        "305-AL1/39-ST1A 110.0",
        "48-AL1/8-ST1A 0.4",
        "48-AL1/8-ST1A 10.0",
        "48-AL1/8-ST1A 20.0",
        "490-AL1/64-ST1A 220.0",
        "490-AL1/64-ST1A 380.0",
        "94-AL1/15-ST1A 0.4",
        "94-AL1/15-ST1A 10.0",
        "94-AL1/15-ST1A 20.0",
    )
    lines_110 = (
        "149-AL1/24-ST1A 110.0",
        "184-AL1/30-ST1A 110.0",
        "243-AL1/39-ST1A 110.0",
        "305-AL1/39-ST1A 110.0",
    )
    for ln in grid.trans_lines
        buses = unique(ln.connecting)
        add_line(pgrid, buses[1], buses[2])
    end
    pp[:create_ext_grid](pgrid, 0, vm_pu=1, max_p_kw=0, min_p_kw=0) # Just so pandapower does
                                                    # not throw an error when doing and OPF.
    return pgrid
end

"""
    to_pandapower(grid::Grid, filename::AbstractString)D

Export a 'grid' to pandapower format and saves it as 'filename'. Pandapower requires
'filename' to be a '.p' file.

Important:
Currently, grid voltages and line properties are ignored. This function places all loads at
100kV and all generation at 380kV. Transmission lines use standard values for 110kV and
380kV lines. Transformers are automatically placed between buses that operate at different
voltage values.
"""
function to_pandapower(grid::Grid, filename::AbstractString)
    pgrid = to_pandapower(grid)
    save_pp_grid(pgrid, filename)
    return pgrid
end

"""
    save_pp_grid(pgrid::PyObject, filename::AbstractString)

Save a pandapower grid as `filename`. Pandapower requires `filename` to be a `.p` file.
"""
function save_pp_grid(pgrid::PyObject, path::AbstractString)
    pp[:to_pickle](pgrid, path)
end

"""
    load_pp_grid(path)

Load a pandapower grid from `path`.
"""
function load_pp_grid(path)
    return pycall(pp[:from_pickle], PyObject,path)
end

# WARNING: exportJLD and loadJLD functions are not working properly yet.
# They should be disregarded.
function save(grid::Grid, outname)
    # Breaking object links in order to avoid StackOverflowError
    dummytuple = (LoadBus(-1, (-1,-1), -1, -1, -1), LoadBus(-1, (-1,-1), -1, -1, -1))
    for ln in grid.trans_lines
        ln.connecting = dummytuple
    end
    for bus in buses(grid)
        bus.connected_to = []
        bus.connections = []
    end
    for sub in substations(grid)
        sub.connected_to = []
        # sub.grouping = []
    end
    #Storing object
    jldopen(outname, "w") do f
        # addrequire(f, SyntheticGrids)
        write(f, "grid", grid)
    end
end

function load_grid(inname)
    # Load grid
    grid = load(inname)
    # Reconnect grid
    connect_buses!(grid)
    count = 1
    for j in 1:(length(buses(grid)) - 1), i in (j + 1):length(buses(grid))
        if grid.bus_conn[i, j] == 1
            grid.trans_lines[count].connecting = (buses(grid)[i], buses(grid)[j])
            push!(buses(grid)[i].connections, grid.trans_lines[count])
            push!(buses(grid)[j].connections, grid.trans_lines[count])
            count += 1
        end
    end
    if length(substations(grid)) > 0
        connect_subs!(grid)
    end
    return grid
end

"""
    save(grid::Grid, outfile::AbstractString)

Save a `grid` in a format that can be recovered later. The format adopted is a JSON file
with a particular structure.
"""
function save(grid::Grid, outfile::AbstractString)
    out_dict = Dict()
    out_dict["seed"] = grid.seed
    bus_vec = []
    for b in buses(grid)
        if isa(b, LoadBus)
            push!(bus_vec, Dict(
                "type" => "LoadBus",
                "coords" => b.coords,
                "load" => b.load,
                "voltage" => b.voltage,
                "population" => b.population
            ))
        elseif isa(b, GenBus)
            push!(bus_vec, Dict(
                "type" => "GenBus",
                "coords" => b.coords,
                "generation" => b.generation,
                "voltage" => b.voltage,
                "tech_type" => b.tech_type,
                "pfactor" => b.pfactor,
                "summgen" => b.summgen,
                "wintgen" => b.wintgen,
                "gens" => [Dict(
                    "coords" => g.coords,
                    "volt" => g.volt,
                    "tech" => g.tech,
                    "cap" => g.cap,
                    "pfactor" => g.pfactor,
                    "minload" => g.minload,
                    "scap" => g.scap,
                    "wcap" => g.wcap,
                    "shut2loadtime" => g.shut2loadtime,
                    "status" => g.status
                ) for g in b.gens]
            ))
        else
            warn("Bus $b is not of a known type. Its type is $(typeof(b))")
        end
    end
    out_dict["buses"] = bus_vec
    out_dict["trans_lines"] = [Dict(
        "impedance" => t.impedance,
        "capacity" => t.capacity,
        "connecting" => [t.connecting[1].id, t.connecting[2].id]
    ) for t in trans_lines(grid)]
    out_dict["substations"] = [Dict(
        "coords" => s.coords,
        "voltages" => s.voltages,
        "load" => s.load,
        "generation" => s.generation,
        "population" => s.population,
        "grouping" => [b.id for b in s.grouping]
    ) for s in substations(grid)]
    out_dict["bus_conn"] = adjacency(grid)

    open(outfile, "w") do f
        JSON.print(f, out_dict)
    end
end

"""
    load_grid(filepath::AbstractString)

Load a synthetic grid from a file previously saved via `save(grid::Grid)`. Due to
conversions during the dumping and the reading of the file, floats may have some noise
added. That should be within machine precision, but may affect operations such as `==`.
"""
function load_grid(filepath::AbstractString)
    in_grid = JSON.parsefile(filepath)
    out_grid = Grid(in_grid["seed"])
    out_grid.bus_conn = hcat(in_grid["bus_conn"]...)
    for id in 1:length(in_grid["buses"])
        if in_grid["buses"][id]["type"] == "LoadBus"
            push!(
                out_grid.buses,
                LoadBus(
                    id,
                    LatLon(
                        in_grid["buses"][id]["coords"]["lat"],
                        in_grid["buses"][id]["coords"]["lon"]
                    ),
                    in_grid["buses"][id]["load"],
                    in_grid["buses"][id]["voltage"],
                    in_grid["buses"][id]["population"],
                )
            )
        elseif in_grid["buses"][id]["type"] == "GenBus"
            push!(
                out_grid.buses,
                GenBus(
                    id,
                    LatLon(
                        in_grid["buses"][id]["coords"]["lat"],
                        in_grid["buses"][id]["coords"]["lon"]
                    ),
                    in_grid["buses"][id]["generation"],
                    in_grid["buses"][id]["voltage"],
                    in_grid["buses"][id]["tech_type"],
                    Set(),
                    Set(),
                    in_grid["buses"][id]["pfactor"],
                    in_grid["buses"][id]["summgen"],
                    in_grid["buses"][id]["wintgen"],
                    [Generator(
                        LatLon(g["coords"]["lat"], g["coords"]["lon"]),
                        g["volt"],
                        g["tech"],
                        g["cap"],
                        g["pfactor"],
                        g["minload"],
                        g["scap"],
                        g["wcap"],
                        g["shut2loadtime"],
                        g["status"]
                    ) for g in in_grid["buses"][id]["gens"]]
                )
            )
        end
    end
    connect_buses!(out_grid)
    for id in 1:length(in_grid["substations"])
        push!(
            out_grid.substations,
            Substation(
                id,
                LatLon(
                    in_grid["substations"][id]["coords"]["lat"],
                    in_grid["substations"][id]["coords"]["lon"]
                ),
                in_grid["substations"][id]["voltages"],
                in_grid["substations"][id]["load"],
                in_grid["substations"][id]["generation"],
                in_grid["substations"][id]["population"],
                Set(),
                [out_grid.buses[i] for i in in_grid["substations"][id]["grouping"]]
            )
        )
    end
    connect_subs!(out_grid)
    for t in in_grid["trans_lines"]
        push!(out_grid.trans_lines, TransLine(
            (out_grid.buses[t["connecting"][1]], out_grid.buses[t["connecting"][2]]),
            t["impedance"],
            t["capacity"]
        ))
        push!(out_grid.buses[t["connecting"][1]].connections, out_grid.trans_lines[end])
        push!(out_grid.buses[t["connecting"][2]].connections, out_grid.trans_lines[end])
    end
    return out_grid
end


text:
import AoC2019
using Test
task_commands = AoC2019.read_lines("input22.txt")
test_commands1 = """deal with increment 7
deal into new stack
deal into new stack
Result: 0 3 6 9 2 5 8 1 4 7"""
test_commands2 = """deal into new stack
cut -2
deal with increment 7
cut 8
cut -4
deal with increment 7
cut 3
deal with increment 9
deal with increment 3
cut -1"""

function do_commands_helper(commands; num_cards, track_card, do_reverse=false, iterations=1)::Array{Int128,1}
  current_num::Int128 = track_card
  operations = Int128[1,0]
  if do_reverse
    commands = reverse(commands)
  end
  for command in commands
    println(command)
    tokens = split(command)
    if tokens[1] == "deal"
      if tokens[2] == "into"
        current_num = (num_cards - 1) - current_num
        #current_num = mod(-current_num, num_cards)
        operations[2] = (num_cards -1) - operations[2]
        operations[1] *= -1
      elseif tokens[2] == "with"
        local N = parse(Int, tokens[4])
        if ! do_reverse
          current_num = mod(current_num * N, num_cards)
          operations .*= N
        else
          # find inverse to N mod num_cards
          N_inv = invmod(N, num_cards)
          current_num = mod(current_num * N_inv, num_cards)
          operations .*= N_inv
        end
      end
    elseif tokens[1] == "cut"
      local N = parse(Int, tokens[2])
      if ! do_reverse
        current_num = mod(current_num - N, num_cards)
        operations[2] -= N
      else
        current_num = mod(current_num + N, num_cards)
        operations[2] += N
      end
    end
    operations = mod.(operations, num_cards)
    println("operations: $operations")
    println("Current_num: $current_num")
    current_num_alt = mod(track_card * operations[1] + operations[2], num_cards)
    println("current_num_alt: $current_num_alt")
  end
  #println(operations)
  return operations
end

function fast_exp(kernel_operation::Array{Int128, 1}; iterations::Int64, modulo::Int64)::Array{Int128, 1}
  binary = bitstring(iterations)
  combined_operation::Array{Int128, 1} = kernel_operation
  new_operation::Array{Int128, 1} = Int128[1,0]
  for entry in reverse(digits(iterations, base=2))[2:end]
    new_operation = copy(combined_operation)
    new_operation .*= combined_operation[1]
    new_operation[2] += combined_operation[2]
    new_operation = mod.(new_operation, modulo)
    if entry == 1
      new_operation .*= kernel_operation[1]
      new_operation[2] += kernel_operation[2]
    end
    new_operation = mod.(new_operation, modulo)
    combined_operation = new_operation
  end
  return map(x -> Int128(x), combined_operation)
end

@test fast_exp(Int128[2,1], iterations=Int64(2), modulo=Int64(50)) == [4,3]
@test fast_exp(Int128[2,1], iterations=Int64(3), modulo=Int64(50)) == [8,7]
@test fast_exp(Int128[2,1], iterations=Int64(4), modulo=Int64(50)) == [16,15]
@test fast_exp(Int128[2,2], iterations=Int64(2), modulo=Int64(50)) == [4, 6]
@test fast_exp(Int128[2,2], iterations=Int64(4), modulo=Int64(50)) == [16, 30]
@test fast_exp(Int128[2,2], iterations=Int64(3), modulo=Int64(50)) == [8,14]

function do_commands(commands; num_cards, track_card, do_reverse=false, iterations=1)::Int128
  operations::Array{Int128, 1} = do_commands_helper(commands, num_cards=num_cards, track_card=track_card, do_reverse=do_reverse, iterations=iterations)
  operations = fast_exp(operations, iterations=iterations, modulo=num_cards)
  println(operations)
  return mod(track_card * operations[1] + operations[2], num_cards)
end

task_result = do_commands(task_commands, num_cards=10007, track_card=2019, do_reverse=false)
@test task_result == 7860

@test do_commands(split(test_commands1, '\n'), num_cards=10, track_card=7) == 9
@test do_commands(split(test_commands1, '\n'), num_cards=10, track_card=9, do_reverse=true) == 7

test_result1 = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=2019, do_reverse=false)
test_result2 = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=test_result1, do_reverse=false)
test_result2_alt = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=2019, do_reverse=false, iterations=2)
@test test_result2 == test_result2_alt

test_result3 = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=test_result2, do_reverse=false)
test_result3_alt = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=2019, do_reverse=false, iterations=3)
@test test_result3 == test_result3_alt

test_result1_back = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=test_result2, do_reverse=true)
test_result0_back = do_commands(split(test_commands2, '\n'), num_cards=10007, track_card=test_result1_back, do_reverse=true)
@test test_result0_back == 2019

iterations= 101741582076661
num_cards = 119315717514047

task_operations_1 = do_commands(task_commands, num_cards=119315717514047, track_card=2020, do_reverse=true, iterations=1)
task_operations_2 = do_commands(task_commands, num_cards=119315717514047, track_card=task_operations_1, do_reverse=true, iterations=1)
task_operations_2_alt = do_commands(task_commands, num_cards=119315717514047, track_card=2020, do_reverse=true, iterations=2)
@test task_operations_2 == task_operations_2_alt

task_operations_3 = do_commands(task_commands, num_cards=119315717514047, track_card=task_operations_2, do_reverse=true, iterations=1)
task_operations_3_alt = do_commands(task_commands, num_cards=119315717514047, track_card=2020, do_reverse=true, iterations=3)
@test task_operations_3 == task_operations_3_alt
task_operations_4 = do_commands(task_commands, num_cards=119315717514047, track_card=task_operations_3, do_reverse=true, iterations=1)
task_operations_4_alt = do_commands(task_commands, num_cards=119315717514047, track_card=2020, do_reverse=true, iterations=4)
@test task_operations_4 == task_operations_4_alt
task_operations_7 = do_commands(task_commands, num_cards=119315717514047, track_card=task_operations_4, do_reverse=true, iterations=3)
task_operations_7_alt = do_commands(task_commands, num_cards=119315717514047, track_card=2020, do_reverse=true, iterations=7)
@test task_operations_7 == task_operations_7_alt
task_operations_0_back = do_commands(task_commands, num_cards=119315717514047, track_card=task_operations_7, do_reverse=false, iterations=7)
@test task_operations_0_back == 2020

task_operations = do_commands(task_commands, num_cards=119315717514047, track_card=2020, do_reverse=true, iterations=iterations)
#fast_exp(task_operations, UInt64(iterations), UInt64(num_cards))
println(UInt64(task_operations))

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc