text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Simulated annealing search with geometrically decreasing temperature schedule.
Propose by randomly choosing an existing node from the parse tree,
and replacing the subtree by a randomly generated one.  Use PSA for parallel version over n_threads.
"""
module SA

export SAESParams, SAESResult, sa_search, exprsearch, SearchParams, SearchResult, get_derivtree
export accept_prob, estimate_temp_params
export psa_search, PSAESParams

using Reexport
using ExprSearch
using RLESUtils, GitUtils, SwapBuffers, CPUTimeUtils, Observers
using DerivationTrees
using GrammaticalEvolution
using Iterators

import DerivationTrees.initialize!
import ..ExprSearch: SearchParams, SearchResult, exprsearch, ExprProblem, initialize!,
    get_grammar, get_fitness
import Base: isless, copy!

type SAESParams <: SearchParams
  #tree params
  maxsteps::Int64

  #SA
  T1::Float64 #initial temperature
  alpha::Float64 #decay rate
  n_epochs::Int64 #iterations
  n_starts::Int64 #number of starts

  observer::Observer
end

type PSAESParams <: SearchParams
  n_threads::Int64
  sa_params::SAESParams
  observer::Observer
end

PSAESParams(n_threads::Int64, sa_params::SAESParams) = PSAESParams(n_threads, sa_params, Observer())

type SAState
  tree::DerivationTree
  fitness::Float64
  expr
end

SAState(tree::DerivationTree) = SAState(tree, realmax(Float64), 0)

type SAESResult <: SearchResult
  tree::DerivationTree
  fitness::Float64
  expr
  best_at_eval::Int64
  totalevals::Int64

  function SAESResult()
    result = new()
    result.fitness = realmax(Float64)
    result.expr = 0
    result.best_at_eval = 0
    result.totalevals = 0
    return result
  end

  function SAESResult(tree::DerivationTree)
    result = SAESResult()
    result.tree = tree
    return result
  end
end

exprsearch(p::SAESParams, problem::ExprProblem, userargs...) = sa_search(p, problem, userargs...)
exprsearch(p::PSAESParams, problem::ExprProblem, userargs...) = psa_search(p, problem, userargs...)

function psa_search(p::PSAESParams, problem::ExprProblem, userargs...)
  @notify_observer(p.observer, "computeinfo", ["starttime", string(now())])
  tic()

  results = pmap(1:p.n_threads) do tid
    sa_search(p.sa_params, problem, userargs...)
  end

  result = minimum(results) #best fitness
  totalevals = sum(map(r -> r.totalevals, results))

  @notify_observer(p.observer, "result", [result.fitness, string(result.expr), totalevals])

  #meta info
  computetime_s = toq()
  @notify_observer(p.observer, "computeinfo", ["computetime_s",  computetime_s])
  @notify_observer(p.observer, "computeinfo", ["endtime",  string(now())])
  @notify_observer(p.observer, "computeinfo", ["hostname", gethostname()])
  @notify_observer(p.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
  @notify_observer(p.observer, "parameters", ["maxsteps", p.sa_params.maxsteps])
  @notify_observer(p.observer, "parameters", ["T1", p.sa_params.T1])
  @notify_observer(p.observer, "parameters", ["alpha", p.sa_params.alpha])
  @notify_observer(p.observer, "parameters", ["n_epochs", p.sa_params.n_epochs])
  @notify_observer(p.observer, "parameters", ["n_starts", p.sa_params.n_starts])
  @notify_observer(p.observer, "parameters", ["n_threads", p.n_threads])

  return result
end

function sa_search(p::SAESParams, problem::ExprProblem, userargs...)
  @notify_observer(p.observer, "verbose1", ["Starting SA search"])
  @notify_observer(p.observer, "computeinfo", ["starttime", string(now())])

  initialize!(problem)
  grammar = get_grammar(problem)
  tree_params = DerivTreeParams(grammar, p.maxsteps)
  s_buffer = SwapBuffer(SAState(DerivationTree(tree_params)),
                        SAState(DerivationTree(tree_params)))

  result = SAESResult()

  tstart = CPUtime_start()
  for j = 1:p.n_starts
    s = active(s_buffer)
    initialize!(s, problem, result) #initialize randomly

    T = p.T1
    for i = 1:p.n_epochs
      @notify_observer(p.observer, "iteration", [j, i])
      @notify_observer(p.observer, "temperature", [j, i, T])

      ###########################################
      # SA algorithm
      s = active(s_buffer)
      sp = inactive(s_buffer)
      neighbor!(sp, s, problem, result)

      #only accept complete trees
      if iscomplete(sp.tree) && accept(s, sp, T)
        swap!(s_buffer)
      end

      #update temperature
      T *= p.alpha
      ###########################################

      @notify_observer(p.observer, "elapsed_cpu_s", [j, i, CPUtime_elapsed_s(tstart)])
      @notify_observer(p.observer, "current_best", [j, i, result.fitness, string(result.expr)])
    end
    gc()
  end

  @notify_observer(p.observer, "result", [result.fitness, string(result.expr), result.best_at_eval, result.totalevals])

  #meta info
  @notify_observer(p.observer, "computeinfo", ["endtime",  string(now())])
  @notify_observer(p.observer, "computeinfo", ["hostname", gethostname()])
  @notify_observer(p.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
  @notify_observer(p.observer, "parameters", ["maxsteps", p.maxsteps])
  @notify_observer(p.observer, "parameters", ["T1", p.T1])
  @notify_observer(p.observer, "parameters", ["alpha", p.alpha])
  @notify_observer(p.observer, "parameters", ["n_epochs", p.n_epochs])
  @notify_observer(p.observer, "parameters", ["n_starts", p.n_starts])

  return result
end

#initialize to random state
function initialize!(s::SAState, problem::ExprProblem, retries::Int64=typemax(Int64))
  rand!(s.tree; retries=retries) #random tree
  s.expr = get_expr(s.tree)
  s.fitness = get_fitness(problem, s.expr)
end

function initialize!(s::SAState, problem::ExprProblem, result::SAESResult,
                     retries::Int64=typemax(Int64))
  initialize!(s, problem, retries)
  update!(result, s)
end

#get a new neighbor and update global trackers
function neighbor!(sp::SAState, s::SAState, problem::ExprProblem, retries::Int64=20)
  while retries > 0
    copy!(sp, s)
    #perturb and evaluate
    if perturb!(sp.tree) #true if successful
      sp.expr = get_expr(sp.tree)
      sp.fitness = get_fitness(problem, sp.expr)
      return sp #keep looping until successful
    end
    retries -= 1
  end
  error("neighbor!(): Max retries exceeded")
end

#update the global best trackers with the current state
function update!(result::SAESResult, s::SAState)
  result.totalevals += 1 #assumes an eval was called prior

  #update globals
  if s.fitness < result.fitness
    result.tree = s.tree
    result.fitness = s.fitness
    result.expr = s.expr
    result.best_at_eval = result.totalevals
  end
end

#get neighbor and update
function neighbor!(sp::SAState, s::SAState, problem::ExprProblem, result::SAESResult)
  neighbor!(sp, s, problem)
  update!(result, sp)
end

#randomly perturb the tree
function perturb!(tree::DerivationTree, retries::Int64=5)
  node = rand(tree) #uniformly select a node
  b_success = rand!(tree, node, retries) #randomly change the subtree
  return b_success
end

#accept rule for simulated annealing, always accept a downhill step or sidestep
function accept(s::SAState, sp::SAState, T::Float64)
  return sp.fitness <= s.fitness || exp((s.fitness - sp.fitness) / T) > rand()
end

############
#for estimating parameters before actual run
function accept_prob(problem::ExprProblem, maxsteps::Int64, P1::Float64, N::Int64)
  sum_uphill = 0.0
  n_uphill = 0

  grammar = get_grammar(problem)
  tree_params = DerivTreeParams(grammar, maxsteps)
  tree = DerivationTree(tree_params)

  s = SAState(tree)
  initialize!(s, problem)

  for i = 1:N
    #println("i=", i)
    sp = neighbor(s, problem)

    #println("s=", s.expr, ", fitness=", s.fitness)
    #println("sp=", sp.expr, ", fitness=", sp.fitness)

    iscomplete(sp.tree) || continue

    if sp.fitness <= s.fitness || rand() < P1
      #track uphill moves
      if sp.fitness > s.fitness
        sum_uphill += sp.fitness - s.fitness
        n_uphill += 1
      end
      s = sp
    end
  end

  avg_uphill = n_uphill != 0.0 ? sum_uphill / n_uphill : 0.0
  P = n_uphill / N

  #P should be close to P1

  return (avg_uphill, P)
end

#Reference: Park and Kim, "A Systematic Procedure for Setting Parameters in Simulated Annealing"
function estimate_temp_params(problem::ExprProblem,
                           P1::Float64,
                           n_epochs::Int64,
                           Tfinal::Float64,
                           maxsteps::Int64,
                           N::Int64,
                           ntrials::Int64=1000)

  avg_uphills = Float64[]
  Ps = Float64[]
  for i = 1:ntrials
    avg_uphill, P = accept_prob(problem, maxsteps, P1, N)
    push!(avg_uphills, avg_uphill)
    push!(Ps, P)
    println("avg_uphill=", avg_uphill, ", P=", P)
  end

  deltabar = mean(avg_uphills)
  #P1 = mean(Ps)
  T1 = -deltabar / log(P1)
  println("deltabar=", deltabar, ", P1=", P1, ", T1=", T1)

  alpha = (Tfinal / T1)^(1 / n_epochs)
  println("T1=", T1, ", alpha=", alpha, ", epochs=", n_epochs)

  return T1, alpha, n_epochs
end

isless(r1::SAESResult, r2::SAESResult) = r1.fitness < r2.fitness

function copy!(dst::SAState, src::SAState)
  copy!(dst.tree, src.tree)
  dst.fitness = src.fitness
  dst.expr = src.expr
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
# update_eta_Met.jl

## To do:
# - construct_Dh doesn't need h as an argument

export βδ_x_h_modify_γ, δ_x_h_modify_γ;

function update_η_h_Met!(model::Model_BNP_WMReg_Joint, h::Int, Λβ0star_ηy::Array{T,1}, βΛβ0star_ηy::T; update_ηy::Bool=true) where T <: Real
    if model.Σx_type == :full
        update_η_h_Met_Σx_full!(model::Model_BNP_WMReg_Joint, h::Int, Λβ0star_ηy::Array{T,1}, βΛβ0star_ηy::T, update_ηy=update_ηy)
    elseif model.Σx_type == :diag
        update_η_h_Met_Σx_diag!(model::Model_BNP_WMReg_Joint, h::Int, Λβ0star_ηy::Array{T,1}, βΛβ0star_ηy::T, update_ηy=update_ηy)
    end
end

# pre_compute Λ0star_ηy*β0star_ηy and β0star_ηy'Λ0star_ηy*β0star_ηy which are not indexed by h
function update_η_h_Met_Σx_full!(model::Model_BNP_WMReg_Joint, h::Int, Λβ0star_ηy::Array{T,1}, βΛβ0star_ηy::T; update_ηy::Bool=true) where T <: Real

    ## if you modify this function, you must modify its sisters.

    indx_h = findall(model.state.S .== h)
    n_h = length(indx_h)

    μ_x_h_old = model.state.μ_x[h,:]
    lδ_x_h_old = log.(model.state.δ_x[h,:])
    β_x_h_old = [ model.state.β_x[k][h,:] for k = 1:(model.K - 1) ]
    ηlδ_x_old = vcat(μ_x_h_old,
        vcat(β_x_h_old...),
        lδ_x_h_old)

    ## Draw candidate for η_x
    ηlδ_x_cand = ηlδ_x_old + rand(MvNormal(model.state.cSig_ηlδx[h]))
    μ_x_h_cand = ηlδ_x_cand[model.indx_ηx[:μ]]
    β_x_h_candvec = ηlδ_x_cand[model.indx_ηx[:β]]
    β_x_h_cand = [ β_x_h_candvec[model.indx_β_x[k]] for k = 1:(model.K - 1) ]
    lδ_x_h_cand = ηlδ_x_cand[model.indx_ηx[:δ]]
    δ_x_h_cand = exp.(lδ_x_h_cand)

    auto_reject = (!all(δ_x_h_cand .> 0.0)) || any(isinf.(δ_x_h_cand))
    not_auto_reject = !auto_reject

    lNX_mat_cand = deepcopy(model.state.lNX)

    ## bookkeeping for variable selection
    if typeof(model.state.γ) == BitArray{1}
        γ_h = deepcopy(model.state.γ)
    elseif typeof(model.state.γ) == BitArray{2}
        γ_h = deepcopy(model.state.γ[h,:])
    end
    γindx = findall(γ_h)
    nγ = length(γindx)

    if not_auto_reject
        if model.state.γδc == Inf # subset method
            if nγ == 0
                lNX_mat_cand[:,h] = zeros(Float64, model.n)
            elseif nγ == 1
                lNX_mat_cand[:,h] = logpdf.(Normal(μ_x_h_cand[γindx[1]], sqrt(δ_x_h_cand[γindx[1]])), vec(model.X[:,γindx]))
            elseif nγ > 1
                βγ_x_h_cand, δγ_x_h_cand = βδ_x_h_modify_γ(β_x_h_cand, δ_x_h_cand,
                    γ_h, model.state.γδc) # either variance-inflated or subset

                tmp = lNX_sqfChol(permutedims(model.X[:,γindx]),
                    μ_x_h_cand[γindx], βγ_x_h_cand, δγ_x_h_cand, false) # false means set to return nothing if not pos.def.

                if isnothing(tmp)
                    auto_reject = true
                    not_auto_reject = !auto_reject
                else
                    lNX_mat_cand[:,h] = deepcopy(tmp)
                end
            end

            lωNX_vec_cand = lωNXvec(model.state.lω, lNX_mat_cand) # n vector

        elseif isnothing(model.state.γδc) # integration method
            if nγ == 0
                lNX_mat_cand[:,h] = zeros(Float64, model.n)
            elseif nγ == 1
                tmp = sqfChol_to_Σ(β_x_h_cand, δ_x_h_cand, false) # false means don't throw error if not PosDef

                if isnothing(tmp)
                    auto_reject = true
                    not_auto_reject = !auto_reject
                else
                    σ2x = tmp.mat[γindx, γindx][1]
                    lNX_mat_cand[:,h] = logpdf.(Normal(μ_x_h_cand[γindx[1]], sqrt(σ2x)), vec(model.X[:,γindx]))
                end
            elseif nγ > 1
                tmp = sqfChol_to_Σ(β_x_h_cand, δ_x_h_cand, false) # false means don't throw error if not PosDef

                if isnothing(tmp)
                    auto_reject = true
                    not_auto_reject = !auto_reject
                else
                    Σx = PDMat_adj(tmp.mat[γindx, γindx])
                    lNX_mat_cand[:,h] = logpdf(MultivariateNormal(μ_x_h_cand[γindx], Σx), permutedims(model.X[:,γindx]))
                end
            end

            lωNX_vec_cand = lωNXvec(model.state.lω, lNX_mat_cand) # n vector

        else # variance-inflation method

            throw("Variance-inflation method no longer supported.")

            # βγ_x_h_cand, δγ_x_h_cand = βδ_x_h_modify_γ(β_x_h_cand, δ_x_h_cand,
            #     γ_h, model.state.γδc) # either variance-inflated or subset

            # tmp = lNX_sqfChol(permutedims(model.X), μ_x_h_cand, βγ_x_h_cand, δγ_x_h_cand, false) # false means don't throw error if not PosDef

            # if isnothing(tmp)
            #     auto_reject = true
            #     not_auto_reject = !auto_reject
            # else
            #     lNX_mat_cand[:,h] = deepcopy(tmp)
            #     lωNX_vec_cand = lωNXvec(model.state.lω, lNX_mat_cand) # n vector
            # end

        end
    end

    if n_h == 0

        ## Metropolis step for η_x

        ## Compute acceptance ratio
        if not_auto_reject
            lar = lG0_ηlδx(μ_x_h_cand, β_x_h_cand, lδ_x_h_cand, model.state) -
                    sum(lωNX_vec_cand) -
                    lG0_ηlδx(μ_x_h_old, β_x_h_old, lδ_x_h_old, model.state) +
                    sum(model.state.lωNX_vec)
        end

        ## Decision and update
        if not_auto_reject && log(rand()) < lar # accept

            model.state.μ_x[h,:] = deepcopy(μ_x_h_cand)
            for k = 1:(model.K - 1)
                model.state.β_x[k][h,:] = deepcopy(β_x_h_cand[k])
            end
            model.state.δ_x[h,:] = deepcopy(δ_x_h_cand)

            model.state.lNX[:,h] = deepcopy(lNX_mat_cand[:,h])
            model.state.lωNX_vec = deepcopy(lωNX_vec_cand)

            ηlδ_x_out = deepcopy(ηlδ_x_cand)
            model.state.accpt[h] += 1

        else # reject
            ηlδ_x_out = deepcopy(ηlδ_x_old)
        end

        ## Full conditional draw (from G0) for η_y

        if update_ηy
            model.state.δ_y[h] = rand(InverseGamma(0.5 * model.state.ν_δy,
                                        0.5 * model.state.ν_δy * model.state.s0_δy))
            βstar_ηy = model.state.β0star_ηy + (model.state.Λ0star_ηy.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
            model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
            model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
        end

    else # some observations are assigned to cluster h (i.e., n_h > 0)

        ## Metropolis step for η_x

            ## Candidate, lNX, lωNXvec, previously computed
            ## Important quantities

        y_h = model.y[indx_h] # doesn't need deepcopy if indexing (result of call to getindex)
        X_h = model.X[indx_h,:]
        D_h_old = construct_Dh(h, X_h, model.state.μ_x[h,:], γ_h)

        Λ1star_ηy_h_old = get_Λ1star_ηy_h(D_h_old, model.state.Λ0star_ηy)
        β1star_ηy_h_old = get_β1star_ηy_h(Λ1star_ηy_h_old, Λβ0star_ηy, D_h_old, y_h)

        a1_δy_old = get_a1_δy_h(model.state.ν_δy, n_h) # posterior IG shape
        b1_δy_old = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, y_h, βΛβ0star_ηy, Λ1star_ηy_h_old, β1star_ηy_h_old) # posterior IG scale

        ## Important quantities for candidate
        if not_auto_reject
            D_h_cand = construct_Dh(h, X_h, μ_x_h_cand, γ_h)

            Λ1star_ηy_h_cand = get_Λ1star_ηy_h(D_h_cand, model.state.Λ0star_ηy)
            β1star_ηy_h_cand = get_β1star_ηy_h(Λ1star_ηy_h_cand, Λβ0star_ηy, D_h_cand, y_h)

            a1_δy_cand = get_a1_δy_h(model.state.ν_δy, n_h) # posterior IG shape
            b1_δy_cand = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, y_h, βΛβ0star_ηy, Λ1star_ηy_h_cand, β1star_ηy_h_cand) # posterior IG scale

            ## Compute acceptance ratio (lcc_ηlδx does not need γ-modified βx and δx)
            lar = lcc_ηlδx(h, indx_h, lNX_mat_cand, lωNX_vec_cand,
                            model.state, μ_x_h_cand, β_x_h_cand, lδ_x_h_cand,
                            Λ1star_ηy_h_cand, a1_δy_cand, b1_δy_cand) -
                lcc_ηlδx(h, indx_h, model.state.lNX, model.state.lωNX_vec,
                            model.state, μ_x_h_old, β_x_h_old, lδ_x_h_old,
                            Λ1star_ηy_h_old, a1_δy_old, b1_δy_old)

        end

        ## Decision and update
        if not_auto_reject && log(rand()) < lar # accept

            model.state.μ_x[h,:] = deepcopy(μ_x_h_cand)
            for k = 1:(model.K - 1)
                model.state.β_x[k][h,:] = deepcopy(β_x_h_cand[k])
            end
            model.state.δ_x[h,:] = deepcopy(δ_x_h_cand)

            model.state.lNX[:,h] = deepcopy(lNX_mat_cand[:,h])
            model.state.lωNX_vec = deepcopy(lωNX_vec_cand)

            ηlδ_x_out = deepcopy(ηlδ_x_cand)
            model.state.accpt[h] += 1

            ## Full conditional draw for η_y ( use prec. prameterization )
            if update_ηy
                model.state.δ_y[h] = rand(InverseGamma(a1_δy_cand, b1_δy_cand))
                βstar_ηy = β1star_ηy_h_cand + (Λ1star_ηy_h_cand.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
                model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
                model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
            end

        else # reject

            ηlδ_x_out = deepcopy(ηlδ_x_old)

            ## Full conditional draw for η_y
            if update_ηy
                model.state.δ_y[h] = rand(InverseGamma(a1_δy_old, b1_δy_old))
                βstar_ηy = β1star_ηy_h_old + (Λ1star_ηy_h_old.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
                model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
                model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
            end
        end
    end

    if model.state.adapt && (model.state.iter % model.state.adapt_thin == 0)
        model.state.adapt_iter += 1
        model.state.runningsum_ηlδx[h,:] += ηlδ_x_out
        runningmean = model.state.runningsum_ηlδx[h,:] / float(model.state.adapt_iter)
        runningdev = ( ηlδ_x_out - runningmean )
        model.state.runningSS_ηlδx[h,:,:] = runningdev * runningdev' # the mean is changing, but this approx. is fine.
    end

    return nothing
end

function update_η_h_Met_Σx_diag!(model::Model_BNP_WMReg_Joint, h::Int, Λβ0star_ηy::Array{T,1}, βΛβ0star_ηy::T; update_ηy::Bool=true) where T <: Real

    ## if you modify this function, you must modify its sisters.

    indx_h = findall(model.state.S .== h)
    n_h = length(indx_h)

    μ_x_h_old = model.state.μ_x[h,:]
    lδ_x_h_old = log.(model.state.δ_x[h,:])
    ηlδ_x_old = vcat(μ_x_h_old,
        # vcat(β_x_h_old...),
        lδ_x_h_old)

    ## Draw candidate for η_x
    ηlδ_x_cand = ηlδ_x_old + rand(MvNormal(model.state.cSig_ηlδx[h]))
    μ_x_h_cand = ηlδ_x_cand[model.indx_ηx[:μ]]
    lδ_x_h_cand = ηlδ_x_cand[model.indx_ηx[:δ]]
    δ_x_h_cand = exp.(lδ_x_h_cand)

    auto_reject = (!all(δ_x_h_cand .> 0.0)) || any(isinf.(δ_x_h_cand))
    not_auto_reject = !auto_reject

    lNX_mat_cand = deepcopy(model.state.lNX)

    ## bookkeeping for variable selection
    if typeof(model.state.γ) == BitArray{1}
        γ_h = deepcopy(model.state.γ)
    elseif typeof(model.state.γ) == BitArray{2}
        γ_h = deepcopy(model.state.γ[h,:])
    end
    γindx = findall(γ_h)
    nγ = length(γindx)

    if not_auto_reject
        if model.state.γδc == Inf || isnothing(model.state.γδc) # subset and integration methods are identical with SigX_diag
            if nγ == 0
                lNX_mat_cand[:,h] = zeros(Float64, model.n)
            elseif nγ == 1
                lNX_mat_cand[:,h] = logpdf.(Normal(μ_x_h_cand[γindx[1]], sqrt(δ_x_h_cand[γindx[1]])), vec(model.X[:,γindx]))
            elseif nγ > 1
                lNX_mat_cand[:,h] = lNX_Σdiag(permutedims(model.X[:,γindx]),
                    μ_x_h_cand[γindx], δ_x_h_cand[γindx])
            end
        else # variance-inflation method
            throw("Variance inflation method not supported when Sigma_x is diagonal.")
        end

        lωNX_vec_cand = lωNXvec(model.state.lω, lNX_mat_cand) # n vector
    end

    if n_h == 0

        ## Metropolis step for η_x

        ## Compute acceptance ratio
        if not_auto_reject
            lar = lG0_ηlδx(μ_x_h_cand, nothing, lδ_x_h_cand, model.state) -
                    sum(lωNX_vec_cand) -
                    lG0_ηlδx(μ_x_h_old, nothing, lδ_x_h_old, model.state) +
                    sum(model.state.lωNX_vec)
        end

        ## Decision and update
        if not_auto_reject && log(rand()) < lar # accept

            model.state.μ_x[h,:] = deepcopy(μ_x_h_cand)
            model.state.δ_x[h,:] = deepcopy(δ_x_h_cand)

            model.state.lNX[:,h] = deepcopy(lNX_mat_cand[:,h])
            model.state.lωNX_vec = deepcopy(lωNX_vec_cand)

            ηlδ_x_out = deepcopy(ηlδ_x_cand)
            model.state.accpt[h] += 1

        else # reject
            ηlδ_x_out = deepcopy(ηlδ_x_old)
        end

        ## Full conditional draw (from G0) for η_y
        if update_ηy
            model.state.δ_y[h] = rand(InverseGamma(0.5 * model.state.ν_δy,
                                        0.5 * model.state.ν_δy * model.state.s0_δy))
            βstar_ηy = model.state.β0star_ηy + (model.state.Λ0star_ηy.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
            model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
            model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
        end
    else # some observations are assigned to cluster h (i.e., n_h > 0)

        ## Metropolis step for η_x

        ## Candidate, lNX, lωNXvec, previously computed
        ## Important quantities

        y_h = model.y[indx_h] # doesn't need deepcopy if indexing (result of call to getindex)
        X_h = model.X[indx_h,:]
        D_h_old = construct_Dh(h, X_h, model.state.μ_x[h,:], γ_h)

        Λ1star_ηy_h_old = get_Λ1star_ηy_h(D_h_old, model.state.Λ0star_ηy)
        β1star_ηy_h_old = get_β1star_ηy_h(Λ1star_ηy_h_old, Λβ0star_ηy, D_h_old, y_h)

        a1_δy_old = get_a1_δy_h(model.state.ν_δy, n_h) # posterior IG shape
        b1_δy_old = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, y_h, βΛβ0star_ηy, Λ1star_ηy_h_old, β1star_ηy_h_old) # posterior IG scale

        ## Important quantities for candidate
        if not_auto_reject
            D_h_cand = construct_Dh(h, X_h, μ_x_h_cand, γ_h)

            Λ1star_ηy_h_cand = get_Λ1star_ηy_h(D_h_cand, model.state.Λ0star_ηy)
            β1star_ηy_h_cand = get_β1star_ηy_h(Λ1star_ηy_h_cand, Λβ0star_ηy, D_h_cand, y_h)

            a1_δy_cand = get_a1_δy_h(model.state.ν_δy, n_h) # posterior IG shape
            b1_δy_cand = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, y_h, βΛβ0star_ηy, Λ1star_ηy_h_cand, β1star_ηy_h_cand) # posterior IG scale

            ## Compute acceptance ratio (lcc_ηlδx does not need γ-modified βx and δx)
            lar = lcc_ηlδx(h, indx_h, lNX_mat_cand, lωNX_vec_cand,
                            model.state, μ_x_h_cand, nothing, lδ_x_h_cand,
                            Λ1star_ηy_h_cand, a1_δy_cand, b1_δy_cand) -
                lcc_ηlδx(h, indx_h, model.state.lNX, model.state.lωNX_vec,
                            model.state, μ_x_h_old, nothing, lδ_x_h_old,
                            Λ1star_ηy_h_old, a1_δy_old, b1_δy_old)

        end

        ## Decision and update
        if not_auto_reject && log(rand()) < lar # accept

            model.state.μ_x[h,:] = deepcopy(μ_x_h_cand)
            model.state.δ_x[h,:] = deepcopy(δ_x_h_cand)

            model.state.lNX[:,h] = deepcopy(lNX_mat_cand[:,h])
            model.state.lωNX_vec = deepcopy(lωNX_vec_cand)

            ηlδ_x_out = deepcopy(ηlδ_x_cand)
            model.state.accpt[h] += 1

            ## Full conditional draw for η_y ( use prec. prameterization )
            if update_ηy
                model.state.δ_y[h] = rand(InverseGamma(a1_δy_cand, b1_δy_cand))
                βstar_ηy = β1star_ηy_h_cand + (Λ1star_ηy_h_cand.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
                model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
                model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
            end
        else # reject

            ηlδ_x_out = deepcopy(ηlδ_x_old)

            ## Full conditional draw for η_y
            if update_ηy
                model.state.δ_y[h] = rand(InverseGamma(a1_δy_old, b1_δy_old))
                βstar_ηy = β1star_ηy_h_old + (Λ1star_ηy_h_old.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
                model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
                model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
            end
        end
    end

    if model.state.adapt && (model.state.iter % model.state.adapt_thin == 0)
        model.state.adapt_iter += 1
        model.state.runningsum_ηlδx[h,:] += ηlδ_x_out
        runningmean = model.state.runningsum_ηlδx[h,:] / float(model.state.adapt_iter)
        runningdev = ( ηlδ_x_out - runningmean )
        model.state.runningSS_ηlδx[h,:,:] = runningdev * runningdev' # the mean is changing, but this approx. is fine.
    end

    return nothing
end

function update_η_h_Met_K1!(model::Model_BNP_WMReg_Joint, h::Int, Λβ0star_ηy::Array{T,1}, βΛβ0star_ηy::T; update_ηy::Bool=true) where T <: Real

    ## if you modify this function, you must modify the originals above.

    indx_h = findall(model.state.S .== h)
    n_h = length(indx_h)

    μ_x_h_old = model.state.μ_x[h,:]
    lδ_x_h_old = log.(model.state.δ_x[h,:])
    ηlδ_x_old = vcat(μ_x_h_old, lδ_x_h_old)

    ## Draw candidate for η_x
    ηlδ_x_cand = ηlδ_x_old + rand(MvNormal(model.state.cSig_ηlδx[h]))
    μ_x_h_cand = ηlδ_x_cand[model.indx_ηx[:μ]]
    lδ_x_h_cand = ηlδ_x_cand[model.indx_ηx[:δ]]
    δ_x_h_cand = exp.(lδ_x_h_cand)

    auto_reject = (!all(δ_x_h_cand .> 0.0)) || any(isinf.(δ_x_h_cand))
    not_auto_reject = !auto_reject

    lNX_mat_cand = deepcopy(model.state.lNX)

    if typeof(model.state.γ) == BitArray{1}
        γ_h = deepcopy(model.state.γ)
    elseif typeof(model.state.γ) == BitArray{2}
        γ_h = deepcopy(model.state.γ[h,:])
    end

    ## bookkeeping for variable selection
    if not_auto_reject
        if model.state.γδc == Inf || isnothing(model.state.γδc) # subset or integration method
            if γ_h[1]
                lNX_mat_cand[:,h] = logpdf.(Normal(μ_x_h_cand[1], sqrt(δ_x_h_cand[1])), vec(model.X))
            else
                lNX_mat_cand[:,h] = zeros(Float64, model.n)
            end
        else # variance-inflation method
            δγ_x_h_cand = δ_x_h_modify_γ(δ_x_h_cand, γ_h, model.state.γδc)
            lNX_mat_cand[:,h] = logpdf.(Normal(μ_x_h_cand[1], sqrt(δγ_x_h_cand[1])), vec(model.X))
        end

        lωNX_vec_cand = lωNXvec(model.state.lω, lNX_mat_cand) # n vector

    end

    if n_h == 0

        ## Metropolis step for η_x

        ## Compute acceptance ratio
        if not_auto_reject
            lar = lG0_ηlδx(μ_x_h_cand[1], lδ_x_h_cand[1], model.state) -
                    sum(lωNX_vec_cand) -
                    lG0_ηlδx(μ_x_h_old[1], lδ_x_h_old[1], model.state) +
                    sum(model.state.lωNX_vec)
        end

        ## Decision and update
        if not_auto_reject && log(rand()) < lar # accept

            model.state.μ_x[h,:] = deepcopy(μ_x_h_cand)
            model.state.δ_x[h,:] = deepcopy(δ_x_h_cand)

            model.state.lNX[:,h] = deepcopy(lNX_mat_cand[:,h])
            model.state.lωNX_vec = deepcopy(lωNX_vec_cand)

            ηlδ_x_out = deepcopy(ηlδ_x_cand)
            model.state.accpt[h] += 1

        else # reject
            ηlδ_x_out = deepcopy(ηlδ_x_old)
        end

        ## Full conditional draw (from G0) for η_y
        if update_ηy
            model.state.δ_y[h] = rand(InverseGamma(0.5 * model.state.ν_δy,
                                        0.5 * model.state.ν_δy * model.state.s0_δy))
            βstar_ηy = model.state.β0star_ηy + (model.state.Λ0star_ηy.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
            model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
            model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
        end
    else # some observations are assigned to cluster h (i.e., n_h > 0)

        ## Metropolis step for η_x

        ## Candidate, lNX, lωNXvec, previously computed
        ## Important quantities

        y_h = model.y[indx_h] # doesn't need to copy if indexing (result of call to getindex)
        X_h = model.X[indx_h,:]
        D_h_old = construct_Dh(h, X_h, model.state.μ_x[h,:], γ_h)

        Λ1star_ηy_h_old = get_Λ1star_ηy_h(D_h_old, model.state.Λ0star_ηy)
        β1star_ηy_h_old = get_β1star_ηy_h(Λ1star_ηy_h_old, Λβ0star_ηy, D_h_old, y_h)

        a1_δy_old = get_a1_δy_h(model.state.ν_δy, n_h) # posterior IG shape
        b1_δy_old = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, y_h, βΛβ0star_ηy, Λ1star_ηy_h_old, β1star_ηy_h_old) # posterior IG scale

        if not_auto_reject
            ## Important quantities for candidate
            D_h_cand = construct_Dh(h, X_h, μ_x_h_cand, γ_h)

            Λ1star_ηy_h_cand = get_Λ1star_ηy_h(D_h_cand, model.state.Λ0star_ηy)
            β1star_ηy_h_cand = get_β1star_ηy_h(Λ1star_ηy_h_cand, Λβ0star_ηy, D_h_cand, y_h)

            a1_δy_cand = get_a1_δy_h(model.state.ν_δy, n_h) # posterior IG shape
            b1_δy_cand = get_b1_δy_h(model.state.ν_δy, model.state.s0_δy, y_h, βΛβ0star_ηy, Λ1star_ηy_h_cand, β1star_ηy_h_cand) # posterior IG scale

            ## Compute acceptance ratio (lcc_ηlδx does not need γ-modified δx)
            lar = lcc_ηlδx(h, indx_h, lNX_mat_cand, lωNX_vec_cand,
                            model.state, μ_x_h_cand[1], lδ_x_h_cand[1],
                            Λ1star_ηy_h_cand, a1_δy_cand, b1_δy_cand) -
                        lcc_ηlδx(h, indx_h, model.state.lNX, model.state.lωNX_vec,
                            model.state, μ_x_h_old[1], lδ_x_h_old[1],
                            Λ1star_ηy_h_old, a1_δy_old, b1_δy_old)
        end

        ## Decision and update
        if not_auto_reject && log(rand()) < lar  # accept

            model.state.μ_x[h,:] = deepcopy(μ_x_h_cand)
            model.state.δ_x[h,:] = deepcopy(δ_x_h_cand)

            model.state.lNX[:,h] = deepcopy(lNX_mat_cand[:,h])
            model.state.lωNX_vec = deepcopy(lωNX_vec_cand)

            ηlδ_x_out = deepcopy(ηlδ_x_cand)
            model.state.accpt[h] += 1

            ## Full conditional draw for η_y ( use prec. prameterization )
            if update_ηy
                model.state.δ_y[h] = rand(InverseGamma(a1_δy_cand, b1_δy_cand))
                βstar_ηy = β1star_ηy_h_cand + (Λ1star_ηy_h_cand.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
                model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
                model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
            end
        else # reject

            ηlδ_x_out = deepcopy(ηlδ_x_old)

            ## Full conditional draw for η_y
            if update_ηy
                model.state.δ_y[h] = rand(InverseGamma(a1_δy_old, b1_δy_old))
                βstar_ηy = β1star_ηy_h_old + (Λ1star_ηy_h_old.chol.U \ randn((model.K + 1)) .* sqrt(model.state.δ_y[h]))
                model.state.μ_y[h] = deepcopy(βstar_ηy[model.indx_ηy[:μ]])
                model.state.β_y[h,:] = deepcopy(βstar_ηy[model.indx_ηy[:β]])
            end
        end
    end

    if model.state.adapt && (model.state.iter % model.state.adapt_thin == 0)
        model.state.adapt_iter += 1
        model.state.runningsum_ηlδx[h,:] += ηlδ_x_out
        runningmean = model.state.runningsum_ηlδx[h,:] / float(model.state.adapt_iter)
        runningdev = ( ηlδ_x_out - runningmean )
        model.state.runningSS_ηlδx[h,:,:] = runningdev * runningdev' # the mean is changing, but this approx. is fine.
    end

    return nothing
end

function construct_Dh(h::Int, Xh::Array{T,2}, μ_x_h::Array{T,1}) where T <: Real
    D0 = broadcast(-, permutedims(μ_x_h), Xh)
    return hcat(ones(size(Xh, 1)), D0)
end
function construct_Dh(h::Int, Xh::Array{T,2}, μ_x_h::Array{T,1},
                      γ::BitArray{1}) where T <: Real
    D0 = broadcast(-, permutedims(μ_x_h), Xh)
    D0[:, findall(.!(γ))] *= 0.0
    return hcat(ones(size(Xh, 1)), D0)
end

function lG0_ηlδx(μ_x::Array{T,1}, β_x::Union{Array{Array{T,1},1}, Nothing},
                  lδ_x::Array{T,1},
                  state::State_BNP_WMReg_Joint) where T <: Real

    K = length(μ_x)
    Q_μ = - 0.5 * PDMats.quad(state.Λ0_μx, (μ_x - state.μ0_μx))
    Q_β = 0.0
    Q_δ = 0.0
    if !isnothing(β_x)
        for k = 1:(K - 1)
            Q_β -= 0.5 * PDMats.quad(state.Λ0_βx[k], (β_x[k] - state.β0_βx[k]))
            # Q_δ -= ( (0.5*state.ν_δx[k] + 1.0)*lδ_x[k] +
            #           0.5*state.ν_δx[k]*state.s0_δx[k]/exp(lδ_x[k]) ) # build Jacobian into it
            # Q_δ -= 0.5*( state.ν_δx[k]*lδ_x[k] + state.ν_δx[k]*state.s0_δx[k]/exp(lδ_x[k]) ) # Jacobian built in
        end
    end
    # Q_δ -= ( (0.5*state.ν_δx[K] + 1.0)*lδ_x[K] +
    #           0.5*state.ν_δx[K]*state.s0_δx[K]/exp(lδ_x[K]) )
    for k = 1:K
        Q_δ -= 0.5 * ( state.ν_δx[k] * lδ_x[k] + state.ν_δx[k] * state.s0_δx[k] / exp(lδ_x[k]) ) # Jacobian built in
    end
    # J = sum(lδ_x) # Jacobian for log(δ) transformation
    # Jacobian now built in

    Q = Q_μ + Q_β + Q_δ # + J
    return Q
end
function lG0_ηlδx(μ_x::T, lδ_x::T,
                  state::State_BNP_WMReg_Joint) where T <: Real
    ## K = 1 case
    Q_μ = - 0.5 * Matrix(state.Λ0_μx)[1] * (μ_x - state.μ0_μx[1])^2.0
    Q_δ = - 0.5 * ( state.ν_δx[1] * lδ_x + state.ν_δx[1] * state.s0_δx[1] / exp(lδ_x) ) # Jacobian built in

    Q = Q_μ + Q_δ
    return Q
end

### log-collapsed conditional
function lcc_ηlδx(h::Int, indx_h::Array{Int,1}, lNX_mat::Array{T,2}, lωNX_vec::Array{T,1},
    state::State_BNP_WMReg_Joint, μ_x_h::Array{T,1}, β_x_h::Union{Array{Array{T,1},1}, Nothing},
    lδ_x_h::Array{T,1},
    Λ1star_ηy_h::PDMat{T}, a1_δy::T, b1_δy::T) where T <: Real

    # Q1
    Q1 = sum(lNX_mat[indx_h, h])

    # Q2
    Q2 = sum(lωNX_vec)

    # Q3 , G0 including Jacobian
    Q3 = lG0_ηlδx(μ_x_h, β_x_h, lδ_x_h, state)

    # Q4
    Q4 = logdet(Λ1star_ηy_h)

    # Q5
    Q5 = log(b1_δy)

    return Q1 - Q2 + Q3 - 0.5 * Q4 - (a1_δy) * Q5

end
function lcc_ηlδx(h::Int, indx_h::Array{Int,1}, lNX_mat::Array{T,2}, lωNX_vec::Array{T,1},
    state::State_BNP_WMReg_Joint, μ_x_h::T, lδ_x_h::T,
    Λ1star_ηy_h::PDMat{T}, a1_δy::T, b1_δy::T) where T <: Real

    ## K = 1 case

    # Q1
    Q1 = sum(lNX_mat[indx_h, h])

    # Q2
    Q2 = sum(lωNX_vec)

    # Q3 , G0 including Jacobian
    Q3 = lG0_ηlδx(μ_x_h, lδ_x_h, state)

    # Q4
    Q4 = logdet(Λ1star_ηy_h)

    # Q5
    Q5 = log(b1_δy)

    return Q1 - Q2 + Q3 - 0.5 * Q4 - (a1_δy) * Q5

end

## variance-inflation method
function βδ_x_h_modify_γ(β_x::Array{Array{T,1},1}, δ_x::Array{T,1},
                        γ::BitArray{1}, γδc::Array{T,1}) where T <: Real
    modify_indx = findall(.!(γ))
    K = length(γ)

    βout = deepcopy(β_x)
    δout = deepcopy(δ_x)

    for k in modify_indx # this works even if no modifications are necessary
        δout[k] += γδc[k]
    end

    for k = 1:(K - 1)
        if !γ[k]
            βout[k] *= 0.0
        else
            modify2 = intersect((k + 1):K, modify_indx)
            βout[k][(modify2 .- k)] *= 0.0
        end
    end

    return βout, δout
end
function δ_x_h_modify_γ(δ_x::Array{T,1},
                      γ::BitArray{1}, γδc::Array{T,1}) where T <: Real
    modify_indx = findall(.!(γ))
    δout = deepcopy(δ_x)

    for k in modify_indx # this works even if no modifications are necessary
        δout[k] += γδc[k]
    end
    return δout
end

## subset method
function βδ_x_h_modify_γ(β_x::Array{Array{T,1},1}, δ_x::Array{T,1},
                        γ::BitArray{1}, γδc::Float64) where T <: Real

    γδc == Inf || throw("A single variance inflation should be equal to Inf")

    γindx = findall(γ)
    nγ = length(γindx)
    nγ > 1 || throw("βδ_x_h_modify_γ requires more than one selected variable.")

    βout = [ deepcopy(β_x[γindx[k]][(γindx[(k + 1):nγ] .- γindx[k]) ])  for k = 1:(nγ - 1) ] # vector of (nγ-1):1 length vectors
    δout = deepcopy(δ_x[γindx]) # H nows and sum(gamma) cols

    return βout, δout
end
function δ_x_h_modify_γ(δ_x::Array{T,1},
                      γ::BitArray{1}, γδc::Float64) where T <: Real
    γδc == Inf || throw("A single variance inflation should be equal to Inf")
    γindx = findall(γ)
    δout = deepcopy(δ_x[γindx]) # H nows and sum(gamma) cols

    return δout
end

## helper functions for marginalization
function get_a1_δy_h(ν_δy::T, n_h::Int) where T <: Real
    return 0.5 * (ν_δy + n_h)
end

## must do no more with y_h than y'y because of its use in update_alloc!()
function get_b1_δy_h(ν_δy::T, s0_δy::T, y_h::Array{T,1}, βΛβ0star_ηy::T, Λ1star_ηy_h::PDMat, β1star_ηy_h::Array{T,1}) where T <: Real
    return 0.5 * (ν_δy * s0_δy + y_h'y_h + βΛβ0star_ηy - PDMats.quad(Λ1star_ηy_h, β1star_ηy_h))
end

function get_Λ1star_ηy_h(D_h::Array{T,2}, Λ0star_ηy::PDMat) where T <: Real
    return PDMat_adj(D_h'D_h + Λ0star_ηy)
end

function get_β1star_ηy_h(Λ1star_ηy_h::PDMat, Λβ0star_ηy::Array{T,1}, D_h::Array{T,2}, y_h::Array{T,1}) where T <: Real
    return Λ1star_ηy_h \ (Λβ0star_ηy + D_h'y_h)
end

ext: jl


text:
"""
    @sync ex

Run expression `ex` and synchronize the GPU afterwards. This is a CPU-friendly
synchronization, i.e. it performs a blocking synchronization without increasing CPU load. As
such, this operation is preferred over implicit synchronization (e.g. when performing a
memory copy) for high-performance applications.

It is also useful for timing code that executes asynchronously.
"""
macro sync(ex)
    quote
        local e = CuEvent(CUDAdrv.EVENT_BLOCKING_SYNC | CUDAdrv.EVENT_DISABLE_TIMING)
        local ret = $(esc(ex))
        CUDAdrv.record(e)
        CUDAdrv.synchronize(e)
        ret
    end
end

"""
    @argout call(..., out(output_arg), ...)

Change the behavior of a function call, returning the value of an argument instead.
A common use case is to pass a newly-created Ref and immediately dereference that output:

    @argout(some_getter(Ref{Int}()))[]

If no output argument is specified, `nothing` will be returned.
Multiple output arguments return a tuple.

"""
macro argout(ex)
    Meta.isexpr(ex, :call) || throw(ArgumentError("@argout macro should be applied to a function call"))

    block = quote end

    # look for output arguments (`out(...)`)
    output_vars = []
    args = ex.args[2:end]
    for (i,arg) in enumerate(args)
        if Meta.isexpr(arg, :call) && arg.args[1] == :out
            # allocate a variable
            @gensym output_val
            push!(block.args, :($output_val = $(ex.args[i+1].args[2]))) # strip `output(...)`
            push!(output_vars, output_val)

            # replace the argument
            ex.args[i+1] = output_val
        end
    end

    # generate a return
    push!(block.args, ex)
    if isempty(output_vars)
        push!(block.args, :(nothing))
    elseif length(output_vars) == 1
        push!(block.args, :($(output_vars[1])))
    else
        push!(block.args, :(tuple($(output_vars...))))
    end

    esc(block)
end

"""
    @workspace size=getWorkspaceSize(args...) [eltyp=UInt8] [fallback=nothing] buffer->begin
      useWorkspace(workspace, sizeof(workspace))
    end

Create a GPU workspace vector with element type `eltyp` and size in number of elements (in
the default case of an UInt8 element type this equals to the amount of bytes) determined by
calling `getWorkspaceSize`, and pass it to the closure for use in calling `useWorkspace`.
A fallback workspace size `fallback` can be specified if the regular size would lead to OOM.
Afterwards, the buffer is put back into the memory pool for reuse.

This helper protects against the rare but real issue of `getWorkspaceSize` returning
different results based on the GPU device memory pressure, which might change _after_
initial allocation of the workspace (which can cause a GC collection).
"""
macro workspace(ex...)
    code = ex[end]
    kwargs = ex[1:end-1]

    sz = nothing
    eltyp = :UInt8
    fallback = nothing
    for kwarg in kwargs
        key,val = kwarg.args
        if key == :size
            sz = val
        elseif key == :eltyp
            eltyp = val
        elseif key == :fallback
            fallback = val
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end

    if sz === nothing
        throw(ArgumentError("@workspace macro needs a size argument"))
    end

    # break down the closure to a let block to prevent JuliaLang/julia#15276
    Meta.isexpr(code, :(->)) || throw(ArgumentError("@workspace macro should be applied to a closure"))
    length(code.args) == 2 || throw(ArgumentError("@workspace closure should take exactly one argument"))
    code_arg = code.args[1]
    code = code.args[2]

    return quote
        sz = $(esc(sz))
        workspace = nothing
        try
          while workspace === nothing || length(workspace) < sz
              workspace = CuArray{$(esc(eltyp))}(undef, sz)
              sz = $(esc(sz))
          end
        catch ex
            $fallback === nothing && rethrow()
            isa(ex, OutOfGPUMemoryError) || rethrow()
            workspace = CuArray{UInt8}(undef, $fallback)
        end

        let $(esc(code_arg)) = workspace
            ret = $(esc(code))
            unsafe_free!(workspace)
            ret
        end
    end
end


text:
struct OperatorTrace end
struct SearchTrace end

@recipe function operator_trace(df::DataFrame, type::Type{OperatorTrace})
    x, y = align(df.moveparams)
    seriestype --> :scatter
    x, y
end

"""
    align(Φ)

Convert parameter-tuples into (x,y)-coordinates, where x indicates the iteration number and y the variable involved.
"""
function align(Φ)
    manipulated_vars = Int[]
    iteration_count = Int[]
    cnt = 1
    for el in Φ # voor elke operator-parameter
        push!(manipulated_vars, el...) # voeg alle variabele-parameters toe
        push!(iteration_count, ntuple(x -> cnt, length(el))...) #
        cnt+=1
    end
    return iteration_count, manipulated_vars
end


text:
############################################################
# joDAparallelToggles - parallel V/MV conversion operators #
############################################################

# includes
include("joAbstractDAparallelToggleOperator/constructors.jl")
include("joAbstractDAparallelToggleOperator/base_functions.jl")
include("joAbstractDAparallelToggleOperator/extra_functions.jl")


text:
abstract Simulation
abstract Galerkin <: Simulation
abstract DG <: Galerkin
abstract CG <: Galerkin

import Base: show

# 1 Dimensional Discontinuous Galerkin Method

type DG1D <: DG
    p::Int64 # The order of ϕ
    K::Int64 # The number of cells to use
    a::Float64 # [a,b] is the interval on which to cumpute
    b::Float64 #
    fflux::Function # The numerical flux function
    Δt::Float64 # Time-step
    nt::Int64 # Number of time steps
    C::Float64 # Diffusion constant
    q₀
    mesh::Any
    boundary::Function
    hackinextraflux::Bool
    useMUSCL::Bool
    DG1D(p, K, a::Float64, b::Float64, fflux::Function, Δt::Float64,
        nt::Int64, C::Float64, q₀; mesh = nothing, boundary = (args...)->zero(typeof(q₀(0))), hackinextraflux=false, useMUSCL = false) =
        new(p,K,a,b,fflux,Δt,nt,C,q₀,mesh,boundary,hackinextraflux, useMUSCL)
end

function show(io::IO, p::DG1D)
    println(io,"1 Dimensional Discontinuous Galerkin Simulation:")
    println(io,"    Ω = [",p.a,",",p.b,"]")
    println(io," mesh = 1D Uniform (K = ",p.K,")")
    println(io,"    p = ",p.p)
    println(io,"   Δt = ",p.Δt)
    println(io,"   nt = ",p.nt)
end

porder(p::DG1D) = p.p
𝖓(p::DG1D) = p.p+1

zerop(p::DG1D) = zero(Float64)
numcells(p::DG1D) = p.K

# 2 Dimensional Discontinuous Galerkin Method

type DG2D <: DG
    mesh::Meshes.Mesh{Vertex2}
    porder::Int64
    fflux::Function # The numerical flux function
    Δt::Float64 # Time-step
    nt::Int64 # Number of time steps
    C::Float64 # Diffusion constant
    q₀
    boundary::Function
end

type CG2D{porder} <: DG
    mesh::Meshes.Mesh{Vertex2}
    dualmesh::Meshes.Mesh{Vertex2}
end

const Galerkin2D = Union(CG2D,DG2D)

function show(io::IO, p::Galerkin2D)
    println(io,"2D Galerkin Simulation:")
    println(io," mesh = 2D Triangular")
    println(io,"    p = ",porder(p))
end

numcells(p::DG2D) = length(p.mesh.faces)
porder(p::DG2D) = p.porder
porder{pp}(::CG2D{pp}) = pp

zerop(p::Galerkin2D) = zero(Vertex2)

# sum_{k=0}^p (k + n - 1 \choose n - 1) for n = 2
𝖓(p::Galerkin2D) = div((porder(p)+1)*(porder(p)+2),2)

# Utilities
const nbf = 𝖓
mesh(p::Galerkin) = p.mesh


text:
@snn_kw struct PoissonParameter{FT=Float32}
    rate::FT = 1Hz
end

@snn_kw mutable struct Poisson{VFT=Vector{Float32},VBT=Vector{Bool}}
    param::PoissonParameter = PoissonParameter()
    N::Int32 = 100
    randcache::VFT = rand(N)
    fire::VBT = zeros(Bool, N)
    records::Dict = Dict()
end

"""
[Poisson Neuron](https://www.cns.nyu.edu/~david/handouts/poisson.pdf)
"""
Poisson

function integrate!(p::Poisson, param::PoissonParameter, dt::Float32)
    @unpack N, randcache, fire = p
    @unpack rate = param
    prob = rate * dt
    rand!(randcache)
    @inbounds for i = 1:N
        fire[i] = randcache[i] < prob
    end
end


text:
struct SingleThreaded end
struct MultiThreaded end
struct GPU end
Processing = Union{SingleThreaded, MultiThreaded, GPU}
struct InvertedIndex end
struct BMatrix end
UpdateMethod = Union{InvertedIndex, BMatrix}
struct Exact end
struct Nesterov end
OptimizationMethod = Union{Exact, Nesterov}

const DEFAULT_CONFIG = (n_codebooks=0,
                        n_centers=8,
                        a=0,
                        b=0,
                        initialise_with_l2_loss=true,
                        T_preclustering=-1,
                        inverted_index=false,
                        max_iter=1000,
                        stopcond=1e-2,
                        verbose=false,
                        max_iter_assignments=10,
                        optimisation="exact",
                        multithreading=false,
                        GPU=false,
                        increment_steps=0,
                        incremental_fitting_iteration=false,
                        reorder                = 0,
                        training_points        = 0)

function generate_data_dependent_defaults(n_dp::Int, n_dims::Int)
    a = floor(Int,sqrt(n_dp))
    training_points = Int(ceil(n_dp/5))
    config_updated  =   (n_codebooks = n_dims ÷ 2,
                        a = a,
                        b = (a ÷ 5) + 1,
                        training_points        = training_points,
                        increment_steps=Int(floor(log10(training_points))-2))
    config = merge(DEFAULT_CONFIG, config_updated)
    return config
end

function check_kwargs(kwargs, n_dp, n_dims, T)
    if (compute_η(T, n_dims) < 1.01) & (compute_η(T, n_dims) > .99)
        @warn("Anisotropic Loss with T=$T for $n_dims dimensions is equivalent to L2, updated loss function to L2 loss.")
        T = 0
    end
    for key in keys(kwargs)
        if key ∉ keys(DEFAULT_CONFIG)
            error("$key is not a valid configuration keyword")
        end
    end
    config = generate_data_dependent_defaults(n_dp, n_dims)
    config = merge(config, kwargs)
    if n_dims % config.n_codebooks != 0
        error("Number of dimensions of data should be divisable by number of codebooks.
               Make sure to set the right n_codebook in the builder.")
    end
    tp = if config.training_points > 0 config.training_points else n_dp end
    if config.increment_steps > 0
        smallest_sample = collect(logrange(tp, config.increment_steps))[1]
        if !(smallest_sample > config.a)
            @warn("Lowest subsample size for incremental fitting is lower than number of centers. Increment_steps have been decreased.")
            config = merge(config, (; increment_steps = Int(ceil(log10(tp))-ceil(log10(config.a)))))
        end
    end
    if tp < config.n_centers
        @error("$tp Training points is lower than $(config.n_centers)")
    end
    if (config.GPU) & (config.optimisation ≠ "Nesterov")
        config = merge(config, (;optimisation="Nesterov"))
        @warn("GPU processing was selected, codebook optimisation has automatically changed to approximate method.")
    end
    if config.T_preclustering < 0
        config = merge(config, (;T_preclustering=T))
    end
    return config, T
end

function generate_config_vars(config)
    thread = if config.multithreading MultiThreaded() else SingleThreaded() end
    update_method = if config.inverted_index InvertedIndex() else BMatrix() end
    optim_method = if lowercase(config.optimisation)[1]=='e' Exact() else Nesterov() end
    processing = if config.GPU GPU() else thread end
    return thread, update_method, optim_method, processing
end

struct AnisotropicWeights{F<:AbstractFloat}
    η::F
    h_orthog::F
    h_par::F
end

struct L2_loss end
Weights = Union{AnisotropicWeights, L2_loss}

compute_η(T::Real, n_dim::Int) = (n_dim-1)*T^2/(1-T^2)

function ComputeWeightsFromT(n_dim::Int, T::AbstractFloat)
    η = compute_η(T, n_dim)
    h_orthog = 1/(1+η)
    h_par    = 1 - h_orthog
    return AnisotropicWeights(η, h_orthog, h_par)
end


text:
## Constants

# Context options
const IO_THREADS = 1
const MAX_SOCKETS = 2
const IPV6 = 42

#Socket Types
const PAIR = 0
const PUB = 1
const SUB = 2
const REQ = 3
const REP = 4
const DEALER = 5
const ROUTER = 6
const PULL = 7
const PUSH = 8
const XPUB = 9
const XSUB = 10
const XREQ = DEALER
const XREP = ROUTER
const UPSTREAM = PULL
const DOWNSTREAM = PUSH

#Message options
const MORE = 1
const SNDMORE = true

#IO Multiplexing
const POLLIN = 1
const POLLOUT = 2
const POLLERR = 4

#Built in devices
const STREAMER = 1
const FORWARDER = 2
const QUEUE = 3

#Send/Recv Options
const ZMQ_DONTWAIT = 1
const ZMQ_SNDMORE = 2


text:
using JuliaTutor
using Documenter

@info "Makeing documentation..."
makedocs(;
    modules=[JuliaTutor],
    authors="caseykneale",
    repo="https://github.com/Humans-of-Julia/JuliaTutor.jl/blob/{commit}{path}#L{line}",
    sitename="JuliaTutor.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://Humans-of-Julia.github.io/JuliaTutor.jl",
        assets=["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Humans-of-Julia/JuliaTutor.jl",
    devbranch="master",
)


text:
# SPDX-License-Identifier: X11
# 2020-11-25
# Traveling Salesman around Lake (300pt)

function main()
  K,N = parse.(Int, (split ∘ readline)())
  A = parse.(Int, (split ∘ readline)())

  dist = -1
  for i ∈ 1:N - 1
    dist = max(dist, abs(A[i] - A[i + 1]))
  end
  dist = max(dist, (K - A[N]) + A[1])

  println(K - dist)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
ztype(::AbstractOp, intype::DataType) = intype

#UnaryOps:
ztype(::UnaryOps.ISINF_T, ::DataType) = Bool
ztype(::UnaryOps.ISNAN_T, ::DataType) = Bool
ztype(::UnaryOps.ISFINITE_T, ::DataType) = Bool

ztype(::UnaryOps.CONJ_T, intype::Type{T}) where {T <: Complex} = intype.parameters[1]
ztype(::UnaryOps.ABS_T, intype::Type{T}) where {T <: Complex} = intype.parameters[1]
ztype(::UnaryOps.CREAL_T, intype::Type{T}) where {T <: Complex} = intype.parameters[1]
ztype(::UnaryOps.CIMAG_T, intype::Type{T}) where {T <: Complex} = intype.parameters[1]
ztype(::UnaryOps.CARG_T, intype::Type{T}) where {T <: Complex} = intype.parameters[1]

ztype(::UnaryOps.POSITIONI_T, ::DataType) = Int64
ztype(::UnaryOps.POSITIONI1_T, ::DataType) = Int64
ztype(::UnaryOps.POSITIONJ_T, ::DataType) = Int64
ztype(::UnaryOps.POSITIONJ1_T, ::DataType) = Int64

#BinaryOps:
ztype(::BinaryOps.EQ_T, ::DataType) = Bool
ztype(::BinaryOps.NE_T, ::DataType) = Bool
ztype(::BinaryOps.GT_T, ::DataType) = Bool
ztype(::BinaryOps.LT_T, ::DataType) = Bool
ztype(::BinaryOps.GE_T, ::DataType) = Bool
ztype(::BinaryOps.LE_T, ::DataType) = Bool
ztype(::BinaryOps.CMPLX_T, intype::Type{T}) where {T <: AbstractFloat} = Complex{T}

ztype(::BinaryOps.FIRSTI_T, ::DataType) = Int64
ztype(::BinaryOps.FIRSTI1_T, ::DataType) = Int64
ztype(::BinaryOps.FIRSTJ_T, ::DataType) = Int64
ztype(::BinaryOps.FIRSTJ1_T, ::DataType) = Int64
ztype(::BinaryOps.SECONDI_T, ::DataType) = Int64
ztype(::BinaryOps.SECONDI1_T, ::DataType) = Int64
ztype(::BinaryOps.SECONDJ_T, ::DataType) = Int64
ztype(::BinaryOps.SECONDJ1_T, ::DataType) = Int64

#Semirings:
ztype(::Semirings.LAND_EQ_T, ::DataType) = Bool
ztype(::Semirings.LOR_EQ_T, ::DataType) = Bool
ztype(::Semirings.LXOR_EQ_T, ::DataType) = Bool
ztype(::Semirings.EQ_EQ_T, ::DataType) = Bool
ztype(::Semirings.ANY_EQ_T, ::DataType) = Bool
ztype(::Semirings.LAND_NE_T, ::DataType) = Bool
ztype(::Semirings.LOR_NE_T, ::DataType) = Bool
ztype(::Semirings.LXOR_NE_T, ::DataType) = Bool
ztype(::Semirings.EQ_NE_T, ::DataType) = Bool
ztype(::Semirings.ANY_NE_T, ::DataType) = Bool
ztype(::Semirings.LAND_GT_T, ::DataType) = Bool
ztype(::Semirings.LOR_GT_T, ::DataType) = Bool
ztype(::Semirings.LXOR_GT_T, ::DataType) = Bool
ztype(::Semirings.EQ_GT_T, ::DataType) = Bool
ztype(::Semirings.ANY_GT_T, ::DataType) = Bool
ztype(::Semirings.LAND_LT_T, ::DataType) = Bool
ztype(::Semirings.LOR_LT_T, ::DataType) = Bool
ztype(::Semirings.LXOR_LT_T, ::DataType) = Bool
ztype(::Semirings.EQ_LT_T, ::DataType) = Bool
ztype(::Semirings.ANY_LT_T, ::DataType) = Bool
ztype(::Semirings.LAND_GE_T, ::DataType) = Bool
ztype(::Semirings.LOR_GE_T, ::DataType) = Bool
ztype(::Semirings.LXOR_GE_T, ::DataType) = Bool
ztype(::Semirings.EQ_GE_T, ::DataType) = Bool
ztype(::Semirings.ANY_GE_T, ::DataType) = Bool
ztype(::Semirings.LAND_LE_T, ::DataType) = Bool
ztype(::Semirings.LOR_LE_T, ::DataType) = Bool
ztype(::Semirings.LXOR_LE_T, ::DataType) = Bool
ztype(::Semirings.EQ_LE_T, ::DataType) = Bool
ztype(::Semirings.ANY_LE_T, ::DataType) = Bool

ztype(::Semirings.MIN_FIRSTI_T, ::DataType) = Int64
ztype(::Semirings.MAX_FIRSTI_T, ::DataType) = Int64
ztype(::Semirings.PLUS_FIRSTI_T, ::DataType) = Int64
ztype(::Semirings.TIMES_FIRSTI_T, ::DataType) = Int64
ztype(::Semirings.ANY_FIRSTI_T, ::DataType) = Int64
ztype(::Semirings.MIN_FIRSTI1_T, ::DataType) = Int64
ztype(::Semirings.MAX_FIRSTI1_T, ::DataType) = Int64
ztype(::Semirings.PLUS_FIRSTI1_T, ::DataType) = Int64
ztype(::Semirings.TIMES_FIRSTI1_T, ::DataType) = Int64
ztype(::Semirings.ANY_FIRSTI1_T, ::DataType) = Int64
ztype(::Semirings.MIN_FIRSTJ_T, ::DataType) = Int64
ztype(::Semirings.MAX_FIRSTJ_T, ::DataType) = Int64
ztype(::Semirings.PLUS_FIRSTJ_T, ::DataType) = Int64
ztype(::Semirings.TIMES_FIRSTJ_T, ::DataType) = Int64
ztype(::Semirings.ANY_FIRSTJ_T, ::DataType) = Int64
ztype(::Semirings.MIN_FIRSTJ1_T, ::DataType) = Int64
ztype(::Semirings.MAX_FIRSTJ1_T, ::DataType) = Int64
ztype(::Semirings.PLUS_FIRSTJ1_T, ::DataType) = Int64
ztype(::Semirings.TIMES_FIRSTJ1_T, ::DataType) = Int64
ztype(::Semirings.ANY_FIRSTJ1_T, ::DataType) = Int64
ztype(::Semirings.MIN_SECONDI_T, ::DataType) = Int64
ztype(::Semirings.MAX_SECONDI_T, ::DataType) = Int64
ztype(::Semirings.PLUS_SECONDI_T, ::DataType) = Int64
ztype(::Semirings.TIMES_SECONDI_T, ::DataType) = Int64
ztype(::Semirings.ANY_SECONDI_T, ::DataType) = Int64
ztype(::Semirings.MIN_SECONDI1_T, ::DataType) = Int64
ztype(::Semirings.MAX_SECONDI1_T, ::DataType) = Int64
ztype(::Semirings.PLUS_SECONDI1_T, ::DataType) = Int64
ztype(::Semirings.TIMES_SECONDI1_T, ::DataType) = Int64
ztype(::Semirings.ANY_SECONDI1_T, ::DataType) = Int64
ztype(::Semirings.MIN_SECONDJ_T, ::DataType) = Int64
ztype(::Semirings.MAX_SECONDJ_T, ::DataType) = Int64
ztype(::Semirings.PLUS_SECONDJ_T, ::DataType) = Int64
ztype(::Semirings.TIMES_SECONDJ_T, ::DataType) = Int64
ztype(::Semirings.ANY_SECONDJ_T, ::DataType) = Int64
ztype(::Semirings.MIN_SECONDJ1_T, ::DataType) = Int64
ztype(::Semirings.MAX_SECONDJ1_T, ::DataType) = Int64
ztype(::Semirings.PLUS_SECONDJ1_T, ::DataType) = Int64
ztype(::Semirings.TIMES_SECONDJ1_T, ::DataType) = Int64
ztype(::Semirings.ANY_SECONDJ1_T, ::DataType) = Int64


text:
using Genie.Router
ENV["GKSwstype"]="nul"
using Analysis,JuliaDB,StatsPlots,GLM

t = Analysis.load_csse_data("public/data")

route("/") do
    print("changed")
    country = @params(:country)
    days = parse(Int,@params(:days))
    chartfilename = "static/$country-$days-chart.png"
    (regress,chart) = Analysis.exp_regress(Analysis.last_n_days(days,Analysis.filter_country(country,t)))
    savefig(chart,"public/$chartfilename")
    "

      In $country, for the last $days days the number of Corona cases has been doubling every $(round(1 / coef(regress)[2],digits=1)) days

"
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module PeriodsTest

using Dates
using Test

@testset "basic arithmetic" begin
    @test -DatesPlus.Year(1) == DatesPlus.Year(-1)
    @test DatesPlus.Year(1) > DatesPlus.Year(0)
    @test (DatesPlus.Year(1) < DatesPlus.Year(0)) == false
    @test DatesPlus.Year(1) == DatesPlus.Year(1)
    @test DatesPlus.Year(1) != 1
    @test DatesPlus.Year(1) + DatesPlus.Year(1) == DatesPlus.Year(2)
    @test DatesPlus.Year(1) - DatesPlus.Year(1) == zero(DatesPlus.Year)
    @test 1 == one(DatesPlus.Year)
    @test_throws MethodError DatesPlus.Year(1) * DatesPlus.Year(1) == DatesPlus.Year(1)
    t = DatesPlus.Year(1)
    t2 = DatesPlus.Year(2)
    @test ([t, t, t, t, t] .+ DatesPlus.Year(1)) == ([t2, t2, t2, t2, t2])
    @test (DatesPlus.Year(1) .+ [t, t, t, t, t]) == ([t2, t2, t2, t2, t2])
    @test ([t2, t2, t2, t2, t2] .- DatesPlus.Year(1)) == ([t, t, t, t, t])
    @test_throws MethodError ([t, t, t, t, t] .* DatesPlus.Year(1)) == ([t, t, t, t, t])
    @test ([t, t, t, t, t] * 1) == ([t, t, t, t, t])
    @test ([t, t, t, t, t] .% t2) == ([t, t, t, t, t])
    @test div.([t, t, t, t, t], DatesPlus.Year(1)) == ([1, 1, 1, 1, 1])
    @test mod.([t, t, t, t, t], DatesPlus.Year(2)) == ([t, t, t, t, t])
    @test [t, t, t] / t2 == [0.5, 0.5, 0.5]
    @test abs(-t) == t
    @test sign(t) == sign(t2) == 1
    @test sign(-t) == sign(-t2) == -1
    @test sign(DatesPlus.Year(0)) == 0
end
@testset "div/mod/gcd/lcm/rem" begin
    @test DatesPlus.Year(10) % DatesPlus.Year(4) == DatesPlus.Year(2)
    @test gcd(DatesPlus.Year(10), DatesPlus.Year(4)) == DatesPlus.Year(2)
    @test lcm(DatesPlus.Year(10), DatesPlus.Year(4)) == DatesPlus.Year(20)
    @test div(DatesPlus.Year(10), DatesPlus.Year(3)) == 3
    @test div(DatesPlus.Year(10), DatesPlus.Year(4)) == 2
    @test div(DatesPlus.Year(10), 4) == DatesPlus.Year(2)
    @test DatesPlus.Year(10) / DatesPlus.Year(4) == 2.5

    @test mod(DatesPlus.Year(10), DatesPlus.Year(4)) == DatesPlus.Year(2)
    @test mod(DatesPlus.Year(-10), DatesPlus.Year(4)) == DatesPlus.Year(2)
    @test mod(DatesPlus.Year(10), 4) == DatesPlus.Year(2)
    @test mod(DatesPlus.Year(-10), 4) == DatesPlus.Year(2)

    @test rem(DatesPlus.Year(10), DatesPlus.Year(4)) == DatesPlus.Year(2)
    @test rem(DatesPlus.Year(-10), DatesPlus.Year(4)) == DatesPlus.Year(-2)
    @test rem(DatesPlus.Year(10), 4) == DatesPlus.Year(2)
    @test rem(DatesPlus.Year(-10), 4) == DatesPlus.Year(-2)
end

y = DatesPlus.Year(1)
q = DatesPlus.Quarter(1)
m = DatesPlus.Month(1)
w = DatesPlus.Week(1)
d = DatesPlus.Day(1)
h = DatesPlus.Hour(1)
mi = DatesPlus.Minute(1)
s = DatesPlus.Second(1)
ms = DatesPlus.Millisecond(1)
us = DatesPlus.Microsecond(1)
ns = DatesPlus.Nanosecond(1)
emptyperiod = ((y + d) - d) - y
@testset "Period arithmetic" begin
    @test DatesPlus.Year(y) == y
    @test DatesPlus.Quarter(q) == q
    @test DatesPlus.Month(m) == m
    @test DatesPlus.Week(w) == w
    @test DatesPlus.Day(d) == d
    @test DatesPlus.Hour(h) == h
    @test DatesPlus.Minute(mi) == mi
    @test DatesPlus.Second(s) == s
    @test DatesPlus.Millisecond(ms) == ms
    @test DatesPlus.Microsecond(us) == us
    @test DatesPlus.Nanosecond(ns) == ns
    @test DatesPlus.Year(convert(Int8, 1)) == y
    @test DatesPlus.Year(convert(UInt8, 1)) == y
    @test DatesPlus.Year(convert(Int16, 1)) == y
    @test DatesPlus.Year(convert(UInt16, 1)) == y
    @test DatesPlus.Year(convert(Int32, 1)) == y
    @test DatesPlus.Year(convert(UInt32, 1)) == y
    @test DatesPlus.Year(convert(Int64, 1)) == y
    @test DatesPlus.Year(convert(UInt64, 1)) == y
    @test DatesPlus.Year(convert(Int128, 1)) == y
    @test DatesPlus.Year(convert(UInt128, 1)) == y
    @test DatesPlus.Year(convert(BigInt, 1)) == y
    @test DatesPlus.Year(convert(BigFloat, 1)) == y
    @test DatesPlus.Year(convert(Complex, 1)) == y
    @test DatesPlus.Year(convert(Rational, 1)) == y
    @test DatesPlus.Year(convert(Float16, 1)) == y
    @test DatesPlus.Year(convert(Float32, 1)) == y
    @test DatesPlus.Year(convert(Float64, 1)) == y
    @test y == y
    @test m == m
    @test w == w
    @test d == d
    @test h == h
    @test mi == mi
    @test s == s
    @test ms == ms
    @test us == us
    @test ns == ns
    y2 = DatesPlus.Year(2)
    @test y < y2
    @test y2 > y
    @test y != y2

    @test DatesPlus.Year(Int8(1)) == y
    @test DatesPlus.Year(UInt8(1)) == y
    @test DatesPlus.Year(Int16(1)) == y
    @test DatesPlus.Year(UInt16(1)) == y
    @test DatesPlus.Year(Int(1)) == y
    @test DatesPlus.Year(UInt(1)) == y
    @test DatesPlus.Year(Int64(1)) == y
    @test DatesPlus.Year(UInt64(1)) == y
    @test DatesPlus.Year(UInt128(1)) == y
    @test DatesPlus.Year(UInt128(1)) == y
    @test DatesPlus.Year(big(1)) == y
    @test DatesPlus.Year(BigFloat(1)) == y
    @test DatesPlus.Year(float(1)) == y
    @test DatesPlus.Year(Float32(1)) == y
    @test DatesPlus.Year(Rational(1)) == y
    @test DatesPlus.Year(complex(1)) == y
    @test_throws InexactError DatesPlus.Year(BigFloat(1.2)) == y
    @test_throws InexactError DatesPlus.Year(1.2) == y
    @test_throws InexactError DatesPlus.Year(Float32(1.2)) == y
    @test_throws InexactError DatesPlus.Year(3//4) == y
    @test_throws InexactError DatesPlus.Year(complex(1.2)) == y
    @test_throws InexactError DatesPlus.Year(Float16(1.2)) == y
    @test DatesPlus.Year(true) == y
    @test DatesPlus.Year(false) != y
    @test_throws MethodError DatesPlus.Year(:hey) == y
    @test DatesPlus.Year(real(1)) == y
    @test_throws InexactError DatesPlus.Year(m) == y
    @test_throws MethodError DatesPlus.Year(w) == y
    @test_throws MethodError DatesPlus.Year(d) == y
    @test_throws MethodError DatesPlus.Year(h) == y
    @test_throws MethodError DatesPlus.Year(mi) == y
    @test_throws MethodError DatesPlus.Year(s) == y
    @test_throws MethodError DatesPlus.Year(ms) == y
    @test DatesPlus.Year(DatesPlus.Date(2013, 1, 1)) == DatesPlus.Year(2013)
    @test DatesPlus.Year(DatesPlus.DateTime(2013, 1, 1)) == DatesPlus.Year(2013)
    @test typeof(y + m) <: DatesPlus.CompoundPeriod
    @test typeof(m + y) <: DatesPlus.CompoundPeriod
    @test typeof(y + w) <: DatesPlus.CompoundPeriod
    @test typeof(y + d) <: DatesPlus.CompoundPeriod
    @test typeof(y + h) <: DatesPlus.CompoundPeriod
    @test typeof(y + mi) <: DatesPlus.CompoundPeriod
    @test typeof(y + s) <: DatesPlus.CompoundPeriod
    @test typeof(y + ms) <: DatesPlus.CompoundPeriod
    @test typeof(y + us) <: DatesPlus.CompoundPeriod
    @test typeof(y + ns) <: DatesPlus.CompoundPeriod
    @test y > m
    @test d < w
    @test mi < h
    @test ms < h
    @test ms < mi
    @test us < ms
    @test ns < ms
    @test ns < us
    @test ns < w
    @test us < w
    @test typemax(DatesPlus.Year) == DatesPlus.Year(typemax(Int64))
    @test typemax(DatesPlus.Year) + y == DatesPlus.Year(-9223372036854775808)
    @test typemin(DatesPlus.Year) == DatesPlus.Year(-9223372036854775808)
end
@testset "Period-Real arithmetic" begin
    @test_throws MethodError y + 1 == DatesPlus.Year(2)
    @test_throws MethodError y + true == DatesPlus.Year(2)
    @test_throws InexactError y + DatesPlus.Year(1.2)
    @test y + DatesPlus.Year(1f0) == DatesPlus.Year(2)
    @test y * 4 == DatesPlus.Year(4)
    @test y * 4f0 == DatesPlus.Year(4)
    @test DatesPlus.Year(2) * 0.5 == y
    @test DatesPlus.Year(2) * 3//2 == DatesPlus.Year(3)
    @test_throws InexactError y * 0.5
    @test_throws InexactError y * 3//4
    @test (1:1:5)*Second(5) === Second(5)*(1:1:5) === Second(5):Second(5):Second(25) === (1:5)*Second(5)
    @test collect(1:1:5)*Second(5) == Second(5)*collect(1:1:5) == (1:5)*Second(5)
    @test (Second(2):Second(2):Second(10))/Second(2) === 1.0:1.0:5.0 == collect(Second(2):Second(2):Second(10))/Second(2)
    @test (Second(2):Second(2):Second(10)) / 2 == Second(1):Second(1):Second(5) == collect(Second(2):Second(2):Second(10)) / 2
    @test DatesPlus.Year(4) / 2 == DatesPlus.Year(2)
    @test DatesPlus.Year(4) / 2f0 == DatesPlus.Year(2)
    @test DatesPlus.Year(4) / 0.5 == DatesPlus.Year(8)
    @test DatesPlus.Year(4) / 2//3 == DatesPlus.Year(6)
    @test_throws InexactError DatesPlus.Year(4) / 3.0
    @test_throws InexactError DatesPlus.Year(4) / 3//2
    @test div(y, 2) == DatesPlus.Year(0)
    @test_throws MethodError div(2, y) == DatesPlus.Year(2)
    @test div(y, y) == 1
    @test y*10 % DatesPlus.Year(5) == DatesPlus.Year(0)
    @test_throws MethodError (y > 3) == false
    @test_throws MethodError (4 < y) == false
    @test 1 != y
    t = [y, y, y, y, y]
    @test t .+ DatesPlus.Year(2) == [DatesPlus.Year(3), DatesPlus.Year(3), DatesPlus.Year(3), DatesPlus.Year(3), DatesPlus.Year(3)]

    let x = DatesPlus.Year(5), y = DatesPlus.Year(2)
        @test div(x, y) * y + rem(x, y) == x
        @test fld(x, y) * y + mod(x, y) == x
    end
end
@testset "Associativity" begin
    dt = DatesPlus.DateTime(2012, 12, 21)
    test = ((((((((dt + y) - m) + w) - d) + h) - mi) + s) - ms)
    @test test == dt + y - m + w - d + h - mi + s - ms
    @test test == y - m + w - d + dt + h - mi + s - ms
    @test test == dt - m + y - d + w - mi + h - ms + s
    @test test == dt + (y - m + w - d + h - mi + s - ms)
    @test test == dt + y - m + w - d + (h - mi + s - ms)
    @test (dt + DatesPlus.Year(4)) + DatesPlus.Day(1) == dt + (DatesPlus.Year(4) + DatesPlus.Day(1))
    @test DatesPlus.Date(2014, 1, 29) + DatesPlus.Month(1) + DatesPlus.Day(1) + DatesPlus.Month(1) + DatesPlus.Day(1) ==
        DatesPlus.Date(2014, 1, 29) + DatesPlus.Day(1) + DatesPlus.Month(1) + DatesPlus.Month(1) + DatesPlus.Day(1)
    @test DatesPlus.Date(2014, 1, 29) + DatesPlus.Month(1) + DatesPlus.Day(1) == DatesPlus.Date(2014, 1, 29) + DatesPlus.Day(1) + DatesPlus.Month(1)
end
@testset "traits" begin
    @test DatesPlus._units(DatesPlus.Year(0)) == " years"
    @test DatesPlus._units(DatesPlus.Year(1)) == " year"
    @test DatesPlus._units(DatesPlus.Year(-1)) == " year"
    @test DatesPlus._units(DatesPlus.Year(2)) == " years"
    @test DatesPlus.string(DatesPlus.Year(0)) == "0 years"
    @test DatesPlus.string(DatesPlus.Year(1)) == "1 year"
    @test DatesPlus.string(DatesPlus.Year(-1)) == "-1 year"
    @test DatesPlus.string(DatesPlus.Year(2)) == "2 years"
    @test isfinite(DatesPlus.Year)
    @test isfinite(DatesPlus.Year(0))
    @test zero(DatesPlus.Year) == DatesPlus.Year(0)
    @test zero(DatesPlus.Year(10)) == DatesPlus.Year(0)
    @test zero(DatesPlus.Month) == DatesPlus.Month(0)
    @test zero(DatesPlus.Month(10)) == DatesPlus.Month(0)
    @test zero(DatesPlus.Day) == DatesPlus.Day(0)
    @test zero(DatesPlus.Day(10)) == DatesPlus.Day(0)
    @test zero(DatesPlus.Hour) == DatesPlus.Hour(0)
    @test zero(DatesPlus.Hour(10)) == DatesPlus.Hour(0)
    @test zero(DatesPlus.Minute) == DatesPlus.Minute(0)
    @test zero(DatesPlus.Minute(10)) == DatesPlus.Minute(0)
    @test zero(DatesPlus.Second) == DatesPlus.Second(0)
    @test zero(DatesPlus.Second(10)) == DatesPlus.Second(0)
    @test zero(DatesPlus.Millisecond) == DatesPlus.Millisecond(0)
    @test zero(DatesPlus.Millisecond(10)) == DatesPlus.Millisecond(0)
    @test DatesPlus.Year(-1) < DatesPlus.Year(1)
    @test !(DatesPlus.Year(-1) > DatesPlus.Year(1))
    @test DatesPlus.Year(1) == DatesPlus.Year(1)
    @test DatesPlus.Year(1) != 1
    @test 1 != DatesPlus.Year(1)
    @test DatesPlus.Month(-1) < DatesPlus.Month(1)
    @test !(DatesPlus.Month(-1) > DatesPlus.Month(1))
    @test DatesPlus.Month(1) == DatesPlus.Month(1)
    @test DatesPlus.Day(-1) < DatesPlus.Day(1)
    @test !(DatesPlus.Day(-1) > DatesPlus.Day(1))
    @test DatesPlus.Day(1) == DatesPlus.Day(1)
    @test DatesPlus.Hour(-1) < DatesPlus.Hour(1)
    @test !(DatesPlus.Hour(-1) > DatesPlus.Hour(1))
    @test DatesPlus.Hour(1) == DatesPlus.Hour(1)
    @test DatesPlus.Minute(-1) < DatesPlus.Minute(1)
    @test !(DatesPlus.Minute(-1) > DatesPlus.Minute(1))
    @test DatesPlus.Minute(1) == DatesPlus.Minute(1)
    @test DatesPlus.Second(-1) < DatesPlus.Second(1)
    @test !(DatesPlus.Second(-1) > DatesPlus.Second(1))
    @test DatesPlus.Second(1) == DatesPlus.Second(1)
    @test DatesPlus.Millisecond(-1) < DatesPlus.Millisecond(1)
    @test !(DatesPlus.Millisecond(-1) > DatesPlus.Millisecond(1))
    @test DatesPlus.Millisecond(1) == DatesPlus.Millisecond(1)
    @test_throws MethodError DatesPlus.Year(1) < DatesPlus.Millisecond(1)
    @test_throws MethodError DatesPlus.Millisecond(1) < DatesPlus.Year(1)

    # issue #27076
    @test DatesPlus.Year(1) != DatesPlus.Millisecond(1)
    @test DatesPlus.Millisecond(1) != DatesPlus.Year(1)
end

struct Beat <: DatesPlus.Period
    value::Int64
end

Beat(p::Period) = Beat(DatesPlus.toms(p) ÷ 86400)

@testset "comparisons with new subtypes of Period" begin
    # https://en.wikipedia.org/wiki/Swatch_Internet_Time
    DatesPlus.value(b::Beat) = b.value
    DatesPlus.toms(b::Beat) = DatesPlus.value(b) * 86400
    DatesPlus._units(b::Beat) = " beat" * (abs(DatesPlus.value(b)) == 1 ? "" : "s")
    Base.promote_rule(::Type{DatesPlus.Day}, ::Type{Beat}) = DatesPlus.Millisecond
    Base.convert(::Type{T}, b::Beat) where {T<:DatesPlus.Millisecond} = T(DatesPlus.toms(b))

    @test Beat(1000) == DatesPlus.Day(1)
    @test Beat(1) < DatesPlus.Day(1)
end

@testset "basic properties" begin
    @test DatesPlus.Year("1") == y
    @test DatesPlus.Quarter("1") == q
    @test DatesPlus.Month("1") == m
    @test DatesPlus.Week("1") == w
    @test DatesPlus.Day("1") == d
    @test DatesPlus.Hour("1") == h
    @test DatesPlus.Minute("1") == mi
    @test DatesPlus.Second("1") == s
    @test DatesPlus.Millisecond("1") == ms
    @test DatesPlus.Microsecond("1") == us
    @test DatesPlus.Nanosecond("1") == ns
    @test_throws ArgumentError DatesPlus.Year("1.0")
    @test DatesPlus.Year(parse(Float64, "1.0")) == y

    dt = DatesPlus.DateTime(2014)
    @test typeof(DatesPlus.Year(dt)) <: DatesPlus.Year
    @test typeof(DatesPlus.Quarter(dt)) <: DatesPlus.Quarter
    @test typeof(DatesPlus.Month(dt)) <: DatesPlus.Month
    @test typeof(DatesPlus.Week(dt)) <: DatesPlus.Week
    @test typeof(DatesPlus.Day(dt)) <: DatesPlus.Day
    @test typeof(DatesPlus.Hour(dt)) <: DatesPlus.Hour
    @test typeof(DatesPlus.Minute(dt)) <: DatesPlus.Minute
    @test typeof(DatesPlus.Second(dt)) <: DatesPlus.Second
    @test typeof(DatesPlus.Millisecond(dt)) <: DatesPlus.Millisecond
end
@testset "Default values" begin
    @test DatesPlus.default(DatesPlus.Year) == y
    @test DatesPlus.default(DatesPlus.Quarter) == q
    @test DatesPlus.default(DatesPlus.Month) == m
    @test DatesPlus.default(DatesPlus.Week) == w
    @test DatesPlus.default(DatesPlus.Day) == d
    @test DatesPlus.default(DatesPlus.Hour) == zero(DatesPlus.Hour)
    @test DatesPlus.default(DatesPlus.Minute) == zero(DatesPlus.Minute)
    @test DatesPlus.default(DatesPlus.Second) == zero(DatesPlus.Second)
    @test DatesPlus.default(DatesPlus.Millisecond) == zero(DatesPlus.Millisecond)
    @test DatesPlus.default(DatesPlus.Microsecond) == zero(DatesPlus.Microsecond)
    @test DatesPlus.default(DatesPlus.Nanosecond) == zero(DatesPlus.Nanosecond)
end
@testset "Conversions" begin
    @test DatesPlus.toms(ms) == DatesPlus.value(DatesPlus.Millisecond(ms)) == 1
    @test DatesPlus.toms(s)  == DatesPlus.value(DatesPlus.Millisecond(s)) == 1000
    @test DatesPlus.toms(mi) == DatesPlus.value(DatesPlus.Millisecond(mi)) == 60000
    @test DatesPlus.toms(h)  == DatesPlus.value(DatesPlus.Millisecond(h)) == 3600000
    @test DatesPlus.toms(d)  == DatesPlus.value(DatesPlus.Millisecond(d)) == 86400000
    @test DatesPlus.toms(w)  == DatesPlus.value(DatesPlus.Millisecond(w)) == 604800000

    @test DatesPlus.days(ms) == DatesPlus.days(s) == DatesPlus.days(mi) == DatesPlus.days(h) == 0
    @test DatesPlus.days(DatesPlus.Millisecond(86400000)) == 1
    @test DatesPlus.days(DatesPlus.Second(86400)) == 1
    @test DatesPlus.days(DatesPlus.Minute(1440)) == 1
    @test DatesPlus.days(DatesPlus.Hour(24)) == 1
    @test DatesPlus.days(d) == 1
    @test DatesPlus.days(w) == 7
end
@testset "issue #9214" begin
    @test 2s + (7ms + 1ms) == (2s + 7ms) + 1ms == 1ms + (2s + 7ms) == 1ms + (1s + 7ms) + 1s == 1ms + (2s + 3d + 7ms) + (-3d) == (1ms + (2s + 3d)) + (7ms - 3d) == (1ms + (2s + 3d)) - (3d - 7ms)
    @test 1ms - (2s + 7ms) == -((2s + 7ms) - 1ms) == (-6ms) - 2s
    @test emptyperiod == ((d + y) - y) - d == ((d + y) - d) - y
    @test emptyperiod == 2y + (m - d) + ms - ((m - d) + 2y + ms)
    @test emptyperiod == 0ms
    @test string(emptyperiod) == "empty period"
    @test string(ms + mi + d + m + y + w + h + s + 2y + m) == "3 years, 2 months, 1 week, 1 day, 1 hour, 1 minute, 1 second, 1 millisecond"
    @test 8d - s == 1w + 23h + 59mi + 59s
    @test h + 3mi == 63mi
    @test y - m == 11m
end
@testset "compound periods and types" begin
    # compound periods should avoid automatically converting period types
    @test (d - h).periods == DatesPlus.Period[d, -h]
    @test d - h == 23h
    @test !isequal(d - h, 23h)
    @test isequal(d - h, 2d - 2h - 1d + 1h)
    @test sprint(show, y + m) == string(y + m)
    @test convert(DatesPlus.CompoundPeriod, y) + m == y + m
    @test DatesPlus.periods(convert(DatesPlus.CompoundPeriod, y)) == convert(DatesPlus.CompoundPeriod, y).periods
end
@testset "compound period simplification" begin
    # reduce compound periods into the most basic form
    @test DatesPlus.canonicalize(h - mi).periods == DatesPlus.Period[59mi]
    @test DatesPlus.canonicalize(-h + mi).periods == DatesPlus.Period[-59mi]
    @test DatesPlus.canonicalize(-y + d).periods == DatesPlus.Period[-y, d]
    @test DatesPlus.canonicalize(-y + m - w + d).periods == DatesPlus.Period[-11m, -6d]
    @test DatesPlus.canonicalize(-y + m - w + ms).periods == DatesPlus.Period[-11m, -6d, -23h, -59mi, -59s, -999ms]
    @test DatesPlus.canonicalize(y - m + w - d + h - mi + s - ms).periods == DatesPlus.Period[11m, 6d, 59mi, 999ms]
    @test DatesPlus.canonicalize(-y + m - w + d - h + mi - s + ms).periods == DatesPlus.Period[-11m, -6d, -59mi, -999ms]

    @test DatesPlus.Date(2009, 2, 1) - (DatesPlus.Month(1) + DatesPlus.Day(1)) == DatesPlus.Date(2008, 12, 31)
    @test_throws MethodError (DatesPlus.Month(1) + DatesPlus.Day(1)) - DatesPlus.Date(2009,2,1)
end

@testset "canonicalize Period" begin
    # reduce individual Period into most basic CompoundPeriod
    @test DatesPlus.canonicalize(DatesPlus.Nanosecond(1000000)) == DatesPlus.canonicalize(DatesPlus.Millisecond(1))
    @test DatesPlus.canonicalize(DatesPlus.Millisecond(1000)) == DatesPlus.canonicalize(DatesPlus.Second(1))
    @test DatesPlus.canonicalize(DatesPlus.Second(60)) == DatesPlus.canonicalize(DatesPlus.Minute(1))
    @test DatesPlus.canonicalize(DatesPlus.Minute(60)) == DatesPlus.canonicalize(DatesPlus.Hour(1))
    @test DatesPlus.canonicalize(DatesPlus.Hour(24)) == DatesPlus.canonicalize(DatesPlus.Day(1))
    @test DatesPlus.canonicalize(DatesPlus.Day(7)) == DatesPlus.canonicalize(DatesPlus.Week(1))
    @test DatesPlus.canonicalize(DatesPlus.Month(12)) == DatesPlus.canonicalize(DatesPlus.Year(1))
    @test DatesPlus.canonicalize(DatesPlus.Minute(24*60*1 + 12*60)) == DatesPlus.canonicalize(DatesPlus.CompoundPeriod([DatesPlus.Day(1),DatesPlus.Hour(12)]))
end
@testset "unary ops and vectorized period arithmetic" begin
    pa = [1y 1m 1w 1d; 1h 1mi 1s 1ms]
    cpa = [1y + 1s 1m + 1s 1w + 1s 1d + 1s; 1h + 1s 1mi + 1s 2m + 1s 1s + 1ms]

    @test +pa == pa == -(-pa)
    @test -pa == map(-, pa)
    @test 1y .+ pa == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test (1y + 1m) .+ pa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]
    @test pa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test pa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]

    @test 1y .+ cpa == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test (1y + 1m) .+ cpa == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]
    @test cpa .+ 1y == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test cpa .+ (1y + 1m) == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]

    @test 1y .+ pa == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test (1y + 1m) .+ pa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]
    @test pa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d; 1y + 1h 1y + 1mi 1y + 1s 1y + 1ms]
    @test pa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d; 1y + 1m + 1h 1y + 1m + 1mi 1y + 1m + 1s 1y + 1m + 1ms]

    @test 1y .+ cpa == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test (1y + 1m) .+ cpa == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]
    @test cpa .+ 1y == [2y + 1s 1y + 1m + 1s 1y + 1w + 1s 1y + 1d + 1s; 1y + 1h + 1s 1y + 1mi + 1s 1y + 2m + 1s 1y + 1ms + 1s]
    @test cpa .+ (1y + 1m) == [2y + 1m + 1s 1y + 2m + 1s 1y + 1m + 1w + 1s 1y + 1m + 1d + 1s; 1y + 1m + 1h + 1s 1y + 1m + 1mi + 1s 1y + 3m + 1s 1y + 1m + 1s + 1ms]

    @test 1y .- pa == [0y 1y-1m 1y-1w 1y-1d; 1y-1h 1y-1mi 1y-1s 1y-1ms]
    @test (1y + 1m) .- pa == [1m 1y 1y + 1m-1w 1y + 1m-1d; 1y + 1m-1h 1y + 1m-1mi 1y + 1m-1s 1y + 1m-1ms]
    @test pa .- (1y + 1m) == [-1m -1y -1y-1m + 1w -1y-1m + 1d; -1y-1m + 1h -1y-1m + 1mi -1y-1m + 1s -1y-1m + 1ms]
    @test pa .- 1y == [0y 1m-1y -1y + 1w -1y + 1d; -1y + 1h -1y + 1mi -1y + 1s -1y + 1ms]

    @test 1y .- cpa == [-1s 1y-1m-1s 1y-1w-1s 1y-1d-1s; 1y-1h-1s 1y-1mi-1s 1y-2m-1s 1y-1ms-1s]
    @test (1y + 1m) .- cpa == [1m-1s 1y-1s 1y + 1m-1w-1s 1y + 1m-1d-1s; 1y + 1m-1h-1s 1y + 1m-1mi-1s 1y-1m-1s 1y + 1m-1s-1ms]
    @test cpa .- 1y == [1s -1y + 1m + 1s -1y + 1w + 1s -1y + 1d + 1s; -1y + 1h + 1s -1y + 1mi + 1s -1y + 2m + 1s -1y + 1ms + 1s]
    @test cpa .- (1y + 1m) == [-1m + 1s -1y + 1s -1y-1m + 1w + 1s -1y-1m + 1d + 1s; -1y-1m + 1h + 1s -1y-1m + 1mi + 1s -1y + 1m + 1s -1y + -1m + 1s + 1ms]

    @test [1y 1m; 1w 1d] + [1h 1mi; 1s 1ms] == [1y + 1h 1m + 1mi; 1w + 1s 1d + 1ms]
    @test [1y 1m; 1w 1d] - [1h 1mi; 1s 1ms] == [1y-1h 1m-1mi; 1w-1s 1d-1ms]
    @test [1y 1m; 1w 1d] - [1h 1mi; 1s 1ms] - [1y-1h 1m-1mi; 1w-1s 1d-1ms] == [emptyperiod emptyperiod; emptyperiod emptyperiod]

    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] + [1h 1mi; 1s 1ms] == [1y + 1h + 1s 1m + 1mi + 1s; 1w + 2s 1d + 1s + 1ms]
    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] - [1h 1mi; 1s 1ms] == [1y-1h + 1s 1m-1mi + 1s; 1w 1d + 1s-1ms]

    @test [1y 1m; 1w 1d] + [1h + 1s 1mi + 1s; 1m + 1s 1s + 1ms] == [1y + 1h + 1s 1m + 1mi + 1s; 1w + 1m + 1s 1d + 1s + 1ms]
    @test [1y 1m; 1w 1d] - [1h + 1s 1mi + 1s; 1m + 1s 1s + 1ms] == [1y-1h-1s 1m-1mi-1s; 1w-1m-1s 1d-1s-1ms]

    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] + [1y + 1h 1y + 1mi; 1y + 1s 1y + 1ms] == [2y + 1h + 1s 1y + 1m + 1mi + 1s; 1y + 1w + 2s 1y + 1d + 1s + 1ms]
    @test [1y + 1s 1m + 1s; 1w + 1s 1d + 1s] - [1y + 1h 1y + 1mi; 1y + 1s 1y + 1ms] == [1s-1h 1m + 1s-1y-1mi; 1w-1y 1d + 1s-1y-1ms]
end
@testset "Equality and hashing between FixedPeriod types" begin
    let types = (DatesPlus.Week, DatesPlus.Day, DatesPlus.Hour, DatesPlus.Minute,
                 DatesPlus.Second, DatesPlus.Millisecond, DatesPlus.Microsecond, DatesPlus.Nanosecond)
        for i in 1:length(types), j in i:length(types), x in (0, 1, 235, -4677, 15250)
            local T, U, y, z
            T = types[i]
            U = types[j]
            y = T(x)
            z = convert(U, y)
            @test y == z
            @test hash(y) == hash(z)
        end
    end
end
@testset "Equality and hashing between OtherPeriod types" begin
    for x in (0, 1, 235, -4677, 15250)
        local x, y, z
        y = DatesPlus.Year(x)
        z = convert(DatesPlus.Month, y)
        @test y == z
        @test hash(y) == hash(z)

        y = DatesPlus.Quarter(x)
        z = convert(DatesPlus.Month, y)
        @test y == z
        @test hash(y) == hash(z)

        y = DatesPlus.Year(x)
        z = convert(DatesPlus.Quarter, y)
        @test y == z
        @test hash(y) == hash(z)
    end
end
@testset "Equality and hashing between FixedPeriod/OtherPeriod/CompoundPeriod (#37459)" begin
    function test_hash_equality(x, y)
        @test x == y
        @test y == x
        @test isequal(x, y)
        @test isequal(y, x)
        @test hash(x) == hash(y)
    end
    for FP = (DatesPlus.Week, DatesPlus.Day, DatesPlus.Hour, DatesPlus.Minute,
              DatesPlus.Second, DatesPlus.Millisecond, DatesPlus.Microsecond, DatesPlus.Nanosecond)
        for OP = (DatesPlus.Year, DatesPlus.Quarter, DatesPlus.Month)
            test_hash_equality(FP(0), OP(0))
        end
    end
end
@testset "Hashing for CompoundPeriod (#37447)" begin
    periods = [DatesPlus.Year(0), DatesPlus.Minute(0), DatesPlus.Second(0), DatesPlus.CompoundPeriod(),
               DatesPlus.Minute(2), DatesPlus.Second(120), DatesPlus.CompoundPeriod(DatesPlus.Minute(2)),
               DatesPlus.CompoundPeriod(DatesPlus.Second(120)), DatesPlus.CompoundPeriod(DatesPlus.Minute(1), DatesPlus.Second(60))]
    for x = periods, y = periods
        @test isequal(x,y) == (hash(x) == hash(y))
    end
end

@testset "#30832" begin
    @test DatesPlus.toms(DatesPlus.Second(1) + DatesPlus.Nanosecond(1)) == 1e3
    @test DatesPlus.tons(DatesPlus.Second(1) + DatesPlus.Nanosecond(1)) == 1e9 + 1
    @test DatesPlus.toms(DatesPlus.Second(1) + DatesPlus.Microsecond(1)) == 1e3
end

@testset "CompoundPeriod and Period isless()" begin
    #tests for allowed comparisons
    #FixedPeriod
    @test (h - ms < h + ns) == true
    @test (h + ns < h -ms) == false
    @test (h  < h -ms) == false
    @test (h-ms  < h) == true
    #OtherPeriod
    @test (2y-m < 25m+1y) == true
    @test (2y < 25m+1y) == true
    @test (25m+1y < 2y) == false
    #Test combined Fixed and Other Periods
    @test (1m + 1d < 1m + 1s) == false
end

@testset "Convert CompoundPeriod to Period" begin
    @test convert(Month, Year(1) + Month(1)) === Month(13)
    @test convert(Second, Minute(1) + Second(30)) === Second(90)
    @test convert(Minute, Minute(1) + Second(60)) === Minute(2)
    @test convert(Millisecond, Minute(1) + Second(30)) === Millisecond(90_000)
    @test_throws InexactError convert(Minute, Minute(1) + Second(30))
    @test_throws MethodError convert(Month, Minute(1) + Second(30))
    @test_throws MethodError convert(Second, Month(1) + Second(30))
    @test_throws MethodError convert(Period, Minute(1) + Second(30))
    @test_throws MethodError convert(DatesPlus.FixedPeriod, Minute(1) + Second(30))
end

end


text:
using Random, Test
using ModelObjectsLH, ModelParams, CollegeStratCollege

mdl = CollegeStratCollege;

function h_shocks_test(switches)
    @testset "H shocks $switches" begin
        # nc = 4;
        iCollege = 2;

        rng = MersenneTwister(43);
        hIdxM = UInt8.(rand(rng, 1:5, 4,3));
        hGridNextV = LinRange(1.0, 3.0, 5);
        hStartM = 1.0 .+ 2.0 .* rand(rng, size(hIdxM)...);

        # for switches in switchV
        #     println(switches);

            objId = ObjectId(:hShocks);
            dhSet = make_h_shock_set(objId, switches);
            @test isa(dhSet, HcShockSet)
            @test validate_hshocks(dhSet)
            println(dhSet);

            hShockS = make_h_shock(dhSet, iCollege);
            @test isa(hShockS, HcShock)
            @test validate_hshocks(hShockS)
            # println(hShockS);

            hOutM = sim_h_shocks(hShockS, hStartM, hIdxM, hGridNextV, rng);
            @test isequal(size(hOutM), size(hIdxM))
            @test eltype(hOutM) == eltype(hIdxM)
            @test all(hOutM .>= 1)

            ihEnd = 6;
            pIdxV, prV = pr_hprime(hShockS, ihEnd,
                hStartM[ihEnd-1], hStartM[ihEnd-1] + 0.5, hGridNextV);
            @test sum(prV) ≈ 1.0
            @test all(0.0 .<= prV .<= 1.0)
            @test length(pIdxV) == length(prV)

            if !has_h_shocks(hShockS)
                @test all(hOutM .== hIdxM)
                @test pIdxV == [ihEnd]
            end
            if isa(hShockS, HcShockDown)
                @test all(pIdxV .>= ihEnd - max_steps(dhSet))
                # Take `Int` of indices. Otherwise negative values roll over into large positive ones.
                lbM = max.(1, Int.(hIdxM) .- max_steps(dhSet));
                @test all(lbM .<= hOutM .<= hIdxM)
                if max_steps(dhSet) == 1
                    @test pIdxV == [ihEnd-1, ihEnd]
                    @test prV[2] ≈ 1.0 - shock_prob(hShockS)

                    # Test this for more than 1 step ++++++
                end
            end
        # end
    end
end

function value_before_hshock_test(switches)
    @testset "Value before shocks: $switches" begin
        rng = MersenneTwister(94);
        ic = 2;
        nh = 5;
        vNextV = collect(range(-0.5, 2.2, length = nh));
        dhSet = mdl.make_h_shock_set(ObjectId(:test), switches);
        dh = make_h_shock(dhSet, ic);

        # Compare with simulation from one state
        nSim = 10_000;
        hStart = 2.3;
        ihEnd = 3;
        hGridNextV = LinRange(2.0, 5.0, ihEnd + 2);
        ihPrimeV = sim_h_shocks(dh, fill(hStart, nSim), fill(ihEnd, nSim),
            hGridNextV, rng);

        vNext = sum(vNextV[ihPrimeV]) / nSim;
        # hEnd = hGridNextV[ihEnd];
        vBefore =
            value_before_h_shock(dh, hStart, ihEnd, hGridNextV, vNextV);
        # @test size(vBeforeV) == size(vNextV);
        @test abs(vNext - vBefore) < 1e-2
    end
end

@testset "H shocks" begin
    nc = 4;
    fracLostV = collect(LinRange(0.5, 0.1, nc));
    for switches in (HcShockSwitchesNone(nc),
        make_hshock_switches_down(nc, 1),
        make_hshock_switches_down(nc, 3),
        make_hshock_switches_lose_learn(nc, fracLostV)
        )
        h_shocks_test(switches);
        value_before_hshock_test(switches);
    end
end

# ---------------------


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

n = 10
plot(x=rand(n), y=rand(n), color=sqrt.(rand(n)), Scale.color_sqrt)


text:
function PiecewiseBicubicSplineInterpolator()
    return PiecewiseBicubicSplineInterpolator(())
end

function interpolate(obj::PiecewiseBicubicSplineInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "interpolate", PiecewiseBicubicSplineInterpolatingFunction, (Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end


text:
"""
Table Schema field constraints
https://github.com/frictionlessdata/TableSchema.jl#constraints
"""
mutable struct Constraints
    required::Bool
    unique::Bool
    minLength::Union{Integer, Nothing}
    maxLength::Union{Integer, Nothing}
    # minimum::Integer
    # maximum::Integer
    # pattern
    # enum

    function Constraints(d::Dict)
        # Defaults
        c = new(
            false, # required::Bool
            false, # unique::Bool
            nothing, # minlength::Integer
            nothing, # maxlength::Integer
            # nothing, # minimum::Integer
            # nothing, # maximum::Integer
                # pattern
                # enum
        )
        # Map from dictionary using all the field names of this type
        for fld in fieldnames(typeof(c))
            if haskey(d, String(fld))
                setfield!(c, fld, d[String(fld)])
            end
        end
        return c
    end

    Constraints(required::Bool, unique::Bool) = new(required, unique, nothing, nothing)
    Constraints() = Constraints(Dict())
end

function build(c::Constraints)
    d = Dict()
    for fld in fieldnames(typeof(c))
        d[fld] = getfield(c, fld)
    end
    d
end

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
# Tests light-travel time delay: planets with
# semi-major axis as follows:
a = [0.232456517,1.162282585,2.32456517]
# with time offset between secondary eclipse
# and transit of:
dt = [0.5027,0.51337,0.52643]-0.5
# Convert to minutes:
dt *= 24.*60.
# for an orbital period of a day.

# Define some constants:

C  = 2.99792458e10
AU = 1.49598e13

# Expected time delay (for circular orbit) is 2a/c:
dt_mod = 2*a*AU/C/60.   # In minutes

using PyPlot
clf()
loglog(a,dt_mod,linewidth=3,label=L"$\Delta t = 2a/c$ ",c="r")
scatter(a,dt,s=50.,label=L"$\Delta t$ from starry")
legend(loc="upper left")
xlabel("Semi-major axis [AU]")
ylabel("Time delay [min]")
axis([.2,3,3.,40.])

println(dt - dt_mod)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark