text:
using DataStructures, DataFrames, StatsFuns, GLM , JuMP,NLopt, StatStack, NLsolve, DBAPI, JavaCall, JDBC

if !isfile("Occasion/Occ_fixed_effects.csv")
    println("\n\n Error: This code can only be run from the Scoring directory under the camapign.\n\n\n")
    quit()
end

"""
-------- app.cfg --------
excludedBreaks =
excludedKeys =
P2_Competitor = true
pvalue_lvl = 0.20
exposed_flag_var = new_exposed_flag
sigLevel = 0.2
random_demos =
random_campaigns =
dropvars =
scoring_vars =
TotalModelsOnly=false
"""

cfgDefaults=OrderedDict( :P2_Competitor => true
                        ,:pvalue_lvl => 0.20  #pvalue_lvl = 0.20
                        ,:excludedBreaks => AbstractString[]    #["estimated_hh_income","hh_age","number_of_children_in_living_un","person_1_gender"]
                        ,:excludedLevels => ["none"]
                        ,:excludedKeys => AbstractString[]
                        ,:exposed_flag_var => :exposed_flag
                        ,:sigLevel => "0.2"
                        ,:random_demos => Symbol[]
                        ,:random_campaigns => Symbol[]
                        ,:dropvars => Symbol[]
                        ,:scoring_vars => Symbol[]
                        ,:TotalModelsOnly=>false
                       )

cfg = loadCFG(cfgDefaults,pwd()*"/"*"app.cfg")
cfg[:counts] = isfile(pwd()*"/hhcounts.csv") ? true : false
if !cfg[:counts] println("\n\n WARNING: HH Count file does not exist. No weighting will be applied to the scoring.\n\n\n") end

# ----------------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------------------------------
# ------------------------------------------------------ SCORING -------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------------------------------------------

df_data = readtable("final_data.csv",header=true); lowercase(df_data)

mocc = MOcc(df_data,cfg)
mdolocc = MDolOcc(df_data,cfg)
mpen = MPen(df_data,cfg)
mdolhh=MDolHH(mocc,mdolocc,mpen)

cnts = Cnts(df_data, mocc,mdolocc,mpen,cfg[:TotalModelsOnly])

if cfg[:counts]
    hhcounts(cnts)
end

function aggregate!(df_data::DataFrame, focc::FOcc, cnts::Cnts)
    sdf = focc.sdf
    #sdf[:M] = cnts.get("Total Campaign")[:M][1]
    #sdf[:Mt] = cnts.get("Total Campaign")[:Mt][1]
    #sdf[:Mc] = cnts.get("Total Campaign")[:Mc][1]

    sdf[:M] = getk(cnts.sdf,"Total Campaign")[:M][1]   #cnts.get("Total Campaign")[:M][1]
    sdf[:Mt] = getk(cnts.sdf,"Total Campaign")[:Mt][1]   #cnts.get("Total Campaign")[:Mt][1]
    sdf[:Mc] = getk(cnts.sdf,"Total Campaign")[:Mc][1]   #cnts.get("Total Campaign")[:Mc][1]

    sdf[:mean_score0]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score0] )
    sdf[:mean_score1]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score1] )

    sdf[:unadj_mean_score0]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0),:occ_score0])
    sdf[:unadj_mean_score1]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1),:occ_score1])

    sdf[:adj_mean_cntrl_grp ] = mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score0] )
    sdf[:adj_mean_expsd_grp ] = mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :occ_score1] )
    # -------  Occ Score
    sdf[:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 0), :trps_pre_p1] )
    sdf[:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 1), :trps_pre_p1] )
    sdf[:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0), :trps_pos_p1] )
    sdf[:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1), :trps_pos_p1] )

    #A hack to fix missing data --
    sdf[isnan(sdf[:unadj_avg_expsd_hh_pre]),:unadj_avg_expsd_hh_pre] = 0.0

end
#aggregate!(df_data, mocc.feff,cnts)

function aggregate!(df_data::DataFrame, fdolocc::FDolOcc, cnts::Cnts)
    sdf = fdolocc.sdf
    #sdf[:M] = cnts.get("Total Campaign")[:M][1]
    #sdf[:Mt] = cnts.get("Total Campaign")[:Mt][1]
    #sdf[:Mc] = cnts.get("Total Campaign")[:Mc][1]

    sdf[:M] = getk(cnts.sdf,"Total Campaign")[:M][1]   #cnts.get("Total Campaign")[:M][1]
    sdf[:Mt] = getk(cnts.sdf,"Total Campaign")[:Mt][1]   #cnts.get("Total Campaign")[:Mt][1]
    sdf[:Mc] = getk(cnts.sdf,"Total Campaign")[:Mc][1]   #cnts.get("Total Campaign")[:Mc][1]

    sdf[:mean_score0]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :dolocc_score0] )
    sdf[:mean_score1]=mean( df_data[ df_data[:buyer_pos_p1] .== 1 ,  :dolocc_score1] )

    sdf[:unadj_mean_score0]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0),:dolocc_score0])
    sdf[:unadj_mean_score1]=mean(df_data[(df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1),:dolocc_score1])

    sdf[:adj_mean_cntrl_grp ] =  mean( df_data[ df_data[:buyer_pos_p1] .== 1 , :dolocc_score0] )
    sdf[:adj_mean_expsd_grp ] =  mean( df_data[ df_data[:buyer_pos_p1] .== 1 ,  :dolocc_score1] )
    # -------  DolOcc Score
    sdf[:unadj_avg_cntrl_hh_pre] =  mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 0), :dol_per_trip_pre_p1] )
    sdf[:unadj_avg_expsd_hh_pre] =  mean(df_data[ (df_data[:buyer_pre_p1] .== 1) & (df_data[:group] .== 1), :dol_per_trip_pre_p1] )
    sdf[:unadj_avg_cntrl_hh_pst] =  mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 0), :dol_per_trip_pos_p1] )
    sdf[:unadj_avg_expsd_hh_pst] =  mean(df_data[ (df_data[:buyer_pos_p1] .== 1) & (df_data[:group] .== 1), :dol_per_trip_pos_p1] )

    #A hack to fix missing data --
    sdf[isnan(sdf[:unadj_avg_expsd_hh_pre]),:unadj_avg_expsd_hh_pre] = 0.0

end
#aggregate!(df_data, mdolocc.feff, cnts)

function aggregate!(df_data::DataFrame, fpen::FPen, cnts::Cnts)
    sdf = fpen.sdf
    #sdf[:M] = cnts.get("Total Campaign")[:N][1]
    #sdf[:Mt] = cnts.get("Total Campaign")[:Nt][1]
    #sdf[:Mc] = cnts.get("Total Campaign")[:Nc][1]

    sdf[:M] = getk(cnts.sdf,"Total Campaign")[:N][1]   #cnts.get("Total Campaign")[:M][1]
    sdf[:Mt] = getk(cnts.sdf,"Total Campaign")[:Nt][1]   #cnts.get("Total Campaign")[:Mt][1]
    sdf[:Mc] = getk(cnts.sdf,"Total Campaign")[:Nc][1]   #cnts.get("Total Campaign")[:Mc][1]

    sdf[:mean_score0]=mean( df_data[ :pen_score0] )
    sdf[:mean_score1]=mean( df_data[ :pen_score1] )
    sdf[:unadj_mean_score0]=mean(df_data[(df_data[:group] .== 0),:pen_score0])
    sdf[:unadj_mean_score1]=mean(df_data[(df_data[:group] .== 1),:pen_score1])

    sdf[:adj_mean_cntrl_grp ] = mean( df_data[ :pen_score0] )
    sdf[:adj_mean_expsd_grp ] = mean( df_data[ :pen_score1] )
    # -------  Pen Score   - buyer_pre_p1
    sdf[:unadj_avg_cntrl_hh_pre ] = mean(df_data[ (df_data[:group] .== 0)  , :buyer_pre_p1] )
    sdf[:unadj_avg_expsd_hh_pre ] = mean(df_data[ (df_data[:group] .== 1)  , :buyer_pre_p1] )
    sdf[:unadj_avg_cntrl_hh_pst ] = mean(df_data[ (df_data[:group] .== 0)  , :buyer_pos_p1] )
    sdf[:unadj_avg_expsd_hh_pst ] = mean(df_data[ (df_data[:group] .== 1)  , :buyer_pos_p1] )

    #A hack to fix missing data --
    sdf[isnan(sdf[:unadj_avg_expsd_hh_pre]),:unadj_avg_expsd_hh_pre] = 0.0

end
#aggregate!(df_data, mpen.feff, cnts)

function aggregateCommon!(fixedeff::FixedEffect)
    sdf = fixedeff.sdf
    #src = fixedeff.src
    sdf[:adj_dod_effct] = ((sdf[:adj_mean_expsd_grp] - sdf[:adj_mean_cntrl_grp])  ./ sdf[:adj_mean_cntrl_grp] ) * 100
    sdf[:twotail_pval] = sdf[:P]   #twotail_pval
    sdf[:onetail_pval] = 1 - (sdf[:twotail_pval] ./ 2)   #onetail_pval
    sdf[:twotail_pval] = 1-sdf[:P]    #redo :twotail_pval
    #For NON-SIG Models, update - adj_mean_expsd_grp (to = adj_mean_cntrl_grp) & adj_dod_effct (= 0)
    #sdf[(sdf[:twotail_pval] .< 0.80) & (sdf[:twotail_pval] .!= 0),:adj_mean_expsd_grp]      =  sdf[(sdf[:twotail_pval] .< 0.80) & (sdf[:twotail_pval] .!= 0),:adj_mean_cntrl_grp]
end
#aggregateCommon!(mocc.feff)
#aggregateCommon!(mdolocc.feff)
#aggregateCommon!(mpen.feff)

# -------- Confidence Intervals - Total --------
function ConfidenceIntervals(feff::FixedEffect)
    sdf = feff.sdf
    for row in eachrow(sdf)
            k = row[:key]
            md = df2dict(row)
            md[:mean_score0] = md[:unadj_mean_score0];
            md[:mean_score1] = md[:unadj_mean_score1];
            cis =  CIs(md, ZDict, 0)
            for (zscore_key,zscore) in  ZDict
                    ubk = symbol(zscore_key*"_ub")
                    lbk = symbol(zscore_key*"_lb")
                    sdf[sdf[:key].==k, lbk] = md[lbk]*100
                    sdf[sdf[:key].==k, ubk] = md[ubk]*100
            end
    end
end

aggregate!(df_data, mocc.feff,cnts)
aggregate!(df_data, mpen.feff, cnts)
aggregate!(df_data, mdolocc.feff, cnts)
aggregateCommon!(mocc.feff)
aggregateCommon!(mdolocc.feff)
aggregateCommon!(mpen.feff)
ConfidenceIntervals(mocc.feff)
ConfidenceIntervals(mdolocc.feff)
ConfidenceIntervals(mpen.feff)

# --------------------------------------------------------
# ---------- BREAKS --------------------------------------
# --------------------------------------------------------

# === Random Effects ====
if !cfg[:TotalModelsOnly]

function init!(cnts::Cnts, m::MModel)
        src=m.reff.src
        # --- GROUP1 -----
        src[:Z1] = src[:B1] ./ src[:SE1]
        src[:P1] = 2.0 * ccdf(Normal(), abs(src[:Z1]))
        src[:SIG1] = map( x -> ((isnan(x))|(x > 0.2)) ?  false : true, src[:P1] )
        # --- GROUP0 -----
        src[:Z0] = src[:B0] ./ src[:SE0]
        src[:P0] = 2.0 * ccdf(Normal(), abs(src[:Z0]))
        src[:SIG0] = map( x -> (isnan(x)|(x > 0.2)) ?  false : true, src[:P0] )

        m.reff.sdf = join(m.reff.sdf,m.reff.src[[:class,:level,:B0,:B1,:SE0,:SE1,:exposed,:P1,:key]], on = :key)
        sdf=m.reff.sdf
        sdf[:B_fixed] = m.feff.sdf[1,:B]
        sdf[:P_fixed] = m.feff.sdf[1,:P]
        sdf[:SE_fixed] = m.feff.sdf[1,:SE]
        #sdf[:B1_combo] = sdf[:B1] + sdf[:B_fixed] #- sdf[:B0]
        #sdf[:SE1_combo] = sqrt(sdf[:SE_fixed].^2+sdf[:SE1].^2)
        sdf[:B1_combo] = sdf[:B1] + sdf[:B_fixed] - sdf[:B0]
        sdf[:SE1_combo] = sqrt(sdf[:SE_fixed].^2+sdf[:SE1].^2+sdf[:SE0].^2)
        sdf[:P1_combo] = 2.0 * ccdf(Normal(), abs( (sdf[:B1_combo]) ./ sdf[:SE1_combo]))

        if m.reff.v_model == "pen"
            for k in sdf[:key]
                sdf[sdf[:key].==k ,:M] = getk(cnts.sdf,k)[:N][1]
                sdf[sdf[:key].==k ,:Mt] = getk(cnts.sdf,k)[:Nt][1]
                sdf[sdf[:key].==k ,:Mc] = getk(cnts.sdf,k)[:Nc][1]
            end
        else
            for k in sdf[:key]
                sdf[sdf[:key].==k ,:M] = getk(cnts.sdf,k)[:M][1]
                sdf[sdf[:key].==k ,:Mt] = getk(cnts.sdf,k)[:Mt][1]
                sdf[sdf[:key].==k ,:Mc] = getk(cnts.sdf,k)[:Mc][1]
            end
        end
end
init!(cnts, mocc)
init!(cnts, mdolocc)
init!(cnts, mpen)

function aggregate!(df_data::DataFrame, mocc::MOcc)
    sdf=mocc.reff.sdf
    for row = eachrow(mocc.reff.sdf)
        k = row[:key]
        ranfx = lowercase(row[:class])
        sranfx = symbol(ranfx)
        v_level = row[:level]
        exposed = row[:exposed]
        dest_colname0=symbol(k*"_occ_0")
        dest_colname1=symbol(k*"_occ_1")
        mean_score0 = mean(df_data[df_data[:buyer_pos_p1] .== 1, dest_colname0])
        mean_score1 = mean(df_data[ df_data[:buyer_pos_p1] .== 1 , dest_colname1])
        #println(mean_score1)
        sdf[sdf[:key].==k ,:mean_score0] = mean_score0
        sdf[sdf[:key].==k ,:mean_score1] = mean_score1
        sdf[sdf[:key].==k ,:adj_mean_cntrl_grp] = mean_score0
        sdf[sdf[:key].==k ,:adj_mean_expsd_grp] = mean_score1

        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pre_p1] )  #(df_data[:buyer_pre_p1] .== 1)
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pos_p1] )  #(df_data[:buyer_pos_p1] .== 1)
        if exposed
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mocc.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mocc.feff.sdf[1,:unadj_avg_cntrl_hh_pst]
        else
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pre_p1] )  #(df_data[:buyer_pre_p1] .== 1)
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :trps_pos_p1] )  #(df_data[:buyer_pos_p1] .== 1)
        end
    end
end
aggregate!(df_data, mocc)

function aggregate!(df_data::DataFrame, mdolocc::MDolOcc)
    sdf=mdolocc.reff.sdf
    for row = eachrow(mdolocc.reff.sdf)
        k = row[:key]
        ranfx = lowercase(row[:class])
        sranfx = symbol(ranfx)
        v_level = row[:level]
        exposed = row[:exposed]
        dest_colname0=symbol(k*"_dolocc_0")
        dest_colname1=symbol(k*"_dolocc_1")
        mean_score0 = mean(df_data[df_data[:buyer_pos_p1] .== 1, dest_colname0])
        mean_score1 = mean(df_data[ df_data[:buyer_pos_p1] .== 1 , dest_colname1])
        sdf[sdf[:key].==k ,:mean_score0] = mean_score0
        sdf[sdf[:key].==k ,:mean_score1] = mean_score1
        sdf[sdf[:key].==k ,:adj_mean_cntrl_grp] = mean_score0
        sdf[sdf[:key].==k ,:adj_mean_expsd_grp] = mean_score1
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )
        #sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
        #sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )

        if exposed
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mdolocc.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mdolocc.feff.sdf[1,:unadj_avg_cntrl_hh_pst]
        else
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pre_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[:buyer_pos_p1] .== 1) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )
            #sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[sranfx] .== v_level) , :dol_per_trip_pre_p1] )
            #sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[sranfx] .== v_level) , :dol_per_trip_pos_p1] )
        end
    end
end
aggregate!(df_data, mdolocc)

function aggregate!(df_data::DataFrame,mpen::MPen)
    sdf=mpen.reff.sdf
    for row = eachrow(mpen.reff.sdf)
        k = row[:key]
        ranfx = lowercase(row[:class])
        sranfx = symbol(ranfx)
        v_level = row[:level]
        exposed = row[:exposed]
        dest_colname0=symbol(k*"_pen_0")
        dest_colname1=symbol(k*"_pen_1")
        mean_score0 = mean(df_data[dest_colname0])
        mean_score1 = mean(df_data[dest_colname1])   #println("..... ",mean_score0," ~ ",mean_score1)
        sdf[sdf[:key].==k ,:mean_score0] = mean_score0
        sdf[sdf[:key].==k ,:mean_score1] = mean_score1
        sdf[sdf[:key].==k ,:adj_mean_cntrl_grp] = mean_score0
        sdf[sdf[:key].==k ,:adj_mean_expsd_grp] = mean_score1
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[symbol(ranfx)] .== v_level) , :buyer_pre_p1] )
        sdf[sdf[:key].==k ,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[symbol(ranfx)] .== v_level) , :buyer_pos_p1] )
        if exposed
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mpen.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mpen.feff.sdf[1,:unadj_avg_cntrl_hh_pst]
        else
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[symbol(ranfx)] .== v_level) , :buyer_pre_p1] )
            sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[symbol(ranfx)] .== v_level) , :buyer_pos_p1] )
        end
    end
end
aggregate!(df_data, mpen)

    #sdf[:B1_combo] = sdf[:B1] + sdf[:B_fixed] - sdf[:B0]
    #sdf[:SE1_combo] = sqrt(sdf[:SE_fixed].^2+sdf[:SE1].^2+sdf[:SE0].^2)
    #sdf[:P1_combo] = 2.0 * ccdf(Normal(), abs( (sdf[:B1_combo]) ./ sdf[:SE1_combo]))

function aggregateCommon!(raneff::RanEffect,cfg::OrderedDict)
    sdf = raneff.sdf
    sdf[:adj_dod_effct] = ((sdf[:adj_mean_expsd_grp] - sdf[:adj_mean_cntrl_grp])  ./ sdf[:adj_mean_cntrl_grp] ) * 100
    sdf[:twotail_pval] = (1-sdf[:P1_combo]) * 100
    sdf[:onetail_pval] = (1-(sdf[:P1_combo] ./ 2)) * 100

    sdf[:Braw] = sdf[:B1] - sdf[:B0]
    sdf[:SEraw] = sqrt(sdf[:SE1].^2+sdf[:SE0].^2)
    sdf[:Zraw] = sdf[:Braw] ./ sdf[:SEraw]
    sdf[:Praw] = 2.0 * ccdf(Normal(), abs(sdf[:Zraw]))
    sdf[:twotail_pval_raw] = (1-sdf[:Praw]) * 100
    sdf[:onetail_pval_raw] = (1-(sdf[:Praw] ./ 2)) * 100

    sigLvl=parse(Float64,get(cfg, :sigLevel, "0.2"))

end
aggregateCommon!(mocc.reff,cfg)
aggregateCommon!(mdolocc.reff,cfg)
aggregateCommon!(mpen.reff,cfg)

#--add 3 random coefs..... combin STD_ERR
# for non-significant - set to 0... or default to total....

function np_val(iDict::OrderedDict)
    M = get(iDict, :M, NA)
    Mt = get(iDict, :Mt, NA)
    Mc = get(iDict, :Mc, NA)
    mean_score0 = get(iDict, :adj_mean_cntrl_grp, NA)
    mean_score1 = get(iDict, :adj_mean_expsd_grp, NA)
    B = get(iDict, :B1_combo, NA)
    SE = get(iDict, :SE1_combo, NA)
    #zscore = get(iDict, :zscore, NA)
    #println("NLsolve : ",iDict,"\n")
    println(get(iDict,:key, NA),"~",M,"~",Mt,"~",Mc,"~",mean_score0,"~",mean_score1,"~",B,"~",SE)

    function f!(x, fvec)
        #fvec[1] = exp(x[1])-45
        Lb_pre = (   (mean_score1*(Mt/M))    -   (mean_score1*exp(-(B-(SE*x[1])))*(Mt/M))   )  +
                 (   (mean_score0*exp((B-(SE*x[1])))*(Mc/M))    -   (mean_score0*(Mc/M))    )
        Lb = Lb_pre/mean_score0
        fvec[1] = Lb
        ### ------------ Lower Bound ---------------
        #Lb_pre = (   (mean_score1*(Mt/M))    -   (mean_score1*exp(-(B-(SE*zscore)))*(Mt/M))   )  +
        #         (   (mean_score0*exp((B-(SE*zscore)))*(Mc/M))    -   (mean_score0*(Mc/M))    )
        #Lb = Lb_pre/mean_score0
        ### ------------ Upper Bound ---------------
        #Ub_pre =  (     ( mean_score1*(Mt/M) )   -   ( mean_score1*exp(-(B+(SE*zscore)))*(Mt/M))   )  +
        #          (     ( mean_score0*exp((B+(SE*zscore)))*(Mc/M))  - (mean_score0*(Mc/M) )   )
        #Ub = Ub_pre/mean_score0
        #md[symbol(zscore_key*"_lb")] = Lb
        #md[symbol(zscore_key*"_ub")] = Ub
    end
    r=nlsolve(f!,[0.1])
    #println(r)
    zvalue=r.zero[1]
    pvalue=2.0 * ccdf(Normal(), abs(zvalue))
    two_tail = 1-pvalue
    one_tail = 1-(pvalue/2)
    return one_tail, two_tail
end

function applyWeights(mocc::MOcc,mdolocc::MDolOcc,mpen::MPen)
    for m in [mocc,mdolocc,mpen]
        println("Processing : ",m.modelName)
        m.reff.sdf = join(m.reff.sdf,cnts.sdf[[:key,:weight]], on = :key)
        sdf=m.reff.sdf
        #sdf[:factor]=0.0  #unadj_mean_score0
        fdf=deepcopy(m.reff.sdf)
        for c in [:unadj_avg_expsd_hh_pre,:unadj_avg_cntrl_hh_pre,:unadj_avg_expsd_hh_pst,:unadj_avg_cntrl_hh_pst,:adj_mean_expsd_grp,:adj_mean_cntrl_grp,:mean_score0,:mean_score1]
            fdf[c] = map((x,y) -> (x*y), fdf[c],fdf[:weight])
        end
        for col in [:unadj_avg_expsd_hh_pre,:unadj_avg_cntrl_hh_pre,:unadj_avg_expsd_hh_pst,:unadj_avg_cntrl_hh_pst,:adj_mean_expsd_grp,:adj_mean_cntrl_grp,:mean_score0,:mean_score1]
            tot=m.feff.sdf[1,col]
            scol=string(col)
            sdf[symbol("f_unweighted_"*scol)] = sdf[col]
            sdf[symbol("f_"*scol)] = 0.0
            lst=by(fdf[!isnan(fdf[col]),:], :class, df -> sum(df[col]))
            rename!(lst, :x1,:creative_weight)
            lst[:factor] = tot ./ lst[:creative_weight]
            for row in eachrow(lst)
                f=row[:factor]
                sdf[sdf[:class].==row[:class], symbol("f_"*scol) ] = f
                sdf[sdf[:class].==row[:class], col ] = sdf[sdf[:class].==row[:class], col ] * f
            end
        end
        sdf[:f_unweighted_onetail_pval] = sdf[:onetail_pval]
        sdf[:f_unweighted_twotail_pval] = sdf[:twotail_pval]
        sdf[:f_unweighted_adj_dod_effct] = sdf[:adj_dod_effct]

        sdf[:adj_dod_effct] = ((sdf[:adj_mean_expsd_grp] - sdf[:adj_mean_cntrl_grp])  ./ sdf[:adj_mean_cntrl_grp] ) * 100
        for i in 1:length(sdf[1])
            d=df2dict(sdf[i,:])
            p1, p2 = np_val(d)
            p1=p1*100
            p2=p2*100
            sdf[i,:onetail_pval] = p1
            sdf[i,:twotail_pval] = p2
        end

       # #RE-Default exposed Break Ctrl cols to total
        ##exposed
        #for col in [:unadj_avg_cntrl_hh_pre,:unadj_avg_cntrl_hh_pst,:adj_mean_cntrl_grp]
        #    totval=m.feff.sdf[1,col]
        #    sdf[sdf[:exposed].==true,col] = totval
        #end
        ##    sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pre] = mpen.feff.sdf[1,:unadj_avg_cntrl_hh_pre]
        ##    sdf[sdf[:key].==k ,:unadj_avg_cntrl_hh_pst] = mpen.feff.sdf[1,:unadj_avg_cntrl_hh_pst]

    end
end
if cfg[:counts]
    applyWeights(mocc,mdolocc,mpen)
end

# -------- Confidence Intervals - Random Effects --------
function ConfidenceIntervals(cnts::Cnts,reff::RanEffect)
    sdf = reff.sdf
    for row in eachrow(sdf)
        k = row[:key]
        md = df2dict(row)    #[:B1_combo,:SE1_combo,:mean_score0,:mean_score1,:M,:Mt,:Mc]
        md[:B] =  md[:B1_combo] # don't really need these...need to clean up
        md[:SE] =  md[:SE1_combo] # ditto

#        md[:M_F] = cnts.sdf[1,:M]
#        md[:Mt_F] = cnts.sdf[1,:Mt]
#        md[:Mc_F] = cnts.sdf[1,:Mc]
        #if row[:exposed]
        #    cis =  CIs(md, ZDict, 1)
        #else
        #    cis =  CIs(md, ZDict, 2)
        #end
        cis =  CIs(md, ZDict)
        for (zscore_key,zscore) in  ZDict
                ubk = symbol(zscore_key*"_ub")
                lbk=symbol(zscore_key*"_lb")
                sdf[sdf[:key].==k, lbk] = md[lbk]*100
                sdf[sdf[:key].==k, ubk] = md[ubk]*100
        end
    end
end
ConfidenceIntervals(cnts, mocc.reff)
ConfidenceIntervals(cnts, mdolocc.reff)
ConfidenceIntervals(cnts, mpen.reff)

#mdolhh.sdf=mdolhh.sdf[mdolhh.sdf[:key].!="creative_groups (none)",:]

end # :TotalModelsOnly

# ==========================================================================================
# =================================== DolHH ================================================
# ==========================================================================================

mdf=MDF(cnts,mocc,mdolocc,mpen)
rdf=RDF(mocc,mdolocc,mpen,mdolhh)

function aggregate!(df_data::DataFrame, mdolhh::MDolHH, smdf::DataFrame)   #mdf::MDF
    sdf=mdolhh.sdf
    xsdf = join(mdolhh.sdf, cnts.sdf[[:key,:class,:level,:exposed]], on = :key)
    for row = eachrow(xsdf)
        k = row[:key]
        #if k != "Total Campaign"
            ranfx = row[:class]
            v_level = row[:level]
            exposed = row[:exposed]
            println("getting : ",k)
            #md = df2dict(mdf.get(k))
            md = df2dict(getk(smdf,k))
            adjctl = md[:o_mean_score0] * md[:y_mean_score0] * md[:p_mean_score0]
            adjexp = md[:o_mean_score1] * md[:y_mean_score1] * md[:p_mean_score1]
            println(k," ~~ ",adjctl,"~",adjexp) #,md)
            sdf[sdf[:key].==k,:adj_mean_cntrl_grp] = adjctl
            sdf[sdf[:key].==k,:adj_mean_expsd_grp] = adjexp
            sdf[sdf[:key].==k,:adj_dod_effct] = ((adjexp - adjctl) ./ adjctl ) * 100
            if k == "Total Campaign"
                sdf[sdf[:key].==k,:unadj_avg_cntrl_hh_pre] =  mean(df_data[ (df_data[:group] .== 0)  , :prd_1_net_pr_pre] )
                sdf[sdf[:key].==k,:unadj_avg_expsd_hh_pre] =  mean(df_data[ (df_data[:group] .== 1)  , :prd_1_net_pr_pre] )
                sdf[sdf[:key].==k,:unadj_avg_cntrl_hh_pst] =  mean(df_data[ (df_data[:group] .== 0)  , :prd_1_net_pr_pos] )
                sdf[sdf[:key].==k,:unadj_avg_expsd_hh_pst] =  mean(df_data[ (df_data[:group] .== 1)  , :prd_1_net_pr_pos] )
            else
                sdf[sdf[:key].== k,:unadj_avg_expsd_hh_pre] = mean(df_data[ (df_data[:group] .== 1) & (df_data[symbol(lowercase(ranfx))] .== v_level) , :prd_1_net_pr_pre] )
                sdf[sdf[:key].== k,:unadj_avg_expsd_hh_pst] = mean(df_data[ (df_data[:group] .== 1) & (df_data[symbol(lowercase(ranfx))] .== v_level) , :prd_1_net_pr_pos] )
                if exposed
                    sdf[sdf[:key].== k,:unadj_avg_cntrl_hh_pre] = sdf[1,:unadj_avg_cntrl_hh_pre]
                    sdf[sdf[:key].== k,:unadj_avg_cntrl_hh_pst] = sdf[1,:unadj_avg_cntrl_hh_pst]
                else
                    sdf[sdf[:key].== k,:unadj_avg_cntrl_hh_pre] = mean(df_data[ (df_data[:group] .== 0) & (df_data[symbol(lowercase(ranfx))] .== v_level) , :prd_1_net_pr_pre] )
                    sdf[sdf[:key].== k,:unadj_avg_cntrl_hh_pst] = mean(df_data[ (df_data[:group] .== 0) & (df_data[symbol(lowercase(ranfx))] .== v_level) , :prd_1_net_pr_pos] )
                end
            end
        #end
    end
end
aggregate!(df_data,mdolhh,genMDF(cnts,mocc,mdolocc,mpen))
#aggregate!(df_data,mdolhh,genMDF(cnts, mocc, mdolocc, mpen))

function calc_rawPvals_Opt(iDict::OrderedDict)
    dout = iDict
    M = get(iDict, :M, NA)
    Mt = get(iDict, :Mt, NA)
    Mc = get(iDict, :Mc, NA)
    N = get(iDict, :N, NA)
    Nt = get(iDict, :Mt, NA)
    Nc = get(iDict, :Nc, NA)
    B1 = get(iDict, :B1, NA)
    B2 = get(iDict, :B2, NA)
    B3 = get(iDict, :B3, NA)
    SE1 = get(iDict, :SE1, NA)
    SE2 = get(iDict, :SE2, NA)
    SE3 = get(iDict, :SE3, NA)
    SEsq=sqrt(SE1^2+SE2^2+SE3^2)
    o_mean_score0 = get(iDict, :o_mean_score0, NA)
    o_mean_score1 = get(iDict, :o_mean_score1, NA)
    y_mean_score0 = get(iDict, :y_mean_score0, NA)
    y_mean_score1 = get(iDict, :y_mean_score1, NA)
    p_mean_score0 = get(iDict, :p_mean_score0, NA)
    p_mean_score1 =get(iDict, :p_mean_score1, NA)
    Bsum=B1+B2+B3
    dout[:Bsum] = Bsum
    ###### PVALUE - ONE & TWO ########
    m=nothing
    m = Model(solver=NLoptSolver(algorithm=:LD_MMA, maxtime=v_ttl))
    if B1 < 0
        @variable(m, Bocc >= B1)
    else
        @variable(m, Bocc <= B1)
    end
    if B2 < 0
        @variable(m, Bdolocc >= B2)
    else
        @variable(m, Bdolocc <= B2)
    end
    if B3 < 0
        @variable(m, Bpen >= B3)
    else
        @variable(m, Bpen <= B3)
    end
    @objective(m, Max, (((Bocc+Bpen+Bdolocc)-Bsum)/SEsq ))   # z = (x-u)/se
    @NLconstraint(m, 0.00000 <= ((((p_mean_score1*(Nt/N))+(p_mean_score0*exp(Bpen)*(Nc/N)))
	                           * ((o_mean_score1*(Mt/M))+(o_mean_score0*exp(Bocc)*(Mc/M)))
	                           * ((y_mean_score1*(Mt/M))+(y_mean_score0*exp(Bdolocc)*(Mc/M)))
	                             )
	                           -(((p_mean_score1*(Nt/N)*exp(-Bpen))+(p_mean_score0*(Nc/N)))
	                           *((o_mean_score1*(Mt/M)*exp(-Bocc))+(o_mean_score0*(Mc/M)))
	                           *((y_mean_score1*(Mt/M)*exp(-Bdolocc))+(y_mean_score0*(Mc/M)))
	                             )
	                          )
	               <= 0.00001
                    )
    #print(m)
    status = solve(m)
    zvalue=getobjectivevalue(m)
    pvalue=2.0 * ccdf(Normal(), abs(zvalue))
    two_tail = 1-pvalue
    one_tail = 1-(pvalue/2)
    dout[:onetail_pval] = one_tail
    dout[:twotail_pval] = two_tail
    println("z-value: ", string(zvalue)," --> p-value: ",string(two_tail))
    return dout
end

function rawPvals(mdolhh::MDolHH,x::DataFrame)
    sdf=mdolhh.sdf
    sdf[:twotail_pval_raw] = 0.0
    sdf[:onetail_pval_raw]= 0.0
    for row = eachrow(x)
        k = row[:key]
        if k != "Total Campaign"
            md = df2dict(x[x[:key].==k,:])
            md[:metakey] = k
            println("Raw PValue : ",k)
            if :B1_orig in keys(md)  #!!!!!!!!!
                md[:B1] = md[:B1_orig]
                md[:SE1] = md[:SE1_orig]
                md[:B2] = md[:B2_orig]
                md[:SE2] = md[:SE2_orig]
                md[:B3] = md[:B3_orig]
                md[:SE3] = md[:SE3_orig]
            println(md)
                calc_rawPvals_Opt(md)
                sdf[sdf[:key].==k, :onetail_pval_raw]=md[:onetail_pval]*100
                sdf[sdf[:key].==k, :twotail_pval_raw]=md[:twotail_pval]*100
            end
        end
    end
end
#rawPvals(mdolhh,genRndMDF(cnts, mocc, mdolocc, mpen, false))
rawPvals(mdolhh,genMDF(cnts, mocc, mdolocc, mpen))
#mdolhh.sdf[[:key,:twotail_pval,:onetail_pval,:twotail_pval_raw,:onetail_pval_raw]]

function ConfidenceIntervals(mdolhh::MDolHH,x::DataFrame)
    sdf=mdolhh.sdf
    for row = eachrow(x)
        k = row[:key]
        md = df2dict(x[x[:key].==k,:])
        md[:metakey] = k
        println("getting : ",k)
        calcPValue_Opt(md)
        CIs_O(md)
        for zk in [ :onetail_pval, :twotail_pval, :onetail_80_pct_intrvl_lb,:onetail_80_pct_intrvl_ub, :onetail_90_pct_intrvl_lb, :onetail_90_pct_intrvl_ub
                    ,:twotail_80_pct_intrvl_lb, :twotail_80_pct_intrvl_ub, :twotail_90_pct_intrvl_lb, :twotail_90_pct_intrvl_ub
                  ]
            sdf[ sdf[:key].==k, zk]= md[zk]*100
        end
    end
end
#ConfidenceIntervals(mdolhh,genRndMDF(cnts, mocc, mdolocc, mpen, false))
#ConfidenceIntervals(mdolhh,genFixedMDF(cnts, mocc, mdolocc, mpen))
ConfidenceIntervals(mdolhh,genMDF(cnts, mocc, mdolocc, mpen))

#function vcatreff(mocc, mdolocc, mpen)
#    #dfo=DataFrame()
#    for m in [mocc,mdolocc, mpen]
#        if m.hasBreaks
#
#    end
#end
#    function rmap(m::MModel)
#        return m.hasBreaks
#    end
#        rvs=vcat(map(x->if x.hasBreaks rmap(x) end,[mocc,mdolocc, mpen]))

function extendRDF(rdf::DataFrame)
    dfo=DataFrame(MODEL_DESC=ASCIIString[], Model=ASCIIString[], TIME_AGG_PERIOD=Int64[], START_WEEK=Int64[], END_WEEK=Int64[], dependent_variable=ASCIIString[],
                  CNT_EXPSD_HH=Int64[], UNADJ_AVG_EXPSD_HH_PRE=Float64[], UNADJ_AVG_CNTRL_HH_PRE=Float64[], UNADJ_AVG_EXPSD_HH_PST=Float64[],
                  UNADJ_AVG_CNTRL_HH_PST=Float64[], UNADJ_DOD_EFFCT=Float64[], UNADJ_DIFF_EFFCT=Float64[], ADJ_MEAN_EXPSD_GRP=Float64[],
                  ADJ_MEAN_CNTRL_GRP=Float64[], ADJ_DOD_EFFCT=Float64[], TWOTAIL_PVAL=Float64[], ONETAIL_PVAL=Float64[], ABS_DIFF=Float64[],
                  DOL_DIFF=Float64[], ONETAIL_80_PCT_INTRVL_UB=Float64[], ONETAIL_80_PCT_INTRVL_LB=Float64[], ONETAIL_90_PCT_INTRVL_UB=Float64[],
                  ONETAIL_90_PCT_INTRVL_LB=Float64[], TWOTAIL_80_PCT_INTRVL_UB=Float64[], TWOTAIL_80_PCT_INTRVL_LB=Float64[], TWOTAIL_90_PCT_INTRVL_UB=Float64[],
                  TWOTAIL_90_PCT_INTRVL_LB=Float64[], CNT_IMPRESSIONS=Int64[], TWOTAIL_PVAL_to_Campaign=Float64[], ONETAIL_PVAL_to_Campaign=Float64[],
                  CNT_Model_HH=Float64[]
                 )

  ###NBNBNB : this needs to come out - only to hack total only
  #          rvs = vcat(mocc.reff.sdf[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]],
#               mdolocc.reff.sdf[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]],
#               mpen.reff.sdf[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]],
##               mdolhh.sdf[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]]
#              )
    if (!cfg[:TotalModelsOnly]) & (length(filter(x-> x.hasBreaks ,[mocc,mdolocc, mpen]))  > 0 )
            rvs=vcat( map(x-> x.reff.sdf, filter(x-> x.hasBreaks ,[mocc,mdolocc, mpen])))[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]]
            rvs = vcat(rvs,mdolhh.sdf[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]])
    else
        rvs = mdolhh.sdf[[:key,:model,:twotail_pval_raw,:onetail_pval_raw]]
    end

    for i in 1:length(rdf[1])
        push!(dfo, [
                     rdf[i,:key],
                     NA, #Model=NA, # code
                     NA, #TIME_AGG_PERIOD=NA,
                     NA, #START_WEEK=NA,
                     NA, #END_WEEK=NA,
                     rdf[i,:model],
                     NA, #CNT_EXPSD_HH=Int64[],
                     rdf[i,:unadj_avg_expsd_hh_pre],
                     rdf[i,:unadj_avg_cntrl_hh_pre],
                     rdf[i,:unadj_avg_expsd_hh_pst],
                     rdf[i,:unadj_avg_cntrl_hh_pst],
                     NA, #UNADJ_DOD_EFFCT=Float64[],
                     NA, #UNADJ_DIFF_EFFCT=Float64[],
                     rdf[i,:adj_mean_expsd_grp],
                     rdf[i,:adj_mean_cntrl_grp],
                     rdf[i,:adj_dod_effct],
                     rdf[i,:twotail_pval],
                     rdf[i,:onetail_pval],
                     NA, #ABS_DIFF=Float64[],
                     NA, #DOL_DIFF=Float64[],
                     rdf[i,:onetail_80_pct_intrvl_ub],
                     rdf[i,:onetail_80_pct_intrvl_lb],
                     rdf[i,:onetail_90_pct_intrvl_ub],
                     rdf[i,:onetail_90_pct_intrvl_lb],
                     rdf[i,:twotail_80_pct_intrvl_ub],
                     rdf[i,:twotail_80_pct_intrvl_lb],
                     rdf[i,:twotail_90_pct_intrvl_ub],
                     rdf[i,:twotail_90_pct_intrvl_lb],
                     NA, #CNT_MPRESSIONS=Float64[],
                     NA, #TWOTAIL_PVAL_to_Campaign=Float64[],
                     NA, #ONETAIL_PVAL_to_Campaign=Float64[],
                     NA #CNT_Model_HH=Float64[]
                   ]
             )
        x=length(dfo[1])
        dfo[x,:UNADJ_DOD_EFFCT] = ( ((dfo[x,:UNADJ_AVG_EXPSD_HH_PST] - dfo[x,:UNADJ_AVG_EXPSD_HH_PRE]) - (dfo[x,:UNADJ_AVG_CNTRL_HH_PST] - dfo[x,:UNADJ_AVG_CNTRL_HH_PRE]))  /  dfo[x,:UNADJ_AVG_CNTRL_HH_PST] ) *100

        dfo[x,:UNADJ_DIFF_EFFCT] = ((dfo[x,:UNADJ_AVG_EXPSD_HH_PST] - dfo[x,:UNADJ_AVG_CNTRL_HH_PST]) / dfo[x,:UNADJ_AVG_CNTRL_HH_PST] )* 100

        if dfo[i,:dependent_variable] == "dolhh"
            k=dfo[x,:MODEL_DESC]
            if cfg[:counts]
                dfo[x,:CNT_EXPSD_HH] =  getk(cnts.sdf,k,:hh)   #  cnts.get(k,:hh)
                dfo[x,:CNT_IMPRESSIONS] =  getk(cnts.sdf,k,:impressions)  #cnts.get(k,:impressions)
            end
            dfo[x,:DOL_DIFF] = dfo[x,:ADJ_MEAN_EXPSD_GRP] - dfo[x,:ADJ_MEAN_CNTRL_GRP]
            dfo[x,:CNT_Model_HH] = cnts.sdf[cnts.sdf[:key].==dfo[x,:MODEL_DESC],:M][1]

        elseif dfo[i,:dependent_variable] == "pen"
            dfo[x,:CNT_Model_HH] = cnts.sdf[cnts.sdf[:key].==dfo[x,:MODEL_DESC],:Nt][1]
        else
            dfo[x,:CNT_Model_HH] = cnts.sdf[cnts.sdf[:key].==dfo[x,:MODEL_DESC],:Mt][1]
        end

        if (dfo[x,:dependent_variable] in ["occ","dolocc","pen","dolhh"]) & (dfo[x,:MODEL_DESC] !== "Total Campaign")

            pv = rvs[(rvs[:key].== dfo[x,:MODEL_DESC])&(rvs[:model].== dfo[x,:dependent_variable]) ,:]
            if length(pv[1]) >=1
                dfo[i,:TWOTAIL_PVAL_to_Campaign] = pv[:twotail_pval_raw][1]
                dfo[i,:ONETAIL_PVAL_to_Campaign] = pv[:onetail_pval_raw][1]
            end
        end

        #format for unify
        dfo[dfo[:dependent_variable].=="occ", :dependent_variable] = "OCC"
        dfo[dfo[:dependent_variable].=="dolocc", :dependent_variable] = "DOL/OCC"
        dfo[dfo[:dependent_variable].=="pen", :dependent_variable] = "PEN"
        dfo[dfo[:dependent_variable].=="dolhh", :dependent_variable] = "DOL/HH"

    end
    dfo[:ABS_DIFF] = dfo[:ADJ_MEAN_EXPSD_GRP] - dfo[:ADJ_MEAN_CNTRL_GRP]

    return dfo
end
xrdf=extendRDF(genRDF(mocc,mdolocc,mpen, mdolhh,cfg) )

#xrdf[[:MODEL_DESC,:dependent_variable,:ADJ_MEAN_EXPSD_GRP,:ADJ_MEAN_CNTRL_GRP]]
#x=genMDF(cnts,mocc,mdolocc,mpen)
#x[[:key,:o_mean_score0,:y_mean_score0,:p_mean_score0,:o_mean_score1,:y_mean_score1,:p_mean_score1]]
#write2disk(xrdf,"/home/rmadmin/g/StatStack/src/xrdf.csv")
write2disk(xrdf,pwd()*"/scored.csv")

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
lines = readlines(ARGS[1])

function solution1(lines)
    hit_map = Dict()
    for line in lines
        coords = []
        step = 1
        for sub_block in split(line, " -> ")
            append!(coords, [parse(Int, x) for x in split(sub_block, ",")])
        end

        x_diff = coords[1] - coords[3]
        y_diff = coords[2] - coords[4]

        if x_diff == 0
            if coords[2] > coords[4]
                step = -1
            end
            for i = coords[2]:step:coords[4]
                entry = string(coords[1], ",", i)
                if entry in keys(hit_map)
                    hit_map[entry] += 1
                else
                    hit_map[entry] = 1
                end
            end
        elseif y_diff == 0
            if coords[1] > coords[3]
                step = -1
            end
            for i = coords[1]:step:coords[3]
                entry = string(i, ",", coords[2])
                if entry in keys(hit_map)
                    hit_map[entry] += 1
                else
                    hit_map[entry] = 1
                end
            end
        else
            println("skip")
        end
    end
    return length([k for (k, v) in hit_map if v >= 2])
end

function solution2(lines)
    hit_map = Dict()
    for line in lines
        coords = []
        step = 1
        for sub_block in split(line, " -> ")
            append!(coords, [parse(Int, x) for x in split(sub_block, ",")])
        end

        x_diff = coords[1] - coords[3]
        y_diff = coords[2] - coords[4]

        if x_diff == 0
            if coords[2] > coords[4]
                step = -1
            end
            for i = coords[2]:step:coords[4]
                entry = string(coords[1], ",", i)
                if entry in keys(hit_map)
                    hit_map[entry] += 1
                else
                    hit_map[entry] = 1
                end
            end
        elseif y_diff == 0
            if coords[1] > coords[3]
                step = -1
            end
            for i = coords[1]:step:coords[3]
                entry = string(i, ",", coords[2])
                if entry in keys(hit_map)
                    hit_map[entry] += 1
                else
                    hit_map[entry] = 1
                end
            end
        elseif abs(x_diff) == abs(y_diff)
            y_step = 1
            if x_diff > 0
                step = -1
            end
            if y_diff > 0
                y_step = -1
            end
            for (x_pos, y_pos) in zip(coords[1]:step:coords[3], coords[2]:y_step:coords[4])
                entry = string(x_pos, ",", y_pos)
                if entry in keys(hit_map)
                    hit_map[entry] += 1
                else
                    hit_map[entry] = 1
                end
            end
        else
            println("skip")
        end
    end
    return length([k for (k, v) in hit_map if v >= 2])
end

println(solution1(lines))
println(solution2(lines))

ext: jl


text:
# ---
# title: 289. Game of Life
# id: problem289
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/game-of-life/description/>
# hidden: true
# ---
#
# According to the [Wikipedia's
# article](https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life): "The **Game
# of Life** , also known simply as **Life** , is a cellular automaton devised by
# the British mathematician John Horton Conway in 1970."
#
# Given a _board_ with _m_ by _n_ cells, each cell has an initial state _live_
# (1) or _dead_ (0). Each cell interacts with its [eight
# neighbors](https://en.wikipedia.org/wiki/Moore_neighborhood) (horizontal,
# vertical, diagonal) using the following four rules (taken from the above
# Wikipedia article):
#
#   1. Any live cell with fewer than two live neighbors dies, as if caused by under-population.
#   2. Any live cell with two or three live neighbors lives on to the next generation.
#   3. Any live cell with more than three live neighbors dies, as if by over-population..
#   4. Any dead cell with exactly three live neighbors becomes a live cell, as if by reproduction.
#
# Write a function to compute the next state (after one update) of the board
# given its current state. The next state is created by applying the above rules
# simultaneously to every cell in the current state, where births and deaths
# occur simultaneously.
#
# **Example:**
#
#
#
#     Input: [
#       [0,1,0],
#       [0,0,1],
#       [1,1,1],
#       [0,0,0]
#     ]
#     Output: [
#       [0,0,0],
#       [1,0,1],
#       [0,1,1],
#       [0,1,0]
#     ]
#
#
# **Follow up** :
#
#   1. Could you solve it in-place? Remember that the board needs to be updated at the same time: You cannot update some cells first and then use their updated values to update other cells.
#   2. In this question, we represent the board using a 2D array. In principle, the board is infinite, which would cause problems when the active area encroaches the border of the array. How would you address these problems?
#
#
## @lc code=start
using LeetCode

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


text:
# ---
# title: 1579. Remove Max Number of Edges to Keep Graph Fully Traversable
# id: problem1579
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Union Find
# link: <https://leetcode.com/problems/remove-max-number-of-edges-to-keep-graph-fully-traversable/description/>
# hidden: true
# ---
#
# Alice and Bob have an undirected graph of `n` nodes and 3 types of edges:
#
#   * Type 1: Can be traversed by Alice only.
#   * Type 2: Can be traversed by Bob only.
#   * Type 3: Can by traversed by both Alice and Bob.
#
# Given an array `edges` where `edges[i] = [typei, ui, vi]` represents a
# bidirectional edge of type `typei` between nodes `ui` and `vi`, find the
# maximum number of edges you can remove so that after removing the edges, the
# graph can still be fully traversed by both Alice and Bob. The graph is fully
# traversed by Alice and Bob if starting from any node, they can reach all other
# nodes.
#
# Return _the maximum number of edges you can remove, or return_ `-1` _if it 's
# impossible for the graph to be fully traversed by Alice and Bob._
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/19/ex1.png)**
#
#
#
#     Input: n = 4, edges = [[3,1,2],[3,2,3],[1,1,3],[1,2,4],[1,1,2],[2,3,4]]
#     Output: 2
#     Explanation: If we remove the 2 edges [1,1,2] and [1,1,3]. The graph will still be fully traversable by Alice and Bob. Removing any additional edge will not make it so. So the maximum number of edges we can remove is 2.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/19/ex2.png)**
#
#
#
#     Input: n = 4, edges = [[3,1,2],[3,2,3],[1,1,4],[2,1,4]]
#     Output: 0
#     Explanation: Notice that removing any edge will not make the graph fully traversable by Alice and Bob.
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/19/ex3.png)**
#
#
#
#     Input: n = 4, edges = [[3,2,3],[1,1,2],[2,3,4]]
#     Output: -1
#     Explanation: In the current graph, Alice cannot reach node 4 from the other nodes. Likewise, Bob cannot reach 1. Therefore it's impossible to make the graph fully traversable.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#   * `1 <= edges.length <= min(10^5, 3 * n * (n-1) / 2)`
#   * `edges[i].length == 3`
#   * `1 <= edges[i][0] <= 3`
#   * `1 <= edges[i][1] < edges[i][2] <= n`
#   * All tuples `(typei, ui, vi)` are distinct.
#
#
## @lc code=start
using LeetCode

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


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

"""
Frictionless 2d small sliding contact without forwarddiff.

true/false flags: finite_sliding, friction, use_forwarddiff
"""
function assemble!(problem::Problem{Contact}, time::Float64,
                   ::Type{Val{1}}, ::Type{Val{false}},
                   ::Type{Val{false}}, ::Type{Val{false}}; debug=false)

    props = problem.properties
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)

    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
    normals, tangents = calculate_normals(slave_elements, time, Val{1};
                        rotate_normals=props.rotate_normals)
    update!(slave_elements, "normal", time => normals)
    update!(slave_elements, "tangent", time => tangents)

    Rn = 0.0

    # 2. loop all slave elements
    for slave_element in slave_elements

        nsl = length(slave_element)
        X1 = slave_element("geometry", time)
        u1 = slave_element("displacement", time)
        la1 = slave_element("reaction force", time)
        n1 = slave_element("normal", time)
        t1 = slave_element("tangent", time)
        x1 = X1 + u1
        Q1_ = [n1[1] t1[1]]
        Q2_ = [n1[2] t1[2]]
        Z = zeros(2, 2)
        Q2 = [Q1_ Z; Z Q2_]
        contact_area = 0.0
        contact_error = 0.0

        if "element area" in props.store_fields
            element_area = 0.0
            for ip in get_integration_points(slave_element)
                detJ = slave_element(ip, time, Val{:detJ})
                w = ip.weight*detJ
                element_area += w
            end
            update!(slave_element, "element area", time => element_area)
        end

        # 3. loop all master elements
        for master_element in slave_element("master elements", time)

            nm = length(master_element)
            X2 = master_element("geometry", time)
            u2 = master_element("displacement", time)
            x2 = X2 + u2

            if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval
                continue
            end

            if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval
                continue
            end

            # 3.1 calculate segmentation
            xi1a = project_from_master_to_slave(slave_element, X2[1], time)
            xi1b = project_from_master_to_slave(slave_element, X2[2], time)
            xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
            l = 1/2*abs(xi1[2]-xi1[1])
            isapprox(l, 0.0) && continue # no contribution in this master element

            # 3.2. bi-orthogonal basis
            De = zeros(nsl, nsl)
            Me = zeros(nsl, nsl)
            Ae = zeros(nsl, nsl)
            if props.dual_basis
                for ip in get_integration_points(slave_element, 3)
                    detJ = slave_element(ip, time, Val{:detJ})
                    w = ip.weight*detJ*l
                    xi = ip.coords[1]
                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    De += w*diagm(N1)
                    Me += w*N1*N1'
                end
                Ae = De*inv(Me)
            else
                Ae = eye(nsl)
            end

            # 3.3. loop integration points of one integration segment and calculate
            # local mortar matrices
            fill!(De, 0.0)
            fill!(Me, 0.0)
            Ne = zeros(nsl, 2*nsl)
            Te = zeros(nsl, 2*nsl)
            He = zeros(nsl, 2*nsl)
            ce = zeros(nsl)
            ge = zeros(nsl)
            for ip in get_integration_points(slave_element, 3)
                detJ = slave_element(ip, time, Val{:detJ})
                w = ip.weight*detJ*l
                xi = ip.coords[1]
                xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                N1 = vec(get_basis(slave_element, xi_s, time))
                Phi = Ae*N1

                # project gauss point from slave element to master element in direction n_s
                X_s = N1*X1 # coordinate in gauss point
                n_s = N1*n1 # normal direction in gauss point
                t_s = N1*t1 # tangent condition in gauss point
                n_s /= norm(n_s)
                t_s /= norm(t_s)
                xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
                N2 = vec(get_basis(master_element, xi_m, time))
                X_m = N2*X2

                u_s = N1*u1
                u_m = N2*u2
                x_s = X_s + u_s
                x_m = X_m + u_m
                la_s = Phi*la1

                # virtual work
                De += w*Phi*N1'
                Me += w*Phi*N2'

                # contact constraints
                Ne += w*reshape(kron(N1, n_s, Phi), 2, 4)
                Te += w*reshape(kron(N2, n_s, Phi), 2, 4)
                He += w*reshape(kron(N1, t_s, Phi), 2, 4)
                ge += w*Phi*dot(n_s, x_m-x_s)
                ce += w*N1*dot(n_s, -la_s)
                Rn += w*dot(n_s, -la_s)

                contact_area += w
                contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
            end

            sdofs = get_gdofs(problem, slave_element)
            mdofs = get_gdofs(problem, master_element)

            # add contribution to contact virtual work
            for i=1:field_dim
                lsdofs = sdofs[i:field_dim:end]
                lmdofs = mdofs[i:field_dim:end]
                add!(problem.assembly.C1, lsdofs, lsdofs, De)
                add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
            end

            # add contribution to contact constraints
            add!(problem.assembly.C2, sdofs[1:field_dim:end], sdofs, Ne)
            add!(problem.assembly.C2, sdofs[1:field_dim:end], mdofs, -Te)
            add!(problem.assembly.D, sdofs[2:field_dim:end], sdofs, He)
            add!(problem.assembly.g, sdofs[1:field_dim:end], ge)
            add!(problem.assembly.c, sdofs[1:field_dim:end], ce)

        end # master elements done

        if "contact area" in props.store_fields
            update!(slave_element, "contact area", time => contact_area)
        end

        if "contact error" in props.store_fields
            update!(slave_element, "contact error", time => contact_error)
        end

    end # slave elements done, contact virtual work ready

    S = sort(collect(keys(normals))) # slave element nodes
    weighted_gap = Dict{Int64, Vector{Float64}}()
    contact_pressure = Dict{Int64, Vector{Float64}}()
    complementarity_condition = Dict{Int64, Vector{Float64}}()
    is_active = Dict{Int64, Int}()
    is_inactive = Dict{Int64, Int}()
    is_slip = Dict{Int64, Int}()
    is_stick = Dict{Int64, Int}()

    la = problem.assembly.la
    ndofs = length(la)
    info("contact ndofs: $ndofs")
    info("Rn = $Rn")

    C1 = sparse(problem.assembly.C1)
    C2 = sparse(problem.assembly.C2, ndofs, ndofs)
    D = sparse(problem.assembly.D, ndofs, ndofs)
    g = full(problem.assembly.g, ndofs, 1)
    c = full(problem.assembly.c, ndofs, 1)

    # active / inactive node detection
    for j in S
        dofs = [2*(j-1)+1, 2*(j-1)+2]
        weighted_gap[j] = g[dofs]

        if length(la) != 0
            p = dot(normals[j], la[dofs])
            t = dot(tangents[j], la[dofs])
            contact_pressure[j] = [p, t]
        else
            contact_pressure[j] = [0.0, 0.0]
        end

#       contact_pressure[j] = c[dofs]
        complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
        if complementarity_condition[j][1] < 0
            is_inactive[j] = 1
            is_active[j] = 0
            is_slip[j] = 0
            is_stick[j] = 0
        else
            is_inactive[j] = 0
            is_active[j] = 1
            is_slip[j] = 1
            is_stick[j] = 0
            _c1 = complementarity_condition[j][1]
            _c2 = c[dofs]
            _c3 = contact_pressure[j][1]
            _c4 = g[dofs]
            info("active $j: c1 = $_c1, c2 = $_c2, c3 = $_c3, c4 = $_c4")
        end
    end

    if "weighted gap" in props.store_fields
        update!(slave_elements, "weighted gap", time => weighted_gap)
    end
    if "contact pressure" in props.store_fields
        update!(slave_elements, "contact pressure", time => contact_pressure)
    end
    if "complementarity condition" in props.store_fields
        update!(slave_elements, "complementarity condition", time => complementarity_condition)
    end
    if "active nodes" in props.store_fields
        update!(slave_elements, "active nodes", time => is_active)
    end
    if "inactive nodes" in props.store_fields
        update!(slave_elements, "inactive nodes", time => is_inactive)
    end
    if "stick nodes" in props.store_fields
        update!(slave_elements, "stick nodes", time => is_stick)
    end
    if "slip nodes" in props.store_fields
        update!(slave_elements, "slip nodes", time => is_slip)
    end

    info("# | active | inactive | stick | slip | gap | pres | comp")
    for j in S
        str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | "
        str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))"
        info(str1 * str2)
    end

    # solve variational inequality

    # constitutive modelling in tangent direction, frictionless contact
    #=
    for j in S
        dofs = [2*(j-1)+1, 2*(j-1)+2]
        if (is_active[j] == 1) && (is_slip[j] == 1)
            info("$j is in active/slip, removing tangential constraint $(dofs[2])")
            C2[dofs[2],:] = 0.0
            g[dofs[2]] = 0.0
            D[dofs[2], dofs] = tangents[j]
        end
    end
    =#

    # remove inactive nodes from assembly
    for j in S
        dofs = [2*(j-1)+1, 2*(j-1)+2]
        if is_inactive[j] == 1
            info("$j is inactive, removing dofs $dofs")
            C1[dofs,:] = 0.0
            C2[dofs,:] = 0.0
            D[dofs,:] = 0.0
            g[dofs,:] = 0.0
        end
    end

    problem.assembly.C1 = C1
    problem.assembly.C2 = C2
    problem.assembly.D = D
    problem.assembly.g = g

end


text:
module SimString

import Base: push!, append!
using DataStructures: DefaultOrderedDict, DefaultDict
using CircularArrays
using OffsetArrays

######### Import modules & utils ################
include("db_collection.jl")
include("dictdb.jl")
include("features.jl")
include("measures.jl")
include("search.jl")

####### Global export of user API #######
export Dice, Jaccard, Cosine, Overlap, ExactMatch,
    AbstractSimStringDB, DictDB, describe_collection,
    CharacterNGrams, WordNGrams,
    search

end


text:
using SeniorProjectAssignments
using Random

rng = MersenneTwister(3)
roles = Set([:hardware, :software, :electronics])

n_students = 100
n_projects = 10

projects = [ProjectData(Symbol(string("project_", j)), (hardware = 2, software = 2, electronics = 2), 9, 12) for j in 1:n_projects]
pnames = [p.id for p in projects]

students = StudentData[]
for i in 1:n_students
    specialty = rand(rng, roles)
    secondary = rand(rng, delete!(copy(roles), specialty))
    sd = StudentData(string(i), (;specialty=>0.7, secondary=>0.3), pnames[randperm(rng, length(pnames))])
    push!(students, sd)
end

group3(i) = [string(j) for j in (i-1)*3+1:i*3]
groups = [group3(i) for i in 1:floor(Int, n_students/3)]

matching = SeniorProjectAssignments.match(students, projects, groups)


text:
################################################################
# This method  computes the set of descriptors through the nonlinear scale space
function Compute_Descriptors(akaze, kpts)

    t1 = time_ns()

    ## Allocate memory for the matrix with the descriptors
    ## We use the full length binary descriptor -> 486 bits
    t = if (akaze.options_.descriptor_size != 0)
        akaze.options_.descriptor_size
    else
        ## We use the random bit selection length binary descriptor
        (6+36+120)*akaze.options_.descriptor_channels
    end
    desc_len = ceil(Int, t/8.0)
    desc = zeros(UInt8, desc_len, length(kpts))

    if akaze.options_.descriptor == MLDB
        Compute_Main_Orientation.([akaze], kpts)
    end

    descriptor_function = if akaze.options_.descriptor == MLDB_UPRIGHT
        if akaze.options_.descriptor_size == 0
            Get_Upright_MLDB_Full_Descriptor
        else
            Get_Upright_MLDB_Descriptor_Subset
        end
    elseif akaze.options_.descriptor == MLDB
        if akaze.options_.descriptor_size == 0
            Get_MLDB_Full_Descriptor
        else
            Get_MLDB_Descriptor_Subset
        end
    end

    for (n,kp) in enumerate(kpts)
        descriptor_function(akaze, kp, @view desc[:,n])
    end

    t2 = time_ns()
    akaze.timing_.descriptor = t2 - t1
    desc
end

mymin(a,b) = ifelse(a<b,a,b)
mymax(a,b) = ifelse(a>b,a,b)

################################################################
function Compute_Main_Orientation(akaze, kpt)

    resX = Vector{Float64}(undef, 109)
    resY = Vector{Float64}(undef, 109)
    Ang = Vector{Float64}(undef, 109)

    ## Get the information from the keypoint
    level = kpt.class_id
    ratio = Float64(1<<akaze.evolution_[level].octave)
    s = round(Int64, 0.5*kpt.size/ratio)
    xf = kpt.pt.x/ratio
    yf = kpt.pt.y/ratio

    idx = 1
    ## Calculate derivatives responses for points within radius of 6*scale
    (limj,limk) = size(akaze.evolution_[level].Lx)
    for i in -6:6
        for j in -6:6
            if i*i + j*j < 36
                iy = mymin(limj, mymax(1, round(Int64, yf + j*s)+1))
                ix = mymin(limk, mymax(1, round(Int64, xf + i*s)+1))

                gweight = gauss25[abs(i)+1, abs(j)+1]
                resX[idx] = gweight * akaze.evolution_[level].Lx[iy, ix]
                resY[idx] = gweight * akaze.evolution_[level].Ly[iy, ix]
                Ang[idx] = mod(atan(resY[idx], resX[idx]), 2π)
                idx += 1
            end
        end
    end

    ## Variables for computing the dominant direction
    maxXY = 0.0
    ## Loop slides pi/3 window around feature point
    for ang1 in 0:0.15:2.0*π
        ang2 = if (ang1 + π/3.0 > 2.0*π) ang1 - 5.0 * π/3.0 else ang1 + π/3.0 end
        sumX = 0.0
        sumY = 0.0

        for k in 1:109
            ## Get angle from the x-axis of the sample point
            ang = Ang[k]

            ## Determine whether the point is within the window
            if (ang1 < ang2 && ang1 < ang && ang < ang2)
                sumX += resX[k]
                sumY += resY[k]
            elseif (ang2 < ang1 && ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2.0*π)))
                sumX += resX[k]
                sumY += resY[k]
            end
        end

        ## if the vector produced from this window is longer than all
        ## previous vectors then this forms the new dominant direction
        if (sumX*sumX + sumY*sumY > maxXY)
            ## store largest orientation
            maxXY = sumX*sumX + sumY*sumY
            kpt.angle = mod(atan(sumY, sumX), 2π)
        end
    end
end

################################################################
function Get_MLDB_Full_Descriptor(akaze, kpt, desc)

    max_channels = 3
    @assert akaze.options_.descriptor_channels <= max_channels
    values = zeros(16*max_channels)
    size_mult = [1.0, 2.0/3.0, 1.0/2.0]

    ratio = 1 << kpt.octave
    scale = round(Int, 0.5 * kpt.size / ratio)
    xf = kpt.pt.x / ratio
    yf = kpt.pt.y / ratio
    co = cos(kpt.angle)
    si = sin(kpt.angle)
    pattern_size = akaze.options_.descriptor_pattern_size

    dpos = Ref(0)
    for lvl in 0:2
        val_count = (lvl + 2) * (lvl + 2)
        sample_step = ceil(Int, pattern_size * size_mult[1+lvl])
        MLDB_Fill_Values(akaze, values, sample_step, kpt.class_id, xf, yf, co, si, scale)
        MLDB_Binary_Comparisons(akaze, values, desc, val_count, dpos)
    end
end

# /* ************************************************************************* */
# void AKAZE::Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const

# const int max_channels = 3
# CV_Assert(options_.descriptor_channels <= max_channels)
# float values[16*max_channels]
# const double size_mult[3] = 1, 2.0/3.0, 1.0/2.0end

# float ratio = (float)(1 << kpt.octave)
# float scale = (float)fRound(0.5f*kpt.size / ratio)
# float xf = kpt.pt.x / ratio
# float yf = kpt.pt.y / ratio
# int pattern_size = options_.descriptor_pattern_size

# int dpos = 0
# for(int lvl = 0; lvl < 3; lvl++)
#     int val_count = (lvl + 2) * (lvl + 2)
#     int sample_step = static_cast<int>(ceil(pattern_size * size_mult[lvl]))
#     MLDB_Fill_Upright_Values(values, sample_step, kpt.class_id, xf, yf, scale)
#     MLDB_Binary_Comparisons(values, desc, val_count, dpos)
# end
# end

################################################################
function MLDB_Fill_Values(akaze, values, sample_step::Int, level::Int,
                          xf::Real, yf::Real, co::Real, si::Real, scale::Real)

    pattern_size = akaze.options_.descriptor_pattern_size
    nr_channels = akaze.options_.descriptor_channels
    valpos = 0

    for i in -pattern_size:sample_step:pattern_size-1
        for j in -pattern_size:sample_step:pattern_size-1

            di = 0.0; dx = 0.0; dy = 0.0
            nsamples = 0

            for k in i:i + sample_step - 1
                for l in j:j + sample_step - 1

                    sample_y = yf + (l*co*scale + k*si*scale)
                    sample_x = xf + (-l*si*scale + k*co*scale)

                    y1 = round(Int, sample_y)
                    x1 = round(Int, sample_x)

                    ri = akaze.evolution_[level].Lt[1+y1,1+x1]
                    di += ri

                    if nr_channels > 1
                        rx = akaze.evolution_[level].Lx[1+y1, 1+x1]
                        ry = akaze.evolution_[level].Ly[1+y1, 1+x1]
                        if (nr_channels == 2)
                            dx += sqrt(rx*rx + ry*ry)
                        else
                            rry = rx*co + ry*si
                            rrx = -rx*si + ry*co
                            dx += rrx
                            dy += rry
                        end
                    end
                    nsamples+=1
                end
            end

            di /= nsamples
            dx /= nsamples
            dy /= nsamples

            values[1+valpos] = di

            if (nr_channels > 1)
                values[1+valpos + 1] = dx
            end

            if (nr_channels > 2)
                values[1+valpos + 2] = dy
            end

            valpos += nr_channels
        end
    end
end

# /* ************************************************************************* */
# void AKAZE::MLDB_Fill_Upright_Values(float* values, int sample_step, int level,
#                                      float xf, float yf, float scale) const

# int pattern_size = options_.descriptor_pattern_size
# int nr_channels = options_.descriptor_channels
# int valpos = 0

# for (int i = -pattern_size; i < pattern_size; i += sample_step)
#     for (int j = -pattern_size; j < pattern_size; j += sample_step)

#         float di = 0.0, dx = 0.0, dy = 0.0
#         int nsamples = 0

#         for (int k = i; k < i + sample_step; k++)
#             for (int l = j; l < j + sample_step; l++)

#                 float sample_y = yf + l*scale
#                 float sample_x = xf + k*scale

#                 int y1 = fRound(sample_y)
#                 int x1 = fRound(sample_x)

#                 float ri = *(evolution_[level].Lt.ptr<float>(y1)+x1)
#                 di += ri

#                 if(nr_channels > 1)
#                     float rx = *(evolution_[level].Lx.ptr<float>(y1)+x1)
#                     float ry = *(evolution_[level].Ly.ptr<float>(y1)+x1)
#                     if (nr_channels == 2)
#                         dx += sqrtf(rx*rx + ry*ry)
#                     end
#                 else
#                     dx += rx
#                     dy += ry
#                 end
#             end
#             nsamples++
#         end
#     end

#     di /= nsamples
#     dx /= nsamples
#     dy /= nsamples

#     values[valpos] = di

#     if (nr_channels > 1)
#         values[valpos + 1] = dx

#         if (nr_channels > 2)
#             values[valpos + 2] = dy

#             valpos += nr_channels
#         end
#     end
# end

################################################################
function MLDB_Binary_Comparisons(akaze, values, desc,
                                 count::Int, dpos::Ref{Int})

    nr_channels = akaze.options_.descriptor_channels

    for pos in 0:nr_channels-1
        for i in 0:count-1
            ival = values[1 + nr_channels * i + pos]
            for j in i + 1:count-1
                res = ival > values[1 + nr_channels * j + pos]
                desc[1+dpos[] >> 3] |= (res << (dpos[] & 7))
                dpos[]+=1
            end
        end
    end
end

# /* ************************************************************************* */
# void AKAZE::Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc)

# float di = 0.f, dx = 0.f, dy = 0.f
# float rx = 0.f, ry = 0.f
# float sample_x = 0.f, sample_y = 0.f
# int x1 = 0, y1 = 0

# ## Get the information from the keypoint
# float ratio = (float)(1<<kpt.octave)
# int scale = fRound(0.5*kpt.size/ratio)
# float angle = kpt.angle
# float level = kpt.class_id
# float yf = kpt.pt.y/ratio
# float xf = kpt.pt.x/ratio
# float co = cos(angle)
# float si = sin(angle)

# ## Allocate memory for the matrix of values
# cv::Mat values = cv::Mat_<float>::zeros((4+9+16)*options_.descriptor_channels, 1)

# ## Sample everything, but only do the comparisons
# vector<int> steps(3)
# steps.at(0) = options_.descriptor_pattern_size
# steps.at(1) = ceil(2.f*options_.descriptor_pattern_size/3.f)
# steps.at(2) = options_.descriptor_pattern_size/2

# for (int i=0; i < descriptorSamples_.rows; i++)
#     int *coords = descriptorSamples_.ptr<int>(i)
#     int sample_step = steps.at(coords[0])
#     di=0.0f
#     dx=0.0f
#     dy=0.0f

#     for (int k = coords[1]; k < coords[1] + sample_step; k++)
#         for (int l = coords[2]; l < coords[2] + sample_step; l++)

#             ## Get the coordinates of the sample point
#             sample_y = yf + (l*scale*co + k*scale*si)
#             sample_x = xf + (-l*scale*si + k*scale*co)

#             y1 = fRound(sample_y)
#             x1 = fRound(sample_x)

#             di += *(evolution_[level].Lt.ptr<float>(y1)+x1)

#             if (options_.descriptor_channels > 1)
#                 rx = *(evolution_[level].Lx.ptr<float>(y1)+x1)
#                 ry = *(evolution_[level].Ly.ptr<float>(y1)+x1)

#                 if (options_.descriptor_channels == 2)
#                     dx += sqrtf(rx*rx + ry*ry)
#                 end
#             else if (options_.descriptor_channels == 3)
#                 ## Get the x and y derivatives on the rotated axis
#                 dx += rx*co + ry*si
#                 dy += -rx*si + ry*co
#             end
#             end
#         end
#     end

#     *(values.ptr<float>(options_.descriptor_channels*i)) = di

#     if (options_.descriptor_channels == 2)
#         *(values.ptr<float>(options_.descriptor_channels*i+1)) = dx
#     else if (options_.descriptor_channels == 3)
#         *(values.ptr<float>(options_.descriptor_channels*i+1)) = dx
#         *(values.ptr<float>(options_.descriptor_channels*i+2)) = dy
#     end
#     end

#     ## Do the comparisons
#     const float *vals = values.ptr<float>(0)
#     const int *comps = descriptorBits_.ptr<int>(0)

#     for (int i=0; i<descriptorBits_.rows; i++)
#         if (vals[comps[2*i]] > vals[comps[2*i +1]])
#             desc[i/8] |= (1<<(i%8))
#         end
#     end
# end

# /* ************************************************************************* */
# void AKAZE::Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc)

# float di = 0.0f, dx = 0.0f, dy = 0.0f
# float rx = 0.0f, ry = 0.0f
# float sample_x = 0.0f, sample_y = 0.0f
# int x1 = 0, y1 = 0

# ## Get the information from the keypoint
# float ratio = (float)(1<<kpt.octave)
# int scale = fRound(0.5*kpt.size/ratio)
# float level = kpt.class_id
# float yf = kpt.pt.y/ratio
# float xf = kpt.pt.x/ratio

# ## Allocate memory for the matrix of values
# cv::Mat values = cv::Mat_<float>::zeros((4+9+16)*options_.descriptor_channels, 1)

# vector<int> steps(3)
# steps.at(0) = options_.descriptor_pattern_size
# steps.at(1) = ceil(2.f*options_.descriptor_pattern_size/3.f)
# steps.at(2) = options_.descriptor_pattern_size/2

# for (int i=0; i < descriptorSamples_.rows; i++)
#     int *coords = descriptorSamples_.ptr<int>(i)
#     int sample_step = steps.at(coords[0])
#     di=0.0f, dx=0.0f, dy=0.0f

#     for (int k = coords[1]; k < coords[1] + sample_step; k++)
#         for (int l = coords[2]; l < coords[2] + sample_step; l++)

#             ## Get the coordinates of the sample point
#             sample_y = yf + l*scale
#             sample_x = xf + k*scale

#             y1 = fRound(sample_y)
#             x1 = fRound(sample_x)
#             di += *(evolution_[level].Lt.ptr<float>(y1)+x1)

#             if (options_.descriptor_channels > 1)
#                 rx = *(evolution_[level].Lx.ptr<float>(y1)+x1)
#                 ry = *(evolution_[level].Ly.ptr<float>(y1)+x1)

#                 if (options_.descriptor_channels == 2)
#                     dx += sqrtf(rx*rx + ry*ry)
#                 end
#             else if (options_.descriptor_channels == 3)
#                 dx += rx
#                 dy += ry
#             end
#             end
#         end
#     end

#     *(values.ptr<float>(options_.descriptor_channels*i)) = di

#     if (options_.descriptor_channels == 2)
#         *(values.ptr<float>(options_.descriptor_channels*i+1)) = dx
#     else if (options_.descriptor_channels == 3)
#         *(values.ptr<float>(options_.descriptor_channels*i+1)) = dx
#         *(values.ptr<float>(options_.descriptor_channels*i+2)) = dy
#     end
#     end

#     ## Do the comparisons
#     const float *vals = values.ptr<float>(0)
#     const int *comps = descriptorBits_.ptr<int>(0)
#     for (int i=0; i<descriptorBits_.rows; i++)
#         if (vals[comps[2*i]] > vals[comps[2*i +1]])
#             desc[i/8] |= (1<<(i%8))
#         end
#     end
# end


text:
MarkovMap(t::IntervalMap) = MarkovMap(t.branches,t.domain,t.rangedomain)
IntervalMap(t::MarkovMap) = IntervalMap(t.branches,t.domain,t.rangedomain)

MarkovMap(t::FwdCircleMap) = forwardmodulomap(t.f,Interval(t.domain),Interval(t.rangedomain),t.dfdx)
MarkovMap(t::RevCircleMap) = reversemodulomap(t.v,Interval(t.domain),Interval(t.rangedomain),t.dvdx)
IntervalMap(t::AbstractCircleMap) = IntervalMap(MarkovMap(t))


text:
using MPI
using MPIClusterManagers
using Distributed
using DistributedArrays
using Test
using LinearAlgebra

let
    manager = MPIManager(np = 6)
    addprocs(manager)

    @everywhere using COSMA

    COSMA.use_manager(manager)

    @testset "Incompatible sizes" begin
        @test_throws DimensionMismatch mul!(drand(6, 6), drand(3, 6), drand(5, 6))
        @test_throws DimensionMismatch mul!(drand(ComplexF64, 6, 6), drand(ComplexF64, 6, 3)', drand(ComplexF64, 5, 6))
        @test_throws DimensionMismatch mul!(drand(6, 6), transpose(drand(6, 3)), drand(5, 6))
        @test_throws DimensionMismatch mul!(drand(4, 3), transpose(drand(3, 3)), drand(3, 3))
    end

    @testset "C = alpha * $f(A) * $g(B) + beta * C ($T)" for f = (identity, transpose, adjoint), g = (identity, transpose, adjoint), T = (Float32, Float64, ComplexF32, ComplexF64)
        # Skip adjoints of real matrices
        if (f == adjoint || g == adjoint) && T <: Real
            continue
        end

        m, n, k = 20, 10, 30

        alpha = T(2.0)
        beta = T(3.0)

        # Distributed arrays
        C_d = drand(T, m, n)
        A_d = f(f != identity ? drand(T, k, m) : drand(T, m, k))
        B_d = g(g != identity ? drand(T, n, k) : drand(T, k, n))

        # Do the same computation locally
        C_h = collect(C_d)
        A_h = collect(A_d)
        B_h = collect(B_d)

        # Reference
        mul!(C_h, A_h, B_h, alpha, beta)
        mul!(C_d, A_d, B_d, alpha, beta)

        # Test by copying to the first node
        @test all(collect(C_d) .≈ C_h)
    end
end


text:
using Test
using TraceGraph
using TraceGraph:tracegraph
using GraphPlot
using LightGraphs

@testset "Basic Tests" begin
    function foo(a)
            b = a+10
            c = b+10
            return (a,b,c)
    end
    gd = tracegraph(foo, 100)
    add_norecurse(reshape)
    show_norecurse()
    rm_norecurse(reshape)
end
@testset "Conditional Statement" begin
    function foo(a::Bool)
        if a == true
            return 100+100
        else
            return 100-100
        end
    end
    gd = tracegraph(foo, false)
    gd = tracegraph(foo, true)
end


text:
struct Grid{T<:Shape,B<:Boundary}
    nrows::Int
    ncols::Int
    cells::Array{Cell{T},2}
end

function Grid(T1::Type{S}, T2::Type{B}, nrows, ncols) where {S<:Shape,B<:Boundary}
    cells = Array{Cell{T1},2}(undef, nrows, ncols)
    for row in 1:nrows
        for col in 1:ncols
            cells[row, col] = Cell(T1, row, col)
        end
    end
    return Grid{T1,T2}(nrows, ncols, cells)
end

function Base.getindex(grid::Grid, i::Int64, j::Int64)
    return grid.cells[i, j]
end

function Base.setindex!(grid::Grid, X, i, j)
    return grid.cells[i, j] = X
end

function size(grid::Grid)
    return grid.ncols * grid.nrows
end


text:
module SEC

include("Calibrate/Calibrate.jl")
#include("Emulate/Emulate.jl")
#include("Sample/Sample.jl")

include("utils.jl")

end # module

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
# Copyright 2019 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

"""
A lightweight (primitive), integer-like handle to a resource.  Useful for dispatch.
i.e. pretty much an integer, but with a unique dispatch signature.
"""
module Handles

export
Handle,
@HandleType

abstract type Handle <: Integer end

"""
    @HandleType(n)

Create a new type named n to be used as a handle to your resource.

Example:
julia> @HandleType(Foo)
Foo

julia> foo(x::Integer) = x*x
foo (generic function with 1 method)

julia> foo(Foo(2))
4

julia> foo(x::Foo) = x-x
foo (generic function with 2 methods)

julia> foo(Foo(2))
0
"""
macro HandleType(tname)
  esc(
      quote
        primitive type $tname <: Handles.Handle sizeof(Int)*8 end
        $tname(x::T) where {T<:Integer} = reinterpret($tname, Int(x))
      end
     )
end

Base.show(io::IO, x::Handle) = print(io, Int(x))

Base.ndigits0zpb( c::C, i::Integer ) where C<:Handle = ndigits0zpb( Int(c), i )
Base.promote_rule( ::Type{C}, ::Type{I} ) where {C<:Handle,I<:Integer} = C
Base.convert( ::Type{C}, number::Number ) where {C<:Handle} = C(Int(number))

Base.Int(x::Handle) = reinterpret(Int,x)
Base.UInt(x::Handle) = reinterpret(UInt,x)
Base.Int32(x::Handle) = Int32(reinterpret(Int,x))
Base.UInt32(x::Handle) = UInt32(reinterpret(UInt,x))

Base.hash(d::Handle, x::UInt64) = hash(reinterpret(UInt64,d),x)

for op in (:+, :-, :*, :/, :mod, :div, :rem, :max, :min)
  @eval Base.$op(a::H, b::H) where {H<:Handle} = H($op(Int(a),Int(b)))
end

for op in (:<, :>, :<=, :>=)
  @eval Base.$op(a::H, b::H) where {H<:Handle} = $op(Int(a),Int(b))
end

for op in (:<<, :>>)
  @eval Base.$op(a::H, b::I) where {H<:Handle, I<:Int} = H($op(Int(a),Int(b)))
end

Base.sub_with_overflow( a::H, b::H ) where {H<:Handle} =
  ((r,f)->(H(r),f))(Base.sub_with_overflow(Int(a), Int(b))...)
Base.add_with_overflow( a::H, b::H ) where {H<:Handle} =
  ((r,f)->(H(r),f))(Base.add_with_overflow(Int(a), Int(b))...)

end


text:
#-------------------------------------------------------------------------------------------------
#--
#-- Create a range of default parameters
#--
using DataFrames, DataFramesMeta

function forecastError(country::String, actual::Array{Float64}, forecast::Array{Float64},
                                        beds::Array{Float64}, icus::Array{Float64})
    # number of forcast points
    l = length(actual)

    # If logarithms are used, avoid errors.
    actual   = max.(actual,   0.0) .+ 1.0  # Should never happen
    forecast = max.(forecast, 0.0) .+ 1.0

    # Log error
    actual = log.(actual)
    forecast = log.(forecast)

    err_deaths = forecast .- actual

    # Penality for being negative (less forecast than actual)
    # To which extent forecast is under actual
    Δ = (100 .* max.(actual .- forecast, 0.0)).^2
    err_deaths = err_deaths .+ Δ

    # Error on bed use
    BED_Max = countryData[country][:hospital_capacity]
    err_beds = exp.(min.(BED_Max .* ones(l) - beds, 0.0))

    # Error on ICU use
    ICU_Max = countryData[country][:ICU_capacity]
    err_icus = exp.(min.(ICU_Max .* ones(l) - icus, 0.0))

    # Prepare the value to never be negative (to avoid log errors)
    return sqrt( sum((err_deaths .+ err_beds .+ err_icus).^2) / l )
end

function singleCountryLoss(country::String, countryparams)
    # finalDate = nothig to force using only the time span of actual recorded deaths
    sol = calculateSolution(country, DiseaseParameters, countryparams; finalDate = nothing)

    # Extract total deaths profile
    actual = convert(Array, countryData[country][:cases][:, :deaths])
    deaths = forecastCompartmentOnActualDates(sol, "D", country)
    beds = forecastVariableOnActualDates(sol, "BED", country)
    icus = forecastVariableOnActualDates(sol, "ICU", country)

    return forecastError(country, actual, deaths, beds, icus)
end

#--
#-- Calculate the sum of all the losses of all the countries to optimise disease params.
#-- Loss per country is sized as if all countries had the same 1m population
#--
function sumCountryLossesDisease(diseaseparams)
    # The parameters passed to the individual loss is created with a mask defined
    totalError = 0.0

    for (country, _) in COUNTRY_LIST
        countryparams = countryData[country][:params]

        totalError += singleCountryLoss(country, countryparams)
    end

    return totalError
end

function fullEpidemyLoss(params)

    # Deconstruct the entire parameter stack
    # First are the disease parameters
    diseaseparams = params[1:DISEASE_N]

    totalError = 0.0

    # Then each country for which the loss is immediately calculated
    for n in 1:COUNTRY_LIST_N
        country, _ = COUNTRY_LIST[n]

        country_start_index = DISEASE_N + (n - 1) * COUNTRY_N + 1
        country_final_index = DISEASE_N + (n - 1) * COUNTRY_N + COUNTRY_N
        countryparams = params[country_start_index:country_final_index]

        # finalDate = nothig to force using only the time span of actual recorded deaths
        sol = calculateSolution(country, diseaseparams, countryparams; finalDate = nothing)

        # Extract total deaths profile
        # Extract total deaths profile
        actual = convert(Array, countryData[country][:cases][:, :deaths])
        deaths = forecastCompartmentOnActualDates(sol, "D", country)
        beds = forecastVariableOnActualDates(sol, "BED", country)
        icus = forecastVariableOnActualDates(sol, "ICU", country)

        # loss =  sum( (log.(actual) .- log.(forecast)).^ 2 ) / length(actual)
        totalError += forecastError(country, actual, deaths, beds, icus)
    end

    return totalError
end

function updateEpidemiologyOnce(;maxtime = 60)
    # Optimise the epidemiology
    println("OPTIMISING EPIDEMIOLOGY---------------------------")

    result = bboptimize(sumCountryLossesDisease,
                        SearchRange = DISEASE_RANGE;
                        Method = :adaptive_de_rand_1_bin,
                        MaxTime = maxtime,
                        NThreads = Threads.nthreads(),
                        TraceMode = :compact)

    global DiseaseParameters = best_candidate(result)
end

function updateCountryOnce(country; maxtime = 60)
    # Make a note of the disease parameters
    println(country)
    print("Before                     ")
    @show countryData[country][:params]

    countryRange = COUNTRY_RANGE
    countryRange[COUNTRY_PARAM_START] = approximateModelStartRange(country)

    # Determine optimal parameters for each countryw
    result = bboptimize(countryData[country][:lossFunction],
                        SearchRange = countryRange;
                        Method = :adaptive_de_rand_1_bin,
                        MaxTime = maxtime,
                        TraceMode = :compact)

    print("After "); @show best_candidate(result)
    println();
    global countryData[country][:params] = best_candidate(result)
end

function sumCountryLossesCountries(params)
    totalError = 0.0

    # Then each country for which the loss is immediately calculated
    for n in 1:COUNTRY_N
        country, _ = COUNTRY_LIST[n]

        country_start_index = (n - 1) * COUNTRY_N + 1
        country_final_index = (n - 1) * COUNTRY_N + COUNTRY_N
        countryparams = params[country_start_index:country_final_index]

        # loss =  sum( (log.(actual) .- log.(forecast)).^ 2 ) / length(actual)
        totalError += singleCountryLoss(country, countryparams)
    end

    return totalError
end


text:
using DataDeps

register(DataDep(
  "breast-cancer-wisconsin-original",
  "http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+%28Original%29",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data",
  "402c585309c399237740f635ef9919dc512cca12cbeb20de5e563a4593f22b64",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=11,
        feature_indices=2:10,
        id_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
    tiles = Tiler(areawidth, areaheight, nrows, ncols, margin=20)

A Tiler is an iterator that, for each iteration, returns a tuple of:

- the `x`/`y` point of the center of each tile in a set of tiles that divide up a
  rectangular space such as a page into rows and columns (relative to current 0/0)

- the number of the tile

`areawidth` and `areaheight` are the dimensions of the area to be tiled, `nrows`/`ncols`
are the number of rows and columns required, and `margin` is applied to all four
edges of the area before the function calculates the tile sizes required.

Tiler and Partition are similar:

- Partition lets you specify the width and height of a cell

- Tiler lets you specify how many rows and columns of cells you want, and a margin:

```
tiles = Tiler(1000, 800, 4, 5, margin=20)
for (pos, n) in tiles
    # the point pos is the center of the tile
end
```

You can access the calculated tile width and height like this:

```
tiles = Tiler(1000, 800, 4, 5, margin=20)
for (pos, n) in tiles
    ellipse(pos.x, pos.y, tiles.tilewidth, tiles.tileheight, :fill)
end
```

It's sometimes useful to know which row and column you're currently on.
`tiles.currentrow` and `tiles.currentcol` should have that information for you.

To use a Tiler to make grid points:

```
first.(collect(Tiler(800, 800, 4, 4)))
```

which returns an array of points that are the center points of the grid.
"""
mutable struct Tiler
    areawidth::Real
    areaheight::Real
    tilewidth::Real
    tileheight::Real
    nrows::Int
    ncols::Int
    currentrow::Int
    currentcol::Int
    margin::Real
    function Tiler(areawidth, areaheight, nrows::Int, ncols::Int; margin=10)
        tilewidth  = (areawidth - 2margin)/ncols
        tileheight = (areaheight - 2margin)/nrows
        currentrow=1
        currentcol=1
        new(areawidth, areaheight, tilewidth, tileheight, nrows, ncols, currentrow, currentcol, margin)
    end
end

function Base.iterate(pt::Tiler)
    x = -(pt.areawidth/2)  + pt.margin + (pt.tilewidth/2)
    y = -(pt.areaheight/2) + pt.margin + (pt.tileheight/2)
    tilenumber = 1
    x1 = x + pt.tilewidth
    y1 = y
    if x1 > (pt.areawidth/2) - pt.margin
        y1 += pt.tileheight
        x1 = -(pt.areawidth/2) + pt.margin + (pt.tilewidth/2)
    end
    pt.currentrow, pt.currentcol = (div(tilenumber-1, pt.ncols)+1, mod1(tilenumber, pt.ncols))
    return ((Point(x, y), tilenumber), (Point(x1, y1), tilenumber + 1))
end

function Base.iterate(pt::Tiler, state)
    if state[2] > (pt.nrows * pt.ncols)
        return
    end
    x = state[1].x
    y = state[1].y
    tilenumber = state[2]
    x1 = x + pt.tilewidth
    y1 = y
    if x1 > (pt.areawidth/2) - pt.margin
        y1 += pt.tileheight
        x1 = -(pt.areawidth/2) + pt.margin + (pt.tilewidth/2)
    end
    pt.currentrow, pt.currentcol = (div(tilenumber-1, pt.ncols)+1, mod1(tilenumber, pt.ncols))
    return ((Point(x, y), tilenumber), (Point(x1, y1), tilenumber + 1))
end

function Base.length(pt::Tiler)
    pt.nrows * pt.ncols
end

function Base.getindex(pt::Tiler, i::Int)
    1 <= i <= pt.ncols *  pt.nrows || throw(BoundsError(pt, i))
    xcoord = -pt.areawidth/2  + pt.margin + mod1(i, pt.ncols) * pt.tilewidth  - pt.tilewidth/2
    ycoord = -pt.areaheight/2 + pt.margin + (div(i - 1,  pt.ncols) * pt.tileheight) + pt.tileheight/2
    return (Point(xcoord, ycoord), i)
end

function Base.size(pt::Tiler)
    return (pt.nrows, pt.ncols)
end

Base.lastindex(pt::Tiler) = length(pt)
Base.firstindex(pt::Tiler) = 1

Base.eltype(::Type{Tiler}) = Tuple

Base.getindex(pt::Tiler, I) = [pt[i] for i in I]

"""
    GridRect(startpoint, xspacing, yspacing, width, height)

Define a rectangular grid, to start at `startpoint` and proceed along the x-axis in
steps of `xspacing`, then along the y-axis in steps of `yspacing`.

    GridRect(startpoint, xspacing=100.0, yspacing=100.0, width=1200.0, height=1200.0)

For a column, set the `xspacing` to 0:

    grid = GridRect(O, 0, 40)

To get points from the grid, use `nextgridpoint(g::Grid)`.

    julia> grid = GridRect(O, 0, 40);
    julia> nextgridpoint(grid)
    Luxor.Point(0.0, 0.0)

    julia> nextgridpoint(grid)
    Luxor.Point(0.0, 40.0)

When you run out of grid points, you'll wrap round and start again.
"""
mutable struct GridRect
    startpoint::Point
    currentpoint::Point
    xspacing::Float64
    yspacing::Float64
    width::Float64
    height::Float64
    rownumber::Int
    colnumber::Int
    function GridRect(startpoint, xspacing, yspacing, width=1200.0, height=1200.0)
        rownumber = 1
        colnumber = 0
        # find the "previous" point, so that the first call gets the real first point
        xspacing == 0 ? sp = startpoint.y - yspacing : sp = startpoint.y
        currentpoint = Point(startpoint.x - xspacing, sp)
        new(startpoint, currentpoint, xspacing, yspacing, width, height, rownumber, colnumber)
    end
end

"""
    GridHex(startpoint, radius, width=1200.0, height=1200.0)

Define a hexagonal grid, to start at `startpoint` and proceed along the x-axis and
then along the y-axis, `radius` is the radius of a circle that encloses each hexagon.
The distance in `x` between the centers of successive hexagons is:

``\\frac{\\sqrt{(3)} radius}{2}``

To get the next point from the grid, use `nextgridpoint(g::Grid)`.

When you run out of grid points, you'll wrap round and start again.
"""
mutable struct GridHex
    startpoint::Point
    radius::Float64
    currentpoint::Point
    width::Float64
    height::Float64
    rownumber::Int
    colnumber::Int
    function GridHex(startpoint, radius, width=1200.0, height=1200.0)
        rownumber = 1
        colnumber = 0
        # find the "previous" point, so that the first call gets the real first point
        currentpoint = Point(startpoint.x - (sqrt(3) * radius), startpoint.y)
        new(startpoint, radius, currentpoint, width, height, rownumber, colnumber)
    end
end

"""
    nextgridpoint(g::GridRect)

Returns the next available (or even the first) grid point of a grid.
"""
function nextgridpoint(g::GridRect)
    temp = Point(g.currentpoint.x + g.xspacing, g.currentpoint.y)
    if g.xspacing == 0
        g.rownumber += 1
        g.colnumber = 1
        g.currentpoint = Point(g.startpoint.x, temp.y + g.yspacing)
    else
        g.colnumber += 1
        g.currentpoint = temp
    end
    if abs(g.startpoint.x - temp.x) >= g.width # next row?
        g.rownumber += 1
        g.colnumber = 1
        g.currentpoint = Point(g.startpoint.x, temp.y + g.yspacing)
    end
    if g.currentpoint.y >= g.height   # finished, start again?
        g.rownumber = 1
        g.colnumber = 1
        g.currentpoint = Point(g.startpoint.x, g.startpoint.y)
    end
    return g.currentpoint
end

"""
    nextgridpoint(g::GridHex)

Returns the next available grid point of a hexagonal grid.
"""
function nextgridpoint(g::GridHex)
    g.currentpoint = Point(g.currentpoint.x + (g.radius * sqrt(3)), g.currentpoint.y)
    g.colnumber += 1
    if g.currentpoint.x >= g.width                 # next row
        g.rownumber += 1
        g.colnumber = 1
        if g.rownumber % 2 == 0
            g.currentpoint = Point(g.startpoint.x + (sqrt(3) * g.radius)/2, g.currentpoint.y + (3/2) * g.radius)
        else
            g.currentpoint = Point(g.startpoint.x, g.currentpoint.y + (3/2) * g.radius)
        end
    end
    if g.currentpoint.y >= g.height     # finished?
        g.rownumber = 1                 # start again...
        g.colnumber = 1
        g.currentpoint = Point(g.startpoint.x, g.startpoint.y)
    end
    return g.currentpoint
end

"""
    p = Partition(areawidth, areaheight, tilewidth, tileheight)

A Partition is an iterator that, for each iteration, returns a tuple of:

- the `x`/`y` point of the center of each tile in a set of tiles that divide up a
rectangular space such as a page into rows and columns (relative to current 0/0)

- the number of the tile

`areawidth` and `areaheight` are the dimensions of the area to be tiled,
`tilewidth`/`tileheight` are the dimensions of the tiles.

Tiler and Partition are similar:

- Partition lets you specify the width and height of a cell

- Tiler lets you specify how many rows and columns of cells you want, and a margin

    tiles = Partition(1200, 1200, 30, 30)
    for (pos, n) in tiles
        # the point pos is the center of the tile
    end

You can access the calculated tile width and height like this:

    tiles = Partition(1200, 1200, 30, 30)
    for (pos, n) in tiles
        ellipse(pos.x, pos.y, tiles.tilewidth, tiles.tileheight, :fill)
    end

It's sometimes useful to know which row and column you're currently on:

    tiles.currentrow
    tiles.currentcol

should have that information for you.

Unless the tilewidth and tileheight are exact multiples of the area width and height, you'll
see a border at the right and bottom where the tiles won't fit.
"""
mutable struct Partition
    areawidth::Real
    areaheight::Real
    tilewidth::Real
    tileheight::Real
    nrows::Int
    ncols::Int
    currentrow::Int
    currentcol::Int
    function Partition(areawidth, areaheight, tilewidth, tileheight)
        ncols = areawidth  / tilewidth  |> floor |> Integer
        nrows = areaheight / tileheight |> floor |> Integer
        if ncols < 1 || nrows < 1
            throw(error("partition(): not enough space for tiles that size"))
        end
        currentrow = 1
        currentcol = 1
        new(areawidth, areaheight, tilewidth, tileheight, nrows, ncols, currentrow, currentcol)
    end
end

function Base.iterate(pt::Partition)
    x = -(pt.areawidth/2)  + (pt.tilewidth/2)
    y = -(pt.areaheight/2) + (pt.tileheight/2)
    tilenumber = 1
    x1 = x + pt.tilewidth
    y1 = y
    if (x1 + pt.tilewidth/2) > (pt.areawidth/2)
        y1 += pt.tileheight
        x1 = -(pt.areawidth/2) + (pt.tilewidth/2)
    end
    pt.currentrow, pt.currentcol = (div(tilenumber-1, pt.ncols)+1, mod1(tilenumber, pt.ncols))
    return ((Point(x, y), tilenumber), (Point(x1, y1), tilenumber + 1))
end

function Base.iterate(pt::Partition, state)
    if state[2] > (pt.nrows * pt.ncols)
        return
    end
    x = state[1].x
    y = state[1].y
    tilenumber = state[2]
    x1 = x + pt.tilewidth
    y1 = y
    if (x1 + pt.tilewidth/2) > (pt.areawidth/2)
        y1 += pt.tileheight
        x1 = -(pt.areawidth/2) + (pt.tilewidth/2)
    end
    pt.currentrow, pt.currentcol = (div(tilenumber-1, pt.ncols)+1, mod1(tilenumber, pt.ncols))
    return ((Point(x, y), tilenumber), (Point(x1, y1), tilenumber + 1))
end

function Base.length(pt::Partition)
    pt.nrows * pt.ncols
end

function Base.getindex(pt::Partition, i::Int)
    1 <= i <= pt.ncols *  pt.nrows || throw(BoundsError(pt, i))
    xcoord = -pt.areawidth/2  + mod1(i, pt.ncols) * pt.tilewidth  - pt.tilewidth/2
    ycoord = -pt.areaheight/2 + (div(i - 1,  pt.ncols) * pt.tileheight) + pt.tileheight/2
    return (Point(xcoord, ycoord), i)
end

"""
    box(tiles::Tiler, n::Integer; action=:none, vertices=false, reversepath=false)
    box(tiles::Tiler, n::Integer, action::Symbol=:none; vertices=false, reversepath=false)

Draw a box in tile `n` of tiles `tiles`.
"""
function box(tiles::Tiler, n::Integer;
        action=:none,
        vertices=false,
        reversepath=false)
    tilew, tileh = tiles.tilewidth, tiles.tileheight
    if vertices == true || action == :none
        box(first(tiles[n]), tilew, tileh, action=:none, vertices=true, reversepath=reversepath)
    else
        box(first(tiles[n]), tilew, tileh, action=action, vertices=false, reversepath=reversepath)
    end
end

box(tiles::Tiler, n::Integer, action::Symbol; vertices=false, reversepath=false) =
        box(tiles, n, action=action, vertices=vertices, reversepath=reversepath)


text:
## Criar uma matrix linha ou coluna
Base.:Matrix(arr::AbstractArray, m::Integer, n::Integer) = begin

    arr_len = length(arr)

    !(m == 1 || n == 1) && throw(DimensionMismatch("Uma das dimensões tem de ser 1"))

    !(m == arr_len || n == arr_len) && throw(DimensionMismatch("Uma das dimensões tem de ter o tamanho do vector"))

    matrix = zeros(Number , m, n)

    for i in 1:m
        for j in 1:n
            matrix[i,j] = m == arr_len ? arr[i] : arr[j]
        end
    end

    matrix
end

verifydimensions(x::T, y::T) where T <: Matrix{Number} = begin
    my, ny = size(y)
    mx, nx = size(x)

    # Verificar se têm o mesmo tamanho
    !(mx in (my, ny) && nx in (my, ny)) && throw(DimensionMismatch("Os vetores têm de ter o mesmo tamanho"))

    (my, ny)
end

max_stars_repo_path: LAFF/Semana 1/helpers.jl

max_issues_repo_path: LAFF/Semana 1/helpers.jl

max_forks_repo_path: LAFF/Semana 1/helpers.jl