text:
"""
    TranscriptionData

A DataType for storing the data mapping an [`InfiniteOpt.InfiniteModel`](@ref)
that has been transcribed to a regular [`JuMP.Model`](@ref) that contains the
transcribed variables. This is stored in the `ext` field of a `JuMP.Model` to
make what is called a `TranscriptionModel` via the [`TranscriptionModel`](@ref)
constructor.

**Fields**
- `infinite_to_vars::Dict{InfiniteOpt.InfiniteVariableRef,
   Vector{JuMP.VariableRef}}`: Infinite variables to their transcribed variables.
- `hold_to_var::Dict{InfiniteOpt.HoldVariableRef, JuMP.VariableRef}`: Hold
  variables to model variables.
- `point_to_var::Dict{InfiniteOpt.PointVariableRef, JuMP.VariableRef}`: Point
  variables to model variables.
- `infvar_to_supports::Dict{InfiniteOpt.InfiniteVariableRef, Dict}`: Infinite
  variables to transcribed supports indexed by their numeric aliases.
- `infinite_to_constrs::Dict{InfiniteOpt.InfiniteConstraintRef,
  Vector{JuMP.ConstraintRef}}`: Infinite constraints to their transcribed
  constraints.
- `measure_to_constrs::Dict{InfiniteOpt.MeasureConstraintRef,
  Vector{JuMP.ConstraintRef}}`: Measure constraints to model constraints.
- `finite_to_constr::Dict{InfiniteOpt.FiniteConstraintRef, JuMP.ConstraintRef}`:
  Finite constraints to model constraints.
- `infconstr_to_supports::Dict{InfiniteOpt.InfiniteConstraintRef, Dict}`: Infinite
  constraints to the transcribed supports indxed by their numeric aliases.
- `measconstr_to_supports::Dict{InfiniteOpt.MeasureConstraintRef, Dict}`:
  Measure constraints to the transcribed supports indxed by their numeric aliases.
- `infconstr_to_params::Dict{InfiniteOpt.InfiniteConstraintRef, Tuple}`: Infinite
  constraints to the parameter tuples associated with each transcribed support.
- `measconstr_to_params::Dict{InfiniteOpt.MeasureConstraintRef, Tuple}`: Measure
  constraints to the parameter tuples associated with each transcribed support.
"""
mutable struct TranscriptionData
    # Variable mapping
    infinite_to_vars::Dict{InfiniteOpt.InfiniteVariableRef,
                           Vector{JuMP.VariableRef}}
    hold_to_var::Dict{InfiniteOpt.HoldVariableRef, JuMP.VariableRef}
    point_to_var::Dict{InfiniteOpt.PointVariableRef, JuMP.VariableRef}

    # Variable support data
    infvar_to_supports::Dict{InfiniteOpt.InfiniteVariableRef, Vector{<:Tuple}}

    # Constraint mapping
    infinite_to_constrs::Dict{InfiniteOpt.InfiniteConstraintRef,
                              Vector{JuMP.ConstraintRef}}
    measure_to_constrs::Dict{InfiniteOpt.MeasureConstraintRef,
                             Vector{JuMP.ConstraintRef}}
    finite_to_constr::Dict{InfiniteOpt.FiniteConstraintRef, JuMP.ConstraintRef}

    # Constraint support data
    infconstr_to_supports::Dict{InfiniteOpt.InfiniteConstraintRef, Vector{<:Tuple}}
    measconstr_to_supports::Dict{InfiniteOpt.MeasureConstraintRef, Vector{<:Tuple}}
    infconstr_to_params::Dict{InfiniteOpt.InfiniteConstraintRef, Tuple}
    measconstr_to_params::Dict{InfiniteOpt.MeasureConstraintRef, Tuple}

    # Default constructor
    function TranscriptionData()
        return new(Dict{InfiniteOpt.InfiniteVariableRef,
                   Vector{JuMP.VariableRef}}(),
                   Dict{InfiniteOpt.HoldVariableRef, JuMP.VariableRef}(),
                   Dict{InfiniteOpt.PointVariableRef, JuMP.VariableRef}(),
                   Dict{InfiniteOpt.InfiniteVariableRef, Vector{Tuple}}(),
                   Dict{InfiniteOpt.InfiniteConstraintRef,
                        Vector{JuMP.ConstraintRef}}(),
                   Dict{InfiniteOpt.MeasureConstraintRef,
                        Vector{JuMP.ConstraintRef}}(),
                   Dict{InfiniteOpt.FiniteConstraintRef, JuMP.ConstraintRef}(),
                   Dict{InfiniteOpt.InfiniteConstraintRef, Vector{Tuple}}(),
                   Dict{InfiniteOpt.MeasureConstraintRef, Vector{Tuple}}(),
                   Dict{InfiniteOpt.InfiniteConstraintRef, Tuple}(),
                   Dict{InfiniteOpt.MeasureConstraintRef, Tuple}())
    end
end

"""
    TranscriptionModel([optimizer_constructor;
                       caching_mode::MOIU.CachingOptimizerMode = MOIU.AUTOMATIC,
                       bridge_constraints::Bool = true])::JuMP.Model

Return a [`JuMP.Model`](@ref) with [`TranscriptionData`](@ref) included in the
`ext` data field. Accepts the same arguments as a typical JuMP `Model`.

**Example**
```julia-repl
julia> TranscriptionModel()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
function TranscriptionModel(; kwargs...)::JuMP.Model
    model = JuMP.Model(; kwargs...)
    model.ext[:TransData] = TranscriptionData()
    return model
end
# Accept optimizer_factorys
function TranscriptionModel(optimizer_constructor;
                            kwargs...)::JuMP.Model
    model = JuMP.Model(optimizer_constructor; kwargs...)
    model.ext[:TransData] = TranscriptionData()
    return model
end

"""
    is_transcription_model(model::JuMP.Model)::Bool

Return true if `model` is a `TranscriptionModel` or false otherwise.

**Example**
```julia-repl
julia> is_transcription_model(model)
true
```
"""
function is_transcription_model(model::JuMP.Model)::Bool
    return haskey(model.ext, :TransData)
end

"""
    transcription_data(model::JuMP.Model)::TranscriptionData

Return the `TranscriptionData` from a `TranscriptionModel`. Errors if it is not
a `TranscriptionModel`.
"""
function transcription_data(model::JuMP.Model)::TranscriptionData
    !is_transcription_model(model) && error("Model is not a transcription model.")
    return model.ext[:TransData]
end

"""
    transcription_variable(model::JuMP.Model,
                           vref::InfiniteOpt.InfOptVariableRef)

Return the transcribed variable reference(s) corresponding to `vref`. Errors
if no transcription variable is found. Also can query via the syntax:
```julia
transcription_variable(vref::InfiniteOpt.InfOptVariableRef)
```
If the infinite model contains a built transcription model.

**Example**
```julia-repl
julia> transcription_variable(trans_model, infvar)
2-element Array{VariableRef,1}:
 infvar(support: 1)
 infvar(support: 2)

julia> transcription_variable(trans_model, hdvar)
hdvar

julia> transcription_variable(infvar)
2-element Array{VariableRef,1}:
 infvar(support: 1)
 infvar(support: 2)

julia> transcription_variable(hdvar)
hdvar
```
"""
function transcription_variable end

## Define the variable mapping functions
# HoldVariableRef
function transcription_variable(model::JuMP.Model,
                                vref::InfiniteOpt.HoldVariableRef)::JuMP.VariableRef
    !haskey(transcription_data(model).hold_to_var, vref) && error("Variable " *
                             "reference $vref not used in transcription model.")
    return transcription_data(model).hold_to_var[vref]
end
# InfiniteVariableRef
function transcription_variable(model::JuMP.Model,
                                vref::InfiniteOpt.InfiniteVariableRef)::Vector
    !haskey(transcription_data(model).infinite_to_vars, vref) && error("Variable " *
                             "reference $vref not used in transcription model.")
    return transcription_data(model).infinite_to_vars[vref]
end
# PointVariableRef
function transcription_variable(model::JuMP.Model,
                                vref::InfiniteOpt.PointVariableRef)::JuMP.VariableRef
    !haskey(transcription_data(model).point_to_var, vref) && error("Variable " *
                             "reference $vref not used in transcription model.")
    return transcription_data(model).point_to_var[vref]
end
# Dispatch for internal models
function transcription_variable(vref::InfiniteOpt.InfOptVariableRef)
    trans_model = InfiniteOpt.optimizer_model(JuMP.owner_model(vref))
    return transcription_variable(trans_model, vref)
end

"""
    InfiniteOpt.optimizer_model_variable(vref::InfiniteOpt.InfOptVariableRef,
                                         ::Val{:TransData})

Proper extension of [`InfiniteOpt.optimizer_model_variable`](@ref) for
`TranscriptionModel`s. This simply dispatches to [`transcription_variable`](@ref).
"""
function InfiniteOpt.optimizer_model_variable(vref::InfiniteOpt.InfOptVariableRef,
                                              ::Val{:TransData})
    return transcription_variable(vref)
end

"""
    InfiniteOpt.variable_supports(model::JuMP.Model,
                                  vref::InfiniteOpt.InfiniteVariableRef,
                                  key::Val{:TransData} = Val(:TransData))::Vector

Return the support alias mapping associated with `vref` in the transcribed model.
Errors if `vref` does not have transcribed variables.
"""
function InfiniteOpt.variable_supports(model::JuMP.Model,
                                       vref::InfiniteOpt.InfiniteVariableRef,
                                       key::Val{:TransData} = Val(:TransData))::Vector
    if !haskey(transcription_data(model).infvar_to_supports, vref)
        error("Variable reference $vref not used in transcription model.")
    end
    return transcription_data(model).infvar_to_supports[vref]
end

"""
    transcription_constraint(model::JuMP.Model,
                             cref::InfiniteOpt.GeneralConstraintRef)

Return the transcribed constraint reference(s) corresponding to `cref`. Errors
if `cref` has not been transcribed. Also can query via the syntax:
```julia
transcription_constraint(cref::InfiniteOpt.GeneralConstraintRef)
```
If the infinite model contains a built transcription model.

**Example**
```julia-repl
julia> transcription_constraint(trans_model, fin_con)
fin_con : x(support: 1) - y <= 3.0

julia> transcription_constraint(fin_con)
fin_con : x(support: 1) - y <= 3.0
```
"""
function transcription_constraint end

## Define the cosntraint mapping functions
# InfiniteConstraintRef
function transcription_constraint(model::JuMP.Model,
                                  cref::InfiniteOpt.InfiniteConstraintRef)::Vector
    if !haskey(transcription_data(model).infinite_to_constrs, cref)
        error("Constraint reference $cref not used in transcription model.")
    end
    return transcription_data(model).infinite_to_constrs[cref]
end

# MeasureConstraintRef
function transcription_constraint(model::JuMP.Model,
                                  cref::InfiniteOpt.MeasureConstraintRef)::Vector
    if !haskey(transcription_data(model).measure_to_constrs, cref)
        error("Constraint reference $cref not used in transcription model.")
    end
    return transcription_data(model).measure_to_constrs[cref]
end

# FiniteConstraintRef
function transcription_constraint(model::JuMP.Model,
                                  cref::InfiniteOpt.FiniteConstraintRef)::JuMP.ConstraintRef
    if !haskey(transcription_data(model).finite_to_constr, cref)
        error("Constraint reference $cref not used in transcription model.")
    end
    return transcription_data(model).finite_to_constr[cref]
end
# Dispatch for internal models
function transcription_constraint(cref::InfiniteOpt.GeneralConstraintRef)
    trans_model = InfiniteOpt.optimizer_model(JuMP.owner_model(cref))
    return transcription_constraint(trans_model, cref)
end

"""
    InfiniteOpt.optimizer_model_constraint(cref::InfiniteOpt.GeneralConstraintRef,
                                           ::Val{:TransData})

Proper extension of [`InfiniteOpt.optimizer_model_constraint`](@ref) for
`TranscriptionModel`s. This simply dispatches to [`transcription_constraint`](@ref).
"""
function InfiniteOpt.optimizer_model_constraint(cref::InfiniteOpt.GeneralConstraintRef,
                                                ::Val{:TransData})
    return transcription_constraint(cref)
end

"""
    InfiniteOpt.constraint_supports(model::JuMP.Model,
                                    cref::InfiniteOpt.GeneralConstraintRef,
                                    key::Val{:TransData} = Val(:TransData))::Vector

Return the support alias mappings associated with `cref`. Errors if `cref` is
not transcribed.
"""
function InfiniteOpt.constraint_supports(model::JuMP.Model,
                                         cref::InfiniteOpt.InfiniteConstraintRef,
                                         key::Val{:TransData} = Val(:TransData))::Vector
    if !haskey(transcription_data(model).infconstr_to_supports, cref)
        error("Constraint reference $cref not used in transcription model.")
    end
    return transcription_data(model).infconstr_to_supports[cref]
end
# MeasureConstraintRef
function InfiniteOpt.constraint_supports(model::JuMP.Model,
                                         cref::InfiniteOpt.MeasureConstraintRef,
                                         key::Val{:TransData} = Val(:TransData))::Vector
    if !haskey(transcription_data(model).measconstr_to_supports, cref)
        error("Constraint reference $cref not used in transcription model " *
              "and/or is finite and doesn't have supports.")
    end
    return transcription_data(model).measconstr_to_supports[cref]
end

"""
    InfiniteOpt.constraint_parameter_refs(model::JuMP.Model,
                                          cref::InfiniteOpt.GeneralConstraintRef,
                                          key::Val{:TransData} = Val(:TransData))::Tuple

Return the a parameter reference tuple of all the parameters that parameterize
`cref` and correspond to the supports. Errors if `cref` has not been transcribed.
"""
function InfiniteOpt.constraint_parameter_refs(model::JuMP.Model,
                                               cref::InfiniteOpt.InfiniteConstraintRef,
                                               key::Val{:TransData} = Val(:TransData))::Tuple
    if !haskey(transcription_data(model).infconstr_to_params, cref)
        error("Constraint reference $cref not used in transcription model " *
              "and/or is finite and doesn't have parameter references.")
    end
    return transcription_data(model).infconstr_to_params[cref]
end
# MeasureConstraintRef
function InfiniteOpt.constraint_parameter_refs(model::JuMP.Model,
                                               cref::InfiniteOpt.MeasureConstraintRef,
                                               key::Val{:TransData} = Val(:TransData))::Tuple
    if !haskey(transcription_data(model).measconstr_to_params, cref)
        error("Constraint reference $cref not used in transcription model " *
              "and/or is finite and doesn't have parameter references.")
    end
    return transcription_data(model).measconstr_to_params[cref]
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The workforce scheduling problem
#
# This model determines a set of workforce levels that will most economically
# meet demands and inventory requirements over time. The formulation is
# motivated by the experiences of a large producer in the United States. The
# data are for three products and 13 periods.
#
# Problem taken from the Appendix C of the expanded version of Fourer, Gay, and
# Kernighan, A Modeling Language for Mathematical Programming
#
# Originally contributed by Louis Luangkesorn, February 26, 2015.

using JuMP
import GLPK
import Test

function example_prod(; verbose = true)
    ## PRODUCTION SETS AND PARAMETERS
    prd = ["18REG"  "24REG" "24PRO"]
    ## Members of the product group
    numprd = length(prd)
    pt =	[1.194,	1.509,	1.509]
    ## Crew-hours to produce 1000 units
    pc =	[2304,	2920,	2910]
    ## Nominal production cost per 1000, used
    ## to compute inventory and shortage costs
    ##
    ## TIME PERIOD SETS AND PARAMETERS
    firstperiod = 1
    ## Index of first production period to be modeled
    lastperiod  = 13
    ## Index of last production period to be modeled
    numperiods = firstperiod:lastperiod
    ## 'planning horizon' := first..last;
    ## EMPLOYMENT PARAMETERS
    ## Workers per crew
    cs = 18
    ## Regular-time hours per shift
    sl =  8
    ## Wage per hour for regular-time labor
    rtr = 16.00
    ## Wage per hour for overtime labor
    otr = 43.85
    ## Crews employed at start of first period
    iw =  8
    ## Regular working days in a production period
    dpp =	 [19.5,	19,	20,	19,	19.5,	19,	19,	20,	19,	20,	20,	18,	18]
    ## Maximum crew-hours of overtime in a period
    ol =	 [96,	96,	96,	96,	96,	96,	96,	96,	96,	96,	96,	96,	96]
    ## Lower limit on average employment in a period
    cmin =	[0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0]
    ## Upper limit on average employment in a period
    cmax =	[8,	8,	8,	8,	8,	8,	8,	8,	8,	8,	8,	8,	8]
    ## Penalty cost of hiring a crew
    hc =	 [7500,	7500,	7500,	7500,	15000,	15000,	15000,	15000,	15000,	15000,	7500,	7500,	7500]
    ## Penalty cost of laying off a crew
    lc =	 [7500,	7500,	7500,	7500,	15000,	15000,	15000,	15000,	15000,	15000,	7500,	7500,	7500]
    ## DEMAND PARAMETERS
    d18REG = [63.8,	76,	88.4,	913.8,	115,	133.8,	79.6,	111,	121.6,	470,	78.4,	99.4,	140.4,	63.8]
    d24REG = [1212,	306.2,	319,	208.4,	298,	328.2,	959.6,	257.6,	335.6,	118,	284.8,	970,	343.8,	1212]
    d24PRO = [0,	0,	0,	0,	0,	0,	0,	0,	0,	1102,	0,	0,	0,	0]
    ## Requirements (in 1000s) to be met from current production and inventory
    dem = Array[d18REG, d24REG, d24PRO]
    ## true if product will be the subject of a special promotion in the period
    pro = Array[
        [0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0],
        [1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	1],
        [0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0],
    ]
    ## INVENTORY AND SHORTAGE PARAMETERS
    ## Proportion of non-promoted demand that must be in inventory the previous
    ## period
    rir = 0.75
    ## Proportion of promoted demand that must be in inventory the previous
    ## period
    pir = 0.80
    ## Upper limit on number of periods that any product may sit in inventory
    life = 2
    ## Inventory cost per 1000 units is cri times nominal production cost
    cri	= [0.015,	0.015,	0.015]
    ## Shortage cost per 1000 units is crs times nominal production cost
    crs	= [1.1,	1.1,	1.1]
    ## Inventory at start of first period; age unknown
    iinv = [82,	792.2,	0]
    ## Initial inventory still available for allocation at end of period t
    iil = [
        [
            max(0, iinv[p] - sum(dem[p][v] for v in firstperiod:t))
            for t in numperiods
        ]
        for p in 1:numprd
    ]
    ## Lower limit on inventory at end of period t
    function checkpro(
        product, timeperiod, production, promotionalrate, regularrate
    )
        if production[product][timeperiod + 1] == 1
            return promotionalrate
        else
            return regularrate
        end
    end
    minv = [
        [dem[p][t + 1] * checkpro(p, t, pro, pir, rir) for t in numperiods]
        for p in 1:numprd
    ]
    ## DEFINE MODEL
    prod = Model(GLPK.Optimizer)
    ## VARIABLES
    ## Average number of crews employed in each period
    @variable(prod, Crews[0:lastperiod] >= 0)
    ## Crews hired from previous to current period
    @variable(prod, Hire[numperiods] >= 0)
    ## Crews laid off from previous to current period
    @variable(prod, Layoff[numperiods]>= 0)
    ## Production using regular-time labor, in 1000s
    @variable(prod, Rprd[1:numprd, numperiods] >= 0)
    ## Production using overtime labor, in 1000s
    @variable(prod, Oprd[1:numprd, numperiods]>= 0)
    ## a numperiods old -- produced in period (t+1)-a --
    ## and still in storage at the end of period t
    @variable(prod, Inv[1:numprd, numperiods, 1:life] >= 0)
    ## Accumulated unsatisfied demand at the end of period t
    @variable(prod, Short[1:numprd, numperiods] >= 0)
    ## CONSTRAINTS
    ## Hours needed to accomplish all regular-time production in a period must
    ## not exceed hours available on all shifts
    @constraint(
        prod,
        [t = numperiods],
        sum(pt[p] * Rprd[p, t] for p in 1:numprd) <= sl * dpp[t] * Crews[t]
    )
    ## Hours needed to accomplish all overtime production in a period must not
    ## exceed the specified overtime limit
    @constraint(
        prod,
        [t = numperiods],
        sum(pt[p] * Oprd[p, t] for p in 1:numprd)  <= ol[t]
    )
    ## Use given initial workforce
    @constraint(prod, Crews[firstperiod - 1] == iw)
    ## Workforce changes by hiring or layoffs
    @constraint(
        prod, [t in numperiods], Crews[t] == Crews[t - 1] + Hire[t] - Layoff[t]
    )
    ## Workforce must remain within specified bounds
    @constraint(prod, [t in numperiods], cmin[t] <= Crews[t])
    @constraint(prod, [t in numperiods], Crews[t] <= cmax[t])
    ## 'first demand requirement
    @constraint(
        prod,
        [p in 1:numprd],
        Rprd[p, firstperiod] + Oprd[p, firstperiod] + Short[p, firstperiod] -
            Inv[p, firstperiod, 1] == max(0, dem[p][firstperiod] - iinv[p])
    )
    ## Production plus increase in shortage plus decrease in inventory must
    ## equal demand
    for t in (firstperiod + 1):lastperiod
        @constraint(
            prod,
            [p in 1:numprd],
            Rprd[p, t] + Oprd[p, t] + Short[p,t] - Short[p,t-1] +
                sum(Inv[p, t - 1, a] - Inv[p, t, a] for a in 1:life) ==
                max(0, dem[p][t] - iil[p][t - 1])
        )
    end
    ## Inventory in storage at end of period t must meet specified minimum
    @constraint(
        prod,
        [p in 1:numprd, t in numperiods],
        sum(Inv[p, t, a] + iil[p][t] for a in 1:life) >= minv[p][t]
    )
    ## In the vth period (starting from first) no inventory may be more than v
    ## numperiods old (initial inventories are handled separately)
    @constraint(
        prod,
        [p in 1:numprd, v in 1:(life - 1), a in (v + 1):life],
        Inv[p, firstperiod + v - 1, a] == 0
    )
    ## New inventory cannot exceed production in the most recent period
    @constraint(
        prod,
        [p in 1:numprd, t in numperiods],
        Inv[p, t, 1] <= Rprd[p, t] + Oprd[p, t]
    )
    ## Inventory left from period (t+1)-p can only decrease as time goes on
    secondperiod = firstperiod + 1
    @constraint(
        prod,
        [p in 1:numprd, t in 2:lastperiod, a in 2:life],
        Inv[p, t, a] <= Inv[p, t - 1, a - 1]
    )
    ## OBJECTIVE
    ## Full regular wages for all crews employed, plus penalties for hiring and
    ## layoffs, plus wages for any overtime worked, plus inventory and shortage
    ## costs. (All other production costs are assumed to depend on initial
    ## inventory and on demands, and so are not included explicitly.)
    @objective(
        prod,
        Min,
        sum(
            rtr * sl * dpp[t] * cs * Crews[t] +
            hc[t] * Hire[t] +
            lc[t] * Layoff[t] +
            sum(
                otr * cs * pt[p] * Oprd[p, t] +
                sum(cri[p] * pc[p] * Inv[p, t, a] for a in 1:life) +
                crs[p] * pc[p] * Short[p, t]
                for p in 1:numprd
            )
           for t in numperiods
        )
    )
    ## Obtain solution
    optimize!(prod)
    Test.@test termination_status(prod) == MOI.OPTIMAL
    Test.@test primal_status(prod) == MOI.FEASIBLE_POINT
    Test.@test objective_value(prod) ≈ 4_426_822.89 atol = 1e-2
    if verbose
        println("RESULTS:")
        println("Crews")
        for t = 0:length(Crews.data) - 1
            print(" $(value(Crews[t])) ")
        end
        println()
        println("Hire")
        for t = 1:length(Hire.data)
            print(" $(value(Hire[t])) ")
        end
        println()
        println("Layoff")
        for t = 1:length(Layoff.data)
            print(" $(value(Layoff[t])) ")
        end
        println()
    end
    return
end

example_prod()

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function showextremes()
  values = [0.0, -0.0, Inf, -Inf, NaN]
  println(1 ./ values)
end

showextremes()

@show Inf + 2.0
@show Inf + Inf
@show Inf - Inf
@show Inf * Inf
@show Inf / Inf
@show Inf * 0
@show 0 == -0
@show NaN == NaN
@show NaN === NaN


text:
# Simple Text File

immutable Simple <: FileFormat end

const _simpleparser_start  = 2
const _simpleparser_first_final  = 2
const _simpleparser_error  = 0
const _simpleparser_en_main  = 2
const __simpleparser_nfa_targs = Int8[ 0, 0 ,  ]
const __simpleparser_nfa_offsets = Int8[ 0, 0, 0, 0 ,  ]
const __simpleparser_nfa_push_actions = Int8[ 0, 0 ,  ]
const __simpleparser_nfa_pop_trans = Int8[ 0, 0 ,  ]
type SimpleParser <: AbstractParser
	state::Ragel.State
	seqbuf::BufferedOutputStream{	BufferedStreams.EmptyStreamSource}

function SimpleParser(input::BufferedInputStream)
	begin
	cs = convert( Int , _simpleparser_start );

end
return new(Ragel.State(cs, input), BufferedOutputStream())
end
	end

Base.eltype(::Type{SimpleParser}) = Line

type Line
data::ASCIIString
Line() = new()
Line(line) = new(line)
end

Base.convert(::Type{ASCIIString}, line::Line) = line.data

function Base.open(input::BufferedInputStream, ::Type{Simple})
return SimpleParser(input)
end

Base.copy(line::Line) = Line(copy(line.data))

Ragel.@generate_read_fuction(
"_simpleparser",
SimpleParser,
Line,
begin
begin
if ( p == pe  )
@goto _test_eof

end
if ( cs  == 2 )
@goto st_case_2
elseif ( cs  == 1 )
@goto st_case_1
elseif ( cs  == 0 )
@goto st_case_0
end
@goto st_out
@label ctr1
begin
output.data = Ragel.@ascii_from_anchor!
yield = true
begin
p+= 1;
cs = 2;
@goto _out

end

end
begin
state.linenum += 1
end
@goto st2
@label st2
p+= 1;
if ( p == pe  )
@goto _test_eof2

end
@label st_case_2
if ( (data[1+(p )])== 10  )
begin
@goto st0

end

end
begin
@goto ctr3

end
@label ctr2
begin
output.data = Ragel.@ascii_from_anchor!
yield = true
begin
p+= 1;
cs = 1;
@goto _out

end

end
@goto st1
@label ctr3
begin
Ragel.@anchor!
end
@goto st1
@label st1
p+= 1;
if ( p == pe  )
@goto _test_eof1

end
@label st_case_1
if ( (data[1+(p )]) == 10 )
begin
@goto ctr1

end
elseif ( (data[1+(p )]) == 13 )
begin
@goto ctr2

end
end
begin
@goto st1

end
@label st_case_0
@label st0
cs = 0;
@goto _out
@label st_out
@label _test_eof2
cs = 2;
@goto _test_eof
@label _test_eof1
cs = 1;
@goto _test_eof
@label _test_eof
begin

end
@label _out
begin

end

end
end
)


text:
# ---
# title: 971. Flip Binary Tree To Match Preorder Traversal
# id: problem971
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/flip-binary-tree-to-match-preorder-traversal/description/>
# hidden: true
# ---
#
# Given a binary tree with `N` nodes, each node has a different value from `{1,
# ..., N}`.
#
# A node in this binary tree can be _flipped_  by swapping the left child and
# the right child of that node.
#
# Consider the sequence of `N` values reported by a preorder traversal starting
# from the root.  Call such a sequence of `N` values the  _voyage_  of the tree.
#
# (Recall that a _preorder traversal_  of a node means we report the current
# node's value, then preorder-traverse the left child, then preorder-traverse
# the right child.)
#
# Our goal is to flip the **least number** of nodes in the tree so that the
# voyage of the tree matches the `voyage` we are given.
#
# If we can do so, then return a list of the values of all nodes flipped.  You
# may return the answer in any order.
#
# If we cannot do so, then return the list `[-1]`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/02/1219-01.png)**
#
#
#
#     Input: root = [1,2], voyage = [2,1]
#     Output: [-1]
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/02/1219-02.png)**
#
#
#
#     Input: root = [1,2,3], voyage = [1,3,2]
#     Output: [1]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/02/1219-02.png)**
#
#
#
#     Input: root = [1,2,3], voyage = [1,2,3]
#     Output: []
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Weber
using Base.Test
include("find_timing.jl")

seq_trial_events,_,rows = find_timing() do
  addtrial(moment(() -> record(:a)),
           moment(() -> record(:b)),
           moment(() -> record(:c)))
  addtrial(moment(() -> record(:a)),
           moment(() -> record(:b)),
           moment(() -> record(:c)))
  addtrial(moment(() -> record(:a)),
           moment(() -> record(:b)),
           moment(() -> record(:c)))
end
seq_trial_index = map(x -> x[:trial],rows)

@testset "Moment Indexing" begin
  @test seq_trial_events == [:a,:b,:c,:a,:b,:c,:a,:b,:c]
  @test seq_trial_index == [1,1,1,2,2,2,3,3,3]
end


text:
module ClpMathProgSolverInterface
using Clp.ClpCInterface
using Compat
using Compat.LinearAlgebra
using Compat.SparseArrays

import MathProgBase
const MPB = MathProgBase

export ClpMathProgModel, ClpSolver

mutable struct ClpMathProgModel <: MPB.AbstractLinearQuadraticModel
    inner::ClpModel
    solveroptions::ClpSolve
end

struct ClpSolver <: MPB.AbstractMathProgSolver
    options
end
ClpSolver(;kwargs...) = ClpSolver(kwargs)

### Options

# map option name to C function
const optionmap = Dict(
   :PrimalTolerance => set_primal_tolerance,
   :DualTolerance => set_dual_tolerance,
   :DualObjectiveLimit => set_dual_objective_limit,
   :MaximumIterations => set_maximum_iterations,
   :MaximumSeconds => set_maximum_seconds,
   :LogLevel => set_log_level,
   :Scaling => scaling,
   :Perturbation => set_perturbation,
   #:Algorithm => set_algorithm
   )
# These options are set by using the ClpSolve object
const solveoptionmap = Dict(
   :PresolveType => set_presolve_type,
   :SolveType => set_solve_type,
   :InfeasibleReturn => set_infeasible_return,
   )

function setoption(m::ClpMathProgModel, name::Symbol, value)
    if haskey(optionmap, name)
        optionmap[name](m.inner,value)
    elseif haskey(solveoptionmap, name)
        solveoptionmap[name](m.solveroptions,value)
    else
        error("Unrecognized option: $name")
    end
end

function ClpMathProgModel(;kwargs...)
    m = ClpMathProgModel(ClpModel(),ClpSolve())
    set_log_level(m.inner,0) # disable output by default
    for (name, value) in kwargs
        setoption(m, name, value)
    end
    return m
end

MPB.LinearQuadraticModel(s::ClpSolver) = ClpMathProgModel(;s.options...)

MPB.ConicModel(s::ClpSolver) = MPB.LPQPtoConicBridge(MPB.LinearQuadraticModel(s))
MPB.supportedcones(s::ClpSolver) = [:Free,:Zero,:NonNeg,:NonPos]

function MPB.loadproblem!(m::ClpMathProgModel, filename::AbstractString)
    if endswith(filename,".mps") || endswith(filename,".mps.gz")
       read_mps(m.inner,filename)
    else
       error("unrecognized input format extension in $filename")
    end
end

function MPB.loadproblem!(m::ClpMathProgModel, A, collb, colub, obj, rowlb,
                          rowub, sense)
    load_problem(m.inner,A,collb,colub,obj,rowlb,rowub)
    MPB.setsense!(m, sense)
end

#writeproblem(m, filename::String)

function replaceInf(x)
    for i in 1:length(x)
        if x[i] > 1e20
            x[i] = Inf
        elseif x[i] < -1e20
            x[i] = -Inf
        end
    end
    return x
end

MPB.getvarLB(m::ClpMathProgModel) = replaceInf(get_col_lower(m.inner))
MPB.setvarLB!(m::ClpMathProgModel, collb) = chg_column_lower(m.inner, collb)

MPB.getvarUB(m::ClpMathProgModel) = replaceInf(get_col_upper(m.inner))
MPB.setvarUB!(m::ClpMathProgModel, colub) = chg_column_upper(m.inner, colub)

MPB.getconstrLB(m::ClpMathProgModel) = replaceInf(get_row_lower(m.inner))
MPB.setconstrLB!(m::ClpMathProgModel, rowlb) = chg_row_lower(m.inner, rowlb)

MPB.getconstrUB(m::ClpMathProgModel) = replaceInf(get_row_upper(m.inner))
MPB.setconstrUB!(m::ClpMathProgModel, rowub) = chg_row_upper(m.inner, rowub)

MPB.getobj(m::ClpMathProgModel) = get_obj_coefficients(m.inner)
MPB.setobj!(m::ClpMathProgModel, obj) = chg_obj_coefficients(m.inner, obj)

function MPB.getconstrmatrix(m::ClpMathProgModel)
    A = get_constraint_matrix(m.inner)
    return convert(SparseMatrixCSC{Float64,Int},A)
end

function MPB.addvar!(m::ClpMathProgModel, rowidx, rowcoef, collb, colub, objcoef)
    @assert length(rowidx) == length(rowcoef)
    colstarts = Int32[0, length(rowcoef)]
    rows = Int32[ i - 1 for i in rowidx ]
    add_columns(m.inner, 1, Float64[collb], Float64[colub], Float64[objcoef],
       colstarts, rows, convert(Vector{Float64},rowcoef))
end

function MPB.addconstr!(m::ClpMathProgModel, colidx, colcoef, rowlb, rowub)
    @assert length(colidx) == length(colcoef)
    rowstarts = Int32[0, length(colcoef)]
    cols = Int32[ i - 1 for i in colidx ]
    add_rows(m.inner, 1, Float64[rowlb], Float64[rowub], rowstarts, cols, convert(Vector{Float64}, colcoef))
end

function MPB.setsense!(m::ClpMathProgModel,sense)
    if sense == :Min
        set_obj_sense(m.inner, 1.0)
    elseif sense == :Max
        set_obj_sense(m.inner, -1.0)
    else
        error("Unrecognized objective sense $sense")
    end
end

function MPB.getsense(m::ClpMathProgModel)
    s = get_obj_sense(m.inner)
    if s == 1.0
        return :Min
    elseif s == -1.0
        return :Max
    else
        error("Internal library error")
    end
end

MPB.numvar(m::ClpMathProgModel) = get_num_cols(m.inner)
MPB.numconstr(m::ClpMathProgModel) = get_num_rows(m.inner)

MPB.optimize!(m::ClpMathProgModel) = initial_solve_with_options(m.inner,m.solveroptions)

function MPB.status(m::ClpMathProgModel)
   s = ClpCInterface.status(m.inner)
   if s == 0
       return :Optimal
   elseif s == 1
       return :Infeasible
   elseif s == 2
       return :Unbounded
   elseif s == 3
       return :UserLimit
   elseif s == 4
       return :Error
   else
       error("Internal library error")
   end
end

MPB.getobjval(m::ClpMathProgModel) = objective_value(m.inner)

MPB.getsolution(m::ClpMathProgModel) = primal_column_solution(m.inner)
MPB.getconstrsolution(m::ClpMathProgModel) = primal_row_solution(m.inner)
MPB.getreducedcosts(m::ClpMathProgModel) = dual_column_solution(m.inner)

MPB.getconstrduals(m::ClpMathProgModel) = dual_row_solution(m.inner)

function MPB.getinfeasibilityray(m::ClpMathProgModel)
    return Compat.rmul!(infeasibility_ray(m.inner),-1.0)
end
MPB.getunboundedray(m::ClpMathProgModel) = unbounded_ray(m.inner)

const statmap = Dict(zip([ 0x00,  0x01,            0x02,            0x03,       0x04,  0x05],
                     [:Free,:Basic,:NonbasicAtUpper,:NonbasicAtLower,:Superbasic,:Fixed]))
function MPB.getbasis(m::ClpMathProgModel)
    num_cols = MPB.numvar(m)
    num_rows = MPB.numconstr(m)
    cbasis = Array{Symbol}(undef, num_cols)
    rbasis = Array{Symbol}(undef, num_rows)
    for i in 1:num_cols
        val = get_column_status(m.inner, i)
        cbasis[i] = statmap[val]
    end
    for i in 1:num_rows
        val = get_row_status(m.inner, i)
        rbasis[i] = statmap[val]
    end
    return cbasis, rbasis
end

MPB.getvartype(m::ClpMathProgModel) = fill(:Cont, get_num_cols(m.inner))
function MPB.setvartype!(m::ClpMathProgModel, typ::Vector{Symbol})
    all(x->isequal(x,:Cont), typ) || error("Clp does not support integer variables")
    return nothing
end

MPB.getrawsolver(m::ClpMathProgModel) = m.inner

end


text:
# INVESTIGATE THE FOLLOWING:
# * Orthogonalize the conjugate vector and the residual vector together, periodically against previous vectors, does it
#   reduce the number of iterations?

mutable struct ConjGradModified{T <: Number}
    A::Function              # Linear opeartor
    b::AbstractArray{T}      # RHS vector
    tensor_size::Tuple       # The size of the tensor the linear operator is acting on.

    x::Vector{T}        # current guess, started with initial guess.
    r::Vector{T}        # previous computed residual
    rnew::Vector{T}     # Current residual
    d::Vector{T}        # Conjugate Direction
    Ad::Vector{T}       # for reducing garbage collector time.

    itr::UInt64              # Iteration count.
    Q::Vector{Vector{T}}     # list of normalized residuals vectors.
    Ds::Vector{Vector{T}}    # ist of normalized Ad vector during the iterations.

    function ConjGradModified(
        A::Function,
        b::AbstractArray,
        x0::Union{AbstractArray,Nothing}=nothing
    )
        x0 = x0===nothing ? b .+ 0.1 : x0
        r = reshape(b - A(x0), :)
        T = typeof(r[1])
        this = new{T}()
        this.r = r
        this.A = A
        this.b = b
        this.tensor_size = size(b)
        this.x = x0
        this.rnew = similar(this.r)
        this.d = this.r
        this.Ad = similar(this.r)
        ComputeVec!(this, this.r, this.Ad)
        this.x = x0
        this.itr = 0

        this.Q = Vector{Vector{T}}()
        this.Ds = Vector{Vector{T}}()
        push!(this.Q, this.r/norm(this.r))
        push!(this.Ds, this.d/norm(this.d))
        return this
    end

    function ConjGradModified(
            A::AbstractMatrix,
            b::AbstractVector,
            x0::Union{AbstractArray,Nothing}=nothing
        )
        return ConjGradModified((x)-> A*x, b, x0)
    end

end

"""
    Apply the linear transformation on the tensor and return the
    residual vector after transformation.
"""
function ComputeResidualVec(this::ConjGradModified, x::AbstractArray)
    resVec = this.b - this.A(reshape(x, this.tensor_size))
    return reshape(resVec, length(resVec))
end

"""
    Get the norm of the current residual of Modified Conjugate
    Gradient.
"""
function GetResidualNorm(this::ConjGradModified)
    return norm(this.r)
end

"""
    Compute the vector using the linear transformation.
"""
function ComputeVec(this::ConjGradModified, x::AbstractArray)
    vec = this.A(reshape(x, this.tensor_size))
    return reshape(vec, length(vec))
end

"""
    Pre-allocates a vector and get the values to the given
    pre-allocated vector.
"""
function ComputeVec!(this::ConjGradModified, x::AbstractArray, vec::AbstractArray)
    vec .= reshape(this.A(reshape(x, this.tensor_size)), :)
    return nothing
end

"""
    Force the orthogonalization happen between the residual vectors against
    all previous residual vectors, and the conjugate vector as well.
"""
function ForceOrthogonalize(this::ConjGradModified)
    for Idx in 1: length(this.Q)
        for Jdx in 1: Idx - 1
            # this.Q[Idx] .-= dot(this.Q[Jdx], this.Q[Idx])*this.Q[Jdx]
            this.Ds[Idx] .-= dot(this.Ds[Jdx], this.Ds[Idx])*this.Ds[Jdx]
        end
    end
    for  Aq in this.Q
        this.d .-= dot(Aq, this.d)*Aq
    end

    return
end

"""
    Performs One step of conjugate Gradient.
"""
function (this::ConjGradModified)()
    r = this.r
    if norm(r) == 0
        return 0 # The problem is solved already.
    end
    d = this.d
    Ad = this.Ad
    a = dot(r, r)/dot(d, Ad)

    if a < 0
        error("CG got a non-definite matrix")
    end

    this.x += a*d
    this.rnew .= r - a*Ad                    # update rnew
    rnewNorm = norm(this.rnew)

    b = dot(this.rnew, this.rnew)/dot(r, r)
    # @assert abs(dot(rnew + β*d, Ad)) < 1e-8 "Not conjugate"
    this.d = this.rnew + b*d
    this.r = copy(this.rnew)                      # Override
    this.itr += 1
    ComputeVec!(this, this.d, Ad)
    return convert(Float64, rnewNorm)
end

"""
    Performs j iterations of conjugate gradient, collect the
    2norm of the residual and return all of them in a vector.
"""
function (this::ConjGradModified)(j::Int64)
    results = Vector()
    for _ in 1: j
        push!(results, this())
    end
    return results
end

### Using stuff.
using Test, LinearAlgebra
@testset begin

    function Test1()
        n = 10
        A = rand(n, n)
        A = A'*A
        display(A)
        b = rand(n)
        cg = ConjGradModified(A, b)
        for I in 1:n + 1
            println("Iter: $I; "*"$(cg())")
        end
        return norm(b - A*cg.x) < 1e-8

    end
    function Test2()
        n = 5
        A = rand(n, n)
        A = A'*A
        display(A)
        b = rand(n)
        cg = ConjGradModified(A, b)
        for I in 1:n + 1
            println("Iter: $I; "*"$(cg())")
            ForceOrthogonalize(cg)
        end
        return norm(b - A*cg.x) < 1e-8
    end

    @test Test1()
    @test Test2()
end


text:
module API

import ..Counter: Domain

using .Domain

export create

"""
    create(name::String, nbr::Int64)::Item

Create an item object.

# Example
```julia
julia> using counter

julia> ar = create("ar", 1000);

julia> ar.created
2020-09-28T10:56:29.997

julia> ar.name
"ar"

julia> ar.nbr
1000
```
"""
function create end

create(name::String, nbr::Int64)::Item = Item(name, nbr)

end


text:
# θ values copied from https://github.com/higham/expmv
# licensed under BSD 2-clause
# Copyright (c) 2010, Awad H. Al-Mohy and Nicholas J. Higham
# All rights reserved.
function coefficients(tol::T) where {T<:Real}
    if tol ≥ 2^(-24)
        θ = T[
            1.1920927998154304e-7,
            0.0005978858893805217,
            0.011233864735286708,
            0.05116619363445085,
            0.13084871645994703,
            0.24952893228466977,
            0.40145824235104804,
            0.5800524627688768,
            0.779511337435803,
            0.9951840790004457,
            1.2234795424241427,
            1.4616615072090335,
            1.7076485296087012,
            1.959850585959898,
            2.2170443949747205,
            2.4782808775219713,
            2.7428171126987797,
            3.0100663628176343,
            3.279561212635997,
            3.5509262147064953,
            3.823857425450966,
            4.098106972191506,
            4.3734713118405,
            4.649782224100758,
            4.926899843755911,
            5.204707228012361,
            5.483106087658634,
            5.762013408447769,
            6.041358758192571,
            6.321082126301961,
            6.601132179501162,
            6.8814648452097185,
            7.162042154487758,
            7.4428312919365975,
            7.723803811553991,
            8.004934986436286,
            8.286203267002167,
            8.567589827662578,
            8.849078185923954,
            9.13065388109011,
            9.412304202219431,
            9.694017956963032,
            9.975785274470699,
            10.257597436797516,
            10.539446734242194,
            10.821326340852181,
            11.103230206980713,
            11.385152966309168,
            11.667089855178881,
            11.949036642429446,
            12.230989568231498,
            12.51294529064492,
            12.7949008388429,
            13.07685357212606,
            13.35880114398722,
            13.64074147050525,
            13.922672702221394,
            14.204593198636742,
            14.486501505025068,
            14.768396332218838,
        ]
    elseif tol ≥ 2^(-53)
        θ = T[
            2.220446049250313e-16,
            2.5809568029946243e-8,
            1.3863478661191185e-5,
            0.00033971688399768305,
            0.002400876357887274,
            0.009065656407595102,
            0.023844555325002736,
            0.049912288711153226,
            0.08957760203223343,
            0.1441829761614378,
            0.21423580684517107,
            0.299615891381158,
            0.3997775336316795,
            0.5139146936124294,
            0.6410835233041199,
            0.7802874256626576,
            0.9305328460786568,
            1.0908637192900361,
            1.2603810606426387,
            1.438252596804337,
            1.6237159502358214,
            1.8160778162150852,
            2.014710780944616,
            2.2190488693650896,
            2.4285825244428265,
            2.6428534574594353,
            2.861449633934264,
            3.084000544989162,
            3.310172839890271,
            3.5396663487436895,
            3.772210495681751,
            4.00756108611804,
            4.245497442579696,
            4.485819859447369,
            4.728347345793539,
            4.972915626191981,
            5.219375371084058,
            5.467590630524544,
            5.717437447572013,
            5.968802630041849,
            6.221582661689891,
            6.4756827360799845,
            6.731015898381024,
            6.98750228213063,
            7.245068429597953,
            7.503646685788864,
            7.763174657377987,
            8.02359472893998,
            8.284853629803917,
            8.546902045684933,
            8.809694269971322,
            9.073187890176145,
            9.337343505612013,
            9.602124472826556,
            9.8674966757534,
            10.133428317897478,
            10.399889734191031,
            10.666853220434106,
            10.934292878475777,
            11.202184475504579,
            11.470505316002537,
            11.739234125080184,
            12.008350942053168,
            12.277837023246892,
            12.547674753126438,
            12.817847562946627,
            13.088339856203294,
            13.359136940242903,
            13.630224963455026,
            13.90159085753186,
            14.17322228433182,
            14.445107586931254,
            14.717235744490083,
            14.989596330594328,
            15.26217947477168,
            15.534975826905702,
            15.807976524300871,
            16.081173161174046,
            16.354557760369328,
            16.628122747112073,
            16.901860924634942,
            17.175765451524093,
            17.449829820647437,
            17.724047839539214,
            17.998413612126303,
            18.27292152169181,
            18.54756621498051,
            18.822342587358953,
            19.09724576895055,
            19.372271111672617,
            19.647414177108576,
            19.92267072515425,
            20.19803670338308,
            20.47350823708314,
            20.74908161993593,
            21.024753305356054,
            21.30051989866348,
            21.576378150721375,
            21.85232495499011,
            22.128357353464594,
        ]
    else
        throw(DomainError(tol, "Tolerance less than 2^(-53) is unsupported."))
    end
    return θ
end


text:
## info command format:
## info time 203 nodes 11111111 score cp 11168 pv 5e9i+ ....

function rememberPV(gs::GameStatus)
    # remember the last PV, and also the 5 previous ones because
    # they usually contain good moves to try
    i::Int = 0

    gs.lastPVLength = gs.triangularLength[1]
    #println("last length = ", gs.triangularLength[1])
    for i = 1:gs.triangularLength[1]
	gs.lastPV[i] = gs.triangularArray[1,i];
    end
end

#using IProfile
#@iprofile off

# qui search
#@iprofile begin
function Qui( gs::GameStatus, ply::Int, alpha::Int, beta::Int)
    movesfound::Int = 0
    gs.pvmovesfound = 0
    gs.triangularLength[ply+1] = ply
    teban::Int = ((ply & 0x1) == 0) ? gs.side : ((gs.side == SENTE) ? GOTE : SENTE)
    ev = (teban == SENTE) ? 1 : -1
    bestValue = -Infinity

    if (gs.inodes) & 1023 == 0
        now = time_ns()
        period::Int = now - gs.nsStart
        if period > gs.maxThinkingTime #MAXTHINKINGTIME
            gs.timedout = true
            return 0
        end
    end

    #x = -1 * ev * Eval( SENTE, gs.board, gs)
    x = ev * EvalBonanza( SENTE, gs.board, gs)
    if x >= beta
	return beta
    end
    if x > alpha
	alpha = x
    end
    if float(ply) > (gs.depth + 4.0) # changed!
	return x
    end

    movesfound = 0
    gs.pvmovesfound = 0
    gs.MoveBeginIndex = gs.moveBufLen[ply+1]
    gs.moveBufLen[ply+1+1] = generateQBB(gs.board, gs.moveBuf, teban, gs.moveBufLen[ply+1], gs)

    for i = gs.moveBufLen[ply+1]+1:gs.moveBufLen[ply+2] # 1 origin
        makeMove( gs.board, i, gs.moveBuf, teban)
        if in_check( teban, gs.board)
            #println("check!")
            takeBack( gs.board, i, gs.moveBuf, teban)
        else
            #if (seeMoveFlag(gs.moveBuf[i]) & (FLAG_NARI|FLAG_TORI)) == 0
            #    takeBack( gs.board, i, gs.moveBuf, teban)
            #    continue
            #end
            gs.inodes += 1
	    movesfound += 1
            val::Int = -Qui( gs, ply+1, -beta, -alpha)
            takeBack( gs.board, i, gs.moveBuf, teban)
            if gs.timedout
                return 0
            end
            if bestValue < val
                # both sides want to maximize from *their* perspective
                gs.pvmovesfound += 1
                gs.triangularArray[ply+1,ply+1] = gs.moveBuf[i]
		# save this move
                for j = (ply +1+1):gs.triangularLength[ply+1+1]
                    gs.triangularArray[ply+1,j] = gs.triangularArray[ply+1+1,j]
                    # and append the latest best PV from deeper plies
                end
                gs.triangularLength[ply+1] = gs.triangularLength[ply + 1+1]
                rememberPV(gs)
            end
            bestValue=max(val,bestValue)
            alpha = max( alpha, val)
            if alpha >= beta
                return beta
            end
        end
    end
    bestValue
end

# alpha-beta search
#@iprofile begin
function AlphaBeta( gs::GameStatus, ply::Int, depth::Float64, alpha::Int, beta::Int)
    movesfound::Int = 0
    gs.pvmovesfound = 0
    gs.triangularLength[ply+1] = ply
    teban::Int = ((ply & 0x1) == 0) ? gs.side : ((gs.side == SENTE) ? GOTE : SENTE)
    bestValue = -Infinity
    tt_flag::Int = TT_ALPHA
    tt_val::Int = 0
    tt_bestMove = Move(0,0,0,0,0,0)
    if depth <= 0.0
        gs.followpv = false
        return Qui( gs, ply, alpha, beta)
    end

    if (gs.inodes) & 1023 == 0
        now = time_ns()
        period::Int = now - gs.nsStart
        if period > gs.maxThinkingTime #MAXTHINKINGTIME
            gs.timedout = true
            return 0
        end
    end

    val, best, fl = tt_probe( depth, alpha, beta, gs)
    # if fl == TT_EXACT
    #     #return val
    # elseif fl == TT_ALPHA
    #     alpha = max(val,alpha)
    # elseif fl == TT_BETA
    #     beta = min(beta,val)
    # end

    if (!gs.followpv) && gs.allownull
        if true
            if in_check( teban, gs.board)
                gs.allownull = false
                gs.inodes += 1
                gs.board.nextMove ⊻= 1
                val = -AlphaBeta( gs, ply, depth - NULLMOVE_REDUCTION, -beta, -beta+1)
                gs.board.nextMove ⊻= 1
                if gs.timedout
                    return 0
                end
                gs.allownull = true
                if val >= beta
                    return val
                end
            else
            end
        end
    end

    # # IID
    # if (!gs.followpv) && gs.allownull
    #     gs.allownull = false
    #     gs.inodes += 1
    #     gs.board.nextMove ⊻= 1
    #     val = -AlphaBeta( gs, ply, depth - NULLMOVE_REDUCTION, -beta, -alpha)
    #     gs.board.nextMove ⊻= 1
    #     if gs.timedout
    #         return 0
    #     end
    #     gs.allownull = true
    #     if val >= beta
    #         return val
    #     end
    # end
    ev = (teban == SENTE) ? 1 : -1
    # tmpeval = ev * EvalBonanza( SENTE, gs.board, gs)

    #if (!gs.followpv) && gs.allownull && depth < 3.0 && ((tmpeval - 100.0 * depth) >= beta)
    # return int(tmpeval - 100.0*depth)
    #end

    movesfound = 0
    gs.pvmovesfound = 0
    gs.MoveBeginIndex = gs.moveBufLen[ply+1]
    gs.moveBufLen[ply+1+1] = generateBB(gs.board, gs.moveBuf, teban, gs.moveBufLen[ply+1], gs)

    beg = gs.moveBufLen[ply+1]+1
    for i = gs.moveBufLen[ply+1]+1:gs.moveBufLen[ply+2] # 1 origin
        if (i == beg)&&(best.move != 0)
            selectHash( gs, ply, i, depth, best)
        else
            selectmove( gs, ply, i, depth)
        end
        makeMove( gs.board, i, gs.moveBuf, teban)
        if in_check( teban, gs.board)
            #println("check!")
            takeBack( gs.board, i, gs.moveBuf, teban)
        else
            gs.inodes += 1
	    movesfound += 1
            ext = 0.0
	    if i > beg
                val = -AlphaBeta( gs, ply+1, depth-1.0+ext, -alpha-1, -alpha)
                if (val > alpha) && (val < beta)
                    val = -AlphaBeta( gs, ply+1, depth-1.0+ext, -beta, -alpha)
                end
            else
                val = -AlphaBeta( gs, ply+1, depth-1.0+ext, -beta, -alpha)
            end

            # original full depth search
            # val::Int = -AlphaBeta( gs, ply+1, depth-1.0, -beta, -alpha)

            takeBack( gs.board, i, gs.moveBuf, teban)
            if gs.timedout
                return 0
            end
            if bestValue < val
                # both sides want to maximize from *their* perspective
                gs.pvmovesfound += 1
                gs.triangularArray[ply+1,ply+1] = gs.moveBuf[i]
		# save this move
                for j = (ply +1+1):gs.triangularLength[ply+1+1]
                    gs.triangularArray[ply+1,j] = gs.triangularArray[ply+1+1,j]
                    # and append the latest best PV from deeper plies
                end
                gs.triangularLength[ply+1] = gs.triangularLength[ply + 1+1]
                rememberPV(gs)
                tt_flag = TT_ALPHA
                tt_val = alpha
                tt_bestMove = gs.moveBuf[i]
                tt_save( depth, tt_val, tt_flag, tt_bestMove, gs)
            end
            bestValue=max(val,bestValue)
            alpha = max( alpha, val)
            if alpha >= beta
                tt_flag = TT_BETA
                tt_val = beta
                tt_bestMove = gs.moveBuf[i]
                tt_save( depth, tt_val, tt_flag, gs.moveBuf[i], gs)
                return beta
            end
        end
    end
    if gs.pvmovesfound > 0
        if gs.board.nextMove == GOTE
            gs.blackHeuristics[seeMoveFrom(gs.triangularArray[ply+1,ply+1])+1,seeMoveTo(gs.triangularArray[ply+1,ply+1])+1] += Int(depth*depth)
        else
            gs.whiteHeuristics[seeMoveFrom(gs.triangularArray[ply+1,ply+1])+1,seeMoveTo(gs.triangularArray[ply+1,ply+1])+1] += Int(depth*depth)
        end
        tt_flag = TT_EXACT
        tt_save( depth, alpha, tt_flag, tt_bestMove, gs)
    end
    if movesfound == 0
        return -Infinity+ply-1
	#if in_check(teban, gs.board)
        #end
    end
    bestValue
end
#end

# search driver
function think( sengo::Int, gs::GameStatus)
    move::Move = Move(0,0,0,0,0,0) # null move
    gs.lastPV = [Move(0,0,0,0,0,0) for x = 1:MaxPly]::Array{Move,1}
    gs.lastPVLength = 0
    gs.whiteHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.blackHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.nsStart = time_ns()
    gs.inodes = 0
    gs.timedout = false
    gs.side = sengo
    gs.board.nextMove = sengo

    for IDdepth = 1.0:1.0:15.0 #MaxPly
        gs.moveBufLen = [0 for x = 1:MaxPly]::Array{Int,1}
        gs.moveBuf = [Move(0,0,0,0,0,0) for x = 1:MaxMoves]
        gs.triangularLength = [0 for x = 1:MaxPly]::Array{Int,1}
        gs.triangularArray  = [Move(0,0,0,0,0,0) for x = 1:MaxPly, y = 1:MaxPly]::Array{Move,2}
        gs.followpv = true
        gs.allownull = true
        gs.depth = IDdepth
        # score::Int = AlphaBeta(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
        score::Int = PVS(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
        # println("Depth = ", IDdepth, ", Score = ", score)
        if gs.timedout
            return gs.lastPV[1]
        end
        rememberPV(gs)
        timeInMSecs::Int = Int(div(time_ns() - gs.nsStart),1000000)
        NPS::Int = Int(div(gs.inodes,div((time_ns() - gs.nsStart),1000000000)))
        ## info time 203 nodes 11111111 score cp 11168 pv 5e9i+ ....
        print("info time ",timeInMSecs," depth ", int(IDdepth), " nodes ", gs.inodes, " score cp ", score, " nps ",NPS," pv")
        for i = 1:gs.triangularLength[1]
	    print(" ",move2USIString(gs.lastPV[i]))
        end
        println()
    end
    move = gs.lastPV[1]
    return move
end

# search driver
function thinkASP( sengo::Int, gs::GameStatus, sock)
    move::Move = Move(0,0,0,0,0,0) # null move
    gs.lastPV = [Move(0,0,0,0,0,0) for x = 1:MaxPly]::Array{Move,1}
    gs.lastPVLength = 0
    gs.whiteHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.blackHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.nsStart = time_ns()
    gs.inodes = 0
    gs.timedout = false
    gs.side = sengo
    gs.board.nextMove = sengo
    gs.tt = Dict{UInt64,TransP}()
    delta::Int = 80
    score::Int = 0
    middle::Int = 0
    smallAlpha::Int = 0
    smallBeta::Int  = 0
    for IDdepth = 1.0:1.0:30.0 #MaxPly
        hitScore::Bool = false

        while hitScore == false
            gs.moveBufLen = [0 for x = 1:MaxPly]::Array{Int,1}
            gs.moveBuf = [Move(0,0,0,0,0,0) for x = 1:MaxMoves]
            gs.triangularLength = [0 for x = 1:MaxPly]::Array{Int,1}
            gs.triangularArray  = [Move(0,0,0,0,0,0) for x = 1:MaxPly, y = 1:MaxPly]::Array{Move,2}
            gs.followpv = true
            gs.allownull = true
            gs.MoveBeginIndex = 0

            # score::Int = AlphaBeta(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
            if IDdepth == 1.0
                #score = PVS(gs, 0, IDdepth,-Infinity,Infinity,true) # ply=0
                score = AlphaBeta(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
                middle = score
                hitScore = true
                smallAlpha = middle - delta
                smallBeta  = middle + delta
            else
                #score = PVS(gs, 0, IDdepth,smallAlpha,smallBeta,true) # ply=0
                score = AlphaBeta(gs, 0, IDdepth,smallAlpha,smallBeta) # ply=0
                middle = score

                if smallAlpha == -Infinity
                    hitScore = true
                end
                if smallBeta  == Infinity
                    hitScore = true
                end

                if smallAlpha < score < smallBeta
                    hitScore = true
                elseif score <= smallAlpha
                    delta *= 2
                    smallAlpha = middle - delta
                    if delta >= 2000
                        smallAlpha = -Infinity
                    end
                elseif score >= smallBeta
                    delta *= 2
                    smallBeta = middle + delta
                    if delta >= 2000
                        smallBeta = Infinity
                    end
                end
            end

            # println("Depth = ", IDdepth, ", Score = ", score)
            if gs.timedout
                return gs.lastPV[1]
            end
            rememberPV(gs)
            timeInMSecs::Int = Int(div((time_ns() - gs.nsStart),1000000))
            frac = div((time_ns() - gs.nsStart),1000000000)
            NPS::Int = Int(div(gs.inodes,(frac == 0) ? 1 : frac ))
            ## info time 203 nodes 11111111 score cp 11168 pv 5e9i+ ....
            print(sock,"info time ",timeInMSecs," depth ", Int(IDdepth), " nodes ", gs.inodes, " score cp ", score, " nps ",NPS," pv")
            for i = 1:gs.triangularLength[1]
	        print(sock," ",move2USIString(gs.lastPV[i]))
            end
            println(sock)
        end
    end

    move = gs.lastPV[1]
    return move
end

function ThinkTest(gs::GameStatus)
    m::Move = think(SENTE,gs)
end

max_stars_repo_name: kimrin/MechajyoForWindows

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
export BinomialGeomRWModel, CRRModel, TianModel, JRModel

"""
    BinomialGeomRWModel

A model for a Binomial Geometric Random Walk (aka Binomial tree).

Arguments:

* `startdate` : start date of process
* `enddate` : start date of process
* `nsteps` : number of steps in the tree
* `S₀` : inital value
* `Δt` : the time-difference between steps, typically `days(startdate - enddate) / (365*nsteps)`
* `iR` : discount rate, `exp(-Δt*interestrate)`
* `u`  : scale factor for up
* `d`  : scale factor for down
* `p`  : up probability
* `q`  : down probability, `1-p`

"""
struct BinomialGeomRWModel{T,V,S} <: AbstractModel
    startdate::Date
    enddate::Date
    nsteps::Int
    S₀::T
    Δt::Float64
    iR::Float64 # exp(-r*Δt)
    logu::V #Float64
    logd::V #Float64
    p::S #Float64
    q::S #Float64
end

abstract type BinomialModel end
struct TianModel <: BinomialModel end
struct CRRModel <: BinomialModel end
struct JRModel <: BinomialModel end
struct JRrnModel <: BinomialModel end

numeraire(m::BinomialGeomRWModel{T}) where {T<:CurrencyQuantity} = unit(m.S₀)
numeraire(m::BinomialGeomRWModel{T}) where {T<:Real} = one(m.S₀)

"""
    CRRModel

Cox-Ross-Rubenstein binomial model.
"""
function CRRModel(startdate::Date, enddate::Date, nsteps::Int,
               startvalue, interestrate::Float64, carryrate::Float64, volatility::T) where T
    Δt = days(enddate - startdate) / (365 * nsteps)
    σ  = volatility
    b  = interestrate - carryrate
    B  = exp(b*Δt)
    iR = exp(-interestrate*Δt)

    logu = σ*√Δt
    logd = -σ*√Δt
    u = exp(logu)
    d = exp(logd)
    p = (B-d)/(u-d)
    q = (u-B)/(u-d)
    BinomialGeomRWModel(startdate,enddate,nsteps,startvalue,Δt,
                                iR,logu,logd,p,q)
end

function JRModel(startdate::Date, enddate::Date, nsteps::Int,
               startvalue, interestrate::Float64, carryrate::Float64, volatility::T) where T

    Δt = days(enddate - startdate) / (365 * nsteps)
    σ  = volatility
    b  = interestrate - carryrate
    iR = exp(-interestrate*Δt)

    logu = (b-σ^2/2)*Δt + σ*√Δt
    logd = (b-σ^2/2)*Δt - σ*√Δt
    p = q = 0.5

    BinomialGeomRWModel(startdate,enddate,nsteps,startvalue,Δt,
                                iR,logu,logd,p,q)
end
function JRrnModel(startdate::Date, enddate::Date, nsteps::Int,
               startvalue, interestrate::Float64, carryrate::Float64, volatility::T) where T
    Δt = days(enddate - startdate) / (365 * nsteps)
    σ  = volatility
    b  = interestrate - carryrate
    B  = exp(b*Δt)
    iR = exp(-interestrate*Δt)

    logu = (b-σ^2/2)*Δt + σ*√Δt
    logd = (b-σ^2/2)*Δt - σ*√Δt
    u = exp(logu)
    d = exp(logd)
    p = (B-d)/(u-d)
    q = (u-B)/(u-d)

    BinomialGeomRWModel(startdate,enddate,nsteps,startvalue,Δt,
                                iR,logu,logd,p,q)
end

# the value at step `n`, index `i` (i.e. `i` ups, `n-i` downs)
@inline valueat(m::BinomialGeomRWModel{T}, s::SingleStock, n, i) where {T} =
    m.S₀*exp(m.logu*i + m.logd*(n-i))

function value(m::BinomialGeomRWModel{T,V}, c::WhenAt{C}) where {C,T,V}
    m.enddate == maturitydate(c) || error("Binomial end date must match maturity of option")
    N = m.nsteps
    S = typeof(valueat(m, c.c, N, N))
    X = S[valueat(m, c.c, N, i) for i = 0:N]
    for n = N-1 : -1 : 0
        for i = 0:n
            X[i+1] = m.iR*(m.q*X[i+1] + m.p*X[i+2])
        end
    end
    X[1]
end
function value(m::BinomialGeomRWModel{T}, c::AnytimeBefore{C}) where {C,T}
    m.enddate == maturitydate(c) || error("Binomial end date must match maturity of option")
    N = m.nsteps
    S = typeof(valueat(m, c.c, N, N))
    X = S[valueat(m, c.c, N, i) for i = 0:N]
    for n = N-1 : -1 : 0
        for i = 0:n
            X[i+1] = max(m.iR*(m.q*X[i+1] + m.p*X[i+2]),
                         valueat(m, c.c, n, i))
        end
    end
    X[1]
end

## Tian
"""
    TianModel

Tian binomial model.
"""
function TianModel(startdate::Date, enddate::Date, nsteps::Int,
               startvalue, interestrate::Float64, carryrate::Float64, volatility::T) where T
    Δt = days(enddate - startdate) / (365 * nsteps)
    σ = volatility
    v = exp(σ^2*Δt)
    r = interestrate - carryrate
    iR = exp(-interestrate*Δt)
    u = 1/2*exp(r*Δt)*v*(v+1+sqrt(v^2+2v-3))
    d = 1/2*exp(r*Δt)*v*(v+1-sqrt(v^2+2v-3))
    logu = log(u)
    logd = log(d)
    p = (exp(r*Δt)-d)/(u-d)
    q = (u-exp(r*Δt))/(u-d)
    BinomialGeomRWModel(startdate,enddate,nsteps,startvalue,Δt,
                                iR,logu,logd,p,q)
end

# From http://www.goddardconsulting.ca/option-pricing-binomial-alts.html

## Cox-Ross-Rubinstein
# u = exp(σ*√Δt)
# d = exp(-σ*√Δt)
# p = (exp(r*Δt)-d)/(u-d)
# q = (u-exp(r*Δt))/(u-d)

## Jarrow-Rudd
# u = exp((r-σ^2/2)*Δt + σ*√Δt)
# d = exp((r-σ^2/2)*Δt - σ*√Δt)
# p = q = 0.5
# NOTE: not risk-neutral

## Jarrow-Rudd risk-neutral
# u = exp((r-σ^2/2)*Δt + σ*√Δt)
# d = exp((r-σ^2/2)*Δt - σ*√Δt)
# p = (exp(r*Δt)-d)/(u-d)
# q = (u-exp(r*Δt))/(u-d)

## Tian
# u = 1/2*exp(r*Δt)*v*(v+1+sqrt(v^2+2v-3))
# d = 1/2*exp(r*Δt)*v*(v+1-sqrt(v^2+2v-3))
#   where v = exp(σ^2*Δt)
# p = (exp(r*Δt)-d)/(u-d)
# q = (u-exp(r*Δt))/(u-d)

## CRR with drift
# u = exp(η*Δt + σ*√Δt)
# d = exp(η*Δt - σ*√Δt)
# p = (exp(r*Δt)-d)/(u-d)
# q = (u-exp(r*Δt))/(u-d)
#   for some η
#   a common choice is η = (log(K)-log(S0))/T

## Leisen-Reimer
# p = h_inv(d1)
# p2 = h_inv(d2)
# u = exp(r*Δt) * p2/p
# d = (exp(r*Δt) - p*u) / (1-p)
#  (ensures standard form for p,q)
# h_inv is a discrete approximation to Φ


text:
export rot, plshift, prshift, arshift

"""
    rot(a, b, θ)

rotate variables `a` and `b` by an angle `θ`
"""
function rot(a, b, θ)
    s, c = sincos(θ)
    a*c-b*s, a*s+b*c
end

"""
    plshift(x, n)

periodic left shift.
"""
plshift(x, n) = (x << n) | (x >> (sizeof(x)*8-n))

"""
    plshift(x, n)

periodic right shift.
"""
prshift(x, n) = (x >> n) | (x << (sizeof(x)*8-n))

"""
    arshift(x, n)

right shift, sign extending.
"""
arshift(x::T, n) where T = (x >> n) | (x & (T(1) << (sizeof(x)*8-1)))


text:
using Trilinos
using Test
using MPI
using LinearAlgebra

if Test.get_testset_depth() == 0
  MPI.Init()
  MPI.finalize_atexit()
end

comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))

const n_mu = 10000
A = 20 * sparse(1.0I, n_mu, n_mu) + sprand(n_mu, n_mu, 20.0 / n_mu)
A = (A + A') / 2

rowmap = Tpetra.Map(n_mu, 0, comm)
A = Tpetra.CrsMatrix(A, rowmap)

muelu_prec = MueLu.CreateTpetraPreconditioner(A, MueLu.parameters_easy())
H = MueLu.GetHierarchy(muelu_prec)
@test MueLu.GetNumLevels(H) == 3


text:
# test Pose2DPoint2D constraint evaluation function

using RoME, IncrementalInference, Distributions
using Base.Test

begin

N = 75
fg = initfg()

initCov = diagm([0.03;0.03;0.001])
odoCov = diagm([3.0;3.0;0.01])

# Some starting position
v1 = addNode!(fg, :x0, zeros(3,1), diagm([1.0;1.0;0.1]), N=N)
initPosePrior = PriorPose2(zeros(3,1), initCov, [1.0])
f1  = addFactor!(fg,[v1], initPosePrior)

# and a second pose
v2 = addNode!(fg, :x1, vectoarr2([50.0;0.0;pi/2]), diagm([1.0;1.0;0.05]), N=N)
ppc = Pose2Pose2(vectoarr2([50.0;0.0;pi/2]), odoCov, [1.0])
f2 = addFactor!(fg, [v1;v2], ppc)

# test evaluation of pose pose constraint
pts = evalFactor2(fg, f2, v2.index)
@test norm(Base.mean(pts,2)[1:2]-[50.0;0.0]) < 10.0
@test abs(Base.mean(pts,2)[3]-pi/2) < 0.5

# @show ls(fg)

tree = wipeBuildNewTree!(fg)
inferOverTreeR!(fg, tree,N=N)
# inferOverTree!(fg, tree, N=N)

# test post evaluation values are correct
pts = getVal(fg, :x0)
@test norm(Base.mean(pts,2)[1:2]-[0.0;0.0]) < 10.0
@test abs(Base.mean(pts,2)[3]) < 0.5

pts = getVal(fg, :x1)
@test norm(Base.mean(pts,2)[1:2]-[50.0;0.0]) < 10.0
@test abs(Base.mean(pts,2)[3]-pi/2) < 0.5

# check that yaw is working
v3 = addNode!(fg, :x2, zeros(3,1), diagm([1.0;1.0;0.05]), N=N)
ppc = Pose2Pose2(vectoarr2([50.0;0.0;0.0]), odoCov, [1.0])
f3 = addFactor!(fg, [v2;v3], ppc)

tree = wipeBuildNewTree!(fg)
[inferOverTree!(fg, tree, N=N) for i in 1:2]

# test post evaluation values are correct
pts = getVal(fg, :x0)
@test norm(Base.mean(pts,2)[1:2]-[0.0;0.0]) < 20.0
@test abs(Base.mean(pts,2)[3]) < 0.5

pts = getVal(fg, :x1)
@test norm(Base.mean(pts,2)[1:2]-[50.0;0.0]) < 20.0
@test abs(Base.mean(pts,2)[3] - pi/2) < 0.5

pts = getVal(fg, :x2)
@test norm(Base.mean(pts,2)[1:2]-[50.0;50.0]) < 20.0
@test abs(Base.mean(pts,2)[3]-pi/2) < 0.5

println("test bearing range evaluations")

# new landmark
l1 = addNode!(fg, :l1, zeros(2,1), diagm([1.0;1.0]), N=N)
# and pose to landmark constraint
rhoZ1 = norm([10.0;0.0])
ppr = Pose2DPoint2DBearingRange{Uniform, Normal}(Uniform(-pi,pi),Normal(rhoZ1,1.0))
f4 = addFactor!(fg, [v1;l1], ppr)

pts = evalFactor2(fg, f4, l1.index)
# @show sum(sqrt(sum(pts.^2, 1 )) .< 5.0)
@test sum(sqrt.(sum(pts.^2, 1 )) .< 5.0) == 0

pts = evalFactor2(fg, f4, v1.index)
# @show sum(sqrt(sum(pts.^2, 1 )) .< 5.0)
# @test sum(sqrt(sum(pts.^2, 1 )) .< 5.0) == 0

# add a prior to landmark
pp2 = PriorPoint2D([10.0;0.0], diagm([1.0;1.0]), [1.0])

f5 = addFactor!(fg,[l1], pp2)
pts = evalFactor2(fg, f5, l1.index)

@test norm(Base.mean(pts,2)[:]-[10.0;0.0]) < 5.0

println("test Pose2D plotting")

# drawPoses(fg);
# drawPosesLandms(fg);

# using KernelDensityEstimate
# using Gadfly
#
# pts = getVal(fg, :l1)
#
# p1= kde!(pts)
# plotKDE( p1 , dimLbls=["x";"y";"z"])

# plotKDE( [marginal(p1c,[1;2]);marginal(p1,[1;2])] , dimLbls=["x";"y";"z"],c=["red";"black"],levels=3)
# p1c = deepcopy(p1)

# plotKDE( marginal(getVertKDE(fg, :x2),[1;2]) , dimLbls=["x";"y";"z"])
#
# axis = [[1.5;3.5]';[-1.25;1.25]';[-1.0;1.0]']
# draw( PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
#       plotKDE( p1, dimLbls=["x";"y";"z"], axis=axis)  )
#

end

#


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

"""
    localvariography(sdata, lpars, var₁, var₂=var₁;
    axis=:X, tol=1e-6, maxratio1=Inf, maxratio2=Inf, [optional parameters])

Computes an unconventional directional (cross-)variogram for the variables `var₁`
and `var₂` of spatial data `sdata`. This uses local anisotropies `lpars` to create
pseudo-partitions along local `axis` direction (`axis` ∈ [:X, :Y, :Z]) and the
EmpiricalVariogram is calculated with these data. The directional tolerance
bandwidth `tol` can be passed. The `maxratio1` and `maxratio2` partially ignore
data if magnitude ratios are above these values (does not ignore anything by
default). Other optional parameters are defined at `GeoStats.EmpiricalVariogram`,
such as `nlags`, `maxlag`, `distance` and `algo`.

Similar (but not equal) to https://github.com/rmcaixeta/Local_variography
"""
function localvariography(obj::SpatialData, lpars::LocalAnisotropy, var₁, var₂=var₁;
                          axis=:X, tol=1e-6, maxratio1=Inf, maxratio2=Inf, kwargs...)
  p = pseudolocalpartition(obj, lpars, axis, tol, maxratio1, maxratio2)
  EmpiricalVariogram(p, var₁, var₂; kwargs...)
end

function pseudolocalpartition(obj, lpars, axis, tol, maxratio1, maxratio2)
  subs = []
  dims = ndims(lpars)
  n = nelements(obj)
  for i in 1:n
    maxratio1 != Inf && ratio1(lpars,i) > maxratio1 && continue
    maxratio2 != Inf && ratio2(lpars,i) > maxratio2 && continue
    x = coordinates(centroid(obj, i))
    v = rotmat(lpars, i)[iaxis(axis),1:dims]
    p = DirectionPartition(Tuple(v), tol=tol)
    s = Int[i]
    for j in 1:n
      i == j && continue
      y = coordinates(centroid(obj, j))
      p(x, y) && push!(s, j)
    end
    push!(subs, s)
  end

  Partition(obj, subs)
end


text:
import MittagLefflerFunctions.H
import MittagLefflerFunctions.Fplus
using PyPlot

α = 0.7
β = 0.95
tol = 1e-10

x = range(0.5, 1.5, length=201)
z = complex.(x)
errH(z) = (z^α-1) - (z-1)*(α+(z-1)*H(α, z, tol))

y = real(errH.(z))
figure(1)
plot(x, y)
grid(true)

G(z) = z^(α-β) * (z-1) / (z^α-1)
errG(z) = G(z) - z^(α-β) / (α+(z-1)*H(α,z,tol))
y = real(errG.(z))
figure(1)
figure(2)
plot(x, y)
grid(true)

y = real(Fplus.(α, β, z, tol))
figure(3)
lo = 1 - 0.1
hi = 1 + 0.1
plot(x, y)
v = axis()
plot([lo, lo, NaN, hi, hi], [v[3], v[4], NaN, v[3], v[4]], ":")
grid(true)


text:
using ProgressBars, Plots, StatsBase, Statistics, LaTeXStrings, JLD

#Diffusion-limited aggregation
function Deposing(;Network,RWpos,cap,ParticlesNumber,H1_0,H1,H2)
    Network[RWpos...,2] += 1
    if Network[RWpos...,2] == cap
        Network[RWpos...,1] = ParticlesNumber
        H1 = max(RWpos[2] + H1_0, H1)
        H2 = H1 + 5
    end
    check = false
    return Network, H1, H2, check
end

function NewRWpos(RWpos,L)
    choices = ((1,0),(0,1),(-1,0),(0,-1))
    retRWpos = RWpos .+ rand(choices)
    if retRWpos[1] == L +1
        retRWpos[1] = 1
    elseif retRWpos[1] == 0
        retRWpos[1] = L
    end
    return retRWpos
end

function DLASimulation(L, H1_0, cap)
    choices = ((1,0),(0,1),(-1,0),(0,-1))
    Frames = []
    Network = zeros(L,floor(Int,L/2),2)
    Network[:,1,1] = ones(L)
    H2_0 = H1_0 + 5
    H1 = H1_0
    H2 = H2_0
    RWpos = [rand(1:L), H1]
    ParticlesNumber = 0
    while H2 < floor(Int,L/2)
        ParticlesNumber += 1
        print("\r$ParticlesNumber   $H2     ")
        push!(Frames,Network[:,:,1])
        RWpos = [rand(1:L), H1]
        check = true
        while RWpos[2] < H2 && check
            RWpos = NewRWpos(RWpos,L)
            if RWpos[2] <= H1 - H1_0 + 2
                Parameters = Dict(:Network => Network, :RWpos => RWpos,
                                    :cap => cap, :ParticlesNumber => ParticlesNumber,
                                        :H1_0 => H1_0, :H1 => H1, :H2 => H2)
                for neighbor in ((1,0),(0,1),(-1,0),(0,-1))
                    if RWpos[1] + neighbor[1] == 0
                        if Network[L,RWpos[2] + neighbor[2],1] != 0
                            Network, H1, H2, check = Deposing(;Parameters...)
                            break
                        end
                    elseif RWpos[1] + neighbor[1] == L + 1
                        if Network[1,RWpos[2] + neighbor[2],1] != 0
                            Network, H1, H2, check = Deposing(;Parameters...)
                            break
                        end
                    elseif Network[RWpos .+ neighbor ...,1] != 0
                        Network, H1, H2, check = Deposing(;Parameters...)
                        break
                    end
                end
            end
        end
    end
    return Frames
end

L = 200
H1_0 = 10
cap = 3

Frames = DLASimulation(L, H1_0, cap)
save("../../Data/Q2/Q2-Frames.jld", "Frames", Frames)

anim = Animation()

for n in ProgressBar(1:100:length(Frames))
    frame(anim, heatmap(rotr90(Frames[n]), c = cgrad(:turbo, rev = false), legend = nothing,
        title = L"Distribution\ of\ %$n\ RW\ in\ SAW"))
end

gif(anim, "../../Figs/Q2/Q2-HM.gif", fps = 100)

savefig("../../Figs/Q2/Q2-HM.pdf")


text:
module QDates

# # Load dependencies
include(joinpath(dirname(@__FILE__), "qref", "QREF.jl"))

import Dates
import Dates:
    year,
    month,
    day,
    dayofweek,
    yearmonth,
    monthday,
    yearmonthday,
    days

const DAYS_OFFSET = 1721425
const FIRST_VALUE = QREF.FIRST_JULIAN - DAYS_OFFSET
# Note: FIRST_VALUE == 162193
const LAST_VALUE = QREF.LAST_JULIAN - DAYS_OFFSET
# Note: LAST_VALUE == 803567
const FIRST_YEAR = QREF.FIRST_YEAR
# Note: FIRST_YEAR == 445
const LAST_YEAR = QREF.LAST_YEAR
# Note: LAST_YEAR == 2200

include("types.jl")
include("periods.jl")
include("accessors.jl")
include("query.jl")
include("arithmetic.jl")
include("conversions.jl")
include("ranges.jl")
include("adjusters.jl")
include("rounding.jl")
include("io.jl")

export
    QDate,
    year,
    month,
    day,
    dayofweek,
    yearmonth,
    monthday,
    yearmonthday,
    days,
    firstdayofyear,
    lastdayofyear,
    firstdayofmonth,
    lastdayofmonth,
    today,
    isleapmonth,
    yearmonthleap,
    monthleapday,
    yearmonthleapday

end # module


text:
function C(n, k)
    if k == 0 || k == n
        return 1
    else
        return C(n - 1, k) + C(n-1, k-1)
    end
end

function get_basic_Bernstein(n, x)
    coeffs = []
    for k in 1:n
        push!(coeffs, C(n, k) * x^k * (1 - x)^(n - k))
    end
    return coeffs
end

function Bernstein_polynom(f::Function, n, x)
    coeffs = get_basic_Bernstein(n,x)
    result = 0
    for k in 1:n
        result += f(k/n) * coeffs[k]
    end
    return result
end

function Bernstein_polynom(f, n, x)
    coeffs = get_basic_Bernstein(n, x)
    result = 0
    for k in 1:n
        result += coeffs[k]*f[k]
    end
    return result
end

function test_Bernstein()
    test1 = []
    test2 = []
    for i in 0:100
        k = i / 100
        push!(test1, 2^k)
        push!(test2, Bernstein_polynom(x->2^x, 30, k))
    end
    plot(test1)
    plot(test2)
end


text:
# Autogenerated wrapper script for cddlib_jll for x86_64-w64-mingw32
export libcdd, libcddgmp

using GMP_jll
JLLWrappers.@generate_wrapper_header("cddlib")
JLLWrappers.@declare_library_product(libcdd, "libcdd-0.dll")
JLLWrappers.@declare_library_product(libcddgmp, "libcddgmp-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libcdd,
        "bin\\libcdd-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcddgmp,
        "bin\\libcddgmp-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Verilog
using Base.Test

#makes life much easier.
macro test_string(string, expr)
  :(@test string($expr) == $string)
end

#test the wire functions
include("test-wires.jl")

#test boolean operations
include("test-boolop.jl")

#test verilog generation
include("test-verigen.jl")

#test multiple outputs
include("test-multioutput.jl")

#test wire arrays
include("test-wirearray.jl")


text:
function CentralDiff(x,y)

    m=size(x,1)

    δ=zeros(m-2)
    for i=2:m-1
       δ[i-1]=(y[i-1]*((x[i]-x[i+1])/((x[i-1]-x[i])*(x[i-1]-x[i+1]))))+
              (y[i]*((2*x[i]-x[i-1]-x[i+1])/((x[i]-x[i-1])*(x[i]-x[i+1]))))+
              (y[i+1]*((x[i]-x[i-1])/((x[i+1]-x[i-1])*(x[i+1]-x[i]))))
    end
    return δ
end


text:
using EDF
using EDF: AnnotationsList, PatientID, RecordAnnotation, RecordingID, Signal
using Dates
using Test

#####
##### Testing utilities
#####

function deep_equal(a::T, b::T) where T
    nfields = fieldcount(T)
    if nfields == 0
        return isequal(a, b)  # Use `isequal` instead of `==` to handle `missing`
    else
        for i = 1:nfields
            isdefined(a, i) || return !isdefined(b, i)  # Call two undefs equal
            deep_equal(getfield(a, i), getfield(b, i)) || return false
        end
    end
    return true
end

function deep_equal(a::T, b::T) where T<:AbstractArray
    length(a) == length(b) || return false
    for (x, y) in zip(a, b)
        deep_equal(x, y) || return false
    end
    return true
end

deep_equal(a::T, b::S) where {T,S} = false

#####
##### Actual tests
#####

const DATADIR = joinpath(@__DIR__, "data")

@testset "Just Do It" begin
    edf = EDF.read(joinpath(DATADIR, "test.edf"))
    @test sprint(show, edf) == "EDF.File with 139 signals"
    @test edf.header.version == "0"
    @test edf.header.patient == PatientID(missing, missing, missing, missing)
    @test edf.header.recording == RecordingID(Date(2014, 4, 29), missing, missing, missing)
    @test edf.header.continuous
    @test edf.header.start == DateTime(2014, 4, 29, 22, 19, 44)
    @test edf.header.n_records == 6
    @test edf.header.duration == 1.0
    @test edf.header.n_signals == 139
    @test edf.signals isa Vector{Signal}
    @test length(edf.signals) == edf.header.n_signals
    for s in edf.signals
        @test length(s.samples) == s.n_samples * edf.header.n_records
    end
    expected = [
        RecordAnnotation(0.0, String[], [AnnotationsList(0.0, nothing, ["start"])], 1024),
        RecordAnnotation(1.0, String[], [AnnotationsList(0.1344, 0.256, ["type A"])], 1024),
        RecordAnnotation(2.0, String[], [AnnotationsList(0.3904, 1.0, ["type A"])], 1024),
        RecordAnnotation(3.0, String[], [AnnotationsList(2.0, nothing, ["type B"])], 1024),
        RecordAnnotation(4.0, String[], [AnnotationsList(2.5, 2.5, ["type A"])], 1024),
        RecordAnnotation(5.0, String[], AnnotationsList[], 1024),
    ]
    @test deep_equal(edf.annotations, expected)

    io = IOBuffer()
    nb = EDF.write_header(io, edf)
    @test nb == edf.header.nb_header
    EDF.write_data(io, edf)
    seekstart(io)
    h, d, i = EDF.read_header(io)
    @test deep_equal(edf.header, h)
    d, a = EDF.read_data!(io, d, h, i)
    @test eof(io)
    @test deep_equal(edf.signals, d)
    @test deep_equal(edf.annotations, a)

    mktempdir() do dir
        file = joinpath(dir, "test2.edf")
        EDF.write(file, edf)
        edf2 = EDF.read(file)
        @test deep_equal(edf, edf2)
    end

    uneven = EDF.read(joinpath(DATADIR, "test_uneven_samp.edf"))
    @test sprint(show, uneven) == "EDF.File with 2 signals"
    @test uneven.header.version == "0"
    @test uneven.header.patient == "A 3Hz sinewave and a 0.2Hz block signal, both starting in their positive phase"
    @test uneven.header.recording == "110 seconds from 13-JUL-2000 12.05.48hr."
    @test uneven.header.continuous
    @test uneven.header.start == DateTime(2000, 7, 13, 12, 5, 48)
    @test uneven.header.n_records == 11
    @test uneven.header.duration == 10.0
    @test uneven.header.n_signals == 2
    @test uneven.signals[1].n_samples != uneven.signals[2].n_samples
    @test uneven.annotations === nothing

    nonint = EDF.read(joinpath(DATADIR, "test_float_extrema.edf"))
    s = first(nonint.signals)
    @test s.physical_min ≈ -29483.1f0
    @test s.physical_max ≈ 29483.12f0
    @test s.digital_min ≈ -32767.0f0
    @test s.digital_max ≈ 32767.0f0

    # Python code for generating the comparison values used here:
    # ```
    # import mne
    # edf = mne.io.read_raw_edf("test/data/test_float_extrema.edf")
    # signal = edf.get_data()[0] * 1e6  # The 1e6 converts volts back to microvolts
    # with open("test/data/mne_values.csv", "w") as f:
    #     for x in signal:
    #         f.write("%s\n" % x)
    # ```
    mne = map(line->parse(Float32, line), eachline(joinpath(DATADIR, "mne_values.csv")))
    for (a, b) in zip(EDF.decode(s), mne)
        @test a ≈ b atol=0.01
    end
end


text:
# Uncomment to debug:
# macro dbgutil(x); esc(:(println(_dbg($x)))); end
macro dbgutil(x); end

### @primitive and @zerograd macros:

# I would like to make these type signatures as specific as possible.
# The following are not allowed yet, see https://github.com/JuliaLang/julia/issues/3766
# f{T<:Number,A<:AbstractArray{T}}(x::Value{A})
# f{T<:Number,A<:AbstractArray}(x::Value{A{T}})

"""

`@primitive fx g1 g2...` can be used to define a new primitive
and (optionally) its gradients.

Julia supports multiple dispatch, i.e. a single function can have
multiple methods with different arg types.  AutoGrad supports
multiple dispatch for primitives and gradients.  Thus fx is a
typed method declaration such as:

* @primitive sin(x::Number)
* @primitive hypot(x1::Array,x2::Array),dy,y

The second example specifies variable names for the output gradient
`dy` and the output `y` after the method declaration which can be used
in gradient expressions.  Untyped, ellipsis and keyword arguments are
ok as in `f(a::Int,b,c...;d=1)`.  Parametric methods such as
`f{T<:Number}(x::T)` cannot be used.

The @primitive macro turns the first example into:

    local sin_r = recorder(sin)
    sin{T<:Number}(x::Value{T}) = sin_r(x)

This will cause any call to `sin` with a Value{T<:Number} argument
to be recorded.  With multiple arguments things are a bit more
complicated.  Here is what happens with the second example:

    local hypot_r = recorder(hypot)
    hypot{T<:Array,S<:Array}(x1::Value{T},x2::Value{S})=hypot_r(x1,x2)
    hypot{T<:Array,S<:Array}(x1::Value{T},x2::S)=hypot_r(x1,x2)
    hypot{T<:Array,S<:Array}(x1::T,x2::Value{S})=hypot_r(x1,x2)

We want the recorder version to be called if any one of the arguments
is a boxed Value.  There is no easy way to specify this in Julia, so
the macro generates all 2^N-1 boxed/unboxed argument combinations.

The method declaration can optionally be followed by gradient
expressions.  Here are the same examples with gradients:

* @primitive sin(x::Number),dy (dy*cos(x))
* @primitive hypot(x1::Array,x2::Array),dy,y  `(dy.*x1./y)`  `(dy.*x2./y)`

Note that the parameters, the return variable and the output gradient
of the original function can be used in the gradient expressions.

In AutoGrad, gradients are defined using gradient methods that have
the following signature:

    f(Grad{i},dy,y,x...) => dx[i]

For the first example here is the generated gradient method:

`sin{T<:Number}(::Type{Grad{1}}, dy, y, x::Value{T})=(dy*cos(x))`

For the second example a different gradient method is generated for
each argument:

`hypot{T<:Array,S<:Array}(::Type{Grad{1}},dy,y,x1::Value{T},x2::Value{S})=(dy.*x1./y)`
`hypot{T<:Array,S<:Array}(::Type{Grad{2}},dy,y,x1::Value{T},x2::Value{S})=(dy.*x2./y)`

In fact @primitive generates four more definitions for the other
boxed/unboxed argument combinations.

Non-differentiable functions such as `sign`, and non-numeric functions
such as `size` should be defined using the @zerograd macro instead.

"""
macro primitive(f,g...)
    isa(f,Expr) || error("'$f' not a method signature")
    if f.head == :tuple # Using f(x),dy,y to indicate return variable for gradients
        if length(f.args) == 3
            (f,dy,y) = f.args
        elseif length(f.args) == 2
            (f,dy) = f.args; y = gensym()
        else
            error("The first arg '$f' should have the format f(x),dy,y")
        end
    else
        dy = gensym(); y = gensym()
    end
    f.head == :call || error("'$f' not a method signature")
    isa(dy,Symbol) || error("Output gradient '$dy' not a symbol")
    isa(y,Symbol) || error("Return variable '$y' not a symbol")
    b = Expr(:block)
    r = gensym()
    push!(b.args, esc(:(local $r = recorder($(fname(f))))))
    rx = rcall(r,f)
    for fx in fsigs(f)
        push!(b.args, esc(:($fx = $rx)))
        for i=1:length(g)
            gx = gsig(fx,dy,y,i)
            push!(b.args, esc(:($gx = $(g[i]))))
        end
    end
    return b
end

"""
`@zerograd f(args...; kwargs...)` allows f to handle its Value inputs
by unboxing them like @primitive, but unlike @primitive it does not
record its actions or return a Value result.  Some functions, like
sign(), have zero gradient.  Others, like length() have discrete or
constant outputs.  These need to handle Value inputs, but do not need
to record anything and can return regular values.  Their output can be
treated like a constant in the program.  Use the @zerograd macro for
those.  Note that kwargs are NOT unboxed.
"""
macro zerograd(f)
    b = Expr(:block)
    f.head == :(::) && (f=f.args[1])
    for fx in fsigs(f)
        zx = zcall(fx)
        push!(b.args, esc(:($fx = $zx)))
    end
    return b
end

function zcall(f)
    z = copy(f)
    z1 = z.args[1]
    isa(z1,Expr) && z1.head==:curly && (z.args[1]=z1.args[1])
    for i=2:length(z.args)
        zi = z.args[i]
        if isa(zi,Symbol)
            # all done
        elseif !isa(zi,Expr)
            error("Unrecognized argtype '$zi'")
        elseif zi.head==:(::)
            (v,t) = zi.args
            if t==:Value || (isa(t,Expr) && t.head==:curly && t.args[1]==:Value)
                z.args[i] = :($v.value)
            else
                z.args[i] = v
            end
        elseif zi.head==:(...)  # done
        elseif zi.head==:parameters # done
        else
            error("Unrecognized argtype '$zi'")
        end
    end
    return z
end

# get name out of function declaration
function fname(f)
    n = f.args[1]
    isa(n,Expr) && n.head==:curly && error("parametric methods not currently supported")
    if isa(n,Symbol)
        return n
    else
        error("$n not a symbol")
    end
end

# create call to r using typeless argument of f
function rcall(r,f)
    rx = notypes(f)
    rx.args[1]=r
    # Need to fix kwargs
    r2 = rx.args[2]
    if isa(r2,Expr) && r2.head == :parameters
        for i in 1:length(r2.args)
            k = r2.args[i]
            if !isa(k,Expr); error("Bad kwarg '$k'")
            elseif k.head == :(...); continue
            elseif k.head != :kw; error("Bad kwarg '$k'")
            elseif !isa(k.args[1],Symbol); error("Bad kwarg '$k'")
            else; k.args[2]=k.args[1]; end
        end
    end
    return rx
end

# eliminate type declarations from a function call
function notypes(ex)
    if isa(ex, Expr)
        if (ex.head == :(::) || ex.head == :curly)
            return notypes(ex.args[1])
        else
            return Expr(ex.head, map(notypes, ex.args)...)
        end
    else
        return ex
    end
end

# create type signatures for f where one or more args are Nodes.
function fsigs(f)
    f1 = copy(f)
    a1 = f1.args[1] = Expr(:curly,fname(f1))
    nargs = 0
    for i=2:length(f1.args)
        ai = f1.args[i]
        if isa(ai,Symbol)
            nargs+=1
            ti = gensym()
            push!(a1.args, Expr(:<:, ti, Any))
            f1.args[i] = Expr(:(::),ai,ti)
        elseif !isa(ai,Expr)
            error("Neither Symbol nor Expr: $ai")
        elseif in(ai.head, (:parameters, :(...)))
            continue
        elseif ai.head == :(::)
            nargs+=1
            ti = gensym()
            push!(a1.args, Expr(:<:,ti,ai.args[2]))
            ai.args[2] = ti
        else
            error("Argtype not supported: '$ai'")
        end
    end
    flist = []
    for nodes=0:(1<<nargs-2)
        fn = copy(f1)
        iargs = 0
        for i=2:length(fn.args)
            ai = fn.args[i]
            in(ai.head, (:parameters, :(...))) && continue
            ai.head == :(::) || error("Bad arg '$ai'")
            if nodes & (1<<iargs) == 0
                ai.args[2] = Expr(:curly,:Value,ai.args[2])
            end
            iargs += 1
        end
        push!(flist, fn)
    end
    return flist
end

function gsig(f,dy,y,i)
    g = copy(f)
    if g.args[2].head == :parameters; a = 3; else; a = 2; end
    insert!(g.args, a, :(::Type{Grad{$i}}))
    insert!(g.args, a+1, dy)
    insert!(g.args, a+2, y)
    return g
end

### Testing Utilities:

if !isdefined(:runtests)
let tests=[]
    global addtest,runtests,alltests
    alltests()=tests
    addtest(t...)=push!(tests,t)
    function runtests(a=tests)
        for fx in a
            try
                tx = fixtest(fx...)
                check_grads(tx...; fname=fx[1])
            catch e
                warn((fx...,"$e"))
            end
        end
    end
end
end

function fixtest(f, x...)
    f = eval(f)
    y = f(x...)
    # detect and prevent testing of zero / undefined grads
    plist = Any[]               # define fnew(plist)
    alist = Any[x...]           # to return f(alist)
    fargs = Any[]               # call fnew(fargs...)
    gargs = (Value(y), Value(y), map(Value,x)...)
    for i=1:length(alist)
        g = nothing
        try
            g = f(Grad{i},gargs...)
        catch e
            if isa(e,MethodError) && e.f === f && e.args[1] === Grad{i}
                continue        # warn("No grad $i for $f: $e")
            else
                error("Error during $f$((Grad{i},gargs...)): $e")
            end
        end
        g == nothing && continue # zero grads
        push!(fargs, alist[i])
        alist[i] = Symbol("x$i")
        push!(plist, alist[i])
    end
    isempty(fargs) && error("$f has no differentiable arguments.")
    f1=f; f = eval(Expr(:->, Expr(:tuple, plist...), Expr(:call, f1, alist...)))
    # if f has non-scalar output, sum it
    isbits(y) || (f2=f; f=(x...)->toscalar(f2(x...)))
    return (f,fargs...)
end

function randin(range, dims...; eps=EPS)
    if isa(range, UnitRange{Int64})
        rand(range, dims...)
    elseif range==(-Inf,Inf)
        randn(dims...)
    elseif range==(0,Inf)
        eps-log(rand(dims...))
    elseif range==(1,Inf)
        eps+1-log(rand(dims...))
    elseif range==(-1,Inf)
        eps-1-log(rand(dims...))
    elseif range==(-1,1)
        (1-eps)*(2rand(dims...)-1)
    elseif range==(0,1)
        eps+(1-2eps)*rand(dims...)
    elseif range==(0,2)
        eps+2*(1-eps)*rand(dims...)
    elseif range==(-Inf,-1,1,Inf)
        x = sec(randn(dims...))
        sign(x)*eps + x
    else
        error("Unknown range $range")
    end
end

function addtest1(f,r)          # unary
    addtest(f,randin(r))
    addtest(f,randin(r,2))
end

function addtest2(f,r1,r2=r1)   # binary
    addtest(f,randin(r1),randin(r2))
    addtest(f,randin(r1),randin(r2,2))
    addtest(f,randin(r1,2),randin(r2))
    addtest(f,randin(r1,2),randin(r2,2))
end

function addtest3(f,r1,r2=r1)   # broadcasting
    addtest2(f,r1,r2)
    addtest(f,randin(r1,2),randin(r2,2,2))
    addtest(f,randin(r1,2,2),randin(r2,2))
    addtest(f,randin(r1,1,2),randin(r2,2,2))
    addtest(f,randin(r1,2,2),randin(r2,1,2))
end

# EPS, RTOL, ATOL = 1e-4, 1e-4, 1e-6
EPS, RTOL, ATOL = 1e-4, 1e-2, 1e-4

# TODO: do sampling or random direction for large args
"""
check_grads(fun, args...) checks the computed gradients for fun(args)
comparing them with numeric approximations.
"""
function check_grads(fun, args...; eps=EPS, rtol=RTOL, atol=ATOL, fname=fun)
    @dbgutil((:check_grads,fname,:args,args...))
    isempty(args) && error("No args given")
    exact = ntuple(i->grad(fun,i)(args...), length(args))
    numeric = nd(fun, args...; eps=eps)
    @dbgutil((:check_grads,fname,:exact,exact,:numeric,numeric))
    same = isequivalent(exact, numeric; rtol=rtol, atol=atol)
    same || warn((:check_grads,fname,:args,args,:exact,exact,:numeric,numeric))
    return same
end

function nd(f, args...; eps=EPS)
    @dbgutil((:nd,f,args..., :eps, eps))
    unary_f = x->f(x...)
    unary_nd(unary_f, args, eps)
end

unary_nd(f, x::Tuple, eps)         = ntuple(i->unary_nd(indexed_function(f, x, i), x[i], eps), length(x))
unary_nd(f, x::Associative, eps)   = (a=similar(x); for(k,v) in x; a[k] = unary_nd(indexed_function(f, x, k), v, eps); end; a)
unary_nd(f, x::AbstractArray, eps) = reshape(eltype(x)[unary_nd(indexed_function(f, x, i), v, eps) for (i,v) in enumerate(x)], size(x))
unary_nd(f, x::Complex, eps)       = ((f(x + eps/2) - f(x - eps/2)) / eps - im*(f(x + im*eps/2) - f(x - im*eps/2)) / eps)
unary_nd(f, x::Real, eps)          = ((f(x + eps/2) - f(x - eps/2)) / eps)

function indexed_function(fun, arg, index)
    function partial_function(x)
        if isa(arg, Tuple)
            local_arg = (arg[1:index-1]..., x, arg[index+1:end]...)
        else
            local_arg = copy(arg); local_arg[index] = x
        end
        return fun(local_arg)
    end
    return partial_function
end

# isequivalent uses isapprox for Number and AbstractArray{T<:Number}
isequivalent(x::Number,y::Number; o...)=isapprox(x,y;o...)
isequivalent{T<:Number,S<:Number}(x::AbstractArray{T},y::AbstractArray{S}; o...)=isapprox(x,y;o...)

# isequivalent extends to Tuple, Associative, and other Arrays, comparing elementwise
isequivalent(x::Tuple, y::Tuple; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))
isequivalent(x::AbstractArray, y::AbstractArray; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))
isequivalent(x::Associative, y::Associative; o...)=all(k->isequivalent(get(x,k,nothing),get(y,k,nothing);o...), unique([keys(x)...,keys(y)...]))

# isequivalent treats `nothing` as equivalent to zero or zero array.
isequivalent(x::Number,z::Void; o...)=isequivalent(z,x;o...)
isequivalent{T<:Number}(x::AbstractArray{T},z::Void; o...)=isequivalent(z,x;o...)
isequivalent(z::Void,x::Number; o...)=isapprox(zero(x),x;o...)
isequivalent{T<:Number}(z::Void,x::AbstractArray{T}; rtol::Real=Base.rtoldefault(T), atol::Real=0, norm::Function=vecnorm) = (norm(x) <= atol/(1-rtol)) # Modified from: linalg/generic.jl:522

# The way broadcasting works in Julia:
# y = f(x...) where f is a broadcasting operation.
# size(y) = broadcast_shape(x...)
# ndims(y) = max ndims(x)
# size(y,i) = max size(x,i)
# size(x,i) = 1 or size(y,i) for all x and i<=ndims(x)
# if ndims(x) < ndims(y) the extra dimensions of x are treated as 1

function unbroadcast(x, dx)
    if size(x)==size(dx)
        return dx
    elseif isa(getval(x),Number)
        return sum(dx)
    else
        d = []
        for i=1:ndims(dx)
            size(x,i) == size(dx,i) && continue
            size(x,i) != 1 && throw(DimensionMismatch())
            push!(d,i)
        end
        length(d)==1 && (d=d[1])
        return sum(dx, d)
    end
end

function toscalar(xv; rng=MersenneTwister())
    x = getval(xv)
    isa(x,Number) && return xv
    isa(x,OneHot) && (x = full(x))
    idx = isa(x,Tuple) ? (1:length(x)) : eachindex(x)
    s = 0
    for i in idx
        s += xv[i] * rand(rng)
    end
    return s
end

# sumvalues sums values of dictionaries, otherwise acts like sum:

sumvalues(x)=sum(x)
sumvalues(x::Associative)=sum(values(x))
@primitive sumvalues(x::Associative),ds fillvalues(ds,x)
fillvalues(v,x)=(y=similar(x);for k in keys(x); y[k]=v; end; y)
@primitive fillvalues(v,x),dxv sumvalues(dxv) nothing
addtest(sumvalues, Dict(1=>1.,2=>2.))
addtest(fillvalues, 0., Dict(1=>1.,2=>2.,3=>3.))

# This needs more work:
# @primitive values(x),dy Dict(map((a,b)->(a=>b), keys(x), dy))

# It gets tiresome to write `Type{Grad{1}}` after a while, here are
# some convenient aliases:

typealias D1 Type{Grad{1}}
typealias D2 Type{Grad{2}}
if !isdefined(:Dn)
typealias Dn{N} Type{Grad{N}}
end

# Pretty print for debugging:
_dbg(x)=x # extend to define short printable representations
_dbg(x::Tuple)=map(_dbg,x)
_dbg(x::Node)="N$(id2(x))_$(id2(x.value))"
_dbg(x::Value)="V$(id2(x))_$(_dbg(x.value))"
_dbg(x::Tape)="T$(join([id2(x),map(id2,x)...],'_'))"
_dbg(x::AbstractArray)="A$(join([id2(x),size(x)...],'_'))"
id2(x)=Int(object_id(x)%1000)

Base.show(io::IO, n::Value) = print(io, _dbg(n))
Base.show(io::IO, n::Node) = print(io, _dbg(n))
Base.show(io::IO, n::Tape) = print(io, _dbg(n))