text:
# Code I needed to add to make conversion from R easier, some cases are just function
# renaming.  Others are major functions.

function factor(v; ordered=false, levels=nothing)
    r = CategoricalArray(v, ordered=ordered)
    levels != nothing && (levels!(r, levels))
    return(r)
end

# No duplicated method?  (add pull request to julialang for this)
# This vectorized method should work on anything that can
# * Fit in an array
# * Can be used as a key in a Dict
# fromLast = true will save last instance of a value
# fromlast is an annoying scenario as we need to flip twice
function duplicated(a::AbstractArray; fromLast = false)
    d = Dict()
    # So we only need one maping line
    @inline getA(x) = fromLast ? reverse(x) : (x)
    getA(map(x -> ((haskey(d, x)) ? true : (d[x] = true; false) ),  getA(a)))
end

#@inline factor=CategoricalArray
# A set of methods for transforming OPM Data for Analysis
# This is a port of R code, some functions are created to replace default R functions.
# Define the factors we will use, these are just pre-generated from the source data.
OPMAgeFactorLevels = vcat("15-19", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60-64", "65-69", "70-74", "75+")
OPMLOSFactorLevels = vcat("< 1", "1-2", "3-4", "5-9", "10-14", "15-19", "20-24", "25-29", "30-34", "35+")

# Changes a vector to a factor based on age ranges
function opm_factor_age(v)
    factor(v, levels = OPMAgeFactorLevels, ordered = true)
end

# Changes a vector to a factor based on los ranges
function opm_factor_los(v)
    factor(v, levels = OPMLOSFactorLevels, ordered = true)
end

# Converts nsftp vector to a logical
function opm_indicator_nsftp(v)
    ifelse.(v .== "2", true, false)
end

#Takes a tibble, assumed to have OPM fields in it, and cleans it up for data analysis
function opm_organize_tibble(filedata, sctdata)
 #= Very long pipeline. Leving this for last since it is hardest
    filedata %>%
    # Put in factors for Lenth of Serivce and Age
    mutate(AGE = opm_factor_age(AGE)) %>%
    mutate(`LOS LEVEL` = opm_factor_los(`LOS LEVEL`)) %>%
    # Change NSFTP to an an indicator
    mutate(NSFTP = opm_indicator_nsftp(`NSFTP IND`)) %>%
    # Hard work!  Change lookup fields into the sctfile into factors
    left_join(opm_sct_get_workschedules(sctdata), by = c("WORK SCHEDULE" = "DATA CODE")) %>%
    mutate(`WORK SCHEDULE` = factor(`WORK SCHEDULE.y`, opm_sct_get_workschedules_values(sctdata))) %>%
    select(-`WORK SCHEDULE.y`) %>%
    left_join(opm_sct_get_appointments(sctdata), by = c("TYPE OF APPOINTMENT" = "DATA CODE")) %>%
    mutate(`TYPE OF APPOINTMENT` = factor(`APPOINTMENT`, opm_sct_get_appointments_values(sctdata))) %>%
    select(-`APPOINTMENT`) %>%
    left_join(opm_sct_get_supervisorystatus(sctdata), by = c("SUPERVISORY STATUS" = "DATA CODE")) %>%
    mutate(`SUPERVISORY STATUS` = factor(`SUPERVISOR`, opm_sct_get_supervisorystatus_values(sctdata))) %>%
    select(-`SUPERVISOR`) %>%
    left_join(opm_sct_get_occupationalcategories(sctdata), by = c("OCCUPATIONAL CATEGORY" = "DATA CODE")) %>%
    mutate(`OCCUPATIONAL CATEGORY` = factor(`OCCUPATIONAL CATEGORY.y`, opm_sct_get_occupationalcategories_values(sctdata))) %>%
    select(-`OCCUPATIONAL CATEGORY.y`) %>%
    left_join(opm_sct_get_occupations(sctdata), by = c("OCCUPATION" = "DATA CODE")) %>%
    mutate(`OCCUPATION` = factor(`OCCUPATION.y`, opm_sct_get_occupations_values(sctdata))) %>%
    select(-`OCCUPATION.y`) %>%
    left_join(opm_sct_get_payplans(sctdata), by = c("PAY PLAN" = "DATA CODE")) %>%
    mutate(`PAY PLAN` = factor(`PAY PLAN.y`, opm_sct_get_payplans_values(sctdata))) %>%
    select(-`PAY PLAN.y`) %>%
    left_join(opm_sct_get_eductionlevels(sctdata), by = c("EDUCATION LEVEL" = "DATA CODE")) %>%
    mutate(`EDUCATION LEVEL` = factor(`EDUCATION LEVEL.y`, opm_sct_get_eductionlevels_values(sctdata))) %>%
    select(-`EDUCATION LEVEL.y`) %>%
    # Duty stations seems to break code parsing, guessing a bad character in one of the names.
    # TODO: try this again, think it was a code bug in another section overwriting a name.
    #left_join(opm_sct_get_dutystations(sctdata), by = c("DUTY STATION" = "DATA CODE")) %>%
    #mutate(`DUTY STATION` = factor(`DUTY STATION.y`, opm_sct_get_dutystations_values(sctdata))) #%>%
    #select(-`DUTY STATION.y`) %>%
    left_join(opm_sct_get_agencies(sctdata), by = c("AGENCY" = "DATA CODE")) %>%
    mutate(`AGENCY` = factor(`AGENCY.y`, opm_sct_get_agencies_values(sctdata))) %>%
    select(-`AGENCY.y`)
    =#
    # Put in factors for Lenth of Serivce and Age
    #mutate(AGE = opm_factor_age(AGE)) %>%
    filedata[:AGE] = opm_factor_age(filedata[:AGE])
    #mutate(`LOS LEVEL` = opm_factor_los(`LOS LEVEL`)) %>%
    filedata[:LOS_LEVEL] = opm_factor_los(filedata[:LOS_LEVEL])
    # Change NSFTP to an an indicator
    #mutate(NSFTP = opm_indicator_nsftp(`NSFTP IND`)) %>%
    filedata[:NSFTP_IND] = opm_indicator_nsftp(filedata[:NSFTP_IND])
    # Hard work!  Change lookup fields into the sctfile into factors
    #left_join(opm_sct_get_workschedules(sctdata), by = c("WORK SCHEDULE" = "DATA CODE")) %>%
    #mutate(`WORK SCHEDULE` = factor(`WORK SCHEDULE.y`, opm_sct_get_workschedules_values(sctdata))) %>%
    #select(-`WORK SCHEDULE.y`) %>%
    # 4julia - these 3 line groups basically performed an in-place left join / replace
    filedata[:WORK_SCHEDULE] =
        factor(join(DataFrame(X1=filedata[:WORK_SCHEDULE]),
            opm_sct_get_workschedules(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:WORK_SCHEDULE],
            levels=opm_sct_get_workschedules_values(sctdata))

    #left_join(opm_sct_get_appointments(sctdata), by = c("TYPE OF APPOINTMENT" = "DATA CODE")) %>%
    #mutate(`TYPE OF APPOINTMENT` = factor(`APPOINTMENT`, opm_sct_get_appointments_values(sctdata))) %>%
    #select(-`APPOINTMENT`) %>%
    filedata[:TYPE_OF_APPOINTMENT] =
        factor(join(DataFrame(X1=filedata[:TYPE_OF_APPOINTMENT]),
            opm_sct_get_appointments(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:APPOINTMENT],
            levels=opm_sct_get_appointments_values(sctdata))
    #left_join(opm_sct_get_supervisorystatus(sctdata), by = c("SUPERVISORY STATUS" = "DATA CODE")) %>%
    #mutate(`SUPERVISORY STATUS` = factor(`SUPERVISOR`, opm_sct_get_supervisorystatus_values(sctdata))) %>%
    #select(-`SUPERVISOR`) %>%
    filedata[:SUPERVISORY_STATUS] =
        factor(join(DataFrame(X1=filedata[:SUPERVISORY_STATUS]),
            opm_sct_get_supervisorystatus(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:SUPERVISOR],
            levels=opm_sct_get_supervisorystatus_values(sctdata))

    #left_join(opm_sct_get_occupationalcategories(sctdata), by = c("OCCUPATIONAL CATEGORY" = "DATA CODE")) %>%
    #mutate(`OCCUPATIONAL CATEGORY` = factor(`OCCUPATIONAL CATEGORY.y`, opm_sct_get_occupationalcategories_values(sctdata))) %>%
    #select(-`OCCUPATIONAL CATEGORY.y`) %>%
    filedata[:OCCUPATIONAL_CATEGORY] =
        factor(join(DataFrame(X1=filedata[:OCCUPATIONAL_CATEGORY]),
            opm_sct_get_occupationalcategories(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:OCCUPATIONAL_CATEGORY],
            levels=opm_sct_get_occupationalcategories_values(sctdata))

    #left_join(opm_sct_get_occupations(sctdata), by = c("OCCUPATION" = "DATA CODE")) %>%
    #mutate(`OCCUPATION` = factor(`OCCUPATION.y`, opm_sct_get_occupations_values(sctdata))) %>%
    #select(-`OCCUPATION.y`) %>%
    filedata[:OCCUPATION] =
        factor(join(DataFrame(X1=filedata[:OCCUPATION]),
            opm_sct_get_occupations(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:OCCUPATION],
            levels=opm_sct_get_occupations_values(sctdata))

    #left_join(opm_sct_get_payplans(sctdata), by = c("PAY PLAN" = "DATA CODE")) %>%
    #mutate(`PAY PLAN` = factor(`PAY PLAN.y`, opm_sct_get_payplans_values(sctdata))) %>%
    #select(-`PAY PLAN.y`) %>%
    filedata[:PAY_PLAN] =
        factor(join(DataFrame(X1=filedata[:PAY_PLAN]),
            opm_sct_get_payplans(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:PAY_PLAN],
            levels=opm_sct_get_payplans_values(sctdata))

    #left_join(opm_sct_get_eductionlevels(sctdata), by = c("EDUCATION LEVEL" = "DATA CODE")) %>%
    #mutate(`EDUCATION LEVEL` = factor(`EDUCATION LEVEL.y`, opm_sct_get_eductionlevels_values(sctdata))) %>%
    #select(-`EDUCATION LEVEL.y`) %>%
    filedata[:EDUCATION_LEVEL] =
        factor(join(DataFrame(X1=filedata[:EDUCATION_LEVEL]),
            opm_sct_get_eductionlevels(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:EDUCATION_LEVEL],
            levels=opm_sct_get_eductionlevels_values(sctdata))

    # Duty stations seems to break code parsing, guessing a bad character in one of the names.
    # TODO: try this again, think it was a code bug in another section overwriting a name.
    #left_join(opm_sct_get_dutystations(sctdata), by = c("DUTY STATION" = "DATA CODE")) %>%
    #mutate(`DUTY STATION` = factor(`DUTY STATION.y`, opm_sct_get_dutystations_values(sctdata))) #%>%
    #select(-`DUTY STATION.y`) %>%

    #left_join(opm_sct_get_agencies(sctdata), by = c("AGENCY" = "DATA CODE")) %>%
    #mutate(`AGENCY` = factor(`AGENCY.y`, opm_sct_get_agencies_values(sctdata))) %>%
    #select(-`AGENCY.y`)
    filedata[:AGENCY] =
        factor(join(DataFrame(X1=filedata[:AGENCY]),
            opm_sct_get_agencies(sctdata), on=(:X1, :DATA_CODE), kind=:left)[:AGENCY],
            levels=opm_sct_get_agencies_values(sctdata))

    return filedata
end

# Will pull all values for a table ID out of tibble parameter sctfile and return them
# sctfile tibble of data to pull from
# tblanme two character code of table to pull
function opm_sct_get_table(sctfile, tblname)
    #sctfile %>%
    #    filter(`SCT TABLE ID` == tblname) %>%
    #    select(c("DATA CODE", "TRANSLATION 1"))
    @from i in sctfile begin
    @where i.SCT_TABLE_ID == tblname
    @select {i.DATA_CODE, i.TRANSLATION_1}
    @collect DataFrame
    end
end

#Pulls appointments "VM" out of the SCT DATA, assums SCT structure of tibble
function opm_sct_get_appointments(sctfile)
    #opm_sct_get_table(sctfile, "VM") %>% rename("APPOINTMENT" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "VM"), :TRANSLATION_1 => :APPOINTMENT)
end

function opm_sct_get_appointments_values(sctfile)
    #opm_sct_get_appointments(sctfile) %>% pull("APPOINTMENT")
    opm_sct_get_appointments(sctfile)[:APPOINTMENT]
end

function opm_sct_get_eductionlevels(sctfile)
    #opm_sct_get_table(sctfile, "EV") %>% rename("EDUCATION LEVEL" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "EV"), :TRANSLATION_1 => :EDUCATION_LEVEL)
end

function opm_sct_get_eductionlevels_values(sctfile)
    #opm_sct_get_eductionlevels(sctfile) %>% pull("EDUCATION LEVEL")
    opm_sct_get_eductionlevels(sctfile)[:EDUCATION_LEVEL]
end

function opm_sct_get_occupationalcategories(sctfile)
    #opm_sct_get_table(sctfile, "GF") %>% rename("OCCUPATIONAL CATEGORY" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "GF"), :TRANSLATION_1 => :OCCUPATIONAL_CATEGORY)
end

function opm_sct_get_occupationalcategories_values(sctfile)
    #opm_sct_get_occupationalcategories(sctfile) %>% pull("OCCUPATIONAL CATEGORY")
    opm_sct_get_occupationalcategories(sctfile)[:OCCUPATIONAL_CATEGORY]
end

function opm_sct_get_payplans(sctfile)
    #opm_sct_get_table(sctfile, "LA") %>% rename("PAY PLAN" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "LA"), :TRANSLATION_1 => :PAY_PLAN)
end

# There are duplicate pay plan names for codes, executive decision: just remove duplicates
function opm_sct_get_payplans_values(sctfile)
    #opm_sct_get_payplans(sctfile) %>% pull("PAY PLAN") %>% unique()
    unique(opm_sct_get_payplans(sctfile)[:PAY_PLAN])
end

function opm_sct_get_supervisorystatus(sctfile)
    #opm_sct_get_table(sctfile, "SU") %>% rename("SUPERVISOR" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "SU"), :TRANSLATION_1 => :SUPERVISOR)
end

function opm_sct_get_supervisorystatus_values(sctfile)
    #opm_sct_get_supervisorystatus(sctfile) %>% pull("SUPERVISOR")
    opm_sct_get_supervisorystatus(sctfile)[:SUPERVISOR]
end

function opm_sct_get_workschedules(sctfile)
    #opm_sct_get_table(sctfile, "WS") %>% rename("WORK SCHEDULE" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "WS"), :TRANSLATION_1 => :WORK_SCHEDULE)
end

function opm_sct_get_workschedules_values(sctfile)
    #opm_sct_get_workschedules(sctfile) %>% pull("WORK SCHEDULE")
    opm_sct_get_workschedules(sctfile)[:WORK_SCHEDULE]
end

function opm_sct_get_occupations(sctfile)
    #opm_sct_get_table(sctfile, "XB") %>% rename("OCCUPATION" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "XB"), :TRANSLATION_1 => :OCCUPATION)
end

# There are duplicate occupation names for codes, executive decision: just remove duplicates
function opm_sct_get_occupations_values(sctfile)
    #opm_sct_get_occupations(sctfile) %>% pull("OCCUPATION") %>% unique()
    unique(opm_sct_get_occupations(sctfile)[:OCCUPATION])
end

function opm_sct_get_dutystations(sctfile)
    #opm_sct_get_table(sctfile, "VX") %>% rename("DUTY STATION" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "VX"), :TRANSLATION_1 => :DUTY_STATION)
end

function opm_sct_get_dutystations_values(sctfile)
    #opm_sct_get_dutystations(sctfile) %>% pull("DUTY STATION")
    opm_sct_get_dutystations(sctfile)[:DUTY_STATION]
end

# Agencies is a bit tricker than others, it has duplicates so we need to de,dupe it in a most likely not correct manner
function opm_sct_get_agencies(sctfile)
    #tmp <- opm_sct_get_table(sctfile, "AH") %>% rename("AGENCY" = `TRANSLATION 1`)
    #tmp[!duplicated(tmp$`DATA CODE`, fromLast=TRUE),]
    tmp = opm_sct_get_table(sctfile, "AH")
    rename!(tmp, :TRANSLATION_1 => :AGENCY)
    tmp[.!duplicated(tmp[:DATA_CODE], fromLast=true), :]
end

function opm_sct_get_agencies_values(sctfile)
    #opm_sct_get_agencies(sctfile) %>% pull("AGENCY") %>% unique()
    unique(opm_sct_get_agencies(sctfile)[:AGENCY])
end

# THESE ARE NOT SALARY GRADES, NOT SURE WHAT THEY ARE
# Figured it out...they are combo codes for combined plan+grade lookup
function opm_sct_get_grades(sctfile)
    #opm_sct_get_table(sctfile, "VK") %>% rename("GRADE" = `TRANSLATION 1`)
    rename!(opm_sct_get_table(sctfile, "VK"), :TRANSLATION_1 => :GRADE)
end

function opm_sct_get_dutystations_values(sctfile)
    #opm_sct_get_grades(sctfile) %>% pull("GRADE")
    opm_sct_get_grades(sctfile)[:GRADE]
end


text:
using DDD, Plots, LinearAlgebra, Statistics, SparseArrays, StaticArrays, BenchmarkTools
plotlyjs()

## Parameters
dlnParams = DislocationParameters(;
    mobility = mobBCC(),
    dragCoeffs = (edge = 1, screw = 1e-1, climb = 1e9, line = 1e-5),
    coreRad = 0.015 * sqrt(2),
    minSegLen = 0.15 * sqrt(2),
    maxSegLen = 1.5 * sqrt(2),
    coreEnergy = 1 / (4 * π) * log(0.015 * sqrt(2) / 3.5e-5),
    coreRadMag = 3.5e-4,
);

matParams = MaterialParameters(; crystalStruct = BCC(), μ = 1, μMag = 80e3, ν = 0.25);

femParamsC = FEMParameters(;
    type = DispatchRegularCuboidMesh(),
    order = LinearElement(),
    model = CantileverLoad(),
    dx = 23.0,
    dy = 17.0,
    dz = 13.0,
    mx = 7,
    my = 5,
    mz = 3,
);

slipSystems = SlipSystem(;
    crystalStruct = BCC(),
    slipPlane = Float64[-1 1; 1 -1; 0 0],
    bVec = Float64[1 1; 1 1; 1 -1],
);

intParams = IntegrationParameters(;
    method = AdaptiveEulerTrapezoid(),
    abstol = dlnParams.collisionDist / 2,
    reltol = dlnParams.collisionDist / 2,
);

intTime = IntegrationTime(; dt = 0.0, time = 0.0);

# See fieldnames
fieldnames(typeof(matParams))

# # Browse docs/get variable info
# ?DislocationParameters

# ?femParamsC

## Build mesh & boundary conditions
meshC = buildMesh(matParams, femParamsC)
# Plot mesh
figFE = plotFEDomain(meshC; camera = (10, -2))
# savefig(figFE, "mesh.svg")

# Generate cantilever boundary conditions.
boundaryC, forceDispC = Boundaries(femParamsC, meshC)
# Plot boundaries
figCBC = plotBoundaries(boundaryC, meshC)
# savefig(figCBC, "cantilever.svg")

# Pillar loading
femParamsP = FEMParameters(;
    type = DispatchRegularCuboidMesh(),
    order = LinearElement(),
    model = PillarLoad(),
    dx = 17.0,
    dy = 13.0,
    dz = 23.0,
    mx = 5,
    my = 3,
    mz = 7,
);
meshP = buildMesh(matParams, femParamsP);
boundaryP, forceDispP = Boundaries(femParamsP, meshP);
figPBC = plotBoundaries(boundaryP, meshP)
# savefig(figPBC, "pillar.svg")

## Dislocations generation
dx, dy, dz = femParamsC.dx, femParamsC.dy, femParamsC.dz
segLen = (dlnParams.minSegLen + dlnParams.maxSegLen) / 2

prismOct = DislocationLoop(;
    loopType = loopPrism(),
    numSides = 8,
    nodeSide = 1,
    numLoops = 5,
    segLen = segLen * ones(8),
    slipSystemIdx = 1,
    slipSystem = slipSystems,
    label = nodeTypeDln.(ones(Int, 8)),
    buffer = 0,
    range = [0 dx; 0 dy; 0 dz],
    dist = Rand(),
)
shearPent = DislocationLoop(;
    loopType = loopShear(),
    numSides = 5,
    nodeSide = 2,
    numLoops = 5,
    segLen = segLen * ones(10),
    slipSystemIdx = 2,
    slipSystem = slipSystems,
    label = nodeTypeDln.(ones(Int, 10)),
    buffer = 0,
    range = [0 dx; 0 dy; 0 dz],
    dist = Rand(),
)
network = DislocationNetwork((prismOct, shearPent))

dlnFig = plotNodes(
    prismOct,
    m = 2,
    l = 3,
    linecolor = :blue,
    markershape = :circle,
    markercolor = :blue,
    legend = false,
)
plotNodes!(
    dlnFig,
    shearPent,
    m = 2,
    l = 3,
    linecolor = :red,
    markershape = :square,
    markercolor = :red,
    legend = false,
    camera = (10, 20),
)
# savefig(dlnFig, "loops.svg")

networkFig = plotNodes(
    network,
    m = 2,
    l = 3,
    linecolor = :blue,
    markershape = :circle,
    markercolor = :blue,
    legend = false,
)

networkFEFig = plotNodes(
    meshC,
    network,
    m = 2,
    l = 3,
    linecolor = :blue,
    markershape = :circle,
    markercolor = :blue,
    legend = false,
    camera = camera = (-15, 25),
)
# savefig(networkFEFig, "networkPreRem.svg")

# Add more to network
network = DislocationNetwork!(network, prismOct);

## Calculations
# Peach-Köhler force. We don't have a function to apply loading yet so we just generate a random applied stress.
forceDispC.uHat .= rand(-50:50, length(forceDispC.uHat)) .* rand(length(forceDispC.uHat))

fPK = calcPKForce(meshC, forceDispC, network)
fPK69 = calcPKForce(meshC, forceDispC, network, 69)
fPK420 = calcPKForce(meshC, forceDispC, network, [4, 20])
fPK[:, 69] ≈ fPK69
fPK[:, [4, 20]] ≈ fPK420

# Self forces on nodes 1 and 2 of the segments.
sf = calcSelfForce(dlnParams, matParams, network)
sf42 = calcSelfForce(dlnParams, matParams, network, 42)
sf619 = calcSelfForce(dlnParams, matParams, network, [6, 1, 9])
sf[1][:, 42] ≈ sf42[1] && sf[2][:, 42] ≈ sf42[2]
sf[1][:, [6, 1, 9]] ≈ sf619[1] && sf[2][:, [6, 1, 9]] ≈ sf619[2]

# Remote forces.
rf = calcSegSegForce(dlnParams, matParams, network)
rf13 = calcSegSegForce(dlnParams, matParams, network, 13)
rf666 = calcSegSegForce(dlnParams, matParams, network, [66, 6])
rf[:, :, 13] ≈ rf13
rf[:, :, [66, 6]] ≈ rf666

# Total forces.
f = calcSegForce(dlnParams, matParams, meshC, forceDispC, network)
f57 = calcSegForce(dlnParams, matParams, meshC, forceDispC, network, 57)
f316 = calcSegForce(dlnParams, matParams, meshC, forceDispC, network, [31, 6])
f[:, :, 57] ≈ f57
f[:, :, [31, 6]] ≈ f316

calcSegForce!(dlnParams, matParams, meshC, forceDispC, network)
f ≈ network.segForce[:, :, 1:network.numSeg[1]]

# Mobility
# ! This mobility law is outdated.
nodeForce, nodeVel = dlnMobility(dlnParams, matParams, network)
dlnMobility!(dlnParams, matParams, network)
nodeVel ≈ network.nodeVel[:, 1:network.numNode[1]]
nodeForce ≈ network.nodeForce[:, 1:network.numNode[1]]

# Remeshing
# Uncomment following 3 lines to add more loops to network.
# network = DislocationNetwork!(network, (prismOct, shearPent))
# calcSegForce!(dlnParams, matParams, meshC, forceDispC, network)
# dlnMobility!(dlnParams, matParams, network)

# Remesh the surface.
network = remeshSurfaceNetwork!(meshC, boundaryC, network)
# Coarsen the internal network.
network = coarsenNetwork!(dlnParams, matParams, meshC, forceDispC, network)
# Refine the internal network.
network = refineNetwork!(dlnParams, matParams, meshC, forceDispC, network)
remFEFig = plotNodes(
    meshC,
    network,
    m = 2,
    l = 3,
    linecolor = :blue,
    markershape = :circle,
    markercolor = :blue,
    legend = false,
    camera = (-15, 25),
)
idx = findall(x -> x ∈ (3, 4), network.label)
scatter!(
    network.coord[1, idx],
    network.coord[2, idx],
    network.coord[3, idx],
    m = 2,
    markercolor = :red,
)
# savefig(remFEFig, "networkPostRem.svg")

#=
    collision == true if a collision was found false otherwise.
    collisionType == :null if no collision was found,
        :hinge if two segments sharing a node collide,
        :twoLine if two unconnected segments collide.
    n1s1 == node 1 segment 1, n2s1 == node 2 segment 1
    n1s2 == node 1 segment 2, n2s2 == node 2 segment 2
    s1 == segment 1, s2 == segment 2
    L1, L2 == normalised position along the segments where collision will occur =#
collision, collisionType, n1s1, n2s1, n1s2, n2s2, s1, s2, L1, L2 =
    detectCollision(dlnParams, network)

# Dislocation displacements on the displacement nodes.
uTilde = calc_uTilde(meshC, boundaryC, matParams, network)
calc_uTilde!(forceDispC, meshC, boundaryC, matParams, network)
uTilde ≈ forceDispC.uTilde[boundaryC.uDofsDln]

reshape(uTilde, :, 3)' == forceDispC.uTilde[reshape(boundaryC.uDofsDln, :, 3)']
meshC.coord[:, boundaryC.uGammaDln] == meshC.coord[reshape(boundaryC.uDofsDln, :, 3)']

coordDisp = copy(meshC.coord[:, boundaryC.uGammaDln])
coordDisp = coordDisp + forceDispC.uTilde[reshape(boundaryC.uDofsDln, :, 3)']

deformedMesh = deepcopy(meshC)
deformedMesh.coord[reshape(boundaryC.uDofsDln, :, 3)'] +=
    20 * forceDispC.uTilde[reshape(boundaryC.uDofsDln, :, 3)']
figDef = plotBoundaries(boundaryC, deformedMesh; camera = (-20, 20))
plotNodes!(
    figDef,
    meshC,
    network,
    m = 2,
    l = 3,
    linecolor = :blue,
    markershape = :circle,
    markercolor = :blue,
    legend = false,
    camera = (-15, 25),
)

# savefig("dispBound.svg")

scatter!(coordDisp[1, :], coordDisp[2, :], coordDisp[3, :])

meshC.coord[reshape(boundaryC.uDofsDln, :, 3)'] - coordDisp
meshC.connectivity

boundaryC.uGammaDln

# Numeric tractions q = 1
N, A = getTGammaDlnNormsArea(boundaryC, meshC)
points = meshC.coord[:, boundaryC.tGammaDln]
σTilde = calc_σTilde(points, dlnParams, matParams, network)
T = calcNumericTractions(σTilde, N, A, Val(1))

# Gauss quadrature

# Field point stresses
numNode = 102;
numSeg = numNode - 1;
len = numNode + 1;
xrange = range(-100, 100, length = len);
yrange = range(-100, 100, length = len);

X = ones(length(yrange)) .* xrange';
Y = ones(length(xrange))' .* yrange;
Z = zeros(length(xrange))' .* zeros(len);
points = [X[:]'; Y[:]'; Z[:]'];

points2 = Iterators.product(xrange, yrange, 0)
reshape(collect(Iterators.flatten(collect.(points2)')), 3, :)
points3 = reshape(collect(Iterators.flatten(collect.(points2)')), 3, :)

points == points3
X, Y, Z = nothing, nothing, nothing;

a = 5.0
# We need some dummy parameters.
matParams = MaterialParameters(; crystalStruct = BCC(), μ = 1.0, μMag = 1.0, ν = 0.28);
dlnParams = DislocationParameters(;
    coreRad = a,
    coreRadMag = 1.0,
    minSegLen = a + 2,
    maxSegLen = a + 3,
    minArea = a + 1,
    maxArea = a + 2,
    mobility = mobBCC(),
);

# Screw dislocation.

l = Float64[0; 0; 1]
b = Float64[0; 0; 1]
n = b × l

links = zeros(Int, 2, numSeg)
slipPlane = zeros(3, numSeg)
bVec = zeros(3, numSeg)
coord = [zeros(len)'; zeros(len)'; xrange']
label = zeros(nodeTypeDln, len)
nodeVel = similar(coord)
nodeForce = similar(coord)
for i in 1:numSeg
    links[:, i] .= (i, i + 1)
    slipPlane[:, i] = n
    bVec[:, i] = b
end

network = DislocationNetwork(;
    links = links,
    slipPlane = slipPlane,
    bVec = bVec,
    coord = coord,
    label = label,
    nodeVel = nodeVel,
    nodeForce = nodeForce,
)
makeConnect!(network)
getSegmentIdx!(network)

# Out of place calculation.
stress = reshape(calc_σTilde(points, dlnParams, matParams, network), 6, len, :)

σ = zeros(6, size(points, 2))
# In-place calculation.
calc_σTilde!(σ, points, dlnParams, matParams, network)
σ = reshape(σ, 6, len, :)

println("σ ≈ stress : ", σ ≈ stress)

figXZ = contourf(xrange, yrange, σ[5, :, :], levels = 30)
figYZ = contourf(xrange, yrange, σ[6, :, :], levels = 30)

# Edge dislocation.

l = Float64[0; 0; 1]
b = Float64[1; 0; 0]
n = b × l

links = zeros(Int, 2, numSeg)
slipPlane = zeros(3, numSeg)
bVec = zeros(3, numSeg)
coord = [zeros(len)'; zeros(len)'; xrange']
label = zeros(nodeTypeDln, len)
nodeVel = similar(coord)
nodeForce = similar(coord)
for i in 1:numSeg
    links[:, i] .= (i, i + 1)
    slipPlane[:, i] = n
    bVec[:, i] = b
end

network = DislocationNetwork(;
    links = links,
    slipPlane = slipPlane,
    bVec = bVec,
    coord = coord,
    label = label,
    nodeVel = nodeVel,
    nodeForce = nodeForce,
)
makeConnect!(network)
getSegmentIdx!(network)

# Out of place calculation.
stress = reshape(calc_σTilde(points, dlnParams, matParams, network), 6, len, :)

σ = zeros(6, size(points, 2))
# In-place calculation.
calc_σTilde!(σ, points, dlnParams, matParams, network)
σ = reshape(σ, 6, len, :)

println("σ ≈ stress : ", σ ≈ stress)

figXX = contourf(xrange, yrange, σ[1, :, :], levels = 30)
figYY = contourf(xrange, yrange, σ[2, :, :], levels = 30)
figZZ = contourf(xrange, yrange, σ[3, :, :], levels = 30)
figXY = contourf(xrange, yrange, σ[4, :, :], levels = 30)

ext: jl

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
#! format: off
# This file is gonna be a hot mess anyway, no reason to format it
using Test, CSV
include("utils.jl")

@testset "fund_args" begin
    @test @testbench fund_args(t) (t::Range{Float64,0,1e7})
end
@testset "norm_ang" begin
    @test @testbench norm_ang(θ) (θ::Range{Float64,-1e5,1e5})
end
@testset "ee_ct" begin
    @test @testbench ee_ct(jd_high, jd_low; accuracy = :full) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1})
    @test @testbench ee_ct(jd_high, jd_low; accuracy = :reduced) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1})
end

@testset "nu2000k" begin
    @test @testbench nu2000k(jd_high, jd_low) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1})
end
@testset "iau2000a" begin
    @test @testbench iau2000a(jd_high, jd_low) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1})
end

@testset "nutation_angles" begin
    @test @testbench nutation_angles(t; accuracy = :full) (t::Range{Float64,-1e3,4e3})
    @test @testbench nutation_angles(t; accuracy = :reduced) (t::Range{Float64,-1e3,4e3})
end
@testset "mean_obliq" begin
    @test @testbench mean_obliq(jd_tdb) (jd_tdb::Range{Float64,0,1e7})
end
@testset "frame_tie" begin
    @test @testbench frame_tie(pos, :dynamic2icrs) (pos::Position)
    @test @testbench frame_tie(pos, :icrs2dynamic) (pos::Position)
end
@testset "nutation" begin
    @test @testbench nutation(jd_tdb, pos; accuracy = :full, direction = :mean2true) (jd_tdb::Range{Float64,-1e7,1e7}) (pos::Position)
    @test @testbench nutation(jd_tdb, pos; accuracy = :reduced, direction = :mean2true) (jd_tdb::Range{Float64,-1e7,1e7}) (pos::Position)
    @test @testbench nutation(jd_tdb, pos; accuracy = :full, direction = :true2mean) (jd_tdb::Range{Float64,-1e7,1e7}) (pos::Position)
    @test @testbench nutation(jd_tdb, pos; accuracy = :reduced, direction = :true2mean) (jd_tdb::Range{Float64,-1e7,1e7}) (pos::Position)
end
@testset "precession" begin
    @test @testbench precession(jd_tdb, pos, NOVAS.T0) (jd_tdb::Range{Float64,-1e7,1e7}) (pos::Position)
end
@testset "e_tilt" begin
    @test @testbench e_tilt(jd_tdb; accuracy = :full) (jd_tdb::Range{Float64,-1e7,1e7})
    @test @testbench e_tilt(jd_tdb; accuracy = :reduced) (jd_tdb::Range{Float64,-1e7,1e7})
end
@testset "ira_equinox" begin
    @test @testbench ira_equinox(jd_tdb; accuracy = :full, equinox = :mean) (jd_tdb::Range{Float64,-1e7,1e7})
    @test @testbench ira_equinox(jd_tdb; accuracy = :full, equinox = :true) (jd_tdb::Range{Float64,-1e7,1e7})
    @test @testbench ira_equinox(jd_tdb; accuracy = :reduced, equinox = :mean) (jd_tdb::Range{Float64,-1e7,1e7})
    @test @testbench ira_equinox(jd_tdb; accuracy = :reduced, equinox = :true) (jd_tdb::Range{Float64,-1e7,1e7})
end
@testset "cio_location" begin
    @test @testbench cio_location(jd_tdb; accuracy = :full) (jd_tdb::Range{Float64,-1e7,1e7})
end
@testset "cio_basis" begin
    @test @testbench cio_basis(jd_tdb, ra_cio; accuracy = :full) (jd_tdb::Range{Float64,-1e7,1e7}) (ra_cio::Range{Float64,0,24})
    @test @testbench cio_basis(jd_tdb, ra_cio; accuracy = :reduced) (jd_tdb::Range{Float64,-1e7,1e7}) (ra_cio::Range{Float64,0,24})
end
@testset "era" begin
    @test @testbench era(jd_high, jd_low) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1})
end
@testset "tdb2tt" begin
    @test @testbench tdb2tt(jdb_jd) (jdb_jd::Range{Float64,-1e7,1e7})
end
@testset "sidereal_time" begin
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :mean, method = :CIO, accuracy = :full) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :mean, method = :CIO, accuracy = :reduced) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :mean, method = :equinox, accuracy = :full) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :mean, method = :equinox, accuracy = :reduced) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :apparent, method = :CIO, accuracy = :full) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :apparent, method = :CIO, accuracy = :reduced) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :apparent, method = :equinox, accuracy = :full) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
    @test @testbench sidereal_time(jd_high, jd_low, delta_t; gst_type = :apparent, method = :equinox, accuracy = :reduced) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70})
end
@testset "wobble" begin
    @test @testbench wobble(tjd, xp, yp, pos; direction = :itrs2terr) (tjd::Range{Float64,-1e7,1e7}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
    @test @testbench wobble(tjd, xp, yp, pos; direction = :terr2itrs) (tjd::Range{Float64,-1e7,1e7}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
end
@testset "spin" begin
    @test @testbench spin(θ, pos) (θ::Range{Float64,0,360}) (pos::Position)
end
@testset "ter2cel" begin
    @test @testbench ter2cel(jd_high, jd_low, delta_t, pos; method = :CIO, accuracy = :full, option = :GCRS, xp = xp, yp = yp) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
    @test @testbench ter2cel(jd_high, jd_low, delta_t, pos; method = :equinox, accuracy = :full, option = :GCRS, xp = xp, yp = yp) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
    @test @testbench ter2cel(jd_high, jd_low, delta_t, pos; method = :CIO, accuracy = :reduced, option = :GCRS, xp = xp, yp = yp) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
    @test @testbench ter2cel(jd_high, jd_low, delta_t, pos; method = :equinox, accuracy = :reduced, option = :GCRS, xp = xp, yp = yp) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
    # Option only matters for method = :equinox
    @test @testbench ter2cel(jd_high, jd_low, delta_t, pos; method = :equinox, accuracy = :full, option = :equinox, xp = xp, yp = yp) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
    @test @testbench ter2cel(jd_high, jd_low, delta_t, pos; method = :equinox, accuracy = :reduced, option = :equinox, xp = xp, yp = yp) (jd_high::Range{Float64,0,1e7}) (jd_low::Range{Float64,0,1}) (delta_t::Range{Float64,-10,70}) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1}) (pos::Position)
end
@testset "refract" begin
    @test @testbench refract(location, zd; ref_option = :standard) (location::Location) (zd::Range{Float64,0,90})
    @test @testbench refract(location, zd; ref_option = :location) (location::Location) (zd::Range{Float64,0,90})
end
@testset "equ2hor" begin
    @test @testbench equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy = :full, xp = xp, yp = yp, ref_option = :none) (jd_ut1::Range{Float64,0,1e7}) (delta_t::Range{Float64,-10,70}) (ra::Range{Float64,0,24}) (dec::Range{Float64,0,360}) (location::Location) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1})
    @test @testbench equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy = :full, xp = xp, yp = yp, ref_option = :standard) (jd_ut1::Range{Float64,0,1e7}) (delta_t::Range{Float64,-10,70}) (ra::Range{Float64,0,24}) (dec::Range{Float64,0,360}) (location::Location) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1})
    @test @testbench equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy = :full, xp = xp, yp = yp, ref_option = :location) (jd_ut1::Range{Float64,0,1e7}) (delta_t::Range{Float64,-10,70}) (ra::Range{Float64,0,24}) (dec::Range{Float64,0,360}) (location::Location) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1})
    @test @testbench equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy = :reduced, xp = xp, yp = yp, ref_option = :none) (jd_ut1::Range{Float64,0,1e7}) (delta_t::Range{Float64,-10,70}) (ra::Range{Float64,0,24}) (dec::Range{Float64,0,360}) (location::Location) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1})
    @test @testbench equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy = :reduced, xp = xp, yp = yp, ref_option = :standard) (jd_ut1::Range{Float64,0,1e7}) (delta_t::Range{Float64,-10,70}) (ra::Range{Float64,0,24}) (dec::Range{Float64,0,360}) (location::Location) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1})
    @test @testbench equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy = :reduced, xp = xp, yp = yp, ref_option = :location) (jd_ut1::Range{Float64,0,1e7}) (delta_t::Range{Float64,-10,70}) (ra::Range{Float64,0,24}) (dec::Range{Float64,0,360}) (location::Location) (xp::Range{Float64,-1,1}) (yp::Range{Float64,-1,1})
end

# Serialize Results
CSV.write("../benchmarks.csv", results)


text:
using Gen
using AdvancedHMC
using Logging
using LoggingExtras

function disable_sample_logging()
    function ignore_sampling_filter(log_args)
        !(occursin("sampling steps",log_args.message) || occursin("adapation steps",log_args.message))
    end

    logger = ActiveFilteredLogger(ignore_sampling_filter, global_logger())

    if !(@isdefined old_logger) #do this only once
        global old_logger = global_logger(logger)
    end
end

function my_nuts(trace, selection,
        transformations,
        n_postadapt_steps = 2,
        n_adapts = 1,
        initial_ϵ_reduce_fac = 10)

    n_NUTS_steps = n_postadapt_steps + n_adapts

    filtered_choices = get_selected(get_choices(trace), selection)
    cur_xy = to_array_transformed(filtered_choices, transformations, 1, Float64)
    dimension = length(cur_xy)
    metric = DiagEuclideanMetric(dimension)

    retval_grad = nothing #accepts_output_grad(get_gen_fn(trace)) ? zero(get_retval(trace)) : nothing

    function update_xy(val)
        extra_constraints, grad_val = from_array_transformed(filtered_choices, transformations, val)
        (new_trace, weight, discard, retdiff) = update(trace, get_args(trace),
                                                    map((_)->NoChange(), get_args(trace)),
                                                    extra_constraints)
        (new_trace, weight, grad_val)
    end

    function val_to_lp_plus_c(val)
        (new_trace, weight, grad_val) = update_xy(val)
        weight
    end

    function val_to_grad(val)
        (new_trace, weight, grad_val) = update_xy(val)
        (retval_grad_out, values_trie, gradient_trie) = choice_gradients(new_trace, selection, retval_grad)
        raw_grad = to_array(gradient_trie, Float64)
        grad = raw_grad .* grad_val #chain rule
        (weight, grad)
    end

    # Define a Hamiltonian system, using metric defined globally above
    hamiltonian = Hamiltonian(metric, val_to_lp_plus_c, val_to_grad)

    # Define a leapfrog solver, with initial step size chosen heuristically
    initial_ϵ = find_good_stepsize(hamiltonian, cur_xy) ./ initial_ϵ_reduce_fac
    integrator = Leapfrog(initial_ϵ)

    # Define an HMC sampler, with the following components
    #   - multinomial sampling scheme,
    #   - generalised No-U-Turn criteria, and
    #   - windowed adaption for step-size and diagonal mass matrix
    proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

    # Run the sampler to draw samples from the specified Gaussian, where
    #   - `samples` will store the samples
    #   - `stats` will store diagnostic statistics for each sample

    samples, stats = sample(hamiltonian, proposal, cur_xy, n_NUTS_steps, adaptor, n_adapts; progress=false)

    #println(samples[3])

    (new_trace, weight, grad_val) = update_xy(samples[n_NUTS_steps])
    new_trace
end

transform_log = (log, exp, exp)
    #(restricted->unrestricted,
        #unrestricted->restricted,
        #unrestricted->d(restricted)/d(unrestricted))

function to_array_transformed(choices::ChoiceMap, transformations::ChoiceMap, whichtrans, ::Type{T}) where {T}
    arr = Vector{T}(undef, 32)
    n = _fill_array_transformed!(choices, transformations, whichtrans, arr, 1)
    @assert n <= length(arr)
    resize!(arr, n)
    arr
end

function _fill_array_transformed!(choices::DynamicChoiceMap, transformations::ChoiceMap, whichtrans,
                arr::Vector{T}, start_idx::Int) where {T}
    leaf_keys_sorted = sort(collect(keys(choices.leaf_nodes)))
    internal_node_keys_sorted = sort(collect(keys(choices.internal_nodes)))
    idx = start_idx
    for key in leaf_keys_sorted
        value = choices.leaf_nodes[key]
        if has_value(transformations,key)
            n_written = _fill_array_transformed!(value, transformations[key][whichtrans], arr, idx)
        else
            n_written = Gen._fill_array!(value, arr, idx)
        end
        idx += n_written
    end
    for key in internal_node_keys_sorted
        if key in keys(transformations.internal_nodes)
            n_written = _fill_array_transformed!(get_submap(choices, key), get_submap(transformations, key), whichtrans, arr, idx)
        else
            n_written = Gen._fill_array!(get_submap(choices, key), arr, idx)
        end
        idx += n_written
    end
    idx - start_idx
end

function _fill_array_transformed!(value::T, transformation, arr::Vector{T}, start_idx::Int) where {T}
    if length(arr) < start_idx
        resize!(arr, 2 * start_idx)
    end
    arr[start_idx] = transformation(value)
    1
end

function _fill_array_transformed!(value::Vector{T}, transformation, arr::Vector{T}, start_idx::Int) where {T}
    if length(arr) < start_idx + length(value)
        resize!(arr, 2 * (start_idx + length(value)))
    end
    arr[start_idx:start_idx+length(value)-1] = transformation.(value)
    length(value)
end

function from_array_transformed(proto_choices::ChoiceMap, transformations::ChoiceMap, arr::Vector)
    grad_arr = copy(arr)
    (n, choices) = _from_array_transformed!(proto_choices, transformations, grad_arr, 1)
    if n != length(arr)
        error("Dimension mismatch: $n, $(length(arr))")
    end
    choices, grad_arr
end

function applyTransformation(transformation, value::T) where {T}
    transformation(value)
end
function applyTransformation(transformation, value::Vector{T}) where {T}
    transformation.(value)
end

function _from_array_transformed!(proto_choices::DynamicChoiceMap, transformations::ChoiceMap, arr::Vector{T}, start_idx::Int) where {T}
    @assert length(arr) >= start_idx
    choices = DynamicChoiceMap()
    leaf_keys_sorted = sort(collect(keys(proto_choices.leaf_nodes)))
    internal_node_keys_sorted = sort(collect(keys(proto_choices.internal_nodes)))
    idx = start_idx
    for key in leaf_keys_sorted
        (n_read, value) = Gen._from_array(proto_choices.leaf_nodes[key], arr, idx)
        if has_value(transformations,key)
            transformedValue = applyTransformation(transformations[key][2],value) #restricted to unrestricted
            choices.leaf_nodes[key] = transformedValue
            if transformations[key][2] != transformations[key][3] #not exp, so grad != inverse
                transformedValue = applyTransformation(transformations[key][3],value) #restricted to grad unrestricted
            end
            if n_read == 1
                arr[idx] = transformedValue
            else
                arr[idx:idx+n_read-1] = transformedValue
            end
        else
            choices.leaf_nodes[key] = value
            arr[idx:idx+n_read-1] = ones(T, n_read) #identity transform has a grad of 1
        end
        idx += n_read
    end
    for key in internal_node_keys_sorted
        if key in keys(transformations.internal_nodes)
            (n_read, node) = _from_array_transformed!(get_submap(proto_choices, key),
                                                get_submap(transformations, key), arr, idx)
        else
            (n_read, node) = Gen._from_array(get_submap(proto_choices, key), arr, idx)
            arr[idx:idx+n_read-1] = ones(T, n_read) #identity transform has a grad of 1
        end
        idx += n_read
        choices.internal_nodes[key] = node
    end
    (idx - start_idx, choices)
end

max_stars_repo_name: jamesonquinn/Gen

max_issues_repo_name: jamesonquinn/Gen

max_forks_repo_name: jamesonquinn/Gen


text:
using Docile, Lexicon, LMDB

const api_directory = "api"
const modules = [LMDB]

cd(dirname(@__FILE__)) do

    # Run the doctests *before* we start to generate *any* documentation.
    for m in modules
        failures = failed(doctest(m))
        if !isempty(failures.results)
            println("\nDoctests failed, aborting commit.\n")
            display(failures)
            exit(1) # Bail when doctests fail.
        end
    end

    # Generate and save the contents of docstrings as markdown files.
    index  = Index()
    for mod in modules
        update!(index, save(joinpath(api_directory, "$(mod).md"), mod))
    end
    save(joinpath(api_directory, "index.md"), index; md_subheader = :category)

    # Add a reminder not to edit the generated files.
    open(joinpath(api_directory, "README.md"), "w") do f
        print(f, """
        Files in this directory are generated using the `build.jl` script. Make
        all changes to the originating docstrings/files rather than these ones.

        Documentation should *only* be build directly on the `master` branch.
        Source links would otherwise become unavailable should a branch be
        deleted from the `origin`. This means potential pull request authors
        *should not* run the build script when filing a PR.
        """)
    end

    info("Adding all documentation changes in $(api_directory) to this commit.")
    success(`git add $(api_directory)`) || exit(1)

end


text:
using Dates

using TeaFiles.Header: TimeSection, _duration_div, field_type, get_section,
    get_primary_time_field, get_time_fields

@testset "core" begin
    @testset "get_section" begin
        metadata = _get_example_metadata()
        time_section = get_section(metadata, TimeSection)
        @test isa(time_section, TimeSection)
        @test time_section.epoch == 719162
    end

    @testset "get_time_fields" begin
        metadata = _get_example_metadata()
        time_fields = get_time_fields(metadata)
        @test length(time_fields) == 1
        time_field = only(time_fields)
        @test time_field == get_primary_time_field(metadata)

        @test field_type(time_field) == Int64
        @test time_field.offset == 0
    end

    @testset "div_duration" begin
        @test _duration_div(Day(1), Hour(1)) == 24
        @test _duration_div(Day(1), Hour(2)) == 12
        @test _duration_div(Day(1), Minute(1)) == 24 * 60
        @test _duration_div(Day(1), Minute(10)) == 24 * 6
        @test _duration_div(Hour(2), Minute(2)) == 60
        @test_throws ArgumentError _duration_div(Minute(3), Minute(4))
    end

    @testset "invalid metadata" begin
        # An invalid item_start.
        @test_throws ArgumentError TeaFileMetadata(1001, 0, [])

        # An invalid item_end.
        @test_throws ArgumentError TeaFileMetadata(1000, 1, [])
    end
end


text:
export SWAP, FLIP
export ROT, IROT
export INC, DEC

"""
    NoGrad{T} <: IWrapper{T}
    NoGrad(x)

A `NoGrad(x)` is equivalent to `GVar^{-1}(x)`, which cancels the `GVar` wrapper.
"""
@pure_wrapper NoGrad

const NullType{T} = Union{NoGrad{T}, Partial{T}}

@selfdual -

@inline FLIP(b::Bool) = !b
@selfdual FLIP

"""
    INC(a!) -> a! + 1
"""
@inline function INC(a!::Number)
    a! + one(a!)
end

"""
    DEC(a!) -> a! - 1
"""
@inline function DEC(a!::Number)
    a! - one(a!)
end
@dual INC DEC

"""
    SWAP(a!, b!) -> b!, a!
"""
@inline function SWAP(a!::T, b!::T) where T
    b!, a!
end
@selfdual SWAP

"""
    ROT(a!, b!, θ) -> a!', b!', θ

```math
\\begin{align}
    {\\rm ROT}(a!, b!, \\theta)  = \\begin{bmatrix}
        \\cos(\\theta) & - \\sin(\\theta)\\\\
        \\sin(\\theta)  & \\cos(\\theta)
    \\end{bmatrix}
    \\begin{bmatrix}
        a!\\\\
        b!
    \\end{bmatrix},
\\end{align}
```
"""
@inline function ROT(i::Real, j::Real, θ::Real)
    a, b = rot(i, j, θ)
    a, b, θ
end

"""
    IROT(a!, b!, θ) -> ROT(a!, b!, -θ)
"""
@inline function IROT(i::Real, j::Real, θ::Real)
    i, j, _ = ROT(i, j, -θ)
    i, j, θ
end
@dual ROT IROT

for F1 in [:(Base.:-)]
    @eval @inline function $F1(a!::NullType)
        @instr $F1(a! |> value)
        a!
    end
end

for F2 in [:SWAP, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    @eval @inline function $F2(a::NullType, b::Real)
        @instr $(NiLangCore.get_argname(F2))(a |> value, b)
        a, b
    end
    @eval @inline function $F2(a::NullType, b::NullType)
        @instr $(NiLangCore.get_argname(F2))(a |> value, b |> value)
        a, b
    end
    @eval @inline function $F2(a::Real, b::NullType)
        @instr $(NiLangCore.get_argname(F2))(a, b |> value)
        a, b
    end
end

function type_except(::Type{TT}, ::Type{T2}) where {TT, T2}
    N = length(TT.parameters)
    setdiff(Base.Iterators.product(zip(TT.parameters, repeat([T2], N))...), [ntuple(x->T2, N)])
end

for F3 in [:ROT, :IROT, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    PS = (:a, :b, :c)
    for PTS in type_except(Tuple{NullType, NullType, NullType}, Real)
        params = map((P,PT)->PT <: NullType ? :($P |> value) : P, PS, PTS)
        params_ts = map((P,PT)->:($P::$PT), PS, PTS)
        @eval @inline function $F3($(params_ts...))
            @instr $F3($(params...))
            ($(PS...),)
        end
    end
end

# patch for fixed point numbers
function (f::PlusEq{typeof(/)})(out!::T, x::Integer, y::Integer) where T<:Fixed
    out!+T(x)/y, x, y
end

function (f::MinusEq{typeof(/)})(out!::T, x::Integer, y::Integer) where T<:Fixed
    out!-T(x)/y, x, y
end

for F in [:exp, :log, :sin, :cos]
    @eval Base.$F(x::Fixed43) = Fixed43($F(Float64(x)))
    @eval (f::PlusEq{typeof($F)})(out!::Fixed43, x::Real) = out! + Fixed43($F(x)), x
    @eval (f::MinusEq{typeof($F)})(out!::Fixed43, x::Real) = out! - Fixed43($F(x)), x
end

Base.:^(x::Integer, y::Fixed43) = Fixed43(x^(Float64(y)))
Base.:^(x::Fixed43, y::Fixed43) = Fixed43(x^(Float64(y)))
Base.:^(x::T, y::Fixed43) where T<:AbstractFloat = x^(T(y))

function (::PlusEq{typeof(convert)})(out!::T, y) where T<:Real
    out! + convert(T, y), y
end

function (::MinusEq{typeof(convert)})(out!::T, y) where T<:Real
    out! - convert(T, y), y
end


text:
function splits(text::AbstractString, obj::HiddenMarkovModel; ntrial = 2, pos = true)
    codes = codeunits(text)
    aca = obj.aca; pos = obj.pos;
    hpr = obj.hpr; h2h = obj.h2h;
    alpha = obj.alpha; h2v = obj.h2v;
    positions = collect(eachmatch(aca, text))
    if (length(positions) == 0) return [String(text)] end
    ncode = length(codes); npos = length(pos);
    dp = Array{Float64, 3}(undef, ntrial, npos, ncode + 1); dp .= 1.0 / 0.0; dp[1, 1, ncode + 1] = 0.0;
    nx = Array{Tuple{Int, Int, Int}, 3}(undef, ntrial, npos, ncode + 1)
    j = length(positions)
    charindexes = collect(eachindex(text)); push!(charindexes, ncode + 1);
    k = nchar = length(charindexes) - 1
    for i = ncode:-1:1
        while j > 0 && i < positions[j].t j -= 1 end
        while j > 0 && i == positions[j].t
            iword = positions[j].i
            s = positions[j].s
            t = positions[j].t
            for ipos = 1:npos
                res = view(dp, :, ipos, s)
                H2V = get(h2v[ipos], iword, alpha[ipos])
                nxt = view(nx, :, ipos, s)
                for ipos2 = 1:npos
                    for itrial2 = 1:ntrial
                        cur = H2V + h2h[ipos2, ipos] + dp[itrial2, ipos2, t + 1]
                        for itrial = 1:ntrial
                            if cur < res[itrial]
                                res[itrial + 1:end] .= res[itrial:end - 1]
                                nxt[itrial + 1:end] .= nxt[itrial:end - 1]
                                res[itrial] = cur
                                nxt[itrial] = (itrial2, ipos2, t)
                                # @show itrial, ipos, s, itrial2, ipos2, t + 1
                                break
                            end
                        end
                    end
                end
            end
            j -= 1
        end
        if i == charindexes[k]
            iword = 0
            s = i
            t = charindexes[k + 1] - 1
            for ipos = 1:npos
                res = view(dp, :, ipos, s)
                H2V = get(h2v[ipos], iword, alpha[ipos])
                nxt = view(nx, :, ipos, s)
                for ipos2 = 1:npos
                    for itrial2 = 1:ntrial
                        cur = H2V + h2h[ipos2, ipos] + dp[itrial2, ipos2, t + 1]
                        for itrial = 1:ntrial
                            if cur < res[itrial]
                                res[itrial + 1:end] .= res[itrial:end - 1]
                                nxt[itrial + 1:end] .= nxt[itrial:end - 1]
                                res[itrial] = cur
                                nxt[itrial] = (itrial2, ipos2, t)
                                break
                            end
                        end
                    end
                end
            end
            k -= 1
        end
    end
    # dp[:, :, 1] .+= hpr
    for i = 1:ntrial dp[i, :, 1] .+= hpr end
    res = Vector()
    dp1 = dp[:, :, 1]
    for i = 1:ntrial
        minind = findmin(dp1)[2]
        dp1[minind] = 1.0 / 0.0
        itrial = minind[1]
        ipos = minind[2]
        pvis = 1
        segs = Vector{Tuple{String, String}}()
        while pvis <= ncode
            # @show nx[ipos, pvis]
            nxt = nx[itrial, ipos, pvis]
            itrial2 = nxt[1]
            ipos2 = nxt[2]
            pvis2 = nxt[3]
            push!(segs, (pos[ipos], String(codes[pvis:pvis2])))
            itrial = itrial2
            ipos = ipos2
            pvis = pvis2 + 1
        end
        push!(res, segs)
    end
    push!(res, dp[:, :, 1])
    push!(res, nx[:, :, 1])
    return res
end


text:
##
using IFGF
using LinearAlgebra
using StaticArrays
using Random
using Plots
plotlyjs()
Random.seed!(1)

includet(joinpath(IFGF.PROJECT_ROOT,"test","testutils.jl"))

# parameters
k        = 4π
λ        = 2π/k
p        = (4,4,4)
splitter = Trees.CardinalitySplitter(;nmax=200)
ds_func  = IFGF.cone_domain_size_func(k,(1.0,1.0,1.0))
K    = HelmholtzKernel(k)
T    = ComplexF64

# specify the function to generate the point clouds
pts_func = sphere_uniform_surface_mesh

dxvec = [0.1/2^n for n in 0:5]
# loop number of points
er   = []
ta   = [] # assemble time
tp   = [] # product time
npts = []
for dx in dxvec
    Xpts = Ypts = pts_func(dx)
    n    = length(Xpts)
    push!(npts,n)
    I   = rand(1:n,100)
    B   = randn(T,n)
    C   = zeros(T,n)
    exa = [sum(K( Xpts[i],Ypts[j])*B[j] for j in 1:n) for i in I]
    tmp = @elapsed A  = IFGFOp(K,Xpts,Ypts;splitter,p,ds_func,threads=false)
    push!(ta,tmp)
    tmp = @elapsed mul!(C,A,B,1,0;threads=false)
    push!(tp,tmp)
    ee = norm(C[I]-exa,2) / norm(exa,2)
    push!(er,ee)
end

fig = plot(npts[2:end],ta[2:end],m=:circle,
    ylabel="time (s)",
    xlabel="#X",
    yscale=:log10,
    xscale=:log10,
    label="precomputation")
plot!(npts[2:end],tp[2:end],m=:cross,label="forward map")
# add a reference slope
p1 = 1
p2 = 1
yy = @. npts[2:end]^p1 * log(4,npts[2:end])^p2*tp[end]/(npts[end]^p1*log(4,npts[end])^p2)
plot!(npts[2:end],yy,ls=:dash,label="loglinear slope")

dir = @__DIR__
savefig(fig,joinpath(dir,"sphere_n_scaling.png"))

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
import Random.rand, PyPlot

function main()
    L::Vector{Int64} = [0 for ind = 1:1000]
    for ind = 1:10000000
        r = rand()
        indice::Int64 = Int64(ceil(1000*r))
        L[indice]+=1
    end
    return L
end

function affichage(L::Vector{Int64})
    PyPlot.plot([ind for ind = 1:length(L)], L, "k-")
    PyPlot.ylim(0, 30000)
    PyPlot.show()
end


text:
importall Base

#abstract Algorithm

type Lbfgs <: Algorithm
  init_alpha::Float64
  tol_obj::Float64
  tol_grad::Float64
  tol_param::Float64
  history_size::Int64
end
Lbfgs(;init_alpha::Number=0.001, tol_obj::Number=1e-8,
  tol_grad::Number=1e-8, tol_param::Number=1e-8, history_size::Number=5) =
    Lbfgs(init_alpha, tol_obj, tol_grad, tol_param, history_size)

type Bfgs <: Algorithm
  init_alpha::Float64
  tol_obj::Float64
  tol_grad::Float64
  tol_param::Float64
end
Bfgs(;init_alpha::Number=0.001, tol_obj::Number=1e-8,
  tol_grad::Number=1e-8, tol_param::Number=1e-8) =
    Bfgs(init_alpha, tol_obj, tol_grad, tol_param)

type Newton <: Algorithm
end

type Optimize <: Methods
  method::Algorithm
  iter::Int64
  save_iterations::Bool
end
Optimize(;method::Algorithm=Lbfgs(), iter::Number=2000, save_iterations::Bool=false) =
  Optimize(method, iter, save_iterations)
Optimize(method::Algorithm) = Optimize(method, 2000, false)

function optimize_show(io::IO, o::Optimize, compact::Bool)
  if compact
    println("Optimize($(o.method), $(o.iter), $(o.save_iterations))")
  else
    println("  method =                  Optimize()")
    if isa(o.method, Lbfgs)
      println("    algorithm =               Lbfgs()")
      println("      init_alpha =              ", o.method.init_alpha)
      println("      tol_obj =                 ", o.method.tol_obj)
      println("      tol_grad =                ", o.method.tol_grad)
      println("      tol_param =               ", o.method.tol_param)
      println("      history_size =            ", o.method.history_size)
    elseif isa(o.method, Bfgs)
      println("    algorithm =               Bfgs()")
      println("      init_alpha =              ", o.method.init_alpha)
      println("      tol_obj =                 ", o.method.tol_obj)
      println("      tol_grad =                ", o.method.tol_grad)
      println("      tol_param =               ", o.method.tol_param)
    else
      println("    algorithm =               Newton()")
    end
    println("    iterations =              ", o.iter)
    println("    save_iterations =         ", o.save_iterations)
  end
end

show(io::IO, o::Optimize) = optimize_show(io, o, false)
showcompact(io::IO, o::Optimize) = optimize_show(io, o, true)


text:
"""
    @latexify expression

Create `LaTeXString` representing `expression`.
Variables and expressions can be interpolated with `\$`.
Keyword arguments can be supplied to `latexify` by appending to the argument.

# Examples
```julia-repl
julia> @latexify x^2 + 3/2
L"\$x^{2} + \\frac{3}{2}\$"

julia> @latexify x^2 + \$(3/2)
L"\$x^{2} + 1.5\$"

julia> @latexify x^2 + 3/2 env=:raw
L"x^{2} + \\frac{3}{2}"
```

See also [`latexify`](@ref), [`@latexrun`](@ref), [`@latexdefine`](@ref).
"""
macro latexify(expr, kwargs...)
    return esc(
        Expr(
            :call, :latexify, Expr(:parameters, _extractparam.(kwargs)...), Meta.quot(expr)
        ),
    )
end

"""
    @latexrun expression

Latexify and evaluate `expression`. Useful for expressions with side effects, like assignments.

# Examples
```julia-repl
julia> @latexrun y = 3/2 + \$(3/2)
L"\$y = \\frac{3}{2} + 1.5\$"

julia> y
3.0
```
See also [`@latexify`](@ref), [`@latexdefine`](@ref).
"""
macro latexrun(expr, kwargs...)
    return esc(
        Expr(
            :block,
            _executable(expr),
            Expr(
                :call,
                :latexify,
                Expr(:parameters, _extractparam.(kwargs)...),
                Meta.quot(expr),
            ),
        ),
    )
end

"""
    @latexdefine expression

Latexify `expression`, followed by an equals sign and the return value of its evaluation.
Any side effects of the expression, like assignments, are evaluated as well.

# Examples
```julia-repl
julia> @latexdefine y = 3/2 + \$(3/2) env=:equation
L"\\begin{equation}
y = \\frac{3}{2} + 1.5 = 3.0
\\end{equation}
"

julia> y
3.0
```
See also [`@latexify`](@ref), [`@latexrun`](@ref).
"""
macro latexdefine(expr, kwargs...)
    return esc(
        Expr(
            :call,
            :latexify,
            Expr(:parameters, _extractparam.(kwargs)...),
            Expr(:call, :Expr, QuoteNode(:(=)), Meta.quot(expr), _executable(expr)),
        ),
    )
end

function _executable(expr)
    return postwalk(expr) do ex
        if ex isa Expr && ex.head == :$
            return ex.args[1]
        end
        return ex
    end
end

_extractparam(arg::Symbol) = arg
_extractparam(arg::Expr) = Expr(:kw, arg.args[1], arg.args[2])


text:
module TentacleEffect

using ..Ahorn, Maple

placements = Maple.TentacleEffect

function Ahorn.canFgBg(effect::Maple.TentacleEffect)
    return true, true
end

function Ahorn.editingOptions(effect::Maple.TentacleEffect)
    return Dict{String, Any}(
        "side" => Maple.tentacle_effect_directions
    )
end

end


text:
"""
    SGLD(n_iters::Int, step_size::Float64)

 Stochastic Gradient Langevin Dynamics sampler.

Usage:

```julia
SGLD(1000, 0.5)
```

Example:

```julia
@model example begin
  ...
end

sample(example, SGLD(1000, 0.5))
```
"""
mutable struct SGLD{T} <: Hamiltonian
    n_iters::Int       # number of samples
    step_size::Float64   # constant scale factor of learning rate
    space::Set{T}    # sampling space, emtpy means all
    gid::Int
end
SGLD(step_size::Float64, space...) = SGLD(1, step_size, space..., 0)
SGLD(n_iters, step_size) = SGLD(n_iters, step_size, Set(), 0)
function SGLD(n_iters, step_size, space...)
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return SGLD(n_iters, step_size, _space, 0)
end
SGLD(alg::SGLD, new_gid::Int) = SGLD(alg.n_iters, alg.step_size, alg.space, new_gid)

function step(model, spl::Sampler{<:SGLD}, vi::VarInfo, is_first::Bool)
    if is_first
        if ~haskey(spl.info, :wum)
            spl.alg.gid != 0 && link!(vi, spl)

            D = length(vi[spl])
            ve = VarEstimator{Float64}(0, zeros(D), zeros(D))
            wum = WarmUpManager(1, Dict(), ve)
            wum[:ϵ] = [spl.alg.step_size]
            wum[:stds] = ones(D)
            spl.info[:wum] = wum

            # Initialize iteration counter
            spl.info[:t] = 0

            spl.alg.gid != 0 && invlink!(vi, spl)
        end

        push!(spl.info[:accept_his], true)
    else
        # Update iteration counter
        spl.info[:t] += 1

        @debug "compute current step size..."
        γ = .35
        ϵ_t = spl.alg.step_size / spl.info[:t]^γ # NOTE: Choose γ=.55 in paper
        push!(spl.info[:wum][:ϵ], ϵ_t)

        @debug "X-> R..."
        if spl.alg.gid != 0
            link!(vi, spl)
            runmodel!(model, vi, spl)
        end

        @debug "recording old variables..."
        θ = vi[spl]
        _, grad = gradient(θ, vi, model, spl)
        verifygrad(grad)

        @debug "update latent variables..."
        θ .-= ϵ_t .* grad ./ 2 .+ rand.(Normal.(zeros(length(θ)), sqrt(ϵ_t)))

        @debug "always accept..."
        push!(spl.info[:accept_his], true)
        vi[spl] = θ

        @debug "R -> X..."
        spl.alg.gid != 0 && invlink!(vi, spl)
    end
    return vi
end


text:
export Scheffer

struct Scheffer <: LakeModel end

@with_kw_noshow mutable struct SchefferParameters <: LakeParameters
    @deftype Float64
    nutrients = 0.7 # current nutrient level
    ib = 2e-5 # g⋅m⁻²⋅day⁻¹ immigration rate of bream (9e-3 g⋅m⁻²⋅year⁻¹)
    ip = 2e-5 # g⋅m⁻²⋅day⁻¹ immigration rate of pike (9e-3 g⋅m⁻²⋅year⁻¹)
    r = 7.5e-3 # day⁻¹ maximum growth rate of bream (2.74 year⁻¹)
    H₁ = 0.5 # half saturation constant
    H₂ = 10 # % half saturation constant
    H₃ = 20 # g⋅m⁻² half saturation constant
    H₄ = 15 # g⋅m⁻² half saturation constant
    cb = 7.5e-5 # m⁻²⋅g⁻¹⋅day⁻¹ intraspecific competition constant for bream (0.0274 m⁻²⋅g⁻¹⋅year⁻¹)
    cp = 2.75e-4 # m⁻²⋅g⁻¹⋅day⁻¹ intraspecific competition constant for pike (0.1 m⁻²⋅g⁻¹⋅year⁻¹)
    prmax = 5e-2 # day⁻¹ maximum predation rate of pike (18.25 year⁻¹)
    ce = 0.1 # pike food conversion efficiency to growth
    mp = 2.25e-3 # day⁻¹ mortality rate of pike (0.82 year⁻¹)
    K = 100 # % maximum vegetation coverage
end

function lake_dynamics!(du, u, p::SchefferParameters, t)
    #TODO: New intervention dynamics not considered for Scheffer model
    B, P = u # g⋅m⁻² bream/pike density

    V = p.K * (p.H₃^2 / (p.H₃^2 + B^2)) # % of lake covered vegetation
    FR = B^2 / (B^2 + p.H₄^2) # fuctional response of pike

    du[1] =
        dB =
            p.ib + p.r * (p.nutrients / (p.nutrients + p.H₁)) * B - p.cb * B^2 -
            p.prmax * FR * P
    du[2] = dP = p.ip + p.ce * p.prmax * FR * P * (V / (V + p.H₂)) - p.mp * P - p.cp * P^2
end

function lake_initial_state(
    nutrients::Float64,
    bream::Float64,
    pike::Float64,
    ::Type{Scheffer};
    kwargs...,
)
    ([bream, pike], SchefferParameters(; nutrients = nutrients, kwargs...))
end

LakeParameters(::Type{Scheffer}, nutr; kwargs...) = SchefferParameters(; nutrients = nutr, kwargs...)

preset_conditions(::Type{Clear}, ::Type{Scheffer}) = (0.7, [20.0, 1.8])
preset_conditions(::Type{Turbid}, ::Type{Scheffer}) = (2.5, [84.0, 0.04])

function Base.show(io::IO, ::MIME"text/plain", p::SchefferParameters)
    println(io, "Parameters for lake dynamics (Scheffer) model):")
    println(io, "Nutrient level: $(p.nutrients)")
    println(io, "Immigration rate (g⋅m⁻²⋅day⁻¹) for bream: $(p.ib), pike: $(p.ip)")
    println(
        io,
        "Growth rate (day⁻¹) for bream: $(p.r), Predation rate (day⁻¹) of pike: $(p.prmax)",
    )
    println(
        io,
        "Half sturation constants: H₁ $(p.H₁), H₂ $(p.H₂)%, H₃ $(p.H₃) g⋅m⁻², H₄ $(p.H₄) g⋅m⁻²",
    )
    println(
        io,
        "Intraspecific competition constant (g⋅m⁻²⋅day⁻¹) for bream: $(p.cb), pike: $(p.cp)",
    )
    println(
        io,
        "Pike food conversion efficiency to growth: $(p.ce), Mortalitiy rate (day⁻¹): $(p.mp)",
    )
    println(io, "Maximum vegetation coverage: $(p.K)%")
end


text:
## Wetted-area method
#======================================================#

# Raymer form factor for wing, pylon, nacelle, fuselage? (Eq. 12.30)
form_factor_wing(x_c, t_c, Λ_m, M) = (1 + 0.6t_c / x_c + 100t_c^4) * (1.34M^0.18 * cos(Λ_m)^0.28)

function form_factor(wing :: HalfWing, M)
    # (x/c)_max, (t/c)_max, (t/c)_max sweep angles
    xcs, tcs, sweeps = max_thickness_to_chord_ratio_sweeps(wing, 60)
    Kf  = form_factor_wing.(xcs, tcs, sweeps, M)
end

form_factor(wing :: Wing, M) = (form_factor(wing.left, M) + form_factor(wing.right, M)) / 2

# Schlichting averaged skin-friction coefficients
cf_lam(Re_c, k_lam = 1.) = 1.328 / √(Re_c * k_lam)
cf_turb(Re_c, M) = 0.455 / log10(Re_c)^2.58 / (1 + 0.144M^2)^0.65
cf_schlichting(Re, Re_xtr, M) = max(cf_lam(Re), cf_turb(Re, M) - (Re_xtr / 320 - 39) / Re)

function skin_friction_coefficients(mean_chords, x_tr, V, ρ, M, μ)
    # Skin-friction coefficients
    Re_c    = reynolds_number.(ρ, V, mean_chords, μ)
    Re_xtr  = reynolds_number.(ρ, V, mean_chords .* x_tr, μ)
    cfs     = cf_schlichting.(Re_c, Re_xtr, M)
end

function wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
    cfs = skin_friction_coefficients(mean_chords, x_tr, V, ρ, M, μ)

    # Profile drag
    Dp_by_q = sum(@. cfs * S_wets * K_fs)
end

function wetted_area_drag(wing :: HalfWing, x_tr, V, ρ, a_ref = 330., μ = 1.5e-5)
    # Chord processing
    mean_chords = (forward_sum ∘ chords)(wing) / 2

    # Wetted areas
    S_wets  = @. mean_chords * wing.spans / cos(wing.dihedrals)

    # Form factors
    M       = V / a_ref
    K_fs    = form_factor(wing, M)

    wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
end

function wetted_area_drag(wing :: WingMesh, x_tr, V, ρ, a_ref = 330., μ = 1.5e-5)
    # Chord processing
    mean_chords = (forward_sum ∘ chords)(wing.surf) / 2

    # Wetted areas
    surf_pans = camber_panels(wing)
    S_wets    = sum(panel_area, surf_pans, dims = 1)

    # Form factors
    M       = V / a_ref
    K_fs    = form_factor(wing.surf, M)

    wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
end

# Sato's local-friction and local-dissipation based on power balance method from Mark Drela, Flight Vehicle Aerodynamics, eq. 4.115.
function local_dissipation_drag(wing :: Wing, wetted_areas, ρ_es, u_es, x_tr, V, ρ, M, μ)
    # Chord processing
    mean_chords = (forward_sum ∘ chords)(wing) / 2

    # Compute weighted wetted areas based on inviscid edge velocity distribution.
    weighted_S_wets = sum(@. ρ_es * norm(u_es)^3 * wetted_areas; dims = 1) ./ (ρ * V^3)

    wetted_area_drag(mean_chords, weighted_S_wets, 1., x_tr, V, ρ, M, μ)
end

function local_dissipation_drag(wing :: WingMesh, ρ_es, u_es, x_tr, V, ρ, M, μ)
    # Chord processing
    mean_chords = (forward_sum ∘ chords)(wing.surf) / 2

    # Compute weighted wetted areas based on inviscid edge velocity distribution.
    S_wets = panel_area.(camber_panels(wing))
    weighted_S_wets = sum(@. ρ_es * u_es^3 * S_wets; dims = 1) ./ (ρ * V^3)

    wetted_area_drag(mean_chords, weighted_S_wets, 1., x_tr, V, ρ, M, μ)
end

abstract type AbstractProfileDrag end

struct FormFactor <: AbstractProfileDrag end
struct Dissipation <: AbstractProfileDrag end

profile_drag_coefficient(wing :: HalfWing, x_tr, V, rho_ref, a_ref, area_ref, μ) = wetted_area_drag(wing, x_tr, V, rho_ref, a_ref, μ) / area_ref
profile_drag_coefficient(wing :: WingMesh, x_tr, V, rho_ref, a_ref, area_ref, μ) = wetted_area_drag(wing, x_tr, V, rho_ref, a_ref, μ) / area_ref
profile_drag_coefficient(wing :: Wing, x_tr, V, rho_ref, a_ref, area_ref, μ) = profile_drag_coefficient(wing.left, x_tr, V, rho_ref, a_ref, area_ref, μ) + profile_drag_coefficient(wing.right, x_tr, V, rho_ref, a_ref, area_ref, μ)
profile_drag_coefficient(wing :: AbstractWing, x_tr, refs :: References) = profile_drag_coefficient(wing, x_tr, refs.speed, refs.density, refs.sound_speed, refs.area, refs.viscosity)

function wave_drag(M, Λ, t_by_c, Cl, κ_A)
    M_drag_divergence = κ_A / cos(Λ) - t_by_c / cos(Λ)^2 - Cl / (10 * cos(Λ)^3) # Drag divergence Mach number
    M_crit   = M_drag_divergnece - (0.1 / 80)^(1/3) # Critical Mach number
    CD_wave = ifelse(M > M_crit, 20 * (Mach - M_crit)^4, 0.)
end


text:
radau_stage_to_order(n::Int64) = 2 * n - 1
radau_rule_to_stage(n::Int64) = 2 * n - 1
radau_rule_to_order(n::Int64) = radau_rule_to_stage(radau_stage_to_order(n))

# get_radau_table_path() = joinpath(first(splitdir(Base.pathof(.Radau))), "table")

function get_radau_table_path()

    path_to_this_file = first(splitdir(Base.functionloc(load_radau_table_from_file)[1]))

    # println(@__MODULE__)
    # the_place = Base.functionloc(load_radau_table_from_file)[1]

    # println("the_place: ", the_place)
    joinpath(path_to_this_file, "table")
end

function load_radau_table_from_file(n_rule::Int64)
    radau_table_path = get_radau_table_path()
    radau_table_path_stage = joinpath(radau_table_path, string(n_rule) * "_rule")

    A      = readdlm(joinpath(radau_table_path_stage, "A.txt"), '\t', Float64, '\n')
    A⁻¹    = readdlm(joinpath(radau_table_path_stage, "inv_A.txt"), '\t', Float64, '\n')
    b      = readdlm(joinpath(radau_table_path_stage, "b.txt"), '\t', Float64, '\n')
    c      = readdlm(joinpath(radau_table_path_stage, "c.txt"), '\t', Float64, '\n')
    λ      = readdlm(joinpath(radau_table_path_stage, "lambda.txt"), '\t', Complex{Float64}, '\n')
    T      = readdlm(joinpath(radau_table_path_stage, "T.txt"), '\t', Complex{Float64}, '\n')
    T⁻¹    = readdlm(joinpath(radau_table_path_stage, "inv_T.txt"), '\t', Complex{Float64}, '\n')
    bi     = readdlm(joinpath(radau_table_path_stage, "bi.txt"), '\t', Float64, '\n')
    b_hat  = readdlm(joinpath(radau_table_path_stage, "b_hat.txt"), '\t', Float64, '\n')

    b = vec(b)
    c = Tuple(vec(c))
    λ = vec(λ)
    b_hat = vec(b_hat)

    b̂ = b_hat[2:end]
    b̂_0 = b_hat[1]

    return A, b, c, λ, T, T⁻¹, bi, b̂, b̂_0
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation_prod_vars::Gar1962.ProdVars,
    formulation_ramping::ArrCon2000.Ramping,
    formulation_status_vars::Gar1962.StatusVars,
)::Nothing
    # TODO: Move upper case constants to model[:instance]
    RESERVES_WHEN_START_UP = true
    RESERVES_WHEN_RAMP_UP = true
    RESERVES_WHEN_RAMP_DOWN = true
    RESERVES_WHEN_SHUT_DOWN = true
    gn = g.name
    RU = g.ramp_up_limit
    RD = g.ramp_down_limit
    SU = g.startup_limit
    SD = g.shutdown_limit
    reserve = model[:reserve]
    eq_ramp_down = _init(model, :eq_ramp_down)
    eq_ramp_up = _init(model, :eq_ramp_up)
    is_initially_on = (g.initial_status > 0)

    # Gar1962.ProdVars
    prod_above = model[:prod_above]

    # Gar1962.StatusVars
    is_on = model[:is_on]
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]

    for t in 1:model[:instance].time
        # Ramp up limit
        if t == 1
            if is_initially_on
                # min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant
                eq_ramp_up[gn, t] = @constraint(
                    model,
                    g.min_power[t] +
                    prod_above[gn, t] +
                    (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <=
                    g.initial_power + RU
                )
            end
        else
            max_prod_this_period =
                g.min_power[t] * is_on[gn, t] +
                prod_above[gn, t] +
                (
                    RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
                    reserve[gn, t] : 0.0
                )
            min_prod_last_period =
                g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]

            # Equation (24) in Kneuven et al. (2020)
            eq_ramp_up[gn, t] = @constraint(
                model,
                max_prod_this_period - min_prod_last_period <=
                RU * is_on[gn, t-1] + SU * switch_on[gn, t]
            )
        end

        # Ramp down limit
        if t == 1
            if is_initially_on
                # TODO If RD < SD, or more specifically if
                #      min_power + RD < initial_power < SD
                #      then the generator should be able to shut down at time t = 1,
                #      but the constraint below will force the unit to produce power
                eq_ramp_down[gn, t] = @constraint(
                    model,
                    g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD
                )
            end
        else
            max_prod_last_period =
                g.min_power[t-1] * is_on[gn, t-1] +
                prod_above[gn, t-1] +
                (
                    RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
                    reserve[gn, t-1] : 0.0
                )
            min_prod_this_period =
                g.min_power[t] * is_on[gn, t] + prod_above[gn, t]

            # Equation (25) in Kneuven et al. (2020)
            eq_ramp_down[gn, t] = @constraint(
                model,
                max_prod_last_period - min_prod_this_period <=
                RD * is_on[gn, t] + SD * switch_off[gn, t]
            )
        end
    end
end


text:
function MeshGrid(x,y)

    X = repeat(x', size(y,1), 1)
    Y = repeat(y, 1, size(x,1))
    return X, Y
    end #End Function


text:
"""
This file simulates the set of firms for step 3 of AFT2017.
"""

## Define module and things to be exported
module DGPsetup

export simulatefirms, vandercorput

## Load packages
using LinearAlgebra, Random, Distributions, Statistics, DataFrames
using CSV, XLSX, StatsBase

## Set seeds
Random.seed!(6)

## Define function for van der Corput sequence (fixed costs are drawn from here)
function vandercorput(n)
    # generate sequence of binary numbers from 1 to n
    d2 = size(digits(n; base = 2),1)
    b = 0.0 * ones(n,d2)
    for i in 1:n
        b[i,:] = digits(i; base = 2, pad = d2)'
    end
    # compute weights
    l = size(b,2)
    w = ones(l,1)
    for i in 1:l
        w[i,1] = i-l-1
    end
    w = reverse(vec(w))
    # get result
    x = b * (2 .^w)
    return x
end

## Define funcion to simulate firms
# ! Remember that bounds_intervals and length_intervals are 13x1 and 12x1 in
# the original codes

function simulatefirms(N; S_fixed = 18000)
    # in original paper S_fixed = 18000 and then S = 2160000; but I do not have so
    # much computational power, so I will reduce the amount of simulated firms

    bounds_intervals = [0;.2;.35;.45;.57;.7;.8;.9;.95;.98;.99;.999;1]
    # oversample more productive firms
    length_intervals = bounds_intervals[2:end] - bounds_intervals[1:end-1]
    num_draws_per_stratum = 10

    # Draw productivity shocks
    prod_draw_uniform = 1.0 * ones(num_draws_per_stratum*size(length_intervals,1),1)
    weights_prod = 1.0 * ones(num_draws_per_stratum*size(length_intervals,1),1)

    for k in 1:size(length_intervals,1)
        lb = (k-1)*num_draws_per_stratum + 1
        ub = k*num_draws_per_stratum
        prod_draw_uniform[lb:ub] = bounds_intervals[k,1] .+ rand(num_draws_per_stratum,1).*length_intervals[k,1]
        weights_prod[lb:ub] = (length_intervals[k,1] ./ num_draws_per_stratum) .* ones(num_draws_per_stratum,1)
    end
    prod_draw_uniform = vec(prod_draw_uniform)
    weights_prod = vec(weights_prod)

    # Draw fixed costs from van der Corput sequence
    corput_seq = vandercorput(S_fixed)
    fc_shock_randn = 0.0 * ones(S_fixed, N)
    for i in 1:N
        fc_shock_randn[:,i] = quantile.(Normal(0,1), shuffle(corput_seq))
        # quantile(d,X) is deprecated; use quantile.(d,X) instead
    end
    fc_shock_randn = repeat(fc_shock_randn, num_draws_per_stratum*size(length_intervals,1),1)

    # Adjust productivity draw and get number of simulated firms
    prod_draw_uniform = vec(kron(prod_draw_uniform, ones(S_fixed,1)))
    weights_prod = vec(kron(weights_prod, ones(S_fixed,1))) ./ S_fixed
    S = size(weights_prod,1)

    ## Add checks in case Jia's algorithm will not work
    num_rand_checks = 100
    rand_check_matrix = rand(num_rand_checks, N)

    return S, prod_draw_uniform, weights_prod, fc_shock_randn, num_rand_checks, rand_check_matrix
end

end     # end module