text:
# Type definitions for test "Purpose: general implementation of getId works for user-defined purposes".
immutable Test1Purpose <: Purpose
  id::Symbol
end
immutable Test2Purpose <: Purpose
  not_id::Symbol
end
getId(p::Test2Purpose) = p.not_id

# Type definitions for test "DamOutput: general implementation of getId works for user-defined releases".
immutable Test1DamOutput <: DamOutput
  id::Symbol
end
immutable Test2DamOutput <: DamOutput
  not_id::Symbol
end
getId(p::Test2DamOutput) = p.not_id

# Actual tests.
facts("Reservoir data structures: basics") do
  function bathymetryTestShortVersions(g::ReservoirBathymetry)
    @fact isPolynomial(g) --> ispolynomial(g)
    @fact isLinear(g) --> islinear(g)
    @fact isQuadratic(g) --> isquadratic(g)
    @fact isPiecewiseLinear(g) --> ispiecewiselinear(g)
    @fact isPWL(g) --> ispiecewiselinear(g)
    @fact isSquareRoot(g) --> issquareroot(g)
  end

  function bathymetryTestValueShortVersions(g::ReservoirBathymetry)
    @fact volume(g, 2.0) --> computeVolume(g, 2.0)
    @fact volume(g, 4.0) --> computeVolume(g, 4.0)
    @fact level(g, 2.0) --> computeLevel(g, 2.0)
    @fact level(g, 4.0) --> computeLevel(g, 4.0)
  end

  const an_id = :a
  const a = 1.
  const a_discharge = a * 1.m^3/s
  const a_power = a * 1.W
  const a_height = a * 1.m
  const b = 2.
  const b_discharge = b * 1.m^3/s
  const b_power = b * 1.W
  const c = 3.
  const xs = [1.1, 2.0, 3.0, 4.0] # xs[1] = 1.1 for QuadraticReservoirBathymetry: avoid root 0 (so xs[1] != 1.0).

  context("bathymetry: NullReservoirBathymetry") do
    g = NullReservoirBathymetry()
    @fact isempty(g) --> true
    for x in xs; @fact_throws level(g, x); end
    for x in xs; @fact_throws volume(g, x); end
    @fact iswhitebox(g) --> false
    # What about isblackbox? No value makes sense...
    @fact ispolynomial(g) --> false
    @fact islinear(g) --> false
    @fact isquadratic(g) --> false
    @fact ispiecewiselinear(g) --> false
    @fact issquareroot(g) --> false
    bathymetryTestShortVersions(g)

    @fact_throws computeLevel(g, 2.0)
    @fact_throws level(g, 2.0)
    @fact_throws computeVolume(g, 2.0)
    @fact_throws volume(g, 2.0)
  end

  context("bathymetry: LinearReservoirBathymetry") do
    g = LinearReservoirBathymetry((a, b))
    @fact isempty(g) --> false
    for x in xs; @fact level(g, x) --> roughly(a + b * x); end
    for x in xs; @fact volume(g, x)  --> roughly((x - a) / b); end
    @fact iswhitebox(g) --> true
    @fact isblackbox(g) --> false
    @fact ispolynomial(g) --> true
    @fact islinear(g) --> true
    @fact isquadratic(g) --> false
    @fact ispiecewiselinear(g) --> false
    @fact issquareroot(g) --> false
    bathymetryTestShortVersions(g)
    bathymetryTestValueShortVersions(g)

    @fact coefficient(g, 1) --> b
    @fact coefficient(g, 0) --> a

    # Other constructors.
    g2 = LinearReservoirBathymetry(a, b)
    @fact g2 --> exactly(g)
  end

  context("bathymetry: QuadraticReservoirBathymetry") do
    g = QuadraticReservoirBathymetry((a, b, c))
    # Roots for this polynomial:
    #     {{v -> 1/3 (-1 - Sqrt[-2 + 3 d])}, {v -> 1/3 (-1 + Sqrt[-2 + 3 d])}}
    @fact isempty(g) --> false
    for x in xs; @fact level(g, x) --> roughly(a + b * x + c * x^2); end
    for x in xs; @fact volume(g, x)  --> roughly((- 1 + sqrt(- 2 + 3 * x))/ 3); end
    @fact iswhitebox(g) --> true
    @fact isblackbox(g) --> false
    @fact ispolynomial(g) --> true
    @fact islinear(g) --> false
    @fact isquadratic(g) --> true
    @fact ispiecewiselinear(g) --> false
    @fact issquareroot(g) --> false
    bathymetryTestShortVersions(g)
    bathymetryTestValueShortVersions(g)

    @fact coefficient(g, 2) --> c
    @fact coefficient(g, 1) --> b
    @fact coefficient(g, 0) --> a

    @fact_throws volume(QuadraticReservoirBathymetry((-1., -2., -3.)), 1.1) # No real root
    @fact_throws volume(QuadraticReservoirBathymetry((1.21, -2.2, 1.)), 1.1) # Two roots.
    @fact_throws volume(QuadraticReservoirBathymetry((1., 2., 3.)), 1.) # Root is zero.

    # Other constructors.
    g2 = QuadraticReservoirBathymetry(a, b, c)
    @fact g2 --> exactly(g)
  end

  context("bathymetry: PiecewiseLinearReservoirBathymetry") do
    @fact_throws PiecewiseLinearReservoirBathymetry([a, b, c], xs[1:3]) # Too many values
    @fact_throws PiecewiseLinearReservoirBathymetry([a, b, c], xs[1:1]) # Too few values

    g = PiecewiseLinearReservoirBathymetry([a, b, c], xs[1:2])
    @fact isempty(g) --> false
    # for x in xs; @fact level(g, x) --> roughly(...); end
    # for x in xs; @fact volume(g, x)  --> roughly(...); end
    @fact iswhitebox(g) --> true
    @fact isblackbox(g) --> false
    @fact ispolynomial(g) --> false
    @fact islinear(g) --> false
    @fact isquadratic(g) --> false
    @fact ispiecewiselinear(g) --> true
    @fact issquareroot(g) --> false
    bathymetryTestShortVersions(g)
    # TODO: bathymetryTestValueShortVersions(g)

    # TODO: values.
  end

  context("bathymetry: SquareRootReservoirBathymetry") do
    g = SquareRootReservoirBathymetry(a, b)
    @fact isempty(g) --> false
    for x in xs; @fact level(g, x) --> roughly(a + sqrt(b * x)); end
    for x in xs; @fact volume(g, x)  --> roughly((x - a)^2 / b); end
    @fact iswhitebox(g) --> true
    @fact isblackbox(g) --> false
    @fact ispolynomial(g) --> false
    @fact islinear(g) --> false
    @fact isquadratic(g) --> false
    @fact ispiecewiselinear(g) --> false
    @fact issquareroot(g) --> true
    bathymetryTestShortVersions(g)
    bathymetryTestValueShortVersions(g)

    # Other constructors.
    g2 = SquareRootReservoirBathymetry((a, b))
    @fact g2 --> exactly(g)
  end

  context("HydropowerUnit") do
    @fact_throws HydropowerUnit(max_discharge=a_discharge, max_power=b_power, efficiency=2.) # Wrong efficiency

    hp = HydropowerUnit(max_discharge=a_discharge, max_power=b_power, efficiency=1./c)
    @fact getId(hp) --> not(getPurposeId(hp))
    @fact getId(hp) --> not(getDamOutputId(hp))
    @fact getPurposeId(hp) --> not(getDamOutputId(hp))
    @fact getMaximumDischarge(hp) --> ReservoirManagement.ReservoirManagement._from_unitful(a_discharge)
    @fact getMaximumPower(hp) --> ReservoirManagement._from_unitful(b_power)
    @fact getEfficiency(hp) --> 1./c

    @fact id(hp) --> getId(hp)
    @fact purposeId(hp) --> getPurposeId(hp)
    @fact damOutputId(hp) --> getDamOutputId(hp)
    @fact maxDischarge(hp) --> getMaximumDischarge(hp)
    @fact maxPower(hp) --> getMaximumPower(hp)
    @fact efficiency(hp) --> getEfficiency(hp)
  end

  context("Purpose: general implementation of methods for user-defined purposes") do
    # Use the existing method, based on the id field.
    t = Test1Purpose(an_id)
    @fact getId(t) --> an_id
    @fact id(t) --> an_id
    @fact isWaterWithdrawal(t) --> false
    @fact isDeterministic(t) --> false

    # Use the user-defined getId.
    t = Test2Purpose(an_id)
    @fact getId(t) --> an_id
    @fact id(t) --> an_id
    @fact isWaterWithdrawal(t) --> false
    @fact isDeterministic(t) --> false
  end

  context("Purpose: DeterministicDrinkingWater and DeterministicEnvironmentalFlow") do
    dw = DeterministicDrinkingWater(a * 1.m^3/s)
    ef = DeterministicEnvironmentalFlow(b * 1.m^3/s)

    @fact getId(dw) --> not(getId(ef))
    @fact getNeed(dw) --> a * 1e-6
    @fact getNeed(ef) --> b * 1e-6
    @fact isHydropower(dw) --> false
    @fact isHydropower(ef) --> false
    @fact isWaterWithdrawal(dw) --> true
    @fact isWaterWithdrawal(ef) --> true
    @fact isDeterministic(dw) --> true
    @fact isDeterministic(ef) --> true

    @fact getNeed(dw, Day(1)) --> roughly(86400 * a * 1e-6)
    @fact getNeed(ef, Day(1)) --> roughly(86400 * b * 1e-6)
    @fact getNeed(dw, Week(1)) --> roughly(86400 * a * 1e-6 * 7)
    @fact getNeed(ef, Week(1)) --> roughly(86400 * b * 1e-6 * 7)

    @fact need(dw) --> getNeed(dw)
    @fact need(ef) --> getNeed(ef)
  end

  context("Purpose: list constructor (DeterministicPurposes)") do
    dw = DeterministicDrinkingWater(a_discharge)
    ef = DeterministicEnvironmentalFlow(b_discharge)
    list = DeterministicPurposes(drinkingWater=a_discharge, environmentalFlow=b_discharge)
    @fact in(dw, list) --> true
    @fact in(ef, list) --> true
    @fact length(list) --> 2

    list = DeterministicPurposes(drinkingWater=a_discharge)
    @fact in(dw, list) --> true
    @fact in(ef, list) --> false
    @fact length(list) --> 1

    list = DeterministicPurposes(environmentalFlow=b_discharge)
    @fact in(dw, list) --> false
    @fact in(ef, list) --> true
    @fact length(list) --> 1
  end

  context("Purpose: DeterministicHydropower") do
    hp = HydropowerUnit(max_discharge=a_discharge, max_power=b_power, efficiency=1./c)
    p = DeterministicHydropower(hp)

    @fact getId(p) --> not(getId(hp))
    @fact getHydropowerUnit(p) --> hp
    @fact getNeed(p) --> 0.
    @fact getNeed(p, Day(1)) --> 0.
    @fact getNeed(p, Week(1)) --> 0.
    @fact isHydropower(p) --> true

    @fact need(p) --> getNeed(p)
    @fact hydropower(p) --> getHydropowerUnit(p)
  end

  function outputTestShortVersions(g::DamOutput)
    @fact isConditional(g) --> hasCondition(g)
    @fact id(g) --> getId(g)
    @fact condition(g) --> getCondition(g)
  end

  function outputTestShortVersions(g::ConveyanceDamOutput)
    @fact isConditional(g) --> hasCondition(g)
    @fact id(g) --> getId(g)
    @fact condition(g) --> getCondition(g)
    @fact conveyance(g) --> getConveyanceFactor(g)
  end

  function outputTestShortVersions(g::ConstantDamOutput)
    @fact isConditional(g) --> hasCondition(g)
    @fact id(g) --> getId(g)
    @fact condition(g) --> getCondition(g)
    @fact discharge(g) --> getDischarge(g)
  end

  function outputTestShortVersions(g::HydropowerDamOutput)
    @fact isConditional(g) --> hasCondition(g)
    @fact id(g) --> getId(g)
    @fact condition(g) --> getCondition(g)
    @fact hydropower(g) --> getHydropowerUnit(g)
  end

  context("DamOutput: ConstantDamOutput") do
    dout1 = ConstantDamOutput(a_discharge)
    @fact hasCondition(dout1) --> false
    @fact getDischarge(dout1) --> ReservoirManagement._from_unitful(a_discharge)
    @fact isConstant(dout1) --> true
    @fact isHydropower(dout1) --> false
    @fact hasConveyanceFactor(dout1) --> false
    outputTestShortVersions(dout1)

    cond = MinimumReservoirLevelCondition(a_height)
    dout2 = ConstantDamOutput(a_discharge, cond)
    @fact hasCondition(dout2) --> true
    @fact condition(dout2) --> exactly(cond)
    @fact getDischarge(dout2) --> ReservoirManagement._from_unitful(a_discharge)
    @fact isConstant(dout2) --> true
    @fact isHydropower(dout2) --> false
    @fact hasConveyanceFactor(dout2) --> false
    outputTestShortVersions(dout2)

    cond = MinimumReservoirLevelCondition(a_height)
    dout3 = ConstantDamOutput(an_id, a_discharge, cond)
    @fact getId(dout3) --> an_id
    @fact hasCondition(dout3) --> true
    @fact condition(dout3) --> exactly(cond)
    @fact getDischarge(dout3) --> ReservoirManagement._from_unitful(a_discharge)
    @fact isConstant(dout3) --> true
    @fact isHydropower(dout3) --> false
    @fact hasConveyanceFactor(dout3) --> false
    outputTestShortVersions(dout3)
  end

  context("DamOutput: HydropowerDamOutput") do
    hp = HydropowerUnit(max_discharge=a_discharge, max_power=b_power, efficiency=1./c)
    dout1 = HydropowerDamOutput(hp)
    @fact hasCondition(dout1) --> false
    @fact getDischarge(dout1) --> getMaximumDischarge(hp)
    @fact isConstant(dout1) --> false
    @fact isHydropower(dout1) --> true
    @fact hasConveyanceFactor(dout1) --> false
    @fact getId(dout1) --> not(getId(hp))
    outputTestShortVersions(dout1)

    cond = MinimumReservoirLevelCondition(a_height)
    dout2 = HydropowerDamOutput(hp, cond)
    @fact hasCondition(dout2) --> true
    @fact condition(dout2) --> exactly(cond)
    @fact getDischarge(dout2) --> ReservoirManagement._from_unitful(a_discharge)
    @fact isConstant(dout2) --> false
    @fact isHydropower(dout2) --> true
    @fact hasConveyanceFactor(dout2) --> false
    @fact getId(dout2) --> not(getId(hp))
    outputTestShortVersions(dout2)
  end

  context("DamOutput: ConveyanceSpillway") do
    dout1 = ConveyanceSpillway(a)
    @fact hasCondition(dout1) --> false
    @fact getConveyanceFactor(dout1) --> a
    @fact isConstant(dout1) --> false
    @fact isHydropower(dout1) --> false
    @fact hasConveyanceFactor(dout1) --> true
    outputTestShortVersions(dout1)

    cond = MinimumReservoirLevelCondition(a_height)
    dout2 = ConveyanceSpillway(a, cond)
    @fact hasCondition(dout2) --> true
    @fact condition(dout2) --> exactly(cond)
    @fact getConveyanceFactor(dout2) --> a
    @fact isConstant(dout2) --> false
    @fact isHydropower(dout2) --> false
    @fact hasConveyanceFactor(dout2) --> true
    outputTestShortVersions(dout2)

    cond = MinimumReservoirLevelCondition(a_height)
    dout3 = ConveyanceSpillway(an_id, a, cond)
    @fact getId(dout3) --> an_id
    @fact hasCondition(dout3) --> true
    @fact condition(dout3) --> exactly(cond)
    @fact getConveyanceFactor(dout3) --> a
    @fact isConstant(dout3) --> false
    @fact isHydropower(dout3) --> false
    @fact hasConveyanceFactor(dout3) --> true
    outputTestShortVersions(dout3)
  end

  context("DamOutput: general implementation of methods for user-defined releases") do
    # Use the existing method, based on the id field.
    t1 = Test1DamOutput(an_id)
    @fact getId(t1) --> an_id
    @fact id(t1) --> an_id

    # Use the user-defined getId.
    t1 = Test2DamOutput(an_id)
    @fact getId(t1) --> an_id
    @fact id(t1) --> an_id
  end
end


text:
"""
The `criterion` module is used to implement various termination (or restart
criteria), such as the number of generations or evaluation passed, or the
amount of an resource consumed, in a flexible manner.
"""
module criterion
  using   state
  export  is_satisfied,
          Criterion,
          EvaluationsCriterion,
          IterationsCriterion

  abstract Criterion

  """
  Limits the number of evaluations an algorithm may perform.
  """
  type EvaluationsCriterion <: Criterion
    limit::Int
    EvaluationsCriterion(limit::Int) = new(limit)
  end

  """
  Limits the number of iterations, or generations, an algorithm may run for.
  """
  type IterationsCriterion <: Criterion
    limit::Int
    IterationsCriterion(limit::Int) = new(limit)
  end

  """
  Limits the algorithm to a fixed number of evaluations. Treated as a hard
  limit, i.e. algorithms should be guaranteed not to exceed this limit, even if
  there are several candidates pending evaluation in the current generation.

  **Parameters:**\n
  `limit::Int`, the evaluation limit.
  """
  evaluations(limit::Int) = EvaluationsCriterion(limit)

  """
  Limits the number of iterations that an algorithm may run for.

  **Parameters:**\n
  `limit::Int`, the iteration limit.
  """
  iterations(limit::Int) = IterationsCriterion(limit)

  """
  Limits the number of generations that an algorithm may run for.

  **Parameters:**\n
  `limit::Int`, the maximum number of generations.
  """
  generations = iterations

  """
  Determines whether a given evaluation limit has been hit by a provided
  algorithm state.
  """
  is_satisfied(c::EvaluationsCriterion, s::State) = s.evaluations >= c.limit

  """
  Determines whether a given iteration limit has been reached by a provided
  algorithm state.
  """
  is_satisfied(c::IterationsCriterion, s::State) = s.iterations >= c.limit
end

ext: jl


text:
######################
### Mutable states ###
######################

mutable struct DAState{TI<:Integer,TF<:Real}
    m     :: TI
    ϵ     :: TF
    μ     :: TF
    x_bar :: TF
    H_bar :: TF
end

function DAState(ϵ::Real)
    μ = computeμ(ϵ) # NOTE: this inital values doesn't affect anything as they will be overwritten
    return DAState(0, ϵ, μ, 0.0, 0.0)
end

function computeμ(ϵ::Real)
    return log(10 * ϵ) # see NUTS paper sec 3.2.1
end

function reset!(dastate::DAState{TI,TF}) where {TI<:Integer,TF<:Real}
    dastate.m = zero(TI)
    dastate.x_bar = zero(TF)
    dastate.H_bar = zero(TF)
end

mutable struct MSSState{T<:Real}
    ϵ :: T
end

################
### Adapters ###
################

abstract type StepSizeAdapter <: AbstractAdapter end

struct FixedStepSize{T<:Real} <: StepSizeAdapter
    ϵ :: T
end

function getss(fss::FixedStepSize)
    return fss.ϵ
end

struct DualAveraging{TI<:Integer,TF<:Real} <: StepSizeAdapter
  γ     :: TF
  t_0   :: TF
  κ     :: TF
  δ     :: TF
  state :: DAState{TI,TF}
end

function DualAveraging(spl::Sampler{<:AdaptiveHamiltonian}, ::Nothing, ϵ::Real)
    return DualAveraging(0.05, 10.0, 0.75, spl.alg.delta, DAState(ϵ))
end

function getss(da::DualAveraging)
    return da.state.ϵ
end

struct ManualSSAdapter{T<:Real} <:StepSizeAdapter
    state :: MSSState{T}
end

function getss(mssa::ManualSSAdapter)
    return mssa.state.ϵ
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/stepsize_adaptation.hpp
function adapt_stepsize!(da::DualAveraging, stats::Real)
    @debug "adapting step size ϵ..."
    @debug "current α = $(stats)"
    da.state.m += 1
    m = da.state.m

    # Clip average MH acceptance probability.
    stats = stats > 1 ? 1 : stats

    γ = da.γ; t_0 = da.t_0; κ = da.κ; δ = da.δ
    μ = da.state.μ; x_bar = da.state.x_bar; H_bar = da.state.H_bar

    η_H = 1.0 / (m + t_0)
    H_bar = (1.0 - η_H) * H_bar + η_H * (δ - stats)

    x = μ - H_bar * sqrt(m) / γ            # x ≡ logϵ
    η_x = m^(-κ)
    x_bar = (1.0 - η_x) * x_bar + η_x * x

    ϵ = exp(x)
    @debug "new ϵ = $(ϵ), old ϵ = $(da.state.ϵ)"

    if isnan(ϵ) || isinf(ϵ)
        @warn "Incorrect ϵ = $ϵ; ϵ_previous = $(da.state.ϵ) is used instead."
    else
        da.state.ϵ = ϵ
    end
    da.state.x_bar = x_bar
    da.state.H_bar = H_bar
end

function adapt!(da::DualAveraging, stats::Real, is_updateμ::Bool)
    adapt_stepsize!(da, stats)
    if is_updateμ
        da.state.μ = computeμ(da.state.ϵ)
        reset!(da.state)
    end
end


text:
"""
    TokenList
Tokenizer accessor.
"""
mutable struct TokenList
    ptr::Ptr{CXToken}
    size::Cuint
    tu::CXTranslationUnit
    function TokenList(tu::CXTranslationUnit, sr::CXSourceRange)
        ptr_ref = Ref{Ptr{CXToken}}(C_NULL)
        num_ref = Ref{Cuint}(0)
        clang_tokenize(tu, sr, ptr_ref, num_ref)
        @assert ptr_ref[] != C_NULL
        obj = new(ptr_ref[], num_ref[], tu)
        finalizer(obj) do x
            if x.ptr != C_NULL
                clang_disposeTokens(x.tu, x.ptr, x.size)
                x.ptr = C_NULL
            end
        end
        return obj
    end
end

Base.firstindex(x::TokenList) = 1
Base.lastindex(x::TokenList) = length(x)
Base.length(x::TokenList) = x.size
Base.iterate(x::TokenList, state=1) = state > x.size ? nothing : (x[state], state + 1)

function Base.getindex(x::TokenList, i::Integer)
    !checkindex(Bool, 1:length(x), i) &&
        throw(ArgumentError("BoundsError: TokenList miss at index $i"))
    GC.@preserve x begin
        token = unsafe_load(x.ptr, i)
        token_kind = kind(token)
        token_content = spelling(x.tu, token)
    end
    return CLTokenMap[token_kind](token, token_kind, token_content)
end

"""
    kind(t::CXToken) -> CXTokenKind
    kind(t::CLToken) -> CXTokenKind
Return the kind of the given token.
"""
kind(t::CXToken) = clang_getTokenKind(t)
kind(t::CLToken) = t.kind

"""
    spelling(tu::TranslationUnit, t::CLToken) -> String
    spelling(tu::TranslationUnit, t::CXToken) -> String
    spelling(tu::CXTranslationUnit, t::CXToken) -> String
Return the spelling of the given token. The spelling of a token is the textual
representation of that token, e.g., the text of an identifier or keyword.
"""
function spelling(tu::CXTranslationUnit, t::CXToken)
    cxstr = clang_getTokenSpelling(tu, t)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end
spelling(tu::TranslationUnit, t::CXToken) = spelling(tu.ptr, t)
spelling(tu::TranslationUnit, t::CLToken) = spelling(tu.ptr, t.token)

"""
    getTokenLocation(tu::Union{CXTranslationUnit,TranslationUnit}, t::Union{CXToken,CLToken}) -> CXSourceLocation
Return the source location of the given token.
"""
function getTokenLocation(tu::Union{CXTranslationUnit,TranslationUnit}, t::Union{CXToken,CLToken})
    return clang_getTokenLocation(tu, t)
end

"""
    getTokenExtent(tu::Union{CXTranslationUnit,TranslationUnit}, t::Union{CXToken,CLToken}) -> CXSourceRange
Return a source range that covers the given token.
"""
function getTokenExtent(tu::Union{CXTranslationUnit,TranslationUnit}, t::Union{CXToken,CLToken})
    return clang_getTokenExtent(tu, t)
end

"""
    tokenize(c::Union{CXCursor,CLCursor}) -> TokenList
Return a TokenList from the given cursor.
"""
function tokenize(c::Union{CXCursor,CLCursor})
    tu = clang_Cursor_getTranslationUnit(c)
    source_range = clang_getCursorExtent(c)
    range_start = clang_getRangeStart(source_range)
    range_end = clang_getRangeEnd(source_range)

    file = Ref{CXFile}(C_NULL)
    start_line = Ref{Cuint}(0)
    end_line = Ref{Cuint}(0)
    start_column = Ref{Cuint}(0)
    end_column = Ref{Cuint}(0)
    offset = Ref{Cuint}(0)
    clang_getExpansionLocation(range_start, file, start_line, start_column, offset)
    clang_getExpansionLocation(range_end, file, end_line, end_column, offset)

    expanded_start = clang_getLocation(tu, file[], start_line[], start_column[])
    expanded_end = clang_getLocation(tu, file[], end_line[], end_column[])
    expanded_range = clang_getRange(expanded_start, expanded_end)

    return TokenList(tu, expanded_range)
end

"""
    annotateTokens(tu::TranslationUnit, tokens, token_num, cursors)
Annotate the given set of tokens by providing cursors for each token that can be mapped to
a specific entity within the abstract syntax tree.
"""
function annotateTokens(tu::Union{CXTranslationUnit,TranslationUnit}, tokens, token_num, cursors)
    return clang_annotateTokens(tu, tokens, Unsigned(token_num), cursors)
end


text:
#-----------------------------------------------------------------------------------
# Testing SGD and SMD, clipped and non-clipped models
# Step size agnostic plots: covergence at range of stepsizes
#-----------------------------------------------------------------------------------

using Random, PyPlot, LaTeXStrings, Printf
include("solve_cov_est.jl");

Random.seed!(321);  # for reproducibility

#-------------------------------------------------------------
#    What two methods are we comparing?
#-------------------------------------------------------------
method1 = "subgradient"
method2 = "mirror"
clipped1 = false
clipped2 = true

#-------------------------------------
#    Set parameters
#-------------------------------------
maxIter = 9000;
r = 10;  # rank
d = 100;
stoch_err = 0.01;  # standard deviation of errors in stochastic measurements b
init_radius = 1.0;  # 2-norm measure of relative deviation between Xtrue and Xinit

stepSizes = 10.0 .^(range(-6,2,length=17))
num_rounds = 10;

#-----------------------------------------------------------
#   Initialize vectors to track final errors
#-----------------------------------------------------------
dist_errors_1 = fill(NaN, length(stepSizes));
fun_errors_1 = fill(NaN, length(stepSizes));
dist_errors_2 = fill(NaN, length(stepSizes));
fun_errors_2 = fill(NaN, length(stepSizes));

#-----------------------------------------------------------------------------------------
#   Run the two methods for each choice of stepsize
#-----------------------------------------------------------------------------------------
for i=1:length(stepSizes)
    η = stepSizes[i];

    dist_errors_tally_1 = 0;
    fun_errors_tally_1 = 0;
    dist_errors_tally_2 = 0;
    fun_errors_tally_2 = 0;

    for j=1:num_rounds
        # Generate true matrix and initialization
        Xtrue = randn(d,r);
        pert = randn(d,r);
        Xinit = Xtrue + init_radius * (norm(Xtrue, 2) / norm(pert, 2)) * pert;

        # Run the two methods and record final errors
        @printf("Running Method 1 for stepsize %1.2e (Round %i of %i)\n", η, 2*num_rounds*(i-1)+2*j-1, 2*length(stepSizes)*num_rounds);
        (err_hist_1, fun_hist_1) = solve_cov_est(Xinit, Xtrue, stoch_err, maxIter,
                                                       fill(η, maxIter), method=method1, clipped=clipped1, verbose=false)
        dist_errors_tally_1 += err_hist_1[end];
        fun_errors_tally_1 += fun_hist_1[end];

        @printf("Running Method 2 for stepsize %1.2e (Round %i of %i)\n", η, 2*num_rounds*(i-1)+2*j, 2*length(stepSizes)*num_rounds);
        (err_hist_2, fun_hist_2)  = solve_cov_est(Xinit, Xtrue, stoch_err, maxIter,
                                                        fill(η, maxIter), method=method2, clipped=clipped2, verbose=false)
        dist_errors_tally_2 += err_hist_2[end];
        fun_errors_tally_2 += fun_hist_2[end];

    end

    dist_errors_1[i] = dist_errors_tally_1 / num_rounds;
    fun_errors_1[i] = fun_errors_tally_1 /num_rounds;

    dist_errors_2[i] = dist_errors_tally_2 / num_rounds;
    fun_errors_2[i] = fun_errors_tally_2 / num_rounds;

end

distance_plot = figure(figsize=[10,6]);
title(@sprintf("Relative Distance Errors (r=%i, d=%i, %i iterations)", r,d,maxIter));
xlabel(L"Stepsize $\eta$");
ylabel(L"Relative distance to solution set $\min_U \; \frac{\| X_k U - X_{true} \|_2^2 }{\| X_{true} \|_2^2}$");
loglog(stepSizes, dist_errors_1, label=string(method1, clipped1 ? " clipped" : "") );
loglog(stepSizes, dist_errors_2, label=string(method2, clipped2 ? " clipped" : ""));
ylim(1e-5,1e3)
legend(loc="upper right")
savefig("plots/agnostic_cov_est_distances.pdf");

fun_plot = figure(figsize=[10,6]);
title(@sprintf("Absolute Function Errors (r=%i, d=%i, %i iterations)", r,d,maxIter));
xlabel(L"Stepsize $\eta$");
ylabel(L"Empirical Function Error $\quad \hat f (X_{final}) - \hat f (X_{true})$");
loglog(stepSizes, fun_errors_1, label=string(method1, clipped1 ? " clipped" : "") );
loglog(stepSizes, fun_errors_2, label=string(method2, clipped2 ? " clipped" : "") );
ylim(1e-5,1e10)
legend(loc="upper right")
savefig("plots/agnostic_cov_est_values.pdf");

max_stars_repo_name: kelliemac/stochastic-model-based

max_issues_repo_name: kelliemac/stochastic-model-based

max_forks_repo_name: kelliemac/stochastic-model-based


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to string processing.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _crop_str(str, crop_size, lstr = -1)

Return a cropped string of `str` with size `crop_size`. Notice that if the last
character before the crop does not fit due to its width, then blank spaces are
added.

The size of the string can be passed to `lstr` to save computational burden. If
`lstr = -1`, then the string length will be computed inside the function.

"""
function _crop_str(str::String, crop_size::Int, lstr::Int = -1)
    lstr < 0 && (lstr = textwidth(str))

    # If the `crop_size` is large than the string length, then just return the
    # current string.
    crop_size ≥ lstr && return str

    # Process every character of the string.
    cstr  = ""
    lcstr = 0

    # Find all ANSI escape sequences and split string there.
    #
    # The regex was obtained from here:
    #
    #     https://stackoverflow.com/questions/14693701/how-can-i-remove-the-ansi-escape-sequences-from-a-string-in-python
    #
    ansi   = collect(eachmatch(r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])", str))
    tokens = split(str, r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])")
    crop_finished = false

    for i = 1:length(tokens)
        for c in tokens[i]
            sc = textwidth(c)

            # Check if we can add the character without passing the crop size.
            if lcstr + sc ≤ crop_size
                cstr  *= c
                lcstr += sc

                # Fill the remaining spaces with blank characters.
            else
                Δ      = crop_size - lcstr
                cstr  *= " "^Δ
                crop_finished = true
                break
            end
        end

        crop_finished && break

        cstr *= ansi[i].match
    end

    return cstr
end

"""
    _render_text(T, v; compact_printing::Bool = true, isstring::Bool = false, linebreaks::Bool = false)

Render the value `v` to strings using the rendered `T` to be displayed in the
text back-end.

`T` can be:

* `Val(:print)`: the function `print` will be used.
* `Val(:show)`: the function `show` will be used.

This function must return a vector of strings in which each element is a line
inside the rendered cell.

If `linebreaks` is `true`, then the rendered should split the created string
into multiple tokens.

In case `show` is used, if `isstring` is `false`, then it means that the
original data is not a string even if `v` is a string. Hence, the surrounding
quotes added by `show` will be removed. This is required to correctly handle
formatters.

"""
@inline function _render_text(::Val{:print}, v;
                              compact_printing::Bool = true,
                              isstring::Bool = false,
                              linebreaks::Bool = false)

    str = sprint(print, v; context = :compact => compact_printing)

    return _render_text(Val(:print), str;
                        compact_printing = compact_printing,
                        isstring = isstring,
                        linebreaks = linebreaks)
end

@inline function _render_text(::Val{:print}, str::AbstractString;
                              compact_printing::Bool = true,
                              isstring::Bool = false,
                              linebreaks::Bool = false)

    if linebreaks
        vstr = string.(split(str, '\n'))
    else
        vstr = [str]
    end

    # NOTE: Here we cannot use `escape_string(str)` because it also adds the
    # character `"` to the list of characters to be escaped.
    return [sprint(escape_string, s, "", sizehint = lastindex(s)) for s in vstr]
end

@inline function _render_text(::Val{:show}, v;
                              compact_printing::Bool = true,
                              linebreaks::Bool = false,
                              isstring::Bool = false)

    if v isa AbstractString
        aux  = linebreaks ? string.(split(v, '\n')) : [v]
        vstr = sprint.(show, aux; context = :compact => compact_printing)

        if !isstring
            for i = 1:length(vstr)
                aux_i   = first(vstr[i], length(vstr[i])-1)
                vstr[i] = last(aux_i, length(aux_i)-1)
            end
        end
    else
        str  = sprint(show, v; context = :compact => compact_printing)
        vstr = linebreaks ? string.(split(str, '\n')) : [str]
    end

    return vstr
end

"""
    _str_aligned(data::String, alignment::Symbol, field_size::Integer, lstr::Integer = -1)

This function returns the string `data` with alignment `alignment` in a field
with size `field_size`. `alignment` can be `:l` or `:L` for left alignment, `:c`
or `:C` for center alignment, or `:r` or `:R` for right alignment. It defaults
to `:r` if `alignment` is any other symbol.

This function also returns the new size of the aligned string.

If the string is larger than `field_size`, then it will be cropped and `⋯` will
be added as the last character.

The size of the string can be passed to `lstr` to save computational burden. If
`lstr = -1`, then the string length will be computed inside the function.

"""
function _str_aligned(data::String, alignment::Symbol,
                      field_size::Integer, lstr::Integer = -1)

    lstr < 0 && (lstr = textwidth(data))
    Δ = field_size - lstr

    # If the length is larger than the field, then we will crop the string.
    if Δ < 0
        data  = _crop_str(data, lstr + Δ - 1)
        data *= "…"

        # In this case, the string has the same size of the field.
        lstr = lstr + Δ - 1
        Δ = 0
    end

    # TODO: If Δ is 0, can we just return the string?

    if alignment == :l || alignment == :L
        return data * " "^Δ, lstr + Δ
    elseif alignment == :c || alignment == :C
        left  = div(Δ,2)
        right = Δ-left
        return " "^left * data * " "^right, lstr + Δ
    else
        return " "^Δ * data, lstr + Δ
    end
end

"""
    _str_autowrap(tokens_raw::Vector{String}, width::Int = 0)

Autowrap the tokens in `tokens_raw` considering a field width of `width`. It
returns a new vector with the new wrapped tokens.

"""
function _str_autowrap(tokens_raw::Vector{String}, width::Int = 0)
    # Check for errors.
    (width <= 0) && error("If `autowrap` is true, then the width must not be positive.")

    tokens = String[]

    for token in tokens_raw
        sub_tokens = String[]
        length_tok = length(token)

        # Get the list of valid indices to handle UTF-8 strings. In this
        # case, the n-th character of the string can be accessed by
        # `token[tok_ids[n]]`.
        tok_ids = collect(eachindex(token))

        if length_tok > width
            # First, let's analyze from the beginning of the token up to the
            # field width.
            #
            # k₀ is the character that will start the sub-token.
            # k₁ is the character that will end the sub-token.
            k₀ = 1
            k₁ = k₀ + width - 1

            while k₀ <= length_tok

                # Check if the remaining string fit in the available space.
                if k₁ == length_tok
                    push!(sub_tokens, token[tok_ids[k₀:k₁]])

                else
                    # If the remaining string does not fit into the
                    # available space, then we search for spaces to crop.
                    Δ = 0
                    for k = k₁:-1:k₀
                        if token[tok_ids[k]] == ' '
                            # If a space is found, then select `k₁` as this
                            # character and use `Δ` to remove it when
                            # printing, so that we hide the space.
                            k₁ = k
                            Δ  = 1

                            break
                        end
                    end

                    push!(sub_tokens, token[tok_ids[k₀:k₁-Δ]])
                end

                # Move to the next analysis window.
                k₀ = k₁+1
                k₁ = clamp(k₀ + width - 1, 0, length_tok)
            end
            push!(tokens, sub_tokens...)
        else
            push!(tokens, token)
        end
    end

    return tokens
end

################################################################################
#                             ANSI Escape Sequence
################################################################################

"""
    _get_composed_ansi_format(ansi::Vector{T}) where T<:AbstractString

Given a vector with a set of ANSI escape sequences, return a composed escape
sequence that leads to the same formatting.

!!! warning

    This function only works with the minimal set used by `Markdown` in
    `stdlib`.

"""
function _get_composed_ansi_format(ansi::Vector{T}) where T<:AbstractString
    bold = false
    underline = false
    color = 39

    for a in ansi
        length(a) ≤ 3  && continue
        a[1]   != '\e' && continue
        a[2]   != '['  && continue
        a[end] != 'm'  && continue

        code = parse(Int, a[3:end-1])

        if code == 0
            color = 39
            bold = false
            underline = false
        elseif code == 1
            bold = true
        elseif code == 4
            underline = true
        elseif code == 22
            bold  = false
            color = 39
        elseif code == 24
            underline = false
        elseif 30 ≤ code ≤ 37
            color = code
        elseif code == 39
            color = 39
        end
    end

    composed_ansi = ""

    bold        && (composed_ansi *= "\e[1m")
    underline   && (composed_ansi *= "\e[4m")
    color != 39 && (composed_ansi *= "\e[" * string(color) * "m")

    return composed_ansi
end

"""
    _reapply_ansi_format!(lines::Vector{T}) where T<:AbstractString

For each line in `lines`, reapply the ANSI format left by the previous line.

"""
function _reapply_ansi_format!(lines::Vector{T}) where T<:AbstractString
    length(lines) ≤ 1 && return nothing

    aux  = collect(eachmatch(r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])", lines[1]))
    ansi = map(x->x.match, aux)

    lines[1] *= "\e[0m"

    @inbounds @views for i = 2:length(lines)
        composed_ansi = _get_composed_ansi_format(ansi)

        # Find the first non-blank character.
        id = findfirst(x->x ≠ ' ', lines[i])

        if (id == nothing) || (id == 1)
            lines[i] = composed_ansi * lines[i]
        else
            lines[i] = lines[i][1:id-1] * composed_ansi * lines[i][id:end]
        end

        aux  = collect(eachmatch(r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])", lines[i]))
        ansi = map(x->x.match, aux)

        lines[i] *= "\e[0m"
    end

    return nothing
end


text:
# include("../src/CurricularAnalyticsSimulationLibrary.jl")
using CurricularAnalyticsSimulationLibrary
using Test

# Testing Students
student = Student(1, Dict())
@test student.id == 1
@test student.gpa == 0.0
@test student.total_credits == 0
@test student.total_points == 0
@test student.termcredits == 0
@test_throws MethodError Student()
@test_throws MethodError Student(1.4, Dict())
@test_throws MethodError Student(1, 8)

# Testing Course
c1 = Course("Course One", 3, Course[], Course[])
c2 = Course("Course Two", 3, Course[], Course[])
c3 = Course("Course Three", 3, [c1], Course[])
c4 = Course("Course Four", 3, [c1], Course[])
@test length(c1.postreqs) == 2
@test length(c3.prereqs) == 1
cruciality(c1)
cruciality(c2)
cruciality(c3)
cruciality(c4)
@test c1.delay == 2
@test c1.blocking == 2
@test c1.cruciality == 4
@test c2.cruciality == 1
@test c3.cruciality == 2
@test c4.cruciality == 2

# Testing Curriculum
curriculum = Curriculum("Test Curriculum", "./test_curriculum.json")
@test curriculum.numCourses == 4
@test curriculum.courses[1].cruciality == 4
@test curriculum.complexity == 9
@test curriculum.delay == 7
@test curriculum.blocking == 2

# Test Simulation with Passrate Model
# sim = Simulation(curriculum)
students = simpleStudents(10)
sim = simulate(curriculum, students, max_credits = 9, duration = 1)
@test sim.gradRate == 0.0
@test c3.enrolled == 0
@test c4.enrolled == 0

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
module SynthControl

using Dates, DataFrames, Optim, LinearAlgebra, RecipesBase, Parameters
import CSV

export TreatmentPanel, SynthControlModel, fit!, isfitted, load_brexit, load_germany, load_basque, load_smoking

abstract type SCM end

"""
    TreatmentPanel

Holds the data to be used in fitting synthetic control models. The object is constructed
from a `DataFrame` which contains all outcome and covariate data. The constructor requires
passing information on which column in the `DataFrame` holds the subject and time period
identifiers, as well as information on the timing of treatment. Treatments are specified
as a `Pair` of either a `String` or `Symbol` identifying the unit treated, and a `Date` or
`Int` value, identifying the treatment period. Where treatments have an end point, the
treatment period os specified as a `Tuple` of either `Date` or `Int`, indicating the first
and last period of treatment. Where individual units have multiple treatment periods,
these are specified as a `Pair{Union{String, Symbol}, Vector{Union{Int, Date}}}` of treatment
unit identifier and vector of timings. Finally, single continuous, single periodic, and
multiple period treatments can be generalized to multiple treated units.

The following table provides an overview of the types of treatment pattern supported:

|                 |  Only starting point        |   Start and end point                     |   Multiple start & end points                       |
|-----------      |------------------------     |-------------------------                  |------------------------------                       |
| **one unit**         |  Pair{String, Date}         |   Pair{String, Tuple{Date, Date}}         |  Pair{String}, Vector{Tuple{Date, Date}}}           |
| **multiple units**  |  Vector{Pair{String, Date}} |   Vector{Pair{String, Tuple{Date, Date}}} |  Vector{Pair{String}, Vector{Tuple{Date, Date}}}}   |

Currently, only single treatment unit and continuous treatment is supported.

"""

@with_kw struct TreatmentPanel
    data::DataFrame
    outcome::Union{Symbol, String}
    id_var::Union{Symbol, String}
    t_var::Union{Symbol, String}
    treatment::Union{Pair{Union{Symbol, String}, Union{Date, Int64}}}
    y_pre::Vector{Float64}
    y_post::Vector{Float64}
    x_pre::Matrix{Float64}
    x_post::Matrix{Float64}
    no_comps::Int64
    no_treatment_periods::Int64
    no_pretreatment_periods::Int64
    comp_labels::Vector{Union{Symbol, String}}
end

function check_tp(outcome, id_var, t_var, treatment, data)

    treatᵢ, treatₜ = treatment

    !isnothing(outcome) || error(ArgumentError(
        "Please specify outcome, the nameof the column in your dataset holding the "*
        "outcome variable of interest in your dataset."
        ))
    !isnothing(id_var) || error(ArgumentError(
            "Please specify id_var, the name of the column in your data set holding the "*
            "identifier of your units of observation (panel dimension)."
            ))
    !isnothing(t_var) || error(ArgumentError(
            "Please specify t_var, the name of the column in your dataset holding the "*
            "time dimension."
        ))
    !isnothing(treatment) || error(ArgumentError(
        "Please specify treatment, the identifier of the treated unit(s) in your "*
        "dataset and associated start(s) of treatment."
        ))
    in(string(outcome), names(data)) || error(DomainError(
        "The specified outcome column, $outcome, was not found in your data set. Please "*
        "double check the name of the outcome column."
        ))
    in(treatᵢ, data[!, id_var]) || error(DomainError(
        "The specified treatment unit, $(treatᵢ), was not found in the specified id "*
        "column, $id_var. Please ensure that treatment unit and id column are correctly "*
        "specified."
        ))
    in(treatₜ, data[!, t_var]) || error(DomainError(
        "The specified treatment period, $(treatₜ), was not found in the specified time "*
        "column, $t_var. Please double check treatment time and time column."
        ))
    treatₜ != maximum(data[!, t_var]) || error(DomainError(
        "The specified treatment period, $(treatₜ), is equal to the last observation period "*
        "of your data set - no estimation of treatment effects is possible. Please double "*
        "check the treatment period."
        ))
    treatₜ < maximum(data[!, t_var]) || error(DomainError(
        "The specified treatment period, $(treatₜ), occurs after the last observation in your "*
        "data set, $(maximum(data[!, t_var])) - no estimation of treatment effects is possible. "*
        "Please double check specification of treatment period."
        ))
    # To add - checks around balance of panel?
end

function TreatmentPanel(treatment::T, data;
                        outcome = nothing, id_var = nothing, t_var = nothing) where T <: Pair

    check_tp(outcome, id_var, t_var, treatment, data)

    treatᵢ, treatₜ = treatment

    y_pre = data[(data[!, id_var] .== treatᵢ) .&
               (data[!, t_var] .< treatₜ), outcome]
    xs_pre = data[(data[!, id_var] .!= treatᵢ) .&
               (data[!, t_var] .< treatₜ), outcome]

    y_post = data[(data[!, id_var] .== treatᵢ) .&
               (data[!, t_var] .>= treatₜ), outcome]
    xs_post = data[(data[!, id_var] .!= treatᵢ) .&
               (data[!, t_var] .>= treatₜ), outcome]

    no_treatment_periods = length(unique(data[data[!, t_var] .>= treatₜ, t_var]))
    no_pretreatment_periods = length(y_pre)
    no_comps = length(unique(data[data[!, id_var] .!= treatᵢ, id_var]))

    comp_labels = unique(data[data[!, id_var] .!= treatᵢ, id_var])

    TreatmentPanel(
        data,
        outcome,
        id_var,
        t_var,
        treatment,
        y_pre,
        y_post,
        reshape(xs_pre, no_comps, no_pretreatment_periods),
        reshape(xs_post, no_comps, no_treatment_periods),
        no_comps,
        no_treatment_periods,
        no_pretreatment_periods,
        comp_labels
    )
end

"""
    SynthControlModel

Takes in a `TreatmentPanel` and holds the results of calling the `fit!` method, optimal
weights for comparator units and the resulting predicted outcome of the synthetic control
unit as well as the estimated treatment effect from comparing synthetic control and observed
outcome of the treated unit.

"""

mutable struct SynthControlModel <: SCM
    treatment_panel::TreatmentPanel
    w::Vector{Float64} # Vector of weights
    ŷ::Vector{Float64} # Vector of post-treatment predicted outcomes for synthetic unit
    δ̂::Vector{Float64} # Vector of estimated treatmet impacts
    p_test_res::Array{Float64, 2} # Matrix of outcomes for untreated under placebo
end

function SynthControlModel(tp::TreatmentPanel)
    @unpack y_pre, y_post, no_comps, no_treatment_periods = tp

    ŷ = zeros(length(y_pre) + length(y_post))
    w = zeros(no_comps)
    δ̂ = zeros(no_treatment_periods)
    p_test_res = zeros(length(ŷ), no_comps)

    SynthControlModel(tp, w, ŷ, δ̂, p_test_res)
end

"""
isfitted(s::SynthControlModel)

Check whether the `SynthControlModel` object `s` has been fitted.

# Examples
```julia-repl
julia> df = load_brexit();

julia> s = SynthControlModel(df, :country, :dateid, :realgdp, 86, "United Kingdom");

julia> isfitted(s)
false

julia> fit!(s);

julia> isfitted(s)
true

```
"""
isfitted(s::SynthControlModel) = sum(s.w) ≈ 0.0 ? false : true

"""
    fit!(s::SynthControlModel; placebo_test = false)

Fit the `SynthControlModel` `s` by finding the weights that minimize the distance between
the pre-treatment outcomes for the observational unit of interest and the weighted average
pre-treatment outcomes for unweighted units.

If `placebo_test = true` is supplied, additional placebo tests will be performed by using
every non-treated unit in the data set as the treated unit in turn and estimating the
treatment impact on this unit. Results are stored in the `p_test_res` field and can be
used as the basis for inference.
"""
function fit!(s::SynthControlModel; placebo_test = false)

    @unpack y_pre, y_post, x_pre, x_post, no_comps, no_pretreatment_periods = s.treatment_panel

    # Objective function
    # Hat tip to Mathieu Tanneau who suggested switching to a quadratic formulation
    obj(w) = dot(y_pre .- vec(w'*x_pre), y_pre .- vec(w'*x_pre)) + 1e6*(1.0 - sum(w))^2

    # Initial condition and bounds
    initial = [1/no_comps for _ in 1:no_comps]
    lower = zeros(no_comps)
    upper = ones(no_comps)

    # Get weights through optimization
    res = optimize(obj, lower, upper, initial)
    s.w = res.minimizer

    # Get estimates
    s.ŷ = vec(s.w'*[x_pre x_post])
    s.δ̂ = ([y_pre; y_post] .- s.ŷ)[no_pretreatment_periods+1:end]

    if placebo_test
        placebos = unique(s.treatment_panel.data[(s.treatment_panel.data[!,s.pid] .!= s.treat_id), s.pid])
        p = deepcopy(s)

        for n ∈ 1:p.no_comps
            # change treatment ID
            p.treat_id = placebos[n]
            # change outcome data
            p.y = p.data[p.data[!, p.pid] .== p.treat_id, p.outcome]
            p.comps = collect(reshape(p.data[(p.data[!,p.pid] .!= p.treat_id) .& (p.data[!,p.tid].<p.treat_t), p.outcome],
                    p.no_comps, p.no_pretreat_t)')
            fit!(p)
            s.p_test_res[:, n] = p.y .- p.ŷ
        end
    end
    return s
end

# Define plotting recipe
@recipe function f(s::SynthControlModel; kind = "overall")
    legend --> :topleft

    if kind == "overall"

        @series begin
            label --> s.treatment_panel.treatment[1]
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), [s.treatment_panel.y_pre; s.treatment_panel.y_post]
        end

        @series begin
            label --> "Control"
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), s.ŷ
        end

        @series begin
            label --> "Impact"
            seriestype := :bar
            seriesalpha := 0.5
            linecolor := "white"
            seriescolor := "darkgreen"
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var]))[findfirst(x -> x >= s.treatment_panel.treatment[2],
                                    sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var]))):end], s.δ̂
        end

    elseif kind == "diffplot"
        @series begin
            label --> s.treatment_panel.treatment[1]
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), [s.treatment_panel.y_pre; s.treatment_panel.y_post] .- s.ŷ
        end

        @series begin
            label --> ""
            seriestype := :scatter
            seriescolor := 1
            markerstrokecolor := "white"
            sort!(unique(s.treatment_panel.data[!, s.treatment_panel.t_var])), [s.treatment_panel.y_pre; s.treatment_panel.y_post] .- s.ŷ
        end
    end
        @series begin
            label --> "Intervention"
            seriestype := :vline
            linestyle := :dash
            seriescolor := "black"
            xguide := "Time"
            yguide := "Outcome"
            [s.treatment_panel.treatment[2]]
        end

     nothing
end

# Data loading for example
"""
    load_brexit()

Returns a `DataFrame` with quarterly GDP data on 30 OECD countries, borrowed from the
analysis of the effect of Brexit on UK GDP undertaken by the [Centre for European Reform](
https://www.cer.eu/insights/cost-brexit-june-2019)
"""
function load_brexit()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","brexit.csv"), DataFrame)
end

function load_germany()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","germany.csv"), DataFrame)
end

function load_basque()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","basque.csv"), DataFrame)
end

function load_smoking()
    CSV.read(joinpath(dirname(@__FILE__),"..","data","smoking.csv"), DataFrame)
end

# Pretty printing
import Base.show

function show(io::IO, ::MIME"text/plain", s::SynthControlModel)

    println(io, "\nSynthetic Control Model")
    println(io, "\tOutcome variable: ", s.treatment_panel.outcome)
    println(io, "\tTime dimension: ",string(s.treatment_panel.t_var)," with ",
            string(length(unique(s.treatment_panel.data[!,s.treatment_panel.t_var]))), " unique values")
    println(io, "\tTreatment period: ", string(s.treatment_panel.treatment[2]))
    println(io, "\tID variable: ",string(s.treatment_panel.id_var), " with ",
            string(length(unique(s.treatment_panel.data[!,s.treatment_panel.id_var]))), " unique values")
    println(io, "\tTreatment ID: ",string(s.treatment_panel.treatment[1]))

    if isfitted(s)
      println(io, "\tModel is fitted")
      println(io, "\tImpact estimates: ",round.(s.δ̂, digits=3))
    else
      println(io, "\n\tModel is not fitted")
    end
end

end # module


text:
using CellularAutomata
using Tk
importall Cairo

function sim_life(init::Array{Int,2},B, S)
    m, n = size(init)
    w = Window("CA", n, m)
    c = Canvas(w)
    pack(c)

#    win = Toplevel("Cellular Automata", 600, 300)
#    f = Frame(win)
#    pack(f, expand = true, fill = "both")
#    c = Canvas(f)
#    pack(c, expand = true, fill = "both")

    done = false
    c.mouse.button1press = (c, x, y)->(done=true)
    cr = getgc(c)

    last = time(); f = 1
    while !done
        img = convert(Array{Uint32, 2}, init) .* 0x00ffffff
        set_source_surface(cr, CairoRGBSurface(img), 0, 0)
        paint(cr)
        reveal(c)

        t = time()
        if (t-last) > 2
            println("$(f/(t-last)) FPS")
            last = t; f = 0
        end
        M = CA2d(B, S, init, 2).cells
        init = int(M[:, :, 2])
        f += 1
        sleep(0.01)
    end
end

#turbine pattern
#init = int(zeros(20,20))
#init[4, 4:9] = 1
#init[5, 4:9] = 1

#init[4:9, 11] = 1
#init[4:9, 12] = 1

#init[11, 7:12] = 1
#init[12, 7:12] = 1

#init[7:12, 4] = 1
#init[7:12, 5] = 1

init = int(zeros(600,600))
init[250:350, 250:350] = int(rand(101,101))
#init[295:305, 295:305] = int(rand(11,11))

#init = int(rand(1000, 1000))

#game of life
sim_life(init,[3], [2,3])

#morley
#sim_life(init,[3,6,8],[2,4,5])

#2x2
#sim_life(init,[3,6],[1,2,5])

#replicator
#sim_life(init,[1,3,5,7], [1,3,5,7])

#dimoeba
#sim_life(init,[3,5,6,7,8], [5,6,7,8])

#no death
#sim_life(init,[3], [0,1,2,3,4,5,6,7,8])

#34 life
#sim_life(init,[3,4],[3,4])


text:
const rinex3_lht = string(bin_dir, "/", "rinex3_lht" ,
                            isapple ? "_macox" :
                                islinux ? "_linux" :
                                    throw("no rinex3_lht in windows yet"),
                    )


text:
using Test
using Logging
using Statistics
using Printf
using JLD2
using CUDA

using Oceananigans
using Oceananigans.Architectures
using JULES

Logging.global_logger(OceananigansLogger())

         Archs = [CPU]
@hascuda Archs = [GPU]

CUDA.allowscalar(true)

@testset "JULES" begin
    include("test_models.jl")
    include("test_lazy_fields.jl")
    include("test_time_stepping.jl")
    include("test_regression.jl")
end


text:
"Generate default set of parameters for Radiative Transfer calculations (from ModelParameters/)"
default_parameters() = parameters_from_yaml(joinpath(dirname(pathof(RadiativeTransfer)), "vSmartMOM", "ModelParameters", "DefaultParameters.yaml"))

# Check that a field exists in yaml file
function check_yaml_field(dict::Dict{Any, Any}, full_keys::Array{String}, curr_keys::Array{String}, final_type::Type, valid_options::Array{String})

    # Should have at least one key
    @assert length(curr_keys) >= 1

    # Only one key (that should exist)
    if length(curr_keys) == 1
        @assert curr_keys[1] in keys(dict) "Missing key in parameters yaml: $(join(full_keys, '/'))"
        @assert dict[curr_keys[1]] isa final_type "Improper type for $(join(full_keys, '/')); must be a $(final_type)"
        if length(valid_options) > 0
            @assert dict[curr_keys[1]] in valid_options "Field $(join(full_keys, '/')) must be one of $(valid_options)"
        end
        return true

    # Recursive case, check sub-dictionary
    else
        @assert curr_keys[1] in keys(dict) "Missing key in parameters yaml: $(join(full_keys, '/'))"
        return check_yaml_field(dict[curr_keys[1]], full_keys, curr_keys[2:end], final_type, valid_options)
    end
end

# Given a list of aerosol parameter dictionaries, validate all of them
function validate_aerosols(aerosols::Array{Dict{Any, Any}})
    fields = [(["τ_ref"], Real),
              (["μ"], Real),
              (["σ"], Real),
              (["nᵣ"], Real),
              (["nᵢ"], Real),
              (["p₀"], Real),
              (["p₀"], Real)]

    for aerosol in aerosols
        for i in 1:length(fields)
            key, elty = fields[i][1:2]
            @assert check_yaml_field(aerosol, key, key, elty, Array{String}([]))
        end
    end
end

"Given a parameter dictionary from a YAML file, validate the dictionary"
function validate_yaml_parameters(params)

    # All the fields that *should* exist and their types
    fields = [(["absorption", "broadening"], String, ["Voigt()", "Doppler()", "Lorentz()"]),
              (["absorption", "CEF"], String, ["HumlicekWeidemann32SDErrorFunction()"]),
              (["absorption", "wing_cutoff"], Integer),
              (["absorption", "spec_bands"], Array{String}),
              (["absorption", "vmr"], Real),
              (["absorption", "molecules"], Array{String}),

              (["scattering", "aerosols"], Array{Dict{Any,Any}}),
              (["scattering", "r_max"], Real),
              (["scattering", "nquad_radius"], Real),
              (["scattering", "λ"], Array{<:Real}),
              (["scattering", "λ_ref"], Real),
              (["scattering", "depol"], Real),
              (["scattering", "polarization_type"], String, ["Stokes_I()", "Stokes_IQ()", "Stokes_IQU()", "Stokes_IQUV()"]),
              (["scattering", "decomp_type"], String, ["NAI2()", "PCW()"]),

              (["geometry", "vza"], Array{<:Real}),
              (["geometry", "vaz"], Array{<:Real}),
              (["geometry", "sza"], Real),
              (["geometry", "obs_alt"], Real),

              (["atmospheric_profile", "file"], String),
              (["atmospheric_profile", "profile_reduction"], Integer),

              (["surface", "BRDFs"], Array{String}),

              (["truncation_type", "l_trunc"], Integer),
              (["truncation_type", "Δ_angle"], Real),

              (["vSmartMOM", "quadrature_type"], String),
              (["vSmartMOM", "max_m"], Integer),
              (["vSmartMOM", "architecture"], String),
              (["vSmartMOM", "float_type"], String)
              ]

    # Check that every field exists and is correctly typed
    for i in 1:length(fields)
        key, elty = fields[i][1:2]
        valid_options = length(fields[i]) == 3 ? fields[i][3] : Array{String}([])
        @assert check_yaml_field(params, key, key, elty, valid_options)
    end

    # Check the aerosols separately
    validate_aerosols(params["scattering"]["aerosols"])
end

"Given a path to a YAML parameters file, load all the parameters into a vSmartMOM_Parameters struct"
function parameters_from_yaml(file_path)

    # Load the YAML file
    params_dict = YAML.load_file(file_path)

    # Validate to make sure it has all necessary fields
    validate_yaml_parameters(params_dict)

    # Evaluate certain fields to get their proper types
    broadening_function = eval(Meta.parse(params_dict["absorption"]["broadening"]))
    CEF = eval(Meta.parse(params_dict["absorption"]["CEF"]))
    polarization_type = eval(Meta.parse(params_dict["scattering"]["polarization_type"]))
    decomp_type = eval(Meta.parse(params_dict["scattering"]["decomp_type"]))
    quadrature_type = eval(Meta.parse(params_dict["vSmartMOM"]["quadrature_type"]))
    architecture = eval(Meta.parse(params_dict["vSmartMOM"]["architecture"]))
    FT = eval(Meta.parse(params_dict["vSmartMOM"]["float_type"]))
    BRDF_per_band = map(x -> eval(Meta.parse(x)), params_dict["surface"]["BRDFs"])
    lengths = convert.(Integer, map(x -> length(collect(eval(Meta.parse(x)))), params_dict["absorption"]["spec_bands"]))

    # Evaluate profile path if it's not a plain string (need to joinpath)
    if occursin("joinpath", params_dict["atmospheric_profile"]["file"])
        profile_path = eval(Meta.parse(params_dict["atmospheric_profile"]["file"]))
    else
        profile_path = params_dict["atmospheric_profile"]["file"]
    end

    # If the specified atmospheric profile is YAML, then set time_index/lat/lon to nothing
    if endswith(profile_path, ".yaml")
        time_index = nothing
        lat = nothing
        lon = nothing

    # Otherwise, store them
    else
        time_index = params_dict["atmospheric_profile"]["time_index"]
        lat = params_dict["atmospheric_profile"]["lat"]
        lon = params_dict["atmospheric_profile"]["lon"]
    end

    return vSmartMOM_Parameters(FT,
                                convert.(FT, params_dict["scattering"]["λ"]),
                                FT(params_dict["scattering"]["λ_ref"]),
                                FT(params_dict["scattering"]["depol"]),
                                params_dict["truncation_type"]["l_trunc"],
                                FT(params_dict["truncation_type"]["Δ_angle"]),
                                params_dict["vSmartMOM"]["max_m"],
                                polarization_type,
                                FT(params_dict["geometry"]["obs_alt"]),
                                FT(params_dict["geometry"]["sza"]),
                                convert.(FT, params_dict["geometry"]["vza"]),
                                convert.(FT, params_dict["geometry"]["vaz"]),
                                length(params_dict["scattering"]["aerosols"]),
                                convert.(FT, map(x -> x["τ_ref"], params_dict["scattering"]["aerosols"])),
                                convert.(FT, map(x -> x["μ"], params_dict["scattering"]["aerosols"])),
                                convert.(FT, map(x -> x["σ"], params_dict["scattering"]["aerosols"])),
                                FT(params_dict["scattering"]["r_max"]),
                                params_dict["scattering"]["nquad_radius"],
                                convert.(FT, map(x -> x["nᵣ"], params_dict["scattering"]["aerosols"])),
                                convert.(FT, map(x -> x["nᵢ"], params_dict["scattering"]["aerosols"])),
                                convert.(FT, map(x -> x["p₀"], params_dict["scattering"]["aerosols"])),
                                convert.(FT, map(x -> x["σp"], params_dict["scattering"]["aerosols"])),
                                profile_path,
                                time_index,
                                lat,
                                lon,
                                params_dict["atmospheric_profile"]["profile_reduction"],
                                convert.(FT, map(x -> first(collect(eval(Meta.parse(x)))), params_dict["absorption"]["spec_bands"])),
                                convert.(FT, map(x -> last(collect(eval(Meta.parse(x)))), params_dict["absorption"]["spec_bands"])),
                                convert(Array{Integer}, convert.(Integer, map(x -> length(collect(eval(Meta.parse(x)))), params_dict["absorption"]["spec_bands"]))),
                                broadening_function,
                                params_dict["absorption"]["wing_cutoff"],
                                CEF,
                                FT(params_dict["absorption"]["vmr"]),
                                BRDF_per_band,
                                decomp_type,
                                quadrature_type,
                                architecture)
end

"Take the parameters specified in the vSmartMOM_Parameters struct, and calculate derived attributes into a vSmartMOM_Model"
function model_from_parameters(params::vSmartMOM_Parameters)

    @unpack λ_band = params
    nBands = length(λ_band)

    # Create observation geometry
    obs_geom = ObsGeometry(params.obs_alt, params.sza, params.vza, params.vaz)

    # Create truncation type
    truncation_type = Scattering.δBGE{params.float_type}(params.l_trunc, params.Δ_angle)

    # Set quadrature points
    quad_points = rt_set_streams(params.quadrature_type, params.l_trunc, obs_geom, params.polarization_type, array_type(params.architecture))

    # Read profile (and generate dry/wet VCDs per layer)
    if isnothing(params.timeIndex)
        profile_hr = vSmartMOM.read_atmos_profile(params.file);
    else
        profile_hr = vSmartMOM.read_atmos_profile(params.file, params.lat, params.lon, params.timeIndex);
    end

    # Reduce the profile to the number of target layers
    profile = vSmartMOM.reduce_profile(params.profile_reduction_n, profile_hr);

    #Rayleigh
    greek_rayleigh = Scattering.get_greek_rayleigh(params.depol)
    τRayl = [zeros(params.float_type, length(profile.p)) for i=1:nBands];

    # Create empty array (can be changed later)
    # spec_grid[iBand][iSpec]
    spec_grid = [zeros(params.float_type,n) for n in params.spec_grid_n];

    # τ_abs[iBand][iSpec,iZ]
    τ_abs     = [zeros(params.float_type,n, length(profile.p)) for n in params.spec_grid_n]

    # Loop over all bands:
    for ib=1:length(λ_band)

        # Compute Rayleigh properties per layer for `ib` band center
        τRayl[ib]   = getRayleighLayerOptProp(profile.psurf / 100, λ_band[ib], params.depol, profile.vcd_dry);

        # Define spectral grid per band
        spec_grid[ib] = range(params.spec_grid_start[ib], params.spec_grid_end[ib], length=params.spec_grid_n[ib]);

        # This has to be more flexible in the future, so that we can have different trace gases per band, Newton's work
        # can be checked for that. (ideally a dictionary that defines gases to be used and profiles per band)
        hitran_data = read_hitran(artifact("O2"), iso=1)
        absorption_model = make_hitran_model(hitran_data,
                                         params.broadening_function,
                                         wing_cutoff = params.wing_cutoff,
                                         CEF = params.CEF,
                                         architecture = params.architecture,
                                         vmr = params.vmr)

        compute_absorption_profile!(τ_abs[ib], absorption_model, range(params.spec_grid_start[ib], params.spec_grid_end[ib], length=params.spec_grid_n[ib]), profile);
    end

    # aerosol_optics[iBand][iAer]
    aerosol_optics = [Array{AerosolOptics}(undef, (params.nAer)) for i=1:nBands];

    FT2 = typeof(params.τAer_ref[1])

    # τAer[iBand][iAer,iProfile]
    τAer = [zeros(FT2, params.nAer, length(profile.p)) for i=1:nBands];

    # Loop over aerosol type
    for iaer=1:params.nAer

        # Create Aerosol size distribution for each aerosol species
        size_distribution = LogNormal(log(params.μ[iaer]), log(params.σ[iaer]))

        # Create a univariate aerosol distribution
        aerosol = make_univariate_aerosol(size_distribution, params.r_max, params.nquad_radius, params.nᵣ[iaer], params.nᵢ[iaer]) #Suniti: why is the refractive index needed here?

        # Create the aerosol extinction cross-section at the reference wavelength:
        mie_model      = make_mie_model(params.decomp_type, aerosol, params.λ_ref, params.polarization_type, truncation_type)
        k_ref          = compute_ref_aerosol_extinction(mie_model, params.float_type)

        # Loop over bands
        for ib=1:nBands

            # Create the aerosols:
            mie_model      = make_mie_model(params.decomp_type, aerosol, params.λ_band[ib], params.polarization_type, truncation_type)

            # Compute raw (not truncated) aerosol optical properties (not needed in RT eventually)
            aerosol_optics_raw = compute_aerosol_optical_properties(mie_model, FT2);

            # Compute truncated aerosol optical properties (phase function and fᵗ), consistent with Ltrunc:
            @show iaer,ib
            aerosol_optics[ib][iaer] = Scattering.truncate_phase(truncation_type, aerosol_optics_raw; reportFit=false)

            # Compute nAer aerosol optical thickness profiles
            τAer[ib][iaer,:] = params.τAer_ref[iaer] * (aerosol_optics[ib][iaer].k/k_ref) * vSmartMOM.getAerosolLayerOptProp(1.0, params.p₀[iaer], params.σp[iaer], profile.p)

        end
    end

    # Just return band 1 for now:
    iBand = 1;
    return vSmartMOM_Model(params, aerosol_optics[iBand],  greek_rayleigh, quad_points, τ_abs[iBand], τRayl[iBand],τAer[iBand], obs_geom, profile, params.brdf[iBand])#, BRDF_type, a_surf[iBand])

end


text:
struct CUDACompilerParams <: AbstractCompilerParams end

CUDACompilerJob = CompilerJob{PTXCompilerTarget,CUDACompilerParams}

GPUCompiler.runtime_module(::CUDACompilerJob) = CUDA

# filter out functions from libdevice and cudadevrt
GPUCompiler.isintrinsic(job::CUDACompilerJob, fn::String) =
    invoke(GPUCompiler.isintrinsic,
           Tuple{CompilerJob{PTXCompilerTarget}, typeof(fn)},
           job, fn) ||
    fn == "__nvvm_reflect" || startswith(fn, "cuda")

function GPUCompiler.finish_module!(job::CUDACompilerJob, mod::LLVM.Module)
    invoke(GPUCompiler.finish_module!,
           Tuple{CompilerJob{PTXCompilerTarget}, typeof(mod)},
           job, mod)
    emit_exception_flag!(mod)
end

function GPUCompiler.link_libraries!(job::CUDACompilerJob, mod::LLVM.Module,
                                     undefined_fns::Vector{String})
    invoke(GPUCompiler.link_libraries!,
           Tuple{CompilerJob{PTXCompilerTarget}, typeof(mod), typeof(undefined_fns)},
           job, mod, undefined_fns)
    link_libdevice!(mod, job.target.cap, undefined_fns)
end


text:
module Othermethods
    using Random,Distributions
    using LinearAlgebra

    using FastGaussQuadrature

    #using DoubleExponentialFormulas
    #using DifferentialEquations
    #using QuadGK

    import ..CGmethods:bicg,cg,shiftedcg,reducedshiftedcg
    import ..Diracoperators:DdagD_operator,DdagD_Staggered_operator,DdagDND_Staggered_operator
    import ..Fermionfields:gauss_distribution_fermi_Z2!
    import ..Fermionfields

    function Tn(n,x)
        return cos(n*acos(x))
    end

    function fitfunc(an,x)
        val = 0
        for n=1:2:length(an)-1
            val += an[n+1]*Tn(n,x)
        end
        return val
    end

    function integcheb(an,a,b)
        wa = 0
        for k=1:length(an)-2
            wa += (Tn(k,b)-Tn(k,a))*(an[k-1+1]-an[k+1+1])/k
        end
        wa /= 2
        return wa
    end

    function calc_vecfunc_cg(func,vec_x,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
        DN = Adiag .+ ANmat

        numdata = length(vec_x)
        data = zeros(Float64,numdata)

        #Aϕ2 = zero(Aϕ)
        #ϕtemp2 = zero(ϕtemp1)

        #x = zero(ϕ)
        #Nx = zero(x)

        #=
        vec_Ax = Array{typeof(x),1}(undef,numdata)
        for l=1:numdata
            vec_Ax[l] = zero(x)
        end
        =#
        vec_tinv = Adiag[1,1]*(1 ./ vec_x  .- 1) #(D/t + A) = ((D(1/t - 1)I + D I + A)

        for k=1:M
            for i=1:N
                ϕ[i] = rand([-1,1])
                #ϕ[i] = randn()
                #ϕ[i] = rand([-1,1,-im,im])
                #r = rand()*2π
                #ϕ[i] = cos(r)+im*sin(r)
            end
            mul!(Aϕ,ANmat,ϕ)

            #cg(ϕtemp1,DtN,ϕ,eps = 1e-15,maxsteps= 3000)
            #shiftedcg(vec_x,vec_β,x,WdagW,univ.η,eps = univ.fparam.eps,maxsteps= univ.fparam.MaxCGstep)
            #@time shiftedcg(vec_Ax,vec_tinv,x,DN,ϕ,eps = 1e-15,maxsteps= 3000)
            ϕtemp1 .= 0
            θ = reducedshiftedcg(Aϕ,vec_tinv,ϕtemp1,DN,ϕ,eps = 1e-15,maxsteps= 3000)
            θ ./= (M .* vec_x)
            data .+= real.(θ)

            #=

            for l=1:numdata
                t = vec_x[l]
                DtN = Adiag .+ t*ANmat
                cg(ϕtemp2,DtN,ϕ,eps = 1e-15,maxsteps= 3000)
                mul!(Aϕ2,ANmat,ϕtemp2)
                val = real(ϕ ⋅ Aϕ2 /M)
                println("$t data,val ",real(θ[l]),"\t",val)
            end
            =#

            #=
            for l=1:numdata
                #mul!(Aϕ,ANmat,vec_Ax[l])
                valshift = real(θ[l]) #real(Nx ⋅ vec_Ax[l] /(M*vec_x[l]))
                #println("val1 ",valshift,"\t theta ",real(θ[l]))

                #=
                mul!(Aϕ,ANmat,vec_Ax[l])
                valshift = real(ϕ ⋅ Aϕ /(M*vec_x[l]))
                println("val2 ",valshift)
                =#

                data[l] += valshift

                t = vec_x[l]
                DtN = Adiag .+ t*ANmat
                cg(ϕtemp1,DtN,ϕ,eps = 1e-15,maxsteps= 3000)
                mul!(Aϕ,ANmat,ϕtemp1)
                val = real(ϕ ⋅ Aϕ /M)
                println("$t data,val ",valshift,"\t",val)
                #val += ϕ ⋅ Aϕ /M
            end
            ϕtemp1 .= 0
            =#
            #=
            for l=1:numdata
                vec_Ax[l] .= 0
            end
            =#

        end
        #exit()

        return data

        exit()

        for (i,x) in enumerate(vec_x)
            #val1 = func(x)
            #println("1 ",val1)
            val1,val2 = calc_fa(x,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
            #println("2 ",val1)
            #exit()
            println("1 ",data[i])
            data[i] = val1
            println("2 ",data[i])
        end
        return data
    end

    function calc_vecfunc(func,vec_x)
        N = length(vec_x)
        data = zeros(Float64,N)
        for (i,x) in enumerate(vec_x)
            data[i] = func(x)
        end
        return data
    end

    function calc_fa(t,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
        DtN = Adiag .+ t*ANmat

        if t < 1e-15
            DtN2 = Adiag
        else
            DtN2 = Adiag ./t .+ ANmat
        end
        #count += 1
        #println(count)
        val = 0
        val2 = 0
        #println("random")
        for k=1:M
            #gauss_distribution_fermi_Z2!(ϕ)

            for i=1:N
                ϕ[i] = rand([-1,1])
            end
            #rand!(d,ϕ)
            #Dinvϕ = DtN \ ϕ
            cg(ϕtemp1,DtN,ϕ,eps = 1e-15,maxsteps= 3000)
            mul!(Aϕ,ANmat,ϕtemp1)
            val += ϕ ⋅ Aϕ /M

            if t < 1e-15
                ϕtemp1 = ϕ/Adiag[1,1]
                mul!(Aϕ,ANmat,ϕtemp1)
                val2 += ϕ ⋅ Aϕ /M
            else
                cg(ϕtemp1,DtN2,ϕ,eps = 1e-15,maxsteps= 3000)
                mul!(Aϕ,ANmat,ϕtemp1)
                val2 += ϕ ⋅ Aϕ /(M*t)
            end

            #exit()

            #val += ϕ'*Aϕ/M#ANmat*Dinvϕ/M
        end
        println("ts: $t ",real(val),"\t",real(val2))
        return real(val),real(val2)
    end

    #calc_f,m,nc,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1
    function chebyshevfit(func,m::T,nc,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1) where T <: Int
        Δt = (π/2)/2m
        vec_x = zeros(Float64,2m+1)

        for i=0:2m
            t = i*Δt
            x = cos(t)
            vec_x[i+1] = x
        end
        @time data = calc_vecfunc_cg(func,vec_x,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)

        #fp = open("hiseki.dat","w")
        #close(fp)

        #=
        data1 = calc_vecfunc(func,vec_x)
        for i=1:length(data)
            println(data[i],"\t",data1[i])
        end
        exit()
        =#
        an = zeros(Float64,nc+1)

        for n=1:2:nc

            for i=1:m
                k = 2i-2
                k1 = 2i-1
                k2 = 2i
                t = k*Δt
                t1 = k1*Δt
                t2 = k2*Δt

                val = (2/π)*2*cos(n*t)*data[k+1]
                val1 = (2/π)*2*cos(n*t1)*data[k1+1]
                val2 = (2/π)*2*cos(n*t2)*data[k2+1]

                an[n+1] += Δt*(val + 4*val1 + val2)/3
            end
        end

        return an

    end

    function chebyshevfit(func,m,nc)
        Δt = (π/2)/2m
        vec_x = zeros(Float64,2m+1)

        for i=0:2m
            t = i*Δt
            x = cos(t)
            vec_x[i+1] = x
        end
        data = calc_vecfunc(vec_x,func)
        an = zeros(Float64,nc+1)

        for n=1:2:nc

            for i=1:m
                k = 2i-2
                k1 = 2i-1
                k2 = 2i
                t = k*Δt
                t1 = k1*Δt
                t2 = k2*Δt

                val = (2/π)*2*cos(n*t)*data[k+1]
                val1 = (2/π)*2*cos(n*t1)*data[k1+1]
                val2 = (2/π)*2*cos(n*t2)*data[k2+1]

                an[n+1] += Δt*(val + 4*val1 + val2)/3
            end
        end

        return an

        data = zeros(Float64,m)
        data1 = zeros(Float64,m)
        data2 = zeros(Float64,m)
        count = 0
        for i=1:m
            t = (2i-2)*Δt
            x = cos(t)

            t1 = (2i-1)*Δt
            x1 = cos(t1)

            t2 = 2i*Δt
            x2 = cos(t2)

            val = func(x)
            val1 = func(x1)
            val2 = func(x2)

            data[i] = val
            data1[i] = val1
            data2[i] = val2
        end
        an = zeros(Float64,nc+1)

        for n=1:2:nc

            for i=1:m
                t = (2i-2)*Δt
                t1 = (2i-1)*Δt
                t2 = 2i*Δt

                val = (2/π)*2*cos(n*t)*data[i]
                val1 = (2/π)*2*cos(n*t1)*data1[i]
                val2 = (2/π)*2*cos(n*t2)*data2[i]

                an[n+1] += Δt*(val + 4*val1 + val2)/3
            end
        end

        return an

        fp = open("hiseki.dat","w")

        for i=1:m
            t = (2i-2)*Δt
            x = cos(t)

            t1 = (2i-1)*Δt
            x1 = cos(t1)

            t2 = 2i*Δt
            x2 = cos(t2)

            val = data[i]
            val1 = data1[i]
            val2 = data2[i]

            valx = 0
            for n=1:2:nc
                valx += an[n+1]*Tn(n,x)
            end

            println("$x $val $valx")
            println(fp,"$x $val $valx")
        end
        close(fp)

        wa2 = 0
        Δx = 1/2m
        for i=1:m
            x = (2i-2)*Δx
            x1 = (2i-1)*Δx
            x2 = 2i*Δx

            val = fitfunc(an,x)
            val1 = fitfunc(an,x1)
            val2 = fitfunc(an,x2)

            wa2 += Δx*(val + 4*val1 + val2)/3
        end

        println(wa2)
        return wa2
        exit()

        return wa

    end

    function tdlogdet(A,M,m,tempvecs;filename=nothing,nc = 20,nonc = false)
        N,_ = size(A)
        Adiag = zero(A)
        d = Normal(0,1)
        ϕ = tempvecs[1]
        Aϕ = tempvecs[2]
        ϕtemp1 = tempvecs[3]

        #=
        for i=1:M
            rand!(d,ϕ)
            mul!(Aϕ,A,ϕ)
            for k=1:N
                Adiag[k,k] += real(conj(ϕ[k])*Aϕ[k]/M)
            end
        end
        =#

        for k=1:N
            Adiag[k,k] = A[1,1]
        end

        #for k=1:N
        #    println("$k $(A[k,k]) $(Adiag[k,k])")
        #end

        ANmat = A .- Adiag
        Δ0 = 0
        for k=1:N
            Δ0 += real(log(Adiag[k,k]))
        end
        Δx = 1/2m

        #@time ldet = real(logdet(Matrix(A)))
        #println("log det A(exact): $ldet")

        count = 0

        function calc_f(t)
            DtN = Adiag .+ t*ANmat

            if t < 1e-15
                DtN2 = Adiag
            else
                DtN2 = Adiag ./t .+ ANmat
            end
            #count += 1
            #println(count)
            val = 0
            val2 = 0
            #println("random")
            for k=1:M
                #gauss_distribution_fermi_Z2!(ϕ)

                for i=1:N
                    ϕ[i] = rand([-1,1])
                end
                #rand!(d,ϕ)
                #Dinvϕ = DtN \ ϕ
                cg(ϕtemp1,DtN,ϕ,eps = 1e-15,maxsteps= 3000)
                mul!(Aϕ,ANmat,ϕtemp1)
                val += ϕ ⋅ Aϕ /M

                if t < 1e-15
                    ϕtemp1 = ϕ/Adiag[1,1]
                    mul!(Aϕ,ANmat,ϕtemp1)
                    val2 += ϕ ⋅ Aϕ /M
                else
                    cg(ϕtemp1,DtN2,ϕ,eps = 1e-15,maxsteps= 3000)
                    mul!(Aϕ,ANmat,ϕtemp1)
                    val2 += ϕ ⋅ Aϕ /(M*t)
                end

                #exit()

                #val += ϕ'*Aϕ/M#ANmat*Dinvϕ/M
            end
            println("t: $t ",real(val),"\t",real(val2))
            return real(val)
        end

        function calc_finv(t)
            if t == 0
                DtN = Adiag
            else
                DtN = Adiag ./t .+ ANmat
            end
            #count += 1
            #println(count)
            val = 0
            #println("random")
            for k=1:M
                #gauss_distribution_fermi_Z2!(ϕ)

                for i=1:N
                    ϕ[i] = rand([-1,1])
                end
                #rand!(d,ϕ)
                #Dinvϕ = DtN \ ϕ
                if t == 0
                    ϕtemp1 = ϕ/Adiag[1,1]
                    mul!(Aϕ,ANmat,ϕtemp1)
                    val += ϕ ⋅ Aϕ /M
                else
                    cg(ϕtemp1,DtN,ϕ,eps = 1e-15,maxsteps= 3000)
                    mul!(Aϕ,ANmat,ϕtemp1)
                    val += ϕ ⋅ Aϕ /(M*t)
                end

                #exit()

                #val += ϕ'*Aϕ/M#ANmat*Dinvϕ/M
            end
            println("t: $t ",val)
            return real(val)
        end

        Δx = 1/2m
        Δ = Δ0
        Δc = Δ0

        #I,E = quadgk(calc_f, 0,1,maxevals=3m)
        #println(I)
        #println(E)
        #exit()

        nodes, weights = gausslegendre( 1000 )
        vec_x = nodes ./ 2 .+ 1/2

        @time data = calc_vecfunc_cg(calc_f,vec_x,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
        Δ += dot(weights,data)/2

        return Δ

        #nc = 20

        if nonc == false
            #@time ldet = real(logdet(Matrix(A)))
            ldet = 0
            println("log det A(exact): $ldet")

            nodes, weights = gausslegendre( 1000 )
            vec_x = nodes ./ 2 .+ 1/2
            #=
            @time data = calc_vecfunc_cg(calc_f,vec_x,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
            Δ += dot(weights,data)/2
            println(" Δ ",Δ)
            exit()
            =#

            #=
            fpk = open("kdep_L8_nc10_complex.dat","w")
            kmax=32
            for k=1:kmax
                M2 = 10
                nc2 = 10
                m2 = 1000
                @time an = chebyshevfit(calc_f,m2,nc2,M2,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
                Δ = integcheb(an,0,1) + Δ0

                #@time data = calc_vecfunc_cg(calc_f,vec_x,M2,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
                #Δ = dot(weights,data)/2 + Δ0
                println("$k $Δ $ldet")
                println(fpk,"$k $Δ $ldet")
            end
            close(fpk)
            exit()
            =#

            fp4 = open("MGL10_complex_k4.dat","w")
            kmax=4
            for M2 in [16,32,64,128,256,512,1024]
                Δ = Δ0
                for k=1:kmax
                    @time data = calc_vecfunc_cg(calc_f,vec_x,M2,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
                    Δ += (dot(weights,data)/2)/kmax
                end
                #Δ = dot(weights,data)/2 + Δ0
                println("M = $M2, Δ = $Δ, exact: $ldet ")
                println(fp4,"$M2 $Δ $ldet ")
            end
            close(fp4)
            exit()

            #=
            vec_x = zeros(Float64,2m+1)

            for i=0:2m
                x = i*Δx
                vec_x[i+1] = x
            end
            @time data = calc_vecfunc_cg(calc_f,vec_x,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)

            for i=1:m
                k = 2i-2
                k1 = 2i-1
                k2 = 2i
                x = k*Δx
                x1 = k1*Δx
                x2 = k2*Δx

                val = data[k+1]
                val1 = data[k1+1]
                val2 = data[k2+1]

                Δ += Δx*(val + 4*val1 + val2)/3
            end
            println(" Δ ",Δ)
            =#

            fp4 = open("Mnc.dat","w")
            an = chebyshevfit(calc_f,m,nc,M,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
            for M2 in [16,32,64,128,256]
                for nc2 in [10,20,50,100,150,200]
                    an = chebyshevfit(calc_f,m,nc2,M2,N,ANmat,Adiag,ϕ,Aϕ,ϕtemp1)
                    Δ3 = integcheb(an,0,1) + Δ0
                    println("M = $M2, nc = $nc2, Δ = $Δ3, exact: $ldet ")
                    println(fp4,"$M2 $nc2 $Δ3 $ldet ")
                end
            end
            close(fp4)
            exit()
            Δ3 = integcheb(an,0,1) + Δ0
            println("an: ",an," Δ ",Δ3)

            #aninv = chebyshevfit(calc_finv,m,nc)
            #println("aninv: ",aninv)
            #Δ3 = integcheb(aninv,0,1) + Δ0
            #println("aninv: ",aninv," Δ ",Δ3)
            exit()
        end

        if filename != nothing
            fp = open(filename,"w")
        end
        for i=1:m
            x = (2i-2)*Δx
            x1 = (2i-1)*Δx
            x2 = 2i*Δx

            if nonc

                val = calc_f(x)
                val1 = calc_f(x1)
                val2 = calc_f(x2)
            else
                valc = fitfunc(an,x)
                val1c = fitfunc(an,x1)
                val2c = fitfunc(an,x2)
            end

            if filename != nothing
                if nonc

                    println(fp,x,"\t",val)
                    println(fp,x1,"\t",val1)
                    println(fp,x2,"\t",val2)
                else

                    println(fp,x,"\t",valc)
                    println(fp,x1,"\t",val1c)
                    println(fp,x2,"\t",val2c)
                end

                flush(fp)
            end

            if nonc
                Δ += Δx*(val + 4*val1 + val2)/3
            else
                Δc += Δx*(valc + 4*val1c + val2c)/3
            end
            #println("i=$i ",Δ,"\t",ldet)

        end

        if nonc
            return real(Δ)
        else
            return real(Δ3)
        end

        println(Δ)
        println("cheb ",Δc)
        println("cheb2 ", Δ3)

        @time ldet = real(logdet(Matrix(A)))
        println("log det A(exact): $ldet")
        exit()

        if filename != nothing
            close(fp)
        end
        #exit()

        return real(Δ)

    end

    function tdlogdet(A::T,M,m;filename=nothing,nc = 20,nonc = false) where T <: DdagD_operator
        ic,ix,iy,iz,it,α = 1,1,1,1,1,1
        temps = A.dirac._temporal_fermi
        xi = temps[7]
        x0 = temps[8]
        Fermionfields.clear!(x0)
        Fermionfields.clear!(xi)
        x0[ic,ix,iy,iz,it,α] = 1
        mul!(xi,A,x0,(ix,iy,iz,it,α))
        diagvalue = xi[ic,ix,iy,iz,it,α]
        #println(diagvalue)

        N,_ = size(A)
        Δ0 = 0
        Δ0 = N*real(log(diagvalue))
        Δx = 1/2m

        ϕ = temps[7]
        ϕtemp1 = temps[8]
        Aϕ = temps[9]

        ANmat = DdagDND_Staggered_operator(A) #non-diagonal part of A

        #@time ldet = real(logdet(Matrix(A)))
        #println("log det A(exact): $ldet")

        count = 0

        function calc_f(t)
            DtN = DdagD_Staggered_operator(A,t)
            #DtN = Adiag .+ t*ANmat
            #count += 1
            #println(count)
            val = 0
            #println("random")
            for k=1:M
                gauss_distribution_fermi_Z2!(ϕ)

                cg(ϕtemp1,DtN,ϕ,eps = 1e-10,maxsteps= 3000)
                mul!(Aϕ,ANmat,ϕtemp1)
                val += ϕ ⋅ Aϕ /M
                println("t: $t ",val)
                exit()
                #val += ϕ'*Aϕ/M#ANmat*Dinvϕ/M
            end
            return real(val)
        end

        Δx = 1/2m
        Δ = Δ0
        Δc = Δ0

        #I,E = quadgk(calc_f, 0,1,maxevals=3m)
        #println(I)
        #println(E)
        #exit()

        #nc = 20

        if nonc == false
            an = chebyshevfit(calc_f,m,nc)
            Δ3 = integcheb(an,0,1) + Δ0
        end

        if filename != nothing
            fp = open(filename,"w")
        end
        for i=1:m
            x = (2i-2)*Δx
            x1 = (2i-1)*Δx
            x2 = 2i*Δx

            if nonc

                val = calc_f(x)
                val1 = calc_f(x1)
                val2 = calc_f(x2)
            else
                valc = fitfunc(an,x)
                val1c = fitfunc(an,x1)
                val2c = fitfunc(an,x2)
            end

            if filename != nothing
                if nonc

                    println(fp,x,"\t",val)
                    println(fp,x1,"\t",val1)
                    println(fp,x2,"\t",val2)
                else

                    println(fp,x,"\t",valc)
                    println(fp,x1,"\t",val1c)
                    println(fp,x2,"\t",val2c)
                end

                flush(fp)
            end

            if nonc
                Δ += Δx*(val + 4*val1 + val2)/3
            else
                Δc += Δx*(valc + 4*val1c + val2c)/3
            end
            #println("i=$i ",Δ,"\t",ldet)

        end

        if nonc
            return real(Δ)
        else
            return real(Δ3)
        end

    end
end


text:
"""
    lis(arr::Array{Int}, ::Val{:dp})

# Arguments:
- `arr`: sequence of integers

# Examples/Tests
```julia
julia> lis([3, 10, 2, 1, 20], Val(:dp))
[3, 10, 20]

julia> lis([2, 2, 3], Val(:dp))
[2, 3]
```

# Brief:
Given a sequence of integers, the function finds the longest, strictly increasing, subsequence in that sequence.

# For more information:
https://cp-algorithms.com/sequences/longest_increasing_subsequence.html

# Contributed by:
- [Igor Malheiros](https://github.com/igormalheiros)
"""

function lis(arr::Array{Int}, ::Val{:dp})
    len = length(arr)
    memo = ones(Int, len)
    p = zeros(Int, len)

    lis_arr = Int[]

    len == 0 && return lis_arr # if arr is empty

    lis_value = 1
    lis_pos = 1

    for i in 2:len
        for j in 1:i
            if arr[j] < arr[i] && memo[i] < memo[j] + 1
                memo[i] = memo[j] + 1
                p[i] = j
            end
        end
        # Updates the longest increasing value and position
        if memo[i] > lis_value
            lis_value = memo[i]
            lis_pos = i
        end
    end

    # Restoring
    while lis_pos != 0
        push!(lis_arr, arr[lis_pos])
        lis_pos = p[lis_pos]
    end

    reverse!(lis_arr)

    return lis_arr
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using Dates

framerate(10)

snake = Snake()
# The game's resolution
rez = 20
# The player's score
score = 0
# Last time the food was eaten
last_eaten = now()
# If the game is paused
paused = false

function createFood()
    # The food will be created in a random location each time
    global food = [rand(0:w-1), rand(0:h-1)]
    # The timer will be reseted
    global last_eaten = now()
end

@use function setup()
    width, height = @width, @height
    global w = width ÷ rez
    global h = height ÷ rez
    createFood()
end

@use function update()
    width, height = @width, @height
    # If the game is paused
    # draw the paused screen
    if paused
        draw(Screen{:Paused})
        return;
    end
    # Save the current transformation state
    gsave()
    # then scale by the resolution variable
    scale(rez)
    background("black")

    # if the player has eaten the food
    if eat(snake, food)
        createFood()
        global score += 1
    end

    # If 10 seconds has passed and
    # the player didn't have eaten the food
    # then a new one will be created
    now() - last_eaten > Second(10) && createFood()

    update(snake)
    show(snake)

    # Draw the food
    # the food is just one yellowgreen rectangle
    sethue("yellowgreen")
    rect(point(food), 1, 1, :fill)

    # If the game ended
    if gameover(snake)
        draw(Screen{:GameOver})
        return;
    end

    grestore()
    # Shows the player score on the top left
    fontsize(16)
    sethue("white")
    text("Score: $score", 15, 30)
end


text:
using Compat
import Compat.String

type BSONError
    _wrap_::Vector{UInt8}

    function BSONError()
        return new(Array(UInt8, 512))
    end
end
export BSONError

const domainDescs = @compat Dict{Int64,AbstractString}(1 => "JSON", 2 => "READER")
const errorDescs = @compat Dict{Int64,Any}(1 => (@compat Dict{Int64,AbstractString}(
        1 => "READ_CORRUPT_JS",
        2 => "READ_INVALID_PARAM",
        3 => "READ_CB_FAILURE"
    )), 2 => @compat Dict{Int64,AbstractString}(1 => "BADFD"))

convert(::Type{String}, bsonError::BSONError) = begin
    uint32s = reinterpret(UInt32, bsonError._wrap_)
    domain = uint32s[1]
    code = uint32s[2]
    return Compat.String(bsonError._wrap_[9:end])
end
export convert

string(bsonError::BSONError) = String(bsonError._wrap_)
export string

error(bsonError::BSONError) = error(string(bsonError))
export error