text:
include("TwoPartCost.jl")
include("ThreePartCost.jl")
include("TechHydro.jl")
include("TechRenewable.jl")
include("TechThermal.jl")
include("Area.jl")
include("LoadZone.jl")
include("Bus.jl")
include("Arc.jl")
include("Line.jl")
include("MonitoredLine.jl")
include("PhaseShiftingTransformer.jl")
include("TapTransformer.jl")
include("Transformer2W.jl")
include("HVDCLine.jl")
include("VSCDCLine.jl")
include("InterruptibleLoad.jl")
include("FixedAdmittance.jl")
include("PowerLoad.jl")
include("HydroEnergyReservoir.jl")
include("HydroDispatch.jl")
include("RenewableDispatch.jl")
include("RenewableFix.jl")
include("ThermalStandard.jl")
include("GenericBattery.jl")
include("StaticReserve.jl")
include("VariableReserve.jl")
include("Transfer.jl")
include("AVRFixed.jl")
include("AVRSimple.jl")
include("AVRTypeI.jl")
include("AVRTypeII.jl")
include("BaseMachine.jl")
include("AndersonFouadMachine.jl")
include("FullMachine.jl")
include("MarconatoMachine.jl")
include("OneDOneQMachine.jl")
include("SimpleAFMachine.jl")
include("SimpleFullMachine.jl")
include("SimpleMarconatoMachine.jl")
include("PSSFixed.jl")
include("PSSSimple.jl")
include("SingleMass.jl")
include("FiveMassShaft.jl")
include("TGFixed.jl")
include("TGTypeI.jl")
include("TGTypeII.jl")
include("AvgCnvFixedDC.jl")
include("FixedDCSource.jl")
include("LCLFilter.jl")
include("LCFilter.jl")
include("PLL.jl")
include("VirtualInertia.jl")
include("ReactivePowerDroop.jl")
include("CombinedVIwithVZ.jl")
include("Source.jl")

export get_Ae
export get_Be
export get_D
export get_D_12
export get_D_23
export get_D_34
export get_D_45
export get_D_ex
export get_D_hp
export get_D_ip
export get_D_lp
export get_Emf
export get_H
export get_H_ex
export get_H_hp
export get_H_ip
export get_H_lp
export get_K0
export get_K_ex
export get_K_hp
export get_K_ip
export get_K_lp
export get_K_p
export get_K_ω
export get_Ka
export get_Ke
export get_Kf
export get_Kv
export get_L_1d
export get_L_1q
export get_L_ad
export get_L_aq
export get_L_d
export get_L_f1d
export get_L_ff
export get_L_q
export get_MVABase
export get_P_max
export get_P_min
export get_R
export get_R_1d
export get_R_1q
export get_R_f
export get_T1
export get_T2
export get_T3
export get_T4
export get_T5
export get_T_AA
export get_Ta
export get_Tc
export get_Td0_p
export get_Td0_pp
export get_Te
export get_Tf
export get_Tq0_p
export get_Tq0_pp
export get_Tr
export get_Ts
export get_V_I
export get_V_R
export get_V_pss
export get_Vr_max
export get_Vr_min
export get_X_th
export get_Xd
export get_Xd_p
export get_Xd_pp
export get_Xq
export get_Xq_p
export get_Xq_pp
export get_Y
export get_activepower
export get_activepower_flow
export get_activepowerlimits
export get_activepowerlimits_from
export get_activepowerlimits_to
export get_angle
export get_anglelimits
export get_arc
export get_area
export get_available
export get_b
export get_basevoltage
export get_bus
export get_bustype
export get_capacity
export get_cf
export get_efficiency
export get_energy
export get_eq_p
export get_ext
export get_fixed
export get_flowlimits
export get_from
export get_fuel
export get_inflow
export get_initial_storage
export get_inputactivepowerlimits
export get_internal
export get_inv_d_fluxlink
export get_inv_q_fluxlink
export get_inverter_firing_angle
export get_inverter_taplimits
export get_inverter_xrc
export get_kad
export get_kd
export get_kffi
export get_kffv
export get_ki_pll
export get_kic
export get_kiv
export get_kp_pll
export get_kpc
export get_kpv
export get_kq
export get_kω
export get_lf
export get_lg
export get_load_zone
export get_loss
export get_lv
export get_maxactivepower
export get_maxreactivepower
export get_model
export get_n_states
export get_number
export get_op_cost
export get_outputactivepowerlimits
export get_powerfactor
export get_primaryshunt
export get_primemover
export get_r
export get_ramplimits
export get_rate
export get_rating
export get_reactivepower
export get_reactivepower_flow
export get_reactivepowerlimits
export get_reactivepowerlimits_from
export get_reactivepowerlimits_to
export get_rectifier_firing_angle
export get_rectifier_taplimits
export get_rectifier_xrc
export get_requirement
export get_rf
export get_rg
export get_rv
export get_s_rated
export get_services
export get_shutdn
export get_startup
export get_states
export get_storage_capacity
export get_tap
export get_tech
export get_timeframe
export get_timelimits
export get_to
export get_v_rated
export get_variable
export get_voltage
export get_voltagelimits
export get_x
export get_α
export get_γd
export get_γq
export get_τ_max
export get_τ_min
export get_ω_lp
export get_ωad
export get_ωb
export get_ωf


text:
"""
This is a sample input file for the CloudScat code.  To use it run julia
and type

julia> include("first_hitting.jl"); run()

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia first_hitting.jl

- - -
2010 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat
using FastGaussQuadrature

const co = CloudScat.constants

function run()
    params = Params(
        # NUmber of simulated photons
        N = 5000000,

        # The wavelength of the scattered light
        λ = 337 * co.nano,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = [0, 0, 10 * co.kilo],
        source_b = [0, 0, 10 * co.kilo],

        # For this example we neglect Rayleigh scattering
        σray = 0.)

    cth = 15 * co.kilo

    # World geometry
    # Describe the cloud geometry here.  It may consist in combinations of
    # geometrical figures (cylinders, spheres, cones, semi-planes).
    # See cloud_gemoetry.jl for a more detailed description
    #
    cloud = Plane(cth, false)

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = Plane(60 * co.kilo, false)

    # Set a homogeneous cloud droplet density and radius.
    # See variable_droplet_radius.jl for how to deal with inhomogeneous clouds.
    # Radius of the scattering particles
    radius = 10e-6

    # Density of the scattering centers.
    nscat = 100 * co.centi^-3

    composition = Homogeneous(params.λ, nscat, radius)

    # Define the full simulation world
    world = World(cloud, domain, composition)

    # To capture all exiting photons we set a few observers around the closest
    # point in th cloud top to the source.  If the radiance would be perfectly
    # lambertian this would not be neccesary.
    μ, w = gausslegendre(5)

    observers = @. observerat(acos(0.5 * (μ + 1)), cth)

    basename = splitext(splitdir(@__FILE__)[2])[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    CloudScat.main(params, world, observers,
                   saveto="$(basename).h5")
end

function observerat(θ, zc)
    R = 400 * co.kilo

    Observer(
        # Location of the observer
        position = [0, R * sin(θ), zc + R * cos(θ)],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)
end

# Check if the code has been 'included' or run from the shell.  In the latter
# case, run the simulation.
isinteractive() || run()

ext: jl


text:
include("sample.jl")
include("../statistics.jl")
"""
Generates initial guess for constrained optimization problem, i.e.
finds y value for which lb < constr(y;s) < ub

OUTPUT:
augmented [y, s] vector of value and
"""
function init_constrained_pt(cdf_fixed, lx, ux; initial = [], quantile = 0.25)
    if length(initial)!=0
        loc = initial #use initial point
    else
        (_, loc) = sample_points( _ -> 1.0, lx[2:end], ux[2:end]; num = 1) #samples random location in box
    end
    #(_, loc) = sample( _ -> 1.0, lx[2:end], ux[2:end]; num = 1)
    cdf_y(y) = (arg = Base.vcat([y], reshape(loc, length(loc))); cdf_fixed(arg))
    #@info "loc" loc
    #@info "type of cdf_y", typeof(cdf_y)
    @info "loc", loc
    @info "lx", lx
    @info "ux", ux
    #try
    (ystar, _, _) = findQuantile(cdf_y, [0, 5000]; p = quantile)
    #catch e
    #    println("findQuantile Error...")
    #    @info "cdf_y(2000)", cdf_y(2000)
    #    @info "cdf_y(0)", cdf_y(0)
    #finally
    println("successfully ran findQuantile...")
    vec = vcat(ystar, reshape(loc, length(loc)))
    @assert typeof(vec)<:Array{T, 1} where T
    return vec
    #end
end

function findQuantile(cdf::Function, support::Array{G,1} where G; pdf = nothing, pdf_deriv=nothing, p=.25)
    bound = support
    #@info "bound", bound
    func = y0 -> cdf(y0) - p
    #@info "func0: ", func(0)
    #@info "func1: ", func(800)
    ystar = fzero(func, bound[1], bound[2])
    relerr = abs(p-cdf(ystar))/p
    @assert relerr < 1 e-1
    return ystar, relerr, bound
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
export
	SolverStatus,
	time_out_status,
	iteration_max_out_status,
	failed_status,
	set_time_out_status!,
	set_iteration_max_out_status!

@with_kw mutable struct SolverStatus
	timed_out::Bool = false
	iterations_maxed_out::Bool = false
end
solver_status(status::SolverStatus) = status
time_out_status(status) 			= solver_status(status).timed_out
iteration_max_out_status(status) 	= solver_status(status).iterations_maxed_out
failed_status(status) = time_out_status(status) || iteration_max_out_status(status)
function set_time_out_status!(status,val=true)
	solver_status(status).timed_out = val
end
function set_iteration_max_out_status!(status,val=true)
	solver_status(status).iterations_maxed_out = val
end

export SolverLogger

"""
    SolverLogger

A logger type for keeping track of thing like runtime, iterations, optimality
gap (including upper and lower bound), etc. The following methods allow for
accessing and modifying the `SolverLogger`'s fields, and are extended to solver
types that wrap a `SolverLogger`:

iterations(logger)
iteration_limit(logger)
max_iterations(logger)
start_time(logger)
runtime_limit(logger)
deadline(logger)
JuMP.lower_bound(logger)
best_cost(logger)
verbosity(logger)
verbosity(val::Int)
debug(logger)
solver_status(logger)

set_iterations!(solver,val)
increment_iteration_count!(logger::SolverLogger)
set_iteration_limit!(solver,val)
set_max_iterations!(solver,val)
set_start_time!(solver,val)
set_runtime_limit!(solver,val)
set_deadline!(solver,val)
set_lower_bound!(solver,val)
set_lower_bound!(logger::SolverLogger{C},val::C) where {C}
set_lower_bound!(logger::SolverLogger{NTuple{N,T}},val::R) where {N,T<:Real,R<:Real}
set_lower_bound!(logger::SolverLogger{T},val::R) where {T<:Tuple,R<:Real}
set_best_cost!(solver,val)
set_best_cost!(logger::SolverLogger,val)
set_best_cost!(logger::SolverLogger{NTuple{N,T}},val::R) where {N,T<:Real,R<:Real}
set_best_cost!(logger::SolverLogger{T},val::R) where {T<:Tuple,R<:Real}
set_verbosity!(solver,val)
set_debug!(solver,val)

The following methods facilitate control flow based on solver status.
check_time(logger)
enforce_time_limit!(logger)
check_iterations(logger)
enforce_iteration_limit!(logger)
"""
@with_kw mutable struct SolverLogger{C}
    iterations      ::Int           = 0
    iteration_limit ::Int           = 100
    max_iterations  ::Int           = 0
    start_time      ::Float64       = time()
    runtime_limit   ::Float64       = 100.0
    deadline        ::Float64       = Inf
    lower_bound     ::C             = typemin(C)
    best_cost       ::C             = typemax(C)
    verbosity       ::Int           = 0
    DEBUG           ::Bool          = false
	status			::SolverStatus	= SolverStatus()
end
cost_type(logger::SolverLogger{C}) where {C} = C
cost_type(solver) = cost_type(get_logger(solver))

export get_logger

get_logger(solver) = solver.logger
get_logger(logger::SolverLogger) = logger

export SolverWrapper

"""
	SolverWrapper

An abstract type whose concrete instances must have a `solver` field.
"""
abstract type SolverWrapper end
get_logger(solver::SolverWrapper) = get_logger(solver.solver)
low_level(solver::SolverWrapper) = low_level(solver.solver)

export
    iterations,
    max_iterations,
    iteration_limit,
    start_time,
    runtime_limit,
    deadline,
    best_cost,
    verbosity,
    debug,
    time_to_deadline

iterations(logger)      = get_logger(logger).iterations
iteration_limit(logger) = get_logger(logger).iteration_limit
max_iterations(logger)  = get_logger(logger).max_iterations
start_time(logger)      = get_logger(logger).start_time
runtime_limit(logger)   = get_logger(logger).runtime_limit
deadline(logger)        = get_logger(logger).deadline
JuMP.lower_bound(logger)= get_logger(logger).lower_bound
best_cost(logger)       = get_logger(logger).best_cost
verbosity(logger)       = get_logger(logger).verbosity
verbosity(val::Int)		= val
verbosity(::Nothing)    = global_verbosity()
debug(logger)           = get_logger(logger).DEBUG
solver_status(logger) 	= get_logger(logger).status

"""
    time_to_deadline(solver)

time to `deadline(solver)` or `runtime_limit(solver)`--whichever is shorter.
"""
time_to_deadline(logger) = min(runtime_limit(logger),deadline(logger)-time())

export
	set_iterations!,
	set_iteration_limit!,
	set_max_iterations!,
	set_start_time!,
	set_runtime_limit!,
	set_deadline!,
	set_lower_bound!,
	set_best_cost!,
	set_verbosity!,
	set_debug!

function set_iterations!(solver,val)
	get_logger(solver).iterations = val
end
function set_iteration_limit!(solver,val)
	get_logger(solver).iteration_limit = val
end
function set_max_iterations!(solver,val)
	get_logger(solver).max_iterations = val
end
function set_start_time!(solver,val)
	get_logger(solver).start_time = val
end
function set_runtime_limit!(solver,val)
	get_logger(solver).runtime_limit = val
end
function set_deadline!(solver,val)
	get_logger(solver).deadline = val
end
function set_lower_bound!(solver,val)
	get_logger(solver).lower_bound = val
end
function set_lower_bound!(logger::SolverLogger{C},val::C) where {C}
    logger.lower_bound = val
end
function set_lower_bound!(logger::SolverLogger{NTuple{N,T}},val::R) where {N,T<:Real,R<:Real}
    logger.lower_bound = NTuple{N,T}((T(val),zeros(T,N-1)...))
end
function set_lower_bound!(logger::SolverLogger{T},val::R) where {T<:Tuple,R<:Real}
    v = cost_type(logger)([val,zeros(length(lower_bound(logger))-1)...])
    set_lower_bound!(logger,v)
end
function set_best_cost!(solver,val)
	set_best_cost!(get_logger(solver), val)
end
function set_best_cost!(logger::SolverLogger,val)
    logger.best_cost = val
end
function set_best_cost!(logger::SolverLogger{NTuple{N,T}},val::R) where {N,T<:Real,R<:Real}
    logger.best_cost = NTuple{N,T}((T(val),zeros(T,N-1)...))
end
function set_best_cost!(logger::SolverLogger{T},val::R) where {T<:Tuple,R<:Real}
    v = cost_type(logger)([val,zeros(length(lower_bound(logger))-1)...])
    set_best_cost!(logger,v)
end
function set_verbosity!(solver,val)
	get_logger(solver).verbosity = val
end
function set_debug!(solver,val)
	get_logger(solver).DEBUG = val
end
for op in [:set_time_out_status!,:set_iteration_max_out_status!,
		:time_out_status,:iteration_max_out_status]
	@eval $op(logger,args...) = $op(get_logger(logger).status,args...)
end

export SolverException

"""
    SolverException

Custom exception type for tracking solver timeouts, etc.
"""
struct SolverException <: Exception
    msg::String
end
SolverException(msg::String,msg2,args...) = SolverException(string(msg,msg2,args...))

export handle_solver_exception

"""
	handle_solver_exception

Takes care of printing `SolverException`s
"""
function handle_solver_exception(solver,e)
    if debug(solver)
        showerror(stdout, e, catch_backtrace())
    else
        printstyled(e,"\n";color=:red)
    end
end

export
	solver_type

solver_type(::SolverLogger) = "Solver"

export
    optimality_gap,
	feasible_status,
	optimal_status

optimality_gap(logger) = best_cost(logger) .- lower_bound(logger)
feasible_status(solver) = best_cost(solver) < typemax(cost_type(solver))
optimal_status(solver) 	= optimality_gap(solver) == 0

export
	check_time,
	check_iterations,
	enforce_time_limit!,
	enforce_iteration_limit!,
    increment_iteration_count!

function check_time(logger)
    t = time()
    if t >= deadline(logger) || t - start_time(logger) >= runtime_limit(logger)
        return true
    end
    return false
end
function enforce_time_limit!(logger)
    if check_time(logger)
		set_time_out_status!(logger,true)
        throw(SolverException("$(solver_type(logger)) time limit exceeded!",
		" deadline was $(deadline(logger)),",
		" runtime_limit was $(runtime_limit(logger))"))
    end
end
function check_iterations(logger)
    if iterations(logger) >= iteration_limit(logger)
		return true
	end
	return false
end
function enforce_iteration_limit!(logger)
    if check_iterations(logger)
		set_iteration_max_out_status!(logger,true)
        throw(SolverException("Solver iterations exceeded! Limit was $(iteration_limit(logger))."))
    end
end

function increment_iteration_count!(logger::SolverLogger)
    logger.iterations += 1
    set_max_iterations!(logger,max(iterations(logger),max_iterations(logger)))
    logger.iterations
end
increment_iteration_count!(solver)  = increment_iteration_count!(get_logger(solver))

export
    reset_solver!,
    soft_reset_solver!,
    hard_reset_solver!

"""
    reset_solver!(solver)

Resets iteration counts and start times, in addition to best cost and lower bound.
"""
function reset_solver!(logger::SolverLogger)
    set_iterations!(logger, 0)
    set_best_cost!(logger,typemax(cost_type(logger)))
    set_lower_bound!(logger,typemin(cost_type(logger)))
    set_start_time!(logger,time())
	set_time_out_status!(logger,false)
	set_iteration_max_out_status!(logger,false)
    logger
end
reset_solver!(solver)               = reset_solver!(get_logger(solver))

"""
    soft_reset_solver!(solver)

Resets iteration counts and start times.
"""
function soft_reset_solver!(logger::SolverLogger)
    set_iterations!(logger, 0)
    set_start_time!(logger,time())
	set_time_out_status!(logger,false)
	set_iteration_max_out_status!(logger,false)
    logger
end
soft_reset_solver!(solver)               = soft_reset_solver!(get_logger(solver))

"""
    hard_reset_solver!(solver)

To be called when no information (other than iteration and time limits) needs to
be stored.
"""
function hard_reset_solver!(logger::SolverLogger)
    reset_solver!(logger)
    set_max_iterations!(logger,0)
end
hard_reset_solver!(solver)          = hard_reset_solver!(get_logger(solver))

# get_infeasible_cost(logger::SolverLogger{C}) where {C} = typemax(C)
# get_infeasible_cost(solver) = get_infeasible_cost(get_logger(solver))

# export
#     global_verbosity,
#     set_global_verbosity!
#
# global VERBOSITY = 0
#
# """
# 	global_verbosity()
#
# Query the global verbosity setting (VERBOSITY ∈ [0,])
# """
# global_verbosity() = copy(VERBOSITY)
#
# """
# 	global_verbosity()
#
# Set the global verbosity to val
# """
# set_global_verbosity!(val::Int) = begin global VERBOSITY = val end
#
# function print_styled_header(header,msg...;bold=true,color=136)
# 	printstyled(header;bold=bold,color=color)
# 	println(msg...)
# end
#
# export @log_info
#
# """
#     @log_info
#
# A helper macro for printing at various verbosity levels.
# Usage:
# 	`@log_info(limit::Int,level::Int,msg...)`
# 	`@log_info(limit::Int,solver,msg...)`
# Args:
# * limit::Int - logging threshold
# * level::Int - the verbosity level
# * msg... - message to print if level > limit
# """
# macro log_info(limit,level,msg...)
# 	filename = split(string(__source__.file),"/")[end]
# 	ex = :(print_styled_header(log_preamble))
#     for x in esc(msg).args[1]
# 		push!(ex.args,esc(x))
#     end
# 	return quote
# 		local log_preamble = string("[ logger ](",$filename,"-",$(__source__.line),"): ")
# 		verbosity($(esc(level))) > $(esc(limit)) || global_verbosity() > $(esc(limit)) ?
# 			$ex : nothing
# 	end
# end


text:
using Test
using LCMCore
using StaticArrays
using Dates: Second, Millisecond
using FileWatching: poll_fd
import LCMCore: encode, decode

include("mymessage.jl")

@testset "close multiple times" begin
    lcm = LCM()
    close(lcm)
    close(lcm)
    close(lcm)
end

@testset "isgood" begin
    lcm = LCM()
    @test isgood(lcm)
    close(lcm)
    @test !isgood(lcm)
end

@testset "publish raw data" begin
    LCM() do lcm
        publish(lcm, "CHANNEL_0", UInt8[1,2,3,4])
    end
end

@testset "publish and subscribe" begin
    lcm = LCM()

    did_check = false
    data = UInt8[1,2,3,4,5]
    channel = "CHANNEL_1"
    function check_data(c, d)
        did_check = true
        @test c == channel
        @test d == data
    end
    subscribe(lcm, channel, check_data)
    publish(lcm, channel, data)
    handle(lcm)
    @test did_check
end

@testset "async publish and subscribe" begin
    lcm = LCM()

    data = UInt8[1,2,3,4,5]
    channel = "CHANNEL_1"
    function check_data(c, d)
        did_check = true
        @test c == channel
        @test d == data
    end
    subscribe(lcm, channel, check_data)
    @async handle(lcm)
    publish(lcm, channel, data)
end

@testset "encode and decode" begin
    lcm = LCM()
    msg = MyMessage(23, 1.234)
    did_check = false
    channel = "CHANNEL_2"
    function check_data(c, d)
        did_check = true
        @test d.field1 == msg.field1
        @test d.field2 == msg.field2
    end
    subscribe(lcm, channel, check_data, MyMessage)
    publish(lcm, channel, msg)
    handle(lcm, Second(1))
    @test did_check

    # Test handle() with a timeout when there's no message available
    did_check = false
    handle(lcm, Second(1))
    @test did_check == false
end

@testset "queue capacity 1" begin
    lcm = LCM()

    did_check = false
    channel = "CHANNEL_1"
    function check_data(c, d)
        did_check = true
    end
    sub = subscribe(lcm, channel, check_data)
    @test set_queue_capacity(sub, 1)
    fd = filedescriptor(lcm)
    publish(lcm, channel, UInt8[1,2,3])
    publish(lcm, channel, UInt8[1,2,3,4])

    # Subscription capacity is 1, so the queue should only have 1 message
    # even though we published twice
    @test LCMCore.get_queue_size(sub) == 1

    # Handle the 1 message in queue and verify that the queue is now empty
    LCMCore.lcm_handle(lcm)
    @test did_check
    @test LCMCore.get_queue_size(sub) == 0

    # We published twice and handled once. Because the queue capacity is only
    # 1, there should not be another message available to read.
    event = poll_fd(fd, 1; readable=true)
    @show event
    @test !event.readable
end

@testset "queue capacity 2" begin
    lcm = LCM()

    did_check = false
    channel = "CHANNEL_1"
    function check_data(c, d)
        did_check = true
    end
    sub = subscribe(lcm, channel, check_data)
    set_queue_capacity(sub, 2)
    fd = filedescriptor(lcm)
    publish(lcm, channel, UInt8[1,2,3,4,5])
    publish(lcm, channel, UInt8[1,2,3,4,5])
    @test LCMCore.get_queue_size(sub) == 2
    LCMCore.lcm_handle(lcm)
    @test did_check

    # We published twice and handled once. The queue capacity is 2, so
    # there should be another message available to read.
    @test LCMCore.get_queue_size(sub) == 1
    event = poll_fd(fd, 1; readable=true)
    @test event.readable

    # Now we handle the second message in the queue
    did_check = false
    LCMCore.lcm_handle(lcm)
    @test did_check

    # Now we've exhausted the queue, so there should be no more messages
    # available to read
    @test LCMCore.get_queue_size(sub) == 0
    event = poll_fd(fd, 1; readable=true)
    @test !event.readable
end

@testset "unsubscribe" begin
    lcm = LCM()

    channel = "FOO"
    did_callback1 = false
    function callback1(channel, data)
        did_callback1 = true
    end
    did_callback2 = false
    function callback2(channel, data)
        did_callback2 = true
    end
    sub1 = subscribe(lcm, channel, callback1)
    unsubscribe(lcm, sub1)
    sub2 = subscribe(lcm, channel, callback2)
    publish(lcm, channel, UInt8[1,2,3])
    handle(lcm)

    @test !did_callback1
    @test did_callback2
end

@testset "UDP ports" begin
    lcm1 = LCM("udpm://239.255.76.67:7667")
    lcm2 = LCM("udpm://239.255.76.67:7668")

    channel = "FOO"
    did_callback1 = false
    function callback1(channel, data)
        @test data == UInt8[1, 2, 3]
        did_callback1 = true
    end
    did_callback2 = false
    function callback2(channel, data)
        @test data == UInt8[1, 2, 3, 4]
        did_callback2 = true
    end
    subscribe(lcm1, channel, callback1)
    subscribe(lcm2, channel, callback2)
    publish(lcm1, channel, UInt8[1,2,3])
    handle(lcm1, Millisecond(100))
    handle(lcm2, Millisecond(100))
    @test did_callback1
    @test !did_callback2

    did_callback1 = false
    did_callback2 = false
    publish(lcm2, channel, UInt8[1,2,3,4])
    handle(lcm1, Millisecond(100))
    handle(lcm2, Millisecond(100))
    @test !did_callback1
    @test did_callback2
end

@testset "lcm_handle allocations" begin
    data = UInt8[1,2,3,4,5]
    channel = "CHANNEL_1"

    # start listening
    sublcm = LCM()
    sub = subscribe(sublcm, channel, (c, d) -> nothing)
    set_queue_capacity(sub, 2)

    # publish two messages
    publcm = LCM()
    for _ = 1 : 2
        publish(publcm, channel, data)
    end

    # check that handling doesn't allocate
    LCMCore.lcm_handle(sublcm)
    allocs = @allocated LCMCore.lcm_handle(sublcm)
    @test allocs == 0
end

include("test_lcmtype.jl")
include("test_readlog.jl")


text:
# Example 2.8 of
# P. Parrilo and A. Jadbabaie
# Approximation of the joint spectral radius using sum of squares
# Linear Algebra and its Applications, Elsevier, 2008, 428, 2385-2402
# Inspired from the construction of:
# Ando, T. and Shih, M.-h.
# Simultaneous Contractibility.
# SIAM Journal on Matrix Analysis & Applications, 1998, 19, 487
# The JSR is 1

@testset "[PJ08] Example 2.8" begin
    expected_lb = [1, 1/2^(1/4)]
    expected_ub = [√2, 1]
    for optimizer_constructor in sdp_factories
        s = discreteswitchedsystem([[1 0; 1 0], [0 1; 0 -1]])
        tol = ismosek(optimizer_constructor) ? 2e-4 : 1e-3
        for d in 1:2
            # We set lb=0 for d=2, otherwise there might be no infeasibile problem and hence no extraction done
            sosdata(s).lb = 0
            lb, ub = soslyapb(s, d, optimizer_constructor=optimizer_constructor, tol=tol)
            @test abs(log(expected_lb[d]) - log(lb)) <= tol
            @test abs(log(expected_ub[d]) - log(ub)) <= tol

            if d == 1
                @test !hassmp(s)
                @test_throws ErrorException getsmp(s)
            else
                @test getlb(s) ≈ 1
                @test hassmp(s)
                @test getsmp(s).growthrate == 1
            end
        end

        for d in 1:2
            seq = sosbuildsequence(s, d, p_0=:Primal)
            psw = findsmp(seq)
            @test psw !== nothing
            @test psw.growthrate == 1
        end
    end
end


text:
using Test
using libcxxwrap_julia_jll

excluded = ["runtests.jl", "testcommon.jl"]

if Sys.iswindows() && Sys.WORD_SIZE == 32
  push!(excluded, "except.jl")
end

@testset "CxxWrap tests" begin
  for f in filter(fname -> fname ∉ excluded, readdir())
    println("Running tests from $f...")
    include(f)
  end
end


text:
# # Generating and manipulating a tetrahedral mesh

# In this tutorial, we will learn
#
#    -  How to generate a simple mesh.
#    -  How to inspect the mesh.
#    -  How to inspect the data stored in the mesh.
#    -  How to export the mesh and visualize it.

# The tutorial comes with a file containing representation of the mesh for
# visualization in the [Paraview](https://www.paraview.org/) format (VTK). One can
# visualize this mesh by loading the file `mymesh.vtu` with `paraview.exe`.
# One of the aims of this tutorial is to work with this mesh: to generate it,
# inspect it, and eventually to store it in the VTK format.

# We will generate the tetrahedral mesh inside a rectangular block.
# The block will have the dimensions shown below:
a, b, c = 2.0, 2.5, 3.0
# The tetrahedra will be generated in a regular pattern, with the number of
# edges per side of the block given as
na, nb, nc = 2, 2, 3

# The mesh will be generated by the package `MeshSteward`.
using MeshSteward: T4block
conn = T4block(a, b, c, na, nb, nc);

# The variable `conn` is an incidence relation. This will become the base
# relation of the mesh.
using MeshSteward: Mesh, attach!
# The mesh is first created.
m = Mesh()
# Then the ``(3, 0)`` incidence relation, which defines the tetrahedral elements in terms of the vertices at their corners, is attached to it.
attach!(m, conn);

# We can now inspect the mesh by printing its summary.
println(summary(m))

# We can see that the relation links the tetrahedral elements and the vertices.
# In addition to its being the only relation in the mesh (at the moment), it is
# also what is called the *base* incidence relation.

# The relations in the mesh are accessed by code. The base relation  is
# described by the code:
using MeshSteward: basecode
@show irc = basecode(m);

# We can retrieve the base incidence relation from the mesh as
using MeshSteward: increl
conn = increl(m, basecode(m));

# We can access the data stored in the incidence relation as follows. For
# instance, the connectivity of the mesh can be accessed with regular array
# indexing. The numbers of the vertices for the first tetrahedron are:
@show conn[1]

# We can check that the incidence relation stores four vertices per tetrahedron as
using MeshCore: nentities
@show nentities(conn, 1) == 4
# We checked it for the first tetrahedron, but for this type of incidence
# relation all relations store the same number of entities.

# The vertex shape collection stores the locations of the vertices as an
# attribute.
using MeshCore: attribute
@show geom = attribute(conn.right, "geom");

# The coordinates of the vertices of the first tetrahedron can be accessed as
for j in 1:nentities(conn, 1)
    k = conn[1, j]
    println("Vertex $(j): global number $(k)")
    println("   $(geom[k])")
end

# The mesh may be exported for viewing with the "Paraview" visualization
# program. Note that a file of this name is provided with the installation.
# When you run this tutorial  you should get the same file. So a quick way of
# checking that things are working as expected is to  visualize what you
# should get by loading the file `mymesh.vtu` with `paraview.exe`, and only
# then running the export in the two lines below and looking at the results again.
using MeshSteward: vtkwrite
vtkwrite("mymesh", conn)

# Start "Paraview", load the file `mymesh.vtu` and select for instance view as
# "Surface with Edges".

# Provided `paraview.exe` is installed and in the path where your operating
# system searches for executables, you may be able to run the following command
# to bring up the visualization program.
@async run(`paraview mymesh.vtu`)

# The executable will load the mesh. The mesh can then be visualized in
# different ways, for instance as a surface with edges.


text:
using Luxor

Drawing(400, 400, joinpath(dirname(@__DIR__), "images", "logo.png"))
origin()
setline(5)
sethue("white")
squircle(O, 185, 185, :fillpreserve)
sethue("black")
strokepath()

@layer begin
    translate(0, 30)
    fontface("AlegreyaSans-BoldItalic")
    fontsize(70)
    text("HUMANS", O - (0, 100), halign=:center)
    fontsize(40)
    text("OF", Point(-65, -50))
    julialogo(centered=true)
end

finish()

max_stars_repo_name: Humans-of-Julia/HoJ-logo-graphics

max_issues_repo_name: Humans-of-Julia/HoJ-logo-graphics

max_forks_repo_name: Humans-of-Julia/HoJ-logo-graphics


text:
using SchroedingersSmoke

vol_size = (4,2,2)# box size
dims = (64, 32, 32) .* 2 # volume resolution
hbar = 0.1f0      # Planck constant
dt = 1f0/48f0     # time step

jet_velocity = (1f0, 0f0, 0f0)
nozzle_cen = Float32.((2-1.7, 1-0.034, 1+0.066))
nozzle_len = 0.5f0
nozzle_rad = 0.5f0
n_particles = 1000   # number of particles

ArrayType = Array

isf2 = ISF{ArrayType, UInt32, Float32}(vol_size, dims, hbar, dt);

# initialize psi
psi = ArrayType([(one(Complex64), one(Complex64) * 0.01f0) for i=1:dims[1], j=1:dims[2], k=1:dims[3]]);

psi .= normalize_psi.(psi);

kvec = jet_velocity ./ hbar;
omega = sum(jet_velocity.^2f0) / (2f0*hbar);
isjetarr = isjet.(isf2.positions, (nozzle_cen,), (nozzle_len,), (nozzle_rad,))
# constrain velocity
for iter = 1:10
    restrict_velocity!(isf2, psi, kvec, isjetarr, 0f0)
    pressure_project!(isf2, psi)
end

# add out of bounds particles
particles = ArrayType(map(x-> (-1f0, -1f0, -1f0), 1:100_000))

add_particles!(particles, 1:n_particles, nozzle_cen, nozzle_rad)

function simloop(
        N, isf, psi, kvec, omega, n_particles, isjetarr,
        nozzle_rad, nozzle_cen, particle
    )
    dt = isf.dt; d = isf.d
    for iter = 1:N
        t = iter * dt
        # incompressible Schroedinger flow
        schroedinger_flow!(isf, psi)
        psi .= normalize_psi.(psi)
        pressure_project!(isf, psi)

        start = mod((iter - 1) * n_particles + 1, length(particles))
        stop = start + n_particles - 1
        add_particles!(particles, start:stop, nozzle_cen, nozzle_rad)

        # constrain velocity
        restrict_velocity!(isf, psi, kvec, isjetarr, omega*t)
        pressure_project!(isf, psi)
        velocity_one_form!(isf, psi, isf.hbar)
        # inplace StaggeredSharp
        dinv = inv.(isf.d)
        broadcast!((x, y)-> x .* y, isf.velocity, isf.velocity, (dinv,))
        staggered_advect!(particle, isf)
    end
    GPUArrays.synchronize(particle)
end

@time simloop(
    3, isf2, psi, kvec, omega, n_particles, isjetarr,
    nozzle_rad, nozzle_cen, particles
)


text:
# Test internals

@testset "Internal tests" begin
    for jtype in ADIOS2.julia_types
        @test ADIOS2.julia_type(ADIOS2.adios_type(jtype)) ≡ jtype
    end
    for atype in ADIOS2.adios_types
        @test ADIOS2.adios_type(ADIOS2.julia_type(atype)) ≡ atype
    end
end


text:
using CLIMA.PlanetParameters
export NoRadiation, DYCOMSRadiation

abstract type RadiationModel end

vars_state(::RadiationModel, FT) = @vars()
vars_aux(::RadiationModel, FT) = @vars()
vars_integrals(::RadiationModel, FT) = @vars()

function atmos_nodal_update_aux!(::RadiationModel, ::AtmosModel, state::Vars,
                                 aux::Vars, t::Real)
end
function preodefun!(::RadiationModel, aux::Vars, state::Vars, t::Real)
end
function integrate_aux!(::RadiationModel, integ::Vars, state::Vars, aux::Vars)
end
function flux_radiation!(::RadiationModel, atmos::AtmosModel, flux::Grad, state::Vars,
                         aux::Vars, t::Real)
end

struct NoRadiation <: RadiationModel
end

"""
  DYCOMSRadiation <: RadiationModel

Stevens et. al (2005) approximation of longwave radiative fluxes in DYCOMS.
Analytical description as a function of the liquid water path and inversion height zᵢ

* Stevens, B. et. al. (2005) "Evaluation of Large-Eddy Simulations via Observations of Nocturnal Marine Stratocumulus". Mon. Wea. Rev., 133, 1443–1462, https://doi.org/10.1175/MWR2930.1

# Fields

$(DocStringExtensions.FIELDS)
"""
struct DYCOMSRadiation{FT} <: RadiationModel
  "mass absorption coefficient `[m^2/kg]`"
  κ::FT
  "Troposphere cooling parameter `[m^(-4/3)]`"
  α_z::FT
  "Inversion height `[m]`"
  z_i::FT
  "Density"
  ρ_i::FT
  "Large scale divergence `[s^(-1)]`"
  D_subsidence::FT
  "Radiative flux parameter `[W/m^2]`"
  F_0::FT
  "Radiative flux parameter `[W/m^2]`"
  F_1::FT
end
vars_integrals(m::DYCOMSRadiation, FT) = @vars(attenuation_coeff::FT)
vars_aux(m::DYCOMSRadiation, FT) = @vars(Rad_flux::FT)
function integrate_aux!(m::DYCOMSRadiation, integrand::Vars, state::Vars, aux::Vars)
  FT = eltype(state)
  integrand.radiation.attenuation_coeff = state.ρ * m.κ * aux.moisture.q_liq
end
function flux_radiation!(m::DYCOMSRadiation, atmos::AtmosModel, flux::Grad, state::Vars,
                         aux::Vars, t::Real)
  FT = eltype(flux)
  z = altitude(atmos.orientation, aux)
  Δz_i = max(z - m.z_i, -zero(FT))
  # Constants
  upward_flux_from_cloud  = m.F_0 * exp(-aux.∫dnz.radiation.attenuation_coeff)
  upward_flux_from_sfc = m.F_1 * exp(-aux.∫dz.radiation.attenuation_coeff)
  free_troposphere_flux = m.ρ_i * FT(cp_d) * m.D_subsidence * m.α_z * cbrt(Δz_i) * (Δz_i/4 + m.z_i)
  F_rad = upward_flux_from_sfc + upward_flux_from_cloud + free_troposphere_flux
  ẑ = vertical_unit_vector(atmos.orientation, aux)
  flux.ρe += F_rad * ẑ
end
function preodefun!(m::DYCOMSRadiation, aux::Vars, state::Vars, t::Real)
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
@symbol_func function phi_res(cur_reactor::AbstractReactor)
  cur_phi = cur_reactor.C_ejima

  cur_phi *= cur_reactor.R_0

  cur_phi *= ( 4e-7 * cur_reactor.pi )

  cur_phi *= f_IN(cur_reactor) * cur_reactor.I_P

  cur_phi *= 1e6

  cur_phi
end


text:
@testset "eltype" begin
    @test eltype_or(HyperCube, nothing) == nothing
    @test eltype_or(HyperCube{2}, nothing) == nothing
    @test eltype_or(HyperCube{2, Float32}, nothing) == Float32
    @test eltype_or(SimpleRectangle, Int) == Int
    @test eltype_or(SimpleRectangle{Float32}, Int) == Float32

    @test eltype(SimpleRectangle(0,0,1,1)) == Int
    @test eltype(SimpleRectangle) == Any
    @test eltype(HyperCube{2}) == Any
    @test eltype(HyperCube{2, Float32}) == Float32
    @test eltype(SimpleRectangle{Float32}) == Float32
end
@testset "ndims" begin
    @test ndims_or(HyperCube, nothing) == nothing
    @test ndims_or(HyperCube{2}, nothing) == 2
    @test ndims_or(HyperCube{2, Float32}, nothing) == 2
    @test ndims_or(SimpleRectangle, 0) == 2

    @test ndims(SimpleRectangle(0,0,1,1)) == 2
    @test ndims(SimpleRectangle) == 2
    @test ndims(HyperCube{2}) == 2
    @test ndims(HyperCube{2, Float32}) == 2
    #@test_throws Exception ndims(HyperCube)
end


text:
@testset "IterationCandidates" begin
    ic = RANSAC.IterationCandidates()
    fp = FittedPlane(SVector(0.5,0.5,0.5), SVector{3}(0,0,1.))
    sc = ConfidenceInterval(0,1)
    inds = [1,2,3,4,5]

    @test length(ic) == 0
    RANSAC.recordscore!(ic, fp, sc, inds)
    @test length(ic) == 1
    @test size(ic.shapes,1) == 1
    @test size(ic.scores,1) == 1
    @test size(ic.inpoints,1) == 1

    deleteat!(ic, 1)
    @test length(ic) == 0
    @test size(ic.shapes,1) == 0
    @test size(ic.scores,1) == 0
    @test size(ic.inpoints,1) == 0
end


text:
export MinMetricGoal

"Goal specification where costs are differences in a state-based metric."
struct MinMetricGoal <: Goal
    terms::Vector{Term} # Goal terms to be satisfied
    metric::Term # Cost metric to be minimized
end

function MinMetricGoal(problem::Problem)
    goals = flatten_conjs(problem.goal)
    sign, metric = problem.metric
    if sign > 0 metric = Compound(:-, [metric]) end
    return MinMetricGoal(goals, metric)
end
MinMetricGoal(goal::Term, metric::Term) =
    MinMetricGoal(flatten_conjs(goal), metric)

Base.hash(spec::MinMetricGoal, h::UInt) =
    hash(spec.metric, hash(Set(spec.terms), h))
Base.isequal(s1::MinMetricGoal, s2::MinMetricGoal) =
    s1.metric == s2.metric && Set(s1.terms) == Set(s2.terms)

is_goal(spec::MinMetricGoal, domain::Domain, state::State) =
    satisfy(spec.terms, state, domain)[1]
is_violated(spec::MinMetricGoal, domain::Domain, state::State) = false
get_cost(spec::MinMetricGoal, domain::Domain, s1::State, ::Term, s2::State) =
    s2[domain, spec.metric] - s1[domain, spec.metric]
get_reward(spec::MinMetricGoal, domain::Domain, s1::State, ::Term, s2::State) =
    s1[domain, spec.metric] - s2[domain, spec.metric]
get_goal_terms(spec::MinMetricGoal) = spec.terms


text:
module SurfaceQuasiGeostrophic

using FFTW

include("grid.jl")
include("sqg.jl")
include("poisson.jl")
include("init_buoyancy.jl")
include("compute_dt.jl")
include("compute_hypercoef.jl")

end


text:
#module Code1

## Pkg.add("Distributions") # If you don't already have it installed.
using Distributions
using Statistics

struct Bounds
    min
    max
end

struct LinearCoefficients
    intercept
    slopes
end

# Example from https://github.com/giob1994/Alistair.jl
n = 100
b = Bounds(0, 20)
β0 = 0.7
β = [0.45]
sd_e = 1

X1 = rand(Uniform(xmin, xmax), n, 1)
X0 = ones(n, 1)
X = [X0 X1]
e = rand(Normal(0, sd_e), n, 1)
Y = X * β + e

function uniform_x(k, b::Bounds, n=n::Int)
   return rand(Uniform(b.min, b.max), n, k)
end

function nrows(X::Array):
    return size(X)[1]

function make_ones_column(X):
    nrows = size(X)[1]
    return ones(nrows, 1)
end

function add_intercept(X::Array)
    X0 = make_ones_column(X)
    return [X0 X]
end

function make_e(;mean=0, sd=1, dist=Normal, n=n)
    return rand(dist(mean, sd), n, 1)
end

function make_x(c::Coefficients, b::Bounds, n::Int):
    k = nrows(c.slopes)
    X = uniform_x(k, b, n)
    if c.intercept != 0:
        X = add_intercept(X)
    end
    return X

function make_beta(c::Coefficients):
    if c.intercept == 0:
        return c.slopes
    else:
        return [c.intercept c.slopes]
    end
end

function make_y(c::Coefficients, X):
    β = make_beta(c)
    return X * β

n = 100
x_bounds = Bounds(0, 20)
lin_coef = Coefficients(0.7, [0.45])

function rnorm(mean=0, sd=1)
    rand(Normal(mean, sd))
end

function make_homoscedastic_normal_distrubance(X; mean=0, sd=1)
    return [rnorm(mean, sd) for i in X[:,1]]
end

# x generator
# x*b
function make_observations_normal(N=100)
    X = make_x(lin_coef, x_bounds, n)
    e = make_homoscedastic_normal_distrubance(X, sd=1)
    Y = make_y(lin_coef, X) + e
    return

# OLS estimation using \
# https://github.com/giob1994/Alistair.jl/blob/3a11c19150169695581b46e4d1895f0641a4c29d/src/linregress.jl#L38-L41
function ols(X, Y)
    # Ax = B
    # x = B \ A
    # https://github.com/JuliaLang/julia/blob/d789231e9985537686052db9b2314c0d51656308/stdlib/LinearAlgebra/src/generic.jl#L827-L855
    return (X' * X) \ (X' * Y)
end

# OLS estimation using inv()
# inv() is a more computationally-intensive procedure than \ operator
# https://github.com/mcreel/Econometrics/blob/b29cffa8ce7f845fa73788522f9c7753d40684f4/Examples/OLS/OlsFit.jl#L13-L14
function ols_(X, Y)
    return inv(X'*X)*X'*Y
end

# the difference between two methods is really small
@assert all(ols(X, Y)-ols_(X,Y) .< 1e-10) # was 14

beta_hat = ols(X, Y)
Y_hat = X * beta_hat
RSS = sum((Y_hat - Y) .^ 2)
TSS = sum((Y .- mean(Y)) .^ 2) # equals var(Y)*n
R2 = 1-(RSS/TSS)
println("R-squared", round(R2, digits=4))

#TODO:
#- https://docs.julialang.org/en/v1/manual/workflow-tips/index.html#REPL-based-workflow-1
# filled an issue: https://github.com/JuliaLang/julia/issues/30190

# DONE:
#-  "\U003B2\U00302" does not make a beta-hat on windows
#   see https://github.com/JuliaLang/julia/issues/30175

struct Observation
    X # enforce n*k
    Y # enforce n*1
end

function add_intercept(X::Array)
    nrows = size(X)[1]
    X0 = ones(nrows, 1)
    return [X0 X]
end

@assert add_intercept(collect(1:3)) == [1 1 1; 1 2 3]'

function add_intercept(obs::Observation)
    X_ = add_intercept(obs.X)
    return Observation(X_, obs.Y)
end

function model(observation, add_intercept=true)
    return LinearModel(observation, add_intercept, beta_hat=nothing)
end
# functions add add_intercept
#X1 = rand(Uniform(xmin, xmax), n, 1)
#X0 = ones(n, 1)
#X = [X0 X1]

mutable struct LinearModel
    data::Observation
    X
    beta_hat # n*(k+1)
end

function evaluate_ols(X, Y, add_intercept=true)
    obs = Observation(X, Y)

    # if add_intercept is True, modify X before evaluation
    beta_hat = ols(X, Y)
    return LinearModel(obs, add_intercept, beta_hat)
end

lm = evaluate_ols(X, Y)

#
# make predicted values based on linear model
#

#function replicate(func, n, k):
#    retrun [func() for i in 1:n, j in 1:k]

#function make_x(n, β::Vector, dist)
#    k = length(β)
#    return replicate(rand(dist), n, k)
#end

function nrows(X::Array)
    return size(X)[1]
end

function make_beta(β, β_0=0)
    if β_0 == 0
        return β
    else
        return [β_0 β]
    end
end


text:
# Interface for using ChainRules.jl

# Base.conj(x::Num) = x   # assuming reals
# Base.complex(x::Num) = x
# Base.float(x::Num) = x

@scalar_rule(^(x::Num, n::Integer), (n==1 ? 1 : n==2 ? 2x : n*x^(n-1), Zero()))

adj(f, z̄::Num, x::Num) = (rrule(f, x)[2](z̄))[2]
adj(f, z̄, x) = adj(f, Num(z̄), Num(x))

Num(x::Real) = x

adj(f, z̄::Num, x, y) = simplify( (rrule(f, x, y)[2])(z̄)[2:end] )
adj(f, z̄, x, y) = adj(f, Num(z̄), Num(x), Num(y))

struct Tangent{T <: Real} <: Real
end

# tangent(var) = Variable(Symbol(value(var), "̇"))
# tangent(var::Num) = Sym{Tangent{Real}}(Symbol(var, "̇"))
tangent(var::Num) = tangent(value(var))
tangent(var::Sym) = Sym{Tangent{Real}}(Symbol(var, "̇"))
tangent(x::Real) = 0  # derivative of a constant

function tangent(eq::Assignment)

    vars = args(eq)
    rhs_tangents = Num.(tangent.(args(eq)))
    lhs_tangent = Num(tangent(lhs(eq)))

    return Assignment(lhs_tangent,
                    frule( (Zero(), rhs_tangents...),
                    op(eq), vars...)[2]
    )
end

# julia> frule((Zero(), ẋ, ẏ), *, x, y)
# (x * y, (ẋ * y) + (x * ẏ))

# julia> dargs = Num.(dotify.(args(eq)))

# julia> a = Num.(args(eq))
# 2-element Vector{Num}:
#  x
#  y

# julia> frule( (Zero(), dargs...), op(eq), a...)
# (x * y, (ẋ * y) + (x * ẏ))


text:
using Test
using MultiScaleTreeGraph
using DataFrames
using MutableNamedTuples
using RecipesBase
using Plots # Add this dependency because else the tests on plot recip return an error (I don't know why)
using Graphs

@testset "read_mtg" begin
    include("test-read_mtg.jl")
end

@testset "descendants()" begin
    include("test-descendants.jl")
end

@testset "ancestors()" begin
    include("test-ancestors.jl")
end

@testset "mutation" begin
    include("test-mutation.jl")
    include("test-insert_node.jl")
    include("test-traverse.jl")
    include("test-transform.jl")
end

@testset "plot recipes" begin
    include("test-plot-recipe.jl")
end

@testset "conversion" begin
    include("test-conversion.jl")
end


text:
# O3DFT.jl
#  - Orbital 3 Density Functional Theory

# Marder is: "Condensed Mattery Physics", Michael Marder, 2000.

# Also very useful is Andrei Postnikov's first DFT lecture, on the Thomas-Fermi method
# http://www.home.uni-osnabrueck.de/apostnik/Lectures/DFT-1.pdf

module O3DFT

println("Orbital 3 Density Functional Theory - Philsophy by Numbers")

include("constants.jl")
include("ThomasFermi.jl")

export ThomasFermi_T, ThomasFermi_T_fnderiv, ThomasFermi_Exc, UAtomic
export ThomasFermi_CoulombPotential
export AtomicTest

end


text:
using julia_pkgdev_example
using Test

@testset "julia_pkgdev_example.jl" begin
    # Write your own tests here.
    @test my_f(2,1) == 5
    @test my_f(2,5) == 9

end


text:
function test_lslq()
  lslq_tol = 1.0e-5

  for npower = 1 : 4
    (b, A, D, HY, HZ, Acond, rnorm) = test(40, 40, 4, npower, 0)  # No regularization.

    (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(A, b)
    resid = norm(A' * (A*x - b)) / norm(b)
    @printf("LSLQ: Relative residual: %8.1e\n", resid)
    @test(resid <= lslq_tol)
    @test(stats.solved)

    λ = 1.0e-3
    (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(A, b, λ=λ)
    resid = norm(A' * (A*x - b) + λ * λ * x) / norm(b)
    @printf("LSLQ: Relative residual: %8.1e\n", resid)
    @test(resid <= lslq_tol)
    @test(stats.solved)
  end

  # Test with smallest singular value estimate
  Σ = diagm(0 => 1:4)
  U, _ = qr(rand(6, 6))
  V, _ = qr(rand(4, 4))
  A = U * [Σ ; zeros(2, 4)] * V'
  b = ones(6)
  (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(A, b, σ=1 - 1.0e-10)
  @test isapprox(err_ubnds_lq[end], 0.0, atol=sqrt(eps(Float64)))
  @test isapprox(err_ubnds_cg[end], 0.0, atol=sqrt(eps(Float64)))
  x_exact = A \ b
  @test norm(x - x_exact) ≤ sqrt(eps(Float64)) * norm(x_exact)
  @test norm(x_cg - x_exact) ≤ sqrt(eps(Float64)) * norm(x_exact)

  # Code coverage.
  (b, A, D, HY, HZ, Acond, rnorm) = test(40, 40, 4, 3, 0)
  (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(Matrix(A), b)
  (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(sparse(Matrix(A)), b)
  show(stats)

  # Test b == 0
  (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [I; rand(1:10, 2, 3)]
  b = A * ones(Int, 3)
  (x, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats) = lslq(A, b)
  @test stats.solved
end

test_lslq()


text:
import SeaMass
using Base.Test
using Plots
pyplot()

cd(ENV["PWD"])

# spectrum ID of spectrum we will be plotting
spectrumID = 1

# load mzMLb input spectrum
mzmlbSpectrumIn = SeaMass.MzmlbSpectrum(
  "data/ubiquitin_5uM_1mTorr_1uscan.mzMLb",
  spectrumID
)

# load SMB input spectrum
smbSpectrumIn = SeaMass.SmbSpectrum(
  "data/out/ubiquitin_5uM_1mTorr_1uscan" *
    "/1.mzmlb2smb/ubiquitin_5uM_1mTorr_1uscan.p-0-005377748.smb",
  spectrumID
)
smbSpectrumInBinWidths = smbSpectrumIn.locations[2:end] -
  smbSpectrumIn.locations[1:end-1]

# load SMB output spectrum
smbSpectrumOut = SeaMass.SmbSpectrum(
  "data/out/ubiquitin_5uM_1mTorr_1uscan" *
    "/4.seamass-td/ubiquitin_5uM_1mTorr_1uscan.p-0-005377748.smb",
  spectrumID
)
smbSpectrumOutBinWidths = smbSpectrumOut.locations[2:end] -
  smbSpectrumOut.locations[1:end-1]

# Compare mzMLb input to 'mzmlb2smb' smb input (binned)
plot(
  mzmlbSpectrumIn.mzs,
  mzmlbSpectrumIn.intensities,
  m = 1,
  label = "mzMLb input",
  title = "Compare mzMLb input (sampled) to 'mzmlb2smb' smb input (binned)",
  xlabel = "m/z (Th)",
  ylabel = "ion count density",
  reuse = false,
)
plot!(
  smbSpectrumIn.locations,
  vcat(
    smbSpectrumIn.counts ./ (smbSpectrumIn.exposure .* smbSpectrumInBinWidths),
    0.0
  ),
  line = :steppost,
  label = "SMB input",
)
gui()

# Compare smb input to 'seamass' ⇨ 'seamass-restore' smb output
plot(
  smbSpectrumIn.locations,
  vcat(
    smbSpectrumIn.counts ./ (smbSpectrumIn.exposure .* smbSpectrumInBinWidths),
    0.0
    ),
  line = :steppost,
  label = "SMB input",
  m = 1,
  title = "Compare smb input to 'seamass' ⇨ 'seamass-restore' smb output)",
  xlabel = "m/z (Th)",
  ylabel = "ion count",
  reuse = false,
)
plot!(
  smbSpectrumOut.locations,
  vcat(
    smbSpectrumOut.counts ./ (smbSpectrumOut.exposure .* smbSpectrumOutBinWidths),
    0.0
  ),
  line = :steppost,
  label = "seaMass-restore SMB output",
)
gui()

println("Press <Enter> to finish")
readline()


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LearnedModel(model, θ)

An object that stores `model` together with learned parameters `θ`.
"""
struct LearnedModel
  model
  θ
end

"""
    learn(task, geodata, model)

Learn the `task` with `geodata` using a learning `model`.
"""
function learn(task::AbstractLearningTask, geodata::AbstractData, model)
  if issupervised(task)
    X = view(geodata, collect(features(task)))
    y = vec(geodata[label(task)])
    θ, _, __ = MI.fit(model, 0, X, y)
  else
    X = view(geodata, collect(features(task)))
    θ, _, __ = MI.fit(model, 0, X)
  end

  LearnedModel(model, θ)
end

function learn(task::CompositeTask, geodata::AbstractData, model)
  @error "not implemented"
end

"""
    perform(task, geodata, lmodel)

Perform the `task` with `geodata` using a *learned* `lmodel`.
"""
function perform(task::AbstractLearningTask, geodata::AbstractData, lmodel::LearnedModel)
  # unpack model and learned parameters
  model, θ = lmodel.model, lmodel.θ

  # apply model to the data
  X = view(geodata, collect(features(task)))
  ŷ = MI.predict(model, θ, X)

  # post-process result
  var = outputvars(task)[1]
  if issupervised(task)
    result = isprobabilistic(model) ? mode.(ŷ) : ŷ
  else
    result = ŷ
  end

  OrderedDict(var => result)
end

function perform(task::CompositeTask, geodata::AbstractData, lmodel::LearnedModel)
  @error "not implemented"
end


text:
module BioPredictions

export pI, weight, gravy, cleave

using BioSequences

include("tables.jl")
include("pI.jl")
include("molecularweight.jl")
include("gravy.jl")
include("peptidecleavage.jl")

end


text:
import WiltonInts84

mutable struct WiltonSERule{P,Q} <: SingularityExtractionRule
    outer_quad_points::P
    regularpart_quadrule::DoubleQuadRule{P,Q}
end

function innerintegrals!(op::MWSingleLayer3DSng, p, g, f, t, s, z,
        strat::WiltonSERule, dx)

    γ = op.gamma
    x = cartesian(p)
    n = cross(s[1]-s[3],s[2]-s[3])
    n /= norm(n)
    ρ = x - ((x-s[1]) ⋅ n) * n

    scal, vec = WiltonInts84.wiltonints(s[1], s[2], s[3], x, Val{1})

    # \int \frac{1}{4 \pi R}
    ∫G = (scal[2] - γ*scal[3] + 0.5*γ^2*scal[4]) / (4π)
    # \int \frac{y}{4 \pi R}
    ∫Gy = SVector((
        (vec[2][1] + scal[2]*ρ[1] - γ*(vec[3][1]+scal[3]*ρ[1]) + 0.5*γ^2*(vec[4][1]+scal[4]*ρ[1]))/(4π),
        (vec[2][2] + scal[2]*ρ[2] - γ*(vec[3][2]+scal[3]*ρ[2]) + 0.5*γ^2*(vec[4][2]+scal[4]*ρ[2]))/(4π),
        (vec[2][3] + scal[2]*ρ[3] - γ*(vec[3][3]+scal[3]*ρ[3]) + 0.5*γ^2*(vec[4][3]+scal[4]*ρ[3]))/(4π),
    ))

    c₁ = op.α
    c₂ = op.β

    α = 1 / volume(t) / volume(s) / 4
    for i in 1 : numfunctions(g)
        a = t[i]
        g = x - a
        dg = 2

        for j in 1 : numfunctions(f)
            b = s[j]

            ∫Gf = SVector(∫Gy[1]-∫G*b[1], ∫Gy[2]-∫G*b[2], ∫Gy[3]-∫G*b[3])
            ∫Gdf = 2 * ∫G
            dg∫Gf = g[1]*∫Gf[1] + g[2]*∫Gf[2] + g[3]*∫Gf[3]
            z[i,j] += ( α*c₁*dg∫Gf + α*c₂*dg*∫Gdf ) * dx

        end # next j
    end #

end

function innerintegrals!(op::MWDoubleLayer3DSng, p, g, f, t, s, z, strat::WiltonSERule, dx)

    γ = op.gamma
    x = cartesian(p)
    n = cross(s[1]-s[3],s[2]-s[3])
    n /= norm(n)
    ρ = x - ((x-s[1]) ⋅ n) * n

    scal, vec, grad = WiltonInts84.wiltonints(s[1], s[2], s[3], x, Val{1})

    # \int \nabla G_s with G_s = \nabla (1/R + 0.5*γ^2*R) / (4\pi)
    ∫∇G = (-grad[1] - 0.5*γ^2*grad[3]) / (4π)

    α = 1 / volume(t) / volume(s) / 4
    for i in 1 : numfunctions(g)
        a = t[i]
        g = (x - a)

        for j in 1 : numfunctions(f)
            b = s[j]

            z[i,j] += ( α * ( (x-b) × g ) ⋅ ∫∇G ) * dx

        end # next j
    end #

end


text:
_isone(x::T) where T = x == one(T)
#
function _ind2sub(a, i)
     i2s = CartesianIndices(a)
     Tuple(i2s[i])
end

function _sub2ind(a,i...)
      s2i = LinearIndices(a)
      s2i[i...]
end

function Base.isless(a::SVector{2}, b::SVector{2})
    if a[1] < b[1]
        return true
    elseif a[1] > b[1]
        return false
    else
        return a[2] < b[2]
    end
end

function central_range(min, max, n)
    dx = (max - min) / n
    range(min + 0.5dx, stop=max - 0.5dx, length=n)
end

function partition_work(N)
    k = Threads.nthreads()

    ls = range(1, stop=N, length=k+1)
    map(1:k) do i
        a = round(Int, ls[i])
        if i > 1
            a += 1
        end
        b = round(Int, ls[i+1])
        a:b
    end
end

function unique_points(points, tol)
    out = Vector{eltype(points)}()
    tol2 = tol * tol
    for p in points
        duplicate = false
        for v in out
            if sum(abs2, p .- v) < tol2
                duplicate = true
                break
            end
        end
        if !duplicate
            push!(out, p)
        end
    end
    out
end
"""
    rangeindex(range, x)

Returns an integer `i` such that `abs(range[i] - x) ≤ 0.5`
"""
rangeindex(range::StepRangeLen, x) = round(Int, (x - range[1]) / step(range)) + 1

"""
    realcoefficients(p::AbstractPolynomial)

Return the coefficient vector of `p` where each coefficient is converted to  a `Float64`.
"""
function realcoefficients(p::MP.AbstractPolynomial)::Vector{Float64}
    map(term -> float(MP.coefficient(term)), MP.terms(p))
end

"""
    mapcoefficients(f, p::AbstractPolynomial)

Construct a new polynomial from `p` where each `f` is applied on each coefficient.
"""
function mapcoefficients(f, p::MP.AbstractPolynomial)
    MP.polynomial(map(t -> f(MP.coefficient(t)) * MP.monomial(t), MP.terms(p)))
end

unzip(xs::Vector{SVector{2,T}}) where {T} = map(x -> x[1], xs), map(x -> x[2], xs)
unzip(xs::Vector{Tuple{T,T}}) where {T} = map(x -> x[1], xs), map(x -> x[2], xs)

logabs(x) = log(abs(x))

inIJulia() = isinteractive() && isdefined(Main, :IJulia) && isdefined(Main.IJulia, :clear_output)
plotsdefined() = isinteractive() && isdefined(Main, :Plots) && isdefined(Main.Plots, :plot)
inAtom() = isinteractive() && isdefined(Main, :Atom)

"""
    intersection_halfray_limits(origin::Point, d::Direction, xmin, xmax, ymin, ymax)

Compute the intersection of the halfray starting at `origin` in direction `d` with the
bounding window defined by `xmin`, `xmax`, `ymin` and `ymax`.
This returns also two (opposite) `Direction`s which are parallel to the window border
"""
function intersection_halfray_limits(origin::Point, d::Direction, xmin, xmax, ymin, ymax)
    x0, y0 = origin
    @assert xmin <= x0 <= xmax
    @assert ymin <= y0 <= ymax

    dx, dy = d

    if dx ≥ 0 && dy ≥ 0
        λx = (xmax - x0) / dx
        λy = (ymax - y0) / dy
    elseif dx < 0 && dy ≥ 0
        λx = (xmin - x0) / dx
        λy = (ymax - y0) / dy
    elseif dx ≥ 0 && dy < 0
        λx = (xmax - x0) / dx
        λy = (ymin - y0) / dy
    else
        λx = (xmin - x0) / dx
        λy = (ymin - y0) / dy
    end

    intersection = origin + min(λx, λy) * d
    # left or right border
    if λx < λy
        borderdir1 = SVector(0.0, 1.0)
        borderdir2 = -SVector(0.0, 1.0)
    # upper or lower border
    elseif λx > λy
        borderdir1 = SVector(1.0, 0.0)
        borderdir2 = -SVector(1.0, 0.0)
    # now we have to handle corners
    # upper right
    elseif dx ≥ 0 && dy ≥ 0
        borderdir1 = -SVector(1.0, 0.0)
        borderdir2 = -SVector(0.0, 1.0)
    # upper left
    elseif dx < 0 && dy ≥ 0
        borderdir1 = SVector(1.0, 0.0)
        borderdir2 = -SVector(0.0, 1.0)
    # lower right
    elseif dx ≥ 0 && dy < 0
        borderdir1 = -SVector(1.0, 0.0)
        borderdir2 = SVector(0.0, 1.0)
    # lower left
    else
        borderdir1 = SVector(1.0, 0.0)
        borderdir2 = SVector(0.0, 1.0)
    end
    intersection, borderdir1, borderdir2
end

"""
    fillrect!(bitmatrix, px, py)

Assigns `true` to the indices of `bitmatrix` which are inside the rectangle defined by `px` and `py`.
"""
function fillrect!(M::BitMatrix, px::SVector{4, Int}, py::SVector{4, Int})
    m, n = size(M)
    left, right = max(1, minimum(px)), min(n, maximum(px))
    @inbounds for x=left:right
        y1 = nothing
        y2 = nothing

        j = 4
        for i=1:4
            if (px[i] <= x <= px[j]) || (px[j] <= x <= px[i])
                # special case: adding the whole cut to ys
                if px[i] == px[j]
                    y1 = py[i]
                    y2 = py[j]
                else
                    y = round(Int, py[i] + (x - px[i]) / (px[j] - px[i]) * (py[j] - py[i]))
                    if y1 === nothing
                        y1 = y
                    elseif y2 === nothing && y != y1
                        y2 = y
                    end
                end
            end
            j = i
        end

        # if we only got a corner only add that value
        if y2 === nothing && 1 <= y1 <= m
            M[y1, x] = true
        elseif y1 !== nothing && y2 !== nothing
            a = max(min(y1, y2), 1)
            b = min(max(y1, y2), m)
            for k in a:b
                M[k, x] = true
            end
        end
    end
    return M
end

"""
    distance(p1, p2, p0)

Computes the distance of `p0` to the line through `p1` and `p2`.
"""
function distance(p1, p2, p0)
    x1, y1 = p1
    x2, y2 = p2
    x0, y0 = p0

    abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / sqrt((y2-y1)^2+(x2-x1)^2)
end

"""
    polyarea(p...)

Compute the area of the (simple) polygon defined by the vertices p.
*Note*: The last vertices does not to be the same as the first.
"""
polyarea(ps...) = polyarea(ps)
function polyarea(ps::NTuple{N, SVector{2, T}}) where {N, T}
    area = zero(T)
    for i=1:N-1
        @inbounds area += ps[i] × ps[i+1]
    end
    area += ps[N] × ps[1]
    0.5area
end

function iscontained(vertices, p)
    on_boundary(vertices, p) || in_interior(vertices, p)
end

function on_boundary(vertices, p)
    n = length(vertices)
    j = n
    i = 1
    while i ≤ n
        a, b = vertices[j], vertices[i]
        if p == a || p == b
            return true
        end
        if a[1] == b[1] == p[1] && 0 ≤ (p[2] - b[2]) / (a[2] - b[2]) ≤ 1
            return true
        elseif a[2] == b[2] == p[2] && 0 ≤ (p[1] - b[1]) / (a[1] - b[1]) ≤ 1
            return true
        end

        t1, t2 = (p - b) ./ (a - b)
        if abs(t1 - t2) < 1e-12 && 0 ≤ t1 ≤ 1
            return true
        end
        j = i
        i += 1
    end
    return false
end

function in_interior(vertices, p)
    n = length(vertices)
    j = n
    i = 1
    x, y = p
    c = false
    while i ≤ n
        vi_x, vi_y = vertices[i]
        vj_x, vj_y = vertices[j]
        if ((vi_y > y) != (vj_y > y)) && (x < (vj_x - vi_x) * (y - vi_y) / (vj_y - vi_y) + vi_x)
            c = !c
        end

        j = i
        i += 1
    end
    c
end

"""
    lineintersection(l1::NTuple{2, SVector{2, T}}, l2::NTuple{2, SVector{2, T}})

Compute the intersection between the lines `l1` and `l2`. `l1` and `l2`
are each defined by two points.
"""
function lineintersection(l1::NTuple{2, SVector{2, T}}, l2::NTuple{2, SVector{2, T}}) where T
    a1, b1 = l1
    a2, b2 = l2

    @inbounds x_a1_b1 = (a1[1] - b1[1])
    @inbounds y_a1_b1 = (a1[2] - b1[2])
    @inbounds x_a2_b2 = a2[1] - b2[1]
    @inbounds y_a2_b2 = (a2[2] - b2[2])

    a1_x_b1 = a1 × b1
    a2_x_b2 = a2 × b2

    denom = inv(x_a1_b1 * y_a2_b2 - y_a1_b1 * x_a2_b2)

    x = (a1_x_b1 * x_a2_b2 - x_a1_b1 * a2_x_b2) * denom
    y = (a1_x_b1 * y_a2_b2 - y_a1_b1 * a2_x_b2) * denom

    SVector{2}(x, y)
end

function equal_aspect_ratio(xmin, xmax, ymin, ymax)
    Δy = ymax - ymin
    Δx = xmax - xmin

    ratio = Δy / Δx

    if ratio ≈ 1.0
        return (xmin, xmax, ymin, ymax)
    elseif ratio > 1.0
        # y is longer than x -> increase x
        δ = 0.5 * (Δy - Δx)
        xmin -= δ
        xmax += δ

        return (xmin, xmax, ymin, ymax)
    else
        δ = 0.5 * (Δy - Δx)
        ymin -= δ
        ymax += δ

        return (xmin, xmax, ymin, ymax)
    end
end


text:
transD_GP.getchi2forall(opt, nxticks=4, fsize=6)
ax = gcf().axes;
ax[2].set_ylim(10, 30)
opt.xall[:] .= zall
transD_GP.plot_posterior(tempest, opt, burninfrac=0.5, figsize=(4,4),
    cmappdf="inferno", qp1=0.01, qp2=0.99)
step(log10.(ρ[2:end]), z[2:end], color="w", alpha=0.4)
##
gcf()


text:
"""
`init_subspec!(m::Model904)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model904)
    if subspec(m) in ["ss9", "realtime1"]
        # Use σ_r_m as SD for all anticipated shocks without normalizing
        # See measurement equation
        return
    elseif subspec(m) in ["ss10", "realtime2"]
        # Normalize s.t. sum of variances of anticipated shocks equals the
        # variance of the contemporaneous shock
        # See measurement equation
        return
    elseif subspec(m) == "realtime3"
        # Normalize s.t. sum of variances of anticipated shocks equals 1/2 the
        # variance of the contemporaneous shock
        # See measurement equation
        return
    elseif subspec(m) == "ss11"
        # Diffuse prior
        ss11!(m)
    elseif subspec(m) == "ss12"
        # Diffuse prior but keep pi_star the same as standard prior
        ss11!(m)
        ss12!(m)
    elseif subspec(m) == "ss13"
        ss13!(m)
    elseif subspec(m) == "ss14"
        return
    else
        error("This subspec is not defined.")
    end
end

function ss11!(m::Model904)
    m <= parameter(:α,      0.1687, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     Normal(0.30, 0.15),         fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's production function (also known as the capital share).",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p,   0.7467, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p,   0.2684, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label="\\iota_p")

    m <= parameter(:δ,      0.025,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    m <= parameter(:Upsilon,  1.000,  (0., 10.),     (1e-5, 0.),      ModelConstructors.Exponential(),    GammaAlt(1., 0.5),          fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\mathcal{\\Upsilon}")

    m <= parameter(:Φ,   1.4933, (1., 10.),     (1.00, 10.00),   ModelConstructors.Exponential(),    Normal(1.25, 0.36),         fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:S′′,       3.3543, (-15., 15.),   (-15., 15.),     Untransformed(),  Normal(4., 4.5),            fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.", tex_label="S\\prime\\prime")

    m <= parameter(:h,        0.4656, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),  Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="h: Consumption habit persistence.", tex_label="h")

    m <= parameter(:ppsi,     0.7614, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ppsi: Utilization costs.", tex_label="\\psi")

    m <= parameter(:ν_l,     1.0647, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    Normal(2, 2.25),            fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor
                   term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:ζ_w,   0.7922, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w,   0.5729, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w,      1.5000,                                                                               fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β,      0.7420, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.25, 0.3),        fixed=false,  scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="100(\\beta^{-1} - 1)")
   m <= parameter(:ψ1,     1.8678, (1e-5, 10.),   (1e-5, 10.00),   ModelConstructors.Exponential(),    Normal(1.5, 0.75),          fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2,     0.0715, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.15),         fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3,     0.2131, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.15),         fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")
    ## The old prior for π* in SW is GammaAlt(0.62, 0.1) with diffuse GammaAlt(0.62, 0.3) and the old prior here was GammaAlt(0.75, 0.4) so we already had a higher variance, but in keeping with the 'multiply by 3' rule, I'm multiplying by 3 anyways.
    m <= parameter(:π_star,   0.6231, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 1.2),        fixed=false,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:σ_c,   1.5073, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    Normal(1.5, 1.11),          fixed=false,
                   description="σ_c: No description available.",
                   tex_label="\\sigma_{c}")

    m <= parameter(:ρ,      0.8519, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),  Distributions.Uniform(0.0, 1.0),        fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho")

    m <= parameter(:ϵ_p,     10.000,                                                                               fixed=true,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    m <= parameter(:ϵ_w,     10.000,                                                                               fixed=true,
                   description="ϵ_w: No description available.",
                   tex_label="\\varepsilon_{w}")

    # financial frictions parameters
    m <= parameter(:Fω,      0.0300, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=true,    scaling = x -> 1 - (1-x)^0.25,
                   description="F(ω): The cumulative distribution function of ω (idiosyncratic iid shock that increases or decreases entrepreneurs' capital).",
                   tex_label="F(\\omega)")

    m <= parameter(:spr,     2.0, (0., 100.),      (1e-5, 0.),    ModelConstructors.Exponential(),   GammaAlt(2., 0.3),           fixed=false,  scaling = x -> (1 + x/100)^0.25,
                   description="spr_*: Steady-state spreads.",
                   tex_label="spr_*")

    m <= parameter(:ζ_spb, 0.05, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), Distributions.Uniform(0.0, 1.0),        fixed=false,
                   description="ζ_spb: The elasticity of the expected exess return on capital (or 'spread') with respect to leverage.",
                   tex_label="\\zeta_{spb}")

    m <= parameter(:γ_star, 0.9900, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), Distributions.Uniform(0.0, 1.0),        fixed=true,
                   description="γ_star: No description available.",
                   tex_label="\\gamma_*")

    # exogenous processes - level
    m <= parameter(:γ,      0.3085, (-5.0, 5.0),     (-5., 5.),     Untransformed(), Normal(0.4, 0.3),            fixed=false, scaling = x -> x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="100\\gamma")

    m <= parameter(:Lmean,  -1.0411, (-1000., 1000.), (-1e3, 1e3),   Untransformed(), Normal(-45, 15),              fixed=false,
                   description="Lmean: No description available.",
                   tex_label="\\bar{L}")

    m <= parameter(:g_star,    0.1800,                                                                               fixed=true,
                   description="g_star: No description available.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_g,      0.9930, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_b,      0.8100, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")

    m <= parameter(:ρ_μ,     0.7790, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_z,      0.9676, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")

    m <= parameter(:ρ_λ_f,    0.8372, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w,    0.9853, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label="\\rho_{\\lambda_w}")

    #monetary policy shock - see eqcond
    m <= parameter(:ρ_rm,     0.9000, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label="\\rho_{rm}")

    ## CHANGED THESE BETAs to UNIFORMS
    m <= parameter(:ρ_σ_w,   0.9, (1e-5, 0.9999), (1e-5, 0.9999),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ρ_σ_w: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with mean ρ_σ_w. Innovations to the process are called _spread shocks_.",
                   tex_label="\\rho_{\\sigma_\\omega}")

    m <= parameter(:ρ_μ_e,    0.7500, (1e-5, 0.99999), (1e-5, 0.99999),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=true,
                   description="ρ_μ_e: Verification costs are a fraction μ_e of the amount the bank extracts from an entrepreneur in case of bankruptcy???? This doesn't seem right because μ_e isn't a process (p12 of PDF)",
                   tex_label="\\rho_{\\mu_e}")

    m <= parameter(:ρ_γ,   0.7500, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=true,  description="ρ_γ: Autocorrelation coefficient on the γ shock.",              tex_label= "\\rho_{\\gamma}")
    m <= parameter(:ρ_π_star,   0.9900, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=true,  description="ρ_π_star: Autocorrelation coefficient on the π_star shock.",  tex_label="\\rho_{\\pi*}")

 # exogenous processes - standard deviation
    m <= parameter(:σ_g,      0.6090, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b,      0.1818, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_b: No description available.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ,     0.4601, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_z,      0.4618, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_z: No description available.",
                   tex_label="\\sigma_{z}")

    m <= parameter(:σ_λ_f,    0.1455, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good.  Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w,    0.2089, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_w: No description available.",
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m,     0.2397, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_r_m: No description available.",
                   tex_label="\\sigma_{rm}")

    m <= parameter(:σ_σ_ω,   0.2/4, (1e-7,100.),     (1e-5, 0.),    ModelConstructors.Exponential(),   RootInverseGamma(4., 0.05),  fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e,    0.0000, (1e-7,100.),     (1e-5, 0.),    ModelConstructors.Exponential(), RootInverseGamma(4., 0.05),  fixed=true,  description="σ_μ_e: Standard deviation of the μ_e shock.",         tex_label="\\sigma_{\\mu_e}")
    m <= parameter(:σ_γ,   0.0000, (1e-7,100.),     (1e-5, 0.),    ModelConstructors.Exponential(),   RootInverseGamma(4., 0.01),  fixed=true,  description="σ_γ: Standard deviation of the γ shock.",              tex_label="\\sigma_{\\gamma}")
    m <= parameter(:σ_π_star,   0.03, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(), RootInverseGamma(6., 0.03),  fixed=false, description="σ_π_star: Standard deviation of the π_star shock.", tex_label="\\sigma_{\\pi^*}")

    m <= parameter(:η_gz,       0.5632, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    m <= parameter(:η_λ_f,      0.7652, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="η_λ_f: Moving average component in the price markup shock.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w,      0.8936, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label="\\eta_{\\lambda_w}")

end

function ss12!(m::Model904)
    m <= parameter(:π_star,   0.6231, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 0.4),        fixed=false,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")
end

function ss13!(m::Model904)
    m <= parameter(:α,      0.1687, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     Normal(0.30, 0.15),         fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's production function (also known as the capital share).",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p,   0.7467, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p,   0.2684, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label="\\iota_p")

    m <= parameter(:δ,      0.025,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    m <= parameter(:Upsilon,  1.000,  (0., 10.),     (1e-5, 0.),      ModelConstructors.Exponential(),    GammaAlt(1., 0.5),          fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\mathcal{\\Upsilon}")

    m <= parameter(:Φ,   1.4933, (1., 10.),     (1.00, 10.00),   ModelConstructors.Exponential(),    Normal(1.25, 0.36),         fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:S′′,       3.3543, (-15., 15.),   (-15., 15.),     Untransformed(),  Normal(4., 4.5),            fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.", tex_label="S\\prime\\prime")

    m <= parameter(:h,        0.4656, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),  Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="h: Consumption habit persistence.", tex_label="h")

    m <= parameter(:ppsi,     0.7614, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ppsi: Utilization costs.", tex_label="\\psi")

    m <= parameter(:ν_l,     1.0647, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    Normal(2, 2.25),            fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor
                   term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:ζ_w,   0.7922, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),          fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w,   0.5729, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w,      1.5000,                                                                               fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β,      0.7420, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.25, 0.3),        fixed=false,  scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="100(\\beta^{-1} - 1)")
   m <= parameter(:ψ1,     1.8678, (1e-5, 10.),   (1e-5, 10.00),   ModelConstructors.Exponential(),    Normal(1.5, 0.75),          fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2,     0.0715, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.15),         fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3,     0.2131, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.15),         fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")
    ## The old prior for π* in SW is GammaAlt(0.62, 0.1) with diffuse GammaAlt(0.62, 0.3) and the old prior here was GammaAlt(0.75, 0.4) so we already had a higher variance, but in keeping with the 'multiply by 3' rule, I'm multiplying by 3 anyways.
    m <= parameter(:π_star,   0.6231, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 0.4),        fixed=false,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:σ_c,   1.5073, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    Normal(1.5, 1.11),          fixed=false,
                   description="σ_c: No description available.",
                   tex_label="\\sigma_{c}")

    m <= parameter(:ρ,      0.8519, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),  Distributions.Uniform(0.0, 1.0),        fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho")

    m <= parameter(:ϵ_p,     10.000,                                                                               fixed=true,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    m <= parameter(:ϵ_w,     10.000,                                                                               fixed=true,
                   description="ϵ_w: No description available.",
                   tex_label="\\varepsilon_{w}")

    # financial frictions parameters
    m <= parameter(:Fω,      0.0300, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), BetaAlt(0.03, 0.01),         fixed=true,    scaling = x -> 1 - (1-x)^0.25,
                   description="F(ω): The cumulative distribution function of ω (idiosyncratic iid shock that increases or decreases entrepreneurs' capital).",
                   tex_label="F(\\omega)")

    m <= parameter(:spr,     2.0, (0., 100.),      (1e-5, 0.),    ModelConstructors.Exponential(),   GammaAlt(2., 0.1),           fixed=false,  scaling = x -> (1 + x/100)^0.25,
                   description="spr_*: Steady-state spreads.",
                   tex_label="spr_*")

    m <= parameter(:ζ_spb, 0.05, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), BetaAlt(0.05, 0.005),        fixed=false,
                   description="ζ_spb: The elasticity of the expected exess return on capital (or 'spread') with respect to leverage.",
                   tex_label="\\zeta_{spb}")

    m <= parameter(:γ_star, 0.9900, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), BetaAlt(0.99, 0.002),        fixed=true,
                   description="γ_star: No description available.",
                   tex_label="\\gamma_*")

    # exogenous processes - level
    m <= parameter(:γ,      0.3085, (-5.0, 5.0),     (-5., 5.),     Untransformed(), Normal(0.4, 0.3),            fixed=false, scaling = x -> x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="100\\gamma")

    m <= parameter(:Lmean,  -1.0411, (-1000., 1000.), (-1e3, 1e3),   Untransformed(), Normal(-45, 15),              fixed=false,
                   description="Lmean: No description available.",
                   tex_label="\\bar{L}")

    m <= parameter(:g_star,    0.1800,                                                                               fixed=true,
                   description="g_star: No description available.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_g,      0.9930, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_b,      0.8100, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")

    m <= parameter(:ρ_μ,     0.7790, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_z,      0.9676, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")

    m <= parameter(:ρ_λ_f,    0.8372, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w,    0.9853, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label="\\rho_{\\lambda_w}")

    #monetary policy shock - see eqcond
    m <= parameter(:ρ_rm,     0.9000, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label="\\rho_{rm}")

    ## CHANGED THESE BETAs to UNIFORMS
    m <= parameter(:ρ_σ_w,   0.9, (1e-5, 0.99), (1e-5, 0.99),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="ρ_σ_w: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with mean ρ_σ_w. Innovations to the process are called _spread shocks_.",
                   tex_label="\\rho_{\\sigma_\\omega}")

    m <= parameter(:ρ_μ_e,    0.7500, (1e-5, 0.99999), (1e-5, 0.99999),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=true,
                   description="ρ_μ_e: Verification costs are a fraction μ_e of the amount the bank extracts from an entrepreneur in case of bankruptcy???? This doesn't seem right because μ_e isn't a process (p12 of PDF)",
                   tex_label="\\rho_{\\mu_e}")

    m <= parameter(:ρ_γ,   0.7500, (1e-5, 0.99999), (1e-5, 0.99),  SquareRoot(), Distributions.Uniform(0.0, 1.0),         fixed=true,  description="ρ_γ: Autocorrelation coefficient on the γ shock.",              tex_label= "\\rho_{\\gamma}")
    m <= parameter(:ρ_π_star,   0.9900, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),           fixed=true,  description="ρ_π_star: Autocorrelation coefficient on the π_star shock.",  tex_label="\\rho_{\\pi*}")

 # exogenous processes - standard deviation
    m <= parameter(:σ_g,      0.6090, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b,      0.1818, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_b: No description available.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ,     0.4601, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_z,      0.4618, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_z: No description available.",
                   tex_label="\\sigma_{z}")

    m <= parameter(:σ_λ_f,    0.1455, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good.  Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w,    0.2089, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_w: No description available.",
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m,     0.2397, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_r_m: No description available.",
                   tex_label="\\sigma_{rm}")

    m <= parameter(:σ_σ_ω,   0.2/4, (1e-7,100.),     (1e-5, 0.),    ModelConstructors.Exponential(),   RootInverseGamma(4., 0.05),  fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e,    0.0000, (1e-7,100.),     (1e-5, 0.),    ModelConstructors.Exponential(), RootInverseGamma(4., 0.05),  fixed=true,  description="σ_μ_e: Standard deviation of the μ_e shock.",         tex_label="\\sigma_{\\mu_e}")
    m <= parameter(:σ_γ,   0.0000, (1e-7,100.),     (1e-5, 0.),    ModelConstructors.Exponential(),   RootInverseGamma(4., 0.01),  fixed=true,  description="σ_γ: Standard deviation of the γ shock.",              tex_label="\\sigma_{\\gamma}")
    m <= parameter(:σ_π_star,   0.03, (1e-8, 5.),      (1e-8, 5.),    ModelConstructors.Exponential(), RootInverseGamma(6., 0.03),  fixed=false, description="σ_π_star: Standard deviation of the π_star shock.", tex_label="\\sigma_{\\pi^*}")

    m <= parameter(:η_gz,       0.5632, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    m <= parameter(:η_λ_f,      0.7652, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="η_λ_f: Moving average component in the price markup shock.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w,      0.8936, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   Distributions.Uniform(0.0, 1.0),         fixed=false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label="\\eta_{\\lambda_w}")

end


text:
module Meshing

using GeometryTypes

abstract type AbstractMeshingAlgorithm end

include("common.jl")
include("marching_tetrahedra.jl")
include("marching_cubes.jl")
include("surface_nets.jl")

export marching_cubes,
       MarchingCubes,
       MarchingTetrahedra,
       NaiveSurfaceNets

end # module


text:
@testset "README examples" begin

sdd = SinkhornDocumentDistance()

# From Kusner's WMD paper they used this illustrative example of two
# quite similar sentences (documents) which are similar semantically
# but not syntactically.
d1 = ["obama", "speaks", "media", "illinois"]
d2 = ["president", "greets", "press", "chicago"]
d12 = evaluate(sdd, d1, d2)

# To compare to something let's take a different document which we expect
# to have a larger distance to the ones above than they have to each other.
d3 = ["lawyer", "tanks", "car", "africa"]

d13 = evaluate(sdd, d1, d3)
d23 = evaluate(sdd, d2, d3)

@test d12 < d13
@test d12 < d23

wdc = WordDistanceCache()
d1 = worddistance(wdc, "robert", "programmer") #
d2 = worddistance(wdc, "robert", "astronaut")
@test d1 < d2

end