text:
using Fermi.Orbitals
function update_amp(T2::Array{T, 4}, newT1::Array{T,2}, newT2::Array{T,4}, foo::Array{T,2}, fov::Array{T,2}, fvv::Array{T,2}, ints::IntegralHelper, alg::A) where { T <: AbstractFloat,
                                                                                A <: CCAlgorithm }

    fill!(newT1, 0.0)
    fill!(newT2, 0.0)

    # Get new amplitudes
    update_T1(T2,newT1,foo,fov,fvv,ints,alg)
    update_T2(T2,newT2,foo,fov,fvv,ints,alg)
end
function BCCD{Ta}(refwfn::RHF, ints::IntegralHelper, newT1::Array{Tb, 2}, newT2::Array{Tc,4}, alg::A; ecT1 = Array{Float64}(undef,0,0), ecT2 = Array{Float64}(undef,0,0,0,0)) where { Ta <: AbstractFloat,
                                                                                Tb <: AbstractFloat,
                                                                                Tc <: AbstractFloat,
                                       A <: CCAlgorithm }

    # Print intro
    Fermi.CoupledCluster.print_header()
    Fermi.CoupledCluster.print_bcc_alg(alg)
    newT1 = convert(Array{Ta},newT1)
    newT2 = convert(Array{Ta},newT2)
    ecT1 = convert(Array{Ta},ecT1)
    ecT2 = convert(Array{Ta},ecT2)

    @output repeat("-",80)*"\n"
    @output "\tBasis: {}\n" ints.bname["primary"]
    @output "\tComputing and Transforming Integrals..."
    tint = @elapsed Fermi.CoupledCluster.compute_integrals(ints,alg)
    @output " done in {} s\n" tint
    for key in keys(ints.cache)
        ints[key] = convert(Array{Ta},ints[key])
    end

    # Process Fock matrix, important for non HF cases
    foo = convert(Array{Ta},ints["FOO"] - Diagonal(ints["FOO"]))
    fvv = convert(Array{Ta},ints["FVV"] - Diagonal(ints["FVV"]))
    fov = convert(Array{Ta},ints["FOV"])

    d = Ta[i - a for i = diag(ints["FOO"]), a = diag(ints["FVV"])]
    D = Ta[i + j - a - b for i = diag(ints["FOO"]), j = diag(ints["FOO"]), a = diag(ints["FVV"]), b = diag(ints["FVV"])]

    # Compute Guess Energy
    oovv = convert(Array{Ta},Fermi.CoupledCluster.compute_oovv(ints,alg))
    Ecc = update_energy(newT1, newT2, fov, ints, alg)
    Eguess = Ecc+refwfn.energy

    #@output "Initial Amplitudes Guess: MP2\n"
    @output "\tGuess Energy:   {:15.10f}\n\n" Ecc
    @output "\tGuess Total Energy:   {:15.10f}\n\n" Ecc+refwfn.energy

    # Start CC iterations

    cc_max_iter = Fermi.CurrentOptions["cc_max_iter"]
    cc_e_conv = Fermi.CurrentOptions["cc_e_conv"]
    cc_max_rms = Fermi.CurrentOptions["cc_max_rms"]
    precision_override = Fermi.CurrentOptions["precision_override"]
    bcc_tol = Fermi.CurrentOptions["bcc_tol"]

    preconv_T1 = Fermi.CurrentOptions["preconv_T1"]
    dp = Fermi.CurrentOptions["cc_damp_ratio"]
    do_diis = Fermi.CurrentOptions["cc_diis"]
    ndiis = Fermi.CurrentOptions["cc_ndiis"]
    diis_start = Fermi.CurrentOptions["diis_start"]
    diis_relax = Fermi.CurrentOptions["cc_diis_relax"]
    prec_selector = Dict{String,Any}(
                                     "half" => Float16,
                                     "single" => Float32,
                                     "double" => Float64
                                    )
    diis_prec = prec_selector[Fermi.CurrentOptions["diis_prec"]]
    do_diis ? DM_T1 = Fermi.DIIS.DIISManager{Ta,diis_prec}(size=ndiis) : nothing
    do_diis ? DM_T2 = Fermi.DIIS.DIISManager{Ta,diis_prec}(size=ndiis) : nothing
    nocc = size(newT2,1)
    nvir = size(newT2,3)

    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]
    @output "\tDropped Occupied Orbitals →  {:3.0f}\n" Int(Fermi.CurrentOptions["drop_occ"])
    @output "\tDropped Virtual Orbitals  →  {:3.0f}\n\n" Int(Fermi.CurrentOptions["drop_vir"])
    @output "\n"*repeat("-",80)*"\n"
    @output "\tOptions:\n"
    @output "\tPrecision? {}\n" Ta
    @output "\tDIIS? {}\n" do_diis
    @output "\tDIIS Vectors? {}\n" (do_diis ? ndiis : 0)
    @output "\tDIIS Precision? {}\n" diis_prec
    @output "\tPreconverge T1 amplitudes? {}\n" preconv_T1
    @output "\tDamping percentage ? {}\n" dp
    @output "\t\tcc_max_iter →  {:3.0d}\n" Int(cc_max_iter)
    @output "\t\tcc_e_conv   →  {:2.0e}\n" cc_e_conv
    @output "\t\tcc_max_rms  →  {:2.0e}\n" cc_max_rms
    if (cc_e_conv < eps(Ta)) && !(precision_override)
        @output "\t⚠ WARNING⚠   cc_e_conv set to less than ϵ ({}) for precision {}\n" eps(Ta) Ta
        @output "\t             CCSD is unlikely to converge to your standards.\n"
        @output "\t             OVERRIDING set convergence criterion.\n"
        @output "\t             Use @set precision_override true to perform the computation as entered.\n"
        cc_e_conv = eps(Ta)
    end
    if (cc_max_rms < eps(Ta)) && !(precision_override)
        @output "\t⚠ WARNING⚠   cc_max_rms set to less than ϵ ({}) for precision {}\n" eps(Ta) Ta
        @output "\t             CCSD is unlikely to converge to your standards.\n"
        @output "\t             OVERRIDING set convergence criterion.\n"
        @output "\t             Use @set precision_override true to perform the computation as entered.\n"
        cc_max_rms = eps(Ta)
    end

    r1 = 1
    r2 = 1
    dE = 1
    rms = 1
    ite = 1
    T1 = deepcopy(newT1)
    T2 = deepcopy(newT2)

    #dfsz,nocc,nvir = size(Bov)

    @output repeat("-",80)*"\n"
    @output "    Starting CC Iterations\n\n"
    preconv_T1 ? T1_time = 0 : nothing
    dummyT2 = zeros(Ta,size(newT2))
    if preconv_T1
        @output "Preconverging T1 amplitudes\n"
        @output "Taking one T2 step\n"
        @output "{:10s}    {: 15s}    {: 12s}    {:12s}    {:10s}\n" "Iteration" "CC Energy" "ΔE" "Max RMS (T1)" "Time (s)"
        t = @elapsed begin
            update_amp(T2, newT1, newT2, foo, fov, fvv, ints, alg)

            #apply external correction
            apply_ec(newT1,dummyT2,ecT1,ecT2)

            # Apply resolvent
            newT1 ./= d
            newT2 ./= D

            # Compute residues
            r1 = sqrt(sum((newT1 .- T1).^2))/length(T1)
            r2 = sqrt(sum((newT2 .- T2).^2))/length(T2)

            if do_diis
                e1 = (newT1 - T1)
                e2 = (newT2 - T2)
                push!(DM_T1,newT1,e1)
                push!(DM_T2,newT2,e2)
            end

            newT1 .= (1-dp)*newT1 .+ dp*T1
            newT2 .= (1-dp)*newT2 .+ dp*T2
        end
        T1_time += t

        rms = max(r1,r2)
        oldE = Ecc
        Ecc = update_energy(newT1, newT2, fov, ints, alg)
        dE = Ecc - oldE
        @output "    {:<5}    {:<15.10f}    {:<12.10f}    {:<12.10f}    {:<10.5f}\n" "pre" Ecc dE rms t

        while abs(dE) > cc_e_conv || rms > cc_max_rms
            if ite > cc_max_iter
                @output "\n⚠  CC Equations did not converge in {:1.0d} iterations.\n" cc_max_iter
                break
            end
            t = @elapsed begin
                T1 .= newT1
                T2 .= newT2
                update_T1(T2,newT1,foo,fov,fvv,ints,alg)
                apply_ec(newT1,dummyT2,ecT1,ecT2)
                newT1 ./= d
                if do_diis
                    e1 = newT1 - T1
                    push!(DM_T1,newT1,e1)
                    newT1 = Fermi.DIIS.extrapolate(DM_T1)
                end

                # Compute residues
                r1 = sqrt(sum((newT1 .- T1).^2))/length(T1)

                newT1 .= (1-dp)*newT1 .+ dp*T1

                rms = r1
                oldE = Ecc
                Ecc = update_energy(newT1, newT2, fov, ints, alg)
                dE = Ecc - oldE
            end
            T1_time += t
            dE > 0 ? sign = "+" : sign = "-"
            @output "    {:<5.0d}    {:<15.10f}    {}{:>12.10f}    {:<12.10f}    {:<10.5f}\n" ite Ecc sign abs(dE) rms t
            ite += 1
        end
        @output "\nT1 pre-convergence took {}s\n" T1_time
    end

    ints.orbs["BCC"] = deepcopy(ints.orbs["FU"])
    Fermi.Orbitals.activate!(ints.orbs,"BCC")
    ints.orbs.frozencore = drop_occ
    ints.orbs.frozenvir = drop_vir
    C = ints.orbs["[FU]"]
    #display(C)

    dE = 1
    rms = 1
    ite = 1

    do_diis ? DM_T1 = Fermi.DIIS.DIISManager{Ta,diis_prec}(size=ndiis) : nothing
    do_diis ? DM_T2 = Fermi.DIIS.DIISManager{Ta,diis_prec}(size=ndiis) : nothing
    if preconv_T1
        @output "Including T2 update\n"
    end

    main_time = 0
    @output "{:10s}    {: 15s}    {: 12s}    {:12s}    {:10s} {:3s}\n" "Iteration" "CC Energy" "ΔE" "Max RMS" "Time (s)" "BCC"
    relax = 1
    function myexp(A::Array,N=6)
        sum = zeros(size(A))
        for n=0:N
            sum += A^n/factorial(n)
        end
        sum
    end

    mT1 = 0
    while (abs(dE) > cc_e_conv || rms > cc_max_rms) || mT1 > bcc_tol
        if ite > cc_max_iter
            @output "\n⚠  CC Equations did not converge in {:1.0d} iterations.\n" cc_max_iter
            break
        end
        t = @elapsed begin

            bcc = false
            T1 .= newT1
            T2 .= newT2
            d = Ta[i - a for i = diag(ints["FOO"]), a = diag(ints["FVV"])]
            D = Ta[i + j - a - b for i = diag(ints["FOO"]), j = diag(ints["FOO"]), a = diag(ints["FVV"]), b = diag(ints["FVV"])]

            update_T2(T2,newT2,foo,fov,fvv,ints,alg)
            update_T1(T2,newT1,foo,fov,fvv,ints,alg)
            newT2 ./= D
            newT1 ./= d
            oldE = Ecc
            Ecc = update_energy(newT1, newT2, fov, ints, alg)

            oldC = ints.orbs["[FU]"]

            #update_amp(T2, newT1, newT2, foo, fov, fvv, ints, alg)
            apply_ec(newT1,newT2,ecT1,ecT2)
            # Apply resolvent
            if do_diis
                relax -= 1
                e1 = (newT1 - T1)
                e2 = (newT2 - T2)
                push!(DM_T1,newT1,e1)
                push!(DM_T2,newT2,e2)
                if relax == 0
                    newT2 = Fermi.DIIS.extrapolate(DM_T2)
                    newT1 = Fermi.DIIS.extrapolate(DM_T1)
                    relax = diis_relax
                end
            end

            newT1 .= (1-dp)*newT1 .+ dp*T1
            newT2 .= (1-dp)*newT2 .+ dp*T2

            #println(maximum(abs.(T1)))

            # Compute residues
            r1 = sqrt(sum((newT1 .- T1).^2))/length(T1)
            r2 = sqrt(sum((newT2 .- T2).^2))/length(T2)

            mT1 = maximum(abs.(newT1))

            #if there's no significant rotation, lets just skip the integral transform
            if mT1 >= bcc_tol #&& max(r1,r2) < 1E-10#&& ite >= 3
                do_diis ? DM_T1 = Fermi.DIIS.DIISManager{Ta,diis_prec}(size=ndiis) : nothing
                do_diis ? DM_T2 = Fermi.DIIS.DIISManager{Ta,diis_prec}(size=ndiis) : nothing
                bcc = true
                X = Array([ zeros(nocc,nocc) zeros(size(T1)); newT1' zeros(nvir,nvir)])
                U = exp(X - X')
                Fermi.Orbitals.rotate!(ints.orbs,1.0*U,fc=drop_occ,fv=drop_vir)
                delete_integrals(ints,alg)
                Fermi.Orbitals.matchphase(ints.orbs,oldC)
                foo = convert(Array{Ta},ints["FOO"] - Diagonal(ints["FOO"]))
                fvv = convert(Array{Ta},ints["FVV"] - Diagonal(ints["FVV"]))
                fov = convert(Array{Ta},ints["FOV"])
            end

        end
        rms = max(r1,r2)
        dE = Ecc - oldE
        main_time += t
        dE > 0 ? sign = "+" : sign = "-"
        @output "    {:<5.0d}    {:<15.10f}    {}{:>12.10f}    {:<12.10f}    {:<10.5f} {:3s}\n" ite Ecc sign abs(dE) rms t bcc
        ite += 1
    end
    @output "\nMain BCC iterations done in {}s\n" main_time

    # Converged?
    if abs(dE) < cc_e_conv && rms < cc_max_rms
        @output "\n 🍾 Equations Converged!\n"
    end
    @output "\n⇒ Final BCCD Energy:     {:15.10f}\n" Ecc+refwfn.energy
    @output "\n  Final T1 norm:         {:15.10f}\n" sqrt(sum(newT1.^2))
    @output repeat("-",80)*"\n"

    return RCCSD{Ta}(Eguess, Ecc+refwfn.energy, Fermi.MemTensor(newT1), Fermi.MemTensor(newT2))
end

function apply_ec(T1,T2,ecT1,ecT2)
    s1,s2,s3,s4 = size(ecT2)
    r1 = 1:s1
    r2 = 1:s2
    r3 = 1:s3
    r4 = 1:s4
    T1[r1,r3] .+= ecT1
    T2[r1,r2,r3,r4] .+= ecT2
end


text:
function movements!(r:: Robot, side::HorizonSide)
    """Move robot to reach the end."""
    c = 0
    while !isborder(r, side)
        move!(r, side)
        c += 1
    end
    return c
end

function movements!(r::Robot, side::HorizonSide, num::Integer)
    """Move robot num steps in side direction."""
    while num > 0
        move!(r, side)
        num -= 1
    end
end

function to_corner(r:: Robot, side1::HorizonSide, side2::HorizonSide)
    """Move robot to the corner of side1 & side2"""
    s1s2 = [0, 0]
    while !isborder(r, side1) || !isborder(r, side2)
        s1s2[1] += movements!(r, side1)
        s1s2[2] += movements!(r, side2)
    end
    return s1s2
end

sw, no = to_corner(r, Sud, West), [0, 0]

for sides in [[Sud, West], [Nord, Ost]]
    for i = 1:2
        movements!(r, HorizonSide((Int(sides[i]) + 2) % 4), abs(sw[i] - no[i]))
        putmarker!(r)
        movements!(r, sides[i], abs(sw[i] - no[i]))
    end
    global no = to_corner(r, Nord, Ost)
end

movements!(r, Sud, abs(no[1] - sw[1]))
movements!(r, West, abs(no[2] - sw[2]))

ext: jl

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
using MOCNeutronTransport

# Model setup
# ----------------------------------------------------------------------------------------
gmsh.initialize()
gmsh.option.set_number("Geometry.OCCParallel", 1)
gmsh.option.set_number("General.NumThreads", 32)

# Import CAD file
gmsh.merge("/home/kcvaughn/Downloads/Core Assembly - Test Layout 2 - 6-30-20 Export.STEP")

# Shift the model to the origin of the xy-plane
dim_tags = gmsh.model.get_entities(3)
# Rotate the entities in dim_tags in the built-in CAD representation by angle radians
# around the axis of revolution defined by the point (x, y, z) and the direction (ax, ay, az).
# gmsh.model.occ.rotate(dim_tags, x, y, z, ax, ay, az, angle)
gmsh.model.occ.rotate(dim_tags, 0, 0, 0, 1, 0, 0, π/2)
# Translate it so that the corner is at the origin
#gmsh.model.occ.translate(dim_tags, 37.3225 - 0.0099 + 10*grid_offset,
#                                   44.09 - 0.00174209 + 10*grid_offset,
#                                   -436.5625)
# Model is in mm. Convert to cm by shrinking everything by 1/10
#gmsh.model.occ.dilate(dim_tags, 0, 0, 0, 1//10, 1//10, 0)

rect_tag = gmsh.model.occ.add_rectangle(-500.0, -500.0, 300, 1000.0, 1000.0)

# Synchronize the CAD and Gmsh models
gmsh.model.occ.synchronize()

#Geometry.OCCParallel
# Verify results
# gmsh.fltk.run()
intersections = gmsh.model.occ.intersect([(Int32(2), rect_tag)], dim_tags)
gmsh.model.occ.synchronize()
gmsh.write("FNR_slice.step")
#gmsh.finalize()

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
""""""
function curl(process, downloadpath)
    mkpath(downloadpath)
    cookiepath = joinpath(downloadpath, "cookie.txt")
    searchpath = joinpath(downloadpath, "search.html")

    cookiecmd = `curl http://awd.mv-regierung.de/$process/index.php -c $cookiepath`
    searchpagecmd = `curl http://awd.mv-regierung.de/$process/anz_abschn.php -b $cookiepath`

    run(cookiecmd)
    run(pipeline(searchpagecmd; stdout=searchpath))

    open(searchpath) do f
        doc = parsehtml(read(f, String))
        options = eachmatch(Selector("select option"), doc.root)
        for option in options
            id = getattr(option, "value")
            id == "" && continue
            pagecmd = `curl http://awd.mv-regierung.de/$process/anz_abschn.php --compressed --data auswahl=$id -b $cookiepath`
            run(pipeline(pagecmd; stdout=joinpath(downloadpath, "$id.html")))
        end
    end
end

""""""
function create(dbname, tablename)
    dbconn = SQLite.DB("$dbname.sqlite")
    SQLite.execute!(dbconn, """
        CREATE TABLE IF NOT EXISTS $tablename (
            id                  INTEGER PRIMARY KEY,
            number              INTEGER,
            author              TEXT,
            ref                 TEXT,
            content             TEXT NOT NULL,
            assessment_content  TEXT
        );
    """)
end

""""""
function transform(path)
    values = []

    open(path) do f
        doc = parsehtml(read(f, String))
        tables = eachmatch(Selector(".table_wrap"), doc.root)
        if length(tables) != 3
            return values
        end

        contributions = eachmatch(Selector(".tblrow"), tables[2])
        for i in 2:length(contributions)
            contribution = contributions[i]
            parts = eachmatch(Selector(".listitem"), contribution)

            participant = nodeText(parts[1])
            ref = nodeText(parts[2])
            content = nodeText(parts[3])
            assessment = nodeText(parts[4])

            m = match(r"lfd\.\sNr\.:\s(\d+)(.+)", participant)
            if isnothing(m)
                continue
            end

            push!(values, m[1], m[2], ref, content, assessment)
        end
    end

    return values
end

""""""
function save(values, dbname, tablename; batchsz=50)
    nvalues = length(values)
    (nvalues > 0 && nvalues % 5 == 0) || return

    ncontribs = div(nvalues, 5)
    niter = div(ncontribs, batchsz)
    restsz = ncontribs % batchsz
    placeholder = "(?, ?, ?, ?, ?)"
    dbconn = SQLite.DB("$dbname.sqlite")
    for i in 1:niter
        SQLite.Query(dbconn, """
            INSERT INTO $tablename (number, author, ref, content, assessment_content)
            VALUES $(((placeholder * ", ") ^ (batchsz - 1)) * placeholder);
        """; values=values[1 + 5 * batchsz * (i - 1):5 * batchsz * i])
    end
    if restsz > 0
        SQLite.Query(dbconn, """
            INSERT INTO $tablename (number, author, ref, content, assessment_content)
            VALUES $(((placeholder * ", ") ^ (restsz - 1)) * placeholder);
        """; values=values[end - 5 * restsz + 1:end])
    end
end

""""""
function transformall(path, dbname, tablename)
    relativepath = relpath(path)
    paths = map(p -> joinpath(relativepath, p), readdir(relativepath))
    filter!(p ->
        isfile(p) && endswith(p, ".html") && !occursin("search", p)
    , paths)
    for p in paths
        trans = transform(p)
        save(trans, dbname, tablename)
    end
end

""""""
function scrape(process, downloadpath, dbname, tablename;
                dodownload=true, docreate=true, dotransform=true)
    dodownload && curl(process, downloadpath)
    docreate && create(dbname, tablename)
    dotransform && transformall(downloadpath, dbname, tablename)
end

""""""
function mwp(download=true, create=true, transform=true)
    processes = ["lep_2009_01", "rrep_2010_wm2", "rrep_2011_wm_3",
            "rrep_2011_wm_4", "rrep_2015_01", "lep_2009_03", "rrep_2011_mm_1",
            "lep_2009_02", "rrep_2010_mv_2", "rrep_2010_ms_3", "rrep_2011_ms_1",
            "lep_2008_01", "rrep_2009_02", "rrep_2013_01", "rrep_2015_02",
            "rrep_2015_03", "lep_2014_01", "lep_2016_01"]
    for p in processes
        scrape(p, "data/mwp/$p", p, "contribution"; dodownload=download,
                docreate=create, dotransform=transform)
    end
end


text:
@testset "5.2.4 (f x)^m (d+e x^2)^p (a+b arccos(c x))^n" begin
    (a, b, c, d, e, x, ) = @variables a b c d e x

    #= ::Package:: =#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(f*x)^m*(d+e*x^2)^p*(a+b*acos(c*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(d-c^2*d*x^2)^p*(a+b*acos(c*x))^1=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^3*(a + b*acos(c*x))/(d - c^2*d*x^2), x, 8, (b*x*sqrt(1 - c^2*x^2))/(4*c^3*d) - (x^2*(a + b*acos(c*x)))/(2*c^2*d) + (im*(a + b*acos(c*x))^2)/(2*b*c^4*d) - (b*asin(c*x))/(4*c^4*d) - ((a + b*acos(c*x))*log(1 - ℯ^(2*im*acos(c*x))))/(c^4*d) + (im*b*PolyLog(2, ℯ^(2*im*acos(c*x))))/(2*c^4*d)]
    @test_int [x^2*(a + b*acos(c*x))/(d - c^2*d*x^2), x, 8, (b*sqrt(1 - c^2*x^2))/(c^3*d) - (x*(a + b*acos(c*x)))/(c^2*d) + (2*(a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/(c^3*d) - (im*b*PolyLog(2, -ℯ^(im*acos(c*x))))/(c^3*d) + (im*b*PolyLog(2, ℯ^(im*acos(c*x))))/(c^3*d)]
    @test_int [x^1*(a + b*acos(c*x))/(d - c^2*d*x^2), x, 5, (im*(a + b*acos(c*x))^2)/(2*b*c^2*d) - ((a + b*acos(c*x))*log(1 - ℯ^(2*im*acos(c*x))))/(c^2*d) + (im*b*PolyLog(2, ℯ^(2*im*acos(c*x))))/(2*c^2*d)]
    @test_int [x^0*(a + b*acos(c*x))/(d - c^2*d*x^2), x, 6, (2*(a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/(c*d) - (im*b*PolyLog(2, -ℯ^(im*acos(c*x))))/(c*d) + (im*b*PolyLog(2, ℯ^(im*acos(c*x))))/(c*d)]
    @test_int [(a + b*acos(c*x))/(x^1*(d - c^2*d*x^2)), x, 7, (2*(a + b*acos(c*x))*atanh(ℯ^(2*im*acos(c*x))))/d - (im*b*PolyLog(2, -ℯ^(2*im*acos(c*x))))/(2*d) + (im*b*PolyLog(2, ℯ^(2*im*acos(c*x))))/(2*d)]
    @test_int [(a + b*acos(c*x))/(x^2*(d - c^2*d*x^2)), x, 10, -((a + b*acos(c*x))/(d*x)) + (2*c*(a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/d + (b*c*atanh(sqrt(1 - c^2*x^2)))/d - (im*b*c*PolyLog(2, -ℯ^(im*acos(c*x))))/d + (im*b*c*PolyLog(2, ℯ^(im*acos(c*x))))/d]
    @test_int [(a + b*acos(c*x))/(x^3*(d - c^2*d*x^2)), x, 9, (b*c*sqrt(1 - c^2*x^2))/(2*d*x) - (a + b*acos(c*x))/(2*d*x^2) + (2*c^2*(a + b*acos(c*x))*atanh(ℯ^(2*im*acos(c*x))))/d - (im*b*c^2*PolyLog(2, -ℯ^(2*im*acos(c*x))))/(2*d) + (im*b*c^2*PolyLog(2, ℯ^(2*im*acos(c*x))))/(2*d)]

    @test_int [x^4*(a + b*acos(c*x))/(d - c^2*d*x^2)^2, x, 12, b/(2*c^5*d^2*sqrt(1 - c^2*x^2)) - (b*sqrt(1 - c^2*x^2))/(c^5*d^2) + (3*x*(a + b*acos(c*x)))/(2*c^4*d^2) + (x^3*(a + b*acos(c*x)))/(2*c^2*d^2*(1 - c^2*x^2)) - (3*(a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/(c^5*d^2) + (3*im*b*PolyLog(2, -ℯ^(im*acos(c*x))))/(2*c^5*d^2) - (3*im*b*PolyLog(2, ℯ^(im*acos(c*x))))/(2*c^5*d^2)]
    @test_int [x^3*(a + b*acos(c*x))/(d - c^2*d*x^2)^2, x, 8, (b*x)/(2*c^3*d^2*sqrt(1 - c^2*x^2)) + (x^2*(a + b*acos(c*x)))/(2*c^2*d^2*(1 - c^2*x^2)) - (im*(a + b*acos(c*x))^2)/(2*b*c^4*d^2) - (b*asin(c*x))/(2*c^4*d^2) + ((a + b*acos(c*x))*log(1 - ℯ^(2*im*acos(c*x))))/(c^4*d^2) - (im*b*PolyLog(2, ℯ^(2*im*acos(c*x))))/(2*c^4*d^2)]
    @test_int [x^2*(a + b*acos(c*x))/(d - c^2*d*x^2)^2, x, 8, b/(2*c^3*d^2*sqrt(1 - c^2*x^2)) + (x*(a + b*acos(c*x)))/(2*c^2*d^2*(1 - c^2*x^2)) - ((a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/(c^3*d^2) + (im*b*PolyLog(2, -ℯ^(im*acos(c*x))))/(2*c^3*d^2) - (im*b*PolyLog(2, ℯ^(im*acos(c*x))))/(2*c^3*d^2)]
    @test_int [x^1*(a + b*acos(c*x))/(d - c^2*d*x^2)^2, x, 2, (b*x)/(2*c*d^2*sqrt(1 - c^2*x^2)) + (a + b*acos(c*x))/(2*c^2*d^2*(1 - c^2*x^2))]
    @test_int [x^0*(a + b*acos(c*x))/(d - c^2*d*x^2)^2, x, 8, b/(2*c*d^2*sqrt(1 - c^2*x^2)) + (x*(a + b*acos(c*x)))/(2*d^2*(1 - c^2*x^2)) + ((a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/(c*d^2) - (im*b*PolyLog(2, -ℯ^(im*acos(c*x))))/(2*c*d^2) + (im*b*PolyLog(2, ℯ^(im*acos(c*x))))/(2*c*d^2)]
    @test_int [(a + b*acos(c*x))/(x^1*(d - c^2*d*x^2)^2), x, 9, (b*c*x)/(2*d^2*sqrt(1 - c^2*x^2)) + (a + b*acos(c*x))/(2*d^2*(1 - c^2*x^2)) + (2*(a + b*acos(c*x))*atanh(ℯ^(2*im*acos(c*x))))/d^2 - (im*b*PolyLog(2, -ℯ^(2*im*acos(c*x))))/(2*d^2) + (im*b*PolyLog(2, ℯ^(2*im*acos(c*x))))/(2*d^2)]
    @test_int [(a + b*acos(c*x))/(x^2*(d - c^2*d*x^2)^2), x, 13, (b*c)/(2*d^2*sqrt(1 - c^2*x^2)) - (a + b*acos(c*x))/(d^2*x*(1 - c^2*x^2)) + (3*c^2*x*(a + b*acos(c*x)))/(2*d^2*(1 - c^2*x^2)) + (3*c*(a + b*acos(c*x))*atanh(ℯ^(im*acos(c*x))))/d^2 + (b*c*atanh(sqrt(1 - c^2*x^2)))/d^2 - (3*im*b*c*PolyLog(2, -ℯ^(im*acos(c*x))))/(2*d^2) + (3*im*b*c*PolyLog(2, ℯ^(im*acos(c*x))))/(2*d^2)]
    @test_int [(a + b*acos(c*x))/(x^3*(d - c^2*d*x^2)^2), x, 12, (b*c)/(2*d^2*x*sqrt(1 - c^2*x^2)) + (c^2*(a + b*acos(c*x)))/(d^2*(1 - c^2*x^2)) - (a + b*acos(c*x))/(2*d^2*x^2*(1 - c^2*x^2)) + (4*c^2*(a + b*acos(c*x))*atanh(ℯ^(2*im*acos(c*x))))/d^2 - (im*b*c^2*PolyLog(2, -ℯ^(2*im*acos(c*x))))/d^2 + (im*b*c^2*PolyLog(2, ℯ^(2*im*acos(c*x))))/d^2]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(d+e*x^2)^p*(a+b*acos(c*x))^1=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(d + e*x^2)*(a + b*acos(c*x))*x^3, x, 6, -((b*(9*c^2*d + 5*e)*x*sqrt(1 - c^2*x^2))/(96*c^5)) - (b*(9*c^2*d + 5*e)*x^3*sqrt(1 - c^2*x^2))/(144*c^3) - (b*e*x^5*sqrt(1 - c^2*x^2))/(36*c) + (1/4)*d*x^4*(a + b*acos(c*x)) + (1/6)*e*x^6*(a + b*acos(c*x)) + (b*(9*c^2*d + 5*e)*asin(c*x))/(96*c^6)]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))*x^2, x, 5, -((b*(5*c^2*d + 3*e)*sqrt(1 - c^2*x^2))/(15*c^5)) + (b*(5*c^2*d + 6*e)*(1 - c^2*x^2)^(3/2))/(45*c^5) - (b*e*(1 - c^2*x^2)^(5/2))/(25*c^5) + (1/3)*d*x^3*(a + b*acos(c*x)) + (1/5)*e*x^5*(a + b*acos(c*x))]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))*x^1, x, 4, -((3*b*(2*c^2*d + e)*x*sqrt(1 - c^2*x^2))/(32*c^3)) - (b*x*sqrt(1 - c^2*x^2)*(d + e*x^2))/(16*c) + ((d + e*x^2)^2*(a + b*acos(c*x)))/(4*e) + (b*(8*c^4*d^2 + 8*c^2*d*e + 3*e^2)*asin(c*x))/(32*c^4*e)]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))*x^0, x, 4, -((b*(3*c^2*d + e)*sqrt(1 - c^2*x^2))/(3*c^3)) + (b*e*(1 - c^2*x^2)^(3/2))/(9*c^3) + d*x*(a + b*acos(c*x)) + (1/3)*e*x^3*(a + b*acos(c*x))]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))/x^1, x, 12, -((b*e*x*sqrt(1 - c^2*x^2))/(4*c)) + (1/2)*e*x^2*(a + b*acos(c*x)) + (b*e*asin(c*x))/(4*c^2) + (1/2)*im*b*d*asin(c*x)^2 - b*d*asin(c*x)*log(1 - ℯ^(2*im*asin(c*x))) + d*(a + b*acos(c*x))*log(x) + b*d*asin(c*x)*log(x) + (1/2)*im*b*d*PolyLog(2, ℯ^(2*im*asin(c*x)))]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))/x^2, x, 5, -((b*e*sqrt(1 - c^2*x^2))/c) - (d*(a + b*acos(c*x)))/x + e*x*(a + b*acos(c*x)) + b*c*d*atanh(sqrt(1 - c^2*x^2))]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))/x^3, x, 10, (b*c*d*sqrt(1 - c^2*x^2))/(2*x) - (d*(a + b*acos(c*x)))/(2*x^2) + (1/2)*im*b*e*asin(c*x)^2 - b*e*asin(c*x)*log(1 - ℯ^(2*im*asin(c*x))) + e*(a + b*acos(c*x))*log(x) + b*e*asin(c*x)*log(x) + (1/2)*im*b*e*PolyLog(2, ℯ^(2*im*asin(c*x)))]
    @test_int [(d + e*x^2)*(a + b*acos(c*x))/x^4, x, 6, (b*c*d*sqrt(1 - c^2*x^2))/(6*x^2) - (d*(a + b*acos(c*x)))/(3*x^3) - (e*(a + b*acos(c*x)))/x + (1/6)*b*c*(c^2*d + 6*e)*atanh(sqrt(1 - c^2*x^2))]

    @test_int [acos(a*x)*(c + d*x^2)^2, x, 5, -(((15*a^4*c^2 + 10*a^2*c*d + 3*d^2)*sqrt(1 - a^2*x^2))/(15*a^5)) + (2*d*(5*a^2*c + 3*d)*(1 - a^2*x^2)^(3/2))/(45*a^5) - (d^2*(1 - a^2*x^2)^(5/2))/(25*a^5) + c^2*x*acos(a*x) + (2/3)*c*d*x^3*acos(a*x) + (1/5)*d^2*x^5*acos(a*x)]

    @test_int [acos(a*x)*(c + d*x^2)^3, x, 5, -(((35*a^6*c^3 + 35*a^4*c^2*d + 21*a^2*c*d^2 + 5*d^3)*sqrt(1 - a^2*x^2))/(35*a^7)) + (d*(35*a^4*c^2 + 42*a^2*c*d + 15*d^2)*(1 - a^2*x^2)^(3/2))/(105*a^7) - (3*d^2*(7*a^2*c + 5*d)*(1 - a^2*x^2)^(5/2))/(175*a^7) + (d^3*(1 - a^2*x^2)^(7/2))/(49*a^7) + c^3*x*acos(a*x) + c^2*d*x^3*acos(a*x) + (3/5)*c*d^2*x^5*acos(a*x) + (1/7)*d^3*x^7*acos(a*x)]

    @test_int [acos(a*x)*(c + d*x^2)^4, x, 5, -(((315*a^8*c^4 + 420*a^6*c^3*d + 378*a^4*c^2*d^2 + 180*a^2*c*d^3 + 35*d^4)*sqrt(1 - a^2*x^2))/(315*a^9)) + (4*d*(105*a^6*c^3 + 189*a^4*c^2*d + 135*a^2*c*d^2 + 35*d^3)*(1 - a^2*x^2)^(3/2))/(945*a^9) - (2*d^2*(63*a^4*c^2 + 90*a^2*c*d + 35*d^2)*(1 - a^2*x^2)^(5/2))/(525*a^9) + (4*d^3*(9*a^2*c + 7*d)*(1 - a^2*x^2)^(7/2))/(441*a^9) - (d^4*(1 - a^2*x^2)^(9/2))/(81*a^9) + c^4*x*acos(a*x) + (4/3)*c^3*d*x^3*acos(a*x) + (6/5)*c^2*d^2*x^5*acos(a*x) + (4/7)*c*d^3*x^7*acos(a*x) + (1/9)*d^4*x^9*acos(a*x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [acos(a*x)/(c + d*x^2)^1, x, 18, (acos(a*x)*log(1 - (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d))))/(2*sqrt(-c)*sqrt(d)) - (acos(a*x)*log(1 + (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d))))/(2*sqrt(-c)*sqrt(d)) + (acos(a*x)*log(1 - (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d))))/(2*sqrt(-c)*sqrt(d)) - (acos(a*x)*log(1 + (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d))))/(2*sqrt(-c)*sqrt(d)) + (im*PolyLog(2, -((sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d)))))/(2*sqrt(-c)*sqrt(d)) - (im*PolyLog(2, (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d))))/(2*sqrt(-c)*sqrt(d)) + (im*PolyLog(2, -((sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d)))))/(2*sqrt(-c)*sqrt(d)) - (im*PolyLog(2, (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d))))/(2*sqrt(-c)*sqrt(d))]

    @test_int [acos(a*x)/(c + d*x^2)^2, x, 26, -(acos(a*x)/(4*c*sqrt(d)*(sqrt(-c) - sqrt(d)*x))) + acos(a*x)/(4*c*sqrt(d)*(sqrt(-c) + sqrt(d)*x)) - (a*atanh((sqrt(d) - a^2*sqrt(-c)*x)/(sqrt(a^2*c + d)*sqrt(1 - a^2*x^2))))/(4*c*sqrt(d)*sqrt(a^2*c + d)) - (a*atanh((sqrt(d) + a^2*sqrt(-c)*x)/(sqrt(a^2*c + d)*sqrt(1 - a^2*x^2))))/(4*c*sqrt(d)*sqrt(a^2*c + d)) - (acos(a*x)*log(1 - (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d))))/(4*(-c)^(3/2)*sqrt(d)) + (acos(a*x)*log(1 + (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d))))/(4*(-c)^(3/2)*sqrt(d)) - (acos(a*x)*log(1 - (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d))))/(4*(-c)^(3/2)*sqrt(d)) + (acos(a*x)*log(1 + (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d))))/(4*(-c)^(3/2)*sqrt(d)) - (im*PolyLog(2, -((sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d)))))/(4*(-c)^(3/2)*sqrt(d)) + (im*PolyLog(2, (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) - im*sqrt(a^2*c + d))))/(4*(-c)^(3/2)*sqrt(d)) - (im*PolyLog(2, -((sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d)))))/(4*(-c)^(3/2)*sqrt(d)) + (im*PolyLog(2, (sqrt(d)*ℯ^(im*acos(a*x)))/(a*sqrt(-c) + im*sqrt(a^2*c + d))))/(4*(-c)^(3/2)*sqrt(d))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d+e*x^2)^(p/2)*(a+b*acos(c*x))^1=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [acos(a*x)*(c + d*x^2)^(1/2), x, 0, Unintegrable(sqrt(c + d*x^2)*acos(a*x), x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [acos(a*x)/(c + d*x^2)^(1/2), x, 0, Unintegrable(acos(a*x)/sqrt(c + d*x^2), x)]

    @test_int [acos(a*x)/(c + d*x^2)^(3/2), x, 6, (x*acos(a*x))/(c*sqrt(c + d*x^2)) - atan((sqrt(d)*sqrt(1 - a^2*x^2))/(a*sqrt(c + d*x^2)))/(c*sqrt(d))]

    @test_int [acos(a*x)/(c + d*x^2)^(5/2), x, 7, -((a*sqrt(1 - a^2*x^2))/(3*c*(a^2*c + d)*sqrt(c + d*x^2))) + (x*acos(a*x))/(3*c*(c + d*x^2)^(3/2)) + (2*x*acos(a*x))/(3*c^2*sqrt(c + d*x^2)) - (2*atan((sqrt(d)*sqrt(1 - a^2*x^2))/(a*sqrt(c + d*x^2))))/(3*c^2*sqrt(d))]

    @test_int [acos(a*x)/(c + d*x^2)^(7/2), x, 8, -((a*sqrt(1 - a^2*x^2))/(15*c*(a^2*c + d)*(c + d*x^2)^(3/2))) - (2*a*(3*a^2*c + 2*d)*sqrt(1 - a^2*x^2))/(15*c^2*(a^2*c + d)^2*sqrt(c + d*x^2)) + (x*acos(a*x))/(5*c*(c + d*x^2)^(5/2)) + (4*x*acos(a*x))/(15*c^2*(c + d*x^2)^(3/2)) + (8*x*acos(a*x))/(15*c^3*sqrt(c + d*x^2)) - (8*atan((sqrt(d)*sqrt(1 - a^2*x^2))/(a*sqrt(c + d*x^2))))/(15*c^3*sqrt(d))]
    end

max_stars_repo_path: test/5 Inverse trig functions/5.2 Inverse cosine/5.2.4 (f x)^m (d+e x^2)^p (a+b arccos(c x))^n.jl


text:
function train_Robot_Conv(model,predict,loss,w; lr=1e-5,epochs=100,batchsize=24,testInterval=1)
    @time begin
        w = map(KnetArray,w)
        m = map(zeros,w)
        v = map(zeros,w)
        beta1 = 0.9; beta2 = 0.999; adamT = 1; eps = 1e-8;
        printlnFlush("lr=$(lr)\tepochs=$(epochs)\tbatchsize=$(batchsize)\tmodel=Conv_$(model)\ttestInterval=$(testInterval)")
        printlnFlush(string("Number of parameters: ",sum(map(length,w))," (",Int(floor(sum(map(length,w))/1e6))," million)"))
        printlnFlush("Loading data...")
        @time xtrain,ytrain = loadData_Robot(model,"train")
        @time xtest,ytest = loadData_Robot(model,"test")
        printlnFlush("Training...")
        grdnt = gradloss(loss)
        trainbatchCount = Int64(ceil(size(ytrain,4)/batchsize))
        testbatchCount = Int64(ceil(size(ytest,4)/batchsize))
        maxacc = 0
        for epoch=1:epochs
            trainloss = testloss = trainncorrect = testncorrect = trainncount = testncount = 0
            for batch=1:trainbatchCount
                printUpdate(@sprintf("Batch: %5d/%-5d",batch,trainbatchCount))
                #Calculate current indices
                startIndex = 1+(batch-1)*batchsize
                endIndex = min(batch*batchsize,size(ytrain,4))
                #Minibatch
                batchx = KnetArray(xtrain[:,:,:,startIndex:endIndex])
                batchy = ytrain[:,:,:,startIndex:endIndex]
                g,l = grdnt(w,batchx,batchy)
                trainloss += l
                #Adam
                for i in 1:length(w)
                    m[i] = beta1*m[i] + (1-beta1)*g[i]
                    v[i] = beta2*v[i] + (1-beta2)*(g[i].*g[i])
                    mhat = m[i] ./ (1-beta1^adamT)
                    vhat = v[i] ./ (1-beta2^adamT)
                    w[i] = w[i] - (lr ./ (sqrt(vhat)+eps)).*mhat
                end
                adamT += 1
            end
            #Testing accuracy
            if epoch%testInterval == 0 || epoch == 1
                for batch=1:testbatchCount
                    printUpdate(@sprintf("Batch: %5d/%-5d",batch,testbatchCount))
                    #Calculate current indices
                    startIndex = 1+(batch-1)*batchsize
                    endIndex = min(batch*batchsize,size(ytest,4))
                    #Minibatch
                    batchx = KnetArray(xtest[:,:,:,startIndex:endIndex])
                    batchy = ytest[:,:,:,startIndex:endIndex]
                    #Accuracy
                    ypred = Array(predict(w,batchx))
                    testloss += pixelwiseSoftloss(ypred,batchy)
                    testAcc = pixelwiseAccuracy(ypred,batchy)
                    testncorrect += testAcc[1]
                    testncount += testAcc[2]
                end
            end
            testaccuracy = testncorrect/testncount*100
            if isnan(testaccuracy)
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f",epoch,trainloss/trainbatchCount))
            else
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f\tTest Loss: %0.10f\tTest Accuracy: %0.4f%%",
                    epoch,trainloss/trainbatchCount,testloss/testbatchCount,testaccuracy))
            end
            if testaccuracy > maxacc
                maxacc = testaccuracy
                saveNetwork(w,"Robot_Conv_$(model)")
            end
            if epoch == epochs
                printlnFlush(@sprintf("Maximum Test Accuracy: %0.4f%%",maxacc))
            end
        end
    end
end

function train_Robot_LSTM(model,predict,loss,w; lr=1e-5,epochs=100,testInterval=1)
    @time begin
        D = 1
        w = map(KnetArray,w)
        m = map(zeros,w)
        v = map(zeros,w)
        h = KnetArray(zeros(Float32,28,28,512,1))
        c = KnetArray(zeros(Float32,28,28,512,1))
        beta1 = 0.9; beta2 = 0.999; adamT = 1; eps = 1e-8;
        printlnFlush("lr=$(lr)\tepochs=$(epochs)\tmodel=LSTM_$(model)\ttestInterval=$(testInterval)")
        printlnFlush(string("Number of parameters: ",sum(map(length,w))," (",Int(floor(sum(map(length,w))/1e6))," million)"))
        printlnFlush("Loading data...")
        @time xtrain,ytrain = loadData_Robot(model,"train")
        @time xtest,ytest = loadData_Robot(model,"test")
        printlnFlush("Training...")
        grdnt = gradloss(loss)
        trainbatchCount = size(ytrain,4)
        testbatchCount = size(ytest,4)
        maxacc = 0
        for epoch=1:epochs
            trainloss = testloss = trainncorrect = testncorrect = trainncount = testncount = 0
            #Train
            h = zeros(h)
            c = zeros(c)
            for t=1:trainbatchCount
                printUpdate(@sprintf("Time: %5d/%-5d",t,trainbatchCount))
                g,l = grdnt(w,h,c,xtrain,ytrain,t,D)
                trainloss += l
                #Adam
                for i in 1:length(w)
                    m[i] = beta1*m[i] + (1-beta1)*g[i]
                    v[i] = beta2*v[i] + (1-beta2)*(g[i].*g[i])
                    mhat = m[i] ./ (1-beta1^adamT)
                    vhat = v[i] ./ (1-beta2^adamT)
                    w[i] = w[i] - (lr ./ (sqrt(vhat)+eps)).*mhat
                end
                adamT += 1
            end
            #Testing accuracy
            if epoch%testInterval == 0 || epoch == 1
                h = zeros(h)
                c = zeros(c)
                for t=1:testbatchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,testbatchCount))
                    ypred = Array(predict(w,h,c,KnetArray(xtest[:,:,:,t:t])))
                    testloss += pixelwiseSoftloss(ypred,ytest[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,ytest[:,:,:,t:t])
                    testncorrect += testAcc[1]
                    testncount += testAcc[2]
                end
            end
            testaccuracy = testncorrect/testncount*100
            if isnan(testaccuracy)
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f",epoch,trainloss/trainbatchCount))
            else
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f\tTest Loss: %0.10f\tTest Accuracy: %0.4f%%",
                    epoch,trainloss/trainbatchCount,testloss/testbatchCount,testaccuracy))
            end
            if testaccuracy > maxacc
                maxacc = testaccuracy
                saveNetwork(w,"Robot_LSTM_$(model)")
            end
            if epoch == epochs
                printlnFlush(@sprintf("Maximum Test Accuracy: %0.4f%%",maxacc))
            end
        end
    end
end

#   Conv
function run_Robot_Conv_D_train()
    train_Robot_Conv("D",Robot_Conv_D,Robot_Conv_D_loss,Robot_Conv_D_w(); lr=1e-5,epochs=100,batchsize=16,testInterval=1)
end
function run_Robot_Conv_RGB_train()
    train_Robot_Conv("RGB",Robot_Conv_RGB,Robot_Conv_RGB_loss,Robot_Conv_RGB_w(); lr=1e-5,epochs=100,batchsize=16,testInterval=1)
end
function run_Robot_Conv_RGBD_train()
    train_Robot_Conv("RGBD",Robot_Conv_RGBD,Robot_Conv_RGBD_loss,Robot_Conv_RGBD_w(); lr=1e-5,epochs=50,batchsize=8,testInterval=1)
end

#   LSTM
function run_Robot_LSTM_D_train()
    train_Robot_LSTM("D",Robot_LSTM_D,Robot_LSTM_D_loss,Robot_LSTM_D_w(); lr=1e-5,epochs=50,testInterval=1)
end
function run_Robot_LSTM_RGB_train()
    train_Robot_LSTM("RGB",Robot_LSTM_RGB,Robot_LSTM_RGB_loss,Robot_LSTM_RGB_w(); lr=1e-5,epochs=50,testInterval=1)
end
function run_Robot_LSTM_RGBD_train()
    train_Robot_LSTM("RGBD",Robot_LSTM_RGBD,Robot_LSTM_RGBD_loss,Robot_LSTM_RGBD_w(); lr=1e-5,epochs=50,testInterval=1)
end

#   Conv
export run_Robot_Conv_D_train
export run_Robot_Conv_RGB_train
export run_Robot_Conv_RGBD_train
#   LSTM
export run_Robot_LSTM_D_train
export run_Robot_LSTM_RGB_train
export run_Robot_LSTM_RGBD_train

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
module Day24

export day24

using CombinedParsers
using CombinedParsers.Regexp: whitespace_maybe

abstract type Direction end

struct SE <: Direction end
struct SW <: Direction end
struct NE <: Direction end
struct NW <: Direction end
struct E <: Direction end
struct W <: Direction end

abstract type Color end

struct White <: Color end
struct Black <: Color end

function parseinput(filename)
    mapdir = Dict{String,Type{<:Direction}}(["se" => SE, "sw" => SW, "ne" => NE, "nw" => NW, "e" => E, "w" => W])
    p = join(Repeat("se" | "sw" | "ne" | "nw" | "e" | "w"),"")

    input = Vector{Vector{Type{<:Direction}}}()
    for line in readlines(filename)
        v = map(s -> mapdir[s], line |> p)
        push!(input, v)
    end
    input
end

function step(d::Type{<: Direction})
    if d == SE
        (0,-1)
    elseif d == SW
        (-1,-1)
    elseif d == NE
        (1,1)
    elseif d == NW
        (0,1)
    elseif d == E
        (1,0)
    else
        (-1,0)
    end
end

import Base.:+

function Base.:+(x::Tuple{<:Integer, <:Integer}, y::Tuple{<:Integer, <:Integer})
    (x[1] + y[1], x[2] + y[2])
end

function findtile(path)
    pos = (0,0)
    for direction in path
        pos += step(direction)
    end
    pos
end

flip(c::Type{<:Color}) = c == Black ? White : Black

function initialize(input)
    floor = Dict{Tuple{Int,Int}, Type{<:Color}}()
    for path in input
        tile = findtile(path)
        floor[tile] = flip(get(floor, tile, White))
    end
    floor
end

function day24(::Val{:a}, filename)
    input = parseinput(filename)
    floor = initialize(input)
    count(c -> c == Black, values(floor))
end

function neighbors(pos)
    (x,y) = pos
    Set{Tuple{Int,Int}}([(x,y+1), (x,y-1), (x-1,y), (x+1,y), (x+1,y+1), (x-1, y-1)])
end

function change(floor::Dict{Tuple{Int,Int}, Type{<:Color}})
    new_floor = copy(floor)
    cells = Set{Tuple{Int,Int}}(keys(floor))
    layout = copy(cells)
    for pos in cells
        union!(layout, neighbors(pos))
    end
    for pos in layout
        blacks = 0
        for nei in neighbors(pos)
            blacks += get(floor, nei, White) == Black ? 1 : 0
        end
        if get!(floor,pos,White) == Black && (blacks == 0 || blacks > 2)
            new_floor[pos] = White
        end
        if get!(floor,pos,White) == White && blacks == 2
            new_floor[pos] = Black
        end
    end
    new_floor
end

function day24(::Val{:b}, filename)
    input = parseinput(filename)
    floor = initialize(input)
    for day in 1:100
        floor = change(floor)
    end
    count(c -> c == Black, values(floor))
end

end


text:
R"""
rm(list=ls())
Start.Time <- Sys.time()

#load libraries
list.of.packages <- c("data.table","bit64","optimx","lme4",'languageR','lmerTest','glmnet','lsmeans','car','RDS','DMwR','stringr','nloptr','ggplot2','reshape2','plyr','viridis','gmodels')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
lapply(list.of.packages,require,character.only=TRUE)

#######################################
#------------- Load Data -------------#
#######################################

# set working directory
setwd("/media/u01/analytics/scoring/CDW5_792")

# set path of data and headers location
input <- c("/media/u01/analytics/scoring/CDW5_792/")

# set output file path
output <- c("/media/u01/analytics/scoring/CDW5_792/Output/")

# set name of initial data file
dataname <- "csv_final_cdw5_792"

# set name of headers file
headers_data <- "Headers_cdw5.csv"

# read initial data and attach headers
finaldata <- fread(paste0(input,dataname),colClasses = "as.character",header=F)
datanames <- fread(paste0(input,headers_data),header=F)
setnames(finaldata,c(datanames$V1))

#check if response and mandatory variables exist in the data
#Put also the sub campaigns as mandatory
all_mandatory_vars <- c("Buyer_Pos_P1","Buyer_Pre_P1","Trps_POS_P1","Trps_PRE_P1","Dol_per_Trip_POS_P1","Dol_per_Trip_PRE_P1","Nonbuyer_Pre_P1","hh_age","estimated_hh_income","number_of_children_in_living_Un","person_1_gender","experian_id")
print(all_mandatory_vars %in% names(finaldata))

# set exposed flag variable
exposed_flag_var <- "exposed_flag_new"

# Put the sub campaigns and the demographic variables that will be used as random effects in the models
random_demos <- c("estimated_hh_income","hh_age","number_of_children_in_living_Un","person_1_gender")
random_campaigns <- c("Publisher_Fct1","Targeting_Fct1")

# set extra exposed flag variable to be excluded (in case you have more than one in dataset) and publisher variables that are not used in the mixed models
dropvars <- c("exposed_flag")

# set True if P2 category products are competitor or False otherwise (need for sign check in the model)
P2_Competitor <- TRUE #or FALSE

# pvalue threshold
pvalue_lvl <- 0.20

# Create folders for outputs
dir.create(paste0(output,"Descriptive Statistics for Publishers and Demographics"))
dir.create(paste0(output,"Final_Modelling_Results"))
dir.create(paste0(output,"Initial_Modelling_Results"))
dir.create(paste0(output,"QCs"))
dir.create(paste0(output,"saved_data"))
dir.create(paste0(output,"Scoring"))
final_modelling_output <- paste0(output,"Final_Modelling_Results/")
initial_modelling_output <- paste0(output,"Initial_Modelling_Results/")
scoring_output <- paste0(output,"Scoring/")
dir.create(paste0(scoring_output,"Occasion/"))
dir.create(paste0(scoring_output,"DollarsOccasion/"))
dir.create(paste0(scoring_output,"Penetration/"))
dir.create(paste0(final_modelling_output,"Occasion/"))
dir.create(paste0(final_modelling_output,"DollarsOccasion/"))
dir.create(paste0(final_modelling_output,"Penetration/"))
dir.create(paste0(initial_modelling_output,"Occasion/"))
dir.create(paste0(initial_modelling_output,"DollarsOccasion/"))
dir.create(paste0(initial_modelling_output,"Penetration/"))

scoring_Occ_output <- paste0(scoring_output,"Occasion/")
scoring_DolOcc_output <- paste0(scoring_output,"DollarsOccasion/")
scoring_Pen_output <- paste0(scoring_output,"Penetration/")
final_Occ_modelling_output <- paste0(final_modelling_output,"Occasion/")
final_DollOcc_modelling_output <- paste0(final_modelling_output,"DollarsOccasion/")
final_Pen_modelling_output <- paste0(final_modelling_output,"Penetration/")
initial_Occ_modelling_output <- paste0(initial_modelling_output,"Occasion/")
initial_DollOcc_modelling_output <- paste0(initial_modelling_output,"DollarsOccasion/")
initial_Pen_modelling_output <- paste0(initial_modelling_output,"Penetration/")

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

# create group, panid columns
finaldata[,`:=`(group=eval(parse(text=exposed_flag_var)),panid=experian_id)]

# remove not needed variables
if(length(dropvars)>0){
  for(i in noquote(dropvars)){
    finaldata <- finaldata[,noquote(i):=NULL]
  }
}

# mandatory variables for scoring
scoring_vars <- c("Prd_1_Net_Pr_PRE","Prd_1_Net_Pr_POS","Buyer_Pos_P0","Buyer_Pre_P0")

# Set the ProScore variable in the dataset
ProScore <- grep("MODEL",colnames(finaldata),value = TRUE)

# aggregate # of children for 4+
finaldata[number_of_children_in_living_Un>=4,number_of_children_in_living_Un:="4+" ]

# create factor variables
allrandoms <- c(random_demos,random_campaigns)
finaldata  <- as.data.frame(finaldata)
for(i in noquote(allrandoms)){
  finaldata[i] <- lapply(finaldata[i],function(x) as.factor(x))
}
finaldata  <- as.data.table(finaldata)

#     for(i in random_campaigns)
#     {
#       finaldata <- as.data.frame(finaldata)
#       finaldata[,i] <- as.factor(finaldata[,i])
#
#       # re-group random intercept levels
#       finaldata[,i] <- str_sub(finaldata[,i], start=1, end=2)
#       finaldata[,i] <- gsub("_","",finaldata[,i])
#       finaldata[,i] <- as.factor(finaldata[,i])
#
#       # create a csv file with the counts of each level of the publisher variable
#       #write.csv(table(finaldata$group,finaldata[,eval(parse(text=random_campaigns[i]))]),file="ExposureFreq5w.csv")
#
#       levelstokeep <- c("14","15","41","3","8","2","6","10","24")
#       all_levels <- unique(finaldata[,i])
#
#       other_levels <- setdiff(all_levels,c(levelstokeep,"\\N"))
#
#       finaldata <- as.data.table(finaldata)
#       for(j in other_levels){
#         finaldata <- finaldata[i==j,i:="Other"]
#       }
#
#     }

# GMT Check

# set variables as.numeric and replace NA's with zero
finaldata <- as.data.frame(finaldata)
indx <- sapply(finaldata, is.character)
options(warn=-1)
finaldata[indx] <- lapply(finaldata[indx],function(x) as.numeric(as.character(x)))
finaldata[,exposed_flag_var] <- NULL
options(warn=0)

# replace NAs with zero for numeric variables
finaldata <- as.data.frame(finaldata)
finaldata[sapply(finaldata, is.numeric)][is.na(finaldata[sapply(finaldata, is.numeric)])] <- 0
finaldata <- as.data.table(finaldata)

# set as zero the //N value of the group variable
finaldata <- finaldata[is.na(group),group:="0"]
finaldata <- finaldata[group=="\\N",group:="0"]
finaldata$group <- factor(finaldata$group)
finaldata$group <- relevel(finaldata$group,"0")

########## QCs ##########

# Create a summary of the dataset
t1 <- as.data.frame(finaldata[,sapply(finaldata,is.numeric)])
finaldata_num <- as.data.frame(finaldata)
finaldata_num <- finaldata_num[,c(which(t1==T))]
finaldata_num <- as.data.table(finaldata_num)
finaldata_num <- finaldata_num[,c(3:(which(names(finaldata_num)=="state")-1),((which(names(finaldata_num)=="mail_responder"))+1):length(names(finaldata_num))),with=F]
finaldata_num <- as.data.frame(finaldata_num)

summary_num <- sapply(finaldata_num[,!names(finaldata_num) %in% c("exposed_flag","experian_id","panid")],function(x) c(min(x,na.rm=T),quantile(x,c(0.25),na.rm=T),mean(x,na.rm=T),median(x,na.rm=T), quantile(x,c(0.75),na.rm=T),max(x,na.rm=T),sum(x,na.rm=T)))
rownames(summary_num) <- c("min","Q1","mean","median","Q3","max","sum")

na_count <-sapply(finaldata_num[,!names(finaldata_num) %in% c("exposed_flag","experian_id","panid")], function(y) sum(length(which(is.na(y)))))
NAs  <- as.vector(na_count)
N <- nrow(finaldata)-na_count

#write csv
write.csv(t(rbind(N,NAs,summary_num)),paste0(output,"QCs/","Summary.csv"),na='')

# Create Crosstables for campaign and pre campaign period
sink(paste0(output,"QCs/","PreCampaign_Summary.txt"))
CrossTable(finaldata$group,finaldata$Buyer_Pre_P1, expected = FALSE)
sink()
sink(paste0(output,"QCs/","Campaign_Summary.txt"))
CrossTable(finaldata$group,finaldata$Buyer_Pos_P1,expected = FALSE)
sink()

# Create check list for factors and Buyers
buy.pre <-  unique(finaldata$Buyer_Pre_P1)
buy.pos <- unique(finaldata$Buyer_Pos_P1)
non.pre <- unique(finaldata$Nonbuyer_Pre_P1)
gr <- unique(finaldata$group)
mos <-  unique(finaldata$Mosaic)
ban <- unique(finaldata$banner)
exp.uniq <- length(unique(finaldata$experian_id))

# check freq for demographic variables
sink(paste0(output,"QCs/","Demographics_Summary.txt"))
xtabs(~estimated_hh_income,finaldata)
cat("\n")
xtabs(~person_1_gender,finaldata)
cat("\n")
xtabs(~number_of_children_in_living_Un,finaldata)
cat("\n")
xtabs(~hh_age,finaldata)

cat("\n")
ftable(xtabs(~group+estimated_hh_income,finaldata))
cat("\n")
ftable(xtabs(~group+person_1_gender,finaldata))
cat("\n")
ftable(xtabs(~group+number_of_children_in_living_Un,finaldata))
cat("\n")
ftable(xtabs(~group+hh_age,finaldata))
sink()

finaldata_buyers <- finaldata[Buyer_Pos_P1=="1",]
sink(paste0(output,"QCs/","Demographics_Summary_Buyers_Campaign_Period.txt"))
cat("\n")
xtabs(~estimated_hh_income,finaldata_buyers)
cat("\n")
xtabs(~person_1_gender,finaldata_buyers)
cat("\n")
xtabs(~number_of_children_in_living_Un,finaldata_buyers)
cat("\n")
xtabs(~hh_age,finaldata_buyers)

cat("\n")
ftable(xtabs(~group+estimated_hh_income,finaldata_buyers))
cat("\n")
ftable(xtabs(~group+person_1_gender,finaldata_buyers))
cat("\n")
ftable(xtabs(~group+number_of_children_in_living_Un,finaldata_buyers))
cat("\n")
ftable(xtabs(~group+hh_age,finaldata_buyers))
sink()

sink(paste0(output,"QCs/","Buyers Data Metrics.txt"))
list(
  Mean_Dollars_Spent_per_Trip_POST_P1_Non_Exposed=mean(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  Mean_Dollars_Spent_per_Trip_PRE_P1_Non_Exposed=mean(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  Mean_Dollars_Spent_per_Trip_POST_P1_Exposed=mean(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  Mean_Dollars_Spent_per_Trip_PRE_P1_Exposed=mean(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  Mean_Trips_POST_P1_Non_Exposed=mean(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  Mean_Trips_PRE_P1_Non_Exposed=mean(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),
  Mean_Trips_POST_P1_Exposed=mean(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  Mean_Trips_PRE_P1_Exposed=mean(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),

  SD_Dol_per_Trip_POST_P1_Non_Exposed=sd(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  SD_Dol_per_Trip_PRE_P1_Non_Exposed=sd(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  SD_Dol_per_Trip_POST_P1_Exposed=sd(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  SD_Dol_per_Trip_PRE_P1_Exposed=sd(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  SD_Trips_POST_P1_Non_Exposed=sd(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  SD_Trips_PRE_P1_Non_Exposed=sd(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),
  SD_Trips_POST_P1_Exposed=sd(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  SD_Trips_PRE_P1_Exposed=sd(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),

  Min_Dol_per_Trip_POST_P1_Non_Exposed=min(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  Min_Dol_per_Trip_PRE_P1_Non_Exposed=min(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  Min_Dol_per_Trip_POST_P1_Exposed=min(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  Min_Dol_per_Trip_PRE_P1_Exposed=min(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  Min_Trips_POST_P1_Non_Exposed=min(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  Min_Trips_PRE_P1_Non_Exposed=min(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),
  Min_Trips_POST_P1_Exposed=min(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  Min_Trips_PRE_P1_Exposed=min(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),

  Max_Dol_per_Trip_POST_P1_Non_Exposed=max(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  Max_Dol_per_Trip_PRE_P1_Non_Exposed=max(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  Max_Dol_per_Trip_POST_P1_Exposed=max(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_POS_P1),
  Max_Dol_per_Trip_PRE_P1_Exposed=max(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Dol_per_Trip_PRE_P1),
  Max_Trips_POST_P1_Non_Exposed=max(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  Max_Trips_PRE_P1_Non_Exposed=max(finaldata[group=="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1),
  Max_Trips_POST_P1_Exposed=max(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_POS_P1),
  Max_Trips_PRE_P1_Exposed=max(finaldata[group!="0" & Buyer_Pos_P1=="1",]$Trps_PRE_P1)
)
sink()

cat("Please check summary.csv and proceed with the modeling process")

# remove HHs with no gender info
finaldata <- finaldata[!person_1_gender=="U",]

# aggregate U and L levels of hh income
finaldata <- finaldata[estimated_hh_income=="U",estimated_hh_income:="L"]

# check and drop exposed HHs with no publisher info or non-exposed HHs with publisher info
panid_out_all <- c()
if(length(random_campaigns)>0){
  for(k in random_campaigns){
    finaldata <- finaldata[eval(parse(text=k))=="\\N"|eval(parse(text=k))=="NULL"|eval(parse(text=k))=="0"|eval(parse(text=k))=="NONE",eval(parse(text=k)):="none"]
    cat(paste0(nrow(finaldata[(group=="0" & eval(parse(text=k))!="none")])," non exposed HHs with publisher info "))
    panid_out <- finaldata[(group=="0" & eval(parse(text=k))!="none"),panid]
    finaldata <- finaldata[!(group=="0" & eval(parse(text=k))!="none"),]
    cat('\n')
    # check and remove HHs that were exposed and we had no publisher info
    cat(paste0(nrow(finaldata[(group!="0" & eval(parse(text=k))=="none")])," exposed HHs with no publisher info "))
    finaldata <- finaldata[!(group!="0" & eval(parse(text=k))=="none")]
    panid_out2 <- finaldata[(group!="0" & eval(parse(text=k))=="none"),panid]
    panid_out_all <- unique(c(panid_out_all,panid_out,panid_out2))
  }
}

# apply matching rules for control group vs exposed (default is TRUE)
Match_control_exposed_groups <- TRUE  #FALSE or TRUE

if(Match_control_exposed_groups==T & length(ProScore)==0){
  finaldata <- as.data.table(finaldata)
  data_exposed <- finaldata[group!="0" ,]
  data_nonexposed1 <- finaldata[group=="0" & Buyer_Pre_P1=="1",]
  data_nonexposed0 <- finaldata[group=="0" & Buyer_Pre_P1=="0",]

  if(nrow(finaldata[group=="0",])>2000000){
    new_data_dim <- 0.3*nrow(finaldata[group=="0",])
  } else {
    if(nrow(finaldata[group=="0",])>1000000){
      new_data_dim <- 0.4*nrow(finaldata[group=="0",])
    } else {
      if(nrow(finaldata[group=="0",])>750000){
        new_data_dim <- 0.6*nrow(finaldata[group=="0",])
      } else {
        new_data_dim <- 0.7*nrow(finaldata[group=="0",])
      }
    }
  }

  dim_sample0 <- round(((new_data_dim-nrow(data_exposed))/(1+(nrow(data_exposed[Buyer_Pre_P1=="1",])/(nrow(data_exposed[Buyer_Pre_P1=="0",]))))))
  dim_sample1 <- new_data_dim-nrow(data_exposed)-dim_sample0

  new_data_nonexposed_11 <- as.data.frame(data_nonexposed1[sample(nrow(data_nonexposed1), dim_sample1),])
  new_data_nonexposed_01 <- as.data.frame(data_nonexposed0[sample(nrow(data_nonexposed0), dim_sample0),])
  data_exposed <- as.data.frame(data_exposed)
  finaldata_sample <- rbind(data_exposed,new_data_nonexposed_11,new_data_nonexposed_01)

  finaldata <- as.data.table(finaldata_sample)
}

if(Match_control_exposed_groups==T & length(ProScore)>0){
  exposed_proportion <- prop.table(xtabs(~eval(parse(text=ProScore)),finaldata[group!="0",]))
  proscore_levels <- names(exposed_proportion)

  if(nrow(finaldata[group=="0",])>2000000){
    new_data_dim <- 0.3*nrow(finaldata[group=="0",])
  } else {
    if(nrow(finaldata[group=="0",])>1000000){
      new_data_dim <- 0.4*nrow(finaldata[group=="0",])
    } else {
      if(nrow(finaldata[group=="0",])>800000){
        new_data_dim <- 0.6*nrow(finaldata[group=="0",])
      } else {
        new_data_dim <- 0.7*nrow(finaldata[group=="0",])
      }
    }
  }

  sample_control_data <- c()

  for (i in proscore_levels){
    sample_i_dim <- round(new_data_dim*exposed_proportion[names(exposed_proportion)==i])
    temp_data <- finaldata[eval(parse(text=ProScore))==i & group=="0",]
    data_i <- temp_data[sample(nrow(temp_data), sample_i_dim),]
    sample_control_data <- rbind(sample_control_data,data_i)
    #    cat(i)
    #    cat("\n")
    #    cat(length(sample_control_data$panid))
    #    cat("\n")
  }

  sample_data <- rbind(sample_control_data,finaldata[group!="0",])

  sample_data <- as.data.table(sample_data)
  data_exposed <- as.data.table(sample_data[group!="0",])

  data_nonexposed1 <- sample_data[group=="0" & Buyer_Pre_P1=="1",]
  data_nonexposed0 <- sample_data[group=="0" & Buyer_Pre_P1=="0",]

  if(nrow(data_nonexposed1)/nrow(data_nonexposed0) > nrow(data_exposed[Buyer_Pre_P1=="1",])/nrow(data_exposed[Buyer_Pre_P1=="0",])){
    new_data_nonexposed_11 <- as.data.frame(data_nonexposed1[sample(nrow(data_nonexposed1), (nrow(data_exposed[Buyer_Pre_P1=="1",])/nrow(data_exposed[Buyer_Pre_P1=="0",]))*nrow(data_nonexposed0)),])
    data_exposed <- as.data.frame(data_exposed)
    finaldata_sample <- rbind(data_exposed,new_data_nonexposed_11,data_nonexposed0)
  } else {
    new_data_nonexposed_01 <- as.data.frame(data_nonexposed0[sample(nrow(data_nonexposed0), (nrow(data_exposed[Buyer_Pre_P1=="0",])/nrow(data_exposed[Buyer_Pre_P1=="1",]))*nrow(data_nonexposed1)),])
    data_exposed <- as.data.frame(data_exposed)
    finaldata_sample <- rbind(data_exposed,new_data_nonexposed_01,data_nonexposed1)
  }
  length(finaldata_sample$panid)

  finaldata <- as.data.table(finaldata_sample)
}

# check proportions of the sample dataset
if(Match_control_exposed_groups==T ){
  finaldata_sample <- as.data.table(finaldata_sample)
  sink(paste0(output,"Matching_data_summary.txt"))
  if(length(ProScore)>0){
    cat("\n")
    xt <- xtabs(~eval(parse(text=ProScore)),finaldata_sample[group!="0",])
    names(dimnames(xt)) <- c("Proportions of ProScore in the exposed group")
    print(prop.table(xt))
    cat("\n\n")
    xt <- xtabs(~eval(parse(text=ProScore)),finaldata_sample[group=="0",])
    names(dimnames(xt)) <- c("Proportions of ProScore in the control group")
    print(prop.table(xt))
  }
  cat("\n\n")
  xt <- xtabs(~Buyer_Pre_P1,finaldata_sample[group!="0",])
  names(dimnames(xt)) <- c("Percentages of target product buyers during the Pre-Campaign Period in the exposed group")
  print(prop.table(xt))
  cat("\n\n")
  xt <- xtabs(~Buyer_Pre_P1,finaldata_sample[group=="0",])
  names(dimnames(xt)) <- c("Percentages of target product buyers during the Pre-Campaign Period in the control group")
  print(prop.table(xt))
  cat("\n\n")
  cat('Average_Trips_PRE_P1_Non_Exposed_matched_data')
  cat("\n")
  cat(mean(finaldata_sample[Buyer_Pre_P1=="1" & group=="0",]$Trps_PRE_P1))
  cat("\n\n")
  cat('Average_Trips_PRE_P1__Exposed_matched_data')
  cat("\n")
  cat(mean(finaldata_sample[Buyer_Pre_P1=="1" & group!="0",]$Trps_PRE_P1))
  cat("\n\n")
  cat('Average_Spent_PRE_P1_Non_Exposed_matched_data')
  cat("\n")
  cat(mean(finaldata_sample[Buyer_Pre_P1=="1" & group=="0",]$Dol_per_Trip_PRE_P1))
  cat("\n\n")
  cat('Average_Spent_PRE_P1_Exposed_matched_data')
  cat("\n")
  cat(mean(finaldata_sample[Buyer_Pre_P1=="1" & group!="0",]$Dol_per_Trip_PRE_P1))
  sink()
}

# remove demographic variables
finaldata <- droplevels(finaldata)
finaldata$panid <- as.character(finaldata$panid)
initial_data <- finaldata[,c(4:(which(names(finaldata)=="state")-1),which(names(finaldata)=="person_1_gender"),which(names(finaldata)=="number_of_children_in_living_Un"),(which(names(finaldata)=="Mosaic")+1):length(names(finaldata))),with=F]
names(initial_data)
options(warn=-1)

#GMT
xVarsDemos <- setdiff(names(finaldata),names(initial_data))

######rm(finaldata,data_nonexposed1,temp_data,data_i,sample_control_data,sample_data,data_exposed,data_nonexposed0,finaldata_sample,finaldata_num,new_data_nonexposed_11,new_data_nonexposed_01)
options(warn=0)

# get number of products in the data
num_products <- length(grep("Buyer_Pre_P",colnames(initial_data)))-1

#exclude POST variables
indexPOS <- grep("POS",colnames(initial_data))
indexPos <- grep("Pos",colnames(initial_data))
indexPre1 <- grep("Buyer_Pre_",colnames(initial_data))
indexPre2 <- grep("Nonbuyer_Pre_",colnames(initial_data))
# check pre vars to exclude (only buyer_pre and nonbuyer_pre should be excluded)
(colnames(initial_data)[c(indexPre1,indexPre2)])
initial_vars_to_exclude <- colnames(initial_data)[c(indexPOS,indexPos,indexPre1,indexPre2)]
indexPREPOS <- grep("PRE_POS",colnames(initial_data))
indexPrePos <- grep("Pre_Pos",colnames(initial_data))
PREPOS_vars_to_include <- colnames(initial_data)[c(indexPREPOS,indexPrePos)]
ALL_vars_to_exclude <- setdiff(initial_vars_to_exclude,c(PREPOS_vars_to_include,all_mandatory_vars,scoring_vars))
initial_data <- as.data.frame(initial_data)
initial_data <- initial_data[ , !names(initial_data) %in% c(ALL_vars_to_exclude)]

#exclude category variables(exclude P0's)
indexP0 <- grep("0",colnames(initial_data))
P0_vars_to_exclude <- colnames(initial_data)[indexP0]
indexMODEL <- grep("MODEL",colnames(initial_data))
ProScore_vars_to_include <- colnames(initial_data)[indexMODEL]
P0_vars_to_exclude <- setdiff(P0_vars_to_exclude,c(ProScore_vars_to_include,all_mandatory_vars,scoring_vars))
initial_data <- initial_data[ , !names(initial_data) %in% c(P0_vars_to_exclude)]

#exclude new variables that are only for reporting
indexPERC <- grep("Perc_",colnames(initial_data))
indexPriPerUn <- grep("Pr_per_",colnames(initial_data))
vars_to_exclude <- colnames(initial_data)[c(indexPERC,indexPriPerUn)]
initial_data <- initial_data[ , !names(initial_data) %in% c(vars_to_exclude)]

# segments for outliers detection
initial_data <- as.data.table(initial_data)
data_NB_NE_B <- initial_data[Buyer_Pre_P1=="0" & group=="0"  & Buyer_Pos_P1=="1",]
#data_B_E_NB <- initial_data[Buyer_Pre_P1=="1" & group!="0" & Buyer_Pos_P1=="0",]
pen_reduction <- initial_data[Buyer_Pre_P1=="1" & group=="0" & Buyer_Pos_P1=="0",]
occ_reduction <- initial_data[group=="0" & Buyer_Pos_P1=="1" & Trps_POS_P1<Trps_PRE_P1,]
dol_reduction <- initial_data[group=="0" & Buyer_Pos_P1=="1" & Dol_per_Trip_POS_P1<Dol_per_Trip_PRE_P1,]

# get variables that need to have negative sign
negative_index <- c()
for(i in 3:num_products){
  indexCompetitorVars <- grep(paste0("P",i),colnames(initial_data))
  indexCompetitorVars2 <- grep(paste0(i,"_"),colnames(initial_data))
  indexGAS <- grep("AVG_PRICE",colnames(initial_data))
  negative_index <- c(negative_index,indexCompetitorVars,indexCompetitorVars2,indexGAS)
}
negativevars <- unique(colnames(initial_data)[negative_index])

if(P2_Competitor==T){
  indexNOCompetitorVars1 <- grep("P2",colnames(initial_data))
  indexNOCompetitorVars2 <- grep(paste0(2,"_"),colnames(initial_data))
  indexNOCompetitorVars <- c(indexNOCompetitorVars1,indexNOCompetitorVars2)
  NOCompetitorVars <- colnames(initial_data)[indexNOCompetitorVars]
  negativevars <- unique(c(negativevars,NOCompetitorVars))
}

# get variables that need to have positive sign
positive_index1 <- grep("P1",colnames(initial_data))
positive_index2 <- grep("1_",colnames(initial_data))
positive_index3 <- grep("MODEL",colnames(initial_data))
positive_index <- c(positive_index1,positive_index2,positive_index3)
positivevars <- colnames(initial_data)[positive_index]

# # kruskal wallis test for pair comparison of publisher's levels keep aggregated levels based on the spent aggregation
# kruskal_var <- 'Dol_per_Trip_POS_P1'
# buyers_data <- initial_data[Buyer_Pos_P1=='1',]
# kruskal_randoms <- c()
# for(i in random_campaigns){
#   kt <- kruskal.test(eval(parse(text=kruskal_var))~eval(parse(text=i)),data=buyers_data)
#   cat(i);print(kt)
#   if(kt$p.value<0.10){
#     kruskal_randoms <- append(kruskal_randoms,i)
#   } else {
#     initial_data[,eval(parse(text=i)):=NULL]
#   }
# }
#
# warning('In case of aggregation based on Kruskal-Wallis test, make sure that the other level of each sub-campaign is named as Other_XXX')
# for(i in kruskal_randoms){
#   lev <- setdiff(levels(buyers_data[,eval(parse(text=i))]),'none')
#   if(length(lev)>2){
#     aggr <- c()
#     for(j in 1:(length(lev)-1)){
#       for( k in (j+1):(length(lev))){
#         df <- subset(buyers_data, eval(parse(text=i))  %in% c(lev[j],lev[k]))
#         l1 <- list(c(lev[j],lev[k]),P_Value=try(round(kruskal.test(df[,eval(parse(text=kruskal_var))]~df[,eval(parse(text=i))])$p.value,5),silent=TRUE))
#         if(l1[2]=="NaN") next
#         vectorElements <- unlist(l1)
#         listnames <- names(vectorElements) <- c(paste0(i,':'),paste0(i,':'),'P_value:')
#         listout <- paste(listnames,vectorElements)
#         write.table(listout,file=paste0(output,'Kruskal.csv'), row.names = FALSE, col.names = FALSE, quote = FALSE,append=TRUE)
#         if(l1[['P_Value']]<0.05){
#           aggr <- append(aggr,c(sapply(l1[1],`[`,1),sapply(l1[1],`[`,2)))
#           aggr <- setdiff(aggr,c('none',paste0("Other_",i)))
#         }
#       }
#     }
#     if(length(unique(aggr))>0){
#       initial_data[!eval(parse(text=i)) %in% c(aggr,'none'),eval(parse(text=i)):=paste0("Other_",i)]
#     } else {
#       initial_data[,eval(parse(text=i)):=NULL]
#     }
#   }
# }

initial_data <- droplevels(initial_data)
saveRDS(initial_data,"initial_data.RDS")
random_campaigns <- intersect(random_campaigns,names(initial_data))

if(length(random_campaigns)>0){
  random_campaigns_factors <- paste0("+(0+group|",random_campaigns,")")
} else {
  random_campaigns_factors <- NULL
}
random_demos_factors <- paste0( "+(0+group|",random_demos,")" )

# Create Descriptives for Campaign breaks
if(length(random_campaigns)>0){
  output_checks <- paste0(output,"Descriptive Statistics for Publishers and Demographics/")
  randoms <- c(random_demos,random_campaigns)
  for(i in randoms){
    write.csv(table(initial_data[,eval(parse(text=i))]),file=paste0(output_checks,i,'.csv'),row.names=F)
    write.csv(table(initial_data$group,initial_data[,eval(parse(text=i))]),file=paste0(output_checks,i,'_group.csv'))

    dolocc_data <- initial_data[Buyer_Pos_P1=="1",]

    write.csv(table(dolocc_data[,eval(parse(text=i))]),file=paste0(output_checks,i,'_buyers.csv'),row.names=F)
    write.csv(table(dolocc_data$group,dolocc_data[,eval(parse(text=i))]),file=paste0(output_checks,i,'_group_buyers.csv'))
  }

  # check exposed and non exposed HHs counts for target brand Buyers in the campaign period
  write.csv(table(initial_data$group,initial_data$Buyer_Pos_P1),file=paste0(output_checks,'buyers_pos_group.csv'))

  # check exposed and non exposed HHs counts for target brand Buyers in the pre-campaign period
  write.csv(table(initial_data$group,initial_data$Buyer_Pre_P1),file=paste0(output_checks,'buyers_pre_group.csv'))

  # calculate correlations between levels of the random intercept variables
  finaldata_cor <- as.data.frame(initial_data)

  for(i in randoms){
    for(level in unique(eval(parse(text=paste0("finaldata_cor$",i))))){
      finaldata_cor[paste(i, level, sep = "_")] <- ifelse(eval(parse(text=paste0("finaldata_cor$",i))) == level, 1, 0)
    }
  }

  cor.data <- finaldata_cor[,c((which(names(finaldata_cor)=="panid")+1):length(finaldata_cor))]
  cor.matrix <- cor(cor.data)
  write.csv(cor.matrix,file=paste0(output_checks,"random_correlation_matrix.csv"))

  mel.cc <- melt(cor.matrix)
  heat_map <- ggplot(mel.cc,aes(Var1, Var2))+geom_tile(aes(fill=value))+ scale_fill_viridis(name="Correlation",option="magma") + theme(axis.text.x = element_text(angle=90, vjust=0.5))+ coord_equal()
  ggsave(filename=paste0(output_checks,"correlation_heat_map.pdf"),heat_map,width = 10, height = 10, units = "in")

  #Make descriptive plots for publishers, creatives etc.
  # create descriptive plots for the DolOcc breaks
  buyers_data <- initial_data[Buyer_Pos_P1=="1",]
  for(i in random_campaigns){
    buyers_data <- buyers_data[eval(parse(text=i))!="none",]
    buyers_data <- droplevels(buyers_data)
    buyers_data <- as.data.frame(buyers_data)
    res <- ldply(levels(buyers_data[,i]),function(x) round(mean(buyers_data[buyers_data[,i]==x,]$Dol_per_Trip_POS_P1),2))
    df <-cbind(res,levels(buyers_data[,i]))
    colnames(df) <- c('Average_Spent_during_the_campaign_period','Levels')

    p <- ggplot(df,aes(x=Levels,y=Average_Spent_during_the_campaign_period))+geom_bar(stat='identity')+ theme(axis.text.x = element_text(angle=90, vjust=0.5))+ylab('Average Spent for the target brand during the campaign period')+
      theme(legend.position="none")+theme(panel.background = element_rect(fill = 'white'))+
      geom_text(stat='identity',aes(label=Average_Spent_during_the_campaign_period,group=Levels), position = position_dodge(0.9), vjust = 0)
    ggsave(filename=paste0(output_checks,i,'_Mean_Spent.pdf'),p,height = 8,width = 15)

    p <- ggplot(buyers_data,aes_string(i,'Dol_per_Trip_POS_P1'))+geom_boxplot()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
      ylab('Spent per Trip for the target bramd during the campaign period')
    ggsave(filename=paste0(output_checks,i,'_Spent_BoxPlot.pdf'),p,height = 8,width = 15)
    buyers_data <- as.data.table(buyers_data)
  }

  # create descriptive plots for the Occ breaks
  options(warn=-1)
  a1 <- buyers_data[Trps_POS_P1>=5,Trps_POS_P1:="5"]
  options(warn=0)
  a1[,Trips_growth:=log(Trps_POS_P1/(Trps_PRE_P1+1))]
  a1$Trps_POS_P1 <- as.factor(a1$Trps_POS_P1)
  a1 <- a1[Trps_POS_P1=="5",Trps_POS_P1:="5+"]

  response_var <- c("Trps_POS_P1")
  for(i in random_campaigns){
    a1 <- as.data.table(a1)
    a1 <- a1[eval(parse(text=i))!="none",]
    a1 <- droplevels(a1)
    rr <- c(response_var,i)
    a1 <- as.data.frame(a1)
    df <- a1[,rr]
    mel <- melt(df,id.var=i)
    p <- ggplot(mel,aes_string(i))+geom_bar(aes(y=..count..,fill=as.factor(value)),position=position_dodge(0.9))+
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+guides(fill=guide_legend(title=response_var))+
      geom_text(stat='count',aes(label=..count..,group=value), position = position_dodge(0.9), vjust = 0)
    ggsave(filename=paste0(output_checks,i,'_Trips.pdf'),p,height = 8,width = 15)

    a2 <- as.data.table(initial_data[Buyer_Pos_P1=="1",])
    a2 <- a2[eval(parse(text=i))!="none",]
    a2 <- as.data.frame(a2)
    a2 <- droplevels(a2)
    res <- ldply(levels(a2[,i]),function(x) round(mean(a2[a2[,i]==x,]$Trps_POS_P1),2))
    df <-cbind(res,levels(a2[,i]))
    colnames(df) <- c('Average_Trips_during_the_campaign_period','Levels')
    p <- ggplot(df,aes(x=Levels,y=Average_Trips_during_the_campaign_period))+geom_bar(stat='identity')+ theme(axis.text.x = element_text(angle=90, vjust=0.5))+ylab('Average Trips for the target brand during the campaign period')+
      theme(legend.position="none")+theme(panel.background = element_rect(fill = 'white'))+
      geom_text(stat='identity',aes(label=Average_Trips_during_the_campaign_period,group=Levels), position = position_dodge(0.9), vjust = 0)
    ggsave(filename=paste0(output_checks,i,'_Mean_Trips.pdf'),p,height = 8,width = 15)
  }

  # create descriptive plots for the Pen breaks
  response_var <- c("Buyer_Pos_P1")
  for(i in random_campaigns){
    data <- as.data.table(copy(initial_data))
    data <- data[eval(parse(text=i))!="none",]
    data <- droplevels(data)
    data <- as.data.frame(data)
    rr <- c(response_var,i)
    df <- data[,rr]
    mel <- melt(df,id.var=i)
    p <- ggplot(mel,aes_string(x=i))+geom_bar(aes(y=..count..,fill=as.factor(value)),position=position_dodge(0.9))+
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+guides(fill=guide_legend(title=response_var))+
      geom_text(stat='count',aes(label=..count..,group=value), position = position_dodge(0.9), vjust = 0)
    ggsave(paste0(output_checks,i,'_Buyers_during_Campaign.pdf'),p,height = 8,width = 15)
  }
}

rm(dolocc_data,mel.cc,finaldata_cor,cor.data,cor.matrix,buyers_data,a1,a2,data)

#Create functions for modelling process and outputs
ggCaterpillar <- function(re,  QQ=FALSE, likeDotplot=FALSE) {
  require(ggplot2)
  f <- function(x) {
    pv   <- attr(x, "postVar")
    cols <- 1:(dim(pv)[1])
    se   <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
    ord  <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x))
    pDf  <- data.frame(y=unlist(x)[ord],
                       ci=1.28*se[ord], #80%CI:1.28, 90%CI:1.645, 95%CI: 1.96
                       nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)),
                       ID=factor(rep(rownames(x), ncol(x))[ord], levels=rownames(x)[ord]),
                       ind=gl(ncol(x), nrow(x), labels=names(x)))

    if(QQ) {  ## normal QQ-plot
      p <- ggplot(pDf, aes(nQQ, y))
      p <- p + facet_wrap(~ ind, scales="free")
      p <- p + xlab("Standard normal quantiles") + ylab("Random effect quantiles")
    } else {  ## caterpillar dotplot
      p <- ggplot(pDf, aes(ID, y)) #+ coord_flip()
      if(likeDotplot) {  ## imitate dotplot() -> same scales for random effects
        p <- p + facet_wrap(~ ind)
      } else {           ## different scales for random effects
        p <- p + facet_wrap(~ ind, scales="free_y")
      }
      p <- p + xlab("Levels") + ylab("Random effects")
    }

    p <- p + theme(legend.position="none")
    p <- p + geom_hline(yintercept=0)
    p <- p + geom_errorbar(aes(ymin=y-ci, ymax=y+ci), width=0, colour="black")
    p <- p + geom_point(aes(size=1.2), colour="blue") + theme(axis.text.x = element_text(angle=90, vjust=0.5))

    return(p)
    #ggsave(paste0(output,"dotplots.pdf"),p)
  }

  lapply(re,f)
}

se.ranef <- function (object) {
  se.bygroup <- ranef(object, condVar = TRUE)
  n.groupings <- length(se.bygroup)
  for (m in 1:n.groupings) {
    vars.m <- attr(se.bygroup[[m]], "postVar")
    K <- dim(vars.m)[1]
    J <- dim(vars.m)[3]
    se.bygroup[[m]] <- array(NA, c(J, K))
    for (j in 1:J) {
      se.bygroup[[m]][j, ] <- sqrt(diag(as.matrix(vars.m[,, j])))  }
    names.full <- dimnames(se.bygroup)
    dimnames(se.bygroup[[m]]) <- list(names.full[[1]], names.full[[2]]) }
  return(se.bygroup) }

# create nloptwrap optimization method
defaultControl <- list(algorithm="NLOPT_LN_BOBYQA",xtol_rel=1e-6,maxeval=1e5)
nloptwrap2 <- function(fn,par,lower,upper,control=list(),...) {
  for (n in names(defaultControl))
    if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
    res <- nloptr(x0=par,eval_f=fn,lb=lower,ub=upper,opts=control,...)
    with(res,list(par=solution,
                  fval=objective,
                  feval=iterations,
                  conv=if (status>0) 0 else status,
                  message=message))
}

augment.ranef.mer <- function(x,
                              ci.level=0.9,
                              reorder=TRUE,
                              order.var=1) {
  tmpf <- function(z) {
    if (is.character(order.var) && !order.var %in% names(z)) {
      order.var <- 1
      warning("order.var not found, resetting to 1")
    }
    ## would use plyr::name_rows, but want levels first
    zz <- data.frame(level=rownames(z),z,check.names=FALSE)
    if (reorder) {
      ## if numeric order var, add 1 to account for level column
      ov <- if (is.numeric(order.var)) order.var+1 else order.var
      zz$level <- reorder(zz$level, zz[,order.var+1], FUN=identity)
    }
    ## Q-Q values, for each column separately
    qq <- c(apply(z,2,function(y) {
      qnorm(ppoints(nrow(z)))[order(order(y))]
    }))
    rownames(zz) <- NULL
    pv   <- attr(z, "postVar")
    cols <- 1:(dim(pv)[1])
    se   <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
    ## n.b.: depends on explicit column-major ordering of se/melt
    zzz <- cbind(melt(zz,id.vars="level",value.name="estimate"),
                 qq=qq,std.error=se)
    ## reorder columns:
    subset(zzz,select=c(variable, level, estimate, qq, std.error))
  }
  dd <- ldply(x,tmpf,.id="grp")
  ci.val <- -qnorm((1-ci.level)/2)
  transform(dd,
            p=2*pnorm(-abs(estimate/std.error)), ## 2-tailed p-val
            lb=estimate-ci.val*std.error,
            ub=estimate+ci.val*std.error)
}

"""

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
#---------------------------------------------------------------------#
#This function performs the time-integration
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
include("ti_SSP.jl")
include("ti_LSRK.jl")

function time_integration!(q0,u,Dhat,Fhat,intma,periodicity,time,ntime,dt,stages,DFloat,ti_type)

    if (ti_type == "SSP")
        q0 = ti_SSP!(q0,u,Dhat,Fhat,intma,periodicity,time,ntime,dt,stages,DFloat)
    elseif (ti_type == "LSRK")
        q0 = ti_LSRK!(q0,u,Dhat,Fhat,intma,periodicity,time,ntime,dt,DFloat)
    end
    return (q0)
end


text:
export teval, @teval, evalfast, @evalfast, tgetfield, @tgetfield, _teval, @_teval

# overloading Meta.parse to support Tokens
Meta.parse(t::T) where {T <: Union{Token, Array{Token}}} = Meta.parse(untokenize(t))

###############################################################
"""
    @teval(t, check_isdefined::Bool=false)
    teval(modul, t::T, check_isdefined::Bool=false)

Parses and evaluates a `Token` (or an experssion or a string) that represents a `Symbol` or `Expr` in Julia. It dispatches on the fast method based on the parsed Token.

Use `@teval`, if you want the token be evaluated in the called scope.

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```jldoctest
julia> t = collect(tokenize("Int64"))

julia> @teval(t)
Int64
```
"""
function teval(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}, String}}
    pt = Meta.parse(t)
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken()
    end
    return evalfast(modul,pt)
end

function teval(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Symbol, Expr}}
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken()
    end
    return evalfast(modul,pt)
end

macro teval(t, check_isdefined::Bool = false) # where {T <: Union{Symbol, Array{Token}, String, Symbol, Expr}}
    macroexpand(__module, t)
    return quote
        TokenizeMeta.teval($__module__, $(Expr(:quote,t)), $check_isdefined)
    end
end
#
# macro teval(mod, ex)
#     :(Core.eval($(esc(mod)), $(Expr(:quote,ex))))
# end

################################################################

"""
    @evalfast(x)
    evalfast(modul, x)

Evaluates `x` fast. Uses `getfield` if `x` is a `Symbol`, and uses `eval` if it is an `Expr`
# Examples
```jldoctest
julia> @evalfast(Int64)

julia> @evalfast([5])

```
"""
evalfast(modul::Module, x::Expr)= modul.eval(x)
evalfast(modul::Module, x::Symbol)= getfield(modul,x)
macro evalfast(x::Expr)
    m = __module__
    return :(Core.eval($m, $x))
end
macro evalfast(x::Symbol)
    m = __module__
    return :(getfield($m,$x))
end

################################################################
"""
    @tgetfield(t, check_isdefined::Bool=false)
    tgetfield(modul, t, check_isdefined::Bool=false)

See [`teval`](@ref) for fast Token evaluation.

Parses and evaluates a `Token` (or an experssion or a string) that represents a `Symbol` in Julia. For `Symbol` it is similar to eval, but much faster.

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```jldoctest
julia> t = collect(tokenize("Int64"))

julia> @tgetfield(t)
Int64
```
"""
function tgetfield(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}, String}}
    pt = Meta.parse(t)
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken()
    end
    return getfield(modul,pt)
end

function tgetfield(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Symbol, Expr}}
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken()
    end
    return getfield(modul,pt)
end

macro tgetfield(t, check_isdefined::Bool = false)
    m = __module__
    return :(tgetfield($m, $t, $check_isdefined))
end

################################################################
"""
    @_teval(t, check_isdefined::Bool=false)
    _teval(modul, t, check_isdefined::Bool=false)

See [`teval`](@ref) for fast Token evaluation.

Parses and evaluates a `Token` (or an expression or a string).

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```jldoctest
julia> t = collect(tokenize("Int64"))

julia> @_teval(t)
Int64
```
"""
function _teval(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}, String}}
    pt = Meta.parse(t)
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken()
    end
    return Core.eval(modul, pt)
end

function _teval(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Symbol, Expr}}
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken()
    end
    return Core.eval(modul, pt)
end

macro _teval(t, check_isdefined::Bool = false)
    m = __module__
    return :(teval($m, $t, $check_isdefined))
end


text: # Citation.jl


text:
# Functions related to sterographic projection to/from a three-dimensional
# unit sphere centered at (0,0) = 0+0im

# function stereo(z::Number)::Vector{Float64}
#     if isinf(z)
#         return [0, 0, 1]
#     end

#     X, Y = reim(z)

#     x = 2X
#     y = 2Y
#     z = X^2 + Y^2 - 1

#     d = 1 / (1 + X^2 + Y^2)

#     return d * [x, y, z]
# end

"""
    stereo
Stereographic projection between the complex plane and a three-dimensional unit sphere
centered at `[0,0,0]`. The north pole, `[0,0,1]`, corresponds to complex infinity and the
south pole, `[0,0,-1]`, corresponds to `0+0im`.

For a complex number `z`, `stereo(z)` maps `z` to the sphere. This may also be invoked as
`stereo(x,y)`.

For a (unit) three-dimensional vector `v`, `stereo(v)` returns the complex number by projecting
`v` to the complex plane. This may also be invoked as `stereo(x,y,z)`.
"""
function stereo(X::Real, Y::Real)::Vector{Float64}
    if isinf(X) || isinf(Y)
        return [0, 0, 1]
    end

    x = 2X
    y = 2Y
    z = X^2 + Y^2 - 1

    d = 1 / (1 + X^2 + Y^2)

    return d * [x, y, z]
end

stereo(z::Number) = stereo(reim(z)...)

function stereo(x::Real, y::Real, z::Real)::Complex
    if z == 1
        return Inf + im * Inf
    end
    X = x / (1 - z)
    Y = y / (1 - z)
    return X + im * Y
end

function stereo(v::Vector{T})::Complex where {T<:Real}
    stereo(v...)
end

"""
    LFTQ
Create a linear fractional transformation from a 3-by-3 unitary matrix.

Given a 3-by-3 real matrix `Q` with `Q*Q'` equal to the identity and `det(Q)` equal to `1`,
create a `LFT` that maps a complex number `v` to
`stereo(Q*stereo(v))`.
"""
function LFTQ(Q::AbstractMatrix)::LFT
    Q = Matrix(Q)
    zz = [0 + 0im, 1 + 0im, 0 + im]

    uu = stereo.(zz)
    uu = [Q * u for u in uu]

    ww = [stereo(u) for u in uu]

    LFT(zz[1], ww[1], zz[2], ww[2], zz[3], ww[3])

end

export stereo, LFTQ


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# RUN: julia --startup-file=no %s | opt -load libjulia%shlibext -AllocOpt -LateLowerGCFrame -S - | FileCheck %s

isz = sizeof(UInt) == 8 ? "i64" : "i32"

println("""
%jl_value_t = type opaque
@tag = external addrspace(10) global %jl_value_t
""")

# CHECK-LABEL: @return_obj
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK: %v = call noalias nonnull %jl_value_t addrspace(10)* @jl_gc_pool_alloc
# CHECK: store %jl_value_t addrspace(10)* @tag, %jl_value_t addrspace(10)* addrspace(10)* {{.*}}, !tbaa !0
println("""
define %jl_value_t addrspace(10)* @return_obj() {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  ret %jl_value_t addrspace(10)* %v
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @return_load
# CHECK: alloca i64
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK-NOT: @jl_gc_pool_alloc
# CHECK: call void @llvm.lifetime.start{{.*}}(i64 8, i8*
# CHECK-NEXT: %v64 = bitcast %jl_value_t* %v to i64*
# CHECK-NOT: @tag
# CHECK-NOT: @llvm.lifetime.end
println("""
define i64 @return_load(i64 %i) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %v64 = bitcast %jl_value_t addrspace(10)* %v to i64 addrspace(10)*
  %v64a11 = addrspacecast i64 addrspace(10)* %v64 to i64 addrspace(11)*
  store i64 %i, i64 addrspace(11)* %v64a11, align 16, !tbaa !4
  call void @external_function()
  %l = load i64, i64 addrspace(11)* %v64a11, align 16, !tbaa !4
  ret i64 %l
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @ccall_obj
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK: @jl_gc_pool_alloc
# CHECK: store %jl_value_t addrspace(10)* @tag, %jl_value_t addrspace(10)* addrspace(10)* {{.*}}, !tbaa !0
println("""
define void @ccall_obj(i8* %fptr) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %f = bitcast i8* %fptr to void (%jl_value_t addrspace(10)*)*
  call void %f(%jl_value_t addrspace(10)* %v)
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @ccall_ptr
# CHECK: alloca i64
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK-NOT: @jl_gc_pool_alloc
# CHECK: call void @llvm.lifetime.start{{.*}}(i64 8, i8*
# CHECK: %f = bitcast i8* %fptr to void (i8*)*
# Currently the GC frame lowering pass strips away all operand bundles
# CHECK-NEXT: call void %f(i8*
# CHECK-NEXT: ret void
println("""
define void @ccall_ptr(i8* %fptr) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %va = addrspacecast %jl_value_t addrspace(10)* %v to %jl_value_t addrspace(11)*
  %ptrj = call %jl_value_t* @julia.pointer_from_objref(%jl_value_t addrspace(11)* %va)
  %ptr = bitcast %jl_value_t* %ptrj to i8*
  %f = bitcast i8* %fptr to void (i8*)*
  call void %f(i8* %ptr) [ "jl_roots"(%jl_value_t addrspace(10)* %v), "unknown_bundle"(i8* %ptr) ]
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @ccall_unknown_bundle
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK: @jl_gc_pool_alloc
# CHECK: store %jl_value_t addrspace(10)* @tag, %jl_value_t addrspace(10)* addrspace(10)* {{.*}}, !tbaa !0
println("""
define void @ccall_unknown_bundle(i8* %fptr) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %va = addrspacecast %jl_value_t addrspace(10)* %v to %jl_value_t addrspace(11)*
  %ptrj = call %jl_value_t* @julia.pointer_from_objref(%jl_value_t addrspace(11)* %va)
  %ptr = bitcast %jl_value_t* %ptrj to i8*
  %f = bitcast i8* %fptr to void (i8*)*
  call void %f(i8* %ptr) [ "jl_not_jl_roots"(%jl_value_t addrspace(10)* %v) ]
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @lifetime_branches
# CHECK: alloca i64
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK: L1:
# CHECK-NEXT: call void @llvm.lifetime.start{{.*}}(i64 8,
# CHECK: %f = bitcast i8* %fptr to void (i8*)*
# CHECK-NEXT: call void %f(i8*
# CHECK-NEXT: br i1 %b2, label %L2, label %L3

# CHECK: L2:
# CHECK-NEXT: %f2 = bitcast i8* %fptr to void (%jl_value_t*)*
# CHECK-NEXT: call void @llvm.lifetime.end{{.*}}(i64 8,
# CHECK-NEXT: call void %f2(%jl_value_t* null)

# CHECK: L3:
# CHECK-NEXT: call void @llvm.lifetime.end{{.*}}(i64 8,
println("""
define void @lifetime_branches(i8* %fptr, i1 %b, i1 %b2) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  br i1 %b, label %L1, label %L3

L1:
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %va = addrspacecast %jl_value_t addrspace(10)* %v to %jl_value_t addrspace(11)*
  %ptrj = call %jl_value_t* @julia.pointer_from_objref(%jl_value_t addrspace(11)* %va)
  %ptr = bitcast %jl_value_t* %ptrj to i8*
  %f = bitcast i8* %fptr to void (i8*)*
  call void %f(i8* %ptr) [ "jl_roots"(%jl_value_t addrspace(10)* %v) ]
  br i1 %b2, label %L2, label %L3

L2:
  %f2 = bitcast i8* %fptr to void (%jl_value_t*)*
  call void %f2(%jl_value_t* null)
  br label %L3

L3:
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @object_field
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK: @jl_gc_pool_alloc
# CHECK: store %jl_value_t addrspace(10)* @tag, %jl_value_t addrspace(10)* addrspace(10)* {{.*}}, !tbaa !0
println("""
define void @object_field(%jl_value_t addrspace(10)* %field) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %va = addrspacecast %jl_value_t addrspace(10)* %v to %jl_value_t addrspace(11)*
  %vab = bitcast %jl_value_t addrspace(11)* %va to %jl_value_t addrspace(10)* addrspace(11)*
  store %jl_value_t addrspace(10)* %field, %jl_value_t addrspace(10)* addrspace(11)* %vab
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @memcpy_opt
# CHECK: alloca i128, align 16
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK-NOT: @jl_gc_pool_alloc
# CHECK: call void @llvm.memcpy.p0i8.p0i8.i64
println("""
define void @memcpy_opt(i8* %v22) {
top:
  %v6 = call %jl_value_t*** @julia.ptls_states()
  %v18 = bitcast %jl_value_t*** %v6 to i8*
  %v19 = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %v18, $isz 16, %jl_value_t addrspace(10)* @tag)
  %v20 = bitcast %jl_value_t addrspace(10)* %v19 to i8 addrspace(10)*
  %v21 = addrspacecast i8 addrspace(10)* %v20 to i8 addrspace(11)*
  call void @llvm.memcpy.p11i8.p0i8.i64(i8 addrspace(11)* %v21, i8* %v22, i64 16, i32 8, i1 false)
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @preserve_opt
# CHECK: alloca i128, align 16
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK-NOT: @julia.gc_alloc_obj
# CHECK-NOT: @jl_gc_pool_alloc
# CHECK: @llvm.lifetime.start
# CHECK-NOT: @llvm.lifetime.end
# CHECK: @external_function
# CHECK-NEXT: @llvm.lifetime.end
# CHECK-NEXT: @external_function
println("""
define void @preserve_opt(i8* %v22) {
top:
  %v6 = call %jl_value_t*** @julia.ptls_states()
  %v18 = bitcast %jl_value_t*** %v6 to i8*
  %v19 = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %v18, $isz 16, %jl_value_t addrspace(10)* @tag)
  %v20 = bitcast %jl_value_t addrspace(10)* %v19 to i8 addrspace(10)*
  %v21 = addrspacecast i8 addrspace(10)* %v20 to i8 addrspace(11)*
  %tok = call token (...) @llvm.julia.gc_preserve_begin(%jl_value_t addrspace(10)* %v19)
  call void @external_function()
  call void @llvm.julia.gc_preserve_end(token %tok)
  call void @external_function()
  ret void
}
""")
# CHECK-LABEL: }

# CHECK-LABEL: @preserve_branches
# CHECK: alloca i64
# CHECK: call %jl_value_t*** @julia.ptls_states()
# CHECK: L1:
# CHECK-NEXT: call void @llvm.lifetime.start{{.*}}(i64 8,
# CHECK-NEXT: @external_function()
# CHECK-NEXT: br i1 %b2, label %L2, label %L3

# CHECK: L2:
# CHECK-NOT: call void @llvm.lifetime.end{{.*}}(i64 8,
# CHECK: @external_function()
# CHECK-NEXT: br label %L3

# CHECK: L3:
# CHECK-NEXT: call void @llvm.lifetime.end{{.*}}(i64 8,
println("""
define void @preserve_branches(i8* %fptr, i1 %b, i1 %b2) {
  %ptls = call %jl_value_t*** @julia.ptls_states()
  %ptls_i8 = bitcast %jl_value_t*** %ptls to i8*
  br i1 %b, label %L1, label %L3

L1:
  %v = call noalias %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, %jl_value_t addrspace(10)* @tag)
  %tok = call token (...) @llvm.julia.gc_preserve_begin(%jl_value_t addrspace(10)* %v)
  call void @external_function()
  br i1 %b2, label %L2, label %L3

L2:
  call void @external_function()
  br label %L3

L3:
  ret void
}
""")
# CHECK-LABEL: }

# CHECK: declare noalias nonnull %jl_value_t addrspace(10)* @jl_gc_pool_alloc(i8*,
# CHECK: declare noalias nonnull %jl_value_t addrspace(10)* @jl_gc_big_alloc(i8*,
println("""
declare void @external_function()
declare %jl_value_t*** @julia.ptls_states()
declare noalias nonnull %jl_value_t addrspace(10)* @julia.gc_alloc_obj(i8*, $isz, %jl_value_t addrspace(10)*)
declare %jl_value_t* @julia.pointer_from_objref(%jl_value_t addrspace(11)*)
declare void @llvm.memcpy.p11i8.p0i8.i64(i8 addrspace(11)* nocapture writeonly, i8* nocapture readonly, i64, i32, i1)
declare token @llvm.julia.gc_preserve_begin(...)
declare void @llvm.julia.gc_preserve_end(token)

!0 = !{!1, !1, i64 0}
!1 = !{!"jtbaa_tag", !2, i64 0}
!2 = !{!"jtbaa_data", !3, i64 0}
!3 = !{!"jtbaa"}
!4 = !{!5, !5, i64 0}
!5 = !{!"jtbaa_mutab", !6, i64 0}
!6 = !{!"jtbaa_value", !2, i64 0}
""")

max_stars_repo_name: ylxdzsw/julia

max_issues_repo_name: ylxdzsw/julia

max_forks_repo_name: ylxdzsw/julia


text:
using DataStructures

slurp_grid(filename) = parse.(Int8, permutedims(hcat(collect.(readlines(filename))...)))

get_neighbors(idx) = [idx + CartesianIndex(Δx, Δy) for Δx = -1:1 for Δy in -1:1 if !(Δx == Δy == 0)]

function in_bounds(grid, idx)
    y_max, x_max = size(grid)
    y, x = Tuple(idx)
    1 <= y <= y_max && 1 <= x <= x_max
end

get_neighbors(grid, idx) = [idx for idx in get_neighbors(idx) if in_bounds(grid, idx)]

get_neighbors(CartesianIndex(1, 1))

function step(grid, num_flashes)
    flashed = []
    to_increment = []
    grid = (i -> i + 1).(grid)
    for idx in CartesianIndices(grid)
        if grid[idx] > 9
            push!(flashed, idx)
            for idx_ in get_neighbors(grid, idx)
                if idx_ ∉ flashed
                    push!(to_increment, idx_)
                end
            end
        end
    end
    while !isempty(to_increment)
        next = pop!(to_increment)
        grid[next] += 1
        if grid[next] > 9 && next ∉ flashed
            push!(flashed, next)
            for idx in get_neighbors(grid, next)
                if idx ∉ flashed
                    push!(to_increment, idx)
                end
            end
        end
    end
    grid = (i -> i > 9 ? 0 : i).(grid)
    grid, num_flashes + length(flashed)
end

function solve1(filename)
    grid = slurp_grid(filename)
    num_flashes = 0
    for _ = 1:100
        grid, num_flashes = step(grid, num_flashes)
    end
    num_flashes
end

solve1("11eg.txt") == 1656
solve1("11.txt")

function step(grid)
    flashed = []
    to_increment = []
    grid = (i -> i + 1).(grid)
    for idx in CartesianIndices(grid)
        if grid[idx] > 9
            push!(flashed, idx)
            for idx_ in get_neighbors(grid, idx)
                if idx_ ∉ flashed
                    push!(to_increment, idx_)
                end
            end
        end
    end
    while !isempty(to_increment)
        next = pop!(to_increment)
        grid[next] += 1
        if grid[next] > 9 && next ∉ flashed
            push!(flashed, next)
            for idx in get_neighbors(grid, next)
                if idx ∉ flashed
                    push!(to_increment, idx)
                end
            end
        end
    end
    (i -> i > 9 ? 0 : i).(grid)
end

function solve2(filename)
    grid = slurp_grid(filename)
    count = 0
    while !all(==(0), grid)
        grid = step(grid)
        count += 1
    end
    count
end

solve2("11eg.txt")
solve2("11.txt")