text:
# ==============================  LOAD DATA =========================================
#save("data.jld", "data", df_data, compress=true)
#load("data.jld")["data"]

#using DataStructures, DataFrames, DataFramesMeta, StatsFuns, GLM, JuMP, NLopt, HDF5, JLD, Distributions, MixedModels, RCall, StatsBase, xCommon
#using StatStack
using DataFrames

"""
    loadDF(dir::String)

Create a DataFrame from the `csv_final` file taking column names from the `Headers` file
"""
function loadDF(dir::String)
    isdir(dir) || throw(ArgumentError("\"$dir\" is not a directory"))
    fnms = readdir(dir)
    headnm = filter(r"^Headers", fnms)
    if length(headnm) != 1
        throw(ArgumentError("\"$dir\" should contain exactly one file starting with \"Headers\""))
    end
    dfrm = readtable(joinpath(dir, filter(r"^final_csv", fnms)[1]), header = false,
        separator = ',',  names = Symbol.(lowercase.(chomp.(readlines(joinpath(dir, headnm[1]))))),
        makefactors = true)
    dfrm = dfrm[convert(Vector{Bool}, map(x -> length(unique(x)) > 1, dfrm.columns))]
#    map!(x -> (isa(x, PooledDataArray) && !anynull(x)) ? pool(map(e -> e.value, x)) : x,
#        dfrm.columns)
    map!(x -> (anyna(x) || isa(x, PooledDataArray)) ? x : Array(x), dfrm.columns)
    dfrm
end

"""
    isValid(df_data::DataFrame, cfg::OrderedDict)

Predicate to check `df_data` against `cfg` for validity
"""
function isValid(df_data::DataFrame, cfg::Dict)
    dfnames = names(df_data)
    if !all(x -> x in dfnames, cfg[:all_mandatory_vars])
        throw(ArgumentError("Not all cfg[:all_mandatory_vars] in dataset"))
    end
    if !all(x -> x in dfnames, cfg[:scoring_vars])
        throw(ArgumentError("Not all cfg[:scoring_vars] in dataset"))
    end
    true
end

Base.lowercase(v::Vector{Symbol}) = Symbol.(lowercase.(string.(v)))
Base.filter(r::Regex, v::Vector{Symbol}) = Symbol.(filter(r, string.(v)))
"""
    getCFG(df_in::DataFrame, basedir::String)

Create the configuration Dict from defaults and the `app.json` file in `basedir`
"""
function getCFG(df_in::DataFrame, basedir::String)
    cfg = Dict(
        :P2_Competitor => true,
        :pvalue_lvl => 0.20,  #pvalue_lvl = 0.20
        :excludedBreaks => String[],    #["estimated_hh_income","hh_age","number_of_children_in_living_un","person_1_gender"]
        :excludedLevels => ["none"],
        :excludedKeys => String[],
        :exposed_flag_var => :exposed_flag,
        :sigLevel => "0.2",  # why a string and not a number.  It is easier to convert a number to a string than vice versa.
        :random_demos => [:estimated_hh_income, :hh_age, :number_of_children_in_living_un, :person_1_gender],
        :random_campaigns => Symbol[],
        :dropvars => [:exposed_flag],
        :scoring_vars => lowercase([:Prd_1_Net_Pr_PRE, :Prd_1_Net_Pr_POS, :Buyer_Pos_P0, :Buyer_Pre_P0]),
        :all_mandatory_vars => lowercase([:Buyer_Pos_P1, :Buyer_Pre_P1, :Trps_POS_P1, :Trps_PRE_P1,
            :Dol_per_Trip_POS_P1, :Dol_per_Trip_PRE_P1, :Nonbuyer_Pre_P1,
            :hh_age, :estimated_hh_income, :number_of_children_in_living_Un, :person_1_gender, :experian_id]),
        )
    override = joinpath(basedir, "app.json")
    if isfile(override)
        for (k, v) in JSON.parsefile(override)
            ksym = Symbol(k)
            if haskey(cfg, ksym)
                current = cfg[ksym]
                if isa(current, Dict{Symbol, Any})
                    cfg[ksym] = union(current, convert(Vector{Symbol}, v))
                else
                    cfg[ksym] = v
                end
            else
                cfg[ksym] = v
            end
        end
    end
    cfg[:exposed_flag_var] = :exposed_flag_new                  # Go In app.cfg
    cfg[:random_campaigns] = [:publisher_fct1, :targeting_fct1]  # Go In app.cfg
    cfg[:allrandoms] = vcat(cfg[:random_demos], cfg[:random_campaigns])
    cfg[:proscore] = filter(r"model", names(df_in))   # Set the ProScore variable in the dataset   #push!(cfg[:all_mandatory_vars], cfg[:ProScore])
    isempty(cfg[:proscore]) && throw(ArgumentError("Missing response variable in data"))
    cfg[:num_products] = length(filter(r"^buyer_pre_p", names(df_in))) - 1  # get number of products in the data
    cfg[:random_campaigns] = intersect(cfg[:random_campaigns], names(df_in))
    cfg
end

function reworkCFG!(df_in::DataFrame, cfg::Dict)
    dnames = names(df_in)
    xflags = [cfg[:exposed_flag_var], :exposed_flag]
    features = setdiff(dnames, xflags)

    cfg[:xVarsDemos] = vcat([:experian_id, :banner, :prd_0_qty_pre],
        features[findfirst(features, :state):findfirst(features, :mosaic)])   # exclude demos between....
    cfg[:xVarsDemos]=setdiff(cfg[:xVarsDemos],[:person_1_gender,:number_of_children_in_living_Un]) # exclude person, # from non-used demos
    cfg[:xVarsPost] = grep(["POS","Pos","Buyer_Pre_","Nonbuyer_Pre_"],features)   #exclude POST variables
    cfg[:iVarsPREPOS] = grep(["PRE_POS","Pre_Pos"],features)
    features=setdiff(features, setdiff(  vcat(cfg[:xVarsDemos],cfg[:xVarsPost]) ,cfg[:iVarsPREPOS])  )

    cfg[:xVarsP0] =  setdiff(grep("0", features ),[cfg[:ProScore]] ) #exclude category variables(exclude P0's)
    features=setdiff(features,cfg[:xVarsP0])

    cfg[:xVarsReports] = grep(["perc_","Pr_per_"],features)  #exclude Reporting vars
    features=setdiff(features,cfg[:xVarsReports])

    cfg[:xvars] = vcat(xflags,cfg[:xVarsDemos], cfg[:xVarsPost],cfg[:xVarsP0], cfg[:xVarsReports],cfg[:dropvars])
    features=setdiff(features,cfg[:dropvars])

    cfg[:ivars] = vcat(cfg[:iVarsPREPOS],cfg[:all_mandatory_vars],cfg[:scoring_vars])

    cfg[:ALL_vars_to_exclude] = setdiff(cfg[:xvars], cfg[:ivars])

    df_in[:group] = df_in[cfg[:exposed_flag_var]]
    df_in[:panid] = df_in[:experian_id]
    features = setdiff(unique(vcat(features,cfg[:iVarsPREPOS],cfg[:all_mandatory_vars],cfg[:scoring_vars],[:group,:panid])  ) ,[:experian_id] )

    cfg[:negativevars] = grep(vec(hcat([["P"*string(i),string(i)*"_"] for i=3:cfg[:num_products]]...)),features)   # get variables that need to have negative sign
    cfg[:positivevars] = grep(["P1", "1_","MODEL"], features)    # get variables that need to have positive sign

    if cfg[:P2_Competitor] == true
        cfg[:negativevars] = unique(vcat(cfg[:negativevars],grep(["P2","2_"],features)))
    else
        cfg[:positivevars] = unique(vcat(cfg[:positivevars],grep(["P2","2_"],features)))
    end

    return df_in[features]
end

#df_data = reworkCFG!(df_in,cfg)

#######################################
#--------- Data Manipulation ---------#
#######################################

function data_Prep(df_data::DataFrame, cfg::OrderedDict)
    df_data[:isO]=false
    df_data[:whyO]=""

    df_data[df_data[:number_of_children_in_living_Un].>=4,:number_of_children_in_living_Un] = 4  # aggregate #of children for 4+ L_114
    #vars = getvars()   # :GROUP
    #vars = DataFrame(names=names(df_data),eltypes=eltypes(df_data))
    if typeof(df_data[:group]) in [DataArray{UTF8String,1},DataArray{ASCIIString,1}]
        df_data[ findin(df_data[:group],["//N","\\N"]), :group] = "0"
        df_data[isna(df_data[:group]), :group]="0"
        df_data[:group] = [parse(Int64,s) for s = df_data[:group]]
    else
        df_data[ isnan(df_data[:group]), :group] = 0
    end

    vars = DataFrame(names=names(df_data),eltypes=eltypes(df_data))
    for c in setdiff(vars[findin(vars[:eltypes],[UTF8String]),:names],cfg[:allrandoms]) # set variables as.numeric and replace NA's with zero
        print("Convert String: String->Numeric: ",c)
        try df_data[c] = map( x -> isna(x) ?  NaN : convert(Float64, x)  , df_data[c]) catch e println("  (failed)") end  #NA to NaN
        try df_data[c] = convert(Array{Float64}, df_data[c]) catch e end   # To Float64
    end
    vars = DataFrame(names=names(df_data),eltypes=eltypes(df_data))
    #vars = getvars()
    for c in setdiff(vars[findin(vars[:eltypes],[Float64]),:names],cfg[:allrandoms])  # replace NaN's with zero for numeric variables
        println("Replace Float64 NaN (0.0): ",c)
        df_data[ isnan(df_data[c]), c] = 0.0
    end
    for c in  setdiff(vars[findin(vars[:eltypes],[Int64]),:names],cfg[:allrandoms])
        println("Replace Int64 NaN (0) : ",c)
        df_data[ isnan(df_data[c]), c] = 0
    end

    # create factor variables -  Need to understand how this impacts the data passed to R - might need to do this in R
    #pool!(df_data, cfg[:allrandoms])
    #pool!(df_data, [,:group])

    # NOTE : Do QCs here
    df_data[df_data[:person_1_gender].=="U",:isO] = true   # remove HHs with no gender info
    df_data[df_data[:person_1_gender].=="U",:whyO] = "person_1_gender=U; remove HHs with no gender info"
    df_data[findin(df_data[:estimated_hh_income],["U","L"]),:estimated_hh_income]="L" # aggregate U and L levels of hh income

    for r in cfg[:random_campaigns]    # check and drop exposed HHs with no publisher info or non-exposed HHs with publisher info
        df_data[findin(df_data[r],["\\N","NULL","0","NONE"])  ,r] ="none"
        df_data[ (df_data[:isO].==false) & (df_data[r].!="none") & (df_data[:group].==0) ,:whyO] = "non exposed HHs with publisher info"
        df_data[ (df_data[:isO].==false) & (df_data[r].!="none") & (df_data[:group].==0) ,:isO] = true
        println(r," non exposed HHs with publisher info : ",nrow(df_data[df_data[:whyO].=="non exposed HHs with publisher info",:] )   )
        df_data[ (df_data[:isO].==false) & (df_data[r].=="none") & (df_data[:group].!=0) ,:whyO] = "exposed HHs with no publisher info"
        df_data[ (df_data[:isO].==false) & (df_data[r].=="none") & (df_data[:group].!=0) ,:isO] = true
        println(r," exposed HHs with no publisher info : ",nrow(df_data[df_data[:whyO].=="exposed HHs with no publisher info",:] )   )
    end

    # segments for outliers detection
    df_data[:data_NB_NE_B] = false
    df_data[ (df_data[:Buyer_Pre_P1].==0 ) & (df_data[:group].==0 ) & (df_data[:Buyer_Pos_P1].==1 ) ,:data_NB_NE_B] = true
    df_data[:data_B_E_NB] = false
    df_data[ (df_data[:Buyer_Pre_P1].==1 ) & (df_data[:group].==0 ) & (df_data[:Buyer_Pos_P1].==0 ) ,:data_B_E_NB] = true
    df_data[:pen_reduction] = false
    df_data[ (df_data[:Buyer_Pre_P1].==1) & (df_data[:group].==0) & (df_data[:Buyer_Pos_P1].==0 )  ,:pen_reduction] = true
    df_data[:occ_reduction] = false
    df_data[ (df_data[:group].==0) & (df_data[:Buyer_Pos_P1].==1) & (df_data[:Trps_POS_P1].< df_data[:Trps_PRE_P1] ) ,:occ_reduction] = true
    df_data[:dolocc_reduction] = false
    df_data[  (df_data[:group].==0) & (df_data[:Buyer_Pos_P1].==1) & (df_data[:Dol_per_Trip_POS_P1].< df_data[:Dol_per_Trip_PRE_P1] )  , :dolocc_reduction] = true

    return df_data[df_data[:isO].==false, : ]   #[setdiff(names(df_data),[:isO,:whyO])]
end

#df_data = data_Prep(df_data, cfg);

function MatchMe(df_data::DataFrame,cfg::OrderedDict)
    df=df_data[df_data[:isO].==false,:]
    df_exp     = df[df[:group].==1,:]
    df_unexp   = df[df[:group].==0,:]
    df_exp_dim   = nrow(df_exp)
    df_unexp_dim = nrow(df_unexp)
    new_unexp_dim = df_unexp_dim*(df_unexp_dim>2000000 ? 0.3 : df_unexp_dim>1000000 ? 0.4 : df_unexp_dim>750000 ? 0.6 : 0.7)
    #df_sample=DataFrame
    if length(string(cfg[:ProScore])) == 0
        df_unexp_1 =  df[(df[:group].==0)&(df[:Buyer_Pre_P1].==1),:]
        df_unexp_0 =  df[(df[:group].==0)&(df[:Buyer_Pre_P1].==0),:]

        df_exp_1_dim = nrow(df[(df[:group].==1)&(df[:Buyer_Pre_P1].==1),:])
        df_exp_0_dim = nrow(df[(df[:group].==1)&(df[:Buyer_Pre_P1].==0),:])
        df_unexp_1_dim = nrow(df_unexp_1)
        df_unexp_0_dim =  nrow(df_unexp_0)
        dim_sample0 = round(Int64, (new_unexp_dim-df_exp_dim  ) / (1+(df_exp_1_dim / df_exp_0_dim)) )
        dim_sample1 = round(Int64, new_unexp_dim - df_exp_dim - dim_sample0)

        new_df_unexp_1 = df_unexp_1[sample(1:size(df_unexp_1,1), dim_sample1 ),:]
        new_df_unexp_0 = df_unexp_0[sample(1:size(df_unexp_0,1), dim_sample0 ),:]
        df_data_sample  = vcat(df_exp,new_df_unexp_1,new_df_unexp_0)
    elseif length(string(cfg[:ProScore])) > 0
        sample_control_data=similar(df_unexp, 0)
        for (key, value) in countmap(df_exp[cfg[:ProScore]])
            sample_dim=round(Int64,new_unexp_dim*(value/df_exp_dim))
            temp_data = df_unexp[df_unexp[cfg[:ProScore]].==key,:]
            samp_data = temp_data[sample(1:size(temp_data,1), sample_dim, replace=false),:]
            sample_control_data = vcat(sample_control_data,    samp_data   )
        end

        sample_data = vcat(sample_control_data,df_exp)
        sample_df_unexp = sample_data[sample_data[:group].==0,:]
        sample_df_exp   = sample_data[sample_data[:group].==1,:]
        sample_df_unexp_1 = sample_data[(sample_data[:group].==0)&(sample_data[:Buyer_Pre_P1].==1),:]
        sample_df_unexp_0 = sample_data[(sample_data[:group].==0)&(sample_data[:Buyer_Pre_P1].==0),:]
        sample_df_exp_1 =  sample_data[(sample_data[:group].==1)&(sample_data[:Buyer_Pre_P1].==1) ,:]
        sample_df_exp_0 =  sample_data[(sample_data[:group].==1)&(sample_data[:Buyer_Pre_P1].==0) ,:]
        sample_df_unexp_1_dim = nrow(sample_df_unexp_1)
        sample_df_unexp_0_dim = nrow(sample_df_unexp_0)
        sample_df_exp_1_dim = nrow(sample_df_exp_1)
        sample_df_exp_0_dim = nrow(sample_df_exp_0)
        dim_sampleA = round(Int64,(sample_df_exp_1_dim/sample_df_exp_0_dim)*sample_df_unexp_0_dim)
        dim_sampleB = round(Int64,(sample_df_exp_0_dim/sample_df_exp_1_dim)*sample_df_unexp_1_dim)

        if sample_df_unexp_1_dim/sample_df_unexp_0_dim > sample_df_exp_1_dim/sample_df_exp_0_dim
            new_df_unexp_1 = sample_df_unexp_1[sample(1:sample_df_unexp_1_dim,dim_sampleA , replace=false),:]
            df_data_sample = vcat(sample_df_exp,new_df_unexp_1,sample_df_unexp_0)
        else
            new_df_unexp_0 = sample_df_unexp_0[sample(1:sample_df_unexp_0_dim, dim_sampleB, replace=false ),:]
            df_data_sample = vcat(sample_df_exp,sample_df_unexp_1,new_df_unexp_0)
        end
    end

    #Generate Report!!!! ---   still need to remove :isO
    function dict2df(d::Dict,kc::AbstractString,c::AbstractString)
        dfo=DataFrame(k=Int64[],v=Int64[])
        for (key, value) in d
            push!(dfo,[key, value])
        end
        rename!(dfo,:k,symbol(kc))
        rename!(dfo,:v,symbol(c))
        return dfo
    end
    ds0=dict2df(countmap(df_data_sample[ df_data_sample[:group].==0,cfg[:ProScore]]),"ProScore","out_unexp")
    ds1=dict2df(countmap(df_data_sample[ df_data_sample[:group].==1,cfg[:ProScore]]),"ProScore","out_exp")
    do0=dict2df(countmap(df[ df[:group].==0,cfg[:ProScore]]),"ProScore","in_unexp")
    do1=dict2df(countmap(df[ df[:group].==1,cfg[:ProScore]]),"ProScore","in_exp")
    dfo=join(join(join(do0, do1, on = :ProScore,kind=:outer), ds0, on = :ProScore,kind=:outer),ds1, on = :ProScore,kind=:outer)
    for c in [:out_unexp, :out_exp, :in_unexp, :in_exp]
        dfo[isna(dfo[c]), c] = 0
    end
    b_ds0=dict2df(countmap(df_data_sample[ df_data_sample[:group].==0,:Buyer_Pre_P1]),"Buyer_Pre_P1","out_unexp")
    b_ds1=dict2df(countmap(df_data_sample[ df_data_sample[:group].==1,:Buyer_Pre_P1]),"Buyer_Pre_P1","out_exp")
    b_do0=dict2df(countmap(df_data[ df_data[:group].==0,:Buyer_Pre_P1]),"Buyer_Pre_P1","in_unexp")
    b_do1=dict2df(countmap(df_data[ df_data[:group].==1,:Buyer_Pre_P1]),"Buyer_Pre_P1","in_exp")
    dfx=DataFrame(Buyer_Pre_P1=[0,1],Buyer_Pre_Period=["No","Yes"])
    b_dfo1=join(b_do0, b_do1, on = :Buyer_Pre_P1,kind=:outer)
    b_dfo1[:in_proportion] = b_dfo1[:in_exp]./sum(b_dfo1[:in_exp])
    b_dfo2=join(b_ds0, b_ds1, on = :Buyer_Pre_P1,kind=:outer)
    b_dfo2[:out_proportion] =b_dfo2[:out_exp]./sum(b_dfo2[:out_exp])
    b_dfo=join(dfx,join(b_dfo1, b_dfo2, on = :Buyer_Pre_P1,kind=:outer), on = :Buyer_Pre_P1,kind=:outer)
    delete!(b_dfo, :Buyer_Pre_P1 )
    #still need to Write QC to file
    #dfo  b_dfo

    rows2remove = setdiff(df_data[df_data[:isO].==false, :panid],df_data_sample[:panid])
    df_data[findin(df_data[:panid],rows2remove),:whyO]="NoMatch"
    df_data[findin(df_data[:panid],rows2remove),:isO]=true

    return df_data[df_data[:isO].==false, : ][setdiff(names(df_data),[:isO,:whyO])]
end

#df_data = MatchMe(df_data,cfg)

#lowercase!(df_data)
#cfg=lowercase(cfg)

abstract MModel

#df_data[:iso]=false
#df_data=df_data[(df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1),:]

######################################
#--------------OCCASION--------------#
######################################

#!!! Multicollinearity, removing :[:fea_or_dis_trps_shr_dpp_p1,:fea_or_dis_trps_shr_dpp_p2,:fea_or_dis_trps_shr_dpp_p3,:fea_or_dis_trps_shr_dpp_p4]
type xGLM_OLD
    vfactors::Vector{Symbol}
    fmula::DataFrames.Formula
    xvars::Vector{Symbol}
    model::Any #DataFrameRegressionModel
    sdf::DataFrame
    wasSuccessful::Bool
    xFromGLM::Vector{Symbol}

    function xGLM(df_data::DataFrame, m::MModel, cfg::OrderedDict)
        this=new()
        this.wasSuccessful=false
        this.xvars=Symbol[]
        this.xFromGLM = vcat([m.y_var, :panid, m.logvar],cfg[:random_demos],cfg[:random_campaigns],cfg[:scoring_vars])
        this.vfactors=unique(setdiff(m.vars, this.xFromGLM))
        for l in 1:30
            this.vfactors=setdiff(this.vfactors,this.xvars)
            this.fmula = genFmula(m.y_var, this.vfactors, m.logvar  )
            try
                #this.model = glm(this.fmula, df_data[  (df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1)  ,vcat(this.vfactors,[m.y_var,m.logvar])], m.dist, LogLink(), offset=log(Array(df_data[(df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1),m.logvar]  + 1)) )
                this.model = glm(this.fmula,  df_data[vcat(m.vars,[m.y_var,m.logvar])], m.dist, LogLink(), offset=log(Array(df_data[m.logvar]  + 1)) )

                this.sdf = DataFrame(vars=vcat([:intercept],this.model.mf.terms.terms), coef=coef(this.model), se=stderr(this.model), zval=coef(this.model)./stderr(this.model) ,pval= ccdf(FDist(1, df_residual(this.model)), abs2(coef(this.model)./stderr(this.model))))

                this.wasSuccessful=true
                break
            catch e
                if isa(e, Base.LinAlg.PosDefException)
                    v=this.vfactors[e.info-1]
                    push!(this.xvars,v)
                    println("!!! Multicollinearity, removing :",v,"~~~",e.info-1, "\n~~~",e)
                else
                    println("....",e)
                    break
                end
            end
        end
        return this
    end

end

type xGLM
    vfactors::Vector{Symbol}
    fmula::DataFrames.Formula
    xvars::Vector{Symbol}
    model::Any #DataFrameRegressionModel
    sdf::DataFrame
    wasSuccessful::Bool
    #xFromGLM::Vector{Symbol}

    function xGLM(dfd::DataFrame, dist::Distribution, y_var::Symbol, logvar::Symbol , vfactors::Array{Symbol} )
        this=new()
        this.wasSuccessful=false
        this.xvars=Symbol[]
        this.vfactors=setdiff(vfactors, [y_var,logvar])
        for l in 1:30
            this.vfactors=setdiff(this.vfactors,this.xvars)
            this.fmula = genFmula(y_var, this.vfactors, logvar  )
            try
                this.model = glm(this.fmula,  dfd[vcat(this.vfactors,[y_var])], dist, LogLink(), offset=log(Array(dfd[logvar]  + 1)) )

                this.sdf = DataFrame(vars=vcat([:intercept],this.model.mf.terms.terms), coef=coef(this.model), se=stderr(this.model), zval=coef(this.model)./stderr(this.model) ,pval= ccdf(FDist(1, df_residual(this.model)), abs2(coef(this.model)./stderr(this.model))))

                this.wasSuccessful=true
                break
            catch e
                if isa(e, Base.LinAlg.PosDefException)
                    v=this.vfactors[e.info-1]
                    push!(this.xvars,v)
                    println("!!! Multicollinearity, removing :",v,"~~~",e.info-1, "\n~~~",e)
                else
                    println("....",e)
                    break
                end
            end
        end
        return this
    end
end

type MOcc <: MModel
    vars::Vector{Symbol}
    y_var::Symbol
    dist::Distribution
    mandatory_vars::Vector{Symbol}
    exclude_vars::Vector{Symbol}
    removed_SingleLevelVars::Vector{Symbol}

    glm1_pvals::xGLM
    glm1_pvals_x::Vector{Symbol}
    glm2_ZnVIF::xGLM
    glm2_ZnVIF_x::Vector{Symbol}
    glm3_PnSigns::xGLM
    glm3_PnSigns_x::Vector{Symbol}
    glm4_PnSignsClean::xGLM
    glm4_PnSignsClean_x::Vector{Symbol}
    glm5::xGLM
    glm5_Z_x::Vector{Symbol}

    Buyer_Pos_P1_is1::Bool

    modelName::AbstractString
    logvar::Symbol

    function MOcc(dfd::DataFrame,cfg::OrderedDict=Dict())
        this=new(); this.modelName="occ"; this.logvar=:trps_pre_p1; this.y_var=:trps_pos_p1; this.dist=Poisson()
        this.mandatory_vars=Symbol[this.y_var]
        this.removed_SingleLevelVars=Symbol[]
        this.glm1_pvals_x=Symbol[]
        this.glm2_ZnVIF_x=Symbol[]
        this.glm3_PnSigns_x=Symbol[]
        this.glm4_PnSignsClean_x=Symbol[]

        this.exclude_vars=Symbol[:dol_per_trip_pos_p1,:nonbuyer_pre_p1,:buyer_pre_p1,:dol_per_trip_pre_p1]
        this.Buyer_Pos_P1_is1=true
        if this.Buyer_Pos_P1_is1  this.exclude_vars = vcat(this.exclude_vars, [:buyer_pos_p1] )  end
        this.vars=setdiff(names(dfd),vcat(this.exclude_vars,[:iso,:whyo,:data_nb_ne_b, :data_b_e_nb ,:pen_reduction,:occ_reduction,:dolocc_reduction]))

        return this
    end

end
#mocc = MOcc(df_data,cfg)

#function rmSingleLevelVar!(dfd::DataFrame,m::MModel)
#    slv=FS_singleLevel(dfd,m.vars)
#    m.removed_SingleLevelVars = slv
#    m.vars=setdiff(m.vars,slv)
#end

function xVIF(g::xGLM)
    include("/home/rmadmin/.julia/v0.4/RegTools/src/diagnostics.jl")
    include("/home/rmadmin/.julia/v0.4/RegTools/src/misc.jl")
    include("/home/rmadmin/.julia/v0.4/RegTools/src/modsel.jl")
    vdf=vif(g.model)
    vdf[:vars] = convert(Array{Symbol}, vdf[:variable])
    g.sdf = join(g.sdf,vdf[[:vars,:vif]], on = :vars, kind=:outer)
end

function modelme(m::MModel,dfd::DataFrame)
    required_vars=vcat([m.y_var,:panid,m.logvar],cfg[:random_demos],cfg[:random_campaigns],cfg[:scoring_vars])
    vars=setdiff(m.vars,required_vars)
    function rmVars(v::Array{Symbol})
        return setdiff(vars,v)
    end

    #SingleValue
    m.removed_SingleLevelVars=FS_singleLevel(dfd,vars)
    vars = rmVars(m.removed_SingleLevelVars)

    #PVals
    m.glm1_pvals = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )
    g1=m.glm1_pvals
    m.glm1_pvals_x=g1.sdf[g1.sdf[:pval].>0.7,:vars]
    vars = rmVars(vcat(m.glm1_pvals_x, g1.xvars ) )

    #Z & Vif
    m.glm2_ZnVIF = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )
    g2=m.glm2_ZnVIF
    xVIF(g2)
    z = g2.sdf[g2.sdf[:zval].<1.96,:vars]
    v = g2.sdf[ !isna(g2.sdf[:vif])&(g2.sdf[:vif].>15),:vars]
    m.glm2_ZnVIF_x =intersect(z,v)
    vars = rmVars(vcat(m.glm2_ZnVIF_x, g2.xvars) )

    # Pvalue & Signs
    m.glm3_PnSigns = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )  # Pvalue & Signs
    g3=m.glm3_PnSigns
    neutralvars = setdiff(vars,vcat(cfg[:negativevars],cfg[:positivevars]))
    neg=intersect(cfg[:negativevars],g3.sdf[ g3.sdf[:coef].<0,:vars])
    pos=intersect(cfg[:positivevars],g3.sdf[ g3.sdf[:coef].>0,:vars])
    varstokeep = intersect(vcat(neutralvars, pos,neg) ,  g3.sdf[ g3.sdf[:pval].<cfg[:pvalue_lvl] ,:vars] )
    m.glm3_PnSigns_x= setdiff(vars, varstokeep )
    vars = rmVars(vcat(m.glm3_PnSigns_x, g3.xvars) )

    # Clean Up Pvalue & Signs
    m.glm4_PnSignsClean = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )
    g4=m.glm4_PnSignsClean
    neutralvars = setdiff(vars,vcat(cfg[:negativevars],cfg[:positivevars]))
    neg=intersect(cfg[:negativevars],g4.sdf[ g4.sdf[:coef].<0,:vars])
    pos=intersect(cfg[:positivevars],g4.sdf[ g4.sdf[:coef].>0,:vars])
    varstokeep = intersect(vcat(neutralvars, pos,neg) ,  g4.sdf[ g4.sdf[:pval].<cfg[:pvalue_lvl] ,:vars] )

    g4.sdf[ g4.sdf[:pval].<cfg[:pvalue_lvl] ,:vars]
#    m.glm4_PnSignsClean_x =
    #m.glm4 = xGLM(dfd, m, cfg)

    #m.vars=vcat(vars,required_vars)
end

#modelme(mocc, df_data[(df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1),:] )

#m=mocc
#dfd=df_data[(df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1),:]

# ==============================================================
# final cleaning of second model based on signs and pvalues
#pvalues_2 <- summary(model_glm_2)$coefficients[,4][-1]
#keep_covariate_2 <- names(pvalues_2[pvalues_2 < pvalue_lvl])
#coeffs_2 <- sign(summary(model_glm_2)$coefficients[,1])
#varstokeep <- c();vars_2 <- c()
#for (i in negativevars) if (i %in% names(coeffs_2[coeffs_2 < 0])) varstokeep <- c(varstokeep,i)
#for (j in positivevars) if (j %in% names(coeffs_2[coeffs_2 > 0])) varstokeep <- c(varstokeep,j)
#for (l in neutralvars) varstokeep <- c(varstokeep,l)
#for (k in varstokeep) if (k %in% keep_covariate_2) vars_2 <- c(vars_2,k)

#modelme(mocc,df_data)
#function glm2df(mdl::DataFrames.DataFrameRegressionModel)
#    return DataFrame(vars=vcat([:intercept],mdl.mf.terms.terms), coef=coef(mdl), se=stderr(mdl), zval=coef(mdl)./stderr(mdl) ,pval= ccdf(FDist(1, df_residual(mdl)), abs2(coef(mdl)./stderr(mdl))))
#end

#g.sdf=glm2df(g.model)

#model_glm_total <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))

#vars = unique(vcat( setdiff(m.vars, m.xFromGLM),[m.y_var, m.logvar]) )
#logvar= m.logvar
#y_var=m.y_var
#runGLM(m.y_var, unique(vcat( setdiff(m.vars, m.xFromGLM),[m.y_var, m.logvar]) ) , m.logvar)
function runGLM(y_var::Symbol, vars::Array{Symbol}, logvar::Symbol=:Missing)
    #*"offset(log("*string(m.logvar)*"+1))"))
    #     vars = setdiff(vars,[:fea_trps_shr_dpp_p1 ,:fea_trps_shr_dpp_p2,:fea_trps_shr_dpp_p3,:fea_trps_shr_dpp_p4
        #,:fea_or_dis_trps_shr_dpp_p1,:fea_or_dis_trps_shr_dpp_p2,:fea_or_dis_trps_shr_dpp_p3,:fea_or_dis_trps_shr_dpp_p4] )
    #cor(df_data[:fea_trps_shr_dpp_p1],df_data[:fea_or_dis_trps_shr_dpp_p1]) #cor(df_data[:fea_trps_shr_dpp_p2],df_data[:fea_or_dis_trps_shr_dpp_p2])
    #cor(df_data[:fea_trps_shr_dpp_p3],df_data[:fea_or_dis_trps_shr_dpp_p3]) #cor(df_data[:fea_trps_shr_dpp_p4],df_data[:fea_or_dis_trps_shr_dpp_p4])

    #vars=setdiff(vars,[:fea_trps_shr_dpp_p3,:fea_or_dis_trps_shr_dpp_p3,:fea_trps_shr_dpp_p4,:fea_or_dis_trps_shr_dpp_p4])
    xvars=Symbol[]
    gmx=nothing
    for l in 1:30
        vars=setdiff(vars,xvars)
        fmla = genFmula(y_var, vars, logvar  )
        try
            if logvar!=:Missing
                #fmla = genFmula(y_var, vars, logvar  )
                gmx = glm(fmla, df_data[vars], Poisson(), LogLink(), offset=log(Array(df_data[logvar]  + 1)) )
            else
                #fmla = genFmula(y_var, vars)
                gmx = glm(fmla, df_data[vars], Poisson() )
            end
            break
        catch e
            if isa(e, Base.LinAlg.PosDefException)
                v=vars[e.info-1]
                push!(xvars,v)
                println("!!! Multicollinearity, removing :",xvars,"~~~",e.info)
            else
                println("....",e)
                break
            end
        end
    end
    return gmx,xvars
end

function modelme(m::MOcc,df_data::DataFrame)
    #Testing
    #vars=unique(vcat( setdiff(m.vars, m.xFromGLM),[m.y_var, m.logvar]) )
    #vars=setdiff(unique(vcat( setdiff(m.vars, m.xFromGLM),[m.y_var, m.logvar]) ) ,  [     :fea_or_dis_trps_shr_dpp_p1,:fea_or_dis_trps_shr_dpp_p2,:fea_or_dis_trps_shr_dpp_p3,:fea_or_dis_trps_shr_dpp_p4 ]  )

    #vars=setdiff(unique(vcat( setdiff(m.vars, m.xFromGLM),[m.y_var, m.logvar]) ) ,  [:fea_trps_shr_dpp_p1 ,:fea_trps_shr_dpp_p2,:fea_trps_shr_dpp_p3,:fea_trps_shr_dpp_p4 ]  )
    #Pkg.clone("https://github.com/joemliang/RegTools.jl")
#    /home/rmadmin/.julia/v0.4/RegTools/src
#    diagnostics.jl  misc.jl  modsel.jl

    #fmla = genFmula(m.y_var, vars, m.logvar  )
    # gmx = glm(fmla, df_data[vars], Poisson(), LogLink(), offset=log(Array(df_data[m.logvar]  + 1)) )
    # df_data[ df_data[:fea_trps_shr_dpp_p1].!=0 ,[:fea_trps_shr_dpp_p1, :fea_or_dis_trps_shr_dpp_p1] ]

    # 1***. drop variables with only one level - cater for buyer_pos_p1
    if m.Buyer_Pos_P1_is1
        m.singleLevel_vars = [c for c in filter(x -> length(unique( df_data[(df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1),x] ))==1, m.vars)]
    else
        m.singleLevel_vars = [c for c in filter(x -> length(unique( df_data[(df_data[:iso].==false),x] ))==1, m.vars)]
    end
    m.vars = setdiff(m.vars,m.singleLevel_vars)
    mdl,xvars = runGLM(m.y_var, unique(vcat( setdiff(m.vars, m.xFromGLM),[m.y_var, m.logvar]) ) , m.logvar)
    m.GLM1_xvars=xvars
    p = glm2df(mdl)
    m.xpvals=vcat(m.xpvals,p[p[:pval].>0.7,:vars])

    include("/home/rmadmin/.julia/v0.4/RegTools/src/diagnostics.jl")
    include("/home/rmadmin/.julia/v0.4/RegTools/src/misc.jl")
    include("/home/rmadmin/.julia/v0.4/RegTools/src/modsel.jl")
    vif(mdl)

end

function xcorr(df::DataFrame,vars::Array{Symbol})
    vx=Symbol[]
    vy=Symbol[]
    vC=Float64[]
    function f(x::Symbol,y::Symbol)
        c=cor(df[x],df[y])
        push!(vx,x)
        push!(vy,y)
        push!(vC,c)
    end
    map(x-> map(y->f(x,y) ,vars)  ,vars)

    for i in 1:length(vC)
        if (vC[i] > 0.8) & (vx[i]!=vy[i])
            println(vx[i]," ~ ",vy[i]," : ",vC[i])
        end
    end
end
#xcorr(df_data,m.vars)

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text: include("SpatialPixelsDataFrame.jl")

ext: jl


text:
using Gnuplot
@gsp "set multiplot title 'Interlocking Tori'"
U = LinRange(-pi, pi, 100) # 50
V = LinRange(-pi, pi, 20)
@gsp :- "unset key" "unset xtics" "unset ytics" "unset ztics"
@gsp :- "set border 0" "set view 60, 30, 1.5, 0.9"
@gsp :- "set style fill transparent solid 0.7"

for loop in 1:2
    if loop == 1
        @gsp :- 1 title="PM3D surface\\nno depth sorting"
        @gsp :- "set origin -0.02,0.0"
        @gsp :- "set size 0.55, 0.9"
        @gsp :- palette(:dense)
        @gsp :- "set colorbox vertical user origin 0.005, 0.15 size 0.02, 0.50"
        @gsp :- "set pm3d scansforward" :-  # scansbackward
    else
        @gsp :- 2 title="PM3D surface\\ndepth sorting"
        @gsp :- "set origin 0.40,0.0"
        @gsp :- "set size 0.55, 0.9"
        @gsp :- palette(:plasma)
        @gsp :- "set colorbox vertical user origin 0.9, 0.15 size 0.02, 0.50"
        @gsp :- "set pm3d depthorder"
    end
    x = [cos(u) + .5 * cos(u) * cos(v)      for u in U, v in V]
    y = [sin(u) + .5 * sin(u) * cos(v)      for u in U, v in V]
    z = [.5 * sin(v)                        for u in U, v in V]
    @gsp :-  x' y' z' "w pm3d"

    x = [1 + cos(u) + .5 * cos(u) * cos(v)  for u in U, v in V]
    y = [.5 * sin(v)                        for u in U, v in V]
    z = [sin(u) + .5 * sin(u) * cos(v)      for u in U, v in V]
    @gsp :- x' y' z' "w pm3d"
end

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
import Base: read, getindex, get, haskey

export JcampdxFile, findfirst_

findfirst_(A, v) = something(findfirst(isequal(v), A), 0)

const JCVAL = Union{AbstractString,Number,Bool,Array,Tuple,Nothing}
const HTSS = Dict{AbstractString,JCVAL}

mutable struct JcampdxFile
  dict::HTSS
end

JcampdxFile() = JcampdxFile(HTSS())

getindex(file::JcampdxFile, key::AbstractString) = file.dict[key]
haskey(file::JcampdxFile, key::AbstractString) = haskey(file.dict, key)
get(file::JcampdxFile, key::AbstractString, default) = get(file.dict, key, default)

function read(file::JcampdxFile, stream::IO, keylist::Vector=String[]; maxEntries=-1)
    finishedReading = true
    currentKey = nothing
    currentIdx = 1
    currentSizes = nothing
    currentEntry = 0
    remainingString = nothing
    tupleReading = false

    skipKeys = ["VisuAcqFrameTime"]

    for line in eachline(stream)

      if line[1] == '#' && line[2] == '#' && line[3] == '$'
        # the last command was not successfully parsed
        finishedReading = true
      end

      if finishedReading
        s = strip(line)
        if s[1] == '#' && s[2] == '#' && s[3] == '$'
          currentEntry += 1
          if maxEntries > 0 && currentEntry == maxEntries
            return file
          end

          i = findfirst_(s, '=')
          key = strip(s[4:i-1])

          # Small HACK
          if in(key, skipKeys)
            return file
          end

          if !isempty(keylist) && findfirst_(keylist,key) == 0
            continue
          end

          val = strip(s[i+1:end])
          if val[1] != '('
            file.dict[key] = val
          else
            if val[2] != ' '
              file.dict[key] = val
            else
              j = findfirst_(val, ')')
              currentSizes = [parse(Int64,s) for s in split(val[2:j-1],",")]
              file.dict[key] = nothing
              currentKey = key
              currentIdx = 1
              finishedReading = false
            end
          end
        end
      else
         if line[1] == '<'
           j = findfirst_(line, '>')
           file.dict[currentKey] = line[2:j-1]
           finishedReading = true
           tupleReading = false
         elseif line[1] == '(' || tupleReading  # skip these for the moment ...
           #finishedReading = true
           tupleReading = true

           if file.dict[currentKey] == nothing
             @debug "Will now allocate memory of size:" currentSizes
             file.dict[currentKey] = Array{Any}(undef, reverse(currentSizes)...)
           end

           totalLine = remainingString == nothing ? line[2:end] : remainingString*line

           parts = split(strip(totalLine), "(")

           for part in parts

             j = findfirst_(part, ')')
             if j != 0
               # Tuple read
               try
                 valsStr = split(strip(part[1:j-1]), ",")
                 vals = Any[]
                 for valStr in valsStr
                   #try
                   #  push!(vals, int(valStr))
                   #catch
                     try
                       push!(vals, parse(Float64,valStr) )
                     catch
                       push!(vals, valStr)
                     end
                   #end
                 end

                 file.dict[currentKey][currentIdx] = vals
                 currentIdx += 1
                 if currentIdx-1 >= length(file.dict[currentKey])
                   finishedReading = true
                   tupleReading = false
                 end
               catch
                 @debug "" currentIdx vals currentSizes size(file.dict[currentKey]) length(keys(file.dict)) currentKey
                 #rethrow()
                 finishedReading = true
                 tupleReading = false
               end

               remainingString = nothing
             else
               remainingString = part
             end
           end

         else
           #print(split(strip(line)," "))
           valsStr = split(strip(line), " ")
           vals = nothing
           try
             vals = parse(Int64,valsStr)
           catch
             try
               vals = parse(Float64,valsStr)
             catch
               vals = valsStr
             end
           end

           if file.dict[currentKey] == nothing
             @debug "Will now allocate memory of size:" currentSizes
             file.dict[currentKey] = Array{eltype(vals)}(undef, reverse(currentSizes)...)
           end

           try
             file.dict[currentKey][currentIdx:currentIdx+length(vals)-1] = vals
             currentIdx += length(vals)
             if currentIdx-1 >= length(file.dict[currentKey])
               finishedReading = true
               tupleReading = false
             end
           catch
             @debug "" currentIdx vals currentSizes size(file.dict[currentKey])
             rethrow()
           end
         end
      end
    end
   file
end

function read(file::JcampdxFile, filename::AbstractString,
              keylist::Vector=String[]; maxEntries=-1)
    open(filename) do f
        read(file, f, keylist, maxEntries=maxEntries)
    end
    file
end


text:
function extend_edge_quadrature_to_surface_quadrature(
    func,
    grad,
    height_dir,
    lo,
    hi,
    x0::V,
    w0,
) where {V<:AbstractVector{T}} where {T}

    extended_func(x) = func(extend(x0, height_dir, x))
    _roots = find_zeros(extended_func, lo, hi)

    num_roots = length(_roots)
    @assert num_roots == 1 || num_roots == 0

    if num_roots == 1
        root = _roots[1]
        p = extend(x0, height_dir, root)
        gradF = grad(p)
        jac = norm(gradF) / (abs(gradF[height_dir]))
        w = w0 * jac
        return p, w
    else
        return nothing, nothing
    end
end

function extend_edge_quadrature_to_surface_quadrature(
    func,
    grad,
    height_dir,
    lo,
    hi,
    x0::M,
    w0::V,
) where {M<:AbstractMatrix{T},V<:AbstractVector} where {T}

    lower_dim, npoints = size(x0)
    nweights = length(w0)
    @assert npoints == nweights

    dim = lower_dim + 1
    quad = TemporaryQuadrature(T, dim)

    for i = 1:npoints
        x0p = view(x0, :, i)
        w0p = w0[i]
        qp, qw = extend_edge_quadrature_to_surface_quadrature(
            func,
            grad,
            height_dir,
            lo,
            hi,
            x0p,
            w0p,
        )
        if !isnothing(qp)
            update!(quad, qp, qw)
        end
    end
    return quad
end

function extend_edge_quadrature_to_surface_quadrature(
    func,
    grad,
    height_dir,
    interval,
    edgequad,
)

    return extend_edge_quadrature_to_surface_quadrature(
        func,
        grad,
        height_dir,
        interval.lo,
        interval.hi,
        edgequad.points,
        edgequad.weights,
    )
end

function cut_surface_quadrature(
    func,
    grad,
    box,
    quad1d,
    recursionlevel,
    maxlevels,
    numsplits,
)

    heightdir = height_direction(grad, box)
    issuitable, gradsign = is_suitable(heightdir, grad, box)
    if !issuitable
        splitboxes = split_box(box, numsplits)
        splitquads = [
            subdivision_surface_quadrature(
                func,
                grad,
                sb,
                quad1d,
                recursionlevel + 1,
                maxlevels,
                numsplits,
            ) for sb in splitboxes
        ]
        return combine_quadratures(splitquads)
    else
        lower(x) = func(extend(x, heightdir, box[heightdir].lo))
        upper(x) = func(extend(x, heightdir, box[heightdir].hi))

        lowersign = sign(gradsign, +1, true, -1)
        uppersign = sign(gradsign, +1, true, +1)

        lowerbox = heightdir == 1 ? box[2] : box[1]

        edgequad = one_dimensional_quadrature(
            [lower, upper],
            [lowersign, uppersign],
            lowerbox,
            quad1d,
        )

        return extend_edge_quadrature_to_surface_quadrature(
            func,
            grad,
            heightdir,
            box[heightdir],
            edgequad,
        )
    end
end

function subdivision_surface_quadrature(
    func,
    grad,
    box,
    quad1d::ReferenceQuadratureRule{T},
    recursionlevel,
    maxlevels,
    numsplits,
) where {T}

    if recursionlevel >= maxlevels
        error("Failed to construct an appropriate surface integration rule after $maxlevels subdivisions")
    else
        s = sign(func, box)
        if s == -1 || s == 0 || s == +1
            if s != 0
                return TemporaryQuadrature(T, 2)
            else
                return cut_surface_quadrature(
                    func,
                    grad,
                    box,
                    quad1d,
                    recursionlevel,
                    maxlevels,
                    numsplits,
                )
            end
        else
            splitboxes = split_box(box, numsplits)
            splitquads = [
                subdivision_surface_quadrature(
                    func,
                    grad,
                    sb,
                    quad1d,
                    recursionlevel + 1,
                    maxlevels,
                    numsplits,
                ) for sb in splitboxes
            ]
            return combine_quadratures(splitquads)
        end
    end
end

function surface_quadrature(func, grad, xL, xR, numqp; maxlevels = 5, numsplits = 2)

    box = IntervalBox(xL, xR)
    quad1d = ReferenceQuadratureRule(numqp)

    tempquad =
        subdivision_surface_quadrature(func, grad, box, quad1d, 0, maxlevels, numsplits)
    return QuadratureRule(tempquad)
end

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
## This kinda follows the algorithm laid out in
## http://www.dgp.toronto.edu/people/mooncake/papers/SIGGRAPH2001_Tupper.pdf
## the paper is much more thorough

## Helper function to find initial partition into squares
function fill_bottom(W,H)
    k = min(floor(Integer,log2(W)), floor(Integer, log2(H)))
    rects = Any[] # Array{Int,1}[] #

    ## fill in
    w,h = 0, 2^k

    while (W - w) >= 2^k
        push!(rects, [w+0, w+2^k, 0, 2^k])
        w = w + 2^k
    end

    j= k - 1
    while j > 0
        if (W - w) >= 2^j
            for i in 1:2^(k-j)
                push!(rects, [w+0, w+2^j, (i-1)*2^j, i*2^j])
            end
            w = w + 2^j
        end
        j = j - 1
    end
    rects, h
end

## Helper function to break a rectangle into squares of size 2^j x 2^j
function break_into_squares(W, H)
    rects = Any[] # Array{Int,1}[] #
    offset = 0

    while H - offset > 0
        nrects, h = fill_bottom(W, H-offset)
        nrects = map(u -> [0,0,offset, offset] + u, nrects)
        append!(rects, nrects)
        offset = offset + h
    end
    rects
end

"""
Main algorithm of Tupper.

Break region into square regions
for each square region, subdivide into 4 regions. For each check if there is no solution, if so add to white. Otherwise, refine the region by repeating until we can't subdivide

At that point, check each pixel-by-pixel region for possible values.

Return red, black and white vectors of Regions.
"""
function GRAPH(r, L, R, B, T, W, H)
    rects = break_into_squares(W, H)

    k = min(floor(Integer,log2(W)), floor(Integer,log2(H))) # largest square is size 2^k x 2^k

    reds = [Region(OInterval(u[1], u[2]), OInterval(u[3], u[4])) for u in rects]
    sizehint!(reds, W)

    red = Region[]         # 1-pixel red, can't decide via check_continuity
    black = Region[]
    white = Region[]

    while (k >= 0) & (length(reds) > 0)
        reds = RefinePixels(r, reds, L, R, B, T, W, H, black, white, red)
        k = k - 1
    end
    red, black, white
end

## Refine the region
function RefinePixels(r, U_k, L, R, B, T, W, H, black, white, red)
    ## Uk_1 a refinement of U_k which hold red regions
    Uk_1 = Region[]
    for u in U_k
        out = compute(r, u,  L, R, B, T, W, H)
        if out == TRUE
            push!(black, u)
        elseif out == FALSE
            push!(white, u)
        else
            ## subdivide and call again if possible,
            x = u.x.val; y = u.y.val
            dx, dy = diam(x), diam(y)
            if (dx > 1) & (dy > 1)
                hx = div(dx,2); hy = div(dy,2)
                for i in 0:1, j in 0:1
                    uij = Region(OInterval(x.lo + i*hx, x.lo + (i+1)*hx),
                                      OInterval(y.lo + j*hy, y.lo + (j+1)*hy))
                    push!(Uk_1, uij)
                end
            else
                ## these are red with diameter 1, could be white or black
                val = check_continuity(r, u, L, R, B, T, W, H)
                if val == TRUE
                    push!(black, u)
                elseif val == FALSE
                    push!(white, u)
                else
                    push!(red, u)
                end
            end
        end
    end
    Uk_1
end

## for 1-pixel squares, check NaN and continuity
## Return TRUE (Black), FALSE (white) or MAYBE (red)
function check_continuity(r::Pred, u, L, R, B, T, W, H)

    fxy = compute_fxy(r, u,  L, R, B, T, W, H)

    ## check for NaN
    if isempty(fxy.val)
        return(FALSE)
    end
    if (fxy.def == FALSE) || (fxy.def == MAYBE)
        return(FALSE)
    end

    ## now check continuity,
    val = FALSE
    if (fxy.cont == TRUE) && ((r.op === ==) || (r.op === <=) || (r.op === >=))
        ## use intermediate value theorem here
        val = val | cross_zero(r, u, L, R, B, T, W, H)

    end

    ## Now check for inequalities
    ineqs = [<, <= , ≤, ≶, ≷, ≥, >=, >]
    if (fxy.def == TRUE) && any([r.op === op for op in ineqs])
        ## just check points
        val = val | check_inequality(r, u, L, R, B, T, W, H)
    end

    ## What to do if fxy.cont !== TRUE...
    default = MAYBE
    if val == MAYBE
        return(default)
    else
        return(val)
    end
end

## Return TRUE, FALSE or MAYBE for predicates
function check_continuity(rs::Preds, u, L, R, B, T, W, H)
    vals = map(r -> check_continuity(r, u, L, R, B, T, W, H), rs.ps)
    val = popfirst!(vals)
    for i in 1:length(rs.ops)
        val = rs.ops[i](val, vals[i])
    end
    val
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

NamespaceA.eval(quote

import ..NamespaceC

FlatBuffers.@with_kw mutable struct SecondTableInA
    refer_to_c::Union{NamespaceC.TableInC, Nothing} = nothing
end
FlatBuffers.@ALIGN(SecondTableInA, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:SecondTableInA} = [
    0x00000004
]

SecondTableInA(buf::AbstractVector{UInt8}) = FlatBuffers.read(SecondTableInA, buf)
SecondTableInA(io::IO) = FlatBuffers.deserialize(io, SecondTableInA)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
using WordErrorRate
using Base.Test

include("test.jl")

x = readtestdata()
y = readsclitepralign()
ncorrect = 0
for (key, val) in y
	ref, hyp = x[key]
	w = WER(ref, hyp)
	alignment = chomp(pralign(String, w))
	ncorrect += alignment == val
end
println(@sprintf("Ran tests: %d/%d = %4.1f%% correct", ncorrect, length(y), 100ncorrect/length(y)))
@test ncorrect == length(y)


text:
#' Boosting the Hodrick-Prescott Filter
#'
#' Coded by Mei Ziwei
#' Documented by Shi Zhentao
#'
#' All in one function of conducting the boosted HP-filter.
#' tesing example is at the end of the script
##############################################################################
#' Parameters:
#'
#' x - a raw time series to be filtered.
#' lam - turning parameter, default value is 1600,
#'   as recommended by Hodrick and Prescott (1997) for quarterly data.
#' iter - logical, True (default) to conduct the boosted HP filter.
#'   False does not iterated, which is exactly the original HP filter.
#' stopping - stopping criterion.  "BIC" (default), or  "adf", or  "nonstop"  means keeping
#'    iteration until the maximum number of iteration, specified by Max_Iter is reached.
#' sig_p - a threshold of the p-value for the ADF test, with default value 0.050.
#'    Only effective when stopping = "adf".
#' Max_Iter - maximal number of iterations. The default is 100.
##############################################################################
#' The function returns a dictionary containing the following items:
#'
#' cycle - The cyclical component in the final iteration.
#' trend - The trend component in the final iteration.
#' trend_hist - The estimated trend in each iteration.
#' iter_num - The total number of iterations when it stops.
#' BIC_hist - The path of the BIC up to the final iterations.
#' adf_p_hist - The path of the ADF test p-value up to the final iteration.
##############################################################################
#' Details:
#'
#' This is the main function of implementing the boosted HP filter (Phillisp and
#' Shi, 2019). The arguments accommendate the orginal HP filter (iter =
#' False), the boosted HP filter with the BIC stopping criterion (stopping = "BIC"),
#' or ADF test stopping criterion (stopping = "adf"),
#' or keep going until the maximum number of iterations is reached (stopping = "nonstop").
#'
#' Either the original HP filter or the bHP filter requires lambda to
#' control the strength of the weak learner for in-sample fitting.
#' The default is lambda = 1600,
#' which is recommended by Hodrick and Prescott (1997) for quarterly data.
#' lambda should be adjusted for different frequencies.
#' For example, lambda = 129600 for monthly data and
#' lambda = 100 or 6.25 for annual data.
#'
#' See the vignette with a brief introduction of the idea of bHP.
#'
##############################################################################
#' References:
#'
#' Phillips, Peter CB, and Zhentao Shi. "Boosting: Why you can use the hp
#' filter." arXiv: 1905.00175, Cowles Foundation Discussion Paper No.2192,
#' (2019).
#'
##############################################################################

using LinearAlgebra
using Statistics
using HypothesisTests
function BoostedHP(x,λ=1600,iter=true,stopping="BIC",sig_p=0.050,Max_Iter=100)

    ## Examine data format
    ## generating trend operator matrix "S"
    raw_x = x;
    n = length(x);

    I_n = Matrix(I,n,n);
    D_temp = [zeros(1,n);Matrix(I,n-1,n)];
    D_temp = (I_n-D_temp)*(I_n-D_temp);
    D = D_temp[3:n, :]';

    S = inv(I_n+λ*(D*D'));
    mS = I_n - S ;

    # the Simple HP Filter
    if ~iter

        # get the trend and cycle
        x_f = S*x;
        x_c = x - x_f;

        # Results
        result = Dict("cycle" => x_c, "trend_hist" => x_f,
                   "stopping" => "nonstop", "trend" => x - x_c, "raw_data" => raw_x);
    end

    # the Boosted HP-filter
    if iter

        ## ADF test as the stopping criterion
        if stopping == "adf"

            r = 1;
            stationary = false;
            x_c = x;

            x_f = zeros(n,Max_Iter);
            adf_p = zeros(Max_Iter);

            while (r <= Max_Iter) & ~stationary

                x_c = mS * x_c; # update
                x_f[:,r] = x - x_c;

                adf_p_r = pvalue(ADFTest(x_c,:trend,Int(floor((n-1)^(1/3)))));
                # x_c is the residual after the mean and linear trend being removed by HP filter
                # we use the critical value for the ADF distribution with
                # the intercept and linear trend specification

                adf_p[r] = adf_p_r;
                stationary = (adf_p_r <= sig_p);
                # Truncate the storage matrix and vectors

                if stationary
                    R = r;
                    x_f = x_f[:,1:R];
                    adf_p = adf_p[1:R];
                    break
                end

                r += 1
            end  # end the while loop

            if r > Max_Iter
                R = Max_Iter;
                warning("The number of iterations exceeds Max_Iter.
                The residual cycle remains non-stationary.")

            end
            result = Dict("cycle" => x_c, "trend_hist" => x_f,  "stopping" => stopping,
                 "signif_p" => sig_p, "adf_p_hist" => adf_p, "iter_num" => R,
                 "trend" => x - x_c, "raw_data" => raw_x);

        else # either

            #assignment
            r = 0;
            x_c_r = x;
            x_f = zeros(n,Max_Iter);
            IC = zeros(Max_Iter);
            # IC_decrease = true;

            I_S_0 = I_n - S;
            c_HP = I_S_0 * x;
            I_S_r = I_S_0;

            while r < Max_Iter

                r = r + 1;

                x_c_r = I_S_r * x;# this is the cyclical component after m iterations
                x_f[:,r] = x - x_c_r;
                B_r = I_n - I_S_r;
                IC[r] = var(x_c_r)/var(c_HP) + log(n)/(n-sum(diag(S)))*sum(diag(B_r));

                I_S_r = I_S_0 * I_S_r; # update for the next round

                if (r>=2) & (stopping == "BIC")
                    if IC[r-1] < IC[r]
                        break
                    end
                end

            end # end the shile loop

            # Final assignment
            R = r - 1;
            x_f = x_f[:, 1:R]
            x_c = x - x_f[:,R]

            if stopping == "BIC"
                # save the path of BIC till iter+1 times to keep the "turning point" of BIC history.
                result = Dict("cycle" => x_c, "trend_hist" => x_f,  "stopping" => stopping,
                       "BIC_hist" => IC[1:(R+1)], "iter_num" => R, "trend" =>  x - x_c, "raw_data" => raw_x);
            end

            if stopping == "nonstop"

                result = Dict("cycle" => x_c, "trend_hist" => x_f,  "stopping" => stopping,
                       "BIC_hist" => IC, "iter_num" => Max_Iter-1, "trend" => x - x_c, "raw_data" => raw_x);
            end

        end

    end  # end the Boosted HP

    return result

end  ## end function
###############################################################################
## Examples

using CSV
using DelimitedFiles
IRE = readdlm("IRE.csv", ',')[:]; # laod the data "IRE"
λ = 100;  # tuning the parameter for annumal data
#
#' # raw HP filter
bx_HP = BoostedHP(IRE, λ, false);
# bx_HP_cycle = bx_HP["cycle"]  # The cyclical component
# bx_HP_trend = bx_HP["trend"]  # The trend component
#'
#' # by BIC
bx_BIC = BoostedHP(IRE, λ, true,"BIC");
#'
#' # by ADF
#' # CAVEAT: Results may be slightly different from other languages
#' #    due to different critical values of ADF test.
bx_ADF = BoostedHP(IRE, λ, true,"adf");
#'
#' # If stopping = "nonstop",
#' # Iterated HP filter until Max_Iter and keep the path of BIC.
bx_nonstop = BoostedHP(IRE, λ, true,"nonstop");


text:
module LocalDescentTests

import ..Algopt # ensure we are using the correct Algopt

using Test

using Algopt.LocalDescent:
search,
find_bracket, zoom_bracket!,
LineSearch, StrongBacktracking,
DescentPointK, MutBracket,
BracketConditions, StrongWolfeConditions

using Algopt.TestFunctions: Rosenbrock

quad0 = x->x'x
∇quad0 = x->2x

quad3 = x->(x - [3, 3])' * (x - [3, 3])
∇quad3 = x->2(x - [3, 3])

ros = Rosenbrock(a = 1, b = 5)

@testset "Algopt.LocalDescent.search(::LineSearch)" begin
    ls = LineSearch()
    @test [0, 0] ≈ search(ls, quad0, [10, 10], [-1, -1])
    @test [3, 3] ≈ search(ls, quad3, [10, 10], [-1, -1])

    # Line search of α on 2D Rosenbrock function
    ls = LineSearch()
    @test [1, 1] ≈ search(ls, ros.f, [10, 10], [-1, -1])
    @test [1, 1] ≈ search(ls, ros.f, [-2, -2], [1, 1])
end

@testset "Algopt.LocalDescent.search(::StrongBacktracking)" begin
    sb = StrongBacktracking()
    @test [0, 0] ≈ search(sb, quad0, ∇quad0, [10, 10], [-1, -1])
    @test [3, 3] ≈ search(sb, quad3, ∇quad3, [10, 10], [-1, -1])

    # Line search of α on 2D Rosenbrock function
    sb = StrongBacktracking()
    @test [-2, -2] ≈ search(sb, ros.f, ros.∇f, [10, 10], [-1, -1])
    @test [1, 1] ≈ search(sb, ros.f, ros.∇f, [-2, -2], [1, 1])
end

@testset "Algopt.LocalDescent.zoom_bracket!" begin
    params = StrongBacktracking()
    point_k = DescentPointK(quad0, ∇quad0, [10,10], [-1,-1])
    wolfe_cond = StrongWolfeConditions(params, quad0, ∇quad0, point_k)
    bracket_cond = BracketConditions(wolfe_cond, quad0, ∇quad0, point_k)

    bracket = MutBracket(0, 20)
    α = zoom_bracket!(bracket, params, wolfe_cond, bracket_cond.third)
    @test α ≈ 10

    bracket = MutBracket(0, 10)
    α = zoom_bracket!(bracket, params, wolfe_cond, bracket_cond.third)
    @test α ≈ 10

    bracket = MutBracket(0, 12)
    α = zoom_bracket!(bracket, params, wolfe_cond, bracket_cond.third)
    @test α > 9
    @test α < 11

    bracket = MutBracket(10, 20)
    α = zoom_bracket!(bracket, params, wolfe_cond, bracket_cond.third)
    @test α > 10
    @test α < 11

    bracket = MutBracket(9, 20)
    α = zoom_bracket!(bracket, params, wolfe_cond, bracket_cond.third)
    @test α > 10
    @test α < 11
end

@testset "Algopt.LocalDescent.find_bracket" begin
    params = StrongBacktracking()
    point_k = DescentPointK(quad0, ∇quad0, [10,10], [-1,-1])
    wolfe_cond = StrongWolfeConditions(params, quad0, ∇quad0, point_k)
    bracket_cond = BracketConditions(wolfe_cond, quad0, ∇quad0, point_k)
    bracket = find_bracket(params, bracket_cond)
    @test bracket.left < 10
    @test bracket.right > 10
end

@testset "Algopt.LocalDescent.DescentPointK" begin
    point_k = DescentPointK(quad0, ∇quad0, [10,10], [-1,-1])
    @test point_k.f == 200
    @test point_k.∇df == -40  # [20, 20]' * [-1, -1]
    @test point_k.x == [10, 10]
    @test point_k.d == [-1, -1]
end

@testset "Algopt.LocalDescent.StrongWolfeConditions" begin
    params = StrongBacktracking()
    point_k = DescentPointK(quad0, ∇quad0, [10,10], [-1,-1])
    wolfe = StrongWolfeConditions(params, quad0, ∇quad0, point_k)

    @test wolfe.first(0)
    @test wolfe.first(1)
    @test wolfe.first(10)

    @test !wolfe.second(1)
    @test !wolfe.second(8)
    @test wolfe.second(9)
    @test wolfe.second(10)
    @test wolfe.second(11)
    @test !wolfe.second(12)
    @test !wolfe.second(20)
end

@testset "Algopt.LocalDescent.BracketConditions" begin
    params = StrongBacktracking(; β = 1e-1)
    point_k = DescentPointK(quad0, ∇quad0, [10,10], [-1,-1])
    wolfe = StrongWolfeConditions(params, quad0, ∇quad0, point_k)
    bracket = BracketConditions(wolfe, quad0, ∇quad0, point_k)

    @test bracket.first(0)
    @test !bracket.first(1)
    @test !bracket.first(10)
    @test !bracket.first(19)
    @test bracket.first(20)
    @test bracket.first(21)

    @test !bracket.second(0)
    @test !bracket.second(10)
    @test !bracket.second(18)
    @test bracket.second(18.000001)
    @test bracket.second(20)

    @test !bracket.third(0)
    @test !bracket.third(9)
    @test bracket.third(10)
    @test bracket.third(11)
    @test bracket.third(20)
end

end # module


text:
## BUILDING CHAINS A FROM HYPERPARAMETERS + INPUT/OUTPUT SHAPE

# We introduce chain builders as a way of exposing neural network
# hyperparameters (describing, architecture, dropout rates, etc) to
# the MLJ user. These parameters generally exclude the input/output
# dimensions/shapes, as the MLJ fit methods will determine these from
# the training data. A `Builder` object stores the parameters and an
# associated `fit` method builds a corresponding chain given the
# input/output dimensions/shape.

# Below n or (n1, n2) etc refers to network inputs, while m or (m1,
# m2) etc refers to outputs.

abstract type Builder <: MLJModelInterface.MLJType end

"""
    Linear(; σ=Flux.relu, rng=Random.GLOBAL_RNG)

MLJFlux builder that constructs a fully connected two layer network
with activation function `σ`. The number of input and output nodes is
determined from the data. The bias and coefficients are initialized
using `Flux.glorot_uniform(rng)`. If `rng` is an integer, it is
instead used as the seed for a `MersenneTwister`.

"""
mutable struct Linear <: Builder
    σ
end
Linear(; σ=Flux.relu) = Linear(σ)
build(builder::Linear, rng, n::Integer, m::Integer) =
    Flux.Chain(Flux.Dense(n, m, builder.σ, init=Flux.glorot_uniform(rng)))

"""
    Short(; n_hidden=0, dropout=0.5, σ=Flux.sigmoid, rng=GLOBAL_RNG)

MLJFlux builder that constructs a full-connected three-layer network
using `n_hidden` nodes in the hidden layer and the specified `dropout`
(defaulting to 0.5). An activation function `σ` is applied between the
hidden and final layers. If `n_hidden=0` (the default) then `n_hidden`
is the geometric mean of the number of input and output nodes.  The
number of input and output nodes is determined from the data.

The each layer is initialized using `Flux.glorot_uniform(rng)`. If
`rng` is an integer, it is instead used as the seed for a
`MersenneTwister`.

"""
mutable struct Short <: Builder
    n_hidden::Int     # if zero use geometric mean of input/output
    dropout::Float64
    σ
end
Short(; n_hidden=0, dropout=0.5, σ=Flux.sigmoid) = Short(n_hidden, dropout, σ)
function build(builder::Short, rng, n, m)
    n_hidden =
        builder.n_hidden == 0 ? round(Int, sqrt(n*m)) : builder.n_hidden
    init=Flux.glorot_uniform(rng)
    Flux.Chain(
        Flux.Dense(n, n_hidden, builder.σ, init=init),
        # TODO: fix next after https://github.com/FluxML/Flux.jl/issues/1617
        Flux.Dropout(builder.dropout),
        Flux.Dense(n_hidden, m, init=init))
end

"""
    MLP(; hidden=(10,), σ=Flux.sigmoid, rng=GLOBAL_RNG)

MLJFlux builder that constructs a Multi-layer perceptron network. The ith element of `hidden` represents the number of neurons in the ith hidden layer. An activation function `σ` is applied between each layer.

The each layer is initialized using `Flux.glorot_uniform(rng)`. If
`rng` is an integer, it is instead used as the seed for a
`MersenneTwister`.

"""
mutable struct MLP{N} <: MLJFlux.Builder
    hidden::NTuple{N, Int}  # count first and last layer
    σ
end
MLP(; hidden=(10,), σ=Flux.sigmoid) = MLP(hidden, σ)
function MLJFlux.build(mlp::MLP, rng, n_in, n_out)
    init=Flux.glorot_uniform(rng)
    hidden = []
    for i ∈ 2:length(mlp.hidden)
        push!(hidden, Flux.Dense(mlp.hidden[i-1], mlp.hidden[i], mlp.σ, init=init))
    end

    return Flux.Chain(
        Flux.Dense(n_in, mlp.hidden[1], mlp.σ, init=init),
        hidden...,
        Flux.Dense(mlp.hidden[end], n_out, init=init)
    )
end

struct GenericBuilder{F} <: Builder
    apply::F
end

"""
    @builder neural_net

Creates a builder for `neural_net`. The variables `rng`, `n_in`, `n_out` and
`n_channels` can be used to create builders for any random number generator `rng`,
input and output sizes `n_in` and `n_out` and number of input channels `n_channels`.

# Examples
```jldoctest
julia> import MLJFlux: @builder;

julia> nn = NeuralNetworkRegressor(builder = @builder(Chain(Dense(n_in, 64, relu),
                                                            Dense(64, 32, relu),
                                                            Dense(32, n_out))));

julia> conv_builder = @builder begin
           front = Chain(Conv((3, 3), n_channels => 16), Flux.flatten)
           d = Flux.outputsize(front, (n_in..., n_channels, 1)) |> first
           Chain(front, Dense(d, n_out));
       end

julia> conv_nn = NeuralNetworkRegressor(builder = conv_builder);
```
"""
macro builder(ex)
    esc(quote
        MLJFlux.GenericBuilder((rng, n_in, n_out, n_channels) -> $ex)
    end)
end

build(b::GenericBuilder, rng, n_in, n_out, n_channels = 1) = b.apply(rng, n_in, n_out, n_channels)


text:
const FLOATING_POINT_REGEX = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?"
const METERS_PER_FOOT = 0.3048

type NGSIMRoadway
    name::Symbol
    boundaries::Vector{Vector{VecE2}}
    centerlines::Vector{Vector{CurvePt}}
end
immutable RoadwayInputParams
    filepath_boundaries::String
    filepath_centerlines::String
end

function read_boundaries(io::IO)

    lines = readlines(io)
    for (i,line) in enumerate(lines)
        lines[i] = strip(line)
    end
    @assert(lines[1] == "BOUNDARIES")

    n_boundaries = parse(Int, lines[2])
    @assert(n_boundaries ≥ 0)

    retval = Array(Vector{VecE2}, n_boundaries)
    line_index = 2
    for i in 1 : n_boundaries

        @assert(lines[line_index+=1] == @sprintf("BOUNDARY %d", i))
        npts = parse(Int, lines[line_index+=1])
        line = Array(VecE2, npts)
        for j in 1 : npts
            matches = matchall(FLOATING_POINT_REGEX, lines[line_index+=1])
            x = parse(Float64, matches[1]) * METERS_PER_FOOT # convert to meters
            y = parse(Float64, matches[2]) * METERS_PER_FOOT
            line[j] = VecE2(x,y)
        end
        retval[i] = line
    end

    retval
end
function read_centerlines(io::IO)

    lines = readlines(io)
    for (i,line) in enumerate(lines)
        lines[i] = strip(line)
    end
    @assert(lines[1] == "CENTERLINES")

    n_centerlines = parse(Int, lines[2])
    @assert(n_centerlines ≥ 0)

    line_index = 2
    retval = Dict{AbstractString, Vector{CurvePt}}()
    for i in 1 : n_centerlines
        @assert(lines[line_index+=1] == "CENTERLINE")
        name = lines[line_index+=1]
        npts = parse(Int, lines[line_index+=1])
        line = Array(VecE2, npts)
        for j in 1 : npts
            matches = matchall(FLOATING_POINT_REGEX, lines[line_index+=1])
            x = parse(Float64, matches[1]) * METERS_PER_FOOT # convert to meters
            y = parse(Float64, matches[2]) * METERS_PER_FOOT # convert to meters
            line[j] = VecE2(x,y)
        end

        # post-process to extract heading and distance along lane
        centerline = Array(CurvePt, npts)
        let
            θ = atan2(line[2]-line[1])
            centerline[1] = CurvePt(VecSE2(line[1],θ), 0.0)
            for i in 2 : npts-1
                s = centerline[i-1].s + hypot(line[i] - line[i-1])
                θ = (atan2(line[i]-line[i-1]) + atan2(line[i+1]-line[i]))/2 # mean angle
                centerline[i] = CurvePt(VecSE2(line[i],θ), s)
            end
            s = centerline[npts-1].s + hypot(line[npts] - line[npts-1])
            θ = atan2(line[npts]-line[npts-1])
            centerline[npts] = CurvePt(VecSE2(line[npts],θ), s)
        end

        retval[name] = centerline

    end

    retval
end
function read_roadway(input_params::RoadwayInputParams)
    boundaries = open(read_boundaries, input_params.filepath_boundaries)
    centerlines = open(read_centerlines, input_params.filepath_centerlines)

    name = symbol(splitext(splitdir(input_params.filepath_boundaries)[2])[1])
    NGSIMRoadway(name, boundaries, collect(values(centerlines)))
end

function write_lwpolyline(io::IO, pts::Vector{CurvePt}, handle_int::Int, is_closed::Bool=false)
    N = length(pts)
    println(io, "  0")
    println(io, "LWPOLYLINE")
    println(io, "  5") # handle (increases)
    @printf(io, "B%d\n", handle_int)
    println(io, "100") # subclass marker
    println(io, "AcDbEntity")
    println(io, "  8") # layer name
    println(io, "150")
    println(io, "  6") # linetype name
    println(io, "ByLayer")
    println(io, " 62") # color number
    println(io, "  256")
    println(io, "370") # lineweight enum
    println(io, "   -1")
    println(io, "100") # subclass marker
    println(io, "AcDbPolyline")
    println(io, " 90") # number of vertices
    @printf(io, "   %d\n", N)
    println(io, " 70") # 0 is default, 1 is closed
    @printf(io, "    %d\n", is_closed ? 1 : 0)
    println(io, " 43") # 0 is constant width
    println(io, "0")

    for i in 1 : N
        println(io, " 10")
        @printf(io, "%.3f\n", pts[i].pos.x)
        println(io, " 20")
        @printf(io, "%.3f\n", pts[i].pos.y)
    end
end

function convert_curves_feet_to_meters!(roadway::Roadway)
    for seg in roadway.segments
        for lane in seg.lanes
            for (i,P) in enumerate(lane.curve)
                lane.curve[i] = CurvePt(
                        VecSE2(P.pos.x*METERS_PER_FOOT, P.pos.y*METERS_PER_FOOT, P.pos.θ),
                        P.s*METERS_PER_FOOT, P.k/METERS_PER_FOOT, P.kd/METERS_PER_FOOT)
            end
        end
    end
    roadway
end

function write_dxf(io::IO, roadway::NGSIMRoadway)

    lines = open(readlines, Pkg.dir("NGSIM", "data", "template.dxf"))
    i = findfirst(lines, "ENTITIES\n")
    i != 0 || error("ENTITIES section not found")

    # write out header
    for j in 1 : i
        print(io, lines[j])
    end

    # write out the lanes
    for (handle_int, lane) in enumerate(roadway.centerlines)
        write_lwpolyline(io, lane, handle_int)
    end

    # write out tail
    for j in i+1 : length(lines)
        print(io, lines[j])
    end
end
function write_roadways_to_dxf()

    roadway_input_80 = RoadwayInputParams(Pkg.dir("NGSIM", "data", "boundaries80.txt"),
                                          Pkg.dir("NGSIM", "data", "centerlines80.txt"))
    roadway_input_101 = RoadwayInputParams(Pkg.dir("NGSIM", "data", "boundaries101.txt"),
                                           Pkg.dir("NGSIM", "data", "centerlines101.txt"))

    ngsimroadway_80 = read_roadway(roadway_input_80)
    ngsimroadway_101 = read_roadway(roadway_input_101)

    open(io->write_dxf(io, ROADWAY_80), Pkg.dir("NGSIM", "data", "ngsim_80.dxf"), "w")
    open(io->write_dxf(io, ROADWAY_101), Pkg.dir("NGSIM", "data", "ngsim_101.dxf"), "w")
end
function write_roadways_from_dxf()

    roadway_80 = open(io->read_dxf(io, Roadway, dist_threshold_lane_connect=2.0), Pkg.dir("NGSIM", "data", "ngsim_80.dxf"), "r")
    roadway_101 = open(io->read_dxf(io, Roadway, dist_threshold_lane_connect=2.0), Pkg.dir("NGSIM", "data", "ngsim_101.dxf"), "r")

    # also converts to meters
    convert_curves_feet_to_meters!(roadway_80)
    convert_curves_feet_to_meters!(roadway_101)

    open(io->write(io, roadway_80), Pkg.dir("NGSIM", "data", "ngsim_80.txt"), "w")
    open(io->write(io, roadway_101), Pkg.dir("NGSIM", "data", "ngsim_101.txt"), "w")
end

const ROADWAY_80 = open(io->read(io, Roadway), Pkg.dir("NGSIM", "data", "ngsim_80.txt"), "r")
const ROADWAY_101 = open(io->read(io, Roadway), Pkg.dir("NGSIM", "data", "ngsim_101.txt"), "r")


text:
using MacroTools

export Parameter
export @gen_convert_parameter, convert_parameter

@doc doc"""
    @gen_convert_parameter(model_typename, (keyname, size_fn, default)...)

Generates `convert_parameter(model::model_typename, param::Parameter)`.

# Example

``` julia
@gen_convert_parameter(A, ("A", numbondtypes, 1.0), ("B", 1, 1))
```

generates a function equivalent to the following:

``` julia
doc\"""
    convert_parameter(model::A, param::Parameter)

# Keynames

- "A": a vector with `numbondtypes(model)` elements (default: 1)
- "B": a scalar (default: 1.0)

\"""
function convert_parameter(model::A, param::Parameter)
    ## if `size_fn` is a `Function`,
    ## result is a vector whose size is `size_fn(model)`.
    ## `param["A"]` can take a scalar or a vector.
    a = get(param, "A", 1.0)
    as = zeros(Float64, numbondtypes(model))
    as .= a

    ## otherwise,
    ## result is a scalar.
    b = convert(Int, get(param, "B", 1))

    return as, b
end
```
"""
macro gen_convert_parameter(model_typename, args...)
    body = Expr(:block)
    str_model = sprint(Base.show_unquoted,model_typename)

    document = "    convert_parameter(model::$(str_model), param::Parameter)\n\n"
    document *= "# Keynames:\n"

    syms = Symbol[]
    for arg in args
        @capture(arg, (name_String, sz_, default_))
        eltyp = typeof(default)
        sym = gensym()
        if isa(eval(sz), Function)
            push!(body.args,
                  esc(:( $(sym) = get(param, $name, $default) ))
                 )
            sym2 = gensym()
            push!(syms, sym2)
            push!(body.args,
                  esc(:( $(sym2) = zeros($eltyp, $sz(model))))
                 )
            push!(body.args,
                  esc(:( $(sym2) .= $(sym) ))
                 )
            document *= "- \"$name\": a vector with `$(sz)(model)` elements (default: $default).\n"
        else
            push!(body.args,
                  esc(:( $(sym) = convert($eltyp, get(param, $name, $default))))
                 )
            push!(syms, sym)
            document *= "- \"$name\": a scalar (default: $default).\n"
        end
    end
    ret = Tuple( s for s in syms )
    push!(body.args,
          esc(:(return tuple($(syms...))  ))
         )

    # res_fn = :(function $(esc(:convert_parameter))($(esc(:(model::$model_typename))),
    res_fn = :(function $(esc(:convert_parameter))($(esc(:model))::$(esc(model_typename)),
                                                   $(esc(:param))::$(esc(:Parameter))
                                                  )
                   $body
               end)

    :( @doc $document $res_fn )
end

@doc doc"""
    convert_parameter(model, param)

Generates arguments of updater and estimator.

# Example
``` julia-repl
julia> model = Ising(chain_lattice(4));

julia> p = convert_parameter(model, Parameter("J"=>1.0))
(1.0, [1.0, 1.0]) # T and Js

julia> p = convert_parameter(model, Parameter("J"=>[1.5, 0.5]))
(1.0, [1.5, 0.5]) # J can take a vector whose size is `numbondtypes(model)`

julia> model.spins
4-element Array{Int64,1}:
  1
  1
  1
 -1

julia> local_update!(model, p...);

julia> model.spins
4-element Array{Int64,1}:
  1
  1
 -1
 -1
```
"""
function convert_parameter end


text:
root_dir = "../."
using Pkg; Pkg.activate(root_dir)
using CSV
using DataFrames
using Dates: Date, DateFormat, today, format
import DotEnv
using Redis
import Tables.namedtupleiterator

cfg = DotEnv.config(path="$root_dir/.env")

curDate = format(today(), DateFormat("yymmdd"))
#srcFile = "$(cfg["OUTPUT_DIR"])/$curDate.trade_info"
srcFile = "$(cfg["OUTPUT_DIR"])/211219.trade_info"
dstFile = "$(cfg["OUTPUT_DIR"])/$curDate.trade_list"
println(dstFile)

df = CSV.File(srcFile; comment="\n", delim=",", types=[String, Int32, Int32, Float64]) |> DataFrame

coef = 0.30 # how many days allowed percentage-wise to have absent trading activity
            # used as an upper bound quantile for the security to be traded

function filterOnColumnNotLower(df::DataFrame, col::String, coef::Float64)
    s_df = copy(df[:, col])
    maxN = maximum(s_df)
    mask = (s_df ./ maxN * 100) .<= coef
    copy(df[mask, :])
end

function filterOnColumnBigger(df::DataFrame, col::String, coef::Float64)
    s_df = copy(df[:, col])
    maxN = maximum(s_df)
    mask = (s_df ./ maxN * 100) .>= coef
    copy(df[mask, :])
end

dfNoIdle = filterOnColumnNotLower(df, "trades_days_idle", 1 - 0.05)
dfUpTrades = filterOnColumnBigger(dfNoIdle, "trades_num_total", 1 - coef)

#conn = RedisConnection(password=cfg["REDIS_PWD"])

open(dstFile, "w") do io
    coin_vec = namedtupleiterator(dfUpTrades[:, "coin"]).x.x
    #l = ReentrantLock()
    Threads.@threads for c in coin_vec
        println(io, c)
        #try
        #finally
        #    unlock(l)
        #end
    end
#    set(conn, "$(curDate)_coins", coin_vec)
end

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
mutable struct Client
    id::Symbol
    stream::Any
    connected_notebook::Union{Notebook,Nothing}
    pendingupdates::Channel
end

Client(id::Symbol, stream) = let
    Client(id, stream, nothing, Channel(1024))
end

struct Initiator
    clientID::Symbol
    requestID::Symbol
end

struct UpdateMessage
    type::Symbol
    message::Any
    notebook::Union{Notebook,Nothing}
    cell::Union{Cell,Nothing}
    initiator::Union{Initiator,Missing}
end

UpdateMessage(type::Symbol, message::Any) = UpdateMessage(type, message, nothing, nothing, missing)
UpdateMessage(type::Symbol, message::Any, notebook::Notebook) = UpdateMessage(type, message, notebook, nothing, missing)

function clientupdate_cell_output(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    payload, mime = cell.output_repr, cell.repr_mime

    return UpdateMessage(:cell_output,
            Dict(:mime => mime,
             :output => payload,
             :errored => cell.errored,
             :runtime => cell.runtime,
             :running => cell.running,
             :rootassignee => cell.rootassignee,
            ),
            notebook, cell, initiator)
end

function clientupdate_cell_input(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_input,
        Dict(:code => cell.code, :folded => cell.code_folded), notebook, cell, initiator)
end

function clientupdate_cell_added(notebook::Notebook, cell::Cell, new_index::Integer; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_added,
        Dict(:index => new_index - 1, # 1-based index (julia) to 0-based index (js)
            ), notebook, cell, initiator)
end

function clientupdate_cell_deleted(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_deleted,
        Dict(), notebook, cell, initiator)
end

function clientupdate_cell_moved(notebook::Notebook, cell::Cell, new_index::Integer; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_moved,
        Dict(:index => new_index - 1, # 1-based index (julia) to 0-based index (js)
            ), notebook, cell, initiator)
end

function clientupdate_cell_running(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_running,
        Dict(), notebook, cell, initiator)
end

function clientupdate_cell_folded(notebook::Notebook, cell::Cell, folded::Bool; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_folded,
        Dict(:folded => folded), notebook, cell, initiator)
end

function clientupdate_notebook_list(notebooks; initiator::Union{Initiator,Missing}=missing)
    update = UpdateMessage(:notebook_list,
        Dict(:notebooks => [Dict(
                :uuid => notebook.uuid,
                :path => notebook.path,
                :shortpath => basename(notebook.path)
                ) for notebook in values(notebooks)]), nothing, nothing, initiator)
end


text:
# Import dependency.
using BenchmarkTools, Expectations, Distributions

#Create benchmark group and benchmarks
benchmarks = BenchmarkGroup()

#Put in specific benchmarks
dist1 = DiscreteUniform(1, 10)
E1 = expectation(dist1)
benchmarks["discreteuniform"] = @benchmarkable $E1($(x -> x))
dist2 = Normal()
E2 = expectation(dist2)
benchmarks["standardnormal"] = @benchmarkable $E2($(x -> x))
dist3 = Exponential(2)
E3 = expectation(dist3)
benchmarks["exponential"] = @benchmarkable $E3($(x -> x))
#...

results = run(benchmarks) # Get results.
results = median(results) # Condense to median.

# To save results, manually call in the REPL: BenchmarkTools.save("benchmarks.json", results)

#Compare to old results
try
  oldresults= BenchmarkTools.load("benchmarks.json")[1]
  judge(oldresults, results)
catch err
  error("Couldn't load file- make sure that you've previously saved results.", err.prefix)
end


text:
# This dictionary maps easy to remember names to Youtube video IDs
# after adding an ID here, you can use the {{youtube <shortname>}}
# syntax in your markdown files to embed the video into the page!
videos = Dict(
              "course-intro" => "ZhC-DIrCe6A",
              "pytorch-tensors" => "BmAS8IH7n3c",
              "automatic-diff" => "Z6H3zakmn6E",
              "loss-functions" => "jReGEZXq4Ac",
              "optim-basics" => "UvM0hK4E2dc",
              "pytorch-module" => "OiyZXdnLHcI",
              "conv" => "HyotcucT-PE",
              "dataloading-emb" => "vm-ZusIUkiY",
              "collab-filt" => "TaEWBB00Els",
              "autoencoder" => "Z1OWiTFafpQ",
              "gan" => "tqhmilLuJRs",
              "rnn-1" => "4G681MnP_OA",
              "rnn-2" => "EsS95AxCnh0",
              "privacy-preserving" => "68aCbiPIgnM",
              "batches-seq" => "Fqx_RCwenfg",
              "intro-julia" => "oL_T_NVoz9w",
              "siamese" => "6x4IPITZ4dw",
              "ben-depth" => "l-J93Onvj70",
              "pbs-depth" => "cBrDJ0KHxCs",
              "dropout" => "gWiJr_y0Tgs",
              "batchnorm" => "cKSCyTODbDI",
              "resnets" => "LdJV4oHgO28",
              "graph1" => "Ch0Iz8BJn98",
              "graph2" => "o1CfrsSTRAU",
              "graph3" => "ZD0zwjEQhz8"

)

max_stars_repo_name: SimonCoste/website

max_issues_repo_name: SimonCoste/website

max_forks_repo_name: SimonCoste/website


text:
# Autogenerated wrapper script for ERFA_jll for x86_64-linux-gnu
export liberfa

JLLWrappers.@generate_wrapper_header("ERFA")
JLLWrappers.@declare_library_product(liberfa, "liberfa.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liberfa,
        "lib/liberfa.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

### fsa2transition
# transition matrix 3 state HMM with self loop in the middle
H = WFST(["a1","a2","a3"],["a"]);

add_arc!(H,1,2,"a1","a");
add_arc!(H,2,2,"a1","<eps>",-log(0.3));
add_arc!(H,2,3,"a2","<eps>",-log(0.7));
add_arc!(H,3,3,"a2","<eps>",-log(0.3));
add_arc!(H,3,4,"a3","<eps>",-log(0.7));
add_arc!(H,4,4,"a3","<eps>",-log(0.3));
add_arc!(H,4,2,"a1","a",-log(0.7));
initial!(H,1);
final!(H,4)
#println(H)

Nt = 10;
time2tr = wfst2tr(H,Nt)
#println(time2tr)

@test length(time2tr) == Nt-1
@test all([ all(keys(tr) .<= length(get_isym(H)) ) for tr in time2tr])

a,A = wfst2tr(H; convert_weight = w -> exp(-get(w)))
@test all(abs.(sum(A,dims=2) .- 1) .< 1e-6)
@test abs(sum(a) - 1) .< 1e-6


text:
t1 = time()
using NiLang

@i function dot(x, y, z)
    for i=1:10
        x += y[i]' * z[i]
    end
end
t2 = time()
println("costs $(t2-t1)s")


text:
using AtomsBase
using Test
using Unitful
using PeriodicTable

@testset "Fast system" begin
    box = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m"
    bcs = [Periodic(), Periodic(), DirichletZero()]
    atoms = Atom[:C => [0.25, 0.25, 0.25]u"m",
                 :C => [0.75, 0.75, 0.75]u"m"]
    system = FastSystem(atoms, box, bcs)

    @test length(system) == 2
    @test size(system)   == (2, )
    @test atomic_mass(system) == [12.011, 12.011]u"u"
    @test boundary_conditions(system) == bcs
    @test bounding_box(system) == box

    # Test AtomView
    for method in (position, atomic_mass, atomic_symbol, atomic_number)
        @test method(system[1]) == method(system, 1)
        @test method(system[2]) == method(system, 2)
    end
    @test ismissing(velocity(system[1]))
    @test n_dimensions(system[1]) == n_dimensions(system)
    @test element(system[1])      == elements[:C]
end


text:
using PaddedViews

to_int(s) = parse(Int, s)

slurp_grid(filename) = to_int.(permutedims(hcat(map(collect, readlines(filename))...)))

function pad(m)
    y, x = size(m)
    PaddedView(10, m, (0:y+1, 0:x+1), (1:y, 1:x))
end

function get_neighbors(idx)
    y, x = Tuple(idx)
    [
        CartesianIndex(y + 1, x),
        CartesianIndex(y - 1, x),
        CartesianIndex(y, x + 1),
        CartesianIndex(y, x - 1),
    ]
end

get_neighbors_values(m, idx) = [m[idx] for idx in get_neighbors(idx)]
function is_low(m, idx)
    current = m[idx]
    neighborhood = get_neighbors_values(m, idx)
    length(filter(<=(current), neighborhood)) == 0
end

function solve1(filename)
    m = slurp_grid(filename)
    p = pad(m)
    sum([1 + p[idx] for idx in CartesianIndices(m) if is_low(p, idx)])
end

solve1("9eg.txt") == 15
solve1("9.txt") == 535

function get_low_points(m)
    p = pad(m)
    [idx for idx in CartesianIndices(m) if is_low(p, idx)]
end

get_basin_neighbors(m, idx) = [idx_ for idx_ in get_neighbors(idx) if m[idx_] < 9]

function basin_size(m, pos)
    m = pad(m)
    explored = Set()
    unexplored = Set([pos])
    while !isempty(unexplored)
        current = pop!(unexplored)
        push!(explored, current)
        neighbors = get_basin_neighbors(m, current)
        for neighbor in neighbors
            if neighbor ∉ explored
                push!(unexplored, neighbor)
            end
        end
    end
    length(explored)
end

function solve2(filename)
    m = slurp_grid(filename)
    sorted = sort([basin_size(m, pos) for pos in get_low_points(m)], rev = true)
    reduce(*, sorted[1:3])
end

solve2("9eg.txt")
solve2("9.txt")


text:
@testset "VectorCacheInt64" begin
    vc = VectorCache{Int64}()
    @test length(vc) == -9999
    empty!(vc)
    @test isempty(vc)
    @test length(vc) == 0
    @test returnNext(vc) isa Int64

    empty!(vc)
    ind_max_1 = vc.ind_max
    for k = 1:ind_max_1
        addCacheItem!(vc, k)
        @test vc[k] == k
    end
    @test !isempty(vc)

    expand!(vc)
    ind_max_2 = vc.ind_max
    @test length(vc) == vc.ind_fill
    @test ind_max_1 < ind_max_2

    vc.ind_fill = vc.ind_max
    i_next = returnNext(vc)
    ind_max_3 = vc.ind_max
    @test ind_max_2 < ind_max_3
end


text:
# Transformations
# ---------------

"""
    push!(seq::BioSequence{A}, x) where {A}

Append a biological symbol `x` to a biological sequence `seq`.
"""
function Base.push!(seq::BioSequence{A}, x) where {A}
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    encoded_setindex!(seq, bin, lastindex(seq))
    return seq
end

"""
    pop!(seq::BioSequence)

Remove the symbol from the end of a biological sequence `seq` and return it.
Returns a variable of `eltype(seq)`.
"""
function Base.pop!(seq::BioSequence)
    if isempty(seq)
        throw(ArgumentError("sequence must be non-empty"))
    end
    x = seq[end]
    deleteat!(seq, lastindex(seq))
    return x
end

"""
    insert!(seq, i, x)

Insert a biological symbol `x` into a biological sequence `seq`, at the given
index `i`.
"""
function Base.insert!(seq::BioSequence{A}, i::Integer, x) where {A}
    checkbounds(seq, i)
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    copyto!(seq, i + 1, seq, i, lastindex(seq) - i)
    encoded_setindex!(seq, bin, i)
    return seq
end

"""
    deleteat!(seq::BioSequence, range::UnitRange{<:Integer})

Deletes a defined `range` from a biological sequence `seq`.

Modifies the input sequence.
"""
function Base.deleteat!(seq::BioSequence{A}, range::UnitRange{<:Integer}) where {A}
    checkbounds(seq, range)
    copyto!(seq, range.start, seq, range.stop + 1, length(seq) - range.stop)
    resize!(seq, length(seq) - length(range))
    return seq
end

"""
    deleteat!(seq::BioSequence, i::Integer)

Delete a biological symbol at a single position `i` in a biological sequence
`seq`.

Modifies the input sequence.
"""
function Base.deleteat!(seq::BioSequence, i::Integer)
    checkbounds(seq, i)
    copyto!(seq, i, seq, i + 1, length(seq) - i)
    resize!(seq, length(seq) - 1)
    return seq
end

"""
    append!(seq, other)

Add a biological sequence `other` onto the end of biological sequence `seq`.
Modifies and returns `seq`.
"""
function Base.append!(seq::BioSequence{A}, other::BioSequence{A}) where {A}
    resize!(seq, length(seq) + length(other))
    copyto!(seq, lastindex(seq) - length(other) + 1, other, 1)
    return seq
end

"""
    popfirst!(seq)

Remove the symbol from the beginning of a biological sequence `seq` and return
it. Returns a variable of `eltype(seq)`.
"""
function Base.popfirst!(seq::BioSequence)
    if isempty(seq)
        throw(ArgumentError("sequence must be non-empty"))
    end
    x = seq[1]
    deleteat!(seq, 1)
    return x
end

"""
    pushfirst!(seq, x)

Insert a biological symbol `x` at the beginning of a biological sequence `seq`.
"""
function Base.pushfirst!(seq::BioSequence{A}, x) where {A}
    bin = enc64(seq, x)
    resize!(seq, length(seq) + 1)
    copyto!(seq, 2, seq, 1, length(seq) - 1)
    encoded_setindex!(seq, bin, 1)
    return seq
end

"""
    resize!(seq, size)

Resize a biological sequence `seq`, to a given `size`.
"""
function Base.resize!(seq::BioSequence{A}, size::Integer) where {A}
    if size < 0
        throw(ArgumentError("size must be non-negative"))
    end
    orphan!(seq, size)
    resize!(seq.data, seq_data_len(A, size + seq.part.start - 1))
    seq.part = seq.part.start:seq.part.start+size-1
    return seq
end

"""
    empty!(seq)

Completely empty a biological sequence `seq` of nucleotides.
"""
Base.empty!(seq::BioSequence) = resize!(seq, 0)

function Base.filter!(f::Function, seq::BioSequence{A}) where {A}
    orphan!(seq)

    len = 0
    next = bitindex(seq, 1)
    j = index(next)
    datum::UInt64 = 0
    for i in 1:lastindex(seq)
        x = inbounds_getindex(seq, i)
        if f(x)
            datum |= enc64(seq, x) << offset(next)
            len += 1
            next += bitsof(A)
            if index(next) != j
                seq.data[j] = datum
                datum = 0
                j = index(next)
            end
        end
    end
    if offset(next) > 0
        seq.data[j] = datum
    end
    resize!(seq, len)

    return seq
end

function Base.filter(f::Function, seq::BioSequence)
    return filter!(f, copy(seq))
end

function Base.map!(f::Function, seq::BioSequence)
    orphan!(seq)
    for i in 1:lastindex(seq)
        unsafe_setindex!(seq, f(inbounds_getindex(seq, i)), i)
    end
    return seq
end

function Base.map(f::Function, seq::BioSequence)
    return map!(f, copy(seq))
end

"""
    reverse!(seq)

Reverse a biological sequence `seq` in place.
"""
function Base.reverse!(seq::BioSequence)
    orphan!(seq)
    for i in 1:div(lastindex(seq), 2)
        x = inbounds_getindex(seq, i)
        unsafe_setindex!(seq, inbounds_getindex(seq, lastindex(seq) - i + 1), i)
        unsafe_setindex!(seq, x, lastindex(seq) - i + 1)
    end
    return seq
end

"""
    reverse(seq)

Create a sequence which is the reverse of the bioloigcal sequence `seq`.
"""
Base.reverse(seq::BioSequence) = reverse!(copy(seq))

@generated function Base.reverse(seq::BioSequence{A}) where {A<:NucAlphs}
    n = bitsof(A)
    if n == 2
        nucrev = :nucrev2
    elseif n == 4
        nucrev = :nucrev4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        data = Vector{UInt64}(undef, seq_data_len(A, length(seq)))
        i = 1
        next = bitindex(seq, lastindex(seq))
        stop = bitindex(seq, 0)
        r = rem(offset(next) + $n, 64)
        if r == 0
            @inbounds while next - stop > 0
                x = seq.data[index(next)]
                data[i] = $nucrev(x)
                i += 1
                next -= 64
            end
        else
            @inbounds while next - stop > 64
                j = index(next)
                x = (seq.data[j] << (64 - r)) | (seq.data[j-1] >> r)
                data[i] = $nucrev(x)
                i += 1
                next -= 64
            end
            if next - stop > 0
                j = index(next)
                x = seq.data[j] << (64 - r)
                if r < next - stop
                    x |= seq.data[j-1] >> r
                end
                data[i] = $nucrev(x)
            end
        end
        return BioSequence{A}(data, 1:length(seq), false)
    end
end

@inline function nucrev2(x::UInt64)
     x = (x & 0x3333333333333333) <<  2 | (x & 0xCCCCCCCCCCCCCCCC) >>  2
     x = (x & 0x0F0F0F0F0F0F0F0F) <<  4 | (x & 0xF0F0F0F0F0F0F0F0) >>  4
     x = (x & 0x00FF00FF00FF00FF) <<  8 | (x & 0xFF00FF00FF00FF00) >>  8
     x = (x & 0x0000FFFF0000FFFF) << 16 | (x & 0xFFFF0000FFFF0000) >> 16
     x = (x & 0x00000000FFFFFFFF) << 32 | (x & 0xFFFFFFFF00000000) >> 32
     return x
end

@inline function nucrev4(x::UInt64)
     x = (x & 0x0F0F0F0F0F0F0F0F) <<  4 | (x & 0xF0F0F0F0F0F0F0F0) >>  4
     x = (x & 0x00FF00FF00FF00FF) <<  8 | (x & 0xFF00FF00FF00FF00) >>  8
     x = (x & 0x0000FFFF0000FFFF) << 16 | (x & 0xFFFF0000FFFF0000) >> 16
     x = (x & 0x00000000FFFFFFFF) << 32 | (x & 0xFFFFFFFF00000000) >> 32
     return x
end

"""
    complement!(seq)

Make a complement sequence of `seq` in place.
"""
function complement!(seq::BioSequence{A}) where {A<:TwoBitNucs}
    orphan!(seq)
    next = bitindex(seq, 1)
    stop = bitindex(seq, lastindex(seq) + 1)
    @inbounds while next < stop
        seq.data[index(next)] = ~seq.data[index(next)]
        next += 64
    end
    return seq
end

"""
    complement!(seq)

Transform `seq` into it's complement.
"""
function complement!(seq::BioSequence{A}) where {A<:FourBitNucs}
    orphan!(seq)
    next = bitindex(seq, 1)
    stop = bitindex(seq, lastindex(seq) + 1)
    @inbounds while next < stop
        x = seq.data[index(next)]
        seq.data[index(next)] = (
            ((x & 0x1111111111111111) << 3) | ((x & 0x8888888888888888) >> 3) |
            ((x & 0x2222222222222222) << 1) | ((x & 0x4444444444444444) >> 1))
        next += 64
    end
    return seq
end

"""
    complement(seq)

Make a complement sequence of `seq`.
"""
function BioSymbols.complement(seq::BioSequence{A}) where {A<:NucAlphs}
    return complement!(copy(seq))
end

"""
    reverse_complement!(seq)

Make a reversed complement sequence of `seq` in place.
"""
function reverse_complement!(seq::BioSequence{A}) where {A<:NucAlphs}
    return complement!(reverse!(seq))
end

"""
    reverse_complement(seq)

Make a reversed complement sequence of `seq`.
"""
function reverse_complement(seq::BioSequence{A}) where {A<:NucAlphs}
    return complement!(reverse(seq))
end

# Shuffle
# -------

function Random.shuffle(seq::BioSequence)
    return shuffle!(copy(seq))
end

function Random.shuffle!(seq::BioSequence)
    orphan!(seq)
    # Fisher-Yates shuffle
    for i in 1:lastindex(seq)-1
        j = rand(i:lastindex(seq))
        seq[i], seq[j] = seq[j], seq[i]
    end
    return seq
end


text:
function scaleandcenter!(x::AbstractArray; center = 0, range = 1)
    minn, maxx = extrema(x)
    if maxx ≈ minn
        x .= center
    else
        @. x = ((x - minn) / (maxx - minn) - 1 // 2) * range + center
    end
    x
end

@inline function symmul!(
    C::Symmetric{T},
    transA::Transpose{<:Any,<:StridedVecOrMat{T}},
    B::StridedVecOrMat{T},
    alpha::Number = true,
    beta::Number = false,
) where {T<:LinearAlgebra.BlasFloat}
    A = transA.parent
    if A === B
        alpha = convert(T, alpha)
        beta = convert(T, beta)
        BLAS.syrk!(C.uplo, 'T', alpha, A, beta, C.data)
        return C
    else
        return LinearAlgebra.gemm_wrapper!(C, 'T', 'N', A, B, alpha, beta)
    end
end

@inline function symmul!(
    C::Symmetric{T},
    A::StridedVecOrMat{T},
    transB::Transpose{<:Any,<:StridedVecOrMat{T}},
    alpha::Number = true,
    beta::Number = false,
) where {T<:LinearAlgebra.BlasFloat}
    B = transB.parent
    if A === B
        alpha = convert(T, alpha)
        beta = convert(T, beta)
        BLAS.syrk!(C.uplo, 'N', alpha, B, beta, C.data)
        return C
    else
        return LinearAlgebra.gemm_wrapper!(C, 'N', 'T', A, B, alpha, beta)
    end
end

"""
    $(SIGNATURES)

Wrap the angle `theta`, specified in radians, to the interval [-pi, pi)
"""
@inline wraptopi(theta::Real) = mod1(theta + pi, 2pi) - pi

mutable struct Converged
    atol::Float64
    rtol::Maybe{Float64}
    nans::Bool
    lastval::Number
    initialized::Bool
    Converged(atol, rtol, nans) = new(atol, rtol, nans, 0, false)
end
Converged(; atol = 0, rtol = nothing, nans = false) = Converged(atol, rtol, nans)

function (x::Converged)(val::Number)
    if x.initialized
        rtol = x.rtol !== nothing ? x.rtol : Base.rtoldefault(x.lastval, val, x.atol)
        converged = isapprox(x.lastval, val, atol = x.atol, rtol = rtol, nans = x.nans)
    else
        x.initialized = true
        converged = false
    end
    x.lastval = val
    return converged
end

reset!(x::Converged) = (x.initialized = false; x.lastval = 0; x)

perturb(rng::AbstractRNG, A::AbstractArray) = perturb!(rng, copy(A))
perturb!(A::AbstractArray) = perturb!(Random.default_rng(), A)
function perturb!(rng::AbstractRNG, A::AbstractArray)
    for i in eachindex(A)
        @inbounds A[i] += rand(rng)
    end
    A
end

perturbn(rng::AbstractRNG, A::AbstractArray) = perturb!(rng, copy(A))
perturbn!(A::AbstractArray) = perturb!(Random.default_rng(), A)
function perturbn!(rng::AbstractRNG, A::AbstractArray)
    for i in eachindex(A)
        @inbounds A[i] += randn(rng)
    end
    A
end

perturb(rng::AbstractRNG, s::Sampleable, A::AbstractArray) = perturb!(rng, s, copy(A))
perturb!(s::Sampleable, A::AbstractArray) = perturb!(Random.default_rng(), s, A)
function perturb!(rng::AbstractRNG, s::Sampleable, A::AbstractArray)
    for i in eachindex(A)
        @inbounds A[i] += rand(rng, s)
    end
    A
end


text:
using SymbolicUtils
using Test

using SymbolicUtils: showraw

const leaf_funcs = [()->100*randn(),
                    ()->rand(-100:100),
                    ()->rand(@syms a b c d e f),
                    ()->rand(@syms a b c d e f)]

const fns = vcat(1 .=> SymbolicUtils.monadic,
                 2 .=> vcat(SymbolicUtils.diadic, fill(+, 5), [-,-], fill(*, 5)),
                 3 .=> [+, *])

function gen_rand_expr(inputs; leaf_prob=0.92, depth=0, min_depth=1, max_depth=5)
    if depth > max_depth  || (min_depth <= depth && rand() < leaf_prob)
        leaf = rand(leaf_funcs)()
        if leaf isa SymbolicUtils.Sym
            push!(inputs, leaf)
        end
        return leaf
    end
    arity, f = rand(fns)
    args = [gen_rand_expr(inputs,
                          leaf_prob=leaf_prob,
                          depth=depth+1,
                          min_depth=min_depth,
                          max_depth=max_depth) for i in 1:arity]
    try
        return f(args...)
    catch err
        if err isa DomainError
            return gen_rand_expr(inputs,
                                 leaf_prob=leaf_prob,
                                 depth=depth,
                                 min_depth=min_depth,
                                 max_depth=max_depth)
        else
            rethrow(err)
        end
    end
end

struct Errored
    err
end

function fuzz_test(ntrials)
    inputs = Set()
    expr = gen_rand_expr(inputs)
    inputs = collect(inputs)
    unsimplifiedstr = """
    function $(tuple(inputs...))
        $(sprint(io->showraw(io, expr)))
    end
    """

    simplifiedstr = """
    function $(tuple(inputs...))
        $(sprint(io->showraw(io, simplify(expr))))
    end
    """
    f = include_string(Main, unsimplifiedstr)
    g = include_string(Main, simplifiedstr)

    for i=1:ntrials
        args = [randn() for j in 1:length(inputs)]
        unsimplified = try
            Base.invokelatest(f, args...)
        catch err
            Errored(err)
        end

        simplified  = try
            Base.invokelatest(g, args...)
        catch err
            Errored(err)
        end

        if unsimplified isa Errored
            @test typeof(simplified.err) == typeof(unsimplified.err)
        else
            try
                @test unsimplified ≈ simplified
                if !(unsimplified ≈ simplified)
                    error("Failed")
                end
            catch err
                println("""Test failed for expression
                               $(sprint(io->showraw(io, expr))) = $unsimplified
                           Simplified to:
                               $(sprint(io->showraw(io, simplify(expr)))) = $simplified
                           On inputs:
                               $inputs = $args
                           """)
            end
        end
    end
end

using Random: seed!

@testset "Fuzz test" begin

    seed!(6174)
    for i=1:2000
        fuzz_test(10)
    end
end