text:
"""
SensitiveFloat objects for expectations involving r_s and c_s.

Args:
vs: A vector of variational parameters

Attributes:
Each matrix has one row for each color and a column for
star / galaxy.  Row 3 is the gamma distribute baseline brightness,
and all other rows are lognormal offsets.
- E_l_a: A NUM_BANDS x NUM_SOURCE_TYPES matrix of expectations and derivatives of
  color terms.  The rows are bands, and the columns
  are star / galaxy.
- E_ll_a: A NUM_BANDS x NUM_SOURCE_TYPES matrix of expectations and derivatives of
  squared color terms.  The rows are bands, and the columns
  are star / galaxy.
"""
struct SourceBrightness{NumType <: Number}
    # [E[l|a=0], E[l]|a=1]]
    E_l_a::Matrix{SensitiveFloat{NumType}}

    # [E[l^2|a=0], E[l^2]|a=1]]
    E_ll_a::Matrix{SensitiveFloat{NumType}}
end

function SourceBrightness{NumType <: Number}(vs::Vector{NumType};
                                             calculate_gradient=true,
                                             calculate_hessian=true)
    flux_loc = vs[ids.flux_loc]
    flux_scale = vs[ids.flux_scale]
    color_mean = vs[ids.color_mean]
    color_var = vs[ids.color_var]

    # E_l_a has a row for each of the five colors and columns
    # for star / galaxy.
    E_l_a  = Matrix{SensitiveFloat{NumType}}(NUM_BANDS, NUM_SOURCE_TYPES)
    E_ll_a = Matrix{SensitiveFloat{NumType}}(NUM_BANDS, NUM_SOURCE_TYPES)

    for i = 1:NUM_SOURCE_TYPES
        for b = 1:NUM_BANDS
            E_l_a[b, i] = SensitiveFloat{NumType}(length(BrightnessParams), 1,
                                       calculate_gradient, calculate_hessian)
        end

        E_l_a[3, i].v[] = exp(flux_loc[i] + 0.5 * flux_scale[i])
        E_l_a[4, i].v[] = exp(color_mean[3, i] + .5 * color_var[3, i])
        E_l_a[5, i].v[] = exp(color_mean[4, i] + .5 * color_var[4, i])
        E_l_a[2, i].v[] = exp(-color_mean[2, i] + .5 * color_var[2, i])
        E_l_a[1, i].v[] = exp(-color_mean[1, i] + .5 * color_var[1, i])

        if calculate_gradient
            # band 3 is the reference band, relative to which the colors are
            # specified.
            # It is denoted r_s and has a lognormal expectation.
            E_l_a[3, i].d[bids.flux_loc] = E_l_a[3, i].v[]
            E_l_a[3, i].d[bids.flux_scale] = E_l_a[3, i].v[] * .5

            if calculate_hessian
                set_hess!(E_l_a[3, i], bids.flux_loc, bids.flux_loc, E_l_a[3, i].v[])
                set_hess!(E_l_a[3, i], bids.flux_loc, bids.flux_scale, E_l_a[3, i].v[] * 0.5)
                set_hess!(E_l_a[3, i], bids.flux_scale, bids.flux_scale, E_l_a[3, i].v[] * 0.25)
            end

            # The remaining indices involve c_s and have lognormal
            # expectations times E_c_3.

            # band 4 = band 3 * color 3.
            E_l_a[4, i].d[bids.color_mean[3]] = E_l_a[4, i].v[]
            E_l_a[4, i].d[bids.color_var[3]] = E_l_a[4, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[4, i], bids.color_mean[3], bids.color_mean[3], E_l_a[4, i].v[])
                set_hess!(E_l_a[4, i], bids.color_mean[3], bids.color_var[3], E_l_a[4, i].v[] * 0.5)
                set_hess!(E_l_a[4, i], bids.color_var[3], bids.color_var[3], E_l_a[4, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[4, i], E_l_a[3, i])

            # Band 5 = band 4 * color 4.
            E_l_a[5, i].d[bids.color_mean[4]] = E_l_a[5, i].v[]
            E_l_a[5, i].d[bids.color_var[4]] = E_l_a[5, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[5, i], bids.color_mean[4], bids.color_mean[4], E_l_a[5, i].v[])
                set_hess!(E_l_a[5, i], bids.color_mean[4], bids.color_var[4], E_l_a[5, i].v[] * 0.5)
                set_hess!(E_l_a[5, i], bids.color_var[4], bids.color_var[4], E_l_a[5, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[5, i], E_l_a[4, i])

            # Band 2 = band 3 * color 2.
            E_l_a[2, i].d[bids.color_mean[2]] = E_l_a[2, i].v[] * -1.
            E_l_a[2, i].d[bids.color_var[2]] = E_l_a[2, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[2, i], bids.color_mean[2], bids.color_mean[2], E_l_a[2, i].v[])
                set_hess!(E_l_a[2, i], bids.color_mean[2], bids.color_var[2], E_l_a[2, i].v[] * -0.5)
                set_hess!(E_l_a[2, i], bids.color_var[2], bids.color_var[2], E_l_a[2, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[2, i], E_l_a[3, i])

            # Band 1 = band 2 * color 1.
            E_l_a[1, i].d[bids.color_mean[1]] = E_l_a[1, i].v[] * -1.
            E_l_a[1, i].d[bids.color_var[1]] = E_l_a[1, i].v[] * .5
            if calculate_hessian
                set_hess!(E_l_a[1, i], bids.color_mean[1], bids.color_mean[1], E_l_a[1, i].v[])
                set_hess!(E_l_a[1, i], bids.color_mean[1], bids.color_var[1], E_l_a[1, i].v[] * -0.5)
                set_hess!(E_l_a[1, i], bids.color_var[1], bids.color_var[1], E_l_a[1, i].v[] * 0.25)
            end
            multiply_sfs!(E_l_a[1, i], E_l_a[2, i])
        else
            # Simply update the values if not calculating derivatives.
            E_l_a[4, i].v[] *= E_l_a[3, i].v[]
            E_l_a[5, i].v[] *= E_l_a[4, i].v[]
            E_l_a[2, i].v[] *= E_l_a[3, i].v[]
            E_l_a[1, i].v[] *= E_l_a[2, i].v[]
        end # Derivs

        ################################
        # Squared terms.

        for b = 1:NUM_BANDS
            E_ll_a[b, i] = SensitiveFloat{NumType}(length(BrightnessParams), 1,
                                           calculate_gradient, calculate_hessian)
        end

        E_ll_a[3, i].v[] = exp(2 * flux_loc[i] + 2 * flux_scale[i])
        E_ll_a[4, i].v[] = exp(2 * color_mean[3, i] + 2 * color_var[3, i])
        E_ll_a[5, i].v[] = exp(2 * color_mean[4, i] + 2 * color_var[4, i])
        E_ll_a[2, i].v[] = exp(-2 * color_mean[2, i] + 2 * color_var[2, i])
        E_ll_a[1, i].v[] = exp(-2 * color_mean[1, i] + 2 * color_var[1, i])

        if calculate_gradient
            # Band 3, the reference band.
            E_ll_a[3, i].d[bids.flux_loc] = 2 * E_ll_a[3, i].v[]
            E_ll_a[3, i].d[bids.flux_scale] = 2 * E_ll_a[3, i].v[]
            if calculate_hessian
                for hess_ids in [(bids.flux_loc, bids.flux_loc),
                                 (bids.flux_loc, bids.flux_scale),
                                 (bids.flux_scale, bids.flux_scale)]
                    set_hess!(E_ll_a[3, i], hess_ids..., 4.0 * E_ll_a[3, i].v[])
                end
            end

            # Band 4 = band 3 * color 3.
            E_ll_a[4, i].d[bids.color_mean[3]] = E_ll_a[4, i].v[] * 2.
            E_ll_a[4, i].d[bids.color_var[3]] = E_ll_a[4, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.color_mean[3], bids.color_mean[3]),
                                 (bids.color_mean[3], bids.color_var[3]),
                                 (bids.color_var[3], bids.color_var[3])]
                    set_hess!(E_ll_a[4, i], hess_ids..., E_ll_a[4, i].v[] * 4.0)
                end
            end
            multiply_sfs!(E_ll_a[4, i], E_ll_a[3, i])

            # Band 5 = band 4 * color 4.
            tmp4 = exp(2 * color_mean[4, i] + 2 * color_var[4, i])
            E_ll_a[5, i].d[bids.color_mean[4]] = E_ll_a[5, i].v[] * 2.
            E_ll_a[5, i].d[bids.color_var[4]] = E_ll_a[5, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.color_mean[4], bids.color_mean[4]),
                                 (bids.color_mean[4], bids.color_var[4]),
                                 (bids.color_var[4], bids.color_var[4])]
                    set_hess!(E_ll_a[5, i], hess_ids..., E_ll_a[5, i].v[] * 4.0)
                end
            end
            multiply_sfs!(E_ll_a[5, i], E_ll_a[4, i])

            # Band 2 = band 3 * color 2
            tmp2 = exp(-2 * color_mean[2, i] + 2 * color_var[2, i])
            E_ll_a[2, i].d[bids.color_mean[2]] = E_ll_a[2, i].v[] * -2.
            E_ll_a[2, i].d[bids.color_var[2]] = E_ll_a[2, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.color_mean[2], bids.color_mean[2]),
                                 (bids.color_var[2], bids.color_var[2])]
                    set_hess!(E_ll_a[2, i], hess_ids..., E_ll_a[2, i].v[] * 4.0)
                end
                set_hess!(E_ll_a[2, i], bids.color_mean[2], bids.color_var[2],
                          E_ll_a[2, i].v[] * -4.0)
            end
            multiply_sfs!(E_ll_a[2, i], E_ll_a[3, i])

            # Band 1 = band 2 * color 1
            E_ll_a[1, i].d[bids.color_mean[1]] = E_ll_a[1, i].v[] * -2.
            E_ll_a[1, i].d[bids.color_var[1]] = E_ll_a[1, i].v[] * 2.
            if calculate_hessian
                for hess_ids in [(bids.color_mean[1], bids.color_mean[1]),
                                 (bids.color_var[1], bids.color_var[1])]
                    set_hess!(E_ll_a[1, i], hess_ids..., E_ll_a[1, i].v[] * 4.0)
                end
                set_hess!(E_ll_a[1, i], bids.color_mean[1], bids.color_var[1],
                          E_ll_a[1, i].v[] * -4.0)
            end
            multiply_sfs!(E_ll_a[1, i], E_ll_a[2, i])
        else
            # Simply update the values if not calculating derivatives.
            E_ll_a[4, i].v[] *= E_ll_a[3, i].v[]
            E_ll_a[5, i].v[] *= E_ll_a[4, i].v[]
            E_ll_a[2, i].v[] *= E_ll_a[3, i].v[]
            E_ll_a[1, i].v[] *= E_ll_a[2, i].v[]
        end # calculate_gradient
    end

    SourceBrightness(E_l_a, E_ll_a)
end

"""
Load the source brightnesses for these model params.  Each SourceBrightness
object has information for all bands and object types.

Returns:
  - An array of SourceBrightness objects for each object in 1:ea.S.  Only
    sources in ea.active_sources will have derivative information.
"""
function load_source_brightnesses{NumType <: Number}(
                    ea::ElboArgs,
                    vp::VariationalParams{NumType};
                    calculate_gradient::Bool=true,
                    calculate_hessian::Bool=true)
    sbs = Vector{SourceBrightness{NumType}}(ea.S)

    for s in 1:ea.S
        this_deriv = (s in ea.active_sources) && calculate_gradient
        this_hess = (s in ea.active_sources) && calculate_hessian
        sbs[s] = SourceBrightness(vp[s];
                                  calculate_gradient=this_deriv,
                                  calculate_hessian=this_hess)
    end

    sbs
end


text:
using Dates

@testset "util.jl" begin
    p1 = DDR2import.util.Point(1.0, 2.0)
    @test p1.lat == 1.0
    @test p1.lon == 2.0

    @test DDR2import.util.extract_lat("S123456") ≈ -12.58222222 atol = 0.0001
    @test DDR2import.util.extract_lat("N123456.0") ≈ 12.58222222 atol = 0.0001
    @test DDR2import.util.extract_lat("123456N") ≈ 12.58222222 atol = 0.0001
    @test DDR2import.util.extract_lat("S1234") ≈ -12.56666666 atol = 0.0001
    @test DDR2import.util.extract_lat("N12") ≈ 12.0 atol = 0.0001
    @test DDR2import.util.extract_lat("1234S") ≈ -12.56666666 atol = 0.0001
    @test DDR2import.util.extract_lat("12N") ≈ 12.0 atol = 0.0001
    @test DDR2import.util.extract_lat("") === NaN

    @test DDR2import.util.extract_lon("W1123456") ≈ -112.58222222 atol = 0.0001
    @test DDR2import.util.extract_lon("E1123456.0") ≈ 112.58222222 atol = 0.0001
    @test DDR2import.util.extract_lon("0123456E") ≈ 12.58222222 atol = 0.0001
    @test DDR2import.util.extract_lon("W01234") ≈ -12.56666666 atol = 0.0001
    @test DDR2import.util.extract_lon("E112") ≈ 112.0 atol = 0.0001
    @test DDR2import.util.extract_lon("11234W") ≈ -112.56666666 atol = 0.0001
    @test DDR2import.util.extract_lon("012E") ≈ 12.0 atol = 0.0001
    @test DDR2import.util.extract_lon("E1") === NaN

    # @test DDR2import.util.latlon("S123456", "012E").lat ≈ -12.58222222 atol = 0.0001
    # @test DDR2import.util.latlon("S123456", "012E").lon ≈ 12.0 atol = 0.0001
    # @test DDR2import.util.latlon("S123456", "E1") === NaN
    #
    # @test DDR2import.util.latlon("S123456E1123456.0").lat ≈ -12.58222222 atol = 0.0001
    # @test DDR2import.util.latlon("S123456E1123456.0").lon ≈ 112.58222222 atol = 0.0001
    # @test DDR2import.util.latlon("S123456W1123456.0").lon ≈ -112.58222222 atol = 0.0001

    @test DDR2import.util.latlon("S123456", "012E").ϕ ≈ -12.58222222 atol = 0.0001
    @test DDR2import.util.latlon("S123456", "012E").λ ≈ 12.0 atol = 0.0001
    @test DDR2import.util.latlon("S123456", "E1") === NaN

    @test DDR2import.util.latlon("123456S1123456.0E").ϕ ≈ -12.58222222 atol = 0.0001
    @test DDR2import.util.latlon("123456S1123456.0E").λ ≈ 112.58222222 atol = 0.0001
    @test DDR2import.util.latlon("123456S1123456.0W").λ ≈ -112.58222222 atol = 0.0001

    yymmdd = DateFormat("YYmmdd")
    testdate = DDR2import.util.format_date("200110", yymmdd, addyear=Year(2000))
    @test Dates.year(testdate) == 2020
    @test Dates.month(testdate) == 1
    @test Dates.day(testdate) == 10

    hhmmss = DateFormat("HHMMSS")
    testtime = DDR2import.util.format_time("102030", hhmmss)
    @test Dates.hour(testtime) == 10
    @test Dates.minute(testtime) == 20
    @test Dates.second(testtime) == 30
    hhmm = DateFormat("HHMM")
    testtime2 = DDR2import.util.format_time("2400", hhmm)
    @test Dates.hour(testtime2) == 23
    @test Dates.minute(testtime2) == 59
    @test Dates.second(testtime2) == 59

    yyyymmddhhmmss = DateFormat("YYYYmmddHHMMSS")
    testdatetime = DDR2import.util.format_datetime("20200110102030", yyyymmddhhmmss)
    @test Dates.year(testdatetime) == 2020
    @test Dates.month(testdatetime) == 1
    @test Dates.day(testdatetime) == 10
    @test Dates.hour(testdatetime) == 10
    @test Dates.minute(testdatetime) == 20
    @test Dates.second(testdatetime) == 30
end

ext: jl


text:
function test10()
  sim = LTspiceSimulation("test10.asc",executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)
  @test sim["m1"] == 0.0
  @test isnan(sim["bad_meas"])
  show(IOBuffer(),sim)
end
test10()


text:
module SokobanDirectedModule

import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase
import ..SokobanUndirectedModule as SUM

#####
##### game logic
#####

const NUM_OBJECTS = SUM.NUM_OBJECTS
const AGENT = SUM.AGENT
const WALL = SUM.WALL
const BOX = SUM.BOX
const TARGET = SUM.TARGET
const NUM_ACTIONS = 4

mutable struct SokobanDirected{R, RNG} <: GW.AbstractGridWorld
    env::SUM.SokobanUndirected{R, RNG}
    agent_direction::Int
end

function SokobanDirected(; R = Float32, file = joinpath(dirname(functionloc(GW.eval)[1]), "envs/sokoban/boxoban-levels/medium/train/000.txt"), rng = Random.GLOBAL_RNG)
    env = SUM.SokobanUndirected(R = R, file = file, rng = rng)
    agent_direction = rand(rng, 0:GW.NUM_DIRECTIONS-1)
    env = SokobanDirected(env, agent_direction)
    GW.reset!(env)
    return env
end

function GW.reset!(env::SokobanDirected)
    GW.reset!(env.env)
    env.agent_direction = rand(env.env.rng, 0:GW.NUM_DIRECTIONS-1)
    return nothing
end

function GW.act!(env::SokobanDirected, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    inner_env = env.env
    tile_map = inner_env.tile_map
    box_positions = inner_env.box_positions
    box_reward = inner_env.box_reward

    r1 = sum(pos -> tile_map[TARGET, pos], box_positions) * box_reward
    agent_position = inner_env.agent_position
    agent_direction = env.agent_direction

    if action in Base.OneTo(2)
        if action == 1
            dest = GW.move_forward(agent_position, agent_direction)
            beyond_dest = GW.move_forward(dest, agent_direction)
        else
            dest = GW.move_backward(agent_position, agent_direction)
            beyond_dest = GW.move_backward(dest, agent_direction)
        end

        if !tile_map[WALL, dest]
            if !tile_map[BOX, dest]
                tile_map[AGENT, agent_position] = false
                inner_env.agent_position = dest
                tile_map[AGENT, dest] = true
            else
                if !tile_map[BOX, beyond_dest] && !tile_map[WALL, beyond_dest]
                    tile_map[BOX, dest] = false
                    tile_map[BOX, beyond_dest] = true

                    box_idx = findfirst(==(dest), box_positions)
                    box_positions[box_idx] = beyond_dest
                    tile_map[AGENT, agent_position] = false
                    inner_env.agent_position = dest
                    tile_map[AGENT, dest] = true
                end
            end
        end
    elseif action == 3
        env.agent_direction = GW.turn_left(agent_direction)
    else
        env.agent_direction = GW.turn_right(agent_direction)
    end

    inner_env.done = all(pos -> tile_map[TARGET, pos], box_positions)
    r2 = sum(pos -> tile_map[TARGET, pos], box_positions) * box_reward
    inner_env.reward = r2 - r1

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::SokobanDirected) = GW.get_height(env.env)
GW.get_width(env::SokobanDirected) = GW.get_width(env.env)

GW.get_action_names(env::SokobanDirected) = (:MOVE_FORWARD, :MOVE_BACKWARD, :TURN_LEFT, :TURN_RIGHT)
GW.get_object_names(env::SokobanDirected) = GW.get_object_names(env.env)

function GW.get_pretty_tile_map(env::SokobanDirected, position::CartesianIndex{2})
    characters = ('☻', '█', '▒', '✖', '→', '↑', '←', '↓', '⋅')

    object = findfirst(@view env.env.tile_map[:, position])
    if isnothing(object)
        return characters[end]
    elseif object == AGENT
        return characters[NUM_OBJECTS + 1 + env.agent_direction]
    else
        return characters[object]
    end
end

function GW.get_pretty_sub_tile_map(env::SokobanDirected, window_size, position::CartesianIndex{2})
    tile_map = env.env.tile_map
    agent_position = env.env.agent_position
    agent_direction = env.agent_direction

    characters = ('☻', '█', '▒', '✖', '→', '↑', '←', '↓', '⋅')

    sub_tile_map = GW.get_sub_tile_map(tile_map, agent_position, window_size, agent_direction)

    object = findfirst(@view sub_tile_map[:, position])
    if isnothing(object)
        return characters[end]
    elseif object == AGENT
        return '↓'
    else
        return characters[object]
    end
end

function Base.show(io::IO, ::MIME"text/plain", env::SokobanDirected)
    str = "tile_map:\n"
    str = str * GW.get_pretty_tile_map(env)
    str = str * "\nsub_tile_map:\n"
    str = str * GW.get_pretty_sub_tile_map(env, GW.get_window_size(env))
    str = str * "\nreward: $(env.env.reward)"
    str = str * "\ndone: $(env.env.done)"
    str = str * "\nlevel_number: $(env.env.level_number)"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::SokobanDirected) = ('w', 's', 'a', 'd')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: SokobanDirected} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: SokobanDirected} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: SokobanDirected} = (env.env.env.tile_map, env.env.agent_direction)

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: SokobanDirected} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: SokobanDirected} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: SokobanDirected} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: SokobanDirected} = env.env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: SokobanDirected} = env.env.env.done

end # module


text:
struct CuVecCodeGen
    eltyp::Type
end

function generate_code(codegen::CuVecCodeGen, g::ExGraph, nd::ExNode)
    # nd = cast_const_type(nd, codegen.eltyp)
    ex = to_expr_kw(nd)
    ex = rewrite_all(ex, CUDA_NATIVE_RULES; phs=[:x, :y, :z, :n])
    ex = subs(ex, Dict(:__FLOAT_TYPE__ => codegen.eltyp))
    return ex
end

function generate_code(codegen::CuVecCodeGen, g::ExGraph)
    g = eliminate_common(g)
    # g = cast_const_type(g, codegen.eltyp)
    ex = to_expr_kw(g)
    ex = rewrite_all(ex, CUDA_NATIVE_RULES; phs=[:x, :y, :z, :n])
    ex = subs(ex, Dict(:__FLOAT_TYPE__ => codegen.eltyp))
    return ex
end

eval_codegen(codegen::CuVecCodeGen) = codegen


text:
export
    m2eul,
    matchi,
    matchw,
    maxd,
    maxi,
    mequ,
    mequg,
    mind,
    mini,
    mtxm,
    mtxmg,
    mtxv,
    mtxvg,
    mxm,
    mxmg,
    mxmt,
    mxmtg,
    mxv,
    mxvg,
    m2q

"""
    m2eul(r, axis3, axis2, axis1)

Factor a rotation matrix as a product of three rotations about specified coordinate axes.

### Arguments ###

- `r`: A rotation matrix to be factored
- `axis3`: Number of the third rotation axis
- `axis2`: Number of the second rotation axis
- `axis1`: Number of the first rotation axis

### Output ###

A tuple consisting of the third, second, and first Euler angles in radians.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/m2eul_c.html)
"""
function m2eul(r, axis3, axis2, axis1)
    @checkdims 3 3 r
    angle3 = Ref{SpiceDouble}()
    angle2 = Ref{SpiceDouble}()
    angle1 = Ref{SpiceDouble}()
    ccall((:m2eul_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceInt,
           Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
           permutedims(r), axis3, axis2, axis1, angle3, angle2, angle1)
    handleerror()
    angle3[], angle2[], angle1[]
end

"""
    matchi(string, templ, wstr, wchar)

Determine whether a string is matched by a template containing wild cards.
The pattern comparison is case-insensitive.

### Arguments ###

- `string`: String to be tested
- `templ`: Template (with wild cards) to test against string
- `wstr`: Wild string token
- `wchr`: Wild character token

### Output ###

Returns `true` if the string matches.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/matchi_c.html)
"""
function matchi(string, templ, wstr, wchar)
    res = Bool(ccall((:matchi_c, libcspice), SpiceBoolean, (Cstring, Cstring, SpiceChar, SpiceChar),
          string, templ, first(wstr), first(wchar)))
    handleerror()
    res
end

"""
    matchw(string, templ, wstr, wchar)

Determine whether a string is matched by a template containing wild cards.

### Arguments ###

- `string`: String to be tested
- `templ`: Template (with wild cards) to test against string
- `wstr`: Wild string token
- `wchr`: Wild character token

### Output ###

Returns `true` if the string matches.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/matchw_c.html)
"""
function matchw(string, templ, wstr, wchar)
    res = Bool(ccall((:matchw_c, libcspice), SpiceBoolean, (Cstring, Cstring, SpiceChar, SpiceChar),
          string, templ, first(wstr), first(wchar)))
    handleerror()
    res
end

function _maxd(args...)
    n = length(args)
    @eval ccall((:maxd_c, libcspice), SpiceDouble, (SpiceInt, SpiceDouble...), $n, $(args...))
end

"""
    maxd(args...)

!!! warning "Deprecated"
    Use `max(args...)` instead.
"""
maxd

@deprecate maxd max

function _maxi(args...)
    n = length(args)
    @eval ccall((:maxi_c, libcspice), SpiceInt, (SpiceInt, SpiceInt...), $n, $(args...))
end

"""
    maxi(args...)

!!! warning "Deprecated"
    Use `max(args...)` instead.
"""
maxi

@deprecate maxi max

function _mequ(a, b)
    ccall((:mequ_c, libcspice), Cvoid, (Ref{SpiceDouble}, Ref{SpiceDouble}), a, b)
end

"""
    mequ(m1, mout)

!!! warning "Deprecated"
    Use `mout .= m1` instead.
"""
mequ

@deprecate mequ(m1, mout) mout .= m1

function _mequg(a, b)
    m, n = size(a)
    ccall((:mequg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, SpiceInt, Ref{SpiceDouble}), a, m, n, b)
end

"""
    mequg(m1, mout)

!!! warning "Deprecated"
    Use `mout .= m1` instead.
"""
mequg

@deprecate mequg(m1, mout) mout .= m1

function _mind(args...)
    n = length(args)
    @eval ccall((:mind_c, libcspice), SpiceDouble, (SpiceInt, SpiceDouble...), $n, $(args...))
end

"""
    mind(args...)

!!! warning "Deprecated"
    Use `min(args...)` instead.
"""
mind

@deprecate mind min

function _mini(args...)
    n = length(args)
    @eval ccall((:mini_c, libcspice), SpiceInt, (SpiceInt, SpiceInt...), $n, $(args...))
end

"""
    mini(args...)

!!! warning "Deprecated"
    Use `min(args...)` instead.
"""
mini

@deprecate mini min

function _mtxm(m1, m2)
    mout = Array{SpiceDouble}(undef, 3, 3)
    ccall((:mtxm_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          permutedims(m1), permutedims(m2), mout)
    permutedims(mout)
end

"""
    mtxm(m1, m2)

!!! warning "Deprecated"
    Use `m1' * m2` instead.
"""
mtxm

@deprecate mtxm(m1, m2) m1' * m2

function _mtxmg(m1, m2)
    l1, n = size(m1)
    l2, k = size(m2)
    @assert l1 == l2
    mout = Array{SpiceDouble}(undef, k, n)
    ccall((:mtxmg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          permutedims(m1), permutedims(m2), n, l1, k, mout)
    permutedims(mout)
end

"""
    mtxmg(m1, m2)

!!! warning "Deprecated"
    Use `m1' * m2` instead.
"""
mtxmg

@deprecate mtxmg(m1, m2) m1' * m2

function _mtxv(m1, v2)
    vout = Array{Float64}(undef, 3)
    ccall((:mtxv_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          permutedims(m1), v2, vout)
    handleerror()
    vout
end

"""
    mtxv(m1,v2)

!!! warning "Deprecated"
    Use `m1' * v2` instead.
"""
mtxv

@deprecate mtxv(m1, v2) m1' * v2

function _mtxvg(m1, v2)
    lm1, lm2 = size(m1)
    lv = length(v2)
    @assert lm1 == lv
    vout = Array{Float64}(undef, lm2)
    ccall((:mtxvg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          permutedims(m1), v2, lm2, lm1, vout)
    handleerror()
    vout
end

"""
    mtxvg(m1,v2)

!!! warning "Deprecated"
    Use `m1' * v2` instead.
"""
mtxvg

@deprecate mtxvg(m1, v2) m1' * v2

function _mxm(m1, m2)
    mout = Array{SpiceDouble}(undef, 3, 3)
    ccall((:mxm_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          permutedims(m1), permutedims(m2), mout)
    permutedims(mout)
end

"""
    mxm(m1, m2)

!!! warning "Deprecated"
    Use `m1 * m2` instead.
"""
mxm

@deprecate mxm(m1, m2) m1 * m2

function _mxmg(m1, m2)
    n, l1 = size(m1)
    l2, k = size(m2)
    @assert l1 == l2
    mout = Array{SpiceDouble}(undef, k, n)
    ccall((:mxmg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          permutedims(m1), permutedims(m2), n, l1, k, mout)
    permutedims(mout)
end

"""
    mxmg(m1, m2)

!!! warning "Deprecated"
    Use `m1 * m2` instead.
"""
mxmg

@deprecate mxmg(m1, m2) m1 * m2

function _mxmt(m1, m2)
    mout = Array{SpiceDouble}(undef, 3, 3)
    ccall((:mxmt_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          permutedims(m1), permutedims(m2), mout)
    permutedims(mout)
end

"""
    mxmt(m1, m2)

!!! warning "Deprecated"
    Use `m1 * m2'` instead.
"""
mxmt

@deprecate mxmt(m1, m2) m1 * m2'

function _mxmtg(m1, m2)
    n, l1 = size(m1)
    k, l2 = size(m2)
    @assert l1 == l2
    mout = Array{SpiceDouble}(undef, k, n)
    ccall((:mxmtg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          permutedims(m1), permutedims(m2), n, l1, k, mout)
    permutedims(mout)
end

"""
    mxmtg(m1, m2)

!!! warning "Deprecated"
    Use `m1 * m2'` instead.
"""
mxmtg

@deprecate mxmtg(m1, m2) m1 * m2'

function _mxv(m1, v2)
    vout = Array{Float64}(undef, 3)
    ccall((:mxv_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          permutedims(m1), v2, vout)
    handleerror()
    vout
end

"""
    mxv(m1,v2)

!!! warning "Deprecated"
    Use `m1 * v2` instead.
"""
mxv

@deprecate mxv(m1, v2) m1 * v2

function _mxvg(m1, v2)
    lm1, lm2 = size(m1)
    lv = length(v2)
    @assert lm2 == lv
    vout = Array{Float64}(undef, lm1)
    ccall((:mxvg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          permutedims(m1), v2, lm1, lm2, vout)
    handleerror()
    vout
end

"""
    mxvg(m1,v2)

!!! warning "Deprecated"
    Use `m1 * v2` instead.
"""
mxvg

@deprecate mxvg(m1, v2) m1 * v2

"""
    m2q(r)

Find a unit quaternion corresponding to a specified rotation matrix.

### Arguments ###

- `r`: A rotation matrix

### Output ###

A unit quaternion representing `r'

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/m2q_c.html)
"""
function m2q(r)
    @checkdims 3 3 r
    q = Array{SpiceDouble}(undef, 4)
    ccall((:m2q_c, libcspice), Cvoid, (Ref{Float64}, Ref{Float64}), permutedims(r), q)
    handleerror()
    q
end


text:
###############################################################################
#
#   nf_elem.jl : Antic number fields
#
###############################################################################

export AnticNumberField, nf_elem, norm, trace, CyclotomicField,
       MaximalRealSubfield, add!, sub!, mul!, signature, sqr_classical

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent_type(::Type{nf_elem}) = AnticNumberField

doc"""
    parent(a::nf_elem)
> Return the parent of the given number field element.
"""
parent(a::nf_elem) = a.parent

elem_type(::Type{AnticNumberField}) = nf_elem

doc"""
    base_ring(a::AnticNumberField)
> Returns `Union{}` since a number field doesn't depend on any ring.
"""
base_ring(a::AnticNumberField) = Union{}

doc"""
    base_ring(a::nf_elem)
> Returns `Union{}` since a number field doesn't depend on any ring.
"""
base_ring(a::nf_elem) = Union{}

isdomain_type(::Type{nf_elem}) = true

doc"""
    var(a::AnticNumberField)
> Returns the identifier (as a symbol, not a string), that is used for printing
> the generator of the given number field.
"""
var(a::AnticNumberField) = a.S

function check_parent(a::nf_elem, b::nf_elem)
   a.parent != b.parent && error("Incompatible number field elements")
end

show_minus_one(::Type{nf_elem}) = false

###############################################################################
#
#   Basic manipulation
#
###############################################################################

const hash_seed = UInt==UInt64 ? 0xc2a44fbe466a1827 : 0xc2a44fb

function hash(a::nf_elem, h::UInt)
   global hash_seed
   b = hash_seed
   d = denominator(a)
   b = hash(d, b)
   for i in 1:degree(parent(a)) + 1
         num_coeff!(d, a, i)
         b = xor(b, xor(hash(d, h), h))
         b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
   end
   return b
end

doc"""
    coeff(x::nf_elem, n::Int)
> Return the $n$-th coefficient of the polynomial representation of the given
> number field element. Coefficients are numbered from $0$, starting with the
> constant coefficient.
"""
function coeff(x::nf_elem, n::Int)
   n < 0 && throw(DomainError())
   z = fmpq()
   ccall((:nf_elem_get_coeff_fmpq, :libflint), Void,
     (Ref{fmpq}, Ref{nf_elem}, Int, Ref{AnticNumberField}), z, x, n, parent(x))
   return z
end

function num_coeff!(z::fmpz, x::nf_elem, n::Int)
   n < 0 && throw(DomainError())
   ccall((:nf_elem_get_coeff_fmpz, :libflint), Void,
     (Ref{fmpz}, Ref{nf_elem}, Int, Ref{AnticNumberField}), z, x, n, parent(x))
   return z
end

doc"""
    gen(a::AnticNumberField)
> Return the generator of the given number field.
"""
function gen(a::AnticNumberField)
   r = nf_elem(a)
   ccall((:nf_elem_gen, :libflint), Void,
         (Ref{nf_elem}, Ref{AnticNumberField}), r, a)
   return r
end

doc"""
    one(a::AnticNumberField)
> Return the multiplicative identity, i.e. one, in the given number field.
"""
function one(a::AnticNumberField)
   r = nf_elem(a)
   ccall((:nf_elem_one, :libflint), Void,
         (Ref{nf_elem}, Ref{AnticNumberField}), r, a)
   return r
end

doc"""
    zero(a::AnticNumberField)
> Return the multiplicative identity, i.e. one, in the given number field.
"""
function zero(a::AnticNumberField)
   r = nf_elem(a)
   ccall((:nf_elem_zero, :libflint), Void,
         (Ref{nf_elem}, Ref{AnticNumberField}), r, a)
   return r
end

doc"""
    isgen(a::nf_elem)
> Return `true` if the given number field element is the generator of the
> number field, otherwise return `false`.
"""
function isgen(a::nf_elem)
   return ccall((:nf_elem_is_gen, :libflint), Bool,
                (Ref{nf_elem}, Ref{AnticNumberField}), a, a.parent)
end

doc"""
    isone(a::nf_elem)
> Return `true` if the given number field element is the multiplicative
> identity of the number field, i.e. one, otherwise return `false`.
"""
function isone(a::nf_elem)
   return ccall((:nf_elem_is_one, :libflint), Bool,
                (Ref{nf_elem}, Ref{AnticNumberField}), a, a.parent)
end

doc"""
    iszero(a::nf_elem)
> Return `true` if the given number field element is the additive
> identity of the number field, i.e. zero, otherwise return `false`.
"""
function iszero(a::nf_elem)
   return ccall((:nf_elem_is_zero, :libflint), Bool,
                (Ref{nf_elem}, Ref{AnticNumberField}), a, a.parent)
end

doc"""
    isunit(a::nf_elem)
> Return `true` if the given number field element is invertible, i.e. nonzero,
> otherwise return `false`.
"""
isunit(a::nf_elem) = a != 0

doc"""
    denominator(a::nf_elem)
> Return the denominator of the polynomial representation of the given number
> field element.
"""
function denominator(a::nf_elem)
   z = fmpz()
   ccall((:nf_elem_get_den, :libflint), Void,
         (Ref{fmpz}, Ref{nf_elem}, Ref{AnticNumberField}),
         z, a, a.parent)
   return z
end

function elem_from_mat_row(a::AnticNumberField, b::fmpz_mat, i::Int, d::fmpz)
   Generic._checkbounds(rows(b), i) || throw(BoundsError())
   cols(b) == degree(a) || error("Wrong number of columns")
   z = a()
   ccall((:nf_elem_set_fmpz_mat_row, :libflint), Void,
        (Ref{nf_elem}, Ref{fmpz_mat}, Int, Ref{fmpz}, Ref{AnticNumberField}),
        z, b, i - 1, d, a)
   return z
end

function elem_to_mat_row!(a::fmpz_mat, i::Int, d::fmpz, b::nf_elem)
   ccall((:nf_elem_get_fmpz_mat_row, :libflint), Void,
         (Ref{fmpz_mat}, Int, Ref{fmpz}, Ref{nf_elem}, Ref{AnticNumberField}),
         a, i - 1, d, b, b.parent)
   nothing
 end

doc"""
    degree(a::AnticNumberField)
> Return the degree of the given number field, i.e. the degree of its
> defining polynomial.
"""
degree(a::AnticNumberField) = degree(a.pol)

doc"""
    signature(a::AnticNumberField)
> Return the signature of the given number field, i.e. a tuple $r, s$
> consisting of $r$, the number of real embeddings and $s$, half the number of
> complex embeddings.
"""
signature(a::AnticNumberField) = signature(a.pol)

function deepcopy_internal(d::nf_elem, dict::ObjectIdDict)
   z = nf_elem(parent(d), d)
   return z
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::AnticNumberField)
   print(io, "Number field over Rational Field")
   print(io, " with defining polynomial ", a.pol)
end

function show(io::IO, x::nf_elem)
   cstr = ccall((:nf_elem_get_str_pretty, :libflint), Ref{UInt8},
                (Ref{nf_elem}, Ref{UInt8}, Ref{AnticNumberField}),
                 x, string(var(parent(x))), parent(x))
   s = unsafe_string(cstr)
   ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)

   s = replace(s, "/", "//")

   print(io, s)

end

needs_parentheses(::Nemo.nf_elem) = true

isnegative(::nf_elem) = false

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(a::nf_elem)
   r = a.parent()
   ccall((:nf_elem_neg, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, a.parent)
   return r
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::nf_elem, b::nf_elem)
   check_parent(a, b)
   r = a.parent()
   ccall((:nf_elem_add, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function -(a::nf_elem, b::nf_elem)
   check_parent(a, b)
   r = a.parent()
   ccall((:nf_elem_sub, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function *(a::nf_elem, b::nf_elem)
   check_parent(a, b)
   r = a.parent()
   ccall((:nf_elem_mul, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function +(a::nf_elem, b::Int)
   r = a.parent()
   ccall((:nf_elem_add_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function +(a::nf_elem, b::fmpz)
   r = a.parent()
   ccall((:nf_elem_add_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function +(a::nf_elem, b::fmpq)
   r = a.parent()
   ccall((:nf_elem_add_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function -(a::nf_elem, b::Int)
   r = a.parent()
   ccall((:nf_elem_sub_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function -(a::nf_elem, b::fmpz)
   r = a.parent()
   ccall((:nf_elem_sub_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function -(a::nf_elem, b::fmpq)
   r = a.parent()
   ccall((:nf_elem_sub_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function -(a::Int, b::nf_elem)
   r = b.parent()
   ccall((:nf_elem_si_sub, :libflint), Void,
         (Ref{nf_elem}, Int, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, b.parent)
   return r
end

function -(a::fmpz, b::nf_elem)
   r = b.parent()
   ccall((:nf_elem_fmpz_sub, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpz}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, b.parent)
   return r
end

function -(a::fmpq, b::nf_elem)
   r = b.parent()
   ccall((:nf_elem_fmpq_sub, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpq}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, b.parent)
   return r
end

+(a::nf_elem, b::Integer) = a + fmpz(b)

-(a::nf_elem, b::Integer) = a - fmpz(b)

-(a::Integer, b::nf_elem) = fmpz(a) - b

+(a::Integer, b::nf_elem) = b + a

+(a::fmpq, b::nf_elem) = b + a

+(a::Rational, b::nf_elem) = fmpq(a) + b

+(a::nf_elem, b::Rational) = b + a

-(a::Rational, b::nf_elem) = fmpq(a) - b

-(a::nf_elem, b::Rational) = a - fmpq(b)

function *(a::nf_elem, b::Int)
   r = a.parent()
   ccall((:nf_elem_scalar_mul_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function *(a::nf_elem, b::fmpz)
   r = a.parent()
   ccall((:nf_elem_scalar_mul_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function *(a::nf_elem, b::fmpq)
   r = a.parent()
   ccall((:nf_elem_scalar_mul_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function *(a::Rational, b::nf_elem)
  return fmpq(a) * b
end

*(a::nf_elem, b::Rational) = b*a

*(a::Integer, b::nf_elem) = b * a

*(a::fmpz, b::nf_elem) = b * a

*(a::fmpq, b::nf_elem) = b * a

//(a::nf_elem, b::Int) = divexact(a, b)

//(a::nf_elem, b::fmpz) = divexact(a, b)

//(a::nf_elem, b::Integer) = a//fmpz(b)

//(a::nf_elem, b::fmpq) = divexact(a, b)

//(a::Integer, b::nf_elem) = divexact(a, b)

//(a::fmpz, b::nf_elem) = divexact(a, b)

//(a::fmpq, b::nf_elem) = divexact(a, b)

//(a::Rational, b::nf_elem) = divexact(fmpq(a), b)

//(a::nf_elem, b::Rational) = divexact(a, fmpq(b))

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::nf_elem, n::Int)
   r = a.parent()
   ccall((:nf_elem_pow, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         r, a, abs(n), a.parent)
   if n < 0
      r = inv(r)
   end
   return r
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::nf_elem, b::nf_elem)
   check_parent(a, b)
   return ccall((:nf_elem_equal, :libflint), Bool,
           (Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}), a, b, a.parent)
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

==(a::nf_elem, b::Integer) = a == parent(a)(b)

==(a::nf_elem, b::fmpz) = a == parent(a)(b)

==(a::nf_elem, b::fmpq) = a == parent(a)(b)

==(a::Integer, b::nf_elem) = parent(b)(a) == b

==(a::fmpz, b::nf_elem) = parent(b)(a) == b

==(a::fmpq, b::nf_elem) = parent(b)(a) == b

###############################################################################
#
#   Inversion
#
###############################################################################

doc"""
    inv(a::nf_elem)
> Return $a^{-1}$. Requires $a \neq 0$.
"""
function inv(a::nf_elem)
   iszero(a) && throw(DivideError())
   r = a.parent()
   ccall((:nf_elem_inv, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, a.parent)
   return r
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(a::nf_elem, b::nf_elem)
   iszero(b) && throw(DivideError())
   check_parent(a, b)
   r = a.parent()
   ccall((:nf_elem_div, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(a::nf_elem, b::Int)
   b == 0 && throw(DivideError())
   r = a.parent()
   ccall((:nf_elem_scalar_div_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

function divexact(a::nf_elem, b::fmpz)
   iszero(b) && throw(DivideError())
   r = a.parent()
   ccall((:nf_elem_scalar_div_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

divexact(a::nf_elem, b::Integer) = divexact(a, fmpz(b))

function divexact(a::nf_elem, b::fmpq)
   iszero(b) && throw(DivideError())
   r = a.parent()
   ccall((:nf_elem_scalar_div_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         r, a, b, a.parent)
   return r
end

divexact(a::Integer, b::nf_elem) = inv(b)*a

divexact(a::fmpz, b::nf_elem) = inv(b)*a

divexact(a::fmpq, b::nf_elem) = inv(b)*a

###############################################################################
#
#   Removal and valuation
#
###############################################################################

doc"""
    divides(f::nf_elem, g::nf_elem)
> Returns a pair consisting of a flag which is set to `true` if $g$ divides
> $f$ and `false` otherwise, and a number field element $h$ such that $f = gh$
> if such exists. If not, the value of $h$ is undetermined.
"""
function divides(a::nf_elem, b::nf_elem)
  iszero(b) && throw(DivideError())
   return true, divexact(a, b)
end

###############################################################################
#
#   Norm and trace
#
###############################################################################

doc"""
    norm(a::nf_elem)
> Return the absolute norm of $a$. The result will be a rational number.
"""
function norm(a::nf_elem)
   z = fmpq()
   ccall((:nf_elem_norm, :libflint), Void,
         (Ref{fmpq}, Ref{nf_elem}, Ref{AnticNumberField}),
         z, a, a.parent)
   return z
end

doc"""
    norm(a::nf_elem)
> Return the absolute trace of $a$. The result will be a rational number.
"""
function trace(a::nf_elem)
   z = fmpq()
   ccall((:nf_elem_trace, :libflint), Void,
         (Ref{fmpq}, Ref{nf_elem}, Ref{AnticNumberField}),
         z, a, a.parent)
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(a::nf_elem)
   ccall((:nf_elem_zero, :libflint), Void,
         (Ref{nf_elem}, Ref{AnticNumberField}), a, parent(a))
   return a
end

function mul!(z::nf_elem, x::nf_elem, y::nf_elem)
   ccall((:nf_elem_mul, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
                                                  z, x, y, parent(x))
   return z
end

function mul_red!(z::nf_elem, x::nf_elem, y::nf_elem, red::Bool)
   ccall((:nf_elem_mul_red, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}, Cint),
                                                z, x, y, parent(x), red)
   return z
end

function addeq!(z::nf_elem, x::nf_elem)
   ccall((:nf_elem_add, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
                                                  z, z, x, parent(x))
   return z
end

function add!(a::nf_elem, b::nf_elem, c::nf_elem)
   ccall((:nf_elem_add, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         a, b, c, a.parent)
  return a
end

function reduce!(x::nf_elem)
   ccall((:nf_elem_reduce, :libflint), Void,
         (Ref{nf_elem}, Ref{AnticNumberField}), x, parent(x))
   return x
end

###############################################################################
#
#   Ad hoc unsafe functions
#
###############################################################################

function add!(c::nf_elem, a::nf_elem, b::fmpq)
   ccall((:nf_elem_add_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function add!(c::nf_elem, a::nf_elem, b::fmpz)
   ccall((:nf_elem_add_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function add!(c::nf_elem, a::nf_elem, b::Int)
   ccall((:nf_elem_add_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

add!(c::nf_elem, a::nf_elem, b::Integer) = add!(c, a, fmpz(b))

function sub!(c::nf_elem, a::nf_elem, b::fmpq)
   ccall((:nf_elem_sub_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function sub!(c::nf_elem, a::nf_elem, b::fmpz)
   ccall((:nf_elem_sub_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function sub!(c::nf_elem, a::nf_elem, b::Int)
   ccall((:nf_elem_sub_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

sub!(c::nf_elem, a::nf_elem, b::Integer) = sub!(c, a, fmpz(b))

function sub!(c::nf_elem, a::fmpq, b::nf_elem)
   ccall((:nf_elem_fmpq_sub, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpq}, Ref{nf_elem}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function sub!(c::nf_elem, a::fmpz, b::nf_elem)
   ccall((:nf_elem_fmpz_sub, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpz}, Ref{nf_elem}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function sub!(c::nf_elem, a::Int, b::nf_elem)
   ccall((:nf_elem_si_sub, :libflint), Void,
         (Ref{nf_elem}, Int, Ref{nf_elem}, Ref{AnticNumberField}),
         c, a, b, b.parent)
   return c
end

sub!(c::nf_elem, a::Integer, b::nf_elem) = sub!(c, fmpz(a), b)

function mul!(c::nf_elem, a::nf_elem, b::fmpq)
   ccall((:nf_elem_scalar_mul_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function mul!(c::nf_elem, a::nf_elem, b::fmpz)
   ccall((:nf_elem_scalar_mul_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

function mul!(c::nf_elem, a::nf_elem, b::Int)
   ccall((:nf_elem_scalar_mul_si, :libflint), Void,
         (Ref{nf_elem}, Ref{nf_elem}, Int, Ref{AnticNumberField}),
         c, a, b, a.parent)
   return c
end

mul!(c::nf_elem, a::nf_elem, b::Integer) = mul!(c, a, fmpz(b))

###############################################################################
#
#   Speedups for polynomials over number fields
#
###############################################################################

function sqr_classical(a::Generic.Poly{nf_elem})
   lena = length(a)

   t = base_ring(a)()

   lenz = 2*lena - 1
   d = Array{nf_elem}(lenz)

   for i = 1:lena - 1
      d[2i - 1] = base_ring(a)()
      d[2i] = base_ring(a)()
      d[2i - 1] = mul_red!(d[2i - 1], coeff(a, i - 1), coeff(a, i - 1), false)
   end
   d[2*lena - 1] = base_ring(a)()
   d[2*lena - 1] = mul_red!(d[2*lena - 1], coeff(a, lena - 1), coeff(a, lena - 1), false)

   for i = 1:lena
      for j = i + 1:lena
         t = mul_red!(t, coeff(a, i - 1), coeff(a, j - 1), false)
         d[i + j - 1] = addeq!(d[i + j - 1], t)
         d[i + j - 1] = addeq!(d[i + j - 1], t)
      end
   end

   for i = 1:lenz
      d[i] = reduce!(d[i])
   end

   z = parent(a)(d)

   set_length!(z, normalise(z, lenz))

   return z
end

function mul_classical(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem})
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   if lena == 0 || lenb == 0
      return parent(a)()
   end

   if a == b
       return sqr_classical(a)
   end

   t = base_ring(a)()

   lenz = lena + lenb - 1
   d = Array{nf_elem}(lenz)

   for i = 1:lena
      d[i] = base_ring(a)()
      d[i] = mul_red!(d[i], coeff(a, i - 1), coeff(b, 0), false)
   end

   for i = 2:lenb
      d[lena + i - 1] = base_ring(a)()
      d[lena + i - 1] = mul_red!(d[lena + i - 1], a.coeffs[lena], coeff(b, i - 1), false)
   end

   for i = 1:lena - 1
      for j = 2:lenb
         t = mul_red!(t, coeff(a, i - 1), b.coeffs[j], false)
         d[i + j - 1] = addeq!(d[i + j - 1], t)
      end
   end

   for i = 1:lenz
      d[i] = reduce!(d[i])
   end

   z = parent(a)(d)

   set_length!(z, normalise(z, lenz))

   return z
end

function *(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem})
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   if min(lena, lenb) < 20
      return mul_classical(a, b)
   end
   lenr = lena + lenb - 1
   r = parent(a)()
   if lena == 0 || lenb == 0
      return r
   end
   pol = base_ring(a).pol
   K = base_ring(a)
   R = parent(pol)
   T = elem_type(R)
   S = Generic.PolyRing{T}(R, :y)
   f = S()
   fit!(f, lena)
   for i = 1:lena
      f = setcoeff!(f, i - 1, R(coeff(a, i - 1)))
   end
   set_length!(f, lena)
   if a !== b
      g = S()
      fit!(g, lenb)
      for i = 1:lenb
         g = setcoeff!(g, i - 1, R(coeff(b, i - 1)))
      end
      set_length!(g, lenb)
   else
      g = f
   end
   p = f*g
   fit!(r, lenr)
   for i = 1:lenr
      r.coeffs[i] = K(p.coeffs[i])
   end
   set_length!(r, normalise(r, lenr))
   return r
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{nf_elem}, ::Type{T}) where {T <: Integer} = nf_elem

promote_rule(::Type{nf_elem}, ::Type{fmpz}) = nf_elem

promote_rule(::Type{nf_elem}, ::Type{fmpq}) = nf_elem

promote_rule(::Type{nf_elem}, ::Type{fmpq_poly}) = nf_elem

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (a::AnticNumberField)()
   z = nf_elem(a)
   ccall((:nf_elem_set_si, :libflint), Void,
         (Ref{nf_elem}, Int, Ref{AnticNumberField}), z, 0, a)
   return z
end

function (a::AnticNumberField)(c::Int)
   z = nf_elem(a)
   ccall((:nf_elem_set_si, :libflint), Void,
         (Ref{nf_elem}, Int, Ref{AnticNumberField}), z, c, a)
   return z
end

(a::AnticNumberField)(c::Integer) = a(fmpz(c))

function (a::AnticNumberField)(c::fmpz)
   z = nf_elem(a)
   ccall((:nf_elem_set_fmpz, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}), z, c, a)
   return z
end

function (a::AnticNumberField)(c::fmpq)
   z = nf_elem(a)
   ccall((:nf_elem_set_fmpq, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpq}, Ref{AnticNumberField}), z, c, a)
   return z
end

function (a::AnticNumberField)(b::nf_elem)
   parent(b) != a && error("Cannot coerce number field element")
   return b
end

function (a::AnticNumberField)(pol::fmpq_poly)
   pol = parent(a.pol)(pol) # check pol has correct parent
   z = nf_elem(a)
   if length(pol) >= length(a.pol)
      pol = mod(pol, a.pol)
   end
   ccall((:nf_elem_set_fmpq_poly, :libflint), Void,
         (Ref{nf_elem}, Ref{fmpq_poly}, Ref{AnticNumberField}), z, pol, a)
   return z
end

function (a::FmpqPolyRing)(b::nf_elem)
   parent(parent(b).pol) != a && error("Cannot coerce from number field to polynomial ring")
   r = a()
   ccall((:nf_elem_get_fmpq_poly, :libflint), Void,
         (Ref{fmpq_poly}, Ref{nf_elem}, Ref{AnticNumberField}), r, b, parent(b))
   return r
end

###############################################################################
#
#   AnticNumberField constructor
#
###############################################################################

doc"""
    AnticNumberField(f::fmpq_poly, s::AbstractString)
> Return a tuple $R, x$ consisting of the parent object $R$ and generator $x$
> of the number field $\mathbb{Q}/(f)$ where $f$ is the supplied polynomial.
> The supplied string `s` specifies how the generator of the number field
> should be printed.
"""
function AnticNumberField(f::fmpq_poly, s::AbstractString; cached = true)
   S = Symbol(s)
   parent_obj = AnticNumberField(f, S, cached)

   return parent_obj, gen(parent_obj)
end

doc"""
    AnticCyclotomicField(n::Int, s::AbstractString, t = "\$")
> Return a tuple $R, x$ consisting of the parent object $R$ and generator $x$
> of the $n$-th cyclotomic field, $\mathbb{Q}(\zeta_n)$. The supplied string
> `s` specifies how the generator of the number field should be printed. If
> provided, the string `t` specifies how the generator of the polynomial ring
> from which the number field is constructed, should be printed. If it is not
> supplied, a default dollar sign will be used to represent the variable.
"""
function AnticCyclotomicField(n::Int, s::AbstractString, t = "\$"; cached = true)
   Zx, x = PolynomialRing(FlintZZ, string(gensym()); cached = cached)
   Qx, = PolynomialRing(FlintQQ, t; cached = cached)
   f = cyclotomic(n, x)
   return AnticNumberField(Qx(f), s; cached = cached)
end

doc"""
    AnticMaximalRealSubfield(n::Int, s::AbstractString, t = "\$")
> Return a tuple $R, x$ consisting of the parent object $R$ and generator $x$
> of the totally real subfield of the $n$-th cyclotomic field,
> $\mathbb{Q}(\zeta_n)$. The supplied string `s` specifies how the generator of
> the number field should be printed. If provided, the string `t` specifies how
> the generator of the polynomial ring from which the number field is
> constructed, should be printed. If it is not supplied, a default dollar sign
> will be used to represent the variable.
"""
function AnticMaximalRealSubfield(n::Int, s::AbstractString, t = "\$"; cached = true)
   Zx, x = PolynomialRing(FlintZZ, string(gensym()); cached = cached)
   Qx, = PolynomialRing(FlintQQ, t; cached = cached)
   f = cos_minpoly(n, x)
   return AnticNumberField(Qx(f), s; cached = cached)
end


text:
using PyPlot
import LsqFit
export mainplot, estimate_embedding_plot

niceloc = matplotlib.ticker.MaxNLocator(nbins=4)

function mainplot(Hs, Cs, eHs, eCs, labels, legendtitle;
        qH = 1, qC = 2, tol = 0.25, offsets = zeros(length(Hs)),
        dimension_fit_H = linear_regression_fit_glm,
        dimension_fit_C = logarithmic_corrected_fit_lsqfit,
    )

    fig, axs = subplots(2,1; sharex = true)

    llines = []
    for j in 1:length(Hs)
        @info labels[j]
        z = offsets[j]
        eH = eHs isa Vector{<:AbstractVector} ? eHs[j] : eHs
        eC = eCs isa Vector{<:AbstractVector} ? eCs[j] : eCs
        color = "C$(j-1)"

        H = Hs[j]
        x, y = log.(eH), -H
        line, = axs[1].plot(x, y .+ z; color, ls = LINESTYLES[j], alpha = 0.9)
        push!(llines, line)
        is, d = linear_region(x, y; tol, warning = false)
        Δ, Δ05, Δ95 = dimension_fit_H(x[is[1]:is[2]], y[is[1]:is[2]])
        Hlabel = "\$$(rdspl(Δ05)), $(rdspl(Δ95))\$"
        axs[1].plot(x[[is[1], is[end]]], y[[is[1], is[end]]] .+ z, label = Hlabel,
        zorder = 2, ms = 10, marker = MARKERS[j], ls = "None", color = color, alpha = 0.75)

        C = Cs[j]
        i = findfirst(z -> z > 0, C)
        x, y = log.(eC)[i:end], log.(C)[i:end]
        is, d = linear_region(x, y; tol, warning = false)
        Δ, Δ05, Δ95 = dimension_fit_C(x[is[1]:is[2]], y[is[1]:is[2]])
        Clabel = "\$$(rdspl(Δ05)), $(rdspl(Δ95))\$"
        axs[2].plot(x, y .+ z, zorder = 1, color, ls = LINESTYLES[j], alpha = 0.9)
        axs[2].plot(x[[is[1], is[end]]], y[[is[1], is[end]]] .+ z, label = Clabel,
        zorder = 2, ms = 10, marker = MARKERS[j], ls = "None", color = color, alpha = 0.75)
        yield()
    end

    axs[1].set_ylabel("\$-H_$(qH)\$")
    axs[2].set_ylabel("\$\\log(C_$(qC))\$")
    axs[2].set_xlabel("\$\\log(\\varepsilon)\$")

    # Make the informative legend
    leg = axs[1].legend(
        handles = llines, labels = labels,
        bbox_to_anchor=(0., 1.02, 1., .102), loc="lower left",
        ncol=3, mode="expand", borderaxespad=0,
        title = legendtitle
    )

    for ax in axs
        ax.legend(;ncol = 3, handlelength = 0.2)
        # This throws ReadOnlyMemoryError for some reason...
        # axs[1].get_yaxis().set_major_locator(niceloc)
    end
    axs[1].add_artist(leg)
    fig.tight_layout(pad = 0.25)
    fig.subplots_adjust(top = 0.8)

    return fig, axs
end

function estimate_embedding_plot(x, t = "")
    fig = figure()
    ax = subplot(2,1,1)
    ax.plot(x)
    ax.set_title(t)

    mi = selfmutualinfo(x, 0:100)
    w = estimate_delay(x, "mi_min")
    ax = subplot(2,2,3)
    ax.plot(0:100, mi; label = "mut. inf.")
    ax.axvline(w; color = "C2")
    ax.legend()

    Es = delay_afnn(x, w, 2:10)
    ax = subplot(2,2,4)
    ax.plot(2:10, Es; label = "Cao") # saturates at d=6
    ax.legend()
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
mutable struct DebugRect <: MakieLayout.LObject
    layoutnodes::MakieLayout.LayoutNodes
    attributes::Attributes
end

function default_attributes(T::Type{DebugRect})
    Attributes(
        height = nothing,
        width = nothing,
        halign = :center,
        valign = :center,
        topprot = 0,
        leftprot = 0,
        rightprot = 0,
        bottomprot = 0
    )
end

function DebugRect(; bbox = nothing, kwargs...)
    attrs = merge!(Attributes(kwargs), default_attributes(DebugRect))

    @extract attrs (valign, halign, topprot, leftprot, rightprot, bottomprot)

    sizeattrs = MakieLayout.sizenode!(attrs.width, attrs.height)
    alignment = lift(tuple, halign, valign)

    suggestedbbox = MakieLayout.create_suggested_bboxnode(bbox)

    autosizenode = Node{NTuple{2, MakieLayout.Optional{Float32}}}((nothing, nothing))

    computedsize = MakieLayout.computedsizenode!(sizeattrs, autosizenode)

    finalbbox = MakieLayout.alignedbboxnode!(suggestedbbox, computedsize, alignment, sizeattrs, autosizenode)

    # no protrusions
    protrusions = lift(leftprot, rightprot, bottomprot, topprot) do l, r, b, t
        MakieLayout.RectSides{Float32}(l, r, b, t)
    end

    layoutnodes = MakieLayout.LayoutNodes{DebugRect, GridLayout}(suggestedbbox, protrusions, computedsize, autosizenode, finalbbox, nothing)

    # trigger bbox
    suggestedbbox[] = suggestedbbox[]

    DebugRect(layoutnodes, attrs)
end


text:
# arithmetic.jl

module ArithmeticTest

using Test
using Dates
using QDates

@testset "basic arithmetic" begin

qdt = QDates.QDate(1999,12,27)
@test qdt + Dates.Year(1) == QDates.QDate(2000,12,27)
@test qdt + Dates.Year(100) == QDates.QDate(2099,12,27)
@test_throws ArgumentError qdt + Dates.Year(1000)
@test qdt - Dates.Year(1) == QDates.QDate(1998,12,27)
@test qdt - Dates.Year(100) == QDates.QDate(1899,12,27)
@test qdt - Dates.Year(1000) == QDates.QDate(999,12,27)
qdt = QDates.QDate(2017,2,30)
@test qdt + Dates.Year(1) == QDates.QDate(2018,2,30)
@test qdt - Dates.Year(1) == QDates.QDate(2016,2,29)
qdt = QDates.QDate(2017,5,true,1)
@test_broken qdt + Dates.Year(1) == QDates.QDate(2018,6,1)
@test qdt + Dates.Year(1) == QDates.QDate(2018,5,1)
@test qdt + Dates.Year(11) == QDates.QDate(2028,5,true,1)
@test_broken qdt - Dates.Year(1) == QDates.QDate(2016,6,1)
@test qdt - Dates.Year(1) == QDates.QDate(2016,5,1)
@test qdt - Dates.Year(8) == QDates.QDate(2009,5,true,1)

qdt = QDates.QDate(1999,12,27)
@test qdt + Dates.Month(1) == QDates.QDate(2000,1,27)
@test qdt + Dates.Month(100) == QDates.QDate(2008,1,27)
@test qdt + Dates.Month(1000) == QDates.QDate(2080,10,27)
@test qdt - Dates.Month(1) == QDates.QDate(1999,11,27)
@test qdt - Dates.Month(100) == QDates.QDate(1991,11,27)
@test qdt - Dates.Month(1000) == QDates.QDate(1919,2,27)
qdt = QDates.QDate(2017,2,30)
@test qdt + Dates.Month(1) == QDates.QDate(2017,3,29)
@test qdt + Dates.Month(4) == QDates.QDate(2017,5,true,29)
@test qdt - Dates.Month(1) == QDates.QDate(2017,1,29)
@test qdt - Dates.Month(3) == QDates.QDate(2016,11,30)
qdt = QDates.QDate(2017,5,true,1)
@test qdt + Dates.Month(1) == QDates.QDate(2017,6,1)
@test qdt - Dates.Month(1) == QDates.QDate(2017,5,1)

# @test_throws MethodError qdt + Dates.Week(1)

qdt = QDates.QDate(1999,12,27)
@test qdt + Dates.Day(1) == QDates.QDate(1999,12,28)
@test qdt + Dates.Day(100) == QDates.QDate(2000,4,9)
@test qdt + Dates.Day(1000) == QDates.QDate(2002,9,24)
@test qdt - Dates.Day(1) == QDates.QDate(1999,12,26)
@test qdt - Dates.Day(100) == QDates.QDate(1999,9,17)
@test qdt - Dates.Day(1000) == QDates.QDate(1997,4,2)

end

# Vectorized arithmetic
@testset "Vectorized arithmetic" begin

a = QDates.QDate(2014,1,1)
dr = [a,a,a,a,a,a,a,a,a,a]
b = a + Dates.Year(1)
@test dr .+ Dates.Year(1) == repeat([b],10)
b = a + Dates.Month(1)
@test dr .+ Dates.Month(1) == repeat([b],10)
b = a + Dates.Day(1)
@test dr .+ Dates.Day(1) == repeat([b],10)
b = a - Dates.Year(1)
@test dr .- Dates.Year(1) == repeat([b],10)
b = a - Dates.Month(1)
@test dr .- Dates.Month(1) == repeat([b],10)
b = a - Dates.Day(1)
@test dr .- Dates.Day(1) == repeat([b],10)

# Month arithmetic minimizes "edit distance", or number of changes
# needed to get a correct answer
# This approach results in a few cases of non-associativity
a = QDates.QDate(2012,1,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,2,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,4,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,5,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,6,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,8,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,9,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,10,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,11,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = QDates.QDate(2012,12,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)

t1 = [QDates.QDate(2009,1,1) QDates.QDate(2009,1,2) QDates.QDate(2009,1,3); QDates.QDate(2009,2,1) QDates.QDate(2009,2,2) QDates.QDate(2009,2,3)]
t2 = [QDates.QDate(2009,1,2) QDates.QDate(2009,2,2) QDates.QDate(2010,1,3); QDates.QDate(2010,2,1) QDates.QDate(2009,3,2) QDates.QDate(2009,2,4)]

# TimeType, Array{TimeType}
@test QDates.QDate(2010,1,1) .- t1 == [Dates.Day(384) Dates.Day(383) Dates.Day(382); Dates.Day(354) Dates.Day(353) Dates.Day(352)]
@test t1 .- QDates.QDate(2010,1,1) == [Dates.Day(-384) Dates.Day(-383) Dates.Day(-382); Dates.Day(-354) Dates.Day(-353) Dates.Day(-352)]

# GeneralPeriod, Array{TimeType}
@test Dates.Day(1) .+ t1 == [QDates.QDate(2009,1,2) QDates.QDate(2009,1,3) QDates.QDate(2009,1,4); QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4)]
@test t1 .+ Dates.Day(1) == [QDates.QDate(2009,1,2) QDates.QDate(2009,1,3) QDates.QDate(2009,1,4); QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4)]
# @test Dates.Day(1) + t1 == [QDates.QDate(2009,1,2) QDates.QDate(2009,1,3) QDates.QDate(2009,1,4); QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4)]
# @test t1 + Dates.Day(1) == [QDates.QDate(2009,1,2) QDates.QDate(2009,1,3) QDates.QDate(2009,1,4); QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4)]

@test (Dates.Month(1) + Dates.Day(1)) .+ t1 == [QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4); QDates.QDate(2009,3,2) QDates.QDate(2009,3,3) QDates.QDate(2009,3,4)]
@test t1 .+ (Dates.Month(1) + Dates.Day(1)) == [QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4); QDates.QDate(2009,3,2) QDates.QDate(2009,3,3) QDates.QDate(2009,3,4)]
# @test (Dates.Month(1) + Dates.Day(1)) + t1 == [QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4); QDates.QDate(2009,3,2) QDates.QDate(2009,3,3) QDates.QDate(2009,3,4)]
# @test t1 + (Dates.Month(1) + Dates.Day(1)) == [QDates.QDate(2009,2,2) QDates.QDate(2009,2,3) QDates.QDate(2009,2,4); QDates.QDate(2009,3,2) QDates.QDate(2009,3,3) QDates.QDate(2009,3,4)]

@test t1 .- Dates.Day(1) == [QDates.QDate(2008,12,30) QDates.QDate(2009,1,1) QDates.QDate(2009,1,2); QDates.QDate(2009,1,30) QDates.QDate(2009,2,1) QDates.QDate(2009,2,2)]
# @test t1 - Dates.Day(1) == [QDates.QDate(2008,12,30) QDates.QDate(2009,1,1) QDates.QDate(2009,1,2); QDates.QDate(2009,1,30) QDates.QDate(2009,2,1) QDates.QDate(2009,2,2)]

@test t1 .- (Dates.Month(1) + Dates.Day(1)) == [QDates.QDate(2008,11,29) QDates.QDate(2008,12,1) QDates.QDate(2008,12,2); QDates.QDate(2008,12,30) QDates.QDate(2009,1,1) QDates.QDate(2009,1,2)]
# @test t1 - (Dates.Month(1) + Dates.Day(1)) == [QDates.QDate(2008,11,29) QDates.QDate(2008,12,1) QDates.QDate(2008,12,2); QDates.QDate(2008,12,30) QDates.QDate(2009,1,1) QDates.QDate(2009,1,2)]

# Array{TimeType}, Array{TimeType}
@test t2 .- t1 == [Dates.Day(1) Dates.Day(30) Dates.Day(384); Dates.Day(384) Dates.Day(30) Dates.Day(1)]
@test t2 - t1 == [Dates.Day(1) Dates.Day(30) Dates.Day(384); Dates.Day(384) Dates.Day(30) Dates.Day(1)]

end

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    RSOCtoPSDBridge{T} <: Bridges.Variable.AbstractBridge

`RSOCtoPSDBridge` implements the following reformulation:

 * ``||x||_2^2 \\le 2tu`` where ``t, u \\ge 0`` into ``Y \\succeq 0``, with the
   substitution rule:
   ``Y = \\left[\\begin{array}{c c}t & x^\\top \\\\ x & 2u \\mathbf{I}\\end{array}\\right].``

Additional bounds are added to ensure the off-diagonals of the ``2uI`` submatrix
are `0`, and linear constraints are added to ennsure the diagonal of ``2uI``
takes the same values.

As a special case, if ``|x|| = 0``, then `RSOCtoPSDBridge` reformulates into
``(t, u) \\in \\mathbb{R}_+``.

## Source node

`RSOCtoPSDBridge` supports:

 * [`MOI.VectorOfVariables`](@ref) in [`MOI.RotatedSecondOrderCone`](@ref)

## Target nodes

`RSOCtoPSDBridge` creates:

 * One variable node that depends on the input dimension:
   * [`MOI.VectorOfVariables`](@ref) in [`MOI.Nonnegatives`](@ref) if dimension
     is `1` or `2`
   * [`MOI.VectorOfVariables`](@ref) in
   [`MOI.PositiveSemidefiniteConeTriangle`](@ref) otherwise
 * The constraint node [`MOI.VariableIndex`](@ref) in [`MOI.EqualTo`](@ref)
 * The constrant node [`MOI.ScalarAffineFunction`](@ref) in [`MOI.EqualTo`](@ref)
"""
struct RSOCtoPSDBridge{T} <: AbstractBridge
    variables::Vector{MOI.VariableIndex}
    psd::Union{
        MOI.ConstraintIndex{
            MOI.VectorOfVariables,
            MOI.PositiveSemidefiniteConeTriangle,
        },
        MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives},
    }
    off_diag::Vector{MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}}}
    diag::Vector{
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
    }
end

const RSOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoPSDBridge{T},OT}

function bridge_constrained_variable(
    ::Type{RSOCtoPSDBridge{T}},
    model::MOI.ModelLike,
    set::MOI.RotatedSecondOrderCone,
) where {T}
    off_diag = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}}[]
    diag = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}[]
    # A special case: if there are no `x` variables, then a
    # RotatedSecondOrderCone is equivalent to [t, u] >= 0.
    if set.dimension <= 2
        variables, psd = MOI.add_constrained_variables(
            model,
            MOI.Nonnegatives(set.dimension),
        )
        return RSOCtoPSDBridge{T}(variables, psd, off_diag, diag)
    end
    psd_set = MOI.PositiveSemidefiniteConeTriangle((set.dimension - 1))
    variables, psd = MOI.add_constrained_variables(model, psd_set)
    zero_set = MOI.EqualTo(zero(T))
    k = 0
    for col in 1:(set.dimension-1), row in 1:col
        k += 1
        if row == 1 || col <= 2
            continue  # These are free variables
        elseif col == row
            # X[2, 2] - X[i, i] == 0 for i = 3:N
            f = MOI.Utilities.operate(-, T, variables[3], variables[k])
            push!(diag, MOI.add_constraint(model, f, zero_set))
        else
            push!(off_diag, MOI.add_constraint(model, variables[k], zero_set))
        end
    end
    return RSOCtoPSDBridge{T}(variables, psd, off_diag, diag)
end

function supports_constrained_variable(
    ::Type{<:RSOCtoPSDBridge},
    ::Type{MOI.RotatedSecondOrderCone},
)
    return true
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:RSOCtoPSDBridge})
    return Tuple{Type}[
        (MOI.PositiveSemidefiniteConeTriangle,),
        (MOI.Nonnegatives,),
    ]
end

function MOI.Bridges.added_constraint_types(
    ::Type{RSOCtoPSDBridge{T}},
) where {T}
    return Tuple{Type,Type}[
        (MOI.VariableIndex, MOI.EqualTo{T}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::RSOCtoPSDBridge, ::MOI.NumberOfVariables)::Int64
    return length(bridge.variables)
end

function MOI.get(bridge::RSOCtoPSDBridge, ::MOI.ListOfVariableIndices)
    return copy(bridge.variables)
end

function MOI.get(
    bridge::RSOCtoPSDBridge,
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,S},
)::Int64 where {S<:Union{MOI.PositiveSemidefiniteConeTriangle,MOI.Nonnegatives}}
    if bridge.psd isa MOI.ConstraintIndex{MOI.VectorOfVariables,S}
        return 1
    end
    return 0
end

function MOI.get(
    bridge::RSOCtoPSDBridge,
    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S},
) where {S<:Union{MOI.PositiveSemidefiniteConeTriangle,MOI.Nonnegatives}}
    if bridge.psd isa MOI.ConstraintIndex{MOI.VectorOfVariables,S}
        return [bridge.psd]
    end
    return MOI.ConstraintIndex{MOI.VectorOfVariables,S}[]
end

function MOI.get(
    bridge::RSOCtoPSDBridge{T},
    ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.EqualTo{T}},
)::Int64 where {T}
    return length(bridge.off_diag)
end

function MOI.get(
    bridge::RSOCtoPSDBridge{T},
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.EqualTo{T}},
) where {T}
    return copy(bridge.off_diag)
end

function MOI.get(
    bridge::RSOCtoPSDBridge{T},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
)::Int64 where {T}
    return length(bridge.diag)
end

function MOI.get(
    bridge::RSOCtoPSDBridge{T},
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
) where {T}
    return copy(bridge.diag)
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::RSOCtoPSDBridge)
    for ci in bridge.diag
        MOI.delete(model, ci)
    end
    MOI.delete(model, bridge.variables)
    return
end

# Attributes, Bridge acting as a constraint

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::RSOCtoPSDBridge)
    CI = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}
    if bridge.psd isa CI
        return MOI.RotatedSecondOrderCone(length(bridge.variables))
    end
    return MOI.RotatedSecondOrderCone(length(bridge.diag) + 3)
end

function _variable_map(bridge::RSOCtoPSDBridge, i::MOI.Bridges.IndexInVector)
    CI = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}
    if bridge.psd isa CI
        return i.value
    elseif i.value == 1
        return 1
    elseif i.value == 2
        return 3
    else
        return MOI.Utilities.trimap(1, i.value - 1)
    end
end

function _variable(bridge::RSOCtoPSDBridge, i::MOI.Bridges.IndexInVector)
    return bridge.variables[_variable_map(bridge, i)]
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::RSOCtoPSDBridge,
)
    values = MOI.get(model, attr, bridge.psd)
    set = MOI.get(model, MOI.ConstraintSet(), bridge)
    n = MOI.dimension(set)
    primal = zeros(n)
    for i in 1:n
        primal[i] = values[_variable_map(bridge, MOI.Bridges.IndexInVector(i))]
    end
    if n >= 2
        primal[2] /= 2
    end
    return primal
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDual,
    bridge::RSOCtoPSDBridge,
)
    values = MOI.get(model, attr, bridge.psd)
    n = MOI.dimension(MOI.get(model, MOI.ConstraintSet(), bridge))
    dual = zeros(n)
    for i in 1:n
        dual[i] = values[_variable_map(bridge, MOI.Bridges.IndexInVector(i))]
    end
    for ci in bridge.diag
        dual[2] += MOI.get(model, attr, ci)
    end
    # For `i = 2`, we multiply by 2 because it is `2u`.
    # For `i > 2`, we multiply by 2 because to account for the difference
    # of scalar product `MOI.Utilities.set_dot`.
    for i in 2:length(dual)
        dual[i] *= 2
    end
    return dual
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimal,
    bridge::RSOCtoPSDBridge,
    i::MOI.Bridges.IndexInVector,
)
    value = MOI.get(model, attr, _variable(bridge, i))
    if i.value == 2
        return value / 2
    end
    return value
end

function MOI.Bridges.bridged_function(
    bridge::RSOCtoPSDBridge{T},
    i::MOI.Bridges.IndexInVector,
) where {T}
    func = _variable(bridge, i)
    if i.value == 2
        return MOI.Utilities.operate(/, T, func, convert(T, 2))
    end
    return convert(MOI.ScalarAffineFunction{T}, func)
end

function unbridged_map(
    bridge::RSOCtoPSDBridge{T},
    vi::MOI.VariableIndex,
    i::MOI.Bridges.IndexInVector,
) where {T}
    if i.value == 2
        func = MOI.Utilities.operate(*, T, convert(T, 2), vi)
        return (_variable(bridge, i) => func,)
    end
    return (_variable(bridge, i) => convert(MOI.ScalarAffineFunction{T}, vi),)
end


text:
module CTProj

using DataFrames, JuMP,  NLopt, Distributions,
  StatsBase, Base.LinAlg.BLAS, Gadfly
#pyplot()

#=markdown path
weave(Pkg.dir("C:\\Users\\Clinton\\Dropbox\\Coursework\\Theory\\PS3\\SMM.jl"),
  informat="script",
  out_path = "C:\\Users\\Clinton\\Dropbox\\Coursework\\Theory\\PS3\\outSMM.html",
  doctype = "md2html")
=#

#start the program
#in/out: nothing
function p5SMM()
  flush(STDOUT)
  #ProcessCompData() #uncomment for data analysis

  const iNSims=1250
  const iBurnIn = 250
  vState=rand(Float64,iNSims,1)
  vStateTemp=rand(Float64,iNSims,1)
  const iSimsForSE = 100

  vMom = ones(4)

  println("Running at θ=0.6:")
  converge!(0.6, vMom = vMom, iNSims=iNSims,iBurnIn=iBurnIn, vState=vState)

  println("\nMean I/K: ", vMom[1],
    "\nσ(I/K): ", vMom[2],
    "\nMean Q: ",vMom[3],
    "\nσ(Q): ",vMom[4])

  mW = eye(4)

  #first stage
  println("Running First Stage")
  getθ(iNSims,iBurnIn,vState,mW,verbose=true)

  #get the standard errors
  mB = zeros(iSimsForSE,4)
  for i in 1:iSimsForSE
    mB[i,:] .= (getθ(iNSims,iBurnIn,rand(Float64,iNSims,1),mW,verbose=false))[:]
  end

  mB .= mB .- ones(iSimsForSE)*
    [mean(mB[:,1]) mean(mB[:,2]) mean(mB[:,3]) mean(mB[:,4])]

  mS = cov(mB)
  dnormConst = mean(abs.(mS))
  mS .= mS/dnormConst
  mW = mS\eye(4)
  println("Normalized Weighting Matrix:\n\n", mW)

  #second stage
  println("Running Second Stage")
  getθ(iNSims,iBurnIn,vState,mW,verbose=true)

end

function getθ(iNSims,iBurnIn,vState::Array{Float64,2},mW::Array{Float64,2};
  verbose=false)
  #the moments from the data

  const dIKMean = 0.26215608182706984
  const dIKσ = 0.18066965524433168
  const dQMean = 1.7303211980649078
  const dQσ = 0.8678322338204875
  const vMomData = [dIKMean, dIKσ, dQMean, dQσ]
  const iM = 4

  const steps = 100
  const dtol = 10.0^-7

  θLow = 0.01
  θHigh = .99
  θMid = (θHigh+θLow)/2.0

  JPrev = 0.0 #for measurement
  JMid = 0.0
  JLow = 0.0
  JHigh = 0.0

  iter = 0

  vMomSim = zeros(iM)
  vMomΔ= zeros(iM)

  #objective function, could be faster if devectorized
  #for notational convenience
  @inline function measureConverge(x::Float64)
    converge!(x, vMom = vMomSim, iNSims=iNSims,
      iBurnIn=iBurnIn, vState=vState)
    vMomΔ .= vMomData.-vMomSim
    return (vMomΔ'*mW*vMomΔ)[1]
  end

  #calculate values for initial J
  JMid = measureConverge(θMid)
  JLow = measureConverge(θLow)
  JHigh = measureConverge(θHigh)
  #println(JMid)

  #main loop
  while ((abs(JLow-JHigh) >= dtol) & (iter < 1000))

    #test points

    JHigh= measureConverge((θMid+θHigh)/2.0)
    JLow = measureConverge((θLow+θMid)/2.0)

    if JHigh<JMid #did we do better? if so, use the upper window
      θLow = θMid
      θMid = (θLow+θHigh)/2.0
      JMid = JHigh
    else
      if JLow<JMid #did we do better? if so, use the lower window
        θHigh = θMid
        θMid = (θLow+θHigh)/2.0
        JMid = JLow
      else #otherwise, just shrink the current window
        θHigh = (θMid + θHigh)/2.0
        θLow= (θMid + θLow)/2.0
      end
    end

    iter+=1

  end

  if verbose
    println("\nIter:", iter,
      "\nθ: ", θMid,
      "\nJ: ", JMid
    )

    println("\nMean I/K: ", vMomSim[1],
      "\nσ(I/K): ", vMomSim[2],
      "\nMean Q: ",vMomSim[3],
      "\nσ(Q): ",vMomSim[4])

    θ=0.01:.01:0.99
    J=map(x->min(measureConverge(x),10.0),θ)

    plOut = plot(melt(DataFrame(θ=θ,
        J=J),:θ),
      x=:θ,y=:value,color=:variable,Geom.line,
      Guide.title("J vs θ"),
      Guide.xlabel("θ"),Guide.ylabel("J"))
    draw(PNG("SMMPlotSecondStage.png", 7inch, 5inch), plOut)
    display(plOut)

  end

  return vMomSim
end

function ProcessCompData()

  #pre-processing data
  dfSMM= readtable("CompData.csv")
  display(summary(dfSMM))

  #drop unneeded columns
  dfSMM = dfSMM[:,[:gvkey,:fyear,:at,:capx,:ceq,:csho,:ppent,:txdb,:prcc_c]]

  #filter out NA values
  dfSMM = dfSMM[(~isna(dfSMM[:,:gvkey])&
    ~isna(dfSMM[:,:ppent])),:]

  dfSMM[:,:txdb] .= map!(((x)->isna(x)?0.0:x),dfSMM[:,:txdb])

  #calculate the market value
  dfSMM[:,:MV] = dfSMM[:,:prcc_c].*dfSMM[:,:csho]
  sort!(dfSMM,cols = (:gvkey,:fyear))
  dfSMM[:,:lppent] = zeros(size(dfSMM,1))

  #make sure it is ok do to do the lagging procedure
  dfSMM[:,:ycnt] = 0
  dcnt = 0

  by(dfSMM, :gvkey) do dfSelect
    dcnt = size(dfSelect,1)
    if maximum(dfSelect[:,:fyear])-minimum(dfSelect[1,:fyear])!=dcnt-1
      dcnt = 0
    end
    dfSelect[:,:ycnt] = dcnt
  end

  dfSMM = dfSMM[dfSMM[:,:ycnt].>1,:]

  by(dfSMM, :gvkey) do dfSelect
    dfSelect[:,:lppent] = [0.0; dfSelect[1:end-1,:ppent]]
  end

  #filter out more NA values
  dfSMM = dfSMM[(~isna(dfSMM[:,:csho])&
    ~isna(dfSMM[:,:prcc_c])&
    ~isna(dfSMM[:,:at])&
    ~isna(dfSMM[:,:capx])&
    ~isna(dfSMM[:,:ceq])),:]

  #filter out non-sensical inputs
  dfSMM = dfSMM[dfSMM[:,:MV].>=100,:]
  dfSMM = dfSMM[dfSMM[:,:at].>.001,:]
  dfSMM = dfSMM[dfSMM[:,:ceq].>=0.001,:]
  dfSMM = dfSMM[dfSMM[:,:lppent].>0.001,:]

  #Calculate a Q
  dfSMM[:,:Q] = (dfSMM[:,:at].+dfSMM[:,:MV].-dfSMM[:,:ceq].-
    dfSMM[:,:txdb])./dfSMM[:,:at]
  dfSMM = dfSMM[dfSMM[:,:Q].>=0,:]

  dfSMM[:,:IOverK] =  dfSMM[:,:capx]./dfSMM[:,:lppent]
  dfSMM = dfSMM[:,[:gvkey,:fyear,:at,:capx,:ceq,:MV,:lppent,:txdb,:Q,:IOverK]]

  #windsorization
  dIKLower = percentile(dfSMM[:,:IOverK],5)
  dIKUpper = percentile(dfSMM[:,:IOverK],95)
  dQLower = percentile(dfSMM[:,:Q],5)
  dQUpper = percentile(dfSMM[:,:Q],95)

  dfSMM = dfSMM[dfSMM[:,:IOverK].>dIKLower,:]
  dfSMM = dfSMM[dfSMM[:,:IOverK].<dIKUpper,:]
  dfSMM = dfSMM[dfSMM[:,:Q].>dQLower,:]
  dfSMM = dfSMM[dfSMM[:,:Q].<dQUpper,:]

  #debug stuff
  println("IOverK mean: ", mean(dfSMM[:,:IOverK]))
  println("Q mean: ", mean(dfSMM[:,:Q]))
  println("IOverK sd: ", var(dfSMM[:,:IOverK])^.5)
  println("Q sd: ", var(dfSMM[:,:Q])^.5)

  display(summary(dfSMM))

end

function converge!(θ::Float64; vMom::Vector{Float64}=[1.0,1.0,1.0,1.0],
  iNSims=1250, iBurnIn = 250,vState::Array{Float64,2}=rand(Float64,iNSims,1),
  iSteps=50)

  θ = (θ==1.0?1.01:θ)
  #println("θ: ", θ)

  const r=0.04
  const β=1/(1+r)
  #const θ=.6
  const δ=0.15

  const vy=[exp(-0.2);exp(0.2)]
  const mpiy=[.6 .4;.4 .6]
  const mpiyβ=β.*mpiy
  const iM = length(vy)

  #the approximate steady state is:
  const dkss = ((r+δ)/(θ*0.5*(exp(-0.2)+exp(0.2))))^(1/(θ-1))
  #println("kss=",dkss)
  #finite dimensional state space
  const dStep=dkss/(iSteps/2.0)
  const dbmax=2.0*dkss

  vb=dStep:dStep:dbmax    #grid steps for bonds are .1*avg. income

  ################begin: My Code
  dTol=10.0^-8.0
  iN = length(vb)

  #first we create our operator
  m3R = ones(iN,iN,iM)
  mI = ones(Int,iN,iM)
  vOneiN = ones(iN)
  mZeroiN = zeros(iN,iN)
  mTPlusV = Array{Float64}(iN,iN)

  for j = 1:iM
    #display(ones(iN,1)*((vb).^θ.*vy[j]+vb*(1-δ))'-vb*ones(1,iN))
    m3R[:,:,j].=ones(iN,1)*(vb.^θ.*vy[j]+vb.*(1-δ))'-vb*ones(1,iN)
  end

  map!((x)->x<0?-10^10:x,m3R) #switch negative consumption to a penalty
  #display(m3R)

  dErr = 1.0
  iIter = 1
  mV = ones(iN,iM)
  mVm1 = ones(iN,iM)

  #BLAS.set_num_threads(12)

  #tic()
  while dErr> dTol && iIter < 10000
      #to get policy and V
      for j = 1:iM
          mTPlusV.=m3R[:,:,j]
          for k = 1:iM
            ger!(mpiyβ[j,k],mVm1[:,k],vOneiN,mTPlusV) #BLAS optimization αxy'+M
          end
          mV[:,j],mI[:,j] = findmax(mTPlusV,1) #functions like matlab max
      end

      dErr = maximum(abs(mVm1.-mV))
      mVm1.=mV

      iIter = iIter + 1;
  end
  #toc()

  mI .= ((mI-1) .% iN)+1  #prior indices were linear, this fixes that
  #display(iIter)

  #Calc our stats by MC

  iState=Int8(2)
  ibCur=Int16(round(iN/2))
  ibNext=Int16(1)
  vStateHist = ones(Int8,iNSims,1)
  vbHist = ones(iNSims,1)

  vInvestment = ones(iNSims,1)
  vIOverK=ones(iNSims,1)
  vQ = ones(iNSims,1)

  #assign the state
  fNextState(x::Float64,s::Int8) =  Int8(x<mpiy[s,1]?1:2)

  for i=1:iNSims
    #update the state

    iState = fNextState(vState[i],iState)

    vStateHist[i] = iState
    ibNext = mI[ibCur,iState]

    vIOverK[i] = (vb[ibNext]-(1-δ)*vb[ibCur])
    vQ[i] = mV[ibCur,iState]/vb[ibCur]

    vbHist[i] = vb[ibNext]
    ibCur = ibNext
  end
  vIOverK[2:end].=vIOverK[2:end]./vbHist[1:(end-1)]

  vbHist = vbHist[iBurnIn:end]
  vStateHist = vStateHist[iBurnIn:end]
  vIOverK = vIOverK[iBurnIn:end]
  vQ = vQ[iBurnIn:end]

  vMom[1] = mean(vIOverK)
  vMom[2] = var(vIOverK)^.5
  vMom[3] = mean(vQ)
  vMom[4] = var(vQ)^.5

end

p5SMM()

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#cgrad(:default)[1]

```
function cMap(x::Number,cSeq::Array)
    n = length(cSeq)
    if x>1 || x<0
        error("x should be normalized to 0-1")
    end
    ind = Int(floor(x*(n-1)))+Int(1)
    if ind ==0 # out of bound caused by 0-eps()
        ind=1
    end
    cSeq[ind]
end
```
"cMap(x::Float,cSeq::Array) Map a value from 0 to one to a color, given a color sequence"
function cMap(x::Number,cSeq)
    cMap(x,cSeq.colors)
end

#cgrad(:default).colors

"cMap_def default c map"
function cMap_def(x::Number)
    cMap(x,cgrad(:default))
end

function cNormalize(x::Union{Array,UnitRange})
    xMin = minimum(x)
    xDiff = maximum(x) - xMin
    [(x[i]-xMin)/xDiff for i in 1:length(x)]
end

"qcMap_def(x) x is an arraym, quick map from a numeric array to a color array"
function qcMap_def(x)
    cMap_def.(cNormalize(x))
end

#qcMap_def(0:.05:1)

"plotInit"
function plotInit()
    return Plots.plot()
end

"plotTrajectory!(x,y,z) for x,y,z as single vector, plotTrajectory(x,y) for 2 vectors, plotTrajectory(x) is a time plot if x is the only input and ncol(x) is 1, else treated as a data-frame"
function plotTrajectory!(p,x::Union{Array,UnitRange},y::Union{Array,UnitRange},z::Union{Array,UnitRange};xlab = "x",ylab = "y",zlab = "z",map_t = true,t = Missing, log_t = false,markeralpha = .5)
    if (length(x)!=length(y)) || (length(y)!=length(z))
        error("Length differ")
    end
    if map_t
        if t==Missing
            t = 1:length(x)
        end
        if log_t == true
            t = log.(t)
        end
        return plot!(p,x,y,z;color = qcMap_def(t),markeralpha = markeralpha)
    end
    if length(t) != length(x)
        error("Check length of t")
    end
    return plot!(p,x,y,z;markeralpha = markeralpha)
end

function plotTrajectory!(p,x::Union{Array,UnitRange},y::Union{Array,UnitRange};xlab = "x",ylab = "y",map_t = true,t = Missing, log_t = false,markeralpha = .5)
    if (length(x)!=length(y))
        error("Length differ")
    end
    if map_t

        if t==Missing

            t = 1:length(x)
        end
        if log_t == true
            t = log.(t)
        end
        if length(t) != length(x)
            error("Check length of t")
        end
        return plot!(p,x,y;color = qcMap_def(t),markeralpha = markeralpha)
    end
    return plot!(p,x,y;markeralpha = markeralpha)
end

function plotTrajectory!(p,x::Union{Array,UnitRange};xlab = "t",ylab = "x",map_t = true,t = Missing,log_t = false,markeralpha = .5)

    if t==Missing           #  need t anyway

        t = 1:length(x)
    end
    if log_t == true
        t = log.(t)
    end
    if map_t
        return plot!(p,x;color = qcMap_def(t),markeralpha = markeralpha)
    end
    return plot!(p,t,x;markeralpha = markeralpha)
end

# plotScatter

"plotScatter!(x,y,z) for x,y,z as single vector, plotScatter(x,y) for 2 vectors, plotScatter(x) is a time plot if x is the only input and ncol(x) is 1, else treated as a data-frame"
function plotScatter!(p,x::Union{Array,UnitRange},y::Union{Array,UnitRange},z::Union{Array,UnitRange};xlab = "x",ylab = "y",zlab = "z",map_t = false,t = Missing,
     log_t = false,alpha = .5)
    if (length(x)!=length(y)) || (length(y)!=length(z))
        error("Length differ")
    end
    if map_t
        if t==Missing
            t = 1:length(x)
        end
        if log_t == true
            t = log.(t)
        end
        return Plots.scatter!(p,x,y,z;color = qcMap_def(t),alpha = alpha)
    end
    if length(t) != length(x)
        error("Check length of t")
    end
    return Plots.scatter!(p,x,y,z;alpha = alpha)
end

function plotScatter!(p,x::Union{Array,UnitRange},y::Union{Array,UnitRange};xlab = "x",ylab = "y",map_t = false,t = Missing, log_t = false,alpha = .5)
    if (length(x)!=length(y))
        error("Length differ")
    end
    if map_t

        if t==Missing

            t = 1:length(x)
        end
        if log_t == true
            t = log.(t)
        end
        if length(t) != length(x)
            error("Check length of t")
        end
        return Plots.scatter!(p,x,y;color = qcMap_def(t),alpha = alpha)
    end
    return Plots.scatter!(p,x,y;alpha = alpha)
end

function plotScatter!(p,x::Union{Array,UnitRange};xlab = "t",ylab = "x",map_t = false,t = Missing,log_t = false,alpha = .5)

    if t==Missing           #  need t anyway

        t = 1:length(x)
    end
    if log_t == true
        t = log.(t)
    end
    if map_t
        return Plots.scatter!(p,x;color = qcMap_def(t),alpha = alpha)
    end
    return Plots.scatter!(p,t,x;alpha = alpha)
end

"plotTrace!(x,y,z) similar to plotScatter, plotTrajectory with predefined plts"
function plotTrace!(p,x,y,z;t = Missing,log_t = false)
    plotTrajectory!(p,x,y,z;t=t,log_t = log_t)
    plotScatter!(p,x,y,z;t=t,log_t = log_t)
end

function plotTrace!(p,x,y;t = Missing,log_t = false)
    plotTrajectory!(p,x,y;t=t,log_t = log_t)
    plotScatter!(p,x,y;t=t,log_t = log_t)
end

function plotTrace!(p,x;t = Missing,log_t = false)
    plotTrajectory!(p,x;t=t,log_t = log_t)
    plotScatter!(p,x;t=t,log_t = log_t)
end


text:
using HorizonSideRobots

function inverse(side::HorizonSide)
    return HorizonSide((Int(side) + 2) % 4)
end

function average(r::Robot)
    nord_border = false
    side_border = false
    side = Ost
    temp_sum::Int = 0
    mark_num::Int = 0
    while !nord_border
        while !side_border
            if ismarker(r)
                temp_sum += temperature(r)
                mark_num += 1
            end
            if !isborder(r, side)
                move!(r, side)
            else
                side_border = true
            end
        end
        if !isborder(r, Nord)
            move!(r, Nord)
            side = inverse(side)
            side_border = false
        else
            nord_border = true
        end
    end
    print(temp_sum / mark_num)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
using Test

@testset "$(file)" for file in readdir(@__DIR__)
    if file == "DeprecatedTest.jl"
        continue
    end
    include(file)
end