text:
const RIGHT_SCALING_COEFFICIENTS = Dict{Int, Vector{Vector{Float64}}}(
2 => Any[
[ 0.8705087534E+00 ; 0.4348969980E+00 ; 0.2303890438E+00 ]
,
[ -0.1942334074E+00 ; 0.1901514184E+00 ; 0.3749553316E+00 ; 0.7675566693E+00 ; 0.4431490496E+00 ]
]
,
3 => Any[
[ 0.9096849943E+00 ; 0.3823606559E+00 ; 0.1509872153E+00 ; 0.5896101069E-01 ]
,
[ -0.2904078511E+00 ; 0.4189992290E+00 ; 0.4969643721E+00 ; 0.4907578307E+00 ; 0.4643627674E+00 ; 0.1914505442E+00 ]
,
[ 0.8183541840E-01 ; -0.1587582156E+00 ; -0.9124735623E-01 ; 0.6042558204E-03 ; 0.7702933610E-01 ; 0.5200601778E+00 ; 0.7642591993E+00 ; 0.3150938230E+00 ]
]
,
4 => Any[
[ 0.9154705188E+00 ; 0.3919142810E+00 ; 0.5947771124E-01 ; -0.2519180851E-01 ; 0.6437934569E-01 ]
,
[ -0.2191626469E+00 ; 0.4488001781E+00 ; 0.7540005084E+00 ; 0.3937758157E+00 ; -0.1581338944E+00 ; -0.1614201190E-01 ; 0.4126840881E-01 ]
,
[ 0.1290078289E-01 ; -0.1390716006E+00 ; 0.2921367950E-01 ; 0.4606168537E+00 ; 0.8164119742E+00 ; 0.2986473346E+00 ; -0.1027663536E+00 ; -0.1257488211E-01 ; 0.3214874197E-01 ]
,
[ -0.6775603652E-02 ; 0.1913244129E-01 ; -0.1770918425E-01 ; -0.6765916174E-01 ; -0.3023588481E-01 ; 0.4977920821E+00 ; 0.8039495996E+00 ; 0.2977111011E+00 ; -0.9910804055E-01 ; -0.1259895190E-01 ; 0.3221027840E-01 ]
]
,
5 => Any[
[ 0.6629994791E+00 ; 0.5590853114E+00 ; -0.3638895712E+00 ; -0.3190292420E+00 ; -0.8577516141E-01 ; 0.7938922949E-01 ]
,
[ -0.1615976602E+00 ; 0.7429922432E+00 ; 0.5771521705E+00 ; 0.2773360273E+00 ; 0.1530007177E-01 ; -0.1063942878E+00 ; -0.1216769629E-01 ; 0.1126647052E-01 ]
,
[ 0.2576681838E+00 ; -0.3995181391E-01 ; 0.1290658921E-01 ; 0.4163708782E+00 ; 0.6469374504E+00 ; 0.5608421202E+00 ; 0.4995700940E-02 ; -0.1569482750E+00 ; -0.2009485730E-01 ; 0.1860648984E-01 ]
,
[ -0.8955158953E-01 ; 0.3895616532E-01 ; 0.4751157158E-02 ; -0.6596188737E-01 ; 0.1520881337E-01 ; 0.2466644601E+00 ; 0.7199998123E+00 ; 0.6131240289E+00 ; 0.1537319456E-01 ; -0.1721947543E+00 ; -0.2083858871E-01 ; 0.1929513523E-01 ]
,
[ 0.1546583496E-01 ; -0.6559247482E-02 ; 0.2728093541E-02 ; 0.1830645747E-01 ; 0.1708296448E-01 ; 0.2156026058E-01 ; -0.3692524465E-01 ; 0.2045334663E+00 ; 0.7233635779E+00 ; 0.6327590443E+00 ; 0.1649999072E-01 ; -0.1751674471E+00 ; -0.2109311508E-01 ; 0.1953080958E-01 ]
]
,
6 => Any[
[ 0.9236675275E+00 ; 0.3793343885E+00 ; 0.4986459778E-01 ; -0.1688947163E-01 ; -0.3311806998E-02 ; 0.2804629640E-02 ; -0.1237553624E-01 ]
,
[ -0.2831955072E+00 ; 0.6209232596E+00 ; 0.6494435367E+00 ; 0.3279016732E+00 ; 0.4613465423E-01  ; -0.2201105042E-01 ; 0.4772996268E-01 ; 0.1811083865E-02 ; -0.7992103956E-02 ]
,
[ 0.1415889950E+00 ; -0.3264117354E+00 ; 0.5898718517E-01 ; 0.5662575494E+00 ; 0.6511436245E+00 ; 0.3434709841E+00 ; -0.6897565912E-01 ; -0.2139852264E-01 ; 0.4514540543E-01 ; 0.1807019829E-02 ; -0.7974169836E-02 ]
,
[ -0.4080929146E-01 ; 0.1251977353E+00 ; -0.9208945297E-01 ; -0.1815441604E+00 ; 0.7492893632E-01 ; 0.4921890166E+00 ; 0.7598885402E+00 ; 0.3298535864E+00 ; -0.6885203537E-01 ; -0.2056224689E-01 ; 0.4357038095E-01 ; 0.1729458762E-02 ; -0.7631901801E-02 ]
,
[ 0.1473338385E-01 ; -0.3879928306E-01 ; 0.3222685441E-01 ; 0.4181351746E-01 ; -0.4470160402E-01 ; -0.1143592801E+00 ; -0.3079264597E-01 ; 0.4953303341E+00 ; 0.7831063037E+00 ; 0.3365261769E+00 ; -0.7161554542E-01 ; -0.2096645980E-01 ; 0.4445995152E-01 ; 0.1762244193E-02 ; -0.7776580122E-02 ]
,
[ -0.2279778462E-02 ; 0.5934267127E-02 ; -0.4653920620E-02 ; -0.4697326198E-02 ; 0.7905498252E-02 ; 0.1478426761E-01 ; 0.6356609850E-03 ; -0.1184137622E+00 ; -0.4719524338E-01 ; 0.4914056494E+00 ; 0.7873827096E+00 ; 0.3378777169E+00 ; -0.7256793624E-01 ; -0.2105690525E-01 ; 0.4471373463E-01 ; 0.1767573251E-02 ; -0.7800096641E-02 ]
]
,
7 => Any[
[ 0.9091449027E+00 ; 0.4122956702E+00 ; 0.8463996964E-02 ; -0.4985469571E-01 ; -0.2885243693E-01 ; -0.8468475604E-02 ; 0.9290995299E-03 ; 0.2377975298E-02 ]
,
[ -0.2699308042E+00 ; 0.6287804833E+00 ; 0.6425124921E+00 ; 0.3295936378E+00 ; 0.9680040842E-01 ; -0.1505090528E-02 ; -0.2576996713E-01 ; -0.1663793375E-01 ; 0.7886583178E-03 ; 0.2019348750E-02 ]
,
[ 0.2157402688E+00 ; -0.3660961628E+00 ; 0.9979560410E-01 ; 0.5343306371E+00 ; 0.5651956498E+00 ; 0.3984651316E+00 ; 0.2006608171E+00 ; 0.3790996947E-01 ; -0.5279896883E-01 ; -0.3509616877E-01 ; 0.1600620645E-02 ; 0.4098367097E-02 ]
,
[ -0.1219999250E+00 ; 0.2497494925E+00 ; -0.1489636838E+00 ; -0.2314941689E+00 ; 0.4306729538E-01 ; 0.3349988297E+00 ; 0.5170449138E+00 ; 0.5314908191E+00 ; 0.3984867666E+00 ; 0.1134176909E+00 ; -0.8664578083E-01 ; -0.6206801690E-01 ; 0.2555164772E-02 ; 0.6542464174E-02 ]
,
[ 0.6643773700E-01 ; -0.1425076588E+00 ; 0.1150002072E+00 ; 0.1110261320E+00 ; -0.5456405522E-01 ; -0.1403852745E+00 ; -0.9002913124E-01 ; 0.9927107272E-01 ; 0.3252695570E+00 ; 0.5955739150E+00 ; 0.6177120462E+00 ; 0.2102093627E+00 ; -0.1207077310E+00 ; -0.9029360085E-01 ; 0.3498469159E-02 ; 0.8957782052E-02 ]
,
[ -0.2472698999E-01 ; 0.5451727819E-01 ; -0.4913268141E-01 ; -0.3543504434E-01 ; 0.2981657384E-01 ; 0.5443081267E-01 ; 0.3080352023E-01 ; -0.2311127518E-01 ; -0.2663146891E-01 ; -0.3986350788E-01 ; 0.1012109980E+00 ; 0.5654368499E+00 ; 0.7419995864E+00 ; 0.2730680885E+00 ; -0.1373359727E+00 ; -0.1050857718E+00 ; 0.3942767529E-02 ; 0.1009540190E-01 ]
,
[ 0.4060965940E-02 ; -0.9082349871E-02 ; 0.8675066430E-02 ; 0.5263347825E-02 ; -0.5426164014E-02 ; -0.8281419894E-02 ; -0.2928559727E-02 ; 0.7407888982E-02 ; 0.1392100926E-02 ; 0.2865053765E-01 ; 0.5527738563E-01 ; -0.5250033174E-01 ; 0.2337385282E-01 ; 0.5390029057E+00 ; 0.7666684946E+00 ; 0.2878132266E+00 ; -0.1399692009E+00 ; -0.1077108034E+00 ; 0.4008607630E-02 ; 0.1026398458E-01 ]
]
,
8 => Any[
[ 0.9280961080E+00 ; 0.3685482772E+00 ; 0.5169019500E-01 ; -0.1072856966E-01 ; -0.3547264354E-02 ; 0.1881724221E-02 ; 0.2253426546E-03 ; -0.4074549423E-03 ; 0.2542121693E-02 ]
,
[ -0.2880715337E+00 ; 0.6422316271E+00 ; 0.6540393284E+00 ; 0.2730170577E+00 ; 0.4300012686E-01 ; -0.1532507378E-01 ; -0.3095082211E-02 ; 0.3081063250E-02 ; -0.1192117285E-01 ; -0.2510321213E-03 ; 0.1566213637E-02 ]
,
[ 0.1554699662E+00 ; -0.3928508188E+00 ; 0.1655006354E+00 ; 0.6186900270E+00 ; 0.5757733437E+00 ; 0.2765573214E+00 ; 0.4215137090E-01 ; -0.2022141803E-01 ; 0.3170527734E-01 ; 0.2969565406E-02 ; -0.1133537854E-01 ; -0.2472629557E-03 ; 0.1542697449E-02 ]
,
[ -0.9273812077E-01 ; 0.2525480154E+00 ; -0.1903081753E+00 ; -0.2325263183E+00 ; 0.1951123754E+00 ; 0.5747752998E+00 ; 0.6031889613E+00 ; 0.3197987067E+00 ; -0.2318110629E-01 ; -0.2335276594E-01 ; 0.3904642048E-01 ; 0.3358324305E-02 ; -0.1296942483E-01 ; -0.2744734004E-03 ; 0.1712466040E-02 ]
,
[ 0.4096472123E-01 ; -0.1240767851E+00 ; 0.1278114403E+00 ; 0.7138489316E-01 ; -0.1584871307E+00 ; -0.1727866355E+00 ; 0.1143640280E+00 ; 0.5129116028E+00 ; 0.7137757767E+00 ; 0.3414317954E+00 ; -0.3881329228E-01 ; -0.2526829049E-01 ; 0.4424255303E-01 ; 0.3575822034E-02 ; -0.1394488775E-01 ; -0.2875903358E-03 ; 0.1794303867E-02 ]
,
[ -0.1790191282E-01 ; 0.5295161530E-01 ; -0.5704547538E-01 ; -0.1894676577E-01 ; 0.6968256270E-01 ; 0.4547075753E-01 ; -0.7337937479E-01 ; -0.1444748191E+00 ; -0.1409191494E-01 ; 0.4948704363E+00 ; 0.7618895951E+00 ; 0.3591100053E+00 ; -0.4810103061E-01 ; -0.2674149513E-01 ; 0.4783092078E-01 ; 0.3754819526E-02 ; -0.1470868404E-01 ; -0.2997260590E-03 ; 0.1870019815E-02 ]
,
[ 0.4842573875E-02 ; -0.1435818114E-01 ; 0.1564773617E-01 ; 0.3635704163E-02 ; -0.1807616145E-01 ; -0.7951972463E-02 ; 0.2177252927E-01 ; 0.3159516732E-01 ; -0.5393073080E-02 ; -0.1454356948E+00 ; -0.5484912621E-01 ; 0.4834662634E+00 ; 0.7753953496E+00 ; 0.3639514649E+00 ; -0.5141318290E-01 ; -0.2717238099E-01 ; 0.4899043649E-01 ; 0.3804055936E-02 ; -0.1492982098E-01 ; -0.3026846009E-03 ; 0.1888478444E-02 ]
,
[ -0.5886912723E-03 ; 0.1750384584E-02 ; -0.1931250902E-02 ; -0.3621632497E-03 ; 0.2085884179E-02 ; 0.6409242844E-03 ; -0.2736289413E-02 ; -0.3368917714E-02 ; 0.1038171472E-02 ; 0.3193936333E-01 ; 0.6823471193E-02 ; -0.1435014917E+00 ; -0.6098305316E-01 ; 0.4814477634E+00 ; 0.7771077176E+00 ; 0.3644295324E+00 ; -0.5192289139E-01 ; -0.2721779258E-01 ; 0.4913249792E-01 ; 0.3808651449E-02 ; -0.1495173150E-01 ; -0.3029170562E-03 ; 0.1889928755E-02 ]
]
)


text:
import Random: rand, seed!
import RandomNumbers: AbstractRNG

# Well1024a pseudo random number generator from "Improved Long-Period Generators Based on Linear Recurrences Modulo 2" f
# by Panneton and L'Ecuyer (2008)

const M1 = 3
const M2 = 24
const M3 = 10
const R = 32
const MASK_MODULO_32 = Int(0x0000001f)
mutable struct Well1024a <: AbstractRNG{UInt32}
    state::Vector{UInt32}
    stateIndex::Int
end

function Well1024a(seed::UInt32 = UInt32(20132109))
    r = Well1024a(Vector{UInt32}(undef, R), 0)
    seed!(r, seed)
end

function seed!(r::Well1024a, seed::UInt32)
    r.stateIndex = 0
    r.state[1] = seed
    @inbounds for j = 2:R
        r.state[j] = 0x6c078965 * (r.state[j-1] ⊻ (r.state[j-1] >> 30)) + (j - 1) % UInt32
    end
    r
end

@inline function rand(r::Well1024a, ::Type{UInt32})
    @inbounds begin
        sM1 = r.state[((r.stateIndex+M1)&MASK_MODULO_32)+1]
        sM2 = r.state[((r.stateIndex+M2)&MASK_MODULO_32)+1]
        sM3 = r.state[((r.stateIndex+M3)&MASK_MODULO_32)+1]
        z2 = sM2 ⊻ sM3 ⊻ (sM2 << 19) ⊻ (sM3 << 14)
        z1 = r.state[r.stateIndex+1] ⊻ sM1 ⊻ (sM1 >> 8)
        z0 = r.state[((r.stateIndex+31)&MASK_MODULO_32)+1]
        z12 = z1 ⊻ z2
        z3 = z12 ⊻ z0 ⊻ (z2 << 13) ⊻ (z0 << 11) ⊻ (z1 << 7)
        r.state[r.stateIndex+1] = z12
        r.stateIndex = (r.stateIndex + 31) & MASK_MODULO_32
        r.state[r.stateIndex+1] = z3
        z3
    end
end

ext: jl


text:
# A (sloppy) example of some actual analysis done using PressureDrop.jl to generate normalized pressure plots for several wells through time.

wells = ["<wellname>",  "<wellname>"]
baseline_pressure = Dict("<wellname>" => 3125, "<wellname>" => 3191)
patternplot = true
max_depth = 8500
start_date = "21-Jul-2019"

import Printf
using DataFrames
using RCall

function initialize_R()

    R"""
    library(DBI)
    library(dplyr)
    library(tidyr)

    print(paste0('Connecting to ', $db))

    connection <- DBI::dbConnect(odbc::odbc(), driver = 'SQL Server',
                                server = '<Aries Server>', database = '<Aries DB>')
    """
end

# note that variable interpolation in @R_str is actual variable interpolation, not string interpolation
function get_production(well)

    R"""

    query <- paste0("
    select d.D_DATE as Date, p.PropNum,
    isnull(case when d.Oil < 0 then 0 else d.Oil end,0) as Oil,
    isnull(case when d.Gas < 0 then 0 else d.Gas end,0) as Gas,
    isnull(case when d.Water < 0 then 0 else d.Water end,0) as Water,
    isnull(d.Press_FTP,0) as TP,
    isnull(d.Press_Csg,0) as CP,
    isnull(d.Press_Inj,0) as GasInj,
    d.REMARKS as Comments

    from ARIES.AC_DAILY as d
    inner join ARIES.AC_PROPERTY as p
    on d.PROPNUM = p.PROPNUM

    WHERE p.AREA = 'oklahoma'
    and p.LEASE like '%",$well,"%'
    and d.D_DATE >= '",$start_date,"'
    ORDER BY D_DATE asc")

    prod <- dbGetQuery(connection, query) %>%
                tidyr::fill(Water, TP, CP, GasInj, .direction = 'down') %>%
                replace_na(list(Oil = 0, Gas = 0))

    """

    @rget prod
    @assert length(unique(prod.PropNum)) == 1 "Query for $well accidentally pulled in $(length(unique(prod.PropNum))) propnums"

    return prod
end

function finish_R()

    R"""
    dbDisconnect(connection)
    """
end

using PressureDrop
using Gadfly

function calculate_pressures(well, prod, default_GOR = 500)

    well_path = read_survey(path = "Surveys/$well.csv", skiplines = 4, maxdepth = max_depth)
    valves = read_valves(path = "GL designs/$well.csv", skiplines = 1)

    model = WellModel(wellbore = well_path, roughness = 0.001,
                    valves = valves, pressurecorrelation = HagedornAndBrown,
                    WHP = 0, CHP = 0, dp_est = 25, temperature_method = "Shiu",
                    BHT = 160, geothermal_gradient = 0.8,
                    q_o = 0, q_w = 500,
                    GLR = 0, naturalGLR = 0,
                    APIoil = 38.2, sg_water = 1.05, sg_gas = 0.65)

    FBHP = Array{Float64, 1}(undef, nrow(prod))

    for day in 1:nrow(prod)
        if prod.Oil[day] + prod.Water[day] == 0
            FBHP[day] = NaN
        else
            model.WHP = prod.TP[day]
            model.CHP = prod.CP[day]
            model.q_o = prod.Oil[day]
            model.q_w = prod.Water[day]
            model.naturalGLR = model.q_o + model.q_w <= 0 ? 0 : max( prod.Gas[day] * 1000 / (model.q_o + model.q_w) , default_GOR * model.q_o / (model.q_o + model.q_w) ) #force minimum GLR
            model.GLR = model.q_o + model.q_w <= 0 ? 0 : max( (prod.Gas[day] + prod.GasInj[day]) * 1000 / (model.q_o + model.q_w) , model.naturalGLR)

            FBHP[day] = gaslift_model!(model, find_injectionpoint = true, dp_min = 100) |> x -> x[1][end]
        end
    end

    return model, FBHP, gaslift_model!(model, find_injectionpoint = true, dp_min = 100)
end

ticks = log10.([0.1:0.1:0.9;1.0:1:10] |> collect)
function plot_normP_data(production, FBHP, well)

    production = production[.!isnan.(FBHP), :]
    FBHP = FBHP[.!isnan.(FBHP)]

    initial_pressure = haskey(baseline_pressure, well) ? baseline_pressure[well] : maximum(FBHP)
    println("Baseline pressure: $initial_pressure")

    norm_rate = (production.Oil .+ production.Water) ./ (initial_pressure .- FBHP)
    println("Initial PI: $(norm_rate[1])")
    println("Max PI: $(maximum(norm_rate))")
    println("Final PI: $(norm_rate[end])")

    return plot(layer(x = 1:length(norm_rate), y = norm_rate, Geom.path),
                Scale.y_log10(minvalue = 10^-1., maxvalue = 10^1, labels = y-> Printf.@sprintf("%0.1f", 10^y)), Scale.x_log10,
                Guide.ylabel("Normalized Total Fluid (bpd / ΔP)"),
                Guide.xlabel("Normalized Time (days)"),
                Guide.yticks(ticks = ticks),
                Guide.title(well))
end

Base.CoreLogging.disable_logging(Base.CoreLogging.Warn) #prevent dumping all of the @infos from normal PressureDrop functions

initialize_R()
GL_plots = Dict()
pressure_plots = Dict()
for well in wells
    println("Evaluating $well...")

    production = get_production(well)[1:end-1, :] #last day is usually 0 prod
    model, FBHPs, GL_data = calculate_pressures(well, production); #GL data is tuple of TP, CP, valve data

    #generate plots 1 by 1
    GL_plots[well] = plot_gaslift(model.wellbore, GL_data[1], GL_data[2], model.temperatureprofile, GL_data[3], well)
    pressure_plots[well] = plot_normP_data(production, FBHPs, well)

    println("$well complete")
end
finish_R()

#set_default_plot_size(6inch,8inch)

using Compose

# pressure plots
if patternplot
    # lay out the wells in question according to gunbarrel view
    set_default_plot_size(12inch,5inch)
    gridstack(Union{Plot,Compose.Context}[pressure_plots["Mayes 1706 2-16MH"]  pressure_plots["Mayes 1706 4-16MH"] ]) |>
                                          SVG("pattern.svg")

else
    set_default_plot_size(6inch,8inch)

    for well in wells
        GL_plots[well] |> SVG("GL plots/$well GL.svg")
    end
end

max_stars_repo_path: examples/actual multi-well example.jl

max_issues_repo_path: examples/actual multi-well example.jl

max_forks_repo_path: examples/actual multi-well example.jl


text:
using StructIO
using COFF
using ObjFileBase
@struct immutable RUNTIME_FUNCTION
    startoff::UInt32
    endoff::UInt32
    unwindoff::UInt32
end

immutable PData <: AbstractArray{RUNTIME_FUNCTION, 1}
    sec::SectionRef
end
import Base: getindex, length
PData(h::COFF.COFFHandle) = PData(
    collect(filter(x->sectionname(x)==ObjFileBase.mangle_sname(h,"pdata"),Sections(h)))[])
getindex(it::PData, n) = (1 <= n <= length(it) || throw(BoundsError());
    seek(it.sec, (n-1)*sizeof(RUNTIME_FUNCTION)); unpack(handle(it.sec), RUNTIME_FUNCTION))
length(it::PData) = Int(div(sectionsize(it.sec), sizeof(RUNTIME_FUNCTION)))

@struct immutable UNWIND_INFO
    # Low 3 bits are version
    # High 5 bits are flags
    VersionFlags::UInt8
    PrologueSize::UInt8
    NumUnwindCodes::UInt8
    FrameRegOff::UInt8
    # Opcodes
    # Exception Handler or chained unwind info
end

immutable UNWIND_CODE
    PrologueOffset::UInt8
    Operation::UInt8
end

@COFF.constants WIN64_UWOP "" begin
    const UWOP_PUSH_NONVOL      =  0
    const UWOP_ALLOC_LARGE      =  1
    const UWOP_ALLOC_SMALL      =  2
    const UWOP_SET_FPREG        =  3
    const UWOP_SAVE_NONVOL      =  4
    const UWOP_SAVE_NONVOL_FAR  =  5
    const UWOP_SAVE_XMM128      =  8
    const UWOP_SAVE_XMM128_FAR  =  9
    const UWOP_PUSH_MACHFRAME   =  0xa
end

_opinfo(Operation) = (Operation & 0xf, (Operation & 0xf0) >> 4)
opinfo(op) = _opinfo(op.Operation)
opinfo(op::UInt16) = _opinfo((op&0xff00)>>8)

function compute_alloc_op_size(ops, i)
    opcode, info = opinfo(ops[i])
    @assert opcode == UWOP_ALLOC_LARGE || opcode == UWOP_ALLOC_SMALL
    (opcode == UWOP_ALLOC_SMALL) && return (8 * (info + 1), 1)
    size = 0; nskip = 1
    if info == 0
        size = 8*ops[i+1]
        nskip = 2
    elseif info == 1
        size = UInt32(ops[i+1]) |
            (UInt32(ops[i+2]) << 8)
        nskip = 3
    else
        error("Invalid operation info")
    end
    (size, nskip)
end

function print_op(io, ops, i)
    opcode, info = opinfo(ops[i])
    print(io, WIN64_UWOP[opcode])
    if opcode == UWOP_PUSH_NONVOL
        println(io, ", register=", X86_64.seh_numbering[info])
        return 1
    elseif opcode == UWOP_ALLOC_LARGE || opcode == UWOP_ALLOC_SMALL
        size, nskip = compute_alloc_op_size(ops, i)
        println(io, ", size=", size)
        return nskip
    else
        println(io)
        return 1
    end
end

const UNW_FLAG_NHANDLER     =  0
const UNW_FLAG_EHANDLER     =  1
const UNW_FLAG_UHANDLER     =  2
const UNW_FLAG_CHAININFO    =  4

immutable UnwindEntry
    start::UInt64
    stop::UInt64
    info::UNWIND_INFO
    opcodes::Vector{UInt16}
end

immutable XPUnwindRef{SR<:ObjFileBase.SectionRef} <: Base.AbstractArray{UnwindEntry,1}
    xdata::SR
    pdata::SR
end
XPUnwindRef(xdata::SectionRef, pdata::SectionRef) = XPUnwindRef(xdata, pdata)
Base.length(xp::XPUnwindRef) = length(PData(xp.pdata))
Base.size(xp::XPUnwindRef) = (length(xp),)
function Base.getindex(xp::XPUnwindRef, idx)
    pentry = PData(xp.pdata)[idx]
    xdataoff = pentry.unwindoff - deref(xp.xdata).VirtualAddress
    seek(xp.xdata, xdataoff)
    info = unpack(handle(xp.xdata), UNWIND_INFO)
    UnwindEntry(pentry.startoff, pentry.endoff, info,
        read(handle(xp.xdata), UInt16, info.NumUnwindCodes))
end


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
np = pyimport("numpy")

"""
1  -  2  -  3  -  4

local node id: left to right
"""

ndofs = 2
nx = 4
nnodes, neles = (nx + 1), nx
ngp=2
T = 200
NT = 200
Δt = T/NT

function BoundaryCondition(tid::Int64)
# PARAMETER: distance between nodes 0-1 and 0-2
l = 1.0
nodes = zeros(Float64, nnodes, ndofs)
for i = 1:nx+1
    nodes[i, 1], nodes[i, 2] = (i - 1)*l, 0.0
end

EBC, g = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)
#fix node 1
EBC[1, :] .= -1

function ggt(t)
    return zeros(sum(EBC.==-2)), zeros(sum(EBC.==-2))
end
gt = ggt

#pull in the x direction
FBC, fext = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)

# todo PARAMETER
FORCE_TYPE = "non-constant"

if FORCE_TYPE == "constant"
    FBC[nx+1, 1]  = -1
    fext[nx+1, 1] = 1.0e3
else
    FBC[nx+1, 1] = -2
end

#force load function
function ft(t)
    return  sin(pi*t/T) * (0.4*tid + 1.6)
end

npoints = nnodes
node_to_point = collect(1:npoints)
return nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point
end


text:
function solve_001()
    sum = 0
    for k = 1:999
        if k % 3 == 0 || k % 5 == 0
            sum += k
        end
    end
    return sum
end


text:
module TexAssistant

    import ArgParse
    import TOML

    include("utils.jl")
    include("regexs.jl")
    include("split_comment.jl")
    include("expand_input.jl")
    include("extract_sections_to_files.jl")
    include("settings.jl")
    include("extractors.jl")
    include("tex_match.jl")

    function __init__()
        _fill_regexs!()
    end

end


text:
using Documenter, MendelPlots

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = Documenter.HTML(),
    sitename = "MendelPlots",
    modules = [MendelPlots]
)

deploydocs(
    repo   = "github.com/OpenMendel/MendelPlots.jl.git",
    target = "build"
)


text:
# here tests should be written for utility functions in src/Utils.jl
using MonteCarloX
using StatsBase

function test_histogram_set_get(;verbose = false)
    pass = true

    # bin uniform values into histogram (each bin has 10 elements)
    list_vals = [i for i = 1:100]
    hist = fit(Histogram, list_vals, 1:10:101)

    # reset the size of the histogram bins (accessing each bin via random
    # elements that are included in the bin)
    target = [i for i = 1:10]
    hist[1] = target[1]
    hist[12] = target[2]
    hist[23] = target[3]
    hist[34] = target[4]
    hist[45] = target[5]
    hist[56] = target[6]
    hist[67] = target[7]
    hist[78] = target[8]
    hist[89] = target[9]
    hist[100] = target[10]

    # check that the histgram entries correspond to the target values
    for i = 1:100
        if verbose
            println("... $(hist[i]) == $(target[1 + floor(Int, (i - 1) / 10)])")
        end
        pass &= hist[i] == target[1 + floor(Int, (i - 1) / 10)]
    end

    return pass
end


text:
###############################################################################
# Affine

affine(f::AffineTransform, μ::AbstractMeasure) = Affine(f, μ)

affine(nt::NamedTuple, μ::AbstractMeasure) = affine(AffineTransform(nt), μ)

affine(f) = μ -> affine(f, μ)

function affine(f::AffineTransform, parent::WeightedMeasure)
    WeightedMeasure(parent.logweight, affine(f, parent.base))
end

function affine(f::AffineTransform, parent::FactoredBase)
    constℓ = parent.constℓ
    varℓ = parent.varℓ
    # Avoid transforming `inbounds`, which is expensive
    base = affine(f, restrict(parent.inbounds, parent.base))
    FactoredBase(Returns(true), constℓ, varℓ, base)
end

###############################################################################
# Half

half(μ::AbstractMeasure) = Half(μ)

###############################################################################
# PointwiseProductMeasure

pointwiseproduct(μ::AbstractMeasure...) = PointwiseProductMeasure(μ)

function pointwiseproduct(
    μ::PointwiseProductMeasure{X},
    ν::PointwiseProductMeasure{Y}
) where {X,Y}
    data = (μ.data..., ν.data...)
    pointwiseproduct(data...)
end

function pointwiseproduct(μ::AbstractMeasure, ν::PointwiseProductMeasure)
    data = (μ, ν.data...)
    return pointwiseproduct(data...)
end

function pointwiseproduct(μ::PointwiseProductMeasure, ν::N) where {N<:AbstractMeasure}
    data = (μ.data..., ν)
    return pointwiseproduct(data...)
end

function pointwiseproduct(μ::AbstractMeasure, ℓ::Likelihood)
    data = (μ, ℓ)
    return PointwiseProductMeasure(data)
end

###############################################################################
# PowerMeaure

function powermeasure(μ::M, dims::NTuple{N,I}) where {M<:AbstractMeasure,N,I}
    productmeasure(Returns(μ), identity, CartesianIndices(dims))
end

function powermeasure(μ::M, dims::Tuple{I}) where {M<:AbstractMeasure,N,I}
    productmeasure(Returns(μ), identity, Base.OneTo(first(dims)))
end

function powermeasure(μ::WeightedMeasure, dims::NTuple{N,I}) where {N,I}
    k = prod(dims) * μ.logweight
    return weightedmeasure(k, μ.base^dims)
end

###############################################################################
# ProductMeasure

productmeasure(f, ops, pars) = ProductMeasure(kernel(f, ops), pars)

productmeasure(f, pars) = productmeasure(f, identity, pars)

productmeasure(μs::Tuple) = TupleProductMeasure(μs)

productmeasure(f::Returns, ops, pars) = ProductMeasure(kernel(f, identity), pars)

productmeasure(k::Kernel, pars) = productmeasure(k.f, k.ops, pars)

productmeasure(nt::NamedTuple) = productmeasure(identity, nt)

function productmeasure(f::Returns{FB}, ops, pars) where {FB<:FactoredBase}
    fb = f.value
    dims = size(pars)
    n = prod(dims)
    inbounds(x) = all(fb.inbounds, x)
    constℓ = n * fb.constℓ
    varℓ() = n * fb.varℓ()
    base = fb.base^dims
    FactoredBase(inbounds, constℓ, varℓ, base)
end

function productmeasure(f::Returns{W}, ::typeof(identity), pars) where {W<:WeightedMeasure}
    ℓ = f.value.logweight
    base = f.value.base
    newbase = productmeasure(Returns(base), identity, pars)
    weightedmeasure(length(pars) * ℓ, newbase)
end

###############################################################################
# RestrictedMeasure

restrict(f, b) = RestrictedMeasure(f, b)

###############################################################################
# SuperpositionMeasure

superpose(a::AbstractArray) = SuperpositionMeasure(a)

superpose(t::Tuple) = SuperpositionMeasure(t)
superpose(nt::NamedTuple) = SuperpositionMeasure(nt)

function superpose(μ::AbstractMeasure, ν::AbstractMeasure)
    components = (μ, ν)
    superpose(components)
end

###############################################################################
# WeightedMeasure

function weightedmeasure(ℓ::R, b::M) where {R,M}
    WeightedMeasure{R,M}(ℓ, b)
end

function weightedmeasure(ℓ, b::WeightedMeasure)
    weightedmeasure(ℓ + b.logweight, b.base)
end

###############################################################################
# Kernel

kernel(μ, ops...) = Kernel(μ, ops)
kernel(μ, op) = Kernel(μ, op)

# kernel(Normal(μ=2))
function kernel(μ::P) where {P<:AbstractMeasure}
    (f, ops) = kernelfactor(μ)
    kernel(f, ops)
end

# kernel(Normal{(:μ,), Tuple{Int64}})
function kernel(::Type{P}) where {P<:AbstractMeasure}
    (f, ops) = kernelfactor(P)
    kernel(f, ops)
end

# kernel(::Type{P}, op::O) where {O, N, P<:ParameterizedMeasure{N}} = kernel{constructorof(P),O}(op)

function kernel(::Type{M}; ops...) where {M}
    nt = NamedTuple(ops)
    kernel(M, nt)
end

kernel(f::Returns, op::typeof(identity)) = Kernel(f, op)
kernel(f::Returns, op) = kernel(f, identity)


text:
"""
$(TYPEDEF)

One sample superiority test for Mean

Constructors
------------
* `OneSampleMeanSuperior(mu1::Real, mu0::Real, delta::Real, stdunknown::Bool)`

Arguments
---------
* `mu1`: True mean

* `mu0`: Mean under null hypothesis

* `delta`: Superiority margin

* `stdunknown`: If the population standard deviation is known, default = false,
when the population parameter is known, the power / sample size are calculated based
on a standard normal distribution, otherwise a t distribution will be used

Fields
------
$(FIELDS)
"""
type OneSampleMeanSuperior <: TrialTest
    mu1::Real
    mu0::Real
    delta::Real
    stdunknown::Bool

    # Validator
    function OneSampleMeanSuperior(mu1, mu0, delta, stdunknown = false)

        if !((-Inf < mu1 < Inf) & (-Inf < mu0 < Inf))

            error("Mean value must be in (-Inf, Inf)")

        end # end if

        if !(0 <= delta < Inf)

            error("The superiority margin δ must be in [0, Inf)")

        end # end if
        new(mu1, mu0, delta, stdunknown)

    end # end function

end # end type

# One-sample superiority test for mean
function hypotheses{T <: OneSampleMeanSuperior}(test::T, n::Real, std::Real, alpha::Real, side::String)

    diff = test.mu1 - test.mu0 - test.delta
    se = sqrt(1 / n) * std
    ts = diff / se
    if test.stdunknown == false

        p = cdf(ZDIST, ts - quantile(ZDIST, 1 - alpha)) + cdf(ZDIST, -ts - quantile(ZDIST, 1 - alpha))

    else

        p = cdf(ZDIST, ts - quantile(TDist(n - 1), 1 - alpha)) + cdf(TDist(n - 1), -ts - quantile(TDist(n - 1), 1 - alpha))

    end # end if
    # power_result(test, p)
    return p

end # end function


text:
using Test
using NamedIndexing

@testset "NamedIndexing" begin

@testset "NamedAxisArray" begin
    A = NamedAxisArray(rand(3, 2), (:foo, :bar))
    @test axisnames(A) == (:foo, :bar)

    Acolons = A[foo=:, bar=:]
    @test axisnames(Acolons) == axisnames(A)
    @test Acolons.data == A.data
end

end


text:
using MAT, JuMP, Ipopt
include("../benchmark.jl")

function read_lp2(dir, name)
    lp_data = matopen("$(dir)/$(name).mat")

    problem = read(lp_data,"Problem")
    A = problem["A"]
    b = problem["b"][:]
    c = problem["aux"]["c"][:]
    lb = problem["aux"]["lo"][:]
    ub = problem["aux"]["hi"][:]
    n = length(c)

    m = Model()
    @variable(m, lb[i] <= x[i=1:n] <= ub[i])
    @NLobjective(m, Min, sum(c[i] * x[i] for i = 1:n))
    @constraint(m, A * x .== b)

    return m
end

function get_file_names(dir::String,ext::String)
  name_list = readdir(dir)
  cleaned_name_list = []
  for name_full in name_list
    try
      if name_full[(end-3):end] == ext
          name = name_full[1:(end-4)]
      end
      cleaned_name_list = [cleaned_name_list; name]
      catch(e)
        @show e
        println("ERROR in running " * name_full)
      end
  end

  return cleaned_name_list
end

function run_infeas_netlib(test_name, my_par, solve_func!)
    summary = Dict{String, problem_summary2}()
    dir = "../data/netlib-infeas"
    problem_list = get_file_names(dir,".mat")

    if_mkdir("../results/$test_name")
    if_mkdir("../results/$test_name/log")
    if_mkdir("../results/$test_name/jld")

    for problem_name in problem_list
        if problem_name != "lpi_cplex2"
          m = read_lp2(dir, problem_name)
          nlp = MathProgNLPModel(m);
          #if nlp.meta.ncon < 2000
          summary[problem_name] = solve_func!(nlp, problem_name, test_name, my_par)
          #end
        end
        write_summary("../results/$(test_name)/summary.txt", summary)
        save("../results/$(test_name)/summary.jld","summary",summary, "pars", my_par)
    end
end

my_par = Class_parameters()
run_infeas_netlib("one_phase/Jan4/infeas-netlib", my_par, one_phase_run_and_store)
run_infeas_netlib("ipopt/Jan4/infeas-netlib", my_par, ipopt_solve_run_and_store)


text:
#=
subplot_2d:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-12
=#

module SubPlot2D
    # packages
    using PyPlot: plt
#     using PyCall
#     @pyimport matplotlib.pyplot as plt

    # flag to switch show or not
    show_plot = true

    # methods
    function set_show_plot(flag)
        global show_plot = flag
    end

    function plot_subplot()
        x = range(-2, 2, length=100)
        fig = plt.figure()
        for i in 1:6
            ax = fig.add_subplot(2, 3, i)
            ax.plot(x, x.^i, "-", c="c", ms=10)
            ax.set_xlabel("X")
            ax.set_ylabel("Y")
            ax.set_title("\$y = x^$(i)\$")
            ax.grid(true)
        end
        fig.tight_layout()
    end

    function main()
        plot_subplot()
        if show_plot == true
            plt.show()
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .SubPlot2D
    SubPlot2D.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
function choose_data(shared)
    d = Dict()
    for s in shared.selectlist
        d[Symbol(s.name)] = [value.name for value in s.values if value.accepted]
    end
    for s in shared.selectvalues
        d[Symbol(s.name)] =
            (t -> (s.values[1].selected_value <= t <= s.values[2].selected_value))
    end
    choose_data(shared.df, d)
end

function choose_data(a,d)
    index = broadcast(t -> true, 1:length(a))
    for key in keys(d)
        combine!(index, d[key], columns(a, key))
    end
    return a[find(index)]
end

function combine!(index, func::Function, values)
    @. index = index & func(values)
end

function combine!{T}(index, els::AbstractArray{T,1}, values)
    func = t -> (string(t) in els)
    combine!(index, func, values)
end

function combine!(index, el, values)
    func = t -> (t == el)
    return combine!(index, func, values)
end


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 3, 2, 4, :x, Quadrilateral(8, 2), false),
  :properties => [1.0e6 0.3;],
  :x_coords => [0.0, 3.0, 6.0],
  :y_coords => [0.0, -3.0, -6.0, -9.0],
  :support => [
    ( 1, [0 1]),
    ( 6, [0 1]),
    ( 9, [0 1]),
    (14, [0 1]),
    (17, [0 1]),
    (22, [0 1]),
    (25, [0 0]),
    (26, [0 0]),
    (27, [0 0]),
    (28, [0 0]),
    ( 5, [0 1]),
    ( 8, [0 1]),
    (13, [0 1]),
    (16, [0 1]),
    (21, [0 1]),
    (24, [0 1]),
    (29, [0 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -0.5]),
    (2, [0.0 -2.0]),
    (3, [0.0 -0.5])
    ]
)

@time m, dis_df, fm_df = p51(data);

@test m.loads[0:36] ≈ [0.0,-5.310672021656476e-6,-4.211153422117119e-7,
  -5.04123153829146e-6,-7.221882550664876e-7,-3.3428571428571473e-6,
  -4.211153422117125e-7,-1.6444827474228391e-6,-1.3750422640578195e-6,
  -4.288483723850779e-6,3.774202388609178e-7,-2.7857142857142914e-6,
  -1.2829448475778024e-6,-3.2431230529276257e-6,2.7081477113201545e-7,
  -2.8731660139934707e-6,3.6699957275470236e-7,-2.2285714285714337e-6,
  2.708147711320155e-7,-1.5839768431493952e-6,-1.2140198042152427e-6,
  -2.217378390101892e-6,2.996313701739242e-7,-1.6714285714285747e-6,
  -1.125478752755258e-6,-1.3787735959184822e-6,1.3704853665438675e-7,
  -1.298707004971293e-6,1.912334941440869e-7,-1.1142857142857151e-6,
  1.3704853665438667e-7,-9.298644236001403e-7,-8.497978326529501e-7,
  -6.453529188458675e-7,1.1219944448620244e-7,-5.571428571428575e-7,
  -4.68932795439848e-7]


text:
module ArchGDAL

using Dates
using GDAL: GDAL
using GeoFormatTypes: GeoFormatTypes
using GeoInterface: GeoInterface
using Tables: Tables
using ImageCore: Normed, N0f8, N0f16, N0f32, ImageCore
using ColorTypes: ColorTypes
using CEnum

const GFT = GeoFormatTypes

include("constants.jl")
include("utils.jl")
include("types.jl")
include("driver.jl")
include("geotransform.jl")
include("spatialref.jl")
include("dataset.jl")
include("raster/rasterband.jl")
include("raster/rasterio.jl")
include("raster/array.jl")
include("raster/rasterattributetable.jl")
include("raster/colortable.jl")
include("raster/images.jl")
include("ogr/geometry.jl")
include("ogr/feature.jl")
include("ogr/featurelayer.jl")
include("ogr/featuredefn.jl")
include("ogr/fielddefn.jl")
include("ogr/styletable.jl")
include("utilities.jl")
include("context.jl")
include("base/iterators.jl")
include("base/display.jl")
include("tables.jl")
include("geointerface.jl")
include("convert.jl")

function __init__()
    return nothing
end

end # module


text:
using StaticArrays

points2D = [SVector(-1.0, -1.0), SVector(1.0, 0.5), SVector(0.5, 1.0) , SVector(0.2, 0.3) ]

bbox = BoundingBox(points2D)

@test bbox.halflength == 1.0
@test bbox.center == SVector(0.0,0.0)

points2D = [SVector(1.0, 2.0), SVector(2.0, 3.0), SVector(1.5, 2.0) , SVector(1.8, 2.8) ]

bbox = BoundingBox(points2D)

@test bbox.halflength == 0.5
@test bbox.center == SVector(1.5, 2.5)

@test getchildbox(bbox, 3) == BoundingBox(0.25, SVector(1.25, 2.25))

@test whichchildbox(bbox, SVector(1.0, 2.0)) == 3
@test whichchildbox(bbox, SVector(2.0, 3.0)) == 1
@test whichchildbox(bbox, SVector(1.5, 2.0)) == 4
@test whichchildbox(bbox, SVector(1.8, 2.8)) == 1

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
using ..Orientations
import ..VariableTemplates: flattened_named_tuple
using ..VariableTemplates

# Sometimes `NodalStack` returns local states
# that is `nothing`. Here, we return `nothing`
# to preserve the keys (e.g., `hyperdiff`)
# when misssing.
flattened_named_tuple(v::Nothing, ft::FlattenType = FlattenArr()) = nothing

"""
    single_stack_diagnostics(
        grid::DiscontinuousSpectralElementGrid,
        bl::BalanceLaw,
        t::Real,
        direction;
        kwargs...,
    )

# Arguments
 - `grid` the grid
 - `bl` the balance law
 - `t` time
 - `direction` direction
 - `kwargs` keyword arguments, passed to [`NodalStack`](@ref).

An array of nested NamedTuples, containing results of
 - `z` - altitude
 - `prog` - the prognostic state
 - `aux` - the auxiliary state
 - `∇flux` - the gradient-flux (diffusive) state
 - `hyperdiff` - the hyperdiffusive state

and all the nested NamedTuples, merged together,
from the `precompute` methods.
"""
function single_stack_diagnostics(
    grid::DiscontinuousSpectralElementGrid,
    bl::BalanceLaw,
    t::Real,
    direction;
    kwargs...,
)
    return [
        begin
            @unpack prog, aux, ∇flux, hyperdiff = local_states
            diffusive = ∇flux
            state = prog
            hyperdiffusive = hyperdiff

            _args_fx1 = (; state, aux, t, direction)
            _args_fx2 = (; state, aux, t, diffusive, hyperdiffusive)
            _args_src = (; state, aux, t, direction, diffusive)

            cache_fx1 = precompute(bl, _args_fx1, Flux{FirstOrder}())
            cache_fx2 = precompute(bl, _args_fx2, Flux{SecondOrder}())
            cache_src = precompute(bl, _args_src, Source())

            z = altitude(bl, aux)

            nt = (;
                z = altitude(bl, aux),
                prog = flattened_named_tuple(prog), # Vars -> flattened NamedTuples
                aux = flattened_named_tuple(aux), # Vars -> flattened NamedTuples
                ∇flux = flattened_named_tuple(∇flux), # Vars -> flattened NamedTuples
                hyperdiff = flattened_named_tuple(hyperdiff), # Vars -> flattened NamedTuples
                cache_fx1,
                cache_fx2,
                cache_src,
            )
            # Flatten top level:
            flattened_named_tuple(nt)
        end for local_states in NodalStack(bl, grid; kwargs...)
    ]
end


text:
@testset "Context" begin

@testset "Create and destroy" begin
    MLIR.IR.Context() do ctx
        @test MLIR.IR.num_registered_dialects(ctx) == 0
        @test MLIR.IR.num_loaded_dialects(ctx) == 1
    end
end

@testset "Create and register the standard dialect" begin
    MLIR.IR.Context() do ctx
        @test MLIR.IR.num_registered_dialects(ctx) == 0
        MLIR.IR.register_standard_dialect!(ctx)
        @test MLIR.IR.num_registered_dialects(ctx) == 1

        prev_loaded = MLIR.IR.num_loaded_dialects(ctx)
        MLIR.IR.load_standard_dialect!(ctx)
        @test MLIR.IR.num_loaded_dialects(ctx) > prev_loaded
    end
end

end


text:
################################################################################
#
#   NfOrd/Ideal/Prime.jl : Prime ideals in orders of absolute number fields
#
# This file is part of Hecke.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann
#  Copyright (C) 2015, 2016, 2017 Claus Fieker
#
################################################################################

export PrimeIdealsSet

doc"""
***
    isramified(O::NfOrd, p::Int) -> Bool

> Returns whether the integer $p$ is ramified in $\mathcal O$.
> It is assumed that $p$ is prime.
"""
function isramified(O::NfOrd, p::Union{Int, fmpz})
  @assert ismaximal_known(O) && ismaximal(O)

  return mod(discriminant(O), p) == 0
end

doc"""
***
    degree(P::NfOrdIdl) -> Int
> The inertia degree of the prime-ideal $P$.
"""
function degree(A::NfOrdIdl)
  @assert isprime(A)
  return A.splitting_type[2]
end

doc"""
***
    ramification_index(P::NfOrdIdl) -> Int
> The ramification index of the prime-ideal $P$.
"""
function ramification_index(A::NfOrdIdl)
  @assert isprime(A)
  return A.splitting_type[1]
end

doc"""
***
    splitting_type(P::NfOrdIdl) -> Int, Int
> The ramification index and inertia degree of the prime ideal $P$.
> First value is the ramificatio index, the second the degree of $P$.
"""
function splitting_type(A::NfOrdIdl)
  @assert isprime(A)
  return A.splitting_type
end

################################################################################
#
#  Prime decomposition
#
################################################################################
doc"""
    intersect_nonindex(f::Map, P::NfOrdIdl) -> NfOrdIdl
> Given a prime ideal $P$ in $K$ and the inclusion map $f:k \to K$
> of number fields, find the unique prime $p$ in $k$ below.
"""
function intersect_nonindex(f::Map, P::NfOrdIdl)
  @assert P.is_prime == 1
  #let g be minpoly of k, G minpoly of K and h in Qt the primitive
  #element of k in K (image of gen(k))
  #then
  #  g(h) = 0 mod G
  k = domain(f)
  K = codomain(f)
  G = K.pol
  Qx = parent(G)
  g = k.pol(gen(Qx))
  h = Qx(f(gen(k)))

  Fp, xp = PolynomialRing(ResidueRing(FlintZZ, Int(minimum(P)), cached=false), cached=false)
  gp = factor(Fp(g))
  hp = Fp(h)
  Gp = gcd(Fp(K(P.gen_two)), Fp(G))
  for (f, e) = gp.fac
    if f(hp) % Gp == 0
      Zk = maximal_order(k)
      p = ideal_from_poly(Zk, Int(minimum(P)), f, 1)
      return p
    end
  end
end

doc"""
    prime_decomposition_nonindex(f::Map, p::NfOrdIdl) -> Array{Tuple{NfOrdIdl, Int}, 1}
> Given a map $f: k\to K$ of number fields defined over $\mathbb Q$ and
> a prime ideal in the maximal order of $k$, find all prime ideals in
> the maximal order of $K$ above.
"""
function prime_decomposition_nonindex(f::Map, p::NfOrdIdl)
  @assert p.is_prime == 1
  k = domain(f)
  K = codomain(f)
  ZK = maximal_order(K)
  G = K.pol
  Qx = parent(G)

  Fp, xp = PolynomialRing(ResidueRing(FlintZZ, Int(minimum(p)), cached=false), cached=false)
  Gp = factor(gcd(Fp(f(K(p.gen_two))), Fp(G)))
  res = []
  Zk = maximal_order(k)
  for (f, e) = Gp.fac
    P = ideal_from_poly(ZK, Int(minimum(p)), f, 1)
    push!(res, (P, e))
  end
  return res
end

doc"""
    lift(K::AnticNumberField, f::nmod_poly) -> nf_elem
> Given a polynomial $f$ over a finite field, lift it to an element of the
> number field $K$. The lift if given by the eleemnt represented by the
> canonical lift of $f$ to a polynomial over the integers.
"""
function lift(K::AnticNumberField, f::nmod_poly)
  if degree(f)>=degree(K)
    f = mod(f, parent(f)(K.pol))
  end
  r = K()
  for i=0:f.length-1
    u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ptr{nmod_poly}, Int), &f, i)
    _num_setcoeff!(r, i, u)
  end
  return r
end

##TODO: make fmpz-safe!!!!
#return <p, lift(O, fi> in 2-element normal presentation given the data
function ideal_from_poly(O::NfOrd, p::Int, fi::nmod_poly, ei::Int)
  b = lift(nf(O), fi)
  idl = ideal(O, fmpz(p), O(b, false))
  idl.is_prime = 1
  idl.splitting_type = ei, degree(fi)
  idl.norm = FlintZZ(p)^degree(fi)
  idl.minimum = FlintZZ(p)

  # We have to do something to get 2-normal presentation:
  # if ramified or valuation val(b,P) == 1, (p,b)
  # is a P(p)-normal presentation
  # otherwise we need to take p+b
  # I SHOULD CHECK THAT THIS WORKS

  if !((mod(norm(b),(idl.norm)^2) != 0) || (ei > 1))
    idl.gen_two = idl.gen_two + O(p)
  end

  idl.gens_normal = p
  idl.gens_weakly_normal = true

  # Find an "anti-uniformizer" in case P is unramified
  # We don't call it anti-unfiformizer anymore

  #if ideal.splitting_type[1] == 1
  #  t = parent(f)(lift(Zx, divexact(fmodp, fi)))
  #  ideal.anti_uniformizer = O(K(t), false)
  #end

  if idl.splitting_type[2] == degree(O)
    # Prime number is inert, in particular principal
    idl.is_principal = 1
    idl.princ_gen = O(p)
  end
  return idl
end

doc"""
***
    prime_decomposition(O::NfOrd,
                        p::Integer,
                        degree_limit::Int = 0,
                        lower_limit::Int = 0) -> Array{Tuple{NfOrdIdl, Int}, 1}

> Returns an array of tuples $(\mathfrak p_i,e_i)$ such that $p \mathcal O$ is the product of
> the $\mathfrak p_i^{e_i}$ and $\mathfrak p_i \neq \mathfrak p_j$ for $i \neq j$.
>
> If `degree_limit` is a nonzero integer $k > 0$, then only those prime ideals
> $\mathfrak p$ with $\deg(\mathfrak p) \leq k$ will be returned.
> Similarly if `\lower_limit` is a nonzero integer $l > 0$, then only those prime ideals
> $\mathfrak p$ with $l \leq \deg(\mathfrak p)$ will be returned.
> Note that in this case it may happen that $p\mathcal O$ is not the product of the
> $\mathfrak p_i^{e_i}$.
"""
function prime_decomposition(O::NfAbsOrd{S, T}, p::Union{Integer, fmpz}, degree_limit::Int = 0, lower_limit::Int = 0, cached::Bool = true) where {S, T}
  if typeof(p) == fmpz && nbits(p) < 64
    return prime_decomposition(O, Int(p), degree_limit, lower_limit)
  end

  if mod(index(O),fmpz(p)) == 0 || !issimple(nf(O))
    if cached
      if haskey(O.index_div, fmpz(p))
        lp = O.index_div[fmpz(p)]
        z = Tuple{NfAbsOrdIdl{S, T}, Int}[]
        for (Q, e) in lp
          if degree(Q) <= degree_limit
            push!(z, (Q, e))
          end
        end
        return z
      end
    end
    return prime_dec_index(O, p, degree_limit, lower_limit)
  end
  return prime_dec_nonindex(O, p, degree_limit, lower_limit)
end

function _fac_and_lift(f::fmpz_poly, p, degree_limit, lower_limit)
  Zx = parent(f)
  Zmodpx = PolynomialRing(ResidueRing(FlintZZ, p, cached = false), "y", cached = false)[1]
  fmodp = Zmodpx(f)
  fac = factor(fmodp)
  lifted_fac = Array{Tuple{fmpz_poly, Int}, 1}()
  for (k, v) in fac
    if degree(k) <= degree_limit && degree(k) >= lower_limit
      push!(lifted_fac, (lift(Zx, k), v))
    end
  end
  return lifted_fac
end

function prime_dec_nonindex(O::NfOrd, p::Union{Integer, fmpz}, degree_limit::Int = 0, lower_limit::Int = 0)
  K = nf(O)
  f = K.pol
  R = parent(f)
  Zx, x = PolynomialRing(FlintIntegerRing(),"x")
  Zf = Zx(f)

  if degree_limit == 0
    degree_limit = degree(K)
  end

  fac = _fac_and_lift(Zf, p, degree_limit, lower_limit)

  result = Array{Tuple{ideal_type(O),Int}}(length(fac))

  for k in 1:length(fac)
    fi = fac[k][1]
    ei = fac[k][2]
    #ideal = ideal_from_poly(O, p, fi, ei)
    t = parent(f)(fi)
    b = K(t)
    ideal = NfAbsOrdIdl(O)
    ideal.gen_one = p
    ideal.gen_two = O(b, false)
    ideal.is_prime = 1
    ideal.splitting_type = ei, degree(fi)
    ideal.norm = FlintZZ(p)^degree(fi)
    ideal.minimum = FlintZZ(p)

    # We have to do something to get 2-normal presentation:
    # if ramified or valuation val(b,P) == 1, (p,b)
    # is a P(p)-normal presentation
    # otherwise we need to take p+b
    # I SHOULD CHECK THAT THIS WORKS

    if !((mod(norm(b),(ideal.norm)^2) != 0) || (ei > 1))
      ideal.gen_two = ideal.gen_two + O(p)
    end

    ideal.gens_normal = p
    ideal.gens_weakly_normal = true

    # Find an "anti-uniformizer" in case P is unramified
    # We don't call it anti-unfiformizer anymore

    #if ideal.splitting_type[1] == 1
    #  t = parent(f)(lift(Zx, divexact(fmodp, fi)))
    #  ideal.anti_uniformizer = O(K(t), false)
    #end

    if length(fac) == 1 && ideal.splitting_type[2] == degree(f)
      # Prime number is inert, in particular principal
      ideal.is_principal = 1
      ideal.princ_gen = O(p)
    end
    result[k] =  (ideal, ei)
    k += 1
  end
  return result
end

function prime_dec_index(O::NfOrd, p::Union{Integer, fmpz}, degree_limit::Int = 0, lower_limit::Int = 0)
  if degree_limit == 0
    degree_limit = degree(O)
  end

  Ip = pradical(O, p)
  A, OtoA = AlgAss(O, Ip, p)
  AtoO = inv(OtoA)
  AA = split(A)

  basisO = basis(O)

  ideals = Vector{NfOrdIdl}()
  m = zero_matrix(FlintZZ, 1, degree(O))
  for (B, BtoA) in AA
    f = dim(B)
    idem = BtoA(B[1]) # Assumes that B == idem*A
    M = representation_matrix(idem)
    ker = left_kernel(M)
    N = basis_mat(Ip)
    for i = 1:length(ker)
      b = elem_in_basis(AtoO(A(ker[i])))
      for j = 1:degree(O)
        m[1, j] = b[j]
      end
      N = vcat(N, m)
    end
    N = sub(_hnf_modular_eldiv(N, fmpz(p), :lowerleft), rows(N) - degree(O) + 1:rows(N), 1:degree(O))
    P = ideal(O, N)
    P.norm = fmpz(p)^f
    P.splitting_type = (0, f)
    #
    fromOtosimplealgebra = Hecke._compose(inv(BtoA), OtoA)
    compute_residue_field_data!(P, fromOtosimplealgebra)
    #primB, minprimB, getcoordpowerbasis = _as_field(B)
    #@assert degree(minprimB) == f
    #P.prim_elem = AtoO(BtoA(primB))
    #P.min_poly_prim_elem = fmpz_poly(fmpz[FlintZZ(coeff(minprimB, i)) for i in 0:degree(minprimB)])
    #P.min_poly_prim_elem.parent = FmpzPolyRing(:$, false)
    #P.basis_in_prim = Vector{fmpz_mat}(degree(O))
    #for i in 1:degree(O)
    #  P.basis_in_prim[i] = zero_matrix(FlintZZ, 1, f)
    #  t = getcoordpowerbasis(BtoA\(OtoA(basisO[i])))
    #  for j in 1:f
    #    P.basis_in_prim[i][1, j] = FlintZZ(t[1, j])
    #  end
    #end
    push!(ideals, P)
  end

  result = Vector{Tuple{NfOrdIdl, Int}}()

  for j in 1:length(ideals)
    P = ideals[j]
    f = P.splitting_type[2]

    if f > degree_limit || f < lower_limit
      continue
    end

    # The following does not work if there is only one prime ideal
    if length(ideals) > 1 && (1-1/BigInt(p))^degree(O) < 0.1
      # This is roughly Algorithm 6.4 of Belabas' "Topics in comptutational algebraic
      # number theory".

      # Compute Vp = P_1 * ... * P_j-1 * P_j+1 * ... P_g

      B, BtoA = AA[j]
      J = ideal(O, AtoO(BtoA(B[1])))
      N = sub(_hnf_modular_eldiv(vcat(basis_mat(Ip), basis_mat(J)), fmpz(p), :lowerleft), degree(O) + 1:2*degree(O), 1:degree(O))
      Vp = ideal(O, N)

      u, v = idempotents(P, Vp)

      x = zero(parent(u))

      if !iszero(mod(norm(u), norm(P)*p))
        x = u
      elseif !iszero(mod(norm(u + p), norm(P)*p))
        x = u + p
      elseif !iszero(mod(norm(u - p), norm(P)*p))
        x = u - p
      else
        for i in 1:degree(O)
          if !iszero(mod(norm(v*basis(P)[i] + u), norm(P)*p))
            x = v*basis(P)[i] + u
          end
        end
      end

      @hassert :NfOrd 1 !iszero(x)
      @hassert :NfOrd 2 O*O(p) + O*x == P

      P.gen_one = p
      P.gen_two = x
      P.gens_normal = p
      P.gens_weakly_normal = 1
    else
      @vprint :NfOrd 1 "Chances for finding second generator: ~$((1-1/p))\n"
      _assure_weakly_normal_presentation(P)
      assure_2_normal(P)
    end

    e = Int(valuation(nf(O)(p), P))
    P.splitting_type = e, f
    P.is_prime = 1
    push!(result, (P, e))
  end
  #=
  if degree_limit >= degree(O)
    O.index_div[fmpz(p)] = result
  end
  =#
  return result
end

function uniformizer(P::NfOrdIdl)
  p = minimum(P)
  if P.gens_normal == p
    return P.gen_two
  else
    if p > 250
      r = 500  # should still have enough elements...
    else
      r = Int(div(p, 2))
    end
    z = rand(P, r)
    while true
      if !iszero(z) && valuation(z, P) == 1
        break
      end
      z = rand(P, r)
    end
    return z
  end
end

# Belabas p. 40
function anti_uniformizer(P::NfOrdIdl)
  if isdefined(P, :anti_uniformizer)
    return P.anti_uniformizer
  else
    p = minimum(P)
    M = representation_matrix(uniformizer(P))
    Mp = MatrixSpace(ResidueRing(FlintZZ, p, cached=false), rows(M), cols(M), false)(M)
    K = kernel(Mp)
    @assert length(K) > 0
    P.anti_uniformizer = elem_in_nf(order(P)(_lift(K[1])))//p
    return P.anti_uniformizer
  end
end

# Don't use the following function. It does not work for index divisors
# TH: Or does it?
function prime_decomposition_type(O::NfOrd, p::Integer)
  if (mod(discriminant(O), p)) != 0 && (mod(fmpz(index(O)), p) != 0)
    K = nf(O)
    f = K.pol
    R = parent(f)
    Zx, x = PolynomialRing(FlintZZ,"x", cached=false)
    Zf = Zx(f)
    fmodp = PolynomialRing(ResidueRing(FlintZZ,p, cached = false), "y", cached = false)[1](Zf)
    fac = factor_shape(fmodp)
    g = sum([ x for x in values(fac)])
    res = Array{Tuple{Int, Int}}(g)
    k = 1
    for (fi, ei) in fac
      for j in 1:ei
        res[k] = (fi, 1)
        k = k + 1
      end
    end
  else
    lp = prime_decomposition(O, p)
    res = Array{Tuple{Int, Int}}(length(lp))
    for i in 1:length(lp)
      res[i] = (lp[i][1].splitting_type[2], lp[i][1].splitting_type[1])
    end
  end
  return res
end

doc"""
***
    prime_ideals_up_to(O::NfOrd,
                       B::Int;
                       degree_limit::Int = 0) -> Array{NfOrdIdl, 1}

> Computes the prime ideals $\mathcal O$ with norm up to $B$.
>
> If `degree_limit` is a nonzero integer $k$, then prime ideals $\mathfrak p$
> with $\deg(\mathfrak p) > k$ will be discarded.
"""
function prime_ideals_up_to(O::NfOrd, B::Int;
                            complete::Bool = false,
                            degree_limit::Int = 0)
  p = 1
  r = NfOrdIdl[]
  while p < B
    p = next_prime(p)
    if p > B
      return r
    end
    if !complete
      deg_lim = Int(floor(log(B)/log(p)))
      if degree_limit >0
        deg_lim = min(deg_lim, degree_limit)
      end
    else
      deg_lim = 0
    end
    @vprint :ClassGroup 2 "decomposing $p ... (bound is $B, deg_lim $deg_lim)\n"
    li = prime_decomposition(O, p, deg_lim)
    for P in li
      push!(r, P[1])
    end
  end
  return r
end

doc"""
***
    prime_ideals_over(O::NfOrd,
                       lp::AbstractArray{Int, 1};
                       degree_limit::Int = 0) -> Array{NfOrdIdl, 1}

> Computes the prime ideals $\mathcal O$ over prime numbers in $lp$.
>
> If `degree_limit` is a nonzero integer $k$, then prime ideals $\mathfrak p$
> with $\deg(\mathfrak p) > k$ will be discarded.
"""
function prime_ideals_over(O::NfOrd,
                           lp::AbstractArray{T};
                           degree_limit::Int = 0) where T <: Union{fmpz, Integer}
  p = 1
  r = NfOrdIdl[]
  for p in lp
    @vprint :ClassGroup 2 "decomposing $p ... (deg_lim $deg_lim)"
    li = prime_decomposition(O, p, degree_limit)
    for P in li
      push!(r, P[1])
    end
  end
  return r
end

doc"""
***
    prime_ideals_up_to(O::NfOrd,
                       B::Int;
                       complete::Bool = false,
                       degree_limit::Int = 0,
                       F::Function,
                       bad::fmpz)

> Computes the prime ideals $\mathcal O$ with norm up to $B$.
>
> If `degree_limit` is a nonzero integer $k$, then prime ideals $\mathfrak p$
> with $\deg(\mathfrak p) > k$ will be discarded.
>
> The function $F$ must be a function on prime numbers not dividing `bad` such that
> $F(p) = \deg(\mathfrak p)$ for all prime ideals $\mathfrak p$ lying above $p$.
"""
function prime_ideals_up_to(O::NfOrd, B::Int, F::Function, bad::fmpz = discriminant(O);
                            complete::Bool = false,
                            degree_limit::Int = 0)
  p = 1
  r = NfOrdIdl[]
  while p < B
    p = next_prime(p)
    if p > B
      return r
    end
    if !complete
      deg_lim = flog(fmpz(B), p) # Int(floor(log(B)/log(p)))
      if degree_limit > 0
        deg_lim = min(deg_lim, degree_limit)
      end
    else
      deg_lim = 0
    end
    @vprint :ClassGroup 2 "decomposing $p ... (bound is $B)"
    if mod(bad, p) == 0
      li = prime_decomposition(O, p, deg_lim)
      for P in li
        push!(r, P[1])
      end
    else
      if F(p) <= deg_lim
        li = prime_decomposition(O, p, deg_lim)
        for P in li
          push!(r, P[1])
        end
      end
    end
  end
  return r
end

################################################################################
#
#  Coprime
#
################################################################################

#TODO: do sth. useful here!!!
function divides(A::NfOrdIdl, B::NfOrdIdl)
  minimum(A) % minimum(B) == 0 || return false
  return valuation(A, B) > 0
end

function coprime_base(A::Array{NfOrdIdl, 1}, p::fmpz)
  #consider A^2 B and A B: if we do gcd with the minimum, we get twice AB
  #so the coprime base is AB
  #however using the p-part of the norm, the coprime basis becomes A, B...
  if iseven(p)
    lp = prime_decomposition(order(A[1]), 2)
    Ap = [x[1] for x = lp if any(y->divides(y, x[1]) > 0, A)]
    a = remove(p, 2)[2]
    if !isone(a)
      Bp = coprime_base(A, a)
      return vcat(Ap, Bp)
    else
      return Ap
    end
  else
    Ap = [gcd(x, p^valuation(norm(x), p)) for x = A if minimum(x) % p == 0]
  end
  return coprime_base_steel(Ap)
end

doc"""
***
    coprime_base(A::Array{NfOrdIdl, 1}) -> Array{NfOrdIdl, 1}
    coprime_base(A::Array{NfOrdElem, 1}) -> Array{NfOrdIdl, 1}
> A coprime base for the (principal) ideals in $A$, ie. the returned array
> generated multiplicatively the same ideals as the input and are pairwise
> coprime.
"""
function coprime_base(A::Array{NfOrdIdl, 1})
  a = collect(Set(map(minimum, A)))
  a = coprime_base(a)
  C = Array{NfOrdIdl, 1}()

  for p = a
    if p == 1
      continue
    end
    cp = coprime_base(A, p)
    append!(C, cp)
  end
  return C
end

function coprime_base(A::Array{NfOrdElem, 1})
  O = parent(A[1])
  return coprime_base([ideal(O, x) for x = A])
end

function integral_split(A::NfOrdIdl)
  return A, ideal(Order(A), fmpz(1))
end

################################################################################
#
#  Factorization into prime ideals
#
################################################################################

#TODO: factoring type??? (with unit)
doc"""
***
    factor(A::NfOrdIdl) -> Dict{NfOrdIdl, Int}

> Computes the prime ideal factorization $A$ as a dictionary, the keys being
> the prime ideal divisors:
> If `lp = factor_dict(A)`, then `keys(lp)` are the prime ideal divisors of A
> and `lp[P]` is the `P`-adic valuation of `A` for all `P` in `keys(lp)`.
"""
factor(A::NfOrdIdl) = factor_dict(A)

function factor_dict(A::NfOrdIdl)
  ## this should be fixed
  lf = factor(minimum(A))
  lF = Dict{NfOrdIdl, Int}()
  n = norm(A)
  O = order(A)
  for (i, (p, v)) in enumerate(lf)
    lP = prime_decomposition(O, p)
    for P in lP
      v = valuation(A, P[1])
      if v != 0
        lF[P[1]] = v
        n = n//norm(P[1])^v
      end
      if n == 1
        return lF
      end
    end
  end
  return lF
end

################################################################################
#
#  Functions for index divisor splitting
#
################################################################################

mutable struct quoringalg{T} <: Ring
  base_order::NfOrd
  ideal::NfOrdIdl
  prime::T
  basis::Array{NfOrdElem, 1}

  function quoringalg(O::NfOrd, I::NfOrdIdl, p::T) where {T}

    z = new{T}()
    z.base_order = O
    z.ideal = I
    z.prime = p

    # compute a basis
    Rp = ResidueRing(FlintZZ, p, cached=false)
    Amodp = MatrixSpace(Rp, degree(O), degree(O), false)(basis_mat(I))
    Amodp = vcat(Amodp, zero_matrix(Rp, 1, degree(O)))
    Amodp[1,1] = 1
    Amodp = sub(Amodp, 1:degree(O), 1:degree(O))

    r, B = _rref(Amodp)
    C = zero_matrix(Rp, degree(O)-r, degree(O))
    BB = Array{NfOrdElem}(degree(O) - r)
    pivots = Array{Int}(0)
#    # get he pivots of B
    for i in 1:r
      for j in 1:degree(O)
        if !iszero(B[i,j])
          push!(pivots, j)
          break
        end
      end
    end
    i = 1
    k = 1
    while i <= degree(O)-r
      for j in k:degree(O)
        if !in(j, pivots)
          BB[i] = basis(O)[j]
          C[i,j] = 1
          k = j + 1
          i = i + 1
          break
        end
      end
    end
    insert!(BB, 1, basis(O)[1])
    z.basis = BB
    return z
  end
end

mutable struct quoelem
  parent::quoringalg
  elem::NfOrdElem
  coord::Array{fmpz, 1}

  function quoelem(R::quoringalg, x::NfOrdElem)
    z = new()
    z.parent = R
    z.elem = x

    return z
  end
end

function _kernel_of_frobenius(R::quoringalg)
  O = R.base_order
  BB = R.basis
  p = R.prime
  Rp = ResidueRing(FlintZZ, R.prime, cached=false)
  C = zero_matrix(Rp, length(BB)+1, degree(O))
  D = zero_matrix(Rp, length(BB), degree(O))

  Q, mQ = quo(O, R.ideal)

  function g(x)
    return mQ\(mQ(x)^p)
  end

  for i in 1:length(BB)
    A = elem_in_basis(mod(g(BB[i]) - BB[i], R.ideal))
    for j in 1:degree(O)
      D[i,j] = A[j]
    end
  end

  DD = NfOrdElem[ dot(BB, _lift(r)) for r in kernel(D) ]

  return [ quoelem(R, r) for r in DD ]
end

function _lift(T::Array{Generic.Res{fmpz}, 1})
  return [ z.data for z in T ]
end

function _lift(T::Array{Nemo.nmod, 1})
  return [ fmpz(z.data) for z in T ]
end

function *(x::quoelem, y::quoelem)
  z = mod(x.elem * y.elem, x.parent.ideal)
  return quoelem(x.parent, z)
end

function ^(x::quoelem, y::Int)
  z = mod(x.elem^y, x.parent.ideal)
  return quoelem(x.parent, z)
end

function ^(x::quoelem, y::Union{Integer, fmpz})
  # Do something stupid
  R, m = quo(x.parent.base_order, x.parent.ideal)
  return quoelem(x.parent, (m\(m(x.elem)^y)))
end

function ==(x::quoelem, y::quoelem)
  z = mod(x.elem - y.elem, x.parent.ideal)
  return zero(parent(z)) == z
end

function minpoly(x::quoelem)
  O = x.parent.base_order
  p = x.parent.prime

  Rp = ResidueRing(FlintZZ, p, cached=false)
  A = zero_matrix(Rp, 0, degree(O))
  B = zero_matrix(Rp, 1, degree(O))

  for i in 0:degree(O)
    ar = elem_in_basis((x^i).elem)
    for j in 1:degree(O)
      B[1, j] = ar[j]
    end
    A = vcat(A, B)
    K = kernel(A)
    if length(K) > 0
      @assert length(K) == 1
      f = PolynomialRing(Rp, "x", cached=false)[1](K[1])
      return f
    end
  end
  error("cannot find minpoly")
end

function split(R::quoringalg)
  if length(R.basis) == 1
    return [ R ]
  end
  K = _kernel_of_frobenius(R)
  O = R.base_order
  p = R.prime

  k = length(K)

  if k == 1
    # the algebra is a field over F_p
    # the ideal Ip is a prime ideal!
    return [ R ]
  end

  maxit = 1

  while true
    maxit = maxit + 1
    r = rand(0:p-1, length(K))

    x = quoelem(R, dot([ x.elem for x in K], r))

    if mod((x^p).elem, R.ideal) != mod(x.elem, R.ideal)
      #println("element^p: $(mod((x^p).elem, R.ideal))")
      #println("element: $(mod(x.elem, R.ideal))")
      #println(R.ideal.basis_mat)
      #println(K)
      error("Strange")
    end

    f = minpoly(x)

    if degree(f) < 2
      continue
    end
    @assert  issquarefree(f)

#    # By theory, all factors should have degree 1 # exploit this if p is small!
    fac = factor(f)
    F = first(keys(fac.fac))
    @assert length(fac) == degree(f)
    H = divexact(f,F)
    E, U, V = gcdx(F, H)
    @assert E == 1
    H = U*F;
    idem = O(coeff(H,0).data)
    for i in 1:degree(H)
      idem = idem + coeff(H,i).data*x.elem^i
    end

    I1 = R.ideal + ideal(O, idem)
    I2 = R.ideal + ideal(O, O(1)-idem)

    return vcat(split(quoringalg(O, I1, p)), split(quoringalg(O, I2, p)))
    break
  end
end

################################################################################
#
#  Primality testing
#
################################################################################

doc"""
***
    isprime_known(A::NfOrdIdl) -> Bool

> Returns whether $A$ knows if it is prime.
"""
function isprime_known(A::NfOrdIdl)
  return A.is_prime != 0
end

doc"""
***
    isprime(A::NfOrdIdl) -> Bool

> Returns whether $A$ is a prime ideal.
"""
function isprime(A::NfOrdIdl)
  if isprime_known(A)
    return A.is_prime == 1
  elseif minimum(A) == 0
    A.is_prime = 2
    return false
  end

  (n, p) = ispower(norm(A, Val{false}))

  if !isprime(p)
    A.is_prime = 2
    return false
  end

  p > 2^62 && error("Not implemented (yet)")

  lp = prime_decomposition(order(A), Int(p))

  for (P, f) in lp
    e = valuation(A, P)
    if e == 1 && n == degree(P)
      A.is_prime = 1
      return true
    elseif e == 0
      continue
    else
      A.is_prime = 2
      return false
    end
  end

  error("Something wrong in isprime")
end

################################################################################
#
#  Valuation
#
################################################################################

# CF:
# Classical algorithm of Cohen, but take a valuation element with smaller (?)
# coefficients. Core idea is that the valuation elementt is, originally, den*gen_two(p)^-1
# where gen_two(p) is "small". Actually, we don't care about gen_two, we
# need gen_two^-1 to be small, hence this version.
function val_func_no_index(p::NfOrdIdl)
  P = p.gen_one
  K = nf(order(p))
  pi = inv(p)
  d = denominator(K(pi.num.gen_two))
  @assert gcd(d, P) == 1
  e = K(pi.num.gen_two)*d
  M = zero_matrix(FlintZZ, 1, degree(K))
  elem_to_mat_row!(M, 1, d, e)
  @assert d == 1
  P2 = P^2
  P22 = div(P2, 2)
  for i=1:degree(K)
    x = M[1,i] % P2
    if x>P22
      x -= P2
    end
    M[1,i] = x
  end
  e = elem_from_mat_row(K, M, 1, P)
  # e is still a valuation element, but with smaller coefficients.
  return function(x::nf_elem)
    v = 0
    d = denominator(x)
    x *= d
    x = x*e
    while denominator(x) % P != 0
      v += 1
      mul!(x, x, e)
    end
    return v-valuation(d, P)*p.splitting_type[1]
  end
end

# CF:
# The idea is that valuations are mostly small, eg. in the class group
# algorithm. So this version computes the completion and the embedding into it
# at small precision and can thus compute (small) valuation at the effective
# cost of an mod(nmod_poly, nmod_poly) operation.
# Isn't it nice?
function valuation(a::UInt, b::UInt)
  return ccall((:n_remove, :libflint), Int, (Ref{UInt}, UInt), a, b)
end

#=
function valuation(a::UInt, b::UInt, bi::Cdouble)
  return ccall((:n_remove2_precomp, :libflint), Int, (Ref{UInt}, UInt, Cdouble), a, b, bi)
end
=#

function val_func_no_index_small(p::NfOrdIdl)
  P = p.gen_one
  @assert P <= typemax(UInt)
  K = nf(order(p))
  Rx = PolynomialRing(ResidueRing(FlintZZ, UInt(P), cached=false), cached=false)[1]
  Zx = PolynomialRing(FlintZZ)[1]
  g = Rx(p.gen_two.elem_in_nf)
  f = Rx(K.pol)
  g = gcd!(g, g, f)
  g = lift(Zx, g)
  k = flog(fmpz(typemax(UInt)), P)
  g = hensel_lift(Zx(K.pol), g, P, k)
  Sx = PolynomialRing(ResidueRing(FlintZZ, UInt(P)^k, cached=false), cached=false)[1]
  g = Sx(g)
  h = Sx()
  uP = UInt(P)
  return function(x::nf_elem)
    d = denominator(x)
    nf_elem_to_nmod_poly!(h, x, false) # ignores the denominator
    h = rem!(h, h, g)
    c = Nemo.coeff_raw(h, 0)
    v = c==0 ? typemax(Int) : valuation(c, uP)
    for i=1:degree(h)
      c = Nemo.coeff_raw(h, i)
      v = min(v, c==0 ? typemax(Int) : valuation(c, uP))
    end
    return v-valuation(d, P)
  end
end

function val_func_index(p::NfOrdIdl)
  # We are in the index divisor case. In larger examples, a lot of
  # time is spent computing denominators of order elements.
  # By using the representation matrix to multiply, we can stay in the order
  # and still be fast (faster even than in field).

  pi = inv(p)
  M = representation_matrix(pi.num.gen_two)
  O = order(p)
  P = p.gen_one
  return function(x::nf_elem)
    v = 0
    d, x_mat = integral_split(x, O)
    Nemo.mul!(x_mat, x_mat, M)
    while gcd(content(x_mat), P) == P  # should divide and test in place
      divexact!(x_mat, x_mat, P)
      Nemo.mul!(x_mat, x_mat, M)
      v += 1
    end
    return v-valuation(d, P)*p.splitting_type[1]
  end
end

doc"""
***
    valuation(a::nf_elem, p::NfOrdIdl) -> fmpz
    valuation(a::NfOrdElem, p::NfOrdIdl) -> fmpz
    valuation(a::fmpz, p::NfOrdIdl) -> fmpz

> Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
> such that $a$ is contained in $\mathfrak p^i$.
"""
function valuation(a::nf_elem, p::NfOrdIdl)
  @hassert :NfOrd 0 !iszero(a)
  #assert(a !=0) # can't handle infinity yet
  if isdefined(p, :valuation)
    return p.valuation(a)::Int
  end
  O = order(p)
  P = p.gen_one

  # for generic ideals
  if p.splitting_type[2] == 0
    #global bad_ideal = p
    p.valuation = function(a::nf_elem)
      d = denominator(a, O)
      x = O(d*a)
      return valuation_naive(O(x), p)::Int - valuation_naive(O(d), p)::Int
    end
    return p.valuation(a)::Int
  end

  if p.splitting_type[1]*p.splitting_type[2] == degree(O)
    p.valuation = function(a::nf_elem)
      return divexact(valuation(norm(a), P)[1], p.splitting_type[2])::Int
    end
  elseif mod(index(O),P) != 0 && p.splitting_type[1] == 1
    if p.gen_one^2 <= typemax(UInt)
      f1 = val_func_no_index_small(p)
      f2 = val_func_no_index(p)
      p.valuation = function(x::nf_elem)
        v = f1(x)
        if v > 100  # can happen ONLY if the precision in the .._small function
                    # was too small.
          return f2(x)::Int
        else
          return v::Int
        end
      end
    else
      p.valuation = val_func_no_index(p)
    end
  else
    p.valuation = val_func_index(p)
  end

  return p.valuation(a)::Int
end

doc"""
***
    valuation(a::nf_elem, p::NfOrdIdl) -> fmpz
    valuation(a::NfOrdElem, p::NfOrdIdl) -> fmpz
    valuation(a::fmpz, p::NfOrdIdl) -> fmpz

> Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
> such that $a$ is contained in $\mathfrak p^i$.
"""
valuation(a::NfOrdElem, p::NfOrdIdl) = valuation(a.elem_in_nf, p)

doc"""
***
    valuation(a::nf_elem, p::NfOrdIdl) -> fmpz
    valuation(a::NfOrdElem, p::NfOrdIdl) -> fmpz
    valuation(a::fmpz, p::NfOrdIdl) -> fmpz

> Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
> such that $a$ is contained in $\mathfrak p^i$.
"""
function valuation(a::fmpz, p::NfOrdIdl)
  if p.splitting_type[1] == 0
    return valuation_naive(order(p)(a), p)
  end
  P = p.gen_one
  return valuation(a, P)* p.splitting_type[1]
end
valuation(a::Integer, p::NfOrdIdl) = valuation(fmpz(a), p)

#TODO: some more intelligence here...
function valuation_naive(A::NfOrdIdl, B::NfOrdIdl)
  @assert !isone(B)
  Bi = inv(B)
  i = 0
  C = simplify(A* Bi)
  while denominator(C) == 1
    C = simplify(Bi*C)
    i += 1
  end
  return i
end

#TODO: some more intelligence here...
#      in non-maximal orders, interesting ideals cannot be inverted
#      maybe this needs to be checked...
function valuation_naive(x::NfOrdElem, B::NfOrdIdl)
  @assert !isone(B)
  i = 0
  C = B
  while x in C
    i += 1
    C *= B
  end
  return i
end

doc"""
***
    valuation(A::NfOrdIdl, p::NfOrdIdl) -> fmpz

> Computes the $\mathfrak p$-adic valuation of $A$, that is, the largest $i$
> such that $A$ is contained in $\mathfrak p^i$.
"""
function valuation(A::NfOrdIdl, p::NfOrdIdl)
  _assure_weakly_normal_presentation(A)
  if !isdefined(p, :splitting_type) || p.splitting_type[1] == 0 #ie. if p is non-prime...
    return valuation_naive(A, p)
  end
  return min(valuation(A.gen_one, p), valuation(elem_in_nf(A.gen_two), p))
end

################################################################################
#
#   Prime ideals iterator
#
################################################################################

mutable struct PrimeIdealsSet
  order::NfOrd
  from::fmpz
  to::fmpz
  primes::PrimesSet{fmpz}
  currentprime::fmpz
  currentindex::Int
  decomposition::Array{Tuple{NfOrdIdl, Int}, 1}
  proof::Bool
  indexdivisors::Bool
  ramified::Bool
  iscoprimeto::Bool
  coprimeto::NfOrdIdl
  degreebound::Int
  unbound::Bool

  function PrimeIdealsSet(O::NfOrd)
    z = new()
    z.order = O
    z.proof = false
    z.indexdivisors = true
    z.ramified = true
    z.unbound = false
    z.degreebound =  degree(O)
    z.iscoprimeto = false
    return z
  end
end

doc"""
***
    PrimeIdealsSet(O::NfOrd, f, t; proof = false,
                                   indexdivisors = true,
                                   ramified = true,
                                   degreebound = degree(O),
                                   coprimeto = false)

Returns an iterable object $S$ representing the prime ideals $\mathfrak p$
of $\mathcal O$ with $f \leq \min(\mathfrak p) \leq t$.

The optional arguments can be used to exclude index divisors, ramified prime
ideals and to include only prime ideals with degree less or equal than
`degreebound` and which are coprime to `coprimeto`.

If $t=-1$, then the upper bound is infinite.

If `coprimeto` is supplied, it must be either an integer, an element of $\mathcal O$,
or a non-zero ideal of $\mathcal O$.
"""
function PrimeIdealsSet(O::NfOrd, from::T, to::S;
                       proof::Bool = false,
                       indexdivisors::Bool = true,
                       ramified::Bool = true,
                       degreebound::Int = degree(O),
                       coprimeto = false) where {T <: Union{fmpz, Integer},
                                                 S <: Union{fmpz, Integer}}
  from < 0 && error("Lower bound must be non-negative")
  to < -1 && error("Upper bound must be non-negative or -1")

  z = PrimeIdealsSet(O)
  z.from = fmpz(from)
  z.to = fmpz(to)
  z.primes = PrimesSet(z.from, z.to)
  if to == -1
    z.unbound = true
  end
  z.proof = proof
  z.indexdivisors = indexdivisors
  z.ramified = ramified
  z.degreebound = degreebound
  if !(coprimeto isa Bool)
    if coprimeto isa NfOrdIdl
      z.coprimeto = coprimeto
    elseif coprimeto isa Union{Integer, fmpz, NfOrdElem}
      z.coprimeto = ideal(O, coprimeto)
    else
      error("Coprime argument of wrong type ($(typeof(coprimeto)))")
    end
    z.iscoprimeto = true
  end
  return z
end

function Base.start(S::PrimeIdealsSet)
  O = S.order
  pstate = start(S.primes)
  found_prime = false
  while !found_prime
    (p, pstate) = next(S.primes, pstate)
    if !S.indexdivisors && isindex_divisor(O, p)
      continue
    end
    lP = prime_decomposition(O, p)
    j = -1
    for i in 1:length(lP)
      e = lP[i][2]
      if !S.ramified && e > 1
        continue
      end
      P = lP[i][1]
      if P.splitting_type[2] > S.degreebound
        continue
      end
      if S.iscoprimeto && !iscoprime(P, S.coprimeto)
        continue
      end
      j = i
      break
    end
    if j != -1
      S.decomposition = lP
      S.currentprime = p
      S.currentindex = j
      return (p, j)
    end
  end
end

function Base.next(S::PrimeIdealsSet, x)
  pstate = x[1]
  j = x[2]
  Q = S.decomposition[j][1] # This we want to return
  newindex = -1
  lP = S.decomposition
  O = S.order

  # Find the next prime ideal in the current decomposition
  for i in (j + 1):length(S.decomposition)
    e = lP[i][2]
    if !S.ramified && e > 1
      continue
    end
    P = lP[i][1]
    if P.splitting_type[2] > S.degreebound
      continue
    end
    newindex = i
    break
  end

  if newindex != -1
    return Q, (pstate, newindex)
  else
    # We have to change the prime
    found_prime = false
    while !found_prime
      (p, pstate) = next(S.primes, pstate)
      if !S.indexdivisors && isindex_divisor(O, pstate)
        continue
      end
      lP = prime_decomposition(O, pstate)
      j = -1
      for i in 1:length(lP)
        e = lP[i][2]
        if !S.ramified && e > 1
          continue
        end
        P = lP[i][1]
        if P.splitting_type[2] > S.degreebound
          continue
        end
        if S.iscoprimeto && !iscoprime(P, S.coprimeto)
          continue
        end
        j = i
        break
      end
      if j != -1
        S.decomposition = lP
        S.currentprime = p
        S.currentindex = j
        return Q, (pstate, j)
      end
    end
  end
end

function Base.done(S::PrimeIdealsSet, x)
  pstate = x[1]
  index = x[2]
  return !S.unbound && pstate > S.to
end

Base.eltype(::PrimeIdealsSet) = NfOrdIdl

Base.iteratorsize(::Type{PrimeIdealsSet}) = Base.SizeUnknown()


text:
using YaoExperiment
using Test

@testset "repl" begin
    include("repl.jl")
end

@testset "quil" begin
    include("quil.jl")
end


text:
#=
run natvsext benchmarks
see natvsext/README.md
=#

using Test
using Printf
import DataFrames
import CSV
import Hypatia
using MosekTools
using ECOS

include(joinpath(@__DIR__, "../../examples/Examples.jl"))
using Main.Examples

using Distributed
include(joinpath(@__DIR__, "spawn.jl"))

# path to write results DataFrame to CSV
results_path = joinpath(mkpath(joinpath(@__DIR__, "raw")), "bench.csv")

# option to keep setting up larger models, only if solver is Hypatia,
# even if last solve was killed
# setup_model_anyway = true
setup_model_anyway = false

verbose = true # make solvers print output
# verbose = false

iter_limit = 250
num_threads = 16 # number of threads for BLAS and Julia processes running instances
free_memory_limit = 8 * 2^30 # keep at least X GB of RAM available
optimizer_time_limit = 1800
setup_time_limit = 1.2 * optimizer_time_limit
check_time_limit = 1.2 * optimizer_time_limit
tol_loose = 1e-7
tol_tight = 1e-3 * tol_loose

hyp_solver = ("Hypatia", Hypatia.Optimizer, (
    verbose = verbose,
    iter_limit = iter_limit,
    time_limit = optimizer_time_limit,
    tol_abs_opt = tol_tight,
    tol_rel_opt = tol_loose,
    tol_feas = tol_loose,
    tol_infeas = tol_tight,
    init_use_indirect = true, # skips dual equalities preprocessing
    use_dense_model = true,
    ))

mosek_solver = ("Mosek", Mosek.Optimizer, (
    QUIET = !verbose,
    MSK_IPAR_INTPNT_MAX_ITERATIONS = iter_limit,
    MSK_IPAR_NUM_THREADS = num_threads,
    MSK_IPAR_OPTIMIZER = Mosek.MSK_OPTIMIZER_CONIC,
    MSK_IPAR_INTPNT_BASIS = Mosek.MSK_BI_NEVER, # no basis identification for LP
    MSK_DPAR_OPTIMIZER_MAX_TIME = optimizer_time_limit,
    MSK_DPAR_INTPNT_CO_TOL_REL_GAP = tol_loose,
    MSK_DPAR_INTPNT_CO_TOL_PFEAS = tol_loose,
    MSK_DPAR_INTPNT_CO_TOL_DFEAS = tol_loose,
    MSK_DPAR_INTPNT_CO_TOL_INFEAS = tol_tight,
    ))

ecos_solver = ("ECOS", ECOS.Optimizer, (
    verbose = 2 * verbose,
    maxit = iter_limit, # no time limit option available
    abstol = tol_tight,
    reltol = tol_loose,
    feastol = tol_loose,
    ))

# instance sets and solvers to run
inst_sets = [
    #= natural formulations paper =#
    # ("nat", hyp_solver),
    # ("ext", hyp_solver),
    # ("ext", mosek_solver),
    # ("extEP", hyp_solver), # ExpPSD extender
    # ("extSEP", hyp_solver), # SOCExpPSD extender
    # ("extEP", mosek_solver), # ExpPSD extender
    # ("extSEP", mosek_solver), # SOCExpPSD extender
    #= spectral function cones paper =#
    ("nat", hyp_solver),
    ("ext", hyp_solver),
    ("ext", mosek_solver),
    # for nonparametricdistr:
    ("vecext", hyp_solver),
    ("vecext", mosek_solver),
    ("vecext", ecos_solver),
    # for covarianceest
    ("logdet", hyp_solver),
    ("sepspec", hyp_solver),
    ("direct", hyp_solver),
    ("direct", mosek_solver),
    ]

# models to run
JuMP_examples = [
    #= natural formulations paper =#
    # "densityest",
    # "doptimaldesign",
    # "matrixcompletion",
    # "matrixregression",
    # "polymin",
    # "portfolio",
    # "shapeconregr",
    #= WSOS cones paper =#
    # "nearestpolymat",
    # "polynorm",
    #= spectral function cones paper =#
    "centralpolymat",
    "classicalquantum",
    "covarianceest",
    "experimentdesign",
    "nonparametricdistr",
    ]

interrupt()
@assert nprocs() == 1
println()

print_memory() = println("free memory (GB): ", Float64(Sys.free_memory()) / 2^30)
print_memory()

extender_name(ext::Nothing) = missing
extender_name(ext::Symbol) = string(ext)

println("\nstarting benchmark runs")
time_all = time()

@testset "examples tests" begin
perf = Examples.setup_benchmark_dataframe()
CSV.write(results_path, perf)

@testset "$ex" for ex in JuMP_examples
(ex_type, ex_insts) = Examples.get_benchmark_instances("JuMP", ex)
@testset "$inst_set, $(solver[1])" for (inst_set, solver) in inst_sets
    haskey(ex_insts, inst_set) || continue
    (extender, inst_subsets) = ex_insts[inst_set]
    isempty(inst_subsets) && continue

    info_perf = (; inst_set, :example => ex, :model_type => "JuMP",
        :real_T => Float64, :solver => solver[1],
        :extender => extender_name(extender))
    str = "$ex $inst_set $(solver[1])"
    println("\nstarting $str")

    for inst_subset in inst_subsets
        solve = true
        # first instance is only used for compilation
        compile_inst = inst_subset[1]
        for (inst_num, inst_data) in enumerate(inst_subset[2:end])
            println("\nstarting $str $inst_num $inst_data")
            flush(stdout); flush(stderr)

            total_time = @elapsed (setup_killed, check_killed, run_perf) =
                spawn_instance(ex, ex_type{Float64}, compile_inst,
                inst_data, extender, solver, solve, num_threads)

            new_perf = (; info_perf..., run_perf..., total_time,
                inst_num, inst_data)
            Examples.write_perf(perf, results_path, new_perf)
            @printf("%8.2e seconds\n", total_time)
            flush(stdout); flush(stderr)

            setup_killed && break
            if check_killed
                if setup_model_anyway && (solver[1] == "Hypatia")
                    solve = false
                else
                    break
                end
            end
        end
    end
end
end

flush(stdout); flush(stderr)
println("\n")
DataFrames.show(perf, allrows = true, allcols = true)
println("\n")
flush(stdout); flush(stderr)
end

interrupt()
@printf("\nbenchmarks total time: %8.2e seconds\n\n", time() - time_all)
;


text:
include("../src/UnicodeExtras.jl")
using UnicodeExtras
using Base.Test

noel1 = UnicodeText("noe\u0308l")
noel2 = UnicodeText("noël")
@test noel1 == noel2
@test length(noel1) == length(noel2) == 4
@test noel1[1:3] == "noe\u0308"
@test noel1[3] == "e\u0308"

set_locale("de")
@test UnicodeText("Köpfe") < UnicodeText("Kypper")
@test uppercase("testingß") == "TESTINGSS"
@test uppercase(UnicodeText("testingß")) == "TESTINGSS"
set_locale("sv")
@test UnicodeText("Kypper") < UnicodeText("Köpfe")
set_locale("tr")
@test uppercase("testingß") == "TESTİNGSS"
@test uppercase(UnicodeText("testingß")) == "TESTİNGSS"

s = "Äpple"
b = encode(s, "iso-8859-1")
s2 = decode(b, "iso-8859-1")
@test s == s2
@test b == [196,112,112,108,101]


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module GMMData
include("load.jl")

export Wishart, GMMInput, GMMOutput, empty_gmm_output, load_gmm_input

struct Wishart
    gamma::Float64
    m::Int
end

struct GMMInput
    alphas::Matrix{Float64}
    means::Matrix{Float64}
    icfs::Matrix{Float64}
    x::Matrix{Float64}
    wishart::Wishart
end

mutable struct GMMOutput
    objective::Float64
    gradient::Vector{Float64}
end

empty_gmm_output() = GMMOutput(0.0, [])

function load_gmm_input(fn::AbstractString, replicate_point::Bool)::GMMInput
    fid = open(fn)
    lines = readlines(fid)
    close(fid)
    line=split(lines[1], " ")
    d = parse(Int, line[1])
    k = parse(Int, line[2])
    n = parse(Int, line[3])
    icf_sz = div(d * (d + 1), 2)
    off = 1

    alphas = zeros(Float64, 1, k)
    for i in 1:k
        alphas[i] = parse(Float64, lines[i + off])
    end
    off += k

    means = zeros(Float64, d, k)
    for ik in 1:k
        line = split(lines[ik + off], " ")
        for id in 1:d
            means[id, ik] = parse(Float64, line[id])
        end
    end
    off += k

    icf = zeros(Float64, icf_sz, k)
    for ik in 1:k
        line = split(lines[ik + off], " ")
        for i in 1:icf_sz
            icf[i, ik] = parse(Float64, line[i])
        end
    end
    off += k

    if replicate_point
        x_ = zeros(Float64, d, 1)
        line = split(lines[1 + off], " ")
        for id in 1:d
            x_[id] = parse(Float64, line[id])
        end
        x = repeat(x_, 1, n)
        off += 1
    else
        x = zeros(Float64, d, n)
        for ix in 1:n
            line = split(lines[ix + off], " ")
            for id in 1:d
                x[id, ix] = parse(Float64, line[id])
            end
        end
        off += n
    end
    line = split(lines[1 + off]," ")
    wishart = Wishart(parse(Float64, line[1]), parse(Int, line[2]))
    GMMInput(alphas, means, icf, x, wishart)
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
## test of basic widgets
using JGUI
using Base.Test

w = window()
g = vbox(w); push!(w, g)

lab = label(w, "label")
btn = button(w, "text")
ed = lineedit(w, "text")
## XXX txt = textedit(w, "more text")
cb = checkbox(w, true, "fred")
rad = radiogroup(w, map(string, [1:3]))
## XXX bgp = buttongroup(w, ["a", "b", "c"])
cmb = combobox(w,  ["a", "b", "c"])
## XXX cmb_ed = combobox(w,  ["a", "b", "c"], editable=true)
sl = slider(w, 1:100)
sp = spinbox(w, 1:100)
im = imageview(w, "/tmp/trash.png")

widgets = [lab, btn, ed, cb, rad, cmb, sl,sp, im]
append!(g, widgets)

raise(w)


text:
module StanModels

const stanmodels_src_path = @__DIR__

"""

# stanmodels_path

Relative path using the StanModels src/ directory.

### Example to get access to the scripts subdirectory of StanModels
```julia
stanmodels_path("..", "scripts")
```
"""
stanmodels_path(parts...) = normpath(joinpath(stanmodels_src_path, parts...))

using StanSample

include("utils/scale.jl")

export
  stanmodels_path

end # module


text:
using DifferentialEquations
using Statistics
using StatsBase
using Plots
using LinearAlgebra
using Distributions
using HDF5
using JLD2
using GaussianMixtures

using Logging
logger = SimpleLogger(stdout, Logging.Error) # Quiet GaussianMixtures output
global_logger(logger);

# Compute total mRNA at each timepoint in each cell (for MMI2-SSB state)
function reduce_mmi(sol)
    popsize = size(sol)[3]
    mT = sum(sol[2:5,:,:], dims=1)
    return reshape(mT[:,:,:], (size(sol)[2], popsize))
end

# Get mRNA 1 concentration at each timepoint in each cell (for repressilator or genetic oscillator state)
function reduce_repressilator(sim)
    return sim[1, :, :]
end
function reduce_geneosc(sim)
    return sim[1, :, :]
end

# Heatmap of state distributions at (final) timepoint across a range of signals
function perturbed_results_heatmap(results; reduction=reduce_mmi, time=nothing, signame="σR", statelabel="[Rtot]", maxrna=5, binwidth=0.25, kwargs...)
    sR_values = keys(results) |> collect |> sort
    bins = 0:binwidth:maxrna
    final_histograms = map(sR_values) do sR
        mTs = reduction(results[sR])
        if isnothing(time)
            states = mTs[end, :]
        else
            states = mTs[time, :]
        end
        fit(Histogram, states, bins).weights
    end
    matrix = hcat(final_histograms...)
    heatmap(sR_values, collect(bins)[begin:(end - 1)] .+ Float64(bins.step / 2), matrix;
        xlabel=signame, ylabel=statelabel, colorbar_title="Cells", xminorticks=10, tick_direction=:out, kwargs...)
end

# Separation (D) of two modes in a one-dimensional dataset
function bimode_distinguishability(data)
    return bimode_distinguishability(GMM(2, reshape(data, length(data), 1)))
end
function bimode_distinguishability(model::GMM)
    # "Separation of means relative to their widths"
    # "D > 2 is required for a clean separation between the modes"
    # --Oleg Y. Gnedin, http://www-personal.umich.edu/~ognedin/gmm/gmm_user_guide.pdf
    μ = means(model)
    return abs(μ[1] - μ[2]) / sqrt(sum(covars(model)) / 2)
end

# Bayesian Information Criterion for a GMM fit to a 1D dataset
function bic(model::GMM, data)
    # https://github.com/scikit-learn/scikit-learn/blob/0d378913b/sklearn/mixture/_gaussian_mixture.py#L821-L823
    average_log_likelihood = avll(model, reshape(data, length(data), 1))
    nparams = length(means(model)) * 3 - 1 # mean, variance, weight (but weights must sum to 1)
    return nparams * log(length(data)) - 2 * average_log_likelihood * length(data)
end

# ΔBIC, the advantage of a 2-Gaussian fit over a 1-Gaussian fit for a 1D dataset
function bic_bimodal_advantage(data)
    data_reshaped = reshape(data, length(data), 1)
    bic1 = bic(GMM(1, data_reshaped), data_reshaped)
    bic2 = bic(GMM(2, data_reshaped), data_reshaped)
    return bic1 - bic2 # Positive if bimodality allows a better model
end

# Time index at which a population timecourse reached equilibrium
function equilibriumindex(results, savedat=0.1, reduction=reduce_mmi; stabletime=100.0, checktime=5.0, maxrna=10, binwidth=0.25, mainspecies=2)
    timesteps = collect(axes(results)[2])
    bins = 0:binwidth:maxrna
    histograms = map(t -> fit(Histogram, results[mainspecies, t, :], bins).weights, timesteps)
    stableindex = stabletime / savedat |> round |> Int
    final = dropdims(mean(hcat(histograms[stableindex:end]...); dims=2); dims=2)
    rmsds = map(t -> sqrt(mean((histograms[t] .- final) .^ 2.0)), timesteps)
    maxrmsd = maximum(rmsds[stableindex:end])
    mTs = reduction(results)
    totalmeans = dropdims(mean(mTs, dims=2); dims=2)
    equilibriummean = mean(totalmeans[stableindex:end])
    equilibriummaxmean = maximum(totalmeans[stableindex:end])
    equilibriumminmean = minimum(totalmeans[stableindex:end])
    ok_rmsd = false
    ok_mean = false
    candidate = missing
    for index in 2:stableindex
        cleared = false
        if rmsds[index] > maxrmsd * 1.1
            ok_rmsd = false
            cleared = true
        end
        if !(equilibriumminmean < totalmeans[index] < equilibriummaxmean)
            ok_mean = false
            cleared = true
        end
        if cleared
            candidate = missing
            continue
        end
        if rmsds[index] < maxrmsd / 2
            ok_rmsd = true
        end
        if sign(totalmeans[index - 1] - equilibriummean) != sign(totalmeans[index] - equilibriummean)
            ok_mean = true
        end
        if ok_rmsd && ok_mean
            if ismissing(candidate)
                candidate = index
            elseif (index - candidate) * savedat >= checktime
                return candidate
            end
        end
    end
    return candidate
end

# Timecourse of population state mean and state distribution for a population
function timecourse_heatmap(results; startindex=first(axes(results)[2]), endindex=last(axes(results)[2]), savedat=0.1, stabletime=100.0, reduction=reduce_mmi, statename="[Rtot]", maxrna=10, binwidth=0.25, kwargs...)
    mTs = reduction(results)
    timesteps = collect(startindex:endindex)
    bins = 0:binwidth:maxrna
    step_histograms = map(timesteps) do time
        fit(Histogram, mTs[time, :], bins).weights
    end
    matrix = hcat(step_histograms...)
    diagram = heatmap(timesteps * savedat, collect(bins)[begin:(end - 1)] .+ Float64(bins.step / 2), matrix; xlabel="Time", ylabel=statename, colorbar_title="Cells", kwargs...)
    means = dropdims(mean(mTs, dims=2); dims=2)
    plot!(diagram, (timesteps .- 1) * savedat, means[timesteps], color=:white, legend=nothing)
    stableindex = stabletime / savedat |> round |> Int
    hline!(diagram, [mean(mTs[stableindex:end, :]), maximum(means[stableindex:end]), minimum(means[stableindex:end])], color=:gray, legend=nothing)
    equilibriumstart = equilibriumindex(results, savedat, reduction; stabletime)
    vline!(diagram, [(equilibriumstart - 1) * savedat], color=:green, legend=nothing)
    xaxis!(diagram, xminorticks=10, tick_direction=:out, xlims=(0, (endindex - 1) * savedat))
    yaxis!(diagram, ylims=(0, maxrna))
end

# Create a wrapper function that runs the original several times and produces the minimum (counters nondeterminism in GMM fit)
function takemin(wrap, runs)
    function wrapper(args...)
        results = map(_ -> wrap(args...), 1:runs)
        return minimum(results)
    end
    return wrapper
end

# ΔBIC and separation line plot(s) for a set of signals (and different noise levels and initial conditions)
function perturbed_reestablishment_subplot(record, title; reduction=reduce_mmi, separation_lim=Inf, linekws=Dict())
    sR_values = keys(record) |> collect |> sort
    endpoints = map(sR_values) do sR
        mTs = reduction(record[sR])
        mTs[end, :]
    end
    distinguishabilities = map(takemin(bimode_distinguishability, 4), endpoints)
    bicadvantages = map(takemin(bic_bimodal_advantage, 4), endpoints)
    titlekws = Dict()
    if ~isnothing(title)
        titlekws = Dict(:title => title, :titlefontsize => 11)
    end
    subplot = plot(sR_values, bicadvantages; legend=nothing, ylabel="ΔBIC", yguidefontcolor=:blue, titlekws..., linekws...)
    hline!(subplot, [0.0], color=:darkblue, legend=nothing)
    twin = twinx(subplot)
    plot!(twin, sR_values, distinguishabilities; color=:orange, xticks=:none, legend=nothing, ylabel="Separation", yguidefontcolor=:orange, linekws...)
    hline!(twin, [2.0], color=:brown, legend=nothing)
    yaxis!(twin, ylims=(0, separation_lim))
    return subplot
end
function perturbed_reestablishment_plot(records; noiselabel="Noise", kwargs...)
    subplots = map(noise -> perturbed_reestablishment_subplot(records[noise], "$noiselabel = $noise"; kwargs...), keys(records) |> collect |> sort)
    plot(subplots..., layout=(length(keys(records)), 1), size=(600, 400), right_margin=15*Plots.mm)
end
function perturbed_reestablishment_plot(records::AbstractVector; noiselabel="Noise", size=(1000, 400), kwargs...)
    columns = map(run -> map(noise ->
        perturbed_reestablishment_subplot(run[noise], isnothing(noiselabel) ? nothing : "$noiselabel = $noise"; kwargs...),
        keys(run) |> collect |> sort), records)
    subplots = hcat(columns...)
    plot(permutedims(subplots)..., layout=grid(length(columns[1]), length(columns)), size=size, left_margin=5*Plots.mm, right_margin=15*Plots.mm)
end

# Plot time to equilibrium for a set of signals (and different noise levels and initial conditions)
function perturbed_equilibrium_plot(records::Dict{T, Array{T, 3}}; savedat=0.1, kwargs...) where T <: Real
    sR_values = keys(records) |> collect |> sort
    equilibriumtimes = map(sR -> (equilibriumindex(records[sR]) - 1) * savedat, sR_values)
    plot(sR_values, equilibriumtimes, xlabel="sR", ylabel="Time", legend=nothing; kwargs...)
end
function perturbed_equilibrium_plot(records::Dict{T, Dict{T, Array{T, 3}}}; kwargs...) where T <: Real
    noises = keys(records) |> collect |> sort
    subplots = map(noise -> perturbed_equilibrium_plot(records[noise]; title="Noise = $noise", titlefontsize=11, kwargs...), noises)
    plot(subplots..., layout=(length(noises), 1))
end
function perturbed_equilibrium_plot(records::AbstractVector; kwargs...)
    columns = map(run -> map(noise -> perturbed_equilibrium_plot(run[noise]; title="Noise = $noise", titlefontsize=11, kwargs...), keys(run) |> collect |> sort), records)
    subplots = hcat(columns...)
    plot(permutedims(subplots)..., layout=grid(length(columns[1]), length(columns)), link=:all, size=(1000, 450), left_margin=5*Plots.mm, bottom_margin=5*Plots.mm)
end

# Print maximum and median equilibrium times for populations in a range of signals (missing or 1000 if ET is too long to determine)
function eqtimestats(desc, records, span)
    eqtimes = map(sR -> (equilibriumindex(records[sR]) - 1) * 0.1, span)
    filled = replace(eqtimes, missing => 1000.0)
    println(desc, "\tMax: ", maximum(eqtimes), "\t", "Avg: ", median(filled))
end

# Load a set of loose Python Gillespie results (filenames under "gillespie" directory, specified by regex) into a noise -> signal -> array dict to match Julia SDE outputs
function loadgillespies(re)
    results = Dict{Float64, Dict{Float64, Array{Float64, 3}}}()
    for filename in filter(f -> match(re, f) !== nothing, readdir("gillespie"))
        file = h5open("gillespie/$filename")
        noise = read(attributes(file)["V"])
        if noise ∉ keys(results)
            results[noise] = Dict{Float64, Array{Float64, 3}}()
        end
        signal = read(attributes(file)["s_m"])
        data = read(file["results"])
        results[noise][signal] = permutedims(data, reverse(1:ndims(data)))
        close(file)
    end
    return results
end