text:
@testset "rl_unplugged_atari" begin
    ds = rl_unplugged_atari_dataset(
        "Pong",
        1,
        [1, 2];
        shuffle_buffer_size = 10_000,
        tf_reader_bufsize = 1*1024*1024,
        tf_reader_sz = 10_000,
        batch_size = 256,
        n_preallocations = Threads.nthreads() * 12
    )

    @test typeof(ds) <: AbstractChannel

    data_1 = take!(ds)

    frame_size = 84
    stack_size = 4
    batch_size = 256

    @test size(data_1.state) == (frame_size, frame_size, stack_size, batch_size)
    @test size(data_1.next_state) == (frame_size, frame_size, stack_size, batch_size)
    @test size(data_1.action) == (batch_size,)
    @test size(data_1.next_action) == (batch_size,)
    @test size(data_1.reward) == (batch_size,)
    @test size(data_1.terminal) == (batch_size,)
    @test size(data_1.episode_id) == (batch_size,)
    @test size(data_1.episode_return) == (batch_size,)

    @test typeof(data_1.state) == Array{UInt8, 4}
    @test typeof(data_1.next_state) == Array{UInt8, 4}
    @test typeof(data_1.action) == Vector{Int64}
    @test typeof(data_1.next_action) == Vector{Int64}
    @test typeof(data_1.reward) == Vector{Float32}
    @test typeof(data_1.terminal) == Vector{Bool}
    @test typeof(data_1.episode_id) == Vector{Int64}
    @test typeof(data_1.episode_return) == Vector{Float32}

    take!(ds)
    take!(ds)
    data_2 = take!(ds)

    @test data_1.state != data_2.state
    @test data_1.next_state != data_2.next_state
    @test data_1.action != data_2.action
    @test data_1.next_action != data_2.next_action
    @test data_1.reward != data_2.reward
    @test data_1.episode_id != data_2.episode_id
    @test data_1.episode_return != data_2.episode_return

end


text:
using JuMP, NLopt, Test

model = Model(NLopt.Optimizer)
set_optimizer_attribute(model, "algorithm", :LD_MMA)

a1 = 2
b1 = 0
a2 = -1
b2 = 1

@variable(model, x1)
@variable(model, x2 >= 0)

@NLobjective(model, Min, sqrt(x2))
@NLconstraint(model, x2 >= (a1*x1+b1)^3)
@NLconstraint(model, x2 >= (a2*x1+b2)^3)

set_start_value(x1, 1.234)
set_start_value(x2, 5.678)

optimize!(model)

println("got ", objective_value(model), " at ", [value(x1), value(x2)])

@test_approx_eq_eps value(x1) 1/3 1e-5
@test_approx_eq_eps value(x2) 8/27 1e-5
@test_approx_eq_eps objective_value(m) sqrt(8/27) 1e-5
@test termination_status(model) == MOI.OPTIMAL

ext: jl


text:
# https://www.hackerrank.com/challenges/s10-basic-statistics/problem

using Statistics
using Printf

data = readlines()
n_data = parse(Float64, data[1])
values = [parse(Float64, x) for x in split(data[2])]

function mode(x)
    a_keys = sort(unique(x))
    a_counts = zeros(length(a_keys))

    for val in x
        ma = a_keys.==val
        a_counts[ma] .+= 1
    end

    max_counts = maximum(a_counts)
    if max_counts == 1.0
        return a_keys[1]
    else
        ma = a_counts.==max_counts
        return minimum(a_keys[ma])
    end
end

mean_val = mean(values)
median_val = median(values)
mode_val = mode(values)

@printf("%.1f\n", mean_val)
@printf("%.1f\n", median_val)
@printf("%.1f\n", mode_val)

max_stars_repo_path: Practice/Tutorials/10 Days of Statistics/Day 00 Mean, Median, and Mode.jl

max_stars_repo_name: mtaufanr/hackerrank-code-julia

max_issues_repo_name: mtaufanr/hackerrank-code-julia

max_forks_repo_name: mtaufanr/hackerrank-code-julia


text:
using Revise
using AdFem

mmesh3 = Mesh3(10,10,10,0.1)
xy = gauss_nodes(mmesh3)

x, y, z = xy[:,1], xy[:,2], xy[:,3]
f = @. (1-x) * x * (1-y) * y * (1+x^2+z^2)

F = compute_fem_source_term1(f, mmesh3)
M = compute_fem_mass_matrix1((@. 1+x^2+z^2), mmesh3)

bd = bcnode(mmesh3)
M, F = impose_Dirichlet_boundary_conditions(M, F, bd, zeros(length(bd)))
s = M\F

sess = Session(); init(sess)
SOL = run(sess, s)

f = vtk(mmesh3, "test")
f["u"] = SOL
save(f)

xy = fem_nodes(mmesh3)
x, y, z = xy[:,1], xy[:,2], xy[:,3]

f = vtk(mmesh3, "test2")
f["u"] =  @. (1-x) * x * (1-y) * y
save(f)


text:
#
# Quadratic fit
#

struct Quadratic
  a :: Float64
  b :: Float64
  c :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

function fitquadratic(X :: AbstractVector, Y :: AbstractVector, options :: Options)
  X, Y = checkdata(X,Y)
  @. model(x,p) = p[1]*x^2 + p[2]*x + p[3]
  fit = find_best_fit(model, X, Y, 3, options)
  R = pearson(X,Y,model,fit)
  x, y, ypred = finexy(X,options.fine,model,fit)
  return Quadratic(fit.param...,R,x,y,ypred,fit.resid)
end
fitquadratic(X :: AbstractVector, Y :: AbstractVector) = fitquadratic(X,Y,Options())
fitquad(X :: AbstractVector, Y :: AbstractVector) = fitquadratic(X,Y,Options())
fitquad(X :: AbstractVector, Y :: AbstractVector, options :: Options) = fitquadratic(X,Y,options)

function Base.show( io :: IO, fit :: Quadratic )
  println("")
  println(" ------------------- Quadratic Fit ------------- ")
  println("")
  println(" Equation: y = ax^2 + bx + c ")
  println("")
  println(" With: a = ", fit.a)
  println("       b = ", fit.b)
  println("       c = ", fit.c)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

export fitquad, fitquadratic


text:
using antras_chor_2018
using Test

@testset "antras_chor_2018.jl" begin
    # Write your tests here.
end


text:
####################
# PullRequest Type #
####################

type PullRequest <: GitHubType
    base::Nullable{Branch}
    head::Nullable{Branch}
    number::Nullable{Int}
    id::Nullable{Int}
    comments::Nullable{Int}
    commits::Nullable{Int}
    additions::Nullable{Int}
    deletions::Nullable{Int}
    changed_files::Nullable{Int}
    state::Nullable{GitHubString}
    title::Nullable{GitHubString}
    body::Nullable{GitHubString}
    merge_commit_sha::Nullable{GitHubString}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    closed_at::Nullable{Dates.DateTime}
    merged_at::Nullable{Dates.DateTime}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    assignee::Nullable{Owner}
    user::Nullable{Owner}
    merged_by::Nullable{Owner}
    milestone::Nullable{Dict}
    _links::Nullable{Dict}
    mergeable::Nullable{Bool}
    merged::Nullable{Bool}
    locked::Nullable{Bool}
end

PullRequest(data::Dict) = json2github(PullRequest, data)
PullRequest(number::Real) = PullRequest(Dict("number" => number))

namefield(pr::PullRequest) = pr.number

###############
# API Methods #
###############

function pull_requests(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/pulls"; options...)
    return map(PullRequest, results), page_data
end

function pull_request(repo, pr; options...)
    result = gh_get_json("/repos/$(name(repo))/pulls/$(name(pr))"; options...)
    return PullRequest(result)
end


text:
type SGsqrtsd <: StepData
	getStep::Function
	stepString::AbstractString
	function SGsqrtsd(stepString::AbstractString)
		new(computeSGsqrtStep,stepString)
	end
	function SGsqrtsd()
		new(computeSGsqrtStep,"SGsqrt")
	end
end

function computeSGsqrtStep(x, k, gnum, sd::SGsqrtsd, problem::Problem; outputLevel =0)

	func = consume(problem.getNextSampleFunction)

	(f, g) = func(x)

	gnum += 1

	(g/sqrt(k+1), gnum)

end


text:
m=mocc
m=mdolocc
#m=mpen

custom_vars=[:dolocc_reduction,:occ_reduction,:pen_reduction,:data_b_e_nb,:data_nb_ne_b,:whyo,:iso]
required_vars=vcat([m.y_var,:panid,m.logvar],cfg[:random_demos],cfg[:random_campaigns],cfg[:scoring_vars])
vars=setdiff(m.vars,vcat(required_vars,custom_vars))

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

    function rmVars(v::Array{Symbol})
        v=setdiff(v,[:group])
        return setdiff(vars,v)
    end

    function getDSname()
        mn = Dict("occ"=>"finaldata_occ","dolocc"=>"finaldata_dolocc","pen"=>"finaldata_pen")
        return get(mn, m.modelName,NA)
    end

    function rv()
        return convert(Array{Symbol},sort(   map(x->lowercase(x),rcopy("  setdiff(names("*getDSname()*"),c(all_mandatory_vars,scoring_vars,random_campaigns,random_demos,'panid')) "))  ))
    end

    function rchk()
        onlyinJ=setdiff(vars,rv())
        onlyinR=setdiff(rv(),vars)
        lj=length(onlyinJ)
        lr=length(onlyinR)
        lj>lr ? DataFrame(onlyinJ=onlyinJ,onlyinR=vcat(onlyinR,fill(:x, lj-lr))) : DataFrame(DataFrame(onlyinJ=vcat(onlyinJ,fill(:x, lr-lj)),onlyinR=onlyinR))
    end

    function isSame()
        println(getDSname()*" : J vals : ",length(vars),"  --- R vals : ",length(rv()))
        sort(vars)==sort(rv())
    end

    function rlst(v::String)
        sort(map(x->Symbol(lowercase(x)),rcopy(v)))
    end

include("/home/rmadmin/g/StatStack/src/RCode.jl")
#dfd=df_data[(df_data[:iso].==false)&(df_data[:buyer_pos_p1].==1),:]

if   m.modelName == "occ"
#OCC
R"""
 y_var <- "Trps_POS_P1"
finaldata_occ <- initial_data[Buyer_Pos_P1=="1",]
# drop variables with only one level
indices <- sapply(finaldata_occ,function(x) length(unique(x))==1)
finaldata_occ <- finaldata_occ[,!indices,with=F]
"""
end
if  m.modelName == "dolocc"
#DOLOCC
R"""
y_var <- "Dol_per_Trip_POS_P1"
finaldata_dolocc <- initial_data[Buyer_Pos_P1=="1",]
# drop variables with only one level
indices <- sapply(finaldata_dolocc,function(x) length(unique(x))==1)
finaldata_dolocc <- finaldata_dolocc[,!indices,with=F]
"""
end

#After Single level reove
isSame()

# GT: Remove cause these will get rmoved in mc test..... so need to QC
if   m.modelName == "occ"
#OCC
R"""
finaldata_occ <- as.data.frame(finaldata_occ)
mc_cols <- c("Fea_or_Dis_Trps_Shr_DPP_P1","Fea_or_Dis_Trps_Shr_DPP_P2","Fea_or_Dis_Trps_Shr_DPP_P3","Fea_or_Dis_Trps_Shr_DPP_P4")
finaldata_occ <- finaldata_occ[, !names(finaldata_occ) %in% mc_cols ]
"""
end
if   m.modelName == "dolocc"
#DOLOCC
R"""
finaldata_dolocc <- as.data.frame(finaldata_dolocc)
mc_cols <- c("Fea_or_Dis_Trps_Shr_DPP_P1","Fea_or_Dis_Trps_Shr_DPP_P2","Fea_or_Dis_Trps_Shr_DPP_P3","Fea_or_Dis_Trps_Shr_DPP_P4")
finaldata_dolocc <- finaldata_dolocc[, !names(finaldata_dolocc) %in% mc_cols ]
"""
end

#Copy R Dataset to avoid sampling
dfd=rcopy(getDSname())
lowercase!(dfd)
mc_cols = [:fea_or_dis_trps_shr_dpp_p1,:fea_or_dis_trps_shr_dpp_p2,:fea_or_dis_trps_shr_dpp_p3,:fea_or_dis_trps_shr_dpp_p4]
vars=setdiff(vars,mc_cols)

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

if  m.modelName == "occ"
#OCC
R"""
finaldata_occ <- as.data.frame(finaldata_occ)
if(exists('dolocc_outliers_panid_to_remove', envir = environment())){
  finaldata_occ <- finaldata_occ[, !(names(finaldata_occ) %in% c("Dol_per_Trip_POS_P1","Nonbuyer_Pre_P1","Buyer_Pre_P1","Dol_per_Trip_PRE_P1"))]
  finaldata_occ <- finaldata_occ[!finaldata_occ$panid %in% dolocc_outliers_panid_to_remove,]
}else{
  finaldata_occ <- finaldata_occ[, !(names(finaldata_occ) %in% c("Dol_per_Trip_POS_P1","Nonbuyer_Pre_P1","Buyer_Pre_P1","Dol_per_Trip_PRE_P1"))]
}

# set initial model formula
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"Trps_PRE_P1","panid",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1))")
options(warn=-1)

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

pval_excluded_vars <- c(rownames(as.data.frame(which(summary(model_glm_total)$coefficient[,4]>0.7))),setdiff(names(model_glm_total$coefficients),rownames(summary(model_glm_total)$coefficients)))
finaldata_occ <- finaldata_occ[, !names(finaldata_occ) %in% c(pval_excluded_vars)]
"""
end
if  m.modelName == "dolocc"
#DOLOCC
R"""
#dolocc_mandatory_covariate <- c("Nonbuyer_Pre_P1")
#finaldata_dolocc <- finaldata_dolocc[, !names(finaldata_dolocc) %in% c(dolocc_mandatory_covariate)]
dolocc_mandatory_covariate <- c()

# set initial model formula
finaldata_dolocc <- as.data.frame(finaldata_dolocc)
finaldata_dolocc <- finaldata_dolocc[, !(names(finaldata_dolocc) %in% c("Trps_POS_P1","Buyer_Pre_P1","Trps_PRE_P1"))]
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_dolocc[,!names(finaldata_dolocc)%in%c(y_var,"panid","Dol_per_Trip_PRE_P1",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," ")," + offset(log(Dol_per_Trip_PRE_P1+1))")
options(warn=-1)

# run glm model
model_glm_total <- glm(Formula_fixed,data=finaldata_dolocc,family=Gamma(link="log"))
summary(model_glm_total)

##### Variable selection - Condition 1 #####
# Remove pvalues greater than 0.7 and rerun the glm model
pval_excluded_vars <- setdiff(c(rownames(as.data.frame(which(summary(model_glm_total)$coefficient[,4]>0.7))),setdiff(names(model_glm_total$coefficients),rownames(summary(model_glm_total)$coefficients))),dolocc_mandatory_covariate)

finaldata_dolocc <- finaldata_dolocc[, !names(finaldata_dolocc) %in% c(pval_excluded_vars)]
"""
end
#filter(x->contains("nonbuyer_pre_p1",string(x)),vars)

#:fea_or_dis_trps_shr_dpp_p1
#:fea_or_dis_trps_shr_dpp_p4
#:fea_or_dis_trps_shr_dpp_p2
#:fea_or_dis_trps_shr_dpp_p3

#cmparr(pval_excluded_vars,j,[:fea_or_dis_trps_shr_dpp_p1,:fea_or_dis_trps_shr_dpp_p2,:fea_or_dis_trps_shr_dpp_p3,:fea_or_dis_trps_shr_dpp_p4])
    #PVals
    m.glm1_pvals = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )
    g1=m.glm1_pvals
    m.glm1_pvals_x=g1.sdf[g1.sdf[:pval].>0.7,:vars]
    vars = rmVars(vcat(m.glm1_pvals_x, g1.xvars ) )

isSame()

#  --------------------QC -------------------------------------
#sort(m.glm1_pvals_x) == rlst("pval_excluded_vars")
#cmparr(sort(m.glm1_pvals_x), rlst("pval_excluded_vars"))

#Cpn_UN_per_Trip_DPP_P4  -0.460885549 1.221531226  -0.37730149 7.059541e-01

# m.glm1_pvals.sdf[m.glm1_pvals.sdf[:vars].==:cpn_un_per_trip_dpp_p4,:]
#length(m.glm1_pvals.sdf[:vars])
#R"""length(summary(model_glm_total)$coefficients[,4])"""

#m.glm1_pvals.model
#R"""summary(model_glm_total)"""
#prd_1_qty_pre

# -0---------- QC end ---------------------------------------------

#writetable("CDW.csv", dfd, header = true)
if  m.modelName == "occ"
# OCC
R"""

# re-run glm model
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"Trps_PRE_P1","panid",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1))")
model_glm_condition1 <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
summary(model_glm_condition1)
rm(model_glm_total)

##### Variable selection - Condition 2 #####
# Remove variables that have z value lower than 2 and VIF greater than 15

# get variables with z value lower than 2
z_excluded_vars <- rownames(as.data.frame(which(abs(summary(model_glm_condition1)$coefficient[,3])<1.96)))

# get the variables with VIF greater than 15
multicol <- alias(model_glm_condition1)
multicols_toremove <- rownames(as.data.frame(multicol$Complete))

if (length(multicols_toremove)>0){
  finaldata_occ <- finaldata_occ[ , !names(finaldata_occ) %in% multicols_toremove]
  Formula_fixed <-  paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"panid",random_demos,random_campaigns)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1)) ")
  model_glm_new <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
} else {
  model_glm_new <- model_glm_condition1
}

vif_output <- (vif(model_glm_new))
vif_to_remove <- names(vif_output[vif_output > 15])

# final vars to exclude under condition 2
vars_condition2 <- intersect(z_excluded_vars,vif_to_remove)

# remove the variables under the above conditions
finaldata_occ <- finaldata_occ[ , !names(finaldata_occ) %in% c(vars_condition2,pval_excluded_vars)]
"""
end
if  m.modelName == "dolocc"
# DOLOCC
R"""
# re-run glm model
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_dolocc[,!names(finaldata_dolocc)%in%c(y_var,"panid","Dol_per_Trip_PRE_P1",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," ")," + offset(log(Dol_per_Trip_PRE_P1+1))")
model_glm_condition1 <- glm(Formula_fixed,data=finaldata_dolocc,family=Gamma(link="log"))
summary(model_glm_condition1)
rm(model_glm_total)

##### Variable selection - Condition 2 #####
# Remove variables that have z value lower than 2 and VIF greater than 15

# get variables with z value lower than 2
z_excluded_vars  <- rownames(as.data.frame(which(abs(summary(model_glm_condition1)$coefficient[,3])<1.96)))

# get the variables with VIF greater than 15
multicol <- alias(model_glm_condition1)
multicols_toremove <- rownames(as.data.frame(multicol$Complete))

if (length(multicols_toremove)>0){
  finaldata_dolocc <- finaldata_dolocc[ , !names(finaldata_dolocc) %in% multicols_toremove]
  Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_dolocc[,!names(finaldata_dolocc)%in%c(y_var,"panid")]),collapse=" + ")," ")," + offset(log(Dol_per_Trip_PRE_P1+1))")
  model_glm_new <- glm(Formula_fixed,data=finaldata_dolocc,family=Gamma(link="log"))
} else {
  model_glm_new <- model_glm_condition1
}

vif_output <- (vif(model_glm_new))
vif_to_remove <- names(vif_output[vif_output > 15])

# final vars to exclude under condition 2
vars_condition2 <- setdiff(intersect(z_excluded_vars,vif_to_remove),dolocc_mandatory_covariate)
finaldata_dolocc <- finaldata_dolocc[ , !names(finaldata_dolocc) %in% vars_condition2]
"""
end

    #Z & Vif
    m.glm2_ZnVIF = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )
    g2=m.glm2_ZnVIF
    vif!(g2)
    #z = g2.sdf[g2.sdf[:zval].<1.96,:vars]
    z = g2.sdf[abs(g2.sdf[:zval]).<1.96,:vars]
    #QC: sort(z) == sort(rlst("z_excluded_vars"))       cmparr(convert(Array{Symbol},z),rlst("z_excluded_vars"))
    v = g2.sdf[ !DataArrays.isna(g2.sdf[:vif])&(g2.sdf[:vif].>15),:vars]
    #QC: cmparr(convert(Array{Symbol},v),rlst("vif_to_remove"))
    m.glm2_ZnVIF_x =intersect(z,v)
    #QC: m.glm2_ZnVIF_x == rlst("vars_condition2")
    #QC: cmparr( m.glm2_ZnVIF_x ,rlst("vars_condition2"))
    vars = rmVars(vcat(m.glm2_ZnVIF_x, g2.xvars) )

#sort(m.glm2_ZnVIF_x)==sort(map(x->Symbol(lowercase(x)) ,rcopy("vars_condition2")))

isSame()

# --------QC ------------------
#cmparr(convert(Array{Symbol},m.glm2_ZnVIF_x), rlst("vars_condition2"))
#cmparr(convert(Array{Symbol},z), rlst("z_excluded_vars")
#g2.model
#R"""summary(model_glm_condition1)$coefficient"""
#R"""vif_output"""
#g2.sdf
#---------QC end --------------

if  m.modelName == "occ"
#OCC
R"""
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_occ[,!names(finaldata_occ)%in%c(y_var,"panid","Trps_PRE_P1",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," "),"+ offset(log(Trps_PRE_P1+1))")

# run model
model_glm <- glm(Formula_fixed,data=finaldata_occ,family=poisson(link="log"))
summary(model_glm)
rm(model_glm_condition1)

# p-values and signs checks
neutralvars <- setdiff(names(finaldata_occ),c(positivevars,negativevars))
pvalues <- summary(model_glm)$coefficients[,4][-1]
keep_covariate <- names(pvalues[pvalues<pvalue_lvl])
coeffs <- sign(summary(model_glm)$coefficients[,1])
varstokeep <- c();initialvars <- c()
for (i in negativevars) if (i %in% names(coeffs[coeffs<0])) varstokeep <- c(varstokeep,i)
neg <- varstokeep
for (j in positivevars) if (j %in% names(coeffs[coeffs>0])) varstokeep <- c(varstokeep,j)
pos <- setdiff(varstokeep,neg)
for (l in neutralvars) varstokeep <- c(varstokeep,l)  #append neut to varstokeep

for (k in varstokeep) if (k %in% keep_covariate) initialvars <- c(initialvars,k)

"""
end
if  m.modelName == "dolocc"
#DOLOCC
R"""
# re-run glm model
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(names(finaldata_dolocc[,!names(finaldata_dolocc)%in%c(y_var,"panid","Dol_per_Trip_PRE_P1",random_demos,random_campaigns,scoring_vars)]),collapse=" + ")," ")," + offset(log(Dol_per_Trip_PRE_P1+1))")
model_glm <- glm(Formula_fixed,data=finaldata_dolocc,family=Gamma(link="log"))
summary(model_glm)
rm(model_glm_condition1)

# p-values and signs checks
neutralvars <- setdiff(names(finaldata_dolocc),c(positivevars,negativevars))
pvalues <- summary(model_glm)$coefficients[,4][-1]
keep_covariate <- names(pvalues[pvalues<pvalue_lvl])
coeffs <- sign(summary(model_glm)$coefficients[,1])
varstokeep <- c();initialvars <- c()
for (i in negativevars) if (i %in% names(coeffs[coeffs<0])) varstokeep <- c(varstokeep,i)
    neg <- varstokeep
for (j in positivevars) if (j %in% names(coeffs[coeffs>0])) varstokeep <- c(varstokeep,j)
    pos <- setdiff(varstokeep,neg)
for (l in neutralvars) varstokeep <- c(varstokeep,l)
for (k in varstokeep) if (k %in% keep_covariate) initialvars <- c(initialvars,k)
#if(dolocc_mandatory_covariate %in% names(finaldata_dolocc) & length(setdiff(dolocc_mandatory_covariate,initialvars))>0){
#  initialvars <- c(dolocc_mandatory_covariate,initialvars)
#}

"""
end

#r=Symbol(lowercase(rcopy("initialvars")))

function chkSigns(m::MModel, vars::Array{Symbol}, dfd::DataFrame, cfg::OrderedDict)  # Pvalue & Signs
    g = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )
    neutralvars = setdiff(vars,vcat(cfg[:negativevars],cfg[:positivevars]))
    #keep_covariate = g.sdf[ (g.sdf[:pval].<cfg[:pvalue_lvl]) & (g.sdf[:vars].!=:intercept),:vars]
    #QC: sort(keep_covariate)==sort(rlst("keep_covariate"))
    neg=intersect(cfg[:negativevars],g.sdf[g.sdf[:coef].<0,:vars])
    #QC: sort(neg)==sort(rlst("neg"))
    #QC: cmparr(neg,rlst("neg"))
    pos=intersect(cfg[:positivevars],g.sdf[g.sdf[:coef].>0,:vars])
    #QC: sort(pos)==sort(rlst("pos"))
    varstokeep = intersect(vcat(neutralvars, pos,neg) ,  g.sdf[ g.sdf[:pval].<cfg[:pvalue_lvl] ,:vars] )
    #m.glm3_PnSigns_x= setdiff(vars, varstokeep )
    #vars = rmVars(vcat(m.glm3_PnSigns_x, g3.xvars) )
    #QC: sort(convert(Array{Symbol},varstokeep)) == sort(rlst("initialvars"))
    return g, varstokeep
end

(m.glm3_PnSigns, initialvars) = chkSigns(m, vars, dfd, cfg)
#sort(initialvars)==rlst("initialvars")

#    # Pvalue & Signs
#    m.glm3_PnSigns = xGLM(dfd, m.dist, m.y_var, m.logvar, vars  )  # Pvalue & Signs
#    g3=m.glm3_PnSigns
#    neutralvars = setdiff(vars,vcat(cfg[:negativevars],cfg[:positivevars]))
#keep_covariate = g3.sdf[ (g3.sdf[:pval].<cfg[:pvalue_lvl]) & (g3.sdf[:vars].!=:intercept),:vars]
##QC: sort(keep_covariate)==sort(rlst("keep_covariate"))
#    neg=intersect(cfg[:negativevars],g3.sdf[ g3.sdf[:coef].<0,:vars])
##QC: sort(neg)==sort(rlst("neg"))
#    pos=intersect(cfg[:positivevars],g3.sdf[ g3.sdf[:coef].>0,:vars])
##QC: sort(pos)==sort(rlst("pos"))
#    varstokeep = intersect(vcat(neutralvars, pos,neg) ,  g3.sdf[ g3.sdf[:pval].<cfg[:pvalue_lvl] ,:vars] )
#    #m.glm3_PnSigns_x= setdiff(vars, varstokeep )
#    #vars = rmVars(vcat(m.glm3_PnSigns_x, g3.xvars) )
##QC: sort(convert(Array{Symbol},varstokeep)) == sort(rlst("initialvars"))

#"""
## QC Varstokeep
#Symbol(lowercase(rcopy("initialvars"))) == varstokeep[1]
##rvars=rlst("varstokeep")
#sort(cfg[:positivevars]) == sort(rlst("positivevars"))
#sort(cfg[:negativevars]) == sort(rlst("negativevars"))
#cmparr(rvs,vcat(cfg[:random_demos],cfg[:random_campaigns],vars))
#sort(neutralvars) == rlst("neutralvars")
#sort(varstokeep) == rlst("varstokeep")
## -------------
#"""

if  m.modelName == "occ"
#OCC
R"""
#set the formula for the second model after cleaning the first model
Formula_final <- paste0(paste0(y_var,"~ ",paste0(initialvars,collapse=" + ")," "),"+ group+ offset(log(Trps_PRE_P1+1))")

model_glm_2 <-  glm(Formula_final,data=finaldata_occ,family=poisson(link="log"))
summary(model_glm_2)
rm(model_glm)

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

"""
end
if  m.modelName == "dolocc"
#DOLOCC
R"""
#set the formula for the second model after cleaning the first model
Formula_fixed <- paste0(paste0(y_var,"~ ",paste0(initialvars,collapse=" + ")," ")," + group + offset(log(Dol_per_Trip_PRE_P1+1))")
model_glm_2 <-  glm(Formula_fixed,data=finaldata_dolocc,family=Gamma(link="log"))
summary(model_glm_2)
rm(model_glm)

# final cleaning of second model based on signs and pvalues
pvalues_2 <- summary(model_glm_2)$coefficients[,4][-1]
keep_covariate_2 <- names(pvalues_2[pvalues_2 < pvalue_lvl])
coeffs_2 <- sign(summary(model_glm_2)$coefficients[,1])
varstokeep <- c();vars_2 <- c()
for (i in negativevars) if (i %in% names(coeffs_2[coeffs_2 < 0])) varstokeep <- c(varstokeep,i)
for (j in positivevars) if (j %in% names(coeffs_2[coeffs_2 > 0])) varstokeep <- c(varstokeep,j)
for (l in neutralvars) varstokeep <- c(varstokeep,l)
for (k in varstokeep) if (k %in% keep_covariate_2) vars_2 <- c(vars_2,k)
"""
end
#r=Symbol(lowercase(rcopy("vars_2")))
(m.glm4_PnSignsClean, vars_2) = chkSigns(m, convert(Array{Symbol},vcat(initialvars,[:group])) , dfd, cfg)

R"""
# final correlation check (find pair vars with correlation greater then 0.8 and drop the one with bigger pvalue) #### Maybe 0.6!
if(length(vars_2)>1){
  finaldata_occ_num <- as.data.table(finaldata_occ[,names(finaldata_occ) %in% vars_2])
  t1 <- as.data.frame(finaldata_occ_num[,sapply(finaldata_occ_num,is.numeric)])
  finaldata_occ_num <- as.data.frame(finaldata_occ_num)
  finaldata_occ_num <- finaldata_occ_num[,c(which(t1==T))]
  temp_cor <- as.matrix(cor(finaldata_occ_num))
  temp_cor[lower.tri(temp_cor, diag = T)] <- NA
  temp_corIndex <- which(temp_cor > 0.8 |temp_cor < (-0.8), arr.ind = TRUE)
  cordata <- data.table(col1=colnames(temp_cor)[temp_corIndex[,1]],col2=colnames(temp_cor)[temp_corIndex[,2]])

  if(nrow(cordata)>0) {
    drop_cor_vars <- c()
    for (i in 1:nrow(cordata)) {
      if(cordata[i][,1,with=F]=="group"|cordata[i][,2,with=F]=="group"){
        if(cordata[i][,1,with=F]=="group") drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,2,with=F]))
        if(cordata[i][,2,with=F]=="group") drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,1,with=F]))
      } else {
        if (summary(model_glm_2)$coefficients[which(rownames(summary(model_glm_2)$coefficients)==as.character(cordata[i][,1,with=F])),4] > summary(model_glm_2)$coefficients[which(rownames(summary(model_glm_2)$coefficients)==as.character(cordata[i][,2,with=F])),4]) {
          drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,1,with=F]))
        } else {
          drop_cor_vars <- c(drop_cor_vars,as.character(cordata[i][,2,with=F]))
        }
      }
    }
    vars_2 <- setdiff(vars_2,drop_cor_vars)
  }
}

"""

#end

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

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

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

"""

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using ArgParse

function parse_commandline(args::AbstractVector{T} where T<:AbstractString)
	s = ArgParseSettings()
	@add_arg_table! s begin
		"CONFIG"
      help = "path to a TOML or JSON file with the parameter settings"
      required = true
      arg_type = String
    "--output-params-dump"
      help = "path where genereated simulation parameters (SimParams) are dumped"
      arg_type = String
		"--output-run-dump-prefix"
      help = "prefix that will be added to a filename were full output for each run should be written. Each subsequent trajectory i will be saved under  \$prefix_\$i.jld2"
      arg_type = String
    "--output-daily"
      help = "path where daily trajectories should be saved"
      arg_type = String
    "--output-summary"
      help = "path where summary should be saved"
      arg_type = String
	end
  parse_args(args, s)
end


text:
freq = collect(linspace(150, 250, 50))THz
neff_data = collect(linspace(3.0, 1.0, 50))
neff = OpticalAttr(freq, neff_data, "n_eff")
corefrac = OpticalAttr(1.0, "Core Fraction")
Aeff = OpticalAttr(1e-12m^2, "Aeff (m²)")
loss = OpticalAttr(0.0m^-1, "Loss (1/m)")
coupling = OpticalAttr(0.6, "Coupling")

mode_TM = Mode(neff, Aeff, loss, corefrac, coupling, :TM, true, true)
mode_TE = Mode(neff, Aeff, loss, corefrac, coupling, :TE, true, true)
mode_TM2 = Mode(neff, Aeff, loss, corefrac, coupling, :TM, true, true)

wg_simple = Waveguide(Silicon)
add_straight!(wg_simple, 1e-3m)
add_turn!(wg_simple, 90°, 1e-3m)
add_mode!(wg_simple, mode_TM)
add_mode!(wg_simple, mode_TE)
add_mode!(wg_simple, mode_TM2)
@test mode_TM.has_raman
@test !(mode_TM.has_interaction)
@test !(mode_TE.has_interaction)
@test !NonlinearOptics.interacts_with_other_mode(mode_TM)
@test NonlinearOptics.get_polarization(mode_TM) == 1
@test NonlinearOptics.get_polarization(mode_TE) == 2

add_interaction!(wg_simple, mode_TM, mode_TE)
@test mode_TM.has_interaction
@test mode_TE.has_interaction
@test NonlinearOptics.interacts_with_other_mode(mode_TM)
@test !NonlinearOptics.interacts_with_other_mode(mode_TM2)

# link waveguides


text:
"""
    mutable struct RoundRotorQuadratic <: Machine
        base_machine::RoundRotorMachine
        saturation_coeffs::Tuple{Float64, Float64}

4-states round-rotor synchronous machine with quadratic saturation:
IEEE Std 1110 §5.3.2 (Model 2.2). GENROU model in PSSE and PSLF.

# Arguments
- `base_machine::RoundRotorMachine`: Round Rotor Machine model.
- `saturation_coeffs::Tuple{Float64, Float64}``: Saturation coefficients for quadratic model.
"""
mutable struct RoundRotorQuadratic <: Machine
    base_machine::RoundRotorMachine
    saturation_coeffs::Tuple{Float64, Float64}
end
IS.@forward((RoundRotorQuadratic, :base_machine), RoundRotorMachine)

function RoundRotorQuadratic(
    R::Float64,
    Td0_p::Float64,
    Td0_pp::Float64,
    Tq0_p::Float64,
    Tq0_pp::Float64,
    Xd::Float64,
    Xq::Float64,
    Xd_p::Float64,
    Xq_p::Float64,
    Xd_pp::Float64,
    Xl::Float64,
    Se::Tuple{Float64, Float64},
)
    saturation_coeffs = get_quadratic_saturation(Se)
    RoundRotorQuadratic(
        RoundRotorMachine(
            R,
            Td0_p,
            Td0_pp,
            Tq0_p,
            Tq0_pp,
            Xd,
            Xq,
            Xd_p,
            Xq_p,
            Xd_pp,
            Xl,
            Se,
        ),
        saturation_coeffs,
    )
end

function RoundRotorQuadratic(;
    R,
    Td0_p,
    Td0_pp,
    Tq0_p,
    Tq0_pp,
    Xd,
    Xq,
    Xd_p,
    Xq_p,
    Xd_pp,
    Xl,
    Se,
)
    RoundRotorQuadratic(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se)
end

function RoundRotorQuadratic(::Nothing)
    RoundRotorQuadratic(;
        R = 0.0,
        Td0_p = 0.0,
        Td0_pp = 0.0,
        Tq0_p = 0.0,
        Tq0_pp = 0.0,
        Xd = 0.0,
        Xq = 0.0,
        Xd_p = 0.0,
        Xq_p = 0.0,
        Xd_pp = 0.0,
        Xl = 0.0,
        Se = (0.0, 0.0),
    )
end

get_base_machine(value::RoundRotorQuadratic) = value.base_machine
get_saturation_coeffs(value::RoundRotorQuadratic) = value.saturation_coeffs
set_base_machine!(value::RoundRotorQuadratic, val::RoundRotorMachine) =
    value.base_machine = val
set_saturation_coeffs!(value::RoundRotorQuadratic, val::Tuple{Float64, Float64}) =
    value.saturation_coeffs = val

function IS.deserialize_struct(::Type{RoundRotorQuadratic}, data::Dict)
    vals = IS.deserialize_to_dict(RoundRotorQuadratic, data)
    return RoundRotorQuadratic(vals[:base_machine], vals[:saturation_coeffs])
end


text:
"""
	ENU(bounds::Bounds{LLA},
		lla_ref::LLA = center(bounds),
		datum::Ellipsoid = WGS84)

Convert LLA Bounds to ENU

there's not an unambiguous conversion, but for now,
returning the minimum bounds that contain all points contained
by the input bounds
"""
function ENU(bounds::Bounds{LLA}, lla_ref::LLA = center(bounds), datum::Ellipsoid = WGS84)

    max_x = max_y = -Inf
    min_x = min_y = Inf

    xs = [bounds.min_x, bounds.max_x]
    ys = [bounds.min_y, bounds.max_y]
    if bounds.min_y < 0.0 < bounds.max_y
        push!(ys, 0.0)
    end
    ref_x = getX(lla_ref)
    if bounds.min_x < ref_x < bounds.max_x ||
       (bounds.min_x > bounds.max_x && !(bounds.min_x >= ref_x >= bounds.max_x))
        push!(xs, ref_x)
    end

    for x_lla in xs, y_lla in ys
        pt = ENU(LLA(y_lla, x_lla), lla_ref, datum)
        x, y = getX(pt), getY(pt)

        min_x, max_x = min(x, min_x), max(x, max_x)
        min_y, max_y = min(y, min_y), max(y, max_y)
    end

    return Bounds{ENU}(min_y, max_y, min_x, max_x)
end

"""
	center(bounds::Bounds{ENU})

Get Center Point of ENU Bounds Region
"""
function center(bounds::Bounds{ENU})
    x_mid = (bounds.min_x + bounds.max_x) / 2
    y_mid = (bounds.min_y + bounds.max_y) / 2

    return ENU(x_mid, y_mid)
end

"""
	center(bounds::Bounds{LLA})

Get Center Point of LLA Bounds Region
"""
function center(bounds::Bounds{LLA})
    x_mid = (bounds.min_x + bounds.max_x) / 2
    y_mid = (bounds.min_y + bounds.max_y) / 2

    if bounds.min_x > bounds.max_x
        x_mid = x_mid > 0 ? x_mid - 180 : x_mid + 180
    end

    return LLA(y_mid, x_mid)
end

"""
    inbounds(loc::ENU, bounds::Bounds{ENU})

Check Whether a location `loc` is within bounds `bounds`
"""
function inbounds(loc::ENU, bounds::Bounds{ENU})
    x, y = getX(loc), getY(loc)
    bounds.min_x <= x <= bounds.max_x &&
    bounds.min_y <= y <= bounds.max_y
end

"""
    inbounds(loc::LLA, bounds::Bounds{LLA})

Check whether a location `loc` is within bounds `bounds`
"""
function inbounds(loc::LLA, bounds::Bounds{LLA})
    x, y = getX(loc), getY(loc)
    min_x, max_x = bounds.min_x, bounds.max_x
    (min_x > max_x ? !(max_x < x < min_x) : min_x <= x <= max_x) &&
    bounds.min_y <= y <= bounds.max_y
end

"""
	onbounds(loc::T, bounds::Bounds{T}) where T<:Union{LLA,ENU}

Check whether a location `loc` is onbounds `bounds`
Works only for points that have passed the inbounds test
"""
function onbounds(loc::T, bounds::Bounds{T}) where T<:Union{LLA,ENU}
    x, y = getX(loc), getY(loc)
    x == bounds.min_x || x == bounds.max_x ||
    y == bounds.min_y || y == bounds.max_y
end

"""
	boundary_point(p1::T, p2::T, bounds::Bounds{T}) where T<:Union{LLA,ENU}

Find the closest point within bounds
Works only for points where inbounds(p1) != inbounds(p2)
"""
function boundary_point(p1::T, p2::T, bounds::Bounds{T}) where T<:Union{LLA,ENU}
    x1, y1 = getX(p1), getY(p1)
    x2, y2 = getX(p2), getY(p2)

    x, y = Inf, Inf

    if bounds.min_x >  bounds.max_x && x1*x2 < 0

        if x1 < bounds.min_x && x2 < bounds.max_x || x2 < bounds.min_x && x1 < bounds.max_x
            x = bounds.min_x
            y = y1 + (y2 - y1) * (bounds.min_x - x1) / (x2 - x1)
        elseif x1 > bounds.max_x && x2 > bounds.min_x || x2 > bounds.max_x && x1 > bounds.min_x
            x = bounds.max_x
            y = y1 + (y2 - y1) * (bounds.max_x - x1) / (x2 - x1)
        end

        p3 = T(XY(x, y))
        inbounds(p3, bounds) && return p3
    end

    # Move x to x bound if segment crosses boundary
    if x1 < bounds.min_x < x2 || x1 > bounds.min_x > x2
        x = bounds.min_x
        y = y1 + (y2 - y1) * (bounds.min_x - x1) / (x2 - x1)
    elseif x1 < bounds.max_x < x2 || x1 > bounds.max_x > x2
        x = bounds.max_x
        y = y1 + (y2 - y1) * (bounds.max_x - x1) / (x2 - x1)
    end

    p3 = T(XY(x, y))
    inbounds(p3, bounds) && return p3

    # Move y to y bound if segment crosses boundary
    if y1 < bounds.min_y < y2 || y1 > bounds.min_y > y2
        x = x1 + (x2 - x1) * (bounds.min_y - y1) / (y2 - y1)
        y = bounds.min_y
    elseif y1 < bounds.max_y < y2 || y1 > bounds.max_y > y2
        x = x1 + (x2 - x1) * (bounds.max_y - y1) / (y2 - y1)
        y = bounds.max_y
    end

    p3 = T(XY(x, y))
    inbounds(p3, bounds) && return p3

    error("Failed to find boundary point.")
end


text:
"""
    _add_exp_score!(ata_model::AbstractModel)

# Description

Add the expected score constraints as specified in the `settings_file`.
If the names of the expected score columns of the item bank are not provided in the `settings_file` they are computed as the item characteristic functions following the IRT model always specified in the `settings_file`.
It requires the [`start_ata`](#ATA.start_ata) step.

# Arguments

- **`ata_model::AbstractModel`** : Required. The model built with `start_ata()` and with settings loaded by [`start_ata`](#ATA.start_ata) function.

"""
function _add_exp_score!(ata_model::AbstractModel)
    T = ata_model.settings.T
    n_items = ata_model.settings.n_items
    n_groups = ata_model.settings.n_groups
    expected_score_var =
        [ata_model.constraints[t].expected_score.var for t = 1:ata_model.settings.T]
    ICF = Vector{Matrix{Float64}}(undef, T)
    for t = 1:T
        if expected_score_var[t] == Symbol("")
            ICF[t] = item_char(
                ata_model.settings.irt.parameters,
                ata_model.constraints[t].expected_score.pts,
                model = ata_model.settings.irt.model,
                parametrization = ata_model.settings.irt.parametrization,
                D = ata_model.settings.irt.D,
            )[1][
                :,
                :,
                1,
            ]# K[t] x I
        else
            ICF[t] = zeros(Float64, 1, n_items)
            ICF[t][1, :] = ata_model.settings.bank[!, expected_score_var[t]]
        end
    end
    for t = 1:T
        ata_model.constraints[t].expected_score.val = ICF[t]
    end
    success!(ata_model, "Expected scores (ICFs) based on IRT parameters computed.")
    # for t = 1:ata_model.settings.T
    # 	DelimitedFiles.writedlm("opt/A_$t.csv", ata_model.constraints[t].constr_A)
    # 	DelimitedFiles.writedlm("opt/b_$t.csv", ata_model.constraints[t].constr_b)
    # end
    JLD2.@save "opt/ICF.jld2" ICF
    success!(ata_model, "Expected Score constrained.")
    return nothing
end


text:
module Distributions

using NumericExtensions
using Stats

export
    # types
    VariateForm,
    ValueSupport,
    Univariate,
    Multivariate,
    Matrixvariate,
    Discrete,
    Continuous,
    Distribution,
    UnivariateDistribution,
    MultivariateDistribution,
    MatrixDistribution,
    NonMatrixDistribution,
    DiscreteDistribution,
    ContinuousDistribution,
    DiscreteUnivariateDistribution,
    DiscreteMultivariateDistribution,
    DiscreteMatrixDistribution,
    ContinuousUnivariateDistribution,
    ContinuousMultivariateDistribution,
    ContinuousMatrixDistribution,
    SufficientStats,
    Arcsine,
    Bernoulli,
    Beta,
    BetaPrime,
    Binomial,
    Categorical,
    Cauchy,
    Chi,
    Chisq,
    Cosine,
    Dirichlet,
    DiscreteUniform,
    DoubleExponential,
    EdgeworthMean,
    EdgeworthSum,
    EdgeworthZ,
    EmpiricalUnivariateDistribution,
    Erlang,
    Exponential,
    FDist,
    Gamma,
    Geometric,
    Gumbel,
    HyperGeometric,
    InverseWishart,
    InverseGaussian,
    InvertedGamma,
    Kolmogorov,
    KSDist,
    KSOneSided,
    Laplace,
    Levy,
    Logistic,
    LogNormal,
    MixtureModel,
    Multinomial,
    MultivariateNormal,
    MvNormal,
    NegativeBinomial,
    NoncentralBeta,
    NoncentralChisq,
    NoncentralF,
    NoncentralT,
    Normal,
    Pareto,
    Poisson,
    Rayleigh,
    Skellam,
    TDist,
    Triangular,
    Truncated,
    Uniform,
    Weibull,
    Wishart,
    QQPair,

    # methods
    binaryentropy, # entropy of distribution in bits
    ccdf,          # complementary cdf, i.e. 1 - cdf
    cdf,           # cumulative distribution function
    cf,            # characteristic function
    cgf,           # cumulant generating function
    cquantile,     # complementary quantile (i.e. using prob in right hand tail)
    cumulant,      # cumulants of distribution
    dim,           # sample dimension of multivariate distribution
    entropy,       # entropy of distribution in nats
    fit,           # fit a distribution to data (using default method)
    fit_mle,       # fit a distribution to data using MLE
    fit_mle!,      # fit a distribution to data using MLE (inplace update to initial guess)
    fit_map,       # fit a distribution to data using MAP
    freecumulant,  # free cumulants of distribution
    insupport,     # predicate, is x in the support of the distribution?
    invlogccdf,    # complementary quantile based on log probability
    invlogcdf,     # quantile based on log probability
    isplatykurtic, # Is excess kurtosis > 0.0?
    isleptokurtic, # Is excess kurtosis < 0.0?
    ismesokurtic,  # Is excess kurtosis = 0.0?
    isprobvec,     # Is a probability vector?
    kde,           # Kernel density estimator
    kurtosis,      # kurtosis of the distribution
    logccdf,       # ccdf returning log-probability
    logcdf,        # cdf returning log-probability
    loglikelihood, # log probability of array of IID draws
    logpdf,        # log probability density
    logpdf!,       # evaluate log pdf to provided storage
    logpmf,        # log probability mass
    logpmf!,       # evaluate log pmf to provided storage
    posterior,     # Bayesian updating
    scale,         # scale parameter of a distribution
    rate,          # rate parameter of a distribution
    sqmahal,       # squared Mahalanobis distance to Gaussian center
    sqmahal!,      # inplace evaluation of sqmahal
    mean,          # mean of distribution
    median,        # median of distribution
    mgf,           # moment generating function
    mode,          # the mode of a unimodal distribution
    modes,         # mode(s) of distribution as vector
    moment,        # moments of distribution
    nsamples,      # get the number of samples in a data array based on distribution types
    pdf,           # probability density function (ContinuousDistribution)
    pmf,           # probability mass function (DiscreteDistribution)
    quantile,      # inverse of cdf (defined for p in (0,1))
    qqbuild,       # build a paired quantiles data structure for qqplots
    rand,          # random sampler
    rand!,         # replacement random sampler
    sample,        # sample from a source array
    sampler,       # create a Sampler object for efficient samples
    skewness,      # skewness of the distribution
    sprand,        # random sampler for sparse matrices
    std,           # standard deviation of distribution
    suffstats,     # compute sufficient statistics
    var,           # variance of distribution
    wsample        # weighted sampling from a source array

import Base.mean, Base.median, Base.quantile, Base.max, Base.min, Base.scale
import Base.Random, Base.rand, Base.rand!, Base.std, Base.var, Base.cor, Base.cov
import Base.show, Base.sprand
import NumericExtensions.dim, NumericExtensions.entropy
import Stats.kurtosis, Stats.skewness, Stats.mode, Stats.modes

#### Distribution type system

abstract ValueSupport
type Discrete <: ValueSupport end
type Continuous <: ValueSupport end

abstract VariateForm
type Univariate <: VariateForm end
type Multivariate <: VariateForm end
type Matrixvariate <: VariateForm end

abstract Distribution{F<:VariateForm,S<:ValueSupport}

typealias UnivariateDistribution{S<:ValueSupport}   Distribution{Univariate,S}
typealias MultivariateDistribution{S<:ValueSupport} Distribution{Multivariate,S}
typealias MatrixDistribution{S<:ValueSupport}       Distribution{Matrixvariate,S}
typealias NonMatrixDistribution Union(UnivariateDistribution, MultivariateDistribution)

typealias DiscreteDistribution{F<:VariateForm}   Distribution{F,Discrete}
typealias ContinuousDistribution{F<:VariateForm} Distribution{F,Continuous}

typealias DiscreteUnivariateDistribution     Distribution{Univariate,    Discrete}
typealias ContinuousUnivariateDistribution   Distribution{Univariate,    Continuous}
typealias DiscreteMultivariateDistribution   Distribution{Multivariate,  Discrete}
typealias ContinuousMultivariateDistribution Distribution{Multivariate,  Continuous}
typealias DiscreteMatrixDistribution         Distribution{Matrixvariate, Discrete}
typealias ContinuousMatrixDistribution       Distribution{Matrixvariate, Continuous}

abstract SufficientStats
abstract GenerativeFormulation

#### Include files

include("constants.jl")

include("fallbacks.jl")
include("rmath.jl")
include("specialfuns.jl")
include("tvpack.jl")
include("utils.jl")

include(joinpath("samplers", "categorical_samplers.jl"))

# Univariate distributions
include(joinpath("univariate", "arcsine.jl"))
include(joinpath("univariate", "bernoulli.jl"))
include(joinpath("univariate", "beta.jl"))
include(joinpath("univariate", "betaprime.jl"))
include(joinpath("univariate", "binomial.jl"))
include(joinpath("univariate", "categorical.jl"))
include(joinpath("univariate", "cauchy.jl"))
include(joinpath("univariate", "chi.jl"))
include(joinpath("univariate", "chisq.jl"))
include(joinpath("univariate", "cosine.jl"))
include(joinpath("univariate", "discreteuniform.jl"))
include(joinpath("univariate", "empirical.jl"))
include(joinpath("univariate", "exponential.jl"))
include(joinpath("univariate", "fdist.jl"))
include(joinpath("univariate", "gamma.jl"))
include(joinpath("univariate", "edgeworth.jl"))
include(joinpath("univariate", "erlang.jl"))
include(joinpath("univariate", "geometric.jl"))
include(joinpath("univariate", "gumbel.jl"))
include(joinpath("univariate", "hypergeometric.jl"))
include(joinpath("univariate", "inversegaussian.jl"))
include(joinpath("univariate", "invertedgamma.jl"))
include(joinpath("univariate", "kolmogorov.jl"))
include(joinpath("univariate", "ksdist.jl"))
include(joinpath("univariate", "ksonesided.jl"))
include(joinpath("univariate", "laplace.jl"))
include(joinpath("univariate", "levy.jl"))
include(joinpath("univariate", "logistic.jl"))
include(joinpath("univariate", "lognormal.jl"))
include(joinpath("univariate", "negativebinomial.jl"))
include(joinpath("univariate", "noncentralbeta.jl"))
include(joinpath("univariate", "noncentralchisq.jl"))
include(joinpath("univariate", "noncentralf.jl"))
include(joinpath("univariate", "noncentralt.jl"))
include(joinpath("univariate", "normal.jl"))
include(joinpath("univariate", "pareto.jl"))
include(joinpath("univariate", "poisson.jl"))
include(joinpath("univariate", "rayleigh.jl"))
include(joinpath("univariate", "skellam.jl"))
include(joinpath("univariate", "tdist.jl"))
include(joinpath("univariate", "triangular.jl"))
include(joinpath("univariate", "uniform.jl"))
include(joinpath("univariate", "weibull.jl"))

# Multivariate distributions
include(joinpath("multivariate", "dirichlet.jl"))
include(joinpath("multivariate", "multinomial.jl"))
include(joinpath("multivariate", "multivariatenormal.jl"))

# Matrix distributions
include(joinpath("matrix", "inversewishart.jl"))
include(joinpath("matrix", "wishart.jl"))

# Truncated distributions
include("truncate.jl")
include(joinpath("univariate", "truncated", "normal.jl"))

# Mixture distributions
include("mixturemodel.jl")

# Sample w/ and w/o replacement
include("sample.jl")

# Link functions for GLM's
#include("glmtools.jl")

# REPL representations
include("show.jl")

# Kernel density estimators
include("kde.jl")

# Expectations, entropy, KL divergence
include("functionals.jl")

include("conjugates.jl")

include("qq.jl")

end # module


text:
"""
kl\\_bicgstab( x0, b, atv, V, eta, ptv = nothing;
      kl_store=nothing, side = "right", lmaxit = 10, pdata = nothing)

BiCGSTAB linear solver. Deals with preconditioning.
Uses bicgstab\\_base with is oblivious to that.

The code works and does what it needs to do, but ...\n
The user interface is unstable and, even worse,
the nonlinear solvers are not hooked up yet.

The way this works it

Input:\n
x0:  initial iterate, this is usually zero for nonlinear solvers

b: right hand side (duh!)

atv:  matrix-vector product which depends on precomputed data pdta
      I expect you to use pdata most or all of the time, so it is not
      an optional argument, even if it's nothing (at least for now).
      If your mat-vec is just A*v, you have to write a function where
      A is the precomputed data.
      API for atv is av=atv(v,pdata)

V: a vector for me to store a Jacobian-vector product. It goes where
   FPS would go in gmres. You are best served if V is Float64.

eta: Termination happens when ||b - Ax|| <= eta || b ||

ptv:  preconditioner-vector product, which will also use pdata. The
      default is nothing, which is no preconditioning at all.
      API for ptv is px=ptv(x,pdat) just like kl\\_gmres

Keyword arguments

kl\\_store: You have the option of giving me some room
         for the vectors bicgstab needs to do its work. These in
         which I will not overwrite and a couple of vectors I use
         in the iteration. If you're only doing a linear solve, it
         does no harm to let me allocate those vectores in kl\\_bicgstab.
         The way to preallocate is ```kl_store=kstore(n,"bicgstab")``` where n
         is the number of unknows. I call this myself in the initialization
         phase if you don't do it ahead of me.

side: left or right preconditioning. The default is "right".

lmaxit: maximum number of linear iterations. The default is 10.

pdata: precomputed data. The default is nothing, but that ain't gonna
        work well for nonlinear equations.

Output:\n
A named tuple (sol, reshist, lits, idid)

where

sol= final result
reshist = residual norm history
lits = number of iterations
idid = status of the iteration
       true -> converged
       false -> failed to converge

### Examples from the docstrings for kl\\_bicgstab

In these examples you have the matrix and use
```
function atv(x, A)
    return A * x
end
```
to compute the matvec.

#### Three dimensional problem.

Will converge in the four iterations (worse than kl_gmres)

```jldoctest
julia> function atv(x, A)
           return A * x
       end
atv (generic function with 1 method)

julia> A = [0.001 0 0; 0 0.0011 0; 0 0 1.e4];

julia> V = zeros(3); b = [1.0; 1.0; 1.0]; x0 = zeros(3);

julia> gout.reshist
5-element Vector{Any}:
 1.73205e+00
 1.41421e+00
 3.21642e-03
 3.20321e-03
 4.98049e-13

julia> norm(b - A*gout.sol,Inf)
3.68594e-13
```

#### Integral equation. Notice that pdata has the kernel of the
operator and we do the matvec directly. Just like the previous example.
We put the grid information and, for this artifical example, the solution
in the precomputed data.

```jldoctest
julia> function integop(u, pdata)
                  K = pdata.K
                  return u - K * u
              end
integop (generic function with 1 method)

julia> function integopinit(n)
                  h = 1 / n
                  X = collect(0.5*h:h:1.0-0.5*h)
                  K = [ker(x, y) for x in X, y in X]
                  K .*= h
                  sol = [usol(x) for x in X]
                  f = sol - K * sol
                  pdata = (K = K, xe = sol, f = f)
                  return pdata
              end
integopinit (generic function with 1 method)

julia> function usol(x)
                  return exp.(x) .* log.(2.0 * x .+ 1.0)
              end
usol (generic function with 1 method)

julia> function ker(x, y)
                  ker = 0.1 * sin(x + exp(y))
              end
ker (generic function with 1 method)

julia> n=100; pdata = integopinit(n); ue = pdata.xe; f=pdata.f;

julia> u0 = zeros(size(f)); V = zeros(size(f));

julia> gout.reshist
4-element Vector{Any}:
 1.48252e+01
 2.90538e-02
 2.07823e-07
 2.17107e-17

julia> norm(gout.sol-ue,Inf)
8.88178e-16
```

"""
function kl_bicgstab( x0, b, atv, V, eta, ptv = nothing;
      kl_store=nothing, side = "right", lmaxit = 10, pdata = nothing)
#
# If you give me too much storage, I will fix it for you.
#
    isa(V,Vector) ? rhs = V : rhs = @view V[:,1]
    rhs .= b
    n=length(b)
#
# If you're playing with both gmres and bicgstab you might have
# lmaxit set to -1. That will break things to I fixed that.
#
    (lmaxit == -1) && (lmaxit = 10)
    if side == "right" || ptv == nothing
        itsleft = false
    else
        itsleft = true
        rhs .= ptv(rhs, pdata)
    end
    n=length(x0)
    kl_store=kstore(n,"bicgstab")
    linsol = kl_store[1]
#    linsol .= b
    y0 = kl_store[2]
    y0 .= x0
#    linsol = copy(b)
#    y0 = copy(x0)
    Kpdata = (pdata = pdata, side = side, ptv = ptv,
              atv = atv, linsol = linsol,kl_store=kl_store)
    bout = bicgstab_base(y0, rhs, Katv, eta;
              lmaxit = lmaxit, pdata = Kpdata)
    if side == "left" || ptv == nothing
        return bout
    else
        sol = y0
        sol .= ptv(sol, pdata)
        return (sol = sol, reshist = bout.reshist, lits=bout.lits,
             idid=bout.idid)
    end
end

"""
bicgstab_base(x0, rhs, atv, eta;
         lmaxit = 10, pdata = nothing)

Base BiCGSTAB. Overwrites initial iterate and right hand side.
"""
function bicgstab_base(x0, rhs, atv, eta; lmaxit = 10, pdata = nothing)
    r = rhs
    x = x0
    (norm(x0) == 0.0) || (r .-= atv(x0, pdata))
    k = 0
    rho = zeros(lmaxit + 2)
    rho[1] = 1.0
    rho[2] = r' * r
    alpha = 1.0
    omega = 1.0
    r0 = copy(r)
    rnorm = norm(r0)
    kl_store = pdata.kl_store
    v = kl_store[3]
    p = kl_store[4]
    s = kl_store[5]
    t = kl_store[6]
#    v = zeros(size(x0))
#    p = zeros(size(x0))
#    s = zeros(size(x0))
#    t = zeros(size(x0))
    tol = eta * norm(rhs)
    k = 0
    reshist = []
    push!(reshist, rnorm)
    idid=true
    while rnorm > tol && k < lmaxit
        k += 1
        abs(omega) > 0 || (println("Breakdown omega = 0"); break)
        beta = (rho[k+1] / rho[k]) * (alpha / omega)
        axpy!(-omega, v, p)
        #        p .= r + beta * (p - omega * v)
        #        p .= r + beta * p
        axpby!(1.0, r, beta, p)
        v .= atv(p, pdata)
        tau=r0'*v
        abs(tau) > 0 || (println("Breakdown r0'*v = 0 "); break)
        alpha = rho[k+1] / tau
        #        s .= r - alpha * v
        copy!(s, r)
        axpy!(-alpha, v, s)
        t .= atv(s, pdata)
        norm(t) > 0 || (println("Breakdown t = 0"); break)
        omega = (t' * s) / (t' * t)
        rho[k+2] = -omega * (r0' * t)
        #        r .= s - omega * t
        copy!(r, s)
        axpy!(-omega, t, r)
        #        x .= x + alpha * p + omega * s
        copy!(t, s)
        axpby!(alpha, p, omega, t)
        #        x .= x + t
        x .+= t
        rnorm = norm(r)
        push!(reshist, rnorm)
    end
    (rnorm <= tol) || (idid=false)
    return (sol = x, reshist = reshist, lits=k, idid=idid)
end


text:
#
# Contains APIs for an application client to interact with the resource manager.
# The underlying methods are implementations of the following protobuf services:
# - applicationclient_protocol.proto
# - application_history_client.proto

@doc doc"""
YarnClient holds a connection to the Yarn Resource Manager and provides
APIs for application clients to interact with Yarn.
""" ->
type YarnClient
    rm_conn::YarnClientProtocol

    function YarnClient(host::AbstractString, port::Integer, ugi::UserGroupInformation=UserGroupInformation())
        new(YarnClientProtocol(host, port, ugi))
    end
end

show(io::IO, client::YarnClient) = show(io, client.rm_conn)

function nodecount(client::YarnClient)
    inp = GetClusterMetricsRequestProto()
    resp = getClusterMetrics(client.rm_conn, inp)
    isfilled(resp, :cluster_metrics) ? resp.cluster_metrics.num_node_managers : 0
end

function nodes(client::YarnClient; all::Bool=false, nodelist::YarnNodes=YarnNodes(client.rm_conn.channel.ugi))
    inp = GetClusterNodesRequestProto()
    if all
        set_field!(inp, :nodeStates, [NodeStateProto.NS_NEW, NodeStateProto.NS_RUNNING, NodeStateProto.NS_UNHEALTHY, NodeStateProto.NS_DECOMMISSIONED, NodeStateProto.NS_LOST, NodeStateProto.NS_REBOOTED])
    else
        set_field!(inp, :nodeStates, [NodeStateProto.NS_RUNNING])
    end
    resp = getClusterNodes(client.rm_conn, inp)
    update(nodelist, resp)
    nodelist
end

@doc doc"""
YarnAppStatus wraps the protobuf type for ease of use
""" ->
type YarnAppStatus
    report::ApplicationReportProto
    function YarnAppStatus(report::ApplicationReportProto)
        new(report)
    end
end

function show(io::IO, status::YarnAppStatus)
    report = status.report

    if isfilled(report, :final_application_status) && report.final_application_status > 0
        final_state = "-$(FINAL_APP_STATES[report.final_application_status])"
    else
        final_state = ""
    end
    if isfilled(report, :progress)
        final_state *= "-$(report.progress)"
    end
    println(io, "YarnApp $(report.applicationType) ($(report.name)/$(report.applicationId.id)): $(APP_STATES[report.yarn_application_state])$(final_state)")
    println(io, "    location: $(report.user)@$(report.host):$(report.rpc_port)/$(report.queue)")
    if report.yarn_application_state > YarnApplicationStateProto.RUNNING
        println(io, "    time: $(report.startTime) to $(report.finishTime)")
        if isfilled(report, :app_resource_Usage)
            rusage = report.app_resource_Usage
            println(io, "    rusage:")
            println(io, "        mem,vcore seconds: $(rusage.memory_seconds), $(rusage.vcore_seconds)")
            println(io, "        containers: used $(rusage.num_used_containers), reserved $(rusage.num_reserved_containers)")
            println(io, "        mem: used $(rusage.used_resources.memory), reserved $(rusage.reserved_resources.memory), needed $(rusage.needed_resources.memory)")
            println(io, "        vcores: used $(rusage.used_resources.virtual_cores), reserved $(rusage.reserved_resources.virtual_cores), needed $(rusage.needed_resources.virtual_cores)")
        end
        if isfilled(report, :diagnostics)
            println(io, "    diagnostics: $(report.diagnostics)")
        end
    elseif report.yarn_application_state == YarnApplicationStateProto.RUNNING
        println(io, "    start time: $(report.startTime)")
    end
    nothing
end

# get am-rm token for an unmanaged app
am_rm_token(status::YarnAppStatus) = status.report.am_rm_token

@doc doc"""
YarnAppAttemptStatus wraps the protobuf type for ease of use
""" ->
type YarnAppAttemptStatus
    report::ApplicationAttemptReportProto
    function YarnAppAttemptStatus(report::ApplicationAttemptReportProto)
        new(report)
    end
end

function show(io::IO, status::YarnAppAttemptStatus)
    report = status.report

    atmpt_id = report.application_attempt_id

    atmpt_str = "$(atmpt_id.application_id.id)"
    if isfilled(report, :am_container_id)
        atmpt_str *= "/$(report.am_container_id.id)"
    else
        atmpt_str *= "/-"
    end
    atmpt_str *= "/$(atmpt_id.attemptId)"

    println(io, "YarnAppAttempt $(atmpt_str): $(ATTEMPT_STATES[report.yarn_application_attempt_state])")
    println(io, "    location: $(report.host):$(report.rpc_port)")
    if isfilled(report, :diagnostics)
        println(io, "    diagnostics: $(report.diagnostics)")
    end
    nothing
end

@doc doc"""
YarnApp represents one instance of application running on the yarn cluster
""" ->
type YarnApp
    client::YarnClient
    appid::ApplicationIdProto
    status::Nullable{YarnAppStatus}
    attempts::Array{YarnAppAttemptStatus}

    function YarnApp(client::YarnClient, appid::ApplicationIdProto)
        new(client, appid, Nullable{YarnAppStatus}(), YarnAppAttemptStatus[])
    end
end

@doc doc"""
APP_STATES: enum value to state map. Used for converting state for display.
""" ->
const APP_STATES = [:new, :new_saving, :submitted, :accepted, :running, :finished, :failed, :killed]

@doc doc"""
FINAL_APP_STATES: enum value to state map. Used for converting state for display.
""" ->
const FINAL_APP_STATES = [:succeeded, :failed, :killed]

@doc doc"""
ATTEMPT_STATES: enum value to state map. Used for converting state for display.
""" ->
const ATTEMPT_STATES = [:new, :submitted, :scheduled, :scheduled, :allocated_saving, :allocated, :launched, :failed, :running, :finishing, :finished, :killed]

function show(io::IO, app::YarnApp)
    if isnull(app.status)
        println(io, "YarnApp: $(app.appid.id)")
    else
        show(io, get(app.status))
    end
    nothing
end

function _new_app(client::YarnClient)
    inp = GetNewApplicationRequestProto()
    resp = getNewApplication(client.rm_conn, inp)
    resp.application_id, resp.maximumCapability.memory, resp.maximumCapability.virtual_cores
end

function submit(client::YarnClient, container_spec::ContainerLaunchContextProto, mem::Integer, cores::Integer;
        priority::Int32=one(Int32), appname::AbstractString="EllyApp", queue::AbstractString="default", apptype::AbstractString="YARN",
        reuse::Bool=false, unmanaged::Bool=false)
    appid, maxmem, maxcores = _new_app(client)

    prio = protobuild(PriorityProto, Dict(:priority => priority))
    res = protobuild(ResourceProto, Dict(:memory => mem, :virtual_cores => cores))
    asc = protobuild(ApplicationSubmissionContextProto, Dict(:application_id => appid,
                :application_name => appname,
                :queue => queue,
                :priority => prio,
                :unmanaged_am => unmanaged,
                :am_container_spec => container_spec,
                :resource => res,
                :maxAppAttempts => 2,
                :applicationType => apptype,
                :keep_containers_across_application_attempts => reuse))

    inp = protobuild(SubmitApplicationRequestProto, Dict(:application_submission_context => asc))

    submitApplication(client.rm_conn, inp)
    YarnApp(client, appid)
end

function kill(app::YarnApp)
    client = app.client
    inp = protobuild(KillApplicationRequestProto, Dict(:application_id => app.appid))

    resp = forceKillApplication(client.rm_conn, inp)
    resp.is_kill_completed
end

function status(app::YarnApp, refresh::Bool=true)
    if refresh || isnull(app.status)
        client = app.client
        inp = protobuild(GetApplicationReportRequestProto, Dict(:application_id => app.appid))

        resp = getApplicationReport(client.rm_conn, inp)
        app.status = isfilled(resp.application_report) ?  Nullable(YarnAppStatus(resp.application_report)) : Nullable{YarnAppStatus}()
    end
    app.status
end

function wait_for_state(app::YarnApp, state::Int32, timeout_secs::Int=60)
    @logmsg("waiting for application to reach $(APP_STATES[state]) ($state) state")
    t1 = time() + timeout_secs
    finalstates = (YarnApplicationStateProto.KILLED, YarnApplicationStateProto.FAILED, YarnApplicationStateProto.FINISHED)
    isfinalstate = state in finalstates
    while time() < t1
        nst = status(app)
        if !isnull(nst)
            appstate = get(nst).report.yarn_application_state
            (appstate == state) && (return true)
            isfinalstate || ((appstate in finalstates) && (return false))
        end
        sleep(1)
    end
    false
end

# get am-rm token for an unmanaged app
function am_rm_token(app::YarnApp)
    wait_for_state(app, YarnApplicationStateProto.ACCEPTED) || throw(YarnException("Application was not accepted"))
    am_rm_token(get(status(app)))
end

function attempts(app::YarnApp, refresh::Bool=true)
    if refresh || isnull(app.attempts)
        client = app.client
        inp = protobuild(GetApplicationAttemptsRequestProto, Dict(:application_id => app.appid))

        resp = getApplicationAttempts(client.rm_conn, inp)
        atmptlist = app.attempts
        empty!(atmptlist)
        if isfilled(resp.application_attempts)
            for atmpt in resp.application_attempts
                push!(atmptlist, YarnAppAttemptStatus(atmpt))
            end
        end
    end
    app.attempts
end

function wait_for_attempt_state(app::YarnApp, attemptid::Int32, state::Int32, timeout_secs::Int=60)
    @logmsg("waiting for application attempt $attemptid to reach $(ATTEMPT_STATES[state]) ($state) state")
    t1 = time() + timeout_secs
    finalstates = (YarnApplicationAttemptStateProto.APP_ATTEMPT_KILLED, YarnApplicationAttemptStateProto.APP_ATTEMPT_FAILED, YarnApplicationAttemptStateProto.APP_ATTEMPT_FINISHED)
    isfinalstate = state in finalstates
    while time() < t1
        atmptlist = attempts(app)
        for atmpt in atmptlist
            report = atmpt.report
            if report.application_attempt_id.attemptId == attemptid
                atmptstate = report.yarn_application_attempt_state
                (atmptstate == state) && (return true)
                isfinalstate || ((atmptstate in finalstates) && (return false))
                @logmsg("application attempt $attemptid is in state $(ATTEMPT_STATES[atmptstate]) ($atmptstate) state")
                break
            end
        end
        sleep(1)
    end
    false
end


text:
module LieAlgebrafields
    using Random
    using LinearAlgebra
    import ..Gaugefields:SU3GaugeFields,SU3GaugeFields_1d,SU2GaugeFields,SU2GaugeFields_1d,GaugeFields,
        GaugeFields_1d,make_staple_double!,substitute!,projlink!,set_wing!, evaluate_wilson_loops!,muladd!,
        SUNGaugeFields_1d,SUNGaugeFields,U1GaugeFields_1d,U1GaugeFields,evaluate_tensor_lines,SU
    import ..SUN_generator:Generator,lie2matrix!,matrix2lie!
    import ..Gaugefields
    import ..Actions:GaugeActionParam_autogenerator,SmearingParam_single,SmearingParam_multi,SmearingParam
    import ..Wilsonloops:make_plaq_staple_prime,make_plaq_staple

    abstract type LieAlgebraFields end

    struct SU3AlgebraFields <: LieAlgebraFields
        a::Array{Float64,5}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NC::Int64
        NumofBasis::Int64

        function SU3AlgebraFields(NC,NX,NY,NZ,NT)
            @assert NC == 3
            NumofBasis = 8
            return new(zeros(Float64,NumofBasis,NX,NY,NZ,NT),NX,NY,NZ,NT,NC,NumofBasis)
        end
    end

    struct SU2AlgebraFields <: LieAlgebraFields
        a::Array{Float64,5}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NC::Int64
        NumofBasis::Int64

        function SU2AlgebraFields(NC,NX,NY,NZ,NT)
            @assert NC == 2
            NumofBasis = 3
            return new(zeros(Float64,NumofBasis,NX,NY,NZ,NT),NX,NY,NZ,NT,NC,NumofBasis)
        end
    end

    struct SUNAlgebraFields <: LieAlgebraFields
        a::Array{Float64,5}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NC::Int64
        NumofBasis::Int64
        generators::Generator

        function SUNAlgebraFields(NC,NX,NY,NZ,NT)
            NumofBasis = NC^2-1
            generators = Generator(NC)
            return new(zeros(Float64,NumofBasis,NX,NY,NZ,NT),NX,NY,NZ,NT,NC,NumofBasis,generators)
        end
    end

    struct U1AlgebraFields <: LieAlgebraFields
        a::Array{Float64,5}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NC::Int64
        NumofBasis::Int64
        #generators::Generator

        function U1AlgebraFields(NC,NX,NY,NZ,NT)
            NumofBasis = 1
            #generators = Generator(NC)
            return new(zeros(Float64,NumofBasis,NX,NY,NZ,NT),NX,NY,NZ,NT,NC,NumofBasis)
        end
    end

    function Base.setindex!(x::LieAlgebraFields,v,i...)
        x.a[i...] = v
    end

    function Base.getindex(x::LieAlgebraFields,i...)
        return x.a[i...]
    end

    function LieAlgebraFields(NC,NX,NY,NZ,NT)
        if NC == 3
            return SU3AlgebraFields(NC,NX,NY,NZ,NT)
        elseif NC == 2
            return SU2AlgebraFields(NC,NX,NY,NZ,NT)
        elseif NC ≥ 4
            return SUNAlgebraFields(NC,NX,NY,NZ,NT)
        elseif NC == 1
            return U1AlgebraFields(NC,NX,NY,NZ,NT)
        end
    end

    function Base.similar(x::LieAlgebraFields)
        return LieAlgebraFields(x.NC,x.NX,x.NY,x.NZ,x.NT)
    end

    function Gaugefields.substitute!(x::LieAlgebraFields,pwork)
        n1,n2,n3,n4,n5 = size(x.a)
        #println(size(pwork))
        x.a[:,:,:,:,:] = reshape(pwork,(n1,n2,n3,n4,n5))
    end

    function Base.display(a::LieAlgebraFields)
        n1,n2,n3,n4,n5 = size(a.a)
        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            println("$i1 $i2 $i3 $i4 $i5 ",a.a[i1,i2,i3,i4,i5])# = α*a.a[i1,i2,i3,i4,i5]
                        end
                    end
                end
            end
        end
        return
    end

    function clear!(a::Array{T,1}) where T <: LieAlgebraFields
        for μ=1:4
            clear!(a[μ])
        end
        return
    end

    function clear!(a::LieAlgebraFields)
        n1,n2,n3,n4,n5 = size(a.a)
        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            a.a[i1,i2,i3,i4,i5] = 0
                        end
                    end
                end
            end
        end
        return
    end

    function LinearAlgebra.mul!(c::LieAlgebraFields,α::Number,a::LieAlgebraFields)
        #n1,n2,n3,n4,n5 = size(a.a)
        @. c.a = α*a.a

        return
        n1,n2,n3,n4,n5 = size(a.a)

        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            c.a[i1,i2,i3,i4,i5] = α*a.a[i1,i2,i3,i4,i5]
                        end
                    end
                end
            end
        end
        return
    end

    function Base.:*(x::Array{T,1},y::Array{T,1})  where T <: LieAlgebraFields
        s = 0
        for μ=1:4
            s += x[μ]*y[μ]
        end

        return s
    end

    function Base.:*(x::T,y::T) where T <: LieAlgebraFields
        n1,n2,n3,n4,n5 = size(x.a)
        s = 0
        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            s += x.a[i1,i2,i3,i4,i5]*y.a[i1,i2,i3,i4,i5]
                        end
                    end
                end
            end
        end
        return s
    end

    function Base.:*(x::SU3AlgebraFields,y::SU3AlgebraFields)
        n1,n2,n3,n4,n5 = size(x.a)
        s = 0
        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            s += x.a[i1,i2,i3,i4,i5]*y.a[i1,i2,i3,i4,i5]
                        end
                    end
                end
            end
        end
        return s
    end

    function Base.:*(x::SU2AlgebraFields,y::SU2AlgebraFields)
        n1,n2,n3,n4,n5 = size(x.a)
        s = 0
        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            s += x.a[i1,i2,i3,i4,i5]*y.a[i1,i2,i3,i4,i5]
                        end
                    end
                end
            end
        end
        return s
    end

    function add!(a::LieAlgebraFields,α,b::LieAlgebraFields)
        n1,n2,n3,n4,n5 = size(a.a)
        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            a.a[i1,i2,i3,i4,i5] += α*b.a[i1,i2,i3,i4,i5]
                        end
                    end
                end
            end
        end
        return
    end

    function gauss_distribution_lie!(p::Array{LieAlgebraFields,1})
        for i=1:length(p)
            gauss_distribution_lie!(p[i])
        end
    end

    function gauss_distribution_lie!(p::LieAlgebraFields)
        n1,n2,n3,n4,n5 = size(p.a)

        for i5=1:n5
            for i4=1:n4
                for i3=1:n3
                    for i2=1:n2
                        for i1=1:n1
                            p.a[i1,i2,i3,i4,i5] = randn()
                        end
                    end
                end
            end
        end

        return
    end

    const sr3 = sqrt(3)
    const sr3i = 1/sr3
    const sr3i2 = 2*sr3i

    function Gauge2Lie!(c::SU3AlgebraFields,x::SU3GaugeFields)
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        x11 = x[1,1,ix,iy,iz,it]
                        x12 = x[1,2,ix,iy,iz,it]
                        x13 = x[1,3,ix,iy,iz,it]
                        x21 = x[2,1,ix,iy,iz,it]
                        x22 = x[2,2,ix,iy,iz,it]
                        x23 = x[2,3,ix,iy,iz,it]
                        x31 = x[3,1,ix,iy,iz,it]
                        x32 = x[3,2,ix,iy,iz,it]
                        x33 = x[3,3,ix,iy,iz,it]

                        c[1,ix,iy,iz,it] = ( imag(x12) + imag(x21) )
                        c[2,ix,iy,iz,it] = ( real(x12) - real(x21) )
                        c[3,ix,iy,iz,it] = ( imag(x11) - imag(x22) )
                        c[4,ix,iy,iz,it] = ( imag(x13) + imag(x31) )
                        c[5,ix,iy,iz,it] = ( real(x13) - real(x31) )

                        c[6,ix,iy,iz,it] = ( imag(x23) + imag(x32) )
                        c[7,ix,iy,iz,it] = ( real(x23) - real(x32) )
                        c[8,ix,iy,iz,it] = sr3i *
                                ( imag(x11) + imag(x22) -
                                        2*imag(x33) )
                    end
                end

            end

        end

    end

    function Gauge2Lie!(c::SU2AlgebraFields,x::SU2GaugeFields)
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX

                        x11 = x[1,1,ix,iy,iz,it]
                        x12 = x[1,2,ix,iy,iz,it]
                        x21 = x[2,1,ix,iy,iz,it]
                        x22 = x[2,2,ix,iy,iz,it]

                        c[1,ix,iy,iz,it] = (imag(x12)+imag(x21))
                        c[2,ix,iy,iz,it] = (real(x12)-real(x21))
                        c[3,ix,iy,iz,it] = (imag(x11)-imag(x22))

                    end
                end

            end

        end

    end

    function Gauge2Lie!(c::SUNAlgebraFields,x::SUNGaugeFields)
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        g = c.generators
        NC = c.NC
        matrix = zeros(ComplexF64,NC,NC)
        a = zeros(ComplexF64,length(g))

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                matrix[k1,k2] =  x[k1,k2,ix,iy,iz,it]
                            end
                        end

                        matrix2lie!(a,g,matrix)
                        for k = 1:length(g)
                            c[k,ix,iy,iz,it] = 2*imag(a[k])
                        end

                    end
                end

            end

        end

    end

    function Gauge2Lie!(c::SU3AlgebraFields,x::SU3GaugeFields_1d)
        NX = c.NX
        NY = c.NY
        NZ = c.NZ
        NT = c.NT

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        i = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix

                        x11 = x[1,1,i]
                        x12 = x[1,2,i]
                        x13 = x[1,3,i]
                        x21 = x[2,1,i]
                        x22 = x[2,2,i]
                        x23 = x[2,3,i]
                        x31 = x[3,1,i]
                        x32 = x[3,2,i]
                        x33 = x[3,3,i]

                        c[1,ix,iy,iz,it] = ( imag(x12) + imag(x21) )
                        c[2,ix,iy,iz,it] = ( real(x12) - real(x21) )
                        c[3,ix,iy,iz,it] = ( imag(x11) - imag(x22) )
                        c[4,ix,iy,iz,it] = ( imag(x13) + imag(x31) )
                        c[5,ix,iy,iz,it] = ( real(x13) - real(x31) )

                        c[6,ix,iy,iz,it] = ( imag(x23) + imag(x32) )
                        c[7,ix,iy,iz,it] = ( real(x23) - real(x32) )
                        c[8,ix,iy,iz,it] = sr3i *
                                ( imag(x11) + imag(x22) -
                                        2*imag(x33) )

                    end
                end

            end

        end

    end

    function Gauge2Lie!(c::SU2AlgebraFields,x::SU2GaugeFields_1d)
        NX = c.NX
        NY = c.NY
        NZ = c.NZ
        NT = c.NT

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        i = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix

                        x11 = x[1,1,i]
                        x12 = x[1,2,i]
                        #x13 = x[1,3,i]
                        x21 = x[2,1,i]
                        x22 = x[2,2,i]

                        c[1,ix,iy,iz,it] = (imag(x12)+imag(x21))
                        c[2,ix,iy,iz,it] = (real(x12)-real(x21))
                        c[3,ix,iy,iz,it] = (imag(x11)-imag(x22))

                    end
                end

            end

        end

    end

    function Gauge2Lie!(c::SUNAlgebraFields,x::SUNGaugeFields_1d)
        NX = c.NX
        NY = c.NY
        NZ = c.NZ
        NT = c.NT

        g = c.generators
        NC = c.NC
        matrix = zeros(ComplexF64,NC,NC)
        a = zeros(ComplexF64,length(g))

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        i = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix

                        for k2=1:NC
                            for k1=1:NC
                                matrix[k1,k2] =  x[k1,k2,i]
                            end
                        end
                        matrix2lie!(a,g,matrix)
                        for k = 1:length(g)
                            c[k,ix,iy,iz,it] = 2*imag(a[k])
                        end

                    end
                end

            end

        end

    end

    function Gauge2Lie!(c::U1AlgebraFields,x::U1GaugeFields_1d)
        NX = c.NX
        NY = c.NY
        NZ = c.NZ
        NT = c.NT

        NC = c.NC

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        i = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix
                        c[1,ix,iy,iz,it] = real(x[1,1,i])
                    end
                end

            end

        end

    end

    function expF_U!(U::Array{T,1},F::Array{N,1},Δτ,temps::Array{T_1d,1},temp_a::N) where {T<: GaugeFields,N<: LieAlgebraFields,T_1d <: GaugeFields_1d}
        @assert Δτ != 0 "Δτ should not be zero in expF_U! function!"
        temp1 = temps[1]
        temp2 = temps[2]
        temp3 = temps[3]
        temp4 = temps[4]

        c = temp_a

        for μ=1:4
            substitute!(temp1,U[μ])

            mul!(c,Δτ,F[μ]) #c = Δτ*F

            expA!(temp2,c,temp3,temp4) #temp2 = exp(c)

            mul!(temp3,temp2,temp1) #temp3 = temp2*temp1

            substitute!(U[μ],temp3)
            set_wing!(U[μ])

        end
    end

    function add_gaugeforce!(F::Array{N,1},U::Array{T,1},temps::Array{T_1d,1},temp_a::N;fac = 1) where {T<: GaugeFields,N <: LieAlgebraFields,T_1d <: GaugeFields_1d}
        temp1 = temps[1]
        temp2 = temps[2]
        temp3 = temps[3]
        staple = temps[4]

        c = temp_a

        for μ=1:4
            #println(1)
            make_staple_double!(staple,U,μ,temp1,temp2,temp3)

            substitute!(temp1,U[μ])

            #=
            c
            c                   ---
            c                  |   |
            c            tmp2 = ---    +   ---
            c                             |   |
            c                              ---
            c
            =#

            mul!(temp2,temp1,staple')

            #.....   Projection onto Lie Algebra   .....

            projlink!(temp3,temp2)

            #display(temp3)
            #exit()

            Gauge2Lie!(c,temp3)

            #display(c)
            #exit()

            #c     .....   p(new) = p(old) + fac * c  .....
            add!(F[μ],fac,c)

        end
        return
    end

    function add_gaugeforce!(F::Array{N,1},U::Array{T,1},temps::Array{T_1d,1},temp_a::N,gparam::GaugeActionParam_autogenerator;fac = 1) where {T<: GaugeFields,N <: LieAlgebraFields,T_1d <: GaugeFields_1d}
        temp1 = temps[1]
        temp2 = temps[2]
        temp3 = temps[3]
        staple = temps[4]
        Utemp = temps[5]
        temp4 = temps[6]

        c = temp_a

        #=
        NX = U[1].NX
        NY = U[1].NY
        NZ = U[1].NZ
        NT = U[1].NT
        NC = U[1].NC
        V = zeros(ComplexF64,NC,NC)
        Vtemp = zeros(ComplexF64,NC,NC)
        =#

        for μ=1:4
            substitute!(Utemp,U[μ])
            Gaugefields.clear!(staple)
            for i=1:gparam.numactions
                evaluate_wilson_loops!(temp4,gparam.staples[i][μ],U,[temp1,temp2,temp3])
                muladd!(staple,gparam.βs[i]/gparam.β,temp4)

                #ix,iy,iz,it = 1,1,1,1
                #icum = (((it-1)*NZ+iz-1)*NY+iy-1)*NX+ix
                #evaluate_wilson_loops!(Vtemp,gparam.staples[i][μ],U,ix,iy,iz,it)
                #println("μ = $μ, i = $i")
                #println("all")
                #display(temp4[:,:,icum])
                #println("\t")
                #println("partial")
                #display(Vtemp)
                #println("\t")
            end
            #exit()

                #substitute!(temp1,U[μ])
                #mul!(temp2,temp1,staple)

            mul!(temp2,Utemp,staple)

            projlink!(temp3,temp2)
            Gauge2Lie!(c,temp3)
            add!(F[μ],fac,c)

        end
        return
    end

    const pi23 = 2pi/3
    const tinyvalue =1e-100

    function expA!(v::SU2GaugeFields_1d ,u::SU2AlgebraFields,temp1,temp2)
        NX=u.NX
        NY=u.NY
        NZ=u.NZ
        NT=u.NT
        NC=u.NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        icum = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix
                        u1 = u[1,ix,iy,iz,it]/2
                        u2 = u[2,ix,iy,iz,it]/2
                        u3 = u[3,ix,iy,iz,it]/2
                        R = sqrt(u1^2+u2^2+u3^2) +  tinyvalue
                        sR = sin(R)/R
                        #sR = ifelse(R == 0,1,sR)
                        a0 = cos(R)
                        a1 = u1*sR
                        a2 = u2*sR
                        a3 = u3*sR

                        v[1,1,icum] = cos(R) + im*a3
                        v[1,2,icum] = im*a1 + a2
                        v[2,1,icum] = im*a1 - a2
                        v[2,2,icum] = cos(R) - im*a3
                    end
                end
            end
        end

    end

    function expA!(v::SUNGaugeFields_1d ,u::SUNAlgebraFields,temp1,temp2)
        g = u.generators

        NX=u.NX
        NY=u.NY
        NZ=u.NZ
        NT=u.NT
        NC=u.NC
        u0 = zeros(ComplexF64,NC,NC)
        a = zeros(Float64,length(g))

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        icum = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix
                        for k=1:length(a)
                            a[k] = u[k,ix,iy,iz,it]
                        end

                        lie2matrix!(u0,g,a)
                        v[:,:,icum] = exp((im/2)*u0)
                    end
                end
            end
        end

    end

    function expA!(v::U1GaugeFields_1d ,u::U1AlgebraFields,temp1,temp2)

        NX=u.NX
        NY=u.NY
        NZ=u.NZ
        NT=u.NT
        NC=u.NC

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        icum = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix
                        v[1,1,icum] = exp(im*u[1,ix,iy,iz,it])
                    end
                end
            end
        end

    end

    """
    The descrition is based on the LTK in fortran.
c----------------------------------------------------------------c
c     Project Lie algebra u to Lie group element g.
c     U = u(1)*lambda(1)/2 + u(2)*lambda(2)/2) + ....
c     g = exp(i*U)
c     ( by exact diagonalization                           )
c     (  X = Adj(V) * D * V                                )
c     (  where D = diag(d1,d2,d3)i                         )
c     (  g = Adj(V) * diag(exp(i*d1),exp(i*d2),exp(i*d3))) )
c----------------------------------------------------------------c
c for i = 1 to n, v(i) is the SU(3) projection of Lie algebra element u(i)
c     !!!!!   This is for SU(3)   !!!!!
c----------------------------------------------------------------c
    """

    function expA!(v::SU3GaugeFields_1d ,u::SU3AlgebraFields,temp1,temp2)
        ww = temp1
        w = temp2

        #ww = similar(v)
        #w =similar(v)

        algemat!(v,u)

        for i=1:v.NV
            v1  = real(v[1,1,i])
            v2  = imag(v[1,1,i])
            v3  = real(v[1,2,i])
            v4  = imag(v[1,2,i])
            v5  = real(v[1,3,i])
            v6  = imag(v[1,3,i])
            v7  = real(v[2,1,i])
            v8  = imag(v[2,1,i])
            v9  = real(v[2,2,i])
            v10 = imag(v[2,2,i])
            v11 = real(v[2,3,i])
            v12 = imag(v[2,3,i])
            v13 = real(v[3,1,i])
            v14 = imag(v[3,1,i])
            v15 = real(v[3,2,i])
            v16 = imag(v[3,2,i])
            v17 = real(v[3,3,i])
            v18 = imag(v[3,3,i])

    #c find eigenvalues of v
            trv3 = (v1 + v9 + v17) / 3.0
            cofac = v1 *  v9  - v3^2  - v4^2 +
                      v1 * v17 -  v5^2  - v6^2 +
                      v9 * v17 - v11^2 - v12^2
            det = v1 * v9 * v17 -
                      v1 * (v11^2 + v12^2) -
                      v9 * (v5^2 + v6^2) -
                      v17 * (v3^2 + v4^2) +
                      (v5 * (v3 * v11 - v4 * v12) +
                       v6 * (v3 * v12 + v4 * v11)) * 2.0
            p3 = cofac / 3.0 - trv3^2
            q = trv3 * cofac - det - 2.0 * trv3^3
            x = sqrt(-4.0 * p3) + tinyvalue
            arg = q / (x * p3)

            arg = min(1, max(-1, arg))
            theta = acos(arg) / 3.0
            e1 = x * cos(theta) + trv3
            theta = theta + pi23
            e2 = x * cos(theta) + trv3
    #       theta = theta + pi23
    #       e3 = x * cos(theta) + trv3
            e3 = 3.0 * trv3 - e1 - e2

    # solve for eigenvectors

            w1 =   v5 * (v9 - e1) - v3 * v11 + v4 * v12
            w2 = - v6 * (v9 - e1) + v4 * v11 + v3 * v12
            w3 =   (v1 - e1) * v11 - v3 * v5 - v4 * v6
            w4 = - (v1 - e1) * v12 - v4 * v5 + v3 * v6
            w5 = - (v1 - e1) * (v9 - e1) +  v3^2 + v4^2
            w6 = 0.0

            coeff = 1.0 / sqrt(w1^2 + w2^2 + w3^2 + w4^2 + w5^2)

            w1 = w1 * coeff
            w2 = w2 * coeff
            w3 = w3 * coeff
            w4 = w4 * coeff
            w5 = w5 * coeff

            w7 =   v5 * (v9 - e2) - v3 * v11 + v4 * v12
            w8 = - v6 * (v9 - e2) + v4 * v11 + v3 * v12
            w9 = (v1 - e2) * v11 - v3 * v5 - v4 * v6
            w10 = - (v1 - e2) * v12 - v4 * v5 + v3 * v6
            w11 = - (v1 - e2) * (v9 - e2) +  v3^2 + v4^2
            w12 = 0.0

            coeff = 1.0 / sqrt(w7^2  + w8^2 + w9^2+ w10^2 + w11^2)
            w7  = w7  * coeff
            w8  = w8  * coeff
            w9  = w9  * coeff
            w10 = w10 * coeff
            w11 = w11 * coeff

            w13 =   v5 * (v9 - e3) - v3 * v11 + v4 * v12
            w14 = - v6 * (v9 - e3) + v4 * v11 + v3 * v12
            w15 =   (v1 - e3) * v11 - v3 * v5 - v4 * v6
            w16 = - (v1 - e3) * v12 - v4 * v5 + v3 * v6
            w17 = - (v1 - e3) * (v9 - e3) +  v3^2 + v4^2
            w18 = 0.0

            coeff = 1.0 / sqrt(w13^2 + w14^2 + w15^2+ w16^2 + w17^2)
            w13 = w13 * coeff
            w14 = w14 * coeff
            w15 = w15 * coeff
            w16 = w16 * coeff
            w17 = w17 * coeff

    # construct the projection v
            c1 = cos(e1)
            s1 = sin(e1)
            ww1  = w1  * c1 - w2 * s1
            ww2  = w2  * c1 + w1 * s1
            ww3  = w3  * c1 - w4 * s1
            ww4  = w4  * c1 + w3 * s1
            ww5  = w5  * c1 - w6 * s1
            ww6  = w6  * c1 + w5 * s1

            c2 = cos(e2)
            s2 = sin(e2)
            ww7  = w7  * c2 - w8 * s2
            ww8  = w8  * c2 + w7 * s2
            ww9  = w9  * c2 - w10 * s2
            ww10 = w10 * c2 + w9 * s2
            ww11 = w11 * c2 - w12 * s2
            ww12 = w12 * c2 + w11 * s2

            c3 = cos(e3)
            s3 = sin(e3)
            ww13 = w13 * c3 - w14 * s3
            ww14 = w14 * c3 + w13 * s3
            ww15 = w15 * c3 - w16 * s3
            ww16 = w16 * c3 + w15 * s3
            ww17 = w17 * c3 - w18 * s3
            ww18 = w18 * c3 + w17 * s3

            w[1,1,i]=   w1+im*w2
            w[1,2,i]=   w3+im*w4
            w[1,3,i]=   w5+im*w6
            w[2,1,i]=   w7+im*w8
            w[2,2,i]=   w9+im*w10
            w[2,3,i]=   w11+im*w12
            w[3,1,i]=   w13+im*w14
            w[3,2,i]=   w15+im*w16
            w[3,3,i]=   w17+im*w18

            ww[1,1,i]=   ww1+im*ww2
            ww[1,2,i]=   ww3+im*ww4
            ww[1,3,i]=   ww5+im*ww6
            ww[2,1,i]=   ww7+im*ww8
            ww[2,2,i]=   ww9+im*ww10
            ww[2,3,i]=   ww11+im*ww12
            ww[3,1,i]=   ww13+im*ww14
            ww[3,2,i]=   ww15+im*ww16
            ww[3,3,i]=   ww17+im*ww18
        end

        mul!(v,w',ww)
        return

    end

    const sr3 = sqrt(3)
    const sr3i = 1/sr3
    const sr3i2 = 2*sr3i

    function algemat!(x::SU3GaugeFields_1d,t::SU3AlgebraFields)
        NX=t.NX
        NY=t.NY
        NZ=t.NZ
        NT=t.NT
        NC=t.NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        icum = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix

                        c1 = t[1,ix,iy,iz,it] * 0.5
                        c2 = t[2,ix,iy,iz,it] * 0.5
                        c3 = t[3,ix,iy,iz,it] * 0.5
                        c4 = t[4,ix,iy,iz,it] * 0.5
                        c5 = t[5,ix,iy,iz,it] * 0.5
                        c6 = t[6,ix,iy,iz,it] * 0.5
                        c7 = t[7,ix,iy,iz,it] * 0.5
                        c8 = t[8,ix,iy,iz,it] * 0.5

                        x[1,1,icum] =  c3+sr3i*c8 +im*(  0.0 )
                        x[1,2,icum] =  c1         +im*( -c2   )
                        x[1,3,icum] =  c4         +im*(-c5   )

                        x[2,1,icum] =  c1         +im*(  c2   )
                        x[2,2,icum] =  -c3+sr3i*c8+im*(  0.0 )
                        x[2,3,icum] =  c6         +im*( -c7   )

                        x[3,1,icum] =  c4         +im*(  c5   )
                        x[3,2,icum] =  c6         +im*(  c7   )
                        x[3,3,icum] =  -sr3i2*c8  +im*(  0.0 )
                    end

                end

            end
        end
    end

    function stoutfource(dSdU,U,smearing::T)  where T <: SmearingParam_single
        dSdUnew = deepcopy(dSdU)
        dSdρs = stoutfource!(dSdUnew,dSdU,U,smearing,smearing.ρs)
        return dSdUnew,dSdρs
    end

    #evaluate_tensor_lines!(V,nu,dSdU::Array{T_1d,1},gparam::GaugeActionParam_autogenerator,U::Array{T,1},umu::Tuple{I,I},ix,iy,iz,it) where {T <: GaugeFields,I <: Int,T_1d <: GaugeFields_1d}
    function stoutfource!(dSdUnew,dSdU,U,smearing::T,ρs)  where T <: SmearingParam
        #nu = 1
        #evaluate_tensor_lines(nu,dSdU,smearing,U,(nu,1),ρs)

        #println("stout force!!")
        #dSdUnew = deepcopy(dSdU)

        dSdρs = zero(ρs)
        NX=U[1].NX
        NY=U[1].NY
        NZ=U[1].NZ
        NT=U[1].NT
        NC=U[1].NC
        V = zeros(ComplexF64,NC,NC)
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        icum = (((it-1)*NZ+iz-1)*NY+iy-1)*NX+ix
                        for nu=1:4

                            #println(ρs)
                            #V,dSdρs_site = evaluate_tensor_lines_2(nu,dSdU,smearing,U,(nu,1),ix,iy,iz,it,ρs)
                            V,dSdρs_site = evaluate_tensor_lines(nu,dSdU,smearing,U,(nu,1),ix,iy,iz,it,ρs)
                            #V = evaluate_tensor_lines_2(nu,dSdU,gparam,U,(nu,1),ρ,ix,iy,iz,it)
                            #V = evaluate_SU2_force(nu,dSdU,gparam,U,(nu,1),ρ,ix,iy,iz,it)
                            #V = evaluate_tensor_lines(nu,dSdU,gparam,U,(nu,1),ρ,ix,iy,iz,it)
                            #println("dSdU\t")
                            #println("VV", V)
                            for i=1:NC
                                for j=1:NC
                                    dSdUnew[nu][j,i,icum] = V[j,i]
                                end
                            end

                            @. dSdρs += dSdρs_site
                            #dSdUnew[nu][:,:,icum] = V2[:,:]

                             #=
                            if nu == 1
                                if (ix,iy,iz,it) == (1,1,1,1) || (ix,iy,iz,it) ==(1,2,2,2) || (ix,iy,iz,it) ==(2,2,2,2) || (ix,iy,iz,it) ==(4,4,4,4)
                                    println("x ",(ix,iy,iz,it))
                                # #=
                                    println("U ", U[nu][:,:,ix,iy,iz,it])

                                    println("dS[Ufat]/dUfat ")
                                    display(dSdU[nu][:,:,icum])
                                    println("\t")
                                    println("dS[Ufat]/dU ")
                                    display(dSdUnew[nu][:,:,icum])
                                    println("\t")
                                    println("new2")
                                    #display(V2)
                                    println("\t")
                                    println("old")
                                    display(V)
                                    println("\t")

                                    #=
                                    V = evaluate_tensor_lines(nu,dSdU,gparam,U,(nu,1),-ρ,ix,iy,iz,it)
                                    #V = evaluate_tensor_lines(nu,dSdU,gparam,U,(nu,1),-ρ,ix,iy,iz,it)
                                    println("new")
                                    display(V)
                                    println("\t")

                                    V2 = evaluate_tensor_lines_2(nu,dSdU,gparam,U,(nu,1),-ρ,ix,iy,iz,it)
                                    println("new2")
                                    display(V2)
                                    println("\t")
                                    =#
                                    exit()

                                 # =#
                                end
                            end
                             =#

                        end

                    end
                end
            end
        end
        #exit()

        return dSdρs
        #exit()
        #exit()
        #return dSdUnew
    end

    function stoutfource(dSdU,U_multi,U,smearing::T)  where T <: SmearingParam_multi
        dSdUnew = deepcopy(dSdU)
        dSdρs = stoutfource!(dSdUnew,dSdU,U_multi,U,smearing)
        return dSdUnew,dSdρs
    end

    function stoutfource!(dSdUnew,dSdU,U_multi,U::Array{GaugeFields{SU{NC}},1},smearing::T) where {NC,T <:  SmearingParam_multi}
        num = length(U_multi)
        dSdU_m = deepcopy(dSdU)
        dSdU_m2 = deepcopy(dSdU)
        dSdρs = deepcopy(smearing.ρs)
        if num == 1
            dSdρs[1] = stoutfource!(dSdUnew,dSdU,U,smearing,smearing.ρs[1])
        else
            for i=num:-1:2
                #println(i)
                #println(smearing.ρs[i])
                dSdρs[i] = stoutfource!(dSdU_m2,dSdU_m,U_multi[i-1],smearing,smearing.ρs[i])
                dSdU_m,dSdU_m2 = dSdU_m2,dSdU_m
            end
            dSdρs[1] =stoutfource!(dSdUnew,dSdU_m,U,smearing,smearing.ρs[1])
        end
        return dSdρs

    end

end


text:
using Test

m = test_page_model()
include("../src/components/LGemissions.jl")

add_comp!(m, LGemissions)

set_param!(m, :LGemissions, :e0_baselineLGemissions, readpagedata(m,"data/e0_baselineLGemissions.csv"))
set_param!(m, :LGemissions, :er_LGemissionsgrowth, readpagedata(m, "data/er_LGemissionsgrowth.csv"))

# run Model
run(m)

emissions= m[:LGemissions,  :e_regionalLGemissions]
emissions_compare=readpagedata(m, "test/validationdata/e_regionalLGemissions.csv")

@test emissions ≈ emissions_compare rtol=1e-3