text:
# ---------------NEW ------------------
using DataFrames
dfd=readtable("orig.csv")
names!(dfd,[Symbol(replace(string(n),"jennie_rf_pos_exposure2_","")) for n in names(dfd)])
rename!(dfd,Symbol("3rdparty"),:thirdparty)

# --------------- NEW END -----------------

# ---------------- REMOVE OTLIERS -------------
using DataFrames
dfd=readtable("out.csv")

isNum{T<:Number}(::AbstractArray{T}) = true
isNum(::Any) = false
isFloat{T<:Float64}(::AbstractArray{T}) = true
isFloat(::Any) = false
isInt{T<:Int64}(::AbstractArray{T}) = true
isInt(::Any) = false
isString{T<:String}(::AbstractArray{T}) = true
isString(::Any) = false
isBool{T<:Bool}(::AbstractArray{T}) = true
isBool(::Any) = false
isSymbol{T<:Symbol}(::AbstractArray{T}) = true
isSymbol(::Any) = false

describe(dfd[:venue_dim_key])
countmap(dfd[:venue_dim_key])

 :venue_dim_key
 :tm_dim_key_day
-- don't change :date1
 :trans_time
???? there is exactly 10000? this should either be unique, or the same/defaulted  :trans_num
 :store_num
 :card_num
 :upc
--- not populated  :ndc_upc_ind
 :quantity
--- not populated :weight
--- overfit :item_list_price
 :net_price
--- not populated :card_discount
--- not populated  :ad_discount
--- not populated  :other_discount
--- not populated :item_type
--- not populated :tender_type
--- not populated :ad_event
--- not populated  :ad_version
:tm_dim_key_week  -- only in control
:operating_company -- only "BJS" in control  & COMSCORE in exp
:item_dim_key - only populated for control
 :upc10 - only populated for control
 :units - only populated for control
 :cents - only populated for control
 :baseline_units  - only populated for control
 :baseline_cents - only populated for control
 :feature - only populated for control
 :display - only populated for control
 :totl_prc_reduc
 :registration_request_id
 :product_id_jennieo6

for c in setdiff(names(dfd[1:33]),[:])
    #b = 2.8*(sqrt(sum((dfd[:quantity].- mean(dfd[:quantity])).^2)/(length(dfd[:quantity])-1)))
    b = 2.8*(sqrt(sum((dfd[(dfd[:gross_imps].==0),:quantity].- mean(dfd[(dfd[:gross_imps].==0),:quantity])).^2)/(length(dfd[(dfd[:gross_imps].==0),:quantity])-1)))

    length(dfd[(dfd[:quantity].> b)|(dfd[:quantity].< (b*-1)),:quantity])
end

function Pre_out(df_in::DataFrame)
     df_cat_pre = df_in[df_in[:Buyer_Pre_P0] .==1 , [:Prd_0_Net_Pr_PRE,:experian_id]]
     df_cat_pos = df_in[(df_in[:Buyer_Pre_P0] .==0) & (df_in[:Buyer_Pos_P0] .==1) , [:experian_id]]
     median_df = median(df_cat_pre[:Prd_0_Net_Pr_PRE])
     df_cat_pre[:Prd_0_Net_Pr_PRE_med1] = abs(df_cat_pre[:Prd_0_Net_Pr_PRE]-median_df)
     MAD=median(df_cat_pre[:Prd_0_Net_Pr_PRE_med1])
     df_cat_pre[:Prd_0_Net_Pr_PRE_med2] = (0.6745*(df_cat_pre[:Prd_0_Net_Pr_PRE]-median_df))/MAD
     df_cat_pre_zsc = df_cat_pre[abs(df_cat_pre[:Prd_0_Net_Pr_PRE_med2]) .< 3.5,:]
     df_cat_pre_zsc_1 = df_cat_pre_zsc[:,[:experian_id]]
     df_cat_pre_zsc_f = vcat(df_cat_pos,df_cat_pre_zsc_1)
     df_in_pout =  join(df_in, df_cat_pre_zsc_f, on =  :experian_id , kind = :inner);
end
df_in = Pre_out(df_in);

# ---------------- END REMOVE OTLIERS -------------

"""
# DataSet
drop table if exists jennie.jennie6_breaks_RF;
CREATE EXTERNAL TABLE IF NOT EXISTS jennie.jennie6_breaks_RF
  ( household_id string,
    iri_week string,
    date1 string,
    creative_id int,
    placement_id int,
    gross_imps bigint,
    val_imps bigint,
    placement_nm string,
    creative_nm string,
    publisher string
   )
ROW FORMAT DELIMITED FIELDS TERMINATED BY ',' LINES TERMINATED BY ‘\n’ STORED AS TEXTFILE LOCATION '/mapr/mapr04p/analytics0001/analytic_users/Models/trees/tables/';

# show create table jennie.jennie6_breaks_RF;

Insert Overwrite Table jennie.jennie6_breaks_RF
select distinct
b.household_id ,
a.iri_week ,
a.date1 ,
a.creative_id ,
a.placement_id ,
a.gross_imps ,
a.val_imps ,
c.placement_nm,
d.creative_nm,
c.publisher
from daily.daily_unioned a
left join jennie.placements6 c on a.placement_id = c.id
left join jennie.creatives6 d on a.creative_id = d.id
inner JOIN WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP b on a.lr_id = b.dependent_id
where b.dependent_type = 'LVRAMP_ID'
      and b.household_type = 'EXPERIAN'
      and b.retailer = 'COMSCORE'
      and b.data_supplier = 'EXPERIAN'
      and b.current_state = 'MATCHED'
      and a.clientid in('21884504')
      and iri_week < 1934
order by household_id, date1;

=============================================================================

select distinct
b.household_id ,
a.iri_week ,
a.date1 ,
a.creative_id ,
a.placement_id ,
a.gross_imps ,
a.val_imps
from daily.daily_unioned a,
WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP b
where a.lr_id = b.dependent_id
      and b.dependent_type = 'LVRAMP_ID'
      and b.household_type = 'EXPERIAN'
      and b.retailer = 'COMSCORE'
      and b.data_supplier = 'COMSCORE'
      and b.current_state = 'MATCHED'
      and a.clientid in('21884504')
      and iri_week < 1934
order by household_id, date1;

select count(*)
from daily.daily_unioned a,
     WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP b
where a.lr_id = b.dependent_id
      and a.clientid in('21884504')
      and iri_week < 1934
      and b.current_state = 'MATCHED'
      and b.retailer = 'COMSCORE'
      and b.data_supplier = 'EXPERIAN'
;

select count(*)
from WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP b
where b.retailer = 'COMSCORE'
      and b.data_supplier = 'COMSCORE'
;

select distinct(retailer) from WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP;

select distinct(data_supplier) from WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP;
"""

# ========================================================== DATASET AGREGATION =============================================

# ================ START HERE ==============================
export JULIA_NUM_THREADS=4

Threads.nthreads()
#addprocs([("10.106.128.212", 1)])
addprocs([("10.106.128.213", 1)])
addprocs([("10.106.128.214", 1)])
addprocs([("10.106.128.215", 1)])
addprocs([("10.106.128.216", 1)])
addprocs([("10.106.128.217", 1)])
addprocs([("10.106.128.218", 1)])
addprocs([("10.106.128.219", 1)])
addprocs([("10.106.128.220", 1)])
addprocs([("10.106.128.221", 1)])
addprocs(1)

#Base.LOAD_CACHE_PATH[1]="/mapr/mapr04p/analytics0001/analytic_users/jpkg/v0.5"
#ENV["JULIA_PKGDIR"]  = "/mapr/mapr04p/analytics0001/analytic_users/jpkg"
#push!(LOAD_PATH, "/mapr/mapr04p/analytics0001/analytic_users/jpkg/v0.5")
#pop!(LOAD_PATH)

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

#dfd = readtable("/mnt/resource/analytics/trees/exposure.csv",header=false);
dfd = readtable("orig.csv",header=true);
#names!(dfdx, [:household_id, :iri_week, :date1, :creative_id, :placement_id, :gross_imps, :val_imps ,:placement, :creative, :publisher] )
#names!(dfd, [:hh_id, :week, :date, :creative_id, :placement_id, :gross_imps, :imps ,:placement, :creative, :publisher] )
names!(dfd,[Symbol(replace(string(n),"jennie_rf_pos_exposure2_","")) for n in names(dfd)])
rename!(dfd,Symbol("3rdparty"), :thirdrdparty )
rename!(dfd,:household_id, :panid )

#dfd[:household_id]
dformat = Dates.DateFormat("y-m-d");
#dfd[:date] = map(x-> x=="NULL" ? NA : Date(replace(x, " 00:00:00",""),dformat), dfd[:date1])
dfd[:date] = map(x-> Date(string(x),Dates.DateFormat("yyyymmdd") ), dfd[:date1])

dfd = dfd[!isna(dfd[:date]),:]
dfd[:date] = convert(Array{Date},dfd[:date])
dfd[:datelag] = dfd[:date] .- Dates.Day(7)
dfd[:datenum] = map(x->Int64(x),dfd[:date])
dfd[:datenumlag] = map(x->Int64(x),dfd[:datelag])

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

#Date("2016-07-06",Dates.DateFormat("y-m-d"))

# ----------------- CHUNKS  ------------------
#Date("2016-06-30",Dates.DateFormat("y-m-d"))-Dates.Day(7)
#Date("2016-06-30",Dates.DateFormat("y-m-d"))
#dfd[:date] = map(x-> Date(replace(x, " 00:00:00",""),dformat), dfd[:date1])
#Date(dfd[:date1],dformat)
#Date(dfd[:date1],dformat)
#Date("2016-06-30",dformat)
#dfd=sort(dfd,cols=[:date])
#dfd[:chunks]=-1

panlst = unique(dfd[:panid])
panchunks=DataFrame( panid=panlst,chunks=-1 )
panlen=length(panlst)
chunksize = Int(floor(length(panlst)/100))
chunknum=1
while (chunknum*chunksize) <= panlen
    s=chunknum*chunksize
    e=s+chunksize
    println("start: ",s,"   end:",e,"   :=  chunk : ",chunknum)
    if e > panlen
        panchunks[s:end,:chunks] = chunknum
    else
        panchunks[s:e,:chunks] = chunknum
    end
    chunknum=chunknum+1
end
panchunks[panchunks[:chunks].==-1,:chunks]=chunknum

countmap(panchunks[:chunks])

dfd2 = join(dfd, panchunks, on = :panid)

root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
writetable(root*"/total.csv", dfd2)
for p in sort(unique(dfd2[:chunks]))
    println(p)
    writetable(root*"/"*string(p)*".csv",  dfd2[dfd2[:chunks].==p,:]    )
end

# ---- END CHUNKS ------

#---  Generate Lags ---

export JULIA_NUM_THREADS=4
Threads.nthreads()

#addprocs([("10.106.128.212", 1)])
addprocs([("10.106.128.213", 1)])
addprocs([("10.106.128.214", 1)])
addprocs([("10.106.128.215", 1)])
addprocs([("10.106.128.216", 1)])
addprocs([("10.106.128.217", 1)])
addprocs([("10.106.128.218", 1)])
addprocs([("10.106.128.219", 1)])
addprocs([("10.106.128.220", 1)])
addprocs([("10.106.128.221", 1)])
addprocs(1)

@everywhere insert!(Base.LOAD_CACHE_PATH, 1, "/mapr/mapr04p/analytics0001/analytic_users/jpkg/lib/v0.5")
@everywhere pop!(Base.LOAD_CACHE_PATH)
@everywhere println((Base.LOAD_CACHE_PATH))

using DataFrames
@everywhere using DataFrames
p=1

@everywhere function genLags(p::Int64)
    brks=[:html5
         ,:static
         ,:native_creative
         ,:video
         ,:content
         ,:contextual
         ,:direct
         ,:prospecting
          ,:retargeting
          ,:behavorial
          ,:native
          ,:predictive
          ,:thirdrdparty
          ,:pmp
          ,:buzzfeed
          ,:hulu
          ,:popsugarus
          ,:rachaelraymag
          ,:shape
          ,:blank
          ,:eatingwell
          ,:meredithcorporation
          ,:realsimple
          ,:turndsp
          ,:cookinglight
          ,:dataxu
          ,:allrecipes
          ,:amazon
          ,:myfitnesspal
          ,:womenshealth
          ,:yummly
         ,:youtube
        ]

    root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
    dfx = readtable(root*"/"*string(p)*".csv",header=true)

    for brk in brks
        lag=Symbol(string(brk)*"_lag")
        println("Generating : ",brk,"  ~ ",lag)
        dfx[lag]=0
        for row in eachrow(dfx)
            #println(r[:video])
            row[lag] = sum(dfx[(dfx[:panid].==row[:panid])&(dfx[:datenum].<=row[:datenum] ) & ( dfx[:datenum].>= row[:datenumlag]  ) ,brk])
        end
    end
    writetable(root*"/"*string(p)*"_out.csv", dfx)
    println("COMPLETE for : ",p,"  in ","/"*string(p)*"_out.csv")
end

"""
@async remotecall_fetch(genLags, 2, 2)
@spawnat 2 genLags(2)

@async remotecall_fetch(genLags, 2, 2)
@async remotecall_fetch(genLags, 3, 3)
@async remotecall_fetch(genLags, 4, 4)
@async remotecall_fetch(genLags, 5, 5)
@async remotecall_fetch(genLags, 6, 6)
@async remotecall_fetch(genLags, 7, 7)
@async remotecall_fetch(genLags, 8, 8)
@async remotecall_fetch(genLags, 9, 9)
@async remotecall_fetch(genLags, 10, 10)
@async remotecall_fetch(genLags, 11, 10)
"""
"""
@async remotecall_fetch(genLags, 2, 11)
@async remotecall_fetch(genLags, 3, 12)
@async remotecall_fetch(genLags, 4, 13)
@async remotecall_fetch(genLags, 5, 14)
@async remotecall_fetch(genLags, 6, 15)
@async remotecall_fetch(genLags, 7, 16)
@async remotecall_fetch(genLags, 8, 17)
@async remotecall_fetch(genLags, 9, 18)
@async remotecall_fetch(genLags, 10, 19)
@async remotecall_fetch(genLags, 11, 20)
"""

x=99
x=59

for p in 2:11
    x=x+1
    @async remotecall_fetch(genLags, p,  x)
end

# ---------------- de-chunk files ----------------------
using DataFrames
root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
#dfx = readtable(root*"/1_out.csv",header=true);

#for i in 31:60   #2:101
#    println("Loading : ",i,"_out.csv")
#    t = readtable(root*"/"*string(i)*"_out.csv",header=true);
#    dfx=vcat(dfx,t)
#end

#writetable(root*"/1_60_tmp.csv", dfx)

using DataFrames
root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
function consolidateChunks(s::Int64,e::Int64)
    println("Loading : ",string(s),"_out.csv")
    dfx = readtable(root*"/"*string(s)*"_out.csv",header=true);
    for i in s+1:e   #2:101
        println("Loading : ",i,"_out.csv")
        t = readtable(root*"/"*string(i)*"_out.csv",header=true);
        dfx=vcat(dfx,t)
    end
    return dfx
end

dfx1 = consolidateChunks(1,30);
dfx2 = consolidateChunks(31,60);
dfx3 = consolidateChunks(61,90);
dfx4 = consolidateChunks(91,101);

dfout=vcat(dfx1,dfx2,dfx3,dfx4)
writetable(root*"/out.csv", dfout)

# ------------------ XGB --------------------------

using BinDeps, DataFrames, XGBoost
#using Dates
#root="/mnt/resource/analytics/models"
root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
#dfd=readtable(root*"/out.csv")
#dfd=readtable(root*"/48_out.csv")
dfd=readtable(root*"/1_60_tmp.csv")
#:datenum, :datenumlag , :chunks

dfd=dfd[setdiff(names(dfd),[:date,:datelag,:date1,:chunks])]
dfd[:trans_time]=pool(dfd[:trans_time])
dfd[:operating_company]=pool(dfd[:operating_company])

#just to get it working - these are strings
dfd=dfd[setdiff(names(dfd),[:trans_time,:ndc_upc_ind,:item_type,:ad_event,:operating_company,:upc10])]

net_price=convert(Array{Float32}, dfd[:net_price])

dfd=dfd[setdiff(names(dfd),[:net_price])]

lags=Symbol[:html5_lag,:static_lag,:native_creative_lag,:video_lag,:content_lag,:contextual_lag,:direct_lag,:prospecting_lag,:retargeting_lag,:behavorial_lag,:native_lag,:predictive_lag,:thirdrdparty_lag,:pmp_lag,:buzzfeed_lag,:hulu_lag,:popsugarus_lag,:rachaelraymag_lag,:shape_lag,:blank_lag,:eatingwell_lag,:meredithcorporation_lag,:realsimple_lag,:turndsp_lag,:cookinglight_lag,:dataxu_lag,:allrecipes_lag,:amazon_lag,:myfitnesspal_lag,:womenshealth_lag,:yummly_lag,:youtube_lag]
dfdc=deepcopy(dfd)
for l in lags
    dfdc[l] = 0.0
end

function DMdata(dfx::DataFrame,lbl::Float32[]=Float32[1.1])
    arr = convert(Array{Float32}, dfx )
    if length(lbl) == 0
        return DMatrix(arr)
    else
        return DMatrix(arr, label = lbl)
    end
end
#Transform data to XGBoost matrices
#trainArray = convert(Array{Float32}, dfd )
#dtrain = DMatrix(trainArray, label = net_price)
#testArray = convert(Array{Float32}, dfd )
#dtest = DMatrix(testArray)
#ctrlArray = convert(Array{Float32}, dfdc )
#dctrl = DMatrix(ctrlArray)
dtrain = DMatrix(convert(Array{Float32},dfd), label = net_price)
dtest = DMatrix(convert(Array{Float32},dfd))
dctrl = DMatrix(convert(Array{Float32},dfdc))

num_round = 250; param = ["eta" => 0.2, "max_depth" => 20, "objective" => "reg:linear", "silent" => 1]
XGBoostModel = xgboost(dtrain, num_round, param = param)

or

#random forrest
param = ["eta" => 0.2, "max_depth" => 5, "objective" => "reg:linear", "silent" => 1,
         "num_parallel_tree" => 1000, "subsample" => 0.5, "colsample_bytree" => 0.5
        ]
XGBoostModel = xgboost(dtrain, 10, param = param)

ptest = XGBoost.predict(XGBoostModel, dtest)
pctrl = XGBoost.predict(XGBoostModel, dctrl)
mean(convert(Array{Float64},ptest))
mean(convert(Array{Float64},pctrl))
mean(ptest) - mean(pctrl)

#importance(XGBoostModel)
im = [ parse(x.fname[2:end]) for x in importance(XGBoostModel)]
icols = [names(dfd)[c] for c in [ parse(x.fname[2:end]) for x in importance(XGBoostModel)]]

# ************************* TEST **** TEST ********************************************
# ************************* TEST **** TEST ********************************************
# ************************* TEST **** TEST ********************************************
# ************************* TEST **** TEST ********************************************
# ************************* TEST **** TEST ********************************************
# ************************* TEST **** TEST ********************************************
# ************************* TEST **** TEST ********************************************

using BinDeps, DataFrames, XGBoost
root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
dfx=readtable(root*"/out.csv")  #dfx=readtable(root*"/1_60_tmp.csv")
net_price=convert(Array{Float32}, dfx[:net_price])
dfxc=deepcopy(dfx)
for l in lags dfxc[l] = 0.0 end
xvars = [:panid, :net_price, :date,:datelag,:date1,:chunks, :trans_time, :trans_num, :operating_company, :ndc_upc_ind,:item_type,:ad_event,   ]

"""
:product_id_jennieo6
 :panid
 :datenum
 :datenumlag
 :trans_num
 :card_num
 :registration_request_id
 :upc

 :venue_dim_key
 :tm_dim_key_day
 :item_list_price

 :tm_dim_key_week

 :item_dim_key
 :units
 :cents
 :baseline_units
 :baseline_cents
 :totl_prc_reduc

 :gross_imps
"""

v_covar=[ :feature,
          :display,
          :ad_version,
          :card_discount,
          :ad_discount,
          :other_discount,
          :tender_type,
          :quantity,
          :weight,
          :store_num
]

v_breaks=[ :html5,
 :static,
 :native_creative,
 :video,
 :content,
 :contextual,
 :direct,
 :prospecting,
 :retargeting,
 :behavorial,
 :native,
 :predictive,
 :thirdrdparty,
 :pmp,
 :buzzfeed,
 :hulu,
 :popsugarus,
 :rachaelraymag,
 :shape,
 :blank,
 :eatingwell,
 :meredithcorporation,
 :realsimple,
 :turndsp,
 :cookinglight,
 :dataxu,
 :allrecipes,
 :amazon,
 :myfitnesspal,
 :womenshealth,
 :yummly,
 :youtube
]

v_lags=[ :html5_lag,
 :static_lag,
 :native_creative_lag,
 :video_lag,
 :content_lag,
 :contextual_lag,
 :direct_lag,
 :prospecting_lag,
 :retargeting_lag,
 :behavorial_lag,
 :native_lag,
 :predictive_lag,
 :thirdrdparty_lag,
 :pmp_lag,
 :buzzfeed_lag,
 :hulu_lag,
 :popsugarus_lag,
 :rachaelraymag_lag,
 :shape_lag,
 :blank_lag,
 :eatingwell_lag,
 :meredithcorporation_lag,
 :realsimple_lag,
 :turndsp_lag,
 :cookinglight_lag,
 :dataxu_lag,
 :allrecipes_lag,
 :amazon_lag,
 :myfitnesspal_lag,
 :womenshealth_lag,
 :yummly_lag,
 :youtube_lag
]

#dfx[v_lags]

#TEST
cols=vcat(v_lags,v_breaks, v_covar)
dtrain = DMatrix(convert(Array{Float32},dfx[cols]), label = net_price)
dtest = DMatrix(convert(Array{Float32},dfx[cols]))
dctrl = DMatrix(convert(Array{Float32},dfxc[cols]))
#param = ["eta" => 0.2, "max_depth" => 20, "objective" => "reg:linear", "silent" => 1]
xgb = xgboost(dtrain, 250, param = ["eta" => 0.2, "max_depth" => 20, "objective" => "reg:linear", "silent" => 1] )
ptest = XGBoost.predict(xgb, dtest)
pctrl = XGBoost.predict(xgb, dctrl)
mean(convert(Array{Float64},ptest))
mean(convert(Array{Float64},pctrl))

mean(ptest) - mean(pctrl)

dfx[dfx[:panid].==1321364138,:]
dfx[dfx[:net_price].!=0.0,:]
julia> mean(dfx[dfx[:gross_imps].==0.0,:net_price])
14.720302357882089
julia> mean(dfx[dfx[:gross_imps].!=0.0,:net_price])
0.0
julia> mean(dfx[dfx[:gross_imps].>0.0,:net_price])
0.0
julia> length(dfx[dfx[:gross_imps].>0.0,:net_price])
6123200
julia> length(dfx[dfx[:gross_imps].==0.0,:net_price])
1091488

# *********************************************************************

"""  """

results[    findin(results[:parameter],map(x->string(x),    setdiff(m[:finalvars],[:group])  ) )         ,[:parameter,:coef]]

findin(results[:parameter],map(x->string(x),    setdiff(m[:finalvars],[:group])  ) )

bst <- xgboost(data = train$data, label = train$label, max.depth = 4,
               num_parallel_tree = 1000, subsample = 0.5, colsample_bytree =0.5,
               nround = 1
               , objective = "binary:logistic")
"""  """

function genData(dfx::DataFrame,ycol::Symbol,z::Symbol[])
    #dfy = deepcopy
    y = convert(Array{Float32},dfx[ycol] )
    x=convert(Array{Float32}, dfx[setdiff(names(dfx),[ycol])])
    return DMatrix(x, label = y), y, x
end
dtrain,y,x = genData(dfd,:net_price)

num_round = 250
param = ["eta" => 0.2, "max_depth" => 20, "objective" => "reg:linear", "silent" => 1]
xgb = xgboost(dtrain, num_round, param=param)
pred = XGBoost.predict(xgb,dtrain)

dtest,ty,tx = genData(dfdc,:net_price)
DMatrix(tx)

pred_ctrl = XGBoost.predict(xgb,dtest)

# ------------------ END XGB --------------------------

==================  TREES WORKING ==================== examples https://github.com/bensadeghi/DecisionTree.jl
using DecisionTree
y=dfd[:dol_per_trip_pre_p1]
x=dfd[setdiff(names(dfd),[:dol_per_trip_pre_p1])]

model = build_tree(Array(y), Array(x), 5)
or
m2 = build_forest(Array(y), Array(x), 2, 10, 5, 0.7)

apply_tree(model, Array(x))

xtest=x[1:1000,:]
xtrain=x[1000:end,:]
ytest=y[1:1000]
ytrain=y[1000:end]
model = build_tree(Array(ytrain), Array(xtrain), 5)
res = apply_tree(model, Array(xtest))
map((x,y)->x-y, res,ytest)

r2 = nfoldCV_tree(Array(ytrain), Array(xtrain), 3, 5)

m2 = build_forest(Array(ytrain), Array(xtrain), 2, 10, 5, 0.7)
res = apply_tree(m2, Array(xtest))

# Symbol in expression
a=:xyz
QuoteNode(a)
:(5 - $t)
:( :a in $( :(:a + :b) ) )

s=:ssss
:( :a in $( :(:($s) + $s) ) )

a=:(:aaaa)
b=:(:bbb)
:( :a in $( :($a + $b) ) )

function t(s::Symbol)
           println(s," ~~ ",typeof(s))
       end
t(:xyz)
:(t(:xyz))
ex = :(t(:xyz))

:(t(:xyz))

v=:xyz
:(t( :($v) ))

ai=:aaa
bi=:bbb
a= QuoteNode(a)
b= QuoteNode(b)
:( :a in $( :($a + $b) ) )

macro tst(v, exp)
    println("v : ",v," ~~ ",typeof(v))
    println("exp : ",exp," ~~ ",typeof(exp))
    return exp
end

vi=:v1
@tst vi :(println("hello"))

#---------

#----
@everywhere function genx(p::Int64)
        brks=[:html5
         ,:static
         ,:native_creative
         ,:video
         ,:content
         ,:contextual
         ,:direct
         ,:prospecting
          ,:retargeting
          ,:behavorial
          ,:native
          ,:predictive
          ,:thirdrdparty
          ,:pmp
          ,:buzzfeed
          ,:hulu
          ,:popsugarus
          ,:rachaelraymag
          ,:shape
          ,:blank
          ,:eatingwell
          ,:meredithcorporation
          ,:realsimple
          ,:turndsp
          ,:cookinglight
          ,:dataxu
          ,:allrecipes
          ,:amazon
          ,:myfitnesspal
          ,:womenshealth
          ,:yummly
         ,:youtube
        ]
    dfx = readtable("/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks/"*string(p)*".csv",header=true)
    println("xyz : ",p)
end
@spawnat 5  genx(64)

p=1
root="/mapr/mapr04p/analytics0001/analytic_users/Models/RF/chunks"
dfx = readtable(root*"/"*string(p)*".csv",header=true)
y=:net_price
#----

  brk=:youtube
  lag=:youtube_lag
  dfx[lag]=0
  for row in eachrow(dfx)
        #println(r[:video])
        row[lag] = sum(dfx[(dfx[:panid].==row[:panid])&(dfx[:datenum].<=row[:datenum] ) & ( dfx[:datenum].>= row[:datenumlag]  ) ,brk])
  end

dfx[(dfx[brk].!=dfx[lag]) ,[:panid,:date,:datelag,:video_lag,:datenum,:datenumlag,brk]]
t = dfx[5634,[:panid,:date,:datelag,lag,:datenum,:datenumlag,brk]]      #...... 1321364407 │ "2016-07-18" │ "2016-07-11" │ 2         │ 736163  │ 736156
t = dfx[(dfx[:datenum].==736174)&(dfx[:panid].==1338983878),[:panid,:date,:datelag,lag,:datenum,:datenumlag,brk]]
dfx[(dfx[:panid].==t[:panid][1]) ,brk]
dfx[(dfx[:panid].==t[:panid][1])&(dfx[:datenum].<=t[:datenum][1] ) & (dfx[:datenum].>= t[:datenumlag][1]  ) ,brk]
sum(dfx[(dfx[:panid].==t[:panid][1])&(dfx[:datenum].<=t[:datenum][1]) & (  dfx[:datenum].>= t[:datenumlag][1]  ) ,brk])

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

dfd[1,:panid]
dfd[1,[:panid,:date,:datelag]]       == 1305262312 │ 2016-07-06 │ 2016-06-29
dfd[dfd[:panid].==1305262312,:]
sum(dfd[dfd[:panid].==1305262312,:video])

dfd[(dfd[:panid].==1305262312)&(dfd[:date].<=Date("2016-07-06",Dates.DateFormat("y-m-d")) ) & (  dfd[:date].>=Date("2016-06-29",Dates.DateFormat("y-m-d"))  ) ,:]

dfd[1,[:panid,:date,:datelag,:datenum,:datenumlag]]
sum(dfd[(dfd[:panid].==1305262312)&(dfd[:datenum].<=736151 ) & (  dfd[:datenum].>= 736144  ) ,:video])

dfd[:video_lag]=0
for row in eachrow(dfd)
    #println(r[:video])
    row[:video_lag] = sum(dfd[(dfd[:panid].==row[:panid])&(dfd[:datenum].<=row[:datenum] ) & (  dfd[:datenum].>= row[:datenumlag]  ) ,:video])
end

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

#hhlen=length(hhlst)
hhlen = Int(floor(length(hhlst)/10000))
for hhcnt in 0:hhlen
    lb=(hhcnt*50000)+1
    if hhcnt < hhlen
        ub=(hhcnt+1)*50000
        println("chuncking : (",lb,":",ub,")")
        #dfd[findin(dfd[:panid],panids) ,:iso] = true
        dfd[ findin(dfd[:panid],hhlst[lb:ub])  ,:chunks] = hhcnt+1
    else
        println("chuncking : (",lb,":end)")
        #dfd[ findin(dfd[:panid],hhlst[lb:end])  ,:chunks] = hhcnt+1
        dfd[dfd[:chunks].<0,:chunks]=hhcnt+1
    end
end

#dfd[:dateX] = map(x->    Date(string(x),Dates.DateFormat("yyyymmdd"))   ,dfd[:date])
#dfd[:dateX] = convert(Array{Date}, dfd[:dateX])

#lag_cnt=Dict(:one=>1,:four=>4,:eight=>8,:sixteen=>16)

for brk in [:placement, :creative, :publisher]
    for lvl in unique(dfd[brk])
        for lag in [1,4] #[:one,:four,:eight]
            col=Symbol(string(brk)*"_"*replace(string(lvl)," ","_")*"_"*string(lag))
            println(col)
            dfd[col]=0
        end
    end
end

for chunk in 1:2   #maximum(dfd[:chunks])
    for row in eachrow(dfd[dfd[:chunks].==chunk,:])
        println(row[:hh_id]," ~~~ ",row[:date])
    end
end

x=by(dfd,[:date,:hh_id], df-> sum(     df[ ()&()  :imps]       ))
x=by(dfd,[:date,:hh_id], df-> sum(  df[ ()&()  :imps]       ))

y=by(dfd, [:date,:hh_id,:publisher], nrow)

by(dfd, [:date,:hh_id]) do df
    DataFrame(m = mean(df[:PetalLength]), s² = var(df[:PetalLength]))
end

gb = groupby(dfd, [:date])

by(iris, :Species) do df
    DataFrame(m = mean(df[:PetalLength]), s² = var(df[:PetalLength]))
end

#=========================== END =============================

x=by(dfd,:hh_id, df-> sum(df[:imps] ))
or : x=by(dfd,:hh_id, df-> sum(df[:gross_imps] ))
x[x[:x1].>0,:]
x[x[:hh_id].==1000120003,:]

function oliers(dfd::DataFrame, k::Symbol, c::Symbol)
    m = median(dfd[c])
    #dfd[c_med1] = abs(dfd[c]-m)
    #MAD=median(dfd[c_med1])
    MAD = median(abs(dfd[c]-m))
    dfd[c_med2] =  ((dfd[c]-m)*0.6745) / MAD

    dfd_zsc = dfd[abs(dfd[c_med2]) .< 3.5 ,:]
    df_in_pout =  join(df_in, dfd_zsc, on = [ k, k ], kind = :inner)
end
oliers(dfd,:household_id, :val_imps)

# by(dfd, :, df -> sum(df[:PetalLength]))
function Pre_out(df_in::DataFram)
    df_cat_pre = df_in[df_in[:Buyer_Pre_P0] .==1 , [:Prd_0_Net_Pr_PRE,:experian_id]]
    median_df = median(df_cat_pre[:Prd_0_Net_Pr_PRE])
    df_cat_pre[:Prd_0_Net_Pr_PRE_med1] = abs(df_cat_pre[:Prd_0_Net_Pr_PRE]-median_df)
    MAD=median(df_cat_pre[:Prd_0_Net_Pr_PRE_med1])
    df_cat_pre[:Prd_0_Net_Pr_PRE_med2] = (0.6745*(df_cat_pre[:Prd_0_Net_Pr_PRE]-median_df))/MAD
    df_cat_pre_zsc = df_cat_pre[abs(df_cat_pre[:Prd_0_Net_Pr_PRE_med2]) .< 3.5 ,:]
    df_in_pout =  join(df_in, df_cat_pre_zsc, on = [ :experian_id, :experian_id ], kind = :inner)
end
Pre_out(df_in)

# --- parallel RF

function build_forest(labels, features, nsubfeatures, ntrees, ncpu=1)
                  if ncpu > nprocs()
                       addprocs(ncpu - nprocs())
                  end
                  Nlabels = length(labels)
                  Nsamples = int(0.7 * Nlabels)
                  forest = @parallel (vcat) for i in [1:ntrees]
                      inds = rand(1:Nlabels, Nsamples)
                      build_tree(labels[inds], features[inds,:], nsubfeatures)
                  end
                  return [forest]
              end

+++++++++Distribute Macro ++++++++++++++++++++++++++++
function sync_add(r)
    spawns = get(task_local_storage(), :SPAWNS, ())
    if spawns !== ()
        push!(spawns[1], r)
        if isa(r, Task)
            tls_r = get_task_tls(r)
            tls_r[:SUPPRESS_EXCEPTION_PRINTING] = true
        end
    end
    r
end

spawnat(p, thunk) = sync_add(remotecall(thunk, p))

macro spawnat(p, expr)
    expr = localize_vars(esc(:(()->($expr))), false)
    :(spawnat($(esc(p)), $expr))
end

macro run(p, expr)
    expr = localize_vars(esc(:(()->($expr))), false)
    :(spawnat($(esc(p)), $expr))
end
+++++++++++++++++++++++++++++++++++++

================== XGBoost ============================================= https://www.kaggle.com/wacaxx/rossmann-store-sales/julia-xgboost-starter-code
using BinDeps, DataFrames, XGBoost
#using Dates

y=convert(Array{Float32},dfd[:dol_per_trip_pre_p1] )
x=dfd[setdiff(names(dfd),[:dol_per_trip_pre_p1])]
#Define target
#y = convert(Array{Float32}, train[:Sales])

#Transform data to XGBoost matrices
#trainArray = convert(Array{Float32},  x[:, vcat(numericalColumns, categoricalColumns)])
#testArray = convert(Array{Float32}, test[:, vcat(numericalColumns, categoricalColumns)])
#dtrain = DMatrix(trainArray, label = costTrainingLog)
#dtest = DMatrix(testArray)

num_round = 250
param = ["eta" => 0.2, "max_depth" => 20, "objective" => "reg:linear", "silent" => 1]

XGBoostModel = xgboost(dtrain, num_round, param = param)

#Predictions using test data
preds = predict(XGBoostModel, dtest)
#Round to zero closed stores
preds[closedStoreIdx] = 0

#Write Results
sampleSubmission = DataFrame(Id = tesdIdx, Sales = preds)

==================  TREES WORKING ==================== examples https://github.com/bensadeghi/DecisionTree.jl
using DecisionTree
y=dfd[:dol_per_trip_pre_p1]
x=dfd[setdiff(names(dfd),[:dol_per_trip_pre_p1])]

model = build_tree(Array(y), Array(x), 5)
or
m2 = build_forest(Array(y), Array(x), 2, 10, 5, 0.7)

apply_tree(model, Array(x))

xtest=x[1:1000,:]
xtrain=x[1000:end,:]
ytest=y[1:1000]
ytrain=y[1000:end]
model = build_tree(Array(ytrain), Array(xtrain), 5)
res = apply_tree(model, Array(xtest))
map((x,y)->x-y, res,ytest)

r2 = nfoldCV_tree(Array(ytrain), Array(xtrain), 3, 5)

m2 = build_forest(Array(ytrain), Array(xtrain), 2, 10, 5, 0.7)
res = apply_tree(m2, Array(xtest))

====================== END ===========================

df_in[:Dol_per_Trip_PRE_P1]

setdiff(names(df_in),:Dol_per_Trip_PRE_P1)

y=df_in[:Dol_per_Trip_PRE_P1]
x=df_in[setdiff(names(df_in),[:Dol_per_Trip_PRE_P1])]
#model = build_forest(y,x, 20, 50, 1.0)
model = build_tree(y, Array(x),20, 50, 1.0);

using DecisionTree

model = build_forest(df_in[:Dol_per_Trip_PRE_P1], df_in[[setdiff(names(df_in),:Dol_per_Trip_PRE_P1)]], 20, 50, 1.0)

x=Array(x[setdiff(names(x),[:core_based_statistical_areas,:person_1_birth_year_and_month])])

model = build_tree(Array(y), Array(x),20, 50, 1.0);

========================================================
#https://github.com/dmlc/XGBoost.jl/blob/master/demo/basic_walkthrough.jl

using XGBoost

function readlibsvm(fname::ASCIIString, shape)
    dmx = zeros(Float32, shape)
    label = Float32[]
    fi = open(fname, "r")
    cnt = 1
    for line in eachline(fi)
        line = split(line, " ")
        push!(label, float(line[1]))
        line = line[2:end]
        for itm in line
            itm = split(itm, ":")
            dmx[cnt, int(itm[1]) + 1] = float(int(itm[2]))
        end
        cnt += 1
    end
    close(fi)
    return (dmx, label)
end

train_X, train_Y = readlibsvm("/home/rmadmin/.julia/v0.4/XGBoost/data/agaricus.txt.train", (6513, 126))
test_X, test_Y = readlibsvm("/home/rmadmin/.julia/v0.4/XGBoost/data/agaricus.txt.test", (1611, 126))

=============================
using DecisionTree

#Train random forest with
#20 for number of features chosen at each random split,
#50 for number of trees,
#and 1.0 for ratio of subsampling.
model = build_forest(yTrain, xTrain, 20, 50, 1.0)

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

function loadDF()
    cd("/media/u01/analytics/scoring/CDW5_792/")
    df_data = readtable("csv_final_cdw5_792.csv",header=false);
    df_h = readtable("Headers_cdw5.csv",header=false);
    names!(df_data, convert(Array{Symbol}, df_h[:x1]) )
end
df_in=loadDF()

# --===========================================================================================================================================
# --===================================== SQL =================================================================================================
# --===========================================================================================================================================

select count(a.household_id), count(distinct a.household_id), b.operating_company
from wh_supplmental.household_dependent_id_map a
join wh_fsp.partitioned_fsp_wkly_fact b
     on a.dependent_id = b.card_num
     and UPPER(trim(b.Operating_company)) = trim(UPPER(case when a.retailer='FOODLION' then 'FOOD LION' else a.retailer end ))
where a.current_state = 'MATCHED'
  and a.household_type = 'EXPERIAN'
  and a.dependent_type = 'FSP_CARD'
  and TM_DIM_KEY_WEEK>=1909
group by b.operating_company
"""
149243  17012   FREDS
196025906       7033205 RITEAID
330798388       994861  WEGMANS
38008898        221081  KEYFOOD
6020735400      15722602        KROGER
16025680        100830  SUPERVALU
302291297       2137143 BJS
"""

hive -e 'SELECT distinct household_id FROM WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP where household_type = "EXPERIAN" and dependent_type = "FSP_CARD" and  dependent_id is not null' | sed 's/[\t]/,/g' > /mapr/mapr04p/analytics0001/analytic_users/mddak/fsp_ids_retailer2.csv

set hive.cli.print.header=true;
SELECT retailer, count(distinct household_id) as ids FROM WH_SUPPLMENTAL.HOUSEHOLD_DEPENDENT_ID_MAP where household_type = "EXPERIAN" and dependent_type = "FSP_CARD" and  dependent_id is not null GROUP BY retailer;

hadoop fs -put /mapr/mapr04p/analytics0001/analytic_users/mddak/Experian/IRI_Key_Food_Shipment.txt /externaldata01/prd/experian/xwalk/raw/

set hive.cli.print.header=true; select * from wh_fsp.partitioned_fsp_wkly_fact limit 10;
hive -e 'set hive.cli.print.header=true; select * from MEDIA_LIFT_CAMPAIGN.lift_product WHERE registration_request_id=931 limit 10' | sed 's/[\t]/,/g' > /mapr/mapr04p/analytics0001/analytic_users/mddak/lift_product_example10.csv
set hive.cli.print.header=true; select * from MEDIA_LIFT_CAMPAIGN.lift_product_item_filter_map WHERE registration_request_id=931 limit 10;
set hive.cli.print.header=true; select * from MEDIA_LIFT_CAMPAIGN.POS_WKLY_FACT POS WHERE registration_request_id=931 limit 10;

set hive.cli.print.header=true; select * from Jennieo6.Kantar_FSP_POS_Jennieo6_931 limit 10;

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
"""
```
KrusellSmith{T} <: AbstractModel{T}
```

The `KrusellSmith` type defines the structure of the simple New Keynesian DSGE
model described in 'Bayesian Estimation of DSGE Models' by Sungbae An and Frank
Schorfheide.

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"an_schorfheide\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.

"""
mutable struct KrusellSmith{T} <: AbstractModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values
    grids::OrderedDict{Symbol,Grid}                        # quadrature grids
    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    state_variables::Vector{Symbol}                        # Vector of symbols of the state variables
    jump_variables::Vector{Symbol}                         # Vector of symbols of the jump variables
    normalized_model_states::Vector{Symbol}                # All of the distributional model
                                                           # state variables that need to be normalized
    endogenous_states_unnormalized::OrderedDict{Symbol,UnitRange} # Vector of unnormalized
                                                           # ranges of indices
    endogenous_states::OrderedDict{Symbol,UnitRange}       # Vector of ranges corresponding
                                                           # to normalized (post Klein solution) indices
    exogenous_shocks::OrderedDict{Symbol,Int}              #
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,UnitRange}  #
    observables::OrderedDict{Symbol,Int}                   #

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
end

description(m::KrusellSmith) = "KrusellSmith, $(m.subspec)"

"""
`init_model_indices!(m::KrusellSmith)`

Arguments:
`m:: KrusellSmith`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::KrusellSmith)
    # Endogenous states
    endogenous_states = collect([
    # These states corresp. to the following in the original notation
    #    LMP,    LELLP,  KKP,   ZP,    MP,    ELLP,
        :μ′_t1, :l′_t1, :K′_t, :z′_t, :μ′_t, :l′_t])

    # Exogenous shocks
    exogenous_shocks = collect([:z_sh])

    # Equilibrium conditions
    equilibrium_conditions = collect([
        :eq_euler, :eq_kolmogorov_fwd, :eq_μ, :eq_l, :eq_K_law_of_motion, :eq_TFP])

    # Observables
    observables = keys(m.observable_mappings)

    ########################################################################################
    # Setting indices of endogenous_states and equilibrium conditions manually for now
    nw = get_setting(m, :nw)
    endo = m.endogenous_states_unnormalized
    eqconds = m.equilibrium_conditions

    # State variables
    endo[:μ′_t1] = 1:nw
    endo[:l′_t1] = nw+1:2*nw
    endo[:K′_t]  = 2*nw+1:2*nw+1
    endo[:z′_t]  = 2*nw+2:2*nw+2

    # Jump variables
    endo[:μ′_t]  = 2*nw+3:3*nw+2
    endo[:l′_t]  = 3*nw+3:4*nw+2

    eqconds[:eq_euler]              = 1:nw
    eqconds[:eq_kolmogorov_fwd]     = nw+1:2*nw
    eqconds[:eq_μ]                  = 2*nw+1:3*nw
    eqconds[:eq_l]                  = 3*nw+1:4*nw
    eqconds[:eq_K_law_of_motion]    = 4*nw+1:4*nw+1
    eqconds[:eq_TFP]                = 4*nw+2:4*nw+2
    ########################################################################################

    m.normalized_model_states = [:μ′_t1, :μ′_t]
    m.endogenous_states = deepcopy(endo)
    m.state_variables = m.endogenous_states.keys[get_setting(m, :state_indices)]
    m.jump_variables = m.endogenous_states.keys[get_setting(m, :jump_indices)]

    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
end

function KrusellSmith(subspec::String="ss0";
                      custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                      testing = false)

    # Model-specific specifications
    spec               = "KrusellSmith"
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = KrusellSmith{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(),
            # grids and keys
            OrderedDict{Symbol,Grid}(), OrderedDict{Symbol,Int}(),

            # normalized_model_states, state_inds, jump_inds
            Vector{Symbol}(), Vector{Symbol}(), Vector{Symbol}(),

            # model indices
            # endogenous states unnormalized, endogenous states normalized
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,UnitRange}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}())

    # Set settings
    model_settings!(m)
    # default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable transformations
    init_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    # Initialize grids
    init_grids!(m)

    init_model_indices!(m)
    init_subspec!(m)

    # Solve for the steady state
    steadystate!(m)

    # So that the indices of m.endogenous_states reflect the normalization
    normalize_model_state_indices!(m)

    return m
end

"""
```
init_parameters!(m::KrusellSmith)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::KrusellSmith)
    # Initialize parameters
    m <= parameter(:β, 0.95, fixed = true,
                   description = "β: Discount rate.", tex_label = "\\beta")
    m <= parameter(:γ, 3.0, fixed = true,
                   description = "γ: CRRA Parameter.", tex_label = "\\gamma")
    m <= parameter(:α, 1.0/3.0, fixed = true,
                   description = "α: Capital Share.", tex_label = "\\alpha")
    m <= parameter(:δ, 0.2, fixed = true,
                   description = "δ: Depreciation.", tex_label = "\\delta")
    m <= parameter(:ρ_z, 0.95, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), BetaAlt(0.5, 0.2), fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:σ_z, sqrt(.007), (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{z}")
    m <= parameter(:μ_s, 0., fixed = true, description = "μ_s: Mu of log normal in income",
                   tex_label = "\\mu_s")
    m <= parameter(:σ_s, 0.5, fixed = true,
                   description = "σ_s: Sigma of log normal in income",
                   tex_label = "\\sigma_s")

    m <= parameter(:e_y, 1e-3, fixed = true, description = "e_y: Measurement error on GDP",
                   tex_label = "e_y")

    # Setting steady-state parameters
    nw = get_setting(m, :nw)

    m <= SteadyStateParameterGrid(:lstar, fill(NaN, nw), description = "Steady-state expected discounted
                                  marginal utility of consumption", tex_label = "l_*")
    m <= SteadyStateParameterGrid(:cstar, fill(NaN, nw), description = "Steady-state consumption",
                                  tex_label = "c_*")
    m <= SteadyStateParameterGrid(:μstar, fill(NaN, nw), description = "Steady-state cross-sectional
                                  density of cash on hand", tex_label = "\\mu_*")
    m <= SteadyStateParameter(:Kstar, NaN, description = "Steady-state capital", tex_label = "K_*")
    m <= SteadyStateParameter(:Lstar, NaN, description = "Steady-state labor", tex_label = "L_*")
    m <= SteadyStateParameter(:Wstar, NaN, description = "Steady-state wages", tex_label =
                              "W_*")
    m <= SteadyStateParameter(:Rstar, NaN, description = "Steady-state net-return on
                              capital", tex_label = "R_*")
    m <= SteadyStateParameterGrid(:KFstar, fill(NaN, (nw, nw)), description = "Steady-state Kolmogorov
                                  Forward Equation", tex_label = "KF_*")
end

"""
```
init_grids!(m::KrusellSmith)
```
"""
function init_grids!(m::KrusellSmith)
    wscale  = get_setting(m, :wscale)
    wlo     = get_setting(m, :wlo)
    whi     = get_setting(m, :whi)
    nw      = get_setting(m, :nw)

    slo     = get_setting(m, :slo)
    shi     = get_setting(m, :shi)
    ns      = get_setting(m, :ns)

    grids = OrderedDict()

    grids[:wgrid] = Grid(uniform_quadrature(wscale), wlo, whi, nw, scale = wscale)
    grids[:sgrid] = Grid(curtis_clenshaw_quadrature(2), slo, shi, ns)

    m.grids = grids
end

function model_settings!(m::KrusellSmith)
    default_settings!(m)

    # Defaults
    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), DSGE_DATE_FORMAT)
    m <= Setting(:data_vintage, vint, true, "vint", "Data vintage")

    saveroot = normpath(joinpath(dirname(@__FILE__), "../../../../","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "../../../../","save","input_data"))

    m <= Setting(:saveroot, saveroot, "Root of data directory structure")
    m <= Setting(:dataroot, datapath, "Input data directory path")

    # Solution method
    m <= Setting(:solution_method, :klein)

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0,
                 "Number of anticipated policy shocks")
    # m <= Setting(:n_anticipated_shocks_padding, 20,
                 # "Padding for anticipated policy shocks")

    # Number of states and jumps
    m <= Setting(:normalize_distr_variables, true, "Whether or not to perform the
                 normalization of the distributional states in the Klein solution step")
    m <= Setting(:n_predetermined_variables, 0, "Number of predetermined variables after
                 removing the densities. Calculated with the Jacobian normalization")

    m <= Setting(:state_indices, 1:4, "Which indices of m.endogenous_states correspond to
                 backward looking state variables")
    m <= Setting(:jump_indices, 5:6, "Which indices of m.endogenous_states correspond to jump
                 variables")

    # Note, these settings assume normalization.
    m <= Setting(:n_degrees_of_freedom_removed, 1, "Number of degrees of freedom from the
                 distributional variables to remove.")
    n_dof_removed = get_setting(m, :n_degrees_of_freedom_removed)
    m <= Setting(:n_backward_looking_distributional_vars, 1, "Number of state variables that are
                 distributional variables.")
    n_backlook_dists = get_setting(m, :n_backward_looking_distributional_vars)
    m <= Setting(:backward_looking_states_normalization_factor,
                 n_dof_removed*n_backlook_dists, "The number of dimensions removed from the
                 backward looking state variables for the normalization.")
    m <= Setting(:n_backward_looking_states, 162 - get_setting(m, :backward_looking_states_normalization_factor),
                 "Number of state variables, in the true sense (fully
                  backward looking) accounting for the discretization across the grid.")
    m <= Setting(:n_jump_distributional_vars, 1, "Number of jump variables that are distributional
                 variables.")
    n_jump_dists = get_setting(m, :n_jump_distributional_vars)
    m <= Setting(:jumps_normalization_factor,
                 n_dof_removed*n_jump_dists, "The number of dimensions removed from the
                 jump variables for the normalization.")
    m <= Setting(:n_jumps, 160 - get_setting(m, :jumps_normalization_factor),
                 "Number of jump variables (forward looking) accounting for
                  the discretization across the grid")

    m <= Setting(:n_model_states, get_setting(m, :n_backward_looking_states) + get_setting(m, :n_jumps),
                 "Number of 'states' in the state space model. Because backward and forward
                 looking variables need to be explicitly tracked for the Klein solution
                 method, we have n_states and n_jumps")

    # Mollifier setting parameters
    m <= Setting(:In, 0.443993816237631, "Normalizing constant for the mollifier")
    m <= Setting(:smoother, 1.0, "Scale of the molifier (how peaked the mollifier is)
                 smoother = 0 you get uniform dist")
    m <= Setting(:zlo, 0.0, "Second income shock to mollify actual income")
    m <= Setting(:zhi, 2.0, "Upper bound on this process")

    # w: Cash on Hand Grid Setup
    m <= Setting(:nw, 80, "Cash on hand distribution grid points")
    m <= Setting(:wlo, 0.5, "Lower bound on cash on hand")
    m <= Setting(:whi, 15.0, "Upper Bound on cash on hand")
    m <= Setting(:wscale, get_setting(m, :whi) - get_setting(m, :wlo), "Size of the wgrid")

    # s: Skill Distribution/ "Units of effective labor" Grid Setup
    m <= Setting(:ns, 50, "Skill distribution grid points")
    m <= Setting(:slo, 0.5, "Lower bound on skill grid")
    m <= Setting(:shi, 1.5, "Upper bound on skill")
    m <= Setting(:sscale, (get_setting(m, :shi) - get_setting(m, :slo))/2.0, "Size of our s
                 grid, relative to size of [-1, 1]")
end

ext: jl


text:
@testset "minimal_sum_test" begin
    # sum_nearest {1.0, 2.0, 3.0} = 6.0
    # sum_nearest {1.0, 2.0, NaN, 3.0} = NaN
    # sum_nearest {1.0, -infinity, 2.0, infinity, 3.0} = NaN
end

@testset "minimal_sum_abs_test" begin
    # sum_abs_nearest {1.0, -2.0, 3.0} = 6.0
    # sum_abs_nearest {1.0, -2.0, NaN, 3.0} = NaN
    # sum_abs_nearest {1.0, -infinity, 2.0, infinity, 3.0} = infinity
end

@testset "minimal_sum_sqr_test" begin
    # sum_sqr_nearest {1.0, 2.0, 3.0} = 14.0
    # sum_sqr_nearest {1.0, 2.0, NaN, 3.0} = NaN
    # sum_sqr_nearest {1.0, -infinity, 2.0, infinity, 3.0} = infinity
end

@testset "minimal_dot_test" begin
    # dot_nearest {1.0, 2.0, 3.0} {1.0, 2.0, 3.0} = 14.0
    # dot_nearest {0x10000000000001p0, 0x1p104} {0x0fffffffffffffp0, -1.0} = -1.0
    # dot_nearest {1.0, 2.0, NaN, 3.0} {1.0, 2.0, 3.0, 4.0} = NaN
    # dot_nearest {1.0, 2.0, 3.0, 4.0} {1.0, 2.0, NaN, 3.0} = NaN
    # dot_nearest {1.0, 2.0, 0.0, 4.0} {1.0, 2.0, infinity, 3.0} = NaN
    # dot_nearest {1.0, 2.0, -infinity, 4.0} {1.0, 2.0, 0.0, 3.0} = NaN
end


text:
"""
Description

        The yahoo() method is a wrapper for downloading historical stock prices from Yahoo.

Usage

        AAPL = yahoo("AAPL)
        SPX = yahoo()

Method Signature(s)

      yahoo(stock::ASCIIString="^GSPC")

Details

      The yahoo method takes a stock name in the form of a string and returns a TimeSeries.TimeArray data structure
      corresponding to the Yahoo Finance ticker. With no argument, the default historical time series is the S&P 500.

References

      http://www.finance.yahoo.com

See Also

      fred() which accesses the St. Louis Federal Reserve financial and economic data sets.

"""

function yahoo(stock::ASCIIString="^GSPC")
    resp = get("http://ichart.yahoo.com/table.csv?s=$stock")
    timearray(resp, m=stock)
end


text:
####################
#File Name: barcodeCounterFuncts.jl
#Author:    Henry C. Wong (wongh@wustl.edu henrycwong@mst.edu cywongx@gmail.com)
#Desc:      Helper file with helper functions for barcodeCounter.jl
####################

#sequence to match
#const MATCH_SEQ = uppercase("AGGCATGGACGAGCTCTATAAGTAATCTAGA");
#const BARCODE_LENGTH = 12;
#const REGEX_TAG = string(".{",BARCODE_LENGTH,"}");

function get_match_seq()
  readArr = ["",""];
  open(ARGS[3]) do f
    line = 1;
    while ! eof(f)
      readArr[line] = readline(f);
      line += 1;
    end
  end
  seq = readArr[1];
  return (seq);
end

##########
#Name: createRegex
#Desc: Forms the Regex object to filter out sequences
#Pre : Tag contains a valid Regex value
#Post: Returns a regex object
##########
function createRegex(seq::String,tag::String)::Regex
  re = Regex(string(seq,tag));
  return re;
end

##########
#Name: readBarcodeFile
#Desc: Reads the barcode tsv or text (in a tsv format) file and converts the input
#      into a managable 2d array
#Pre : barcodeFile is a tsv or tsv-like file that contains two columns. The first
#      column being a list of barcodes, the second being a list of names(strings)
#Post: Returns a 2d array with 2 columns
##########
function readBarcodeFile(barcodeFile::String)::Array{String}
  grid = readdlm(barcodeFile,String);
  return grid;
end

##########
#Name: readFastqFile
#Desc: Reads the passed in fastq file and converts it into a managable FASTQ
#      format provided by the FASTX library
#Pre : file in fileName is a valid fastq file
#Post: Returns the contesnts of the passed in fastq file of type FASTQ
#      Specifically FASTX.FASTQ.Reader{TranscodingStreams.TranscodingStream
#        {TranscodingStreams.Noop,IOStream}}
##########
function readFastqFile(fileName)
  reader = FASTQ.Reader(open(fileName,"r"));
  return reader;
end

##########
#Name: fillDict
#Desc: Fills count and name dictionary passed in with the fastq file passed in
#Pre : countDict and nameDict are the same length and both have the same keys.
#      fileName is a valid fastq file
#Post: Both countDict will be filled with a number of times a barcode has
#      appeared. Both countDict and nameDict will still have equal number
#      of entries. However, the entry count for both may chance.
##########
function fillDict(fileName::String,countDict::Dict,nameDict::Dict)::Nothing
  i = 0;
  j = 0;
  reader = readFastqFile(fileName);
  seq = get_match_seq();
  barcodeLength = length(collect(countDict)[1][1]);
  re = createRegex(seq,string(".{",barcodeLength,"}"));
  for record in reader
    i += 1;
    seq2 = sequence(record);
    m = match(re,string(seq2));
    if !(isnothing(m))
      mString = m.match;
      barcode = mString[end-11:end];
      if barcode in keys(countDict)
        countDict[barcode] += 1;
        j+=1;
       else
        #uncomment if you want non matching barcodes to be recorded
        #WARNING: Uncommenting may break script used to merge all counts together
        #push!(countDict,barcode => 1);
        #push!(nameDict,barcode => "NA");
      end
    end
  end
  print("File: ");
  println(fileName);
  print("Size: " );
  println(i);
  print("Hits: ");
  println(j);
  print("Hit Ratio(hits/size): ");
  println(j/i);
end

##########
#Name: syncDictsToArray
#Desc: Synchronizes both countDict and namesDict into a 2d 3-column array.
#Pre : countDict and namesDict have the same number of entries.
#Post: Returns a 3-column array.
##########
function syncDictsToArray(countDict::Dict,namesDict::Dict)::Array
  ks = collect(keys(countDict));
  vals = collect(values(countDict));
  newNames = fill(" ", length(ks));
  for i=1:length(ks)
    newNames[i] = nameDict[ks[i]];
  end
  return hcat(ks,hcat(newNames,vals));
end

##########
#Name: createOutputFileName
#Desc: creates a filename(string) based on given fastq file
#Pre : N/A
#Post: Returns a string
##########
function createOutputFileName(fileName::String)::String
  fm = fileName[1:end-6]; #remove .fastq extention
  fm = string(fm,".counts");
  return fm;
end

##########
#Name: writeResultsToFile
#Desc: Writes contentArray to a file named fileName
#Pre : N/A
#Post: File fileName contains the content of contentArray
##########
function writeResultsToFile(contentArray::Array,fileName::String)::Nothing
  writedlm(fileName,contentArray,"\t");
end

#grid = readBarcodeFile(ARGS[1]);
#barcodeArray = grid[:,1];
#names = grid[:,2];
#counts = fill(0,length(names));
#countDict = Dict(barcodeArray .=> counts);
#nameDict = Dict(barcodeArray .=> names);
#fillDict(ARGS[2],countDict,nameDict); #fill dictionaries with data from fastq file
#contentArray = syncDictsToArray(countDict,nameDict);
#sortedContentArr = sortslices(contentArray,dims=1,by=x->x[3],rev=true); #sort by count
#fileName = createOutputFileName(ARGS[2]);
#writeResultsToFile(sortedContentArr,fileName);

max_stars_repo_name: HenryCWong/genomic-barcode-counter

max_issues_repo_name: HenryCWong/genomic-barcode-counter

max_forks_repo_name: HenryCWong/genomic-barcode-counter


text: @test Foo.say_hello() == "hello"


text:
include(joinpath(dirname(@__FILE__), "initialization.jl"))
include(joinpath(dirname(@__FILE__), "updateP!.jl"))
include(joinpath(dirname(@__FILE__), "updateD!.jl"))
include(joinpath(dirname(@__FILE__), "updateA!.jl"))

function TrainDPL(Data::Matrix{Float64}, Label::Matrix{Int64}, DictSize::Int64, τ::Float64, λ::Float64, γ::Float64, DEMO=false)
    # This is the DPL training function

    # Initialize D and P, compute the inverse matrix used in Eq. (10), update A for one time
    DataMat, D, P, DataInvMat, A = initialization(Data, Label, DictSize, τ, λ, γ, DEMO)

    # Alternatively update P, D and A
    for i = 1:20
        updateP!(P, A, DataInvMat, DataMat, τ)
        updateD!(D, A, DataMat)
        updateA!(A, D, DataMat, P, τ, DictSize)
        if DEMO
            println("Training ... iteration $i finished.")
        end
    end

    EncoderMat = zeros(size(P[1], 1) * size(P, 1), size(P[1], 2))
    # Reorganize the P matrix to make the classification fast
    for i = 1:size(P, 1)
        @inbounds EncoderMat[(i - 1)*DictSize+1 : i*DictSize, :] = P[i]
    end

    D, EncoderMat
end


text:
# Test that Mocking can interact with functions defined at different scopes or externally

# Global scope
global_scope() = "org"

# The @mock macro is essentially a no-op
@test (@mock global_scope()) == global_scope()

# Create a patched version of function
global_patch = (@patch global_scope() = "alt")
apply(global_patch) do
    @test (@mock global_scope()) != global_scope()
end

# The @mock macro should return to the original behaviour
@test (@mock global_scope()) == global_scope()

# Local scope within a function
function test_function_scope()
    function_scope() = "org"
    @test (@mock function_scope()) == function_scope()

    patch = @patch function_scope() = "alt"
    apply(patch) do
        @test (@mock function_scope()) != function_scope()
    end

    @test (@mock function_scope()) == function_scope()
end

test_function_scope()

# Local scope within a let-block
let let_scope
    let_scope() = "org"
    @test (@mock let_scope()) == let_scope()

    patch = @patch let_scope() = "alt"
    apply(patch) do
        @test (@mock let_scope()) != let_scope()
    end

    @test (@mock let_scope()) == let_scope()
end

@testset "qualified function" begin
    # Define a new generic function named zero
    zero() = 0

    patch = @patch Base.zero(T::Type{<:Integer}) = one(T)
    apply(patch) do
        # Call alternative
        @test (@mock Base.zero(Int)) == 1

        # Call original function
        @test (@mock Base.zero(1)) == 0

        # Call new function named zero
        @test (@mock zero()) == 0
        @test_throws MethodError (@mock zero(1))
    end
end

@testset "constructor function" begin
    @test Vector() == []

    p = @patch Vector() = [1,2,3]
    apply(p) do
        # Call alternative
        @test (@mock Vector()) == [1,2,3]
        @test (@mock Base.Vector()) == [1,2,3]

        # Call original function
        @test (@mock Base.Vector{Any}()) == Any[]
    end
end

struct Coefficient
    x
end

(c::Coefficient)(y) = c.x * y

# https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects-1
@testset "function-like objects" begin
    c = Coefficient(2)

    p = @patch (c::Coefficient)(y::Integer) = c.x + y
    apply(p) do
        # Call alternative
        @test (@mock c(3)) == 5

        # Call original function
        @test (@mock c(3.0)) == 6.0
    end
end

@testset "missing target" begin
    # The function `f` is not defined so we cannot target it.
    @test !isdefined(@__MODULE__, :f)
    @test_throws UndefVarError (@patch f() = 1)
end


text:
using Random

# for bug with display
ENV["GKSwstype"] = "100"

example_files = filter(readdir(@__DIR__; join=true)) do f
    endswith(f, ".jl") && !occursin("large", f) && !occursin("result", f) && f != "activate.jl"
end

example_shuffle = randperm(length(example_files))

activate_file = joinpath(@__DIR__, "activate.jl")

for file in example_files[example_shuffle[1:2]]
    @info "Including example $file"
    instruction = """include("$activate_file"); include("$file")"""
    run(`julia -e $instruction`)
end


text:
using LinearAlgebra
using FermiCG
using Printf
using Test
using LinearMaps
using Arpack

using Profile

atoms = []
push!(atoms,Atom(1,"H",[0,0,0]))
push!(atoms,Atom(2,"H",[0,0,1]))
push!(atoms,Atom(3,"H",[0,0,2]))
push!(atoms,Atom(4,"H",[0,0,3]))
push!(atoms,Atom(5,"H",[0,0,4]))
push!(atoms,Atom(6,"H",[0,0,5]))
push!(atoms,Atom(7,"H",[0,0,6]))
push!(atoms,Atom(8,"H",[0,0,7]))
#push!(atoms,Atom(9,"H",[0,0,8]))
#push!(atoms,Atom(10,"H",[0,0,9]))
#push!(atoms,Atom(11,"H",[0,0,10]))
#push!(atoms,Atom(12,"H",[0,0,11]))
#basis = "6-31g"
basis = "sto-3g"

mol     = Molecule(0,1,atoms)
mf = FermiCG.pyscf_do_scf(mol,basis)
ints = FermiCG.pyscf_build_ints(mf.mol,mf.mo_coeff);
@time e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,4,4)
# @printf(" FCI Energy: %12.8f\n", e_fci)

norbs = size(ints.h1)[1]

problem = StringCI.FCIProblem(norbs, 4, 4)

display(problem)

v0 = zeros(problem.dim,1)
v0[1] = 1

e, v = StringCI.run_fci(ints, problem, v0=v0, precompute_ss=true)
#e, v = StringCI.run_fci(ints, problem, v0=v)
#
#ket_a = StringCI.DeterminantString(problem.no, problem.na)
#ket_b = StringCI.DeterminantString(problem.no, problem.nb)
#
#v = reshape(v,ket_a.max,ket_b.max)
#F = svd(v)
#N = 50
#v0 = F.U[:,1:N] * Diagonal(F.S[1:N]) * F.Vt[1:N,:]
#v0 = reshape(v0,ket_a.max*ket_b.max)
#println(v0'*v0)
#v0 = v0/norm(v0)
#e, v = StringCI.run_fci(ints, problem, v0=v0)

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

symmetric-indefinite linear system solver
solves linear system in naive.jl by first eliminating s and kap via the method
in naiveelim.jl and then eliminating tau via a procedure similar to that
described by S7.4 of
http://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf

3x3 nonsymmetric system in (x, y, z):
A'*y + G'*z = [xrhs, -c]
-A*x = [yrhs, -b]
(pr bar) -mu*H_k*G_k*x + z_k = [mu*H_k*zrhs_k + srhs_k, -mu*H_k*h_k]
(du bar) -G_k*x + mu*H_k*z_k = [zrhs_k + srhs_k, -h_k]

multiply pr bar constraint by (mu*H_k)^-1 to get 3x3 symmetric indefinite system
A'*y + G'*z = [xrhs, -c]
A*x = [-yrhs, b]
(pr bar) G_k*x - (mu*H_k)\z_k = [-zrhs_k - (mu*H_k)\srhs_k, h_k]
(du bar) G_k*x - mu*H_k*z_k = [-zrhs_k - srhs_k, h_k]

TODO not implemented: to avoid inverse hessian products, let for pr bar
w_k = (mu*H_k)\z_k (later recover z_k = mu*H_k*w_k)
to get 3x3 symmetric indefinite system
A'*y + sum_{pr bar} G_k'*mu*H_k*w_k + sum_{du bar} G_k'*z_k = [xrhs, -c]
A*x = [-yrhs, b]
(pr bar) mu*H_k*G_k*x - mu*H_k*w_k = [-mu*H_k*zrhs_k - srhs_k, mu*H_k*h_k]
(du bar) G_k*x - mu*H_k*z_k = [-zrhs_k - srhs_k, h_k]
=#

abstract type SymIndefSystemSolver{T <: Real} <: SystemSolver{T} end

function setup_rhs3(
    ::SymIndefSystemSolver{T},
    model::Models.Model{T},
    rhs::Point{T},
    sol::Point{T},
    rhs_sub::Point{T},
    ) where {T <: Real}
    @inbounds for (k, cone_k) in enumerate(model.cones)
        rhs_z_k = rhs.z_views[k]
        rhs_s_k = rhs.s_views[k]
        rhs_sub_z_k = rhs_sub.z_views[k]
        if Cones.use_dual_barrier(cone_k)
            # G_k*x - mu*H_k*z_k = [-zrhs_k - srhs_k, h_k]
            @. rhs_sub_z_k = -rhs_z_k - rhs_s_k
        else
            # G_k*x - (mu*H_k)\z_k = [-zrhs_k - (mu*H_k)\srhs_k, h_k]
            Cones.inv_hess_prod!(rhs_sub_z_k, rhs_s_k, cone_k)
            axpby!(-1, rhs_z_k, -1, rhs_sub_z_k)
        end
    end
    return
end

#=
direct sparse
=#

mutable struct SymIndefSparseSystemSolver{T <: Real} <: SymIndefSystemSolver{T}
    lhs_sub::SparseMatrixCSC{T}
    fact_cache::SparseSymCache{T}
    hess_idxs::Vector
    rhs_sub::Point{T}
    sol_sub::Point{T}
    sol_const::Point{T}
    rhs_const::Point{T}
    function SymIndefSparseSystemSolver{T}(;
        fact_cache::SparseSymCache{T} = SparseSymCache{T}(),
        ) where {T <: Real}
        syssolver = new{T}()
        syssolver.fact_cache = fact_cache
        return syssolver
    end
end

function load(
    syssolver::SymIndefSparseSystemSolver{T},
    solver::Solver{T},
    ) where {T <: Real}
    syssolver.fact_cache.analyzed = false
    model = solver.model
    (n, p, q) = (model.n, model.p, model.q)
    cones = model.cones
    cone_idxs = model.cone_idxs

    # form lower triangle of symmetric sparse LHS
    # without Hessians and inverse Hessians in z/z block
    # TODO check for inefficiency, may need to implement more manually
    spz(a, b) = spzeros(T, a, b)
    lhs_sub = hvcat((1, 2, 2),
        spz(n, n + p + q),
        model.A, spz(p, p + q),
        model.G, spz(q, p + q))
    @assert lhs_sub isa SparseMatrixCSC{T}
    dropzeros!(lhs_sub)
    (Is, Js, Vs) = findnz(lhs_sub)

    # add I, J, V for Hessians and inverse Hessians
    if isempty(cones)
        hess_nz_total = 0
    else
        hess_nz_total = sum(Cones.use_dual_barrier(cone_k) ?
            Cones.hess_nz_count_tril(cone_k) :
            Cones.inv_hess_nz_count_tril(cone_k) for cone_k in cones)
    end
    H_Is = Vector{Int}(undef, hess_nz_total)
    H_Js = Vector{Int}(undef, hess_nz_total)
    offset = 1
    y_start = n + p - 1
    for (cone_k, idxs_k) in zip(cones, cone_idxs)
        z_start_k = y_start + first(idxs_k)
        for j in 1:Cones.dimension(cone_k)
            nz_rows_kj = z_start_k .+ (Cones.use_dual_barrier(cone_k) ?
                Cones.hess_nz_idxs_col_tril(cone_k, j) :
                Cones.inv_hess_nz_idxs_col_tril(cone_k, j))
            len_kj = length(nz_rows_kj)
            IJV_idxs = offset:(offset + len_kj - 1)
            offset += len_kj
            @. H_Is[IJV_idxs] = nz_rows_kj
            @. H_Js[IJV_idxs] = z_start_k + j
        end
    end
    append!(Is, H_Is)
    append!(Js, H_Js)
    append!(Vs, ones(T, hess_nz_total))

    diag_pert = T(diag_min(syssolver.fact_cache))
    append!(Is, 1:(n + p))
    append!(Js, 1:(n + p))
    append!(Vs, fill(diag_pert, n))
    append!(Vs, fill(-diag_pert, p))

    # integer type supported by the sparse system solver library to be used
    Ti = int_type(syssolver.fact_cache)
    # convert integer types here rather than inside external wrappers
    Is = convert(Vector{Ti}, Is)
    Js = convert(Vector{Ti}, Js)
    dim = size(lhs_sub, 1)
    lhs_sub = syssolver.lhs_sub = sparse(Is, Js, Vs, dim, dim)

    # cache indices of nonzeros of Hessians in sparse LHS nonzeros vector
    syssolver.hess_idxs = [Vector{Union{UnitRange, Vector{Int}}}(
        undef, Cones.dimension(cone_k)) for cone_k in cones]
    for (k, cone_k) in enumerate(cones)
        cone_idxs_k = cone_idxs[k]
        z_start_k = y_start + first(cone_idxs_k)
        for j in 1:Cones.dimension(cone_k)
            col = z_start_k + j
            # get nonzero rows in the current column of the LHS
            col_idx_start = lhs_sub.colptr[col]
            nz_rows = lhs_sub.rowval[col_idx_start:(lhs_sub.colptr[col + 1] - 1)]
            # get nonzero rows in column j of the Hessian or inverse Hessian
            nz_hess_indices = (Cones.use_dual_barrier(cone_k) ?
                Cones.hess_nz_idxs_col_tril(cone_k, j) :
                Cones.inv_hess_nz_idxs_col_tril(cone_k, j))
            # get index corresponding to first nonzero element of the current col
            first_H = findfirst(isequal(z_start_k + first(nz_hess_indices)),
                nz_rows)
            # indices of nonzero values for cone k column j
            syssolver.hess_idxs[k][j] = (col_idx_start + first_H - 2) .+
                (1:length(nz_hess_indices))
        end
    end

    setup_point_sub(syssolver, model)

    return syssolver
end

function update_lhs(syssolver::SymIndefSparseSystemSolver, solver::Solver)
    for (k, cone_k) in enumerate(solver.model.cones)
        H_k = (Cones.use_dual_barrier(cone_k) ? Cones.hess(cone_k) :
            Cones.inv_hess(cone_k))
        for j in 1:Cones.dimension(cone_k)
            nz_rows = (Cones.use_dual_barrier(cone_k) ?
                Cones.hess_nz_idxs_col_tril(cone_k, j) :
                Cones.inv_hess_nz_idxs_col_tril(cone_k, j))
            @. @views syssolver.lhs_sub.nzval[syssolver.hess_idxs[k][j]] =
                -H_k[nz_rows, j]
        end
    end

    solver.time_upfact += @elapsed update_fact(syssolver.fact_cache,
        syssolver.lhs_sub)
    solve_subsystem3(syssolver, solver, syssolver.sol_const, syssolver.rhs_const)

    return syssolver
end

function solve_subsystem3(
    syssolver::SymIndefSparseSystemSolver,
    ::Solver,
    sol::Point,
    rhs::Point,
    )
    inv_prod(syssolver.fact_cache, sol.vec, syssolver.lhs_sub, rhs.vec)
    return sol
end

#=
direct dense
=#

mutable struct SymIndefDenseSystemSolver{T <: Real} <: SymIndefSystemSolver{T}
    lhs_sub::Symmetric{T, Matrix{T}}
    lhs_sub_fact::Symmetric{T, Matrix{T}}
    fact::Factorization{T}
    rhs_sub::Point{T}
    sol_sub::Point{T}
    sol_const::Point{T}
    rhs_const::Point{T}
    function SymIndefDenseSystemSolver{T}(;
        ) where {T <: Real}
        syssolver = new{T}()
        return syssolver
    end
end

function load(
    syssolver::SymIndefDenseSystemSolver{T},
    solver::Solver{T},
    ) where {T <: Real}
    model = solver.model
    (n, p, q) = (model.n, model.p, model.q)
    npq = n + p + q

    # fill symmetric lower triangle
    lhs_sub = zeros(T, npq, npq)
    @views copyto!(lhs_sub[n .+ (1:p), 1:n], model.A)
    @views copyto!(lhs_sub[n + p .+ (1:q), 1:n], model.G)
    @assert lhs_sub isa Matrix{T}
    syssolver.lhs_sub = Symmetric(lhs_sub, :L)
    syssolver.lhs_sub_fact = zero(syssolver.lhs_sub)

    setup_point_sub(syssolver, model)

    return syssolver
end

function update_lhs(syssolver::SymIndefDenseSystemSolver, solver::Solver)
    model = solver.model
    z_start = model.n + model.p
    lhs_sub = syssolver.lhs_sub.data

    for (cone_k, idxs_k) in zip(model.cones, model.cone_idxs)
        z_rows_k = z_start .+ idxs_k
        H_k = (Cones.use_dual_barrier(cone_k) ? Cones.hess :
            Cones.inv_hess)(cone_k)
        @. lhs_sub[z_rows_k, z_rows_k] = -H_k
    end

    solver.time_upfact += @elapsed syssolver.fact =
        symm_fact_copy!(syssolver.lhs_sub_fact, syssolver.lhs_sub)

    if !issuccess(syssolver.fact)
        println("symmetric linear system factorization failed")
    end

    solve_subsystem3(syssolver, solver, syssolver.sol_const, syssolver.rhs_const)

    return syssolver
end

function solve_subsystem3(
    syssolver::SymIndefDenseSystemSolver,
    ::Solver,
    sol::Point,
    rhs::Point,
    )
    ldiv!(sol.vec, syssolver.fact, rhs.vec)
    return sol
end

#=
indirect (using LinearMaps and IterativeSolvers)
TODO
- precondition
- optimize operations
- tune tolerances etc
- try to make initial point in sol a good guess (currently zeros)
=#

mutable struct SymIndefIndirectSystemSolver{T <: Real} <: SymIndefSystemSolver{T}
    lhs::LinearMaps.LinearMap{T}
    rhs_sub::Point{T}
    sol_sub::Point{T}
    sol_const::Point{T}
    rhs_const::Point{T}
    SymIndefIndirectSystemSolver{T}() where {T <: Real} = new{T}()
end

function load(
    syssolver::SymIndefIndirectSystemSolver{T},
    solver::Solver{T},
    ) where {T <: Real}
    model = solver.model
    (n, p, q) = (model.n, model.p, model.q)
    x_idxs = 1:n
    y_idxs = n .+ (1:p)
    z_start = n + p
    z_idxs = z_start .+ (1:q)

    function symindef_mul(b::AbstractVector, a::AbstractVector)
        # x part
        @views mul!(b[x_idxs], model.A', a[y_idxs])
        @views mul!(b[x_idxs], model.G', a[z_idxs], true, true)
        # y part
        @views mul!(b[y_idxs], model.A, a[x_idxs])
        # z part
        for (cone_k, idxs_k) in zip(model.cones, model.cone_idxs)
            z_rows_k = z_start .+ idxs_k
            prod_fun = (Cones.use_dual_barrier(cone_k) ? Cones.hess_prod! :
                Cones.inv_hess_prod!)
            @views prod_fun(b[z_rows_k], a[z_rows_k], cone_k)
        end
        @views mul!(b[z_idxs], model.G, a[x_idxs], true, -one(T))
        return b
    end

    syssolver.lhs = LinearMaps.LinearMap{T}(symindef_mul, z_start + q,
        ismutating = true, issymmetric = true, isposdef = false)

    setup_point_sub(syssolver, model)

    return syssolver
end

function update_lhs(syssolver::SymIndefIndirectSystemSolver, solver::Solver)
    solve_subsystem3(syssolver, solver, syssolver.sol_const, syssolver.rhs_const)
    return syssolver
end

function solve_subsystem3(
    syssolver::SymIndefIndirectSystemSolver,
    ::Solver,
    sol::Point,
    rhs::Point,
    )
    sol.vec .= 0 # initially_zero = true
    # TODO tune options, initial guess?
    IterativeSolvers.minres!(sol.vec, syssolver.lhs, rhs.vec,
        initially_zero = true) #, maxiter = 2 * size(sol.vec, 1))
    return sol
end


text:
using Test
using Cthulhu
using StaticArrays
using VariableTemplates

abstract type OneLayerModel end

struct EmptyModel <: OneLayerModel end
vars_state(m::EmptyModel, T) = @vars()

struct ScalarModel <: OneLayerModel end
vars_state(m::ScalarModel, T) = @vars(x::T)

struct VectorModel{N} <: OneLayerModel end
vars_state(m::VectorModel{N}, T) where {N} = @vars(x::SVector{N, T})

struct MatrixModel{N,M} <: OneLayerModel end
vars_state(m::MatrixModel{N,M}, T) where {N,M} = @vars(x::SHermitianCompact{N, T, M})

abstract type TwoLayerModel end

Base.@kwdef struct CompositModel{Nv,N,M} <: TwoLayerModel
    empty_model = EmptyModel()
    scalar_model = ScalarModel()
    vector_model = VectorModel{Nv}()
    matrix_model = MatrixModel{N,M}()
end
function vars_state(m::CompositModel, T)
    @vars begin
        empty_model::vars_state(m.empty_model, T)
        scalar_model::vars_state(m.scalar_model, T)
        vector_model::vars_state(m.vector_model, T)
        matrix_model::vars_state(m.matrix_model, T)
    end
end

Base.@kwdef struct NTupleModel <: OneLayerModel
    scalar_model = ScalarModel()
end

function vars_state(m::NTupleModel, T)
    @vars begin
        scalar_model::vars_state(m.scalar_model, T)
    end
end

vars_state(m::NTuple{N, NTupleModel}, FT) where {N} =
    Tuple{ntuple(i -> vars_state(m[i], FT), N)...}

Base.@kwdef struct NTupleContainingModel{N,Nv} <: TwoLayerModel
    ntuple_model = ntuple(i -> NTupleModel(), N)
    vector_model = VectorModel{Nv}()
    scalar_model = ScalarModel()
end

function vars_state(m::NTupleContainingModel, T)
    @vars begin
        ntuple_model::vars_state(m.ntuple_model, T)
        vector_model::vars_state(m.vector_model, T)
        scalar_model::vars_state(m.scalar_model, T)
    end
end

FT = Float32;

Nv = 3
N = 3
m = NTupleContainingModel{N,Nv}()
st = vars_state(m, FT)

@descend concretize(st)


text:
#
# semaphores.jl --
#
# Implements semaphores for Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of InterProcessCommunication.jl released under the MIT
# "expat" license.
#
# Copyright (C) 2016-2019, Éric Thiébaut
# (https://github.com/emmt/InterProcessCommunication.jl).
#

"""
## Named Semaphores

```julia
Semaphore(name, value; perms=0o600, volatile=true) -> sem
```

creates a new named semaphore identified by the string `name` of the form
`"/somename"` and initial value set to `value`.  An instance of
`Semaphore{String}` is returned.  Keyword `perms` can be used to specify access
permissions.  Keyword `volatile` specify whether the semaphore should be
unlinked when the returned object is finalized.

```julia
Semaphore(name) -> sem
```

opens an existing named semaphore and returns an instance of
`Semaphore{String}`.

To unlink (remove) a persistent named semaphore, simply do:

```julia
rm(Semaphore, name)
```

If the semaphore does not exists, the error is ignored.  A `SystemError` is
however thrown for other errors.

For maximum flexibility, an instance of a named semaphore may also be created
by:

```julia
open(Semaphore, name, flags, mode, value, volatile) -> sem
```

where `flags` may have the bits `IPC.O_CREAT` and `IPC.O_EXCL` set, `mode`
specifies the granted access permissions, `value` is the initial semaphore
value and `volatile` is a boolean indicating whether the semaphore should be
unlinked when the returned object `sem` is finalized.  The values of `mode` and
`value` are ignored if an existing named semaphore is open.

## Anonymous Semaphores

Anonymous semaphores are backed by *memory* objects providing the necessary
storage.

```julia
Semaphore(mem, value; offset=0, volatile=true) -> sem
```

initializes an anonymous semaphore backed by memory object `mem` with initial
value set to `value` and returns an instance of `Semaphore{typeof(mem)}`.
Keyword `offset` can be used to specify the address (in bytes) of the semaphore
data relative to `pointer(mem)`.  Keyword `volatile` specify whether the
semaphore should be destroyed when the returned object is finalized.

```julia
Semaphore(mem; offset=0) -> sem
```

yields an an instance of `Semaphore{typeof(mem)}` associated with an
initialized anonymous semaphore and backed by memory object `mem` at relative
position (in bytes) specified by keyword `offset`.

The number of bytes needed to store an anonymous semaphore is given by
`sizeof(Semaphore)` and anonymous semaphore must be aligned in memory at
multiples of the word size (that is `Sys.WORD_SIZE >> 3` in bytes).  Memory
objects used to store an anonymous semaphore must implement two methods:
`pointer(mem)` and `sizeof(mem)` to yield respectively the base address and the
size (in bytes) of the associated memory.

See also: [`post`](@ref), [`wait`](@ref), [`timedwait`](@ref),
          [`trywait`](@ref).

"""
function Semaphore(name::AbstractString, value::Integer;
                   perms::Integer = S_IRUSR | S_IWUSR,
                   volatile::Bool = true)
    val = _check_semaphore_value(value)
    mode = maskmode(perms)
    flags = O_CREAT | O_EXCL
    return open(Semaphore, name, flags, mode, val, volatile)
end

function Semaphore(name::AbstractString)
    flags = zero(Cint)
    mode = zero(_typeof_mode_t)
    value = zero(Cuint)
    return open(Semaphore, name, flags, mode, value, false)
end

function Base.open(::Type{Semaphore}, name::AbstractString, flags::Integer,
                   mode::Unsigned, value::Unsigned, volatile::Bool)
    ptr = _sem_open(name, flags, mode, value)
    systemerror("sem_open", ptr == SEM_FAILED)
    sem = Semaphore{String}(ptr, name)
    if volatile
        finalizer(_close_and_unlink, sem,)
    else
        finalizer(_close, sem)
    end
    return sem
end

# Initialize and anonymous semaphore.
function Semaphore(mem::M, value::Integer;
                   offset::Integer = 0,
                   volatile::Bool = true) where {M}
    val = _check_semaphore_value(value)
    ptr = _get_semaphore_address(mem, offset)
    systemerror("sem_init",
                _sem_init(ptr, true, val) != SUCCESS)
    sem = Semaphore{M}(ptr, mem)
    if volatile
        finalizer(_destroy, sem)
    end
    return sem
end

# Connect to an initialized anonymous semaphore.
function Semaphore(mem::M; offset::Integer = 0) where {M}
    ptr = _get_semaphore_address(mem, offset)
    return Semaphore{M}(ptr, mem)
end

# Unlink (remove) a named semaphore.
function Base.rm(::Type{Semaphore}, name::AbstractString)
    if _sem_unlink(name) != SUCCESS
        errno = Libc.errno()
        # On OSX, errno is EINVAL when semaphore does not exists.
        if errno != (@static Sys.isapple() ? Libc.EINVAL : Libc.ENOENT)
            throw_system_error("sem_unlink", errno)
        end
    end
end

Base.sizeof(::Type{<:Semaphore}) = _sizeof_sem_t
Base.typemin(::Type{<:Semaphore}) = zero(Cuint)
Base.typemax(::Type{<:Semaphore}) = SEM_VALUE_MAX

function _check_semaphore_value(value::Integer)
    typemin(Semaphore) ≤ value ≤ typemax(Semaphore) ||
        throw_argument_error("invalid semaphore value (", value, ")")
    return convert(Cuint, value)
end

function _get_semaphore_address(mem, off::Integer)::Ptr{Cvoid}
    off ≥ 0 || throw_argument_error("offset must be nonnegative (", off, ")")
    ptr, siz = get_memory_parameters(mem)
    siz ≥ off + _sizeof_sem_t ||
        throw_argument_error("not enough memory at given offset")
    ptr += off
    align = (Sys.WORD_SIZE >> 3) # assume alignment is word size (in bytes)
    rem(convert(Int, ptr), align) == 0 ||
        throw_argument_error("address of semaphore must be a multiple of ",
                             align, " bytes")
    return ptr
end

# Finalize a named semaphore.
function _close(obj::Semaphore{String})
    _sem_close(obj.ptr)
end

# Finalize a named semaphore.
function _close_and_unlink(obj::Semaphore{String})
    _sem_close(obj.ptr)
    _sem_unlink(obj.lnk)
end

# Finalize an anonymous semaphore.
function _destroy(obj::Semaphore)
    _sem_destroy(obj.ptr)
end

function Base.getindex(sem::Semaphore)
    val = Ref{Cint}()
    systemerror("sem_getvalue", _sem_getvalue(sem.ptr, val) != SUCCESS)
    return val[]
end

"""
```julia
post(sem)
```

increments (unlocks) the semaphore `sem`.  If the semaphore's value
consequently becomes greater than zero, then another process or thread blocked
in a [`wait`](@ref) call on this semaphore will be woken up.

See also: [`Semaphore`](@ref), [`wait`](@ref), [`timedwait`](@ref),
          [`trywait`](@ref).

"""
post(sem::Semaphore) =
    systemerror("sem_post", _sem_post(sem.ptr) != SUCCESS)

"""
```julia
wait(sem)
```

decrements (locks) the semaphore `sem`.  If the semaphore's value is greater
than zero, then the decrement proceeds and the function returns immediately.
If the semaphore currently has the value zero, then the call blocks until
either it becomes possible to perform the decrement (i.e., the semaphore value
rises above zero), or a signal handler interrupts the call (in which case an
instance of `InterruptException` is thrown).  A `SystemError` may be thrown if
an unexpected error occurs.

See also: [`Semaphore`](@ref), [`post`](@ref), [`timedwait`](@ref),
          [`trywait`](@ref).

"""
function Base.wait(sem::Semaphore)
    if _sem_wait(sem.ptr) != SUCCESS
        code = Libc.errno()
        if code == Libc.EINTR
            throw(InterruptException())
        else
            throw_system_error("sem_wait", code)
        end
    end
    nothing
end

"""
```julia
timedwait(sem, secs)
```

decrements (locks) the semaphore `sem`.  If the semaphore's value is greater
than zero, then the decrement proceeds and the function returns immediately.
If the semaphore currently has the value zero, then the call blocks until
either it becomes possible to perform the decrement (i.e., the semaphore value
rises above zero), or the limit of `secs` seconds expires (in which case an
instance of [`TimeoutError`](@ref) is thrown), or a signal handler interrupts
the call (in which case an instance of `InterruptException` is thrown).

See also: [`Semaphore`](@ref), [`post`](@ref), [`wait`](@ref),
          [`trywait`](@ref).

"""
Base.timedwait(sem::Semaphore, secs::Real) =
    timedwait(sem::Semaphore, convert(Float64, secs))

function Base.timedwait(sem::Semaphore, secs::Float64)
    tsref = Ref{TimeSpec}(time() + secs)
    if _sem_timedwait(sem.ptr, tsref) != SUCCESS
        code = Libc.errno()
        if code == Libc.EINTR
            throw(InterruptException())
        elseif code == Libc.ETIMEDOUT
            throw(TimeoutError())
        else
            throw_system_error("sem_timedwait", code)
        end
    end
    nothing
end

"""
```julia
trywait(sem) -> boolean
```

attempts to immediately decrement (lock) the semaphore `sem` returning `true`
if successful.  If the decrement cannot be immediately performed, then the call
returns `false`.  If an interruption is received or if an unexpected error
occurs, an exception is thrown (`InterruptException` or `SystemError`
repectively).

See also: [`Semaphore`](@ref), [`post`](@ref), [`wait`](@ref),
          [`timedwait`](@ref).

"""
function trywait(sem::Semaphore)
    if _sem_trywait(sem.ptr) == SUCCESS
        return true
    end
    code = Libc.errno()
    if code == Libc.EAGAIN
        return false
    end
    if code == Libc.EINTR
        throw(InterruptException())
    else
        throw_system_error("sem_trywait", code)
    end
end


text:
#=
P04 (*) Find the number of elements of a list.
=#

function my_length(coll)
    # With built-in `length` function
    length(coll)
end

function len_with_count(coll, count)
    # With tail recursion
    if isempty(coll)
        count
    else
        len_with_count(coll[2:end], count + 1)
    end
end

function my_length_2(coll)
    # Counting from 0
    len_with_count(coll, 0)
end

@assert my_length(()) == 0
@assert my_length((:a, :b, :c)) == 3
@assert my_length((1, 2, 3, 4)) == 4
@assert my_length(((:a, :b), (:c, :d), (:e, :f))) == 3

@assert my_length_2(()) == 0
@assert my_length_2((:a, :b, :c)) == 3
@assert my_length_2((1, 2, 3, 4)) == 4
@assert my_length_2(((:a, :b), (:c, :d), (:e, :f))) == 3

println("Tests passed: JL-04.jl")


text:
export  println_time_flush,
    strderr!, tstats!, pvals!,
    strderr, tstats, pvals,
    Fstat!,
    _model,
    coeftable,
    evaluate_likelihood

function println_time_flush(str)
    println(Dates.format(now(), "HH:MM:SS   ") * str)
    flush(stdout)
end

coef(  leo::LocalEstObj) = theta1(leo)
stderr(leo::LocalEstObj) = sqrt.(diag(invhess(leo)))
tstats(leo::LocalEstObj) = coef(leo) ./ stderr(leo)

function pvals(leo::LocalEstObj)
    t = tstats(leo)
    p = 2 .* ccdf.(Normal(), abs.(t))
    return p
end

function coef_and_se(leo::LocalEstObj)
    se = stderr(leo)
    t = coef(leo) ./ se
    p = pvals(leo)
    return hcat(coef(leo), se, t, p)
end

tstats!(leo)      = (update_invhess!(leo); stderr(leo))
stderr!(leo)      = (update_invhess!(leo); tstats(leo))
pvals!(leo)       = (update_invhess!(leo); pvals(leo))
coef_and_se!(leo) = (update_invhess!(leo); coef_and_se(leo))

@deprecate coef_and_se(args...) coeftable(args...)
@deprecate coef_and_se!(args...) coeftable(args...)

function Fstat!(leo::LocalEstObj; H0 = theta0(leo), alpha=0.05)
    k = _nparm(leo)
    err = coef(leo) .- H0
    waldtest = err'*invhess!(leo)*err
    p = ccdf(Chisq(k), waldtest)
    reject = p < alpha
    return waldtest, p, reject
end

function critval(alpha::Real=0.05, twosided::Bool=true)
    0 < alpha < 1 || throw(DomainError(alpha))
    alpha >= 0.5 && @warn "α = $alpha > 0.05"
    a = twosided ? alpha/2 : alpha
    return quantile(Normal(), a)
end

function coeftable(leo, alpha::Real=0.05)
    # see defn at
    # https://github.com/JuliaStats/StatsBase.jl/blob/da42557642046116097ddaca39fd5dc2c41402cc/src/statmodels.jl#L376-L401
    update_invhess!(leo)
    cc = coef(leo)
    se = stderr(leo)
    tt = tstats(leo)
    p = pvals(leo)
    ci = se .* critval(alpha, true)

    omlev100 = (1 - alpha)*100
    levstr = isinteger(omlev100) ? string(Integer(omlev100)) : string(omlev100)

    tbl = hcat(cc, se, tt, p, cc+ci, cc-ci)
    cols = ["Estimate","Std. Error","t value","Pr(>|t|)","Lower $levstr%","Upper $levstr%"]
    rows = coefnames(data(leo))
    pvalcol = 4
    return CoefTable(tbl, cols, rows, pvalcol)
end

function solve_model(d::DataSetofSets, theta, M, maxtime, OptimMethod=bfgs)
    check_finite(theta)
    leo = LocalEstObj(d, theta)
    reo = RemoteEstObj(leo, M)
    ew = EstimationWrapper(leo, reo)
    leograd = grad(leo)
    @eval @everywhere set_g_RemoteEstObj($reo)

    resetcount!()
    startcount!([100, 500, 100000,], [1, 5, 100,])
    opts = Optim.Options(show_trace=true, time_limit=maxtime, allow_f_increases=true)

    res = solve_model(ew, theta; OptimOpts=opts, OptimMethod=OptimMethod(ew))

    println(res)
    println("Recomputing final gradient / hessian")
    let dograd=true, theta=Optim.minimizer(res)
        parallel_simloglik!(ew, theta, dograd)
        update!(ew, theta, dograd)
    end
    @assert Optim.minimizer(res) == theta1(leo)
    println(coeftable(leo))
    print("Parameter estimates are\n\t")
    print(round.(theta1(leo); digits=3))
    print("\n\t")
    print(sprintf_binary(theta1(leo)))
    print("\n\n")
    return res, ew
end

function evaluate_likelihood(d::DataSetofSets, theta, M; compute=true)
    check_finite(theta)
    leo = LocalEstObj(d, theta)
    reo = RemoteEstObj(leo, M)
    ew = EstimationWrapper(leo, reo)
    leograd = grad(leo)

    hess(leo) .= Matrix(I, length(leograd), length(leograd))
    theta1(leo) .= theta

    if compute
        @eval @everywhere set_g_RemoteEstObj($reo)
        dograd=true
        parallel_simloglik!(ew, theta, dograd)
        update!(ew, theta, dograd)
    end

    println(coeftable(leo))
    print("Parameter estimates are\n\t")
    print(round.(theta1(leo); digits=3))
    print("\n\t")
    print(sprintf_binary(theta1(leo)))
    print("\n\n")
    return ew
end


text:
#=
fft:
- Julia version:
- Author: ymocquar
- Date: 2019-11-15
=#
# function that reverse the order of the pos lowest bits
using FFTW
function _reverse_num(num, pos)
    result = 0
    pos_m1 = pos - 1
    for i = 0:pos_m1
        if (num & (1 << i)) != 0
            result |= 1 << (pos_m1 - i)
        end
    end
    return result
end
"""
    PrepareFftBig( size_fft::Unsigned, [T=BigFloat])

Immutable structure to operate fft transform,
x is the type of non transformed data also called signal.

# Arguments :
- `size_fft::Integer` : Number of values, must be a power of two
- `[T=BigFloat | x::T ]` : type of the values

# Implementation
- size_fft : size of the signal
- tab_permut : permutation
- root_one : size order roots of one
- root_one_conj : conjugate of root_one

"""
struct PrepareFftBig
    size_fft::Any
    tab_permut::Any
    root_one::Any
    root_one_conj::Any
    type::DataType
    function PrepareFftBig(size_fft::Integer, x::T) where {T<:AbstractFloat}
        @assert prevpow(2, size_fft) == size_fft "size_fft=$size_fft is not a power of 2"
        power = convert(Int64, log2(size_fft))
        tab_permut = zeros(Int64, size_fft)
        root_one = zeros(Complex{T}, size_fft)
        root_one_conj = zeros(Complex{T}, size_fft)
        prec = precision(T)
        setprecision(prec + 32) do
            for i = 1:size_fft
                tab_permut[i] = _reverse_num(i - 1, power) + 1
                root_one[i] = round(
                    exp(one(T) * 2big(pi) * im * (i - 1) / size_fft),
                    digits = prec + 16,
                    base = 2,
                )
                root_one_conj[i] = round(
                    exp(-one(T) * 2big(pi) * im * (i - 1) / size_fft),
                    digits = prec + 16,
                    base = 2,
                )
            end
        end
        return new(
            size_fft,
            tab_permut,
            Complex{T}.(root_one),
            Complex{T}.(root_one_conj),
            typeof(x),
        )
    end
end
PrepareFftBig(size_fft::Integer, type::DataType) = PrepareFftBig(size_fft, one(type))
PrepareFftBig(size_fft::Integer) = PrepareFftBig(size_fft, one(BigFloat))
function fftbig!(par::PrepareFftBig, signal; flag_inv = false)
    s = size(signal, 2)
    @assert prevpow(2, s) == s "size_fft(signal)=$s is not a power of 2"
    s_div2 = div(s, 2)
    len = s
    n_len = len >> 1
    nb_r = 1
    rootO = flag_inv ? par.root_one : par.root_one_conj
    #    prec= precision(real(rootO[1]))
    #    setprecision(prec+32) do
    while n_len != 0
        start = 1
        suite = start + n_len
        for i = 1:nb_r
            deb = 1
            for j = deb:nb_r:s_div2
                signal[:, start], signal[:, suite] = (signal[:, start] + signal[:, suite]),
                (signal[:, start] - signal[:, suite]) * rootO[j]
                start += 1
                suite += 1
            end
            start = suite
            suite = start + n_len
        end
        len = n_len
        n_len >>= 1
        nb_r <<= 1
    end
    #    end
    signal .= signal[:, par.tab_permut]

    if flag_inv
        signal ./= s
    end
    return signal
end
function fftbig(par::PrepareFftBig, signal; flag_inv = false)
    fl = flag_inv
    return fftbig!(
        par::PrepareFftBig,
        copy(convert(Array{Complex{par.type}}, signal)),
        flag_inv = fl,
    )
end
fftgen(_::Any, t::Array{Complex{Float64}}) = fft(t, (2,))
fftgen(_::Any, t::Array{Float64}) = fft(t, (2,))
fftgen(p::PrepareFftBig, t::Array{T}) where {T<:AbstractFloat} = fftbig(p, t)
fftgen(p::PrepareFftBig, t::Array{Complex{T}}) where {T<:AbstractFloat} = fftbig(p, t)
ifftgen(_::Any, t::Array{Complex{Float64}}) = ifft(t, (2,))
ifftgen(_::Any, t::Array{Float64}) = ifft(t, (2,))
ifftgen(p::PrepareFftBig, t::Array{T}) where {T<:AbstractFloat} =
    fftbig(p, t, flag_inv = true)
ifftgen(p::PrepareFftBig, t::Array{Complex{T}}) where {T<:AbstractFloat} =
    fftbig(p, t, flag_inv = true)