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

# Missing, missing and ismissing are defined in essentials.jl

show(io::IO, x::Missing) = print(io, "missing")

"""
    MissingException(msg)

Exception thrown when a [`missing`](@ref) value is encountered in a situation
where it is not supported. The error message, in the `msg` field
may provide more specific details.
"""
struct MissingException <: Exception
    msg::AbstractString
end

showerror(io::IO, ex::MissingException) =
    print(io, "MissingException: ", ex.msg)

promote_rule(::Type{Missing}, ::Type{T}) where {T} = Union{T, Missing}
promote_rule(::Type{Union{S,Missing}}, ::Type{T}) where {T,S} = Union{promote_type(T, S), Missing}
promote_rule(::Type{Any}, ::Type{T}) where {T} = Any
promote_rule(::Type{Any}, ::Type{Missing}) = Any
promote_rule(::Type{Missing}, ::Type{Any}) = Any
promote_rule(::Type{Missing}, ::Type{Missing}) = Missing

convert(::Type{Union{T, Missing}}, x) where {T} = convert(T, x)
# To print more appropriate message than "T not defined"
convert(::Type{Missing}, x) = throw(MethodError(convert, (Missing, x)))
convert(::Type{Missing}, ::Missing) = missing

# Comparison operators
==(::Missing, ::Missing) = missing
==(::Missing, ::Any) = missing
==(::Any, ::Missing) = missing
# To fix ambiguity
==(::Missing, ::WeakRef) = missing
==(::WeakRef, ::Missing) = missing
isequal(::Missing, ::Missing) = true
isequal(::Missing, ::Any) = false
isequal(::Any, ::Missing) = false
<(::Missing, ::Missing) = missing
<(::Missing, ::Any) = missing
<(::Any, ::Missing) = missing
isless(::Missing, ::Missing) = false
isless(::Missing, ::Any) = false
isless(::Any, ::Missing) = true

# Unary operators/functions
for f in (:(!), :(+), :(-), :(identity), :(zero), :(one), :(oneunit),
          :(abs), :(abs2), :(sign),
          :(acos), :(acosh), :(asin), :(asinh), :(atan), :(atanh),
          :(sin), :(sinh), :(cos), :(cosh), :(tan), :(tanh),
          :(exp), :(exp2), :(expm1), :(log), :(log10), :(log1p),
          :(log2), :(exponent), :(sqrt), :(gamma), :(lgamma),
          :(iseven), :(ispow2), :(isfinite), :(isinf), :(isodd),
          :(isinteger), :(isreal), :(isnan), :(isempty),
          :(iszero), :(transpose), :(float))
    @eval Math.$(f)(::Missing) = missing
end

for f in (:(Base.zero), :(Base.one), :(Base.oneunit))
    @eval function $(f)(::Type{Union{T, Missing}}) where T
        T === Any && throw(MethodError($f, (Any,)))  # To prevent StackOverflowError
        $f(T)
    end
end

# Binary operators/functions
for f in (:(+), :(-), :(*), :(/), :(^),
          :(div), :(mod), :(fld), :(rem), :(min), :(max))
    @eval begin
        # Scalar with missing
        ($f)(::Missing, ::Missing) = missing
        ($f)(d::Missing, x::Number) = missing
        ($f)(d::Number, x::Missing) = missing
    end
end

# Rounding and related functions
for f in (:(ceil), :(floor), :(round), :(trunc))
    @eval begin
        ($f)(::Missing, digits::Integer=0, base::Integer=0) = missing
        ($f)(::Type{>:Missing}, ::Missing) = missing
        ($f)(::Type{T}, ::Missing) where {T} =
            throw(MissingException("cannot convert a missing value to type $T"))
    end
end

# to avoid ambiguity warnings
(^)(::Missing, ::Integer) = missing

# Bit operators
(&)(::Missing, ::Missing) = missing
(&)(a::Missing, b::Bool) = ifelse(b, missing, false)
(&)(b::Bool, a::Missing) = ifelse(b, missing, false)
(&)(::Missing, ::Integer) = missing
(&)(::Integer, ::Missing) = missing
(|)(::Missing, ::Missing) = missing
(|)(a::Missing, b::Bool) = ifelse(b, true, missing)
(|)(b::Bool, a::Missing) = ifelse(b, true, missing)
(|)(::Missing, ::Integer) = missing
(|)(::Integer, ::Missing) = missing
xor(::Missing, ::Missing) = missing
xor(a::Missing, b::Bool) = missing
xor(b::Bool, a::Missing) = missing
xor(::Missing, ::Integer) = missing
xor(::Integer, ::Missing) = missing

*(d::Missing, x::AbstractString) = missing
*(d::AbstractString, x::Missing) = missing

function float(A::AbstractArray{Union{T, Missing}}) where {T}
    U = typeof(float(zero(T)))
    convert(AbstractArray{Union{U, Missing}}, A)
end
float(A::AbstractArray{Missing}) = A

max_stars_repo_name: JuliaLangEs/julieta

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using CategoricalArrays, XLSX

export shadedmap, myhist, GISturbines_density, aggregate, exportGISturbinedata, landtypes,
        makePixelDataframe, savePixelData, grouppopulationdensity, groupwindspeeds,
        scatter, plot, plot!, Point2f0, RGBA, FRect, plotfix!, heatmap, colorlegend, vbox

const landtypes = ["Evergreen Needleleaf", "Evergreen Broadleaf",
                "Deciduous Needleleaf", "Deciduous Broadleaf", "Mixed Forests",
                "Closed Shrublands", "Open Shrublands", "Woody Savannas", "Savannas", "Grasslands",
                "Permanent Wetlands", "Croplands", "Urban", "Cropland/Natural", "Snow/Ice", "Barren"]

# const landtypes = ["Forest", "Forest", "Forest", "Forest", "Forest", "Shrubland", "Shrubland", "Forest",
#                 "Savanna", "Grassland", "Wetland", "Cropland", "Urban", "Cropland", "Snow/Ice", "Barren"]

function map_protected(; downsample=1, resolutionscale=1, textscale=1, optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, res, protected_codes = options
    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land,
        protected, lonrange, latrange = read_datasets(options)

    n2000_full, ex_n2000 = readraster(in_datafolder("natura2000.tif"), :getextent)
    latrange_n2000, lonrange_n2000 = bbox2ranges(extent2bbox(ex_n2000), 100)
    lons_n2000 = (lonrange[1]:lonrange[end]) .- lonrange_n2000[1] .+ 1
    lats_n2000 = latrange .- latrange_n2000[1] .+ 1
    n2000 = n2000_full[lons_n2000, lats_n2000]

    wdpa = zeros(Bool, size(protected))
    for i in protected_codes
        wdpa[protected .== i] .= true
    end

    reg = (regions.>0) .& (regions .< NOREGION)

    mask = zeros(Int16, size(regions))
    mask[reg] .= 1
    mask[wdpa .& reg] .= 2
    mask[reg .& (n2000 .> 0)] .= 3
    mask[reg .& wdpa .& (n2000 .> 0)] .= 4
    mask[regions.==NOREGION] .= NOREGION

    legendtext = ["unprotected", "WDPA", "Natura2000", "both", "wind turbines"]

    mask = mask[1:downsample:end, 1:downsample:end]
    lonrange = lonrange[1:downsample:end]
    latrange = latrange[1:downsample:end]
    nreg = length(legendtext)

    colors = RGBA.([RGB(0.2,0.3,0.4), RGB(0.7,0.7,0.7), RGB(0,.7,0), RGB(0,0,.7), RGB(.9,.9,0), RGB(1,0,0), RGB(0.4,0.4,0.4)], 0.8)

    println("Mapping colors to regions (avoid same color in adjacent regions)...")
    connected = connectedregions(mask, nreg)
    # colorindices = (nreg > 7) ? greedycolor(connected, 1:7, 1:nreg, randseed=randseed) : collect(1:nreg)
    # colors = colorschemes[Symbol("Set2_$nreg")].colors[colorindices]

    println("\nProjecting coordinates (Mollweide)...")
    res = 0.01
    res2 = res/2
    lons = (-180+res2:res:180-res2)[lonrange]         # longitude values (pixel center)
    lats = (90-res2:-res:-90+res2)[latrange]          # latitude values (pixel center)
    source = Projection("+proj=longlat +datum=WGS84")
    dest = Projection("+proj=moll +lon_0=$(mean(lons)) +ellps=WGS84")
    xs, ys = xygrid(lons, lats)
    Proj4.transform!(source, dest, vec(xs), vec(ys))

    cols = [:lon, :lat, :capac, :year, :onshore, :elec2018]
    df = DataFrame(CSV.File(in_datafolder("turbines_DE.csv")))[:, cols[1:4]]
    tx, ty = df.lon, df.lat
    Proj4.transform!(source, dest, vec(tx), vec(ty))

    println("\nOnshore map...")
    createmap("$(gisregion)_protected", mask, legendtext, lons, lats, colors, source, dest, xs, ys,
        [], [], connected, connected; lines=false, labels=false, resolutionscale=resolutionscale,
        textscale=textscale, legend=true, dots=(tx,ty))
end

function makePixelDataframe(; optionlist...)
    cols = [:lon, :lat, :capac, :year, :onshore, :elec2018]
    df_DK = DataFrame(CSV.File(in_datafolder("turbines_DK.csv")))[:, cols]
    df_SE = DataFrame(CSV.File(in_datafolder("turbines_SE.csv")))[:, cols[1:5]]
    df_DE = DataFrame(CSV.File(in_datafolder("turbines_DE.csv")))[:, cols[1:4]]
    df_USA = DataFrame(CSV.File(in_datafolder("turbines_USA.csv")))[:, cols[1:4]]
    df_SE[:, :elec2018] .= missing    # MWh/year
    df_DE[:, :elec2018] .= missing    # MWh/year
    df_DE[:, :onshore] .= missing    # MWh/year
    df_USA[:, :elec2018] .= missing    # MWh/year
    df_USA[:, :onshore] .= missing    # MWh/year
    df_SE[!,:year] = [ismissing(d) ? 2000 : year(d) for d in df_SE[!,:year]]
    turbines = vcat(df_DK, df_SE, df_DE, df_USA)

    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, res = options
    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land,
        protected, lonrange, latrange = read_datasets(options)

    windatlas = getwindatlas()[lonrange,latrange]
    miuu = zeros(Float32, size(windatlas))

    if gisregion == "SwedenGADM3"
        miuu_raw, ex_miuu = readraster(in_datafolder("miuu_windatlas.tif"), :getextent)
        latrange_miuu, lonrange_miuu = bbox2ranges(extent2bbox(ex_miuu), 100)
        # MIUU has a smaller extent than SwedenGADM3
        lons_miuu = lonrange_miuu .- lonrange[1] .+ 1
        lats_miuu = latrange_miuu .- latrange[1] .+ 1
        miuu[lons_miuu, lats_miuu] = miuu_raw
    end

    n2000, ex_n2000 = readraster(in_datafolder("natura2000.tif"), :getextent)
    latrange_n2000, lonrange_n2000 = bbox2ranges(extent2bbox(ex_n2000), 100)
    lons_n2000 = (lonrange[1]:lonrange[end]) .- lonrange_n2000[1] .+ 1
    lats_n2000 = latrange .- latrange_n2000[1] .+ 1
    n2000 = occursin("USA", gisregion) ? zeros(UInt8, size(windatlas)) : n2000[lons_n2000, lats_n2000]

    nreg = length(regionlist)
    nlon, nlat = size(windatlas)
    turbinecapac = zeros(Int, nlon, nlat)
    nturbines = zeros(Int, nlon, nlat)
    turbineyear = zeros(Int, nlon, nlat)

    for row in eachrow(turbines)
        reg, flon, flat = getregion_and_index(row[:lon], row[:lat], regions, lonrange, latrange)
        if (reg == 0 || reg == NOREGION)
            reg, flon, flat = getregion_and_index(row[:lon], row[:lat], offshoreregions, lonrange, latrange)
            (reg == 0 || reg == NOREGION) && continue
        end
        nturbines[flon,flat] += 1
        turbineyear[flon,flat] = max(turbineyear[flon,flat], row[:year])
        turbinecapac[flon,flat] += row[:capac]
    end

    region = max.(regions, offshoreregions)
    lons = (-180+res/2:res:180-res/2)[lonrange]   # longitude values (pixel center)
    lats = (90-res/2:-res:-90+res/2)[latrange]    # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)
    lon = repeat(lons, 1, nlat)
    lat = repeat(lats', nlon)
    area = repeat(cellarea', nlon)

    r = (region .> 0) .& (region .!= NOREGION)

    return DataFrame(lat=lat[r], lon=lon[r], munic=region[r],
        area=round.(area[r], sigdigits=4), landtype=land[r], protected=protected[r], natura2000=Int.(n2000[r]),
        popdens=round.(popdens[r], sigdigits=4), windspeed=round.(windatlas[r], sigdigits=4), miuu=round.(miuu[r], sigdigits=4),
        nturbines=nturbines[r], turbinecapac=turbinecapac[r], turbineyear=turbineyear[r])
end

function savePixelData()
    gisregions = ["SwedenGADM3", "Denmark83", "GermanyGADM3", "USA_Texas", "USA_Iowa", "USA_Oklahoma",
        "USA_Kansas", "USA_Illinois", "USA_South Dakota", "USA_North Dakota", "USA_California",
        "USA_Minnesota", "USA_Colorado", "USA_Oregon"]
    countries = ["Sweden", "Denmark", "Germany", "Texas", "Iowa", "Oklahoma", "Kansas", "Illinois",
        "South Dakota", "North Dakota", "California", "Minnesota", "Colorado", "Oregon"]
    println("Building municipality list...")
    munic, nmunic = Dict(), Dict()
    for (gisregion, country) in zip(gisregions, countries)
        _, _, regionlist, _, _ = loadregions(gisregion)
        munic[country] = regionlist
        nmunic[gisregion] = length(regionlist)
    end
    dfmunic = vcat((DataFrame(country=country, municipality=munic[country])
                for country in countries)...)
    insertcols!(dfmunic, 1, :index => 1:nrow(dfmunic))
    CSV.write("D:/GISdata/municipalities.csv", dfmunic)
    println("Building main DataFrame...")
    println(countries[1])
    dfall = makePixelDataframe(gisregion=gisregions[1])
    countrynum = 1
    insertcols!(dfall, 3, :country => countrynum)
    len = nmunic[gisregions[1]]
    for gisregion in gisregions[2:end]
        println(gisregion)
        df = makePixelDataframe(gisregion=gisregion)
        df[!,:munic] .+= len
        countrynum += 1
        insertcols!(df, 3, :country => countrynum)
        dfall = vcat(dfall, df)
        len += nmunic[gisregion]
    end
    println("Saving DataFrame...")
    CSV.write("D:/GISdata/windpixeldata.csv", dfall)
end

extent2bbox(ex) = [ex[2] ex[1]; ex[4] ex[3]]

function protected_vs_natura2000(; optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land,
        protected, lonrange, latrange = read_datasets(options)

    n2000, ex_n2000 = readraster(in_datafolder("natura2000.tif"), :getextent)
    latrange_n2000, lonrange_n2000 = bbox2ranges(extent2bbox(ex_n2000), 100)
    lons_n2000 = (lonrange[1]:lonrange[end]) .- lonrange_n2000[1] .+ 1
    lats_n2000 = latrange .- latrange_n2000[1] .+ 1
    n2000 = n2000[lons_n2000, lats_n2000]

    onshore = (regions .> 0) .& (regions .!= NOREGION) .& (land .> 0)
    offshore = (offshoreregions .> 0) .& (offshoreregions .!= NOREGION) .& (land .== 0)

    count_onshore = zeros(Int, 11, 4)
    count_offshore = zeros(Int, 11, 4)
    for i = 1:length(protected)
        if onshore[i]
            count_onshore[protected[i]+1, n2000[i]+1] += 1
        elseif offshore[i]
            count_offshore[protected[i]+1, n2000[i]+1] += 1
        end
    end
    display(count_onshore)
    display(count_offshore)

    return protected, n2000
end

function miuu_vs_windatlas()
    options = windoptions()
    options[:gisregion] = "SwedenGADM3"
    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land,
        protected, lonrange, latrange = read_datasets(options)

    windatlas = getwindatlas()[lonrange,latrange]
    onshore = (regions .> 0) .& (regions .!= NOREGION) .& (land .> 0)
    offshore = (offshoreregions .> 0) .& (offshoreregions .!= NOREGION) .& (land .== 0)

    miuu, ex_miuu = readraster(in_datafolder("miuu_windatlas.tif"), :getextent)
    latrange_miuu, lonrange_miuu = bbox2ranges(extent2bbox(ex_miuu), 100)
    # MIUU has a smaller extent than SwedenGADM3
    lons_miuu = lonrange_miuu .- lonrange[1] .+ 1
    lats_miuu = latrange_miuu .- latrange[1] .+ 1

    windatlas = windatlas[lons_miuu, lats_miuu]
    onshore = onshore[lons_miuu, lats_miuu] .& (miuu .> 0)
    offshore = offshore[lons_miuu, lats_miuu] .& (miuu .> 0)

    return miuu, windatlas, onshore, offshore
end

function plotfix!(scene)
    xlabel!(scene, "Global Wind Atlas [m/s]")
    ylabel!(scene, "MIUU [m/s]")
    plot!(scene, [0,14],[0,14], color=:red)
end

#=
mm,ww,on,off = GE.miuu_vs_windatlas();
p = heatmap(reverse(ww.*(mm.>0),dims=2), resolution=(950,950), colorrange=(3,10), highclip=:red, lowclip=:black)
p = heatmap(reverse(mm,dims=2), resolution=(950,950), colorrange=(3,10), highclip=:red, lowclip=:black)
        # cl = colorlegend(p[end], resolution=(950,950)); vbox(p,vl)
opt = (color=RGBA(0,0,0,0.2), markersize=0.01, limits=FRect(0,0,14,14), textsize=10)
p = scatter(Point2f0.(ww[on], mm[on]+randn(sum(on))*.05); opt..., color=RGBA(0,0,0,1)); plotfix!(p)
p = scatter(Point2f0.(ww[off], mm[off]+randn(sum(off))*.05); opt..., color=RGBA(0,0,0,1)); plotfix!(p)

pp,nn = GE.protected_vs_natura2000(gisregion="SwedenGADM3");
p = heatmap(reverse(nn,dims=2), resolution=(950,950), highclip=:red, lowclip=:black)
p = heatmap(reverse(pp,dims=2), resolution=(950,950), highclip=:red, lowclip=:black)

=#

function analyze_protected(; firstyear=1978, lastyear=2021)
    df = CSV.File("D:/GISdata/windpixeldata.csv") |> DataFrame

    df.inyears = (df.turbineyear.>=firstyear) .& (df.turbineyear.<=lastyear)
    countries = ["Sweden", "Denmark", "Germany", "Texas", "Iowa", "Oklahoma", "Kansas", "Illinois",
        "South Dakota", "North Dakota", "California", "Minnesota", "Colorado", "Oregon"]
    df.countryname = countries[df.country]
    df.onshore = df.landtype .> 0
    protected_names = [
       "(Ia) Strict Nature Reserve"
       "(Ib) Wilderness Area"
       "(II) National Park"
       "(III) Natural Monument"
       "(IV) Habitat/Species Management"
       "(V) Protected Landscape/Seascape"
       "(VI) Managed Resource Protected Area"
       "Not Reported"
       "Not Applicable"
       "Not Assigned"
    ]
    protected_type = Dict(i => n for (i,n) in enumerate(protected_names))
    protected_type[0] = "UNPROTECTED"

    gdf = groupby(df, [:countryname, :onshore, :protected])
    cdf = combine(gdf,
            :area => (a -> sum(a)/1000) => :area,
            [:turbinecapac, :inyears] => ((c,y) -> sum(c.*y)/1000) => :capac
        )
    gdf_tot = groupby(cdf, [:countryname, :onshore])
    cdf_tot = combine(gdf_tot, :area => sum, :capac => sum)

    cdf.turbinedensity = cdf.capac ./ cdf.area
    cdf.protected_type = getindex.(Ref(protected_type), cdf.protected)
    cdf_tot.turbinedensity = cdf_tot.capac_sum ./ cdf_tot.area_sum

    sort!(cdf, [order(:countryname, by = x -> findfirst(countries.==x)), :protected, order(:onshore, rev=true)])
    sort!(cdf_tot, [order(:countryname, by = x -> findfirst(countries.==x)), order(:onshore, rev=true)])
    # CSV.write("protectedinfo.csv", cdf)
    # CSV.write("protectedinfo_tot.csv", cdf_tot)

    filename = "protected WDPA.xlsx"
    XLSX.openxlsx(filename, mode=isfile(filename) ? "rw" : "w") do xf
        sheetname = "$firstyear-$lastyear"
        sheet = sheetname in XLSX.sheetnames(xf) ? xf[sheetname] : XLSX.addsheet!(xf, sheetname)
        XLSX.writetable!(sheet, [cdf[:,i] for i=1:size(cdf,2)], names(cdf), anchor_cell=XLSX.CellRef("A4"))
        XLSX.writetable!(sheet, [cdf_tot[:,i] for i=1:size(cdf_tot,2)], names(cdf_tot), anchor_cell=XLSX.CellRef("L4"))
    end

    return cdf, cdf_tot
end

function analyze_natura2000(; firstyear=1978, lastyear=2021)
    df = CSV.File("D:/GISdata/windpixeldata.csv") |> DataFrame

    df.inyears = (df.turbineyear.>=firstyear) .& (df.turbineyear.<=lastyear)
    countries = ["Sweden", "Denmark", "Germany", "Texas", "Iowa", "Oklahoma", "Kansas", "Illinois",
        "South Dakota", "North Dakota", "California", "Minnesota", "Colorado", "Oregon"]
    df.countryname = countries[df.country]
    df.onshore = df.landtype .> 0
    protected_names = [
       "A: SPAs"
       "B: SCIs and SACs"
       "C: both categories A + B"
    ]
    protected_type = Dict(i => n for (i,n) in enumerate(protected_names))
    protected_type[0] = "UNPROTECTED"

    gdf = groupby(df, [:countryname, :onshore, :natura2000])
    cdf = combine(gdf,
            :area => (a -> sum(a)/1000) => :area,
            [:turbinecapac, :inyears] => ((c,y) -> sum(c.*y)/1000) => :capac
        )
    gdf_tot = groupby(cdf, [:countryname, :onshore])
    cdf_tot = combine(gdf_tot, :area => sum, :capac => sum)

    cdf.turbinedensity = cdf.capac ./ cdf.area
    cdf.protected_type = getindex.(Ref(protected_type), cdf.natura2000)
    cdf_tot.turbinedensity = cdf_tot.capac_sum ./ cdf_tot.area_sum

    sort!(cdf, [order(:countryname, by = x -> findfirst(countries.==x)), :natura2000, order(:onshore, rev=true)])
    sort!(cdf_tot, [order(:countryname, by = x -> findfirst(countries.==x)), order(:onshore, rev=true)])

    filename = "protected Natura2000.xlsx"
    XLSX.openxlsx(filename, mode=isfile(filename) ? "rw" : "w") do xf
        sheetname = "$firstyear-$lastyear"
        sheet = sheetname in XLSX.sheetnames(xf) ? xf[sheetname] : XLSX.addsheet!(xf, sheetname)
        XLSX.writetable!(sheet, [cdf[:,i] for i=1:size(cdf,2)], names(cdf), anchor_cell=XLSX.CellRef("A4"))
        XLSX.writetable!(sheet, [cdf_tot[:,i] for i=1:size(cdf_tot,2)], names(cdf_tot), anchor_cell=XLSX.CellRef("L4"))
    end
    # CSV.write("protectedinfo.csv", cdf)
    # CSV.write("protectedinfo_tot.csv", cdf_tot)

    return cdf, cdf_tot
end

function analyze_landtype(; firstyear=1978, lastyear=2021, minspeed=0)
    df = CSV.File("D:/GISdata/windpixeldata.csv") |> DataFrame

    df.inyears = (df.turbineyear.>=firstyear) .& (df.turbineyear.<=lastyear) .& (df.windspeed .>= minspeed)
    countries = ["Sweden", "Denmark", "Germany", "Texas", "Iowa", "Oklahoma", "Kansas", "Illinois",
        "South Dakota", "North Dakota", "California", "Minnesota", "Colorado", "Oregon"]
    df.countryname = countries[df.country]
    landtype_names = [
        "Water"
        "Evergreen Needleleaf Forests"
        "Evergreen Broadleaf Forests"
        "Deciduous Needleleaf Forests"
        "Deciduous Broadleaf Forests"
        "Mixed Forests"
        "Closed Shrublands"
        "Open Shrublands"
        "Woody Savannas"
        "Savannas"
        "Grasslands"
        "Permanent Wetlands"
        "Croplands"
        "Urban"
        "Cropland/Natural"
        "Snow/Ice"
        "Barren"
    ]
    landtype = Dict(i-1 => n for (i,n) in enumerate(landtype_names))

    gdf = groupby(df, [:countryname, :landtype])
    cdf = combine(gdf,
            :area => (a -> sum(a)/1000) => :area,
            [:turbinecapac, :inyears] => ((c,y) -> sum(c.*y)/1000) => :capac
        )
    gdf_tot = groupby(cdf, :countryname)
    cdf_tot = combine(gdf_tot, :area => sum, :capac => sum)
    cdf.turbinedensity = cdf.capac ./ cdf.area
    cdf.landtype_name = getindex.(Ref(landtype), cdf.landtype)
    cdf_tot.turbinedensity = cdf_tot.capac_sum ./ cdf_tot.area_sum

    sort!(cdf, [order(:countryname, by = x -> findfirst(countries.==x)), :landtype])
    sort!(cdf_tot, [order(:countryname, by = x -> findfirst(countries.==x))])

    filename = "landtypes minspeed $minspeed.xlsx"
    XLSX.openxlsx(filename, mode=isfile(filename) ? "rw" : "w") do xf
        sheetname = "$firstyear-$lastyear"
        sheet = sheetname in XLSX.sheetnames(xf) ? xf[sheetname] : XLSX.addsheet!(xf, sheetname)
        XLSX.writetable!(sheet, [cdf[:,i] for i=1:size(cdf,2)], names(cdf), anchor_cell=XLSX.CellRef("A4"))
        XLSX.writetable!(sheet, [cdf_tot[:,i] for i=1:size(cdf_tot,2)], names(cdf_tot), anchor_cell=XLSX.CellRef("L4"))
    end

    return cdf, cdf_tot
end

function grouppopulationdensity(country; firstyear=1978, lastyear=2021)
    df = CSV.File("D:/GISdata/windpixeldata.csv") |> DataFrame
    df = df[df.country .== country, :]
    df.inyears = (df.turbineyear.>=firstyear) .& (df.turbineyear.<=lastyear)
    df.nturbines = df.nturbines .* df.inyears
    df.capac = df.turbinecapac/1000 .* df.inyears
    pops = [[i*10.0^j for j = -1:3 for i in [1,2,5]]; 100_000]
    popmin = [0; pops[1:end-1]]
    df.poprange = CategoricalArrays.cut(df.popdens, pops, extend=true)
    gdf = groupby(df, :poprange)
    cdf = combine(gdf, :nturbines .=> sum, :capac .=> sum, nrow => :pixels, :area .=> sum,
        renamecols = false)
    allranges = CategoricalArrays.cut(popmin, pops, extend=true)
    df_out = similar(cdf, length(allranges))
    df_out.poprange = allranges
    df_out[:,2:end] .= 0
    insertcols!(df_out, 2, :popmin => popmin, :popmax => pops)
    indexes = [findfirst(allranges .== r) for r in cdf.poprange]
    df_out[indexes, 4:end] = cdf[:, 2:end]
    countries = ["Sweden", "Denmark", "Germany", "Texas", "Iowa", "Oklahoma", "Kansas", "Illinois",
        "South Dakota", "North Dakota", "California", "Minnesota", "Colorado", "Oregon"]

    filename = "wind_popdens 1978-2021.xlsx"
    XLSX.openxlsx(filename, mode=isfile(filename) ? "rw" : "w") do xf
        sheetname = "$(countries[country]) $firstyear-$lastyear"
        sheet = sheetname in XLSX.sheetnames(xf) ? xf[sheetname] : XLSX.addsheet!(xf, sheetname)
        XLSX.writetable!(sheet, [df_out[:,i] for i=2:size(df_out,2)], names(df_out)[2:end], anchor_cell=XLSX.CellRef("A1"))
    end
    # CSV.write(filename, df_out[:, 2:end])
end

# just for the histogram data Fredrik wanted
function groupwindspeeds(country; firstyear=1978, lastyear=2021, usemiuu=false)
    df = CSV.File("D:/GISdata/windpixeldata.csv") |> DataFrame
    df = df[df.country .== country, :]
    df.inyears = (df.turbineyear.>=firstyear) .& (df.turbineyear.<=lastyear)
    df.is_onshore = (df.landtype .> 0)
    df.is_offshore = (df.landtype .== 0)
    df.nturbines = df.nturbines .* df.inyears
    df.nturbines_on = df.nturbines .* df.is_onshore
    df.nturbines_off = df.nturbines .* df.is_offshore
    df.capac = df.turbinecapac/1000 .* df.inyears
    df.capac_on = df.capac .* df.is_onshore
    df.capac_off = df.capac .* df.is_offshore
    df.area_on = df.area .* df.is_onshore
    df.area_off = df.area .* df.is_offshore
    df.speed = usemiuu ? df.miuu : df.windspeed
    df.windspeed_range = CategoricalArrays.cut(df.speed, 0:.25:21, extend=true)
    gdf = groupby(df, :windspeed_range)
    cdf = combine(gdf, [:nturbines, :nturbines_on, :nturbines_off] .=> sum,
        [:capac, :capac_on, :capac_off] .=> sum,
        nrow => :pixels, [:is_onshore, :is_offshore] .=> sum .=> [:pixels_on, :pixels_off],
        [:area, :area_on, :area_off] .=> sum,
        renamecols = false)
    allranges = CategoricalArrays.cut(0:.25:20.75, 0:.25:21, extend=true)
    df_out = similar(cdf, length(allranges))
    df_out.windspeed_range = allranges
    df_out[:,2:end] .= 0
    indexes = [findfirst(allranges .== r) for r in cdf.windspeed_range]
    df_out[indexes, :] = cdf
    insertcols!(df_out, 2, :windspeed_low => 0:.25:20.75, :windspeed_high => .25:.25:21)

    countries = ["Sweden", "Denmark", "Germany", "Texas", "Iowa", "Oklahoma", "Kansas", "Illinois",
        "South Dakota", "North Dakota", "California", "Minnesota", "Colorado", "Oregon"]

    filename = "winddata 1978-2021.xlsx"
    XLSX.openxlsx(filename, mode=isfile(filename) ? "rw" : "w") do xf
        sheetname = "$(countries[country]) $(usemiuu ? "MIUU " : "")$firstyear-$lastyear"
        sheet = sheetname in XLSX.sheetnames(xf) ? xf[sheetname] : XLSX.addsheet!(xf, sheetname)
        XLSX.writetable!(sheet, [df_out[:,i] for i=2:size(df_out,2)], names(df_out)[2:end], anchor_cell=XLSX.CellRef("A1"))
    end
    # CSV.write(filename, df_out)
end

#=
# run in REPL with using Plots & plotly()
function plothists(region, windspeed, onshoreturbine, onshore, offshore, years)
    region = replace(region, "GADM" => "")
    region = replace(region, r"\d" => "")
    title = "$region $(years[1]) - $(years[2])"
    histopt = (legend=:none, size=(1200,550), tickfont=12, bins=2.8:.1:11.2,
                    xticks=3:1:11, normalize=true)
    barhist(windspeed[onshoreturbine], title=title, c=:lightgray; histopt...)
    stephist!(windspeed[onshore .& (windspeed.>0)], lw=1.5, c=:blue; histopt...)
    stephist!(windspeed[offshore .& (windspeed.>0)], lw=1.5, c=:red; histopt...)
end

    # groupedhist([windspeed[onshoreturbine]; windspeed[offshoreturbine]];
    #     group=[ones(sum(onshoreturbine)); 2*ones(sum(offshoreturbine))],
    #     bar_position = :stack, c=[:lightgray :red], histopt...)
=#

function exportGISturbinedata(; mincapac=0, minclass=0, firstyear=1978, plotmasks=false, optionlist...)
    cols = [:lon, :lat, :capac, :onshore, :elec2018]
    df_DK = DataFrame(CSV.File(in_datafolder("turbines_DK.csv")))[:, cols]
    df_SE = DataFrame(CSV.File(in_datafolder("turbines_SE.csv")))[:, cols[1:4]]
    df_DE = DataFrame(CSV.File(in_datafolder("turbines_DE.csv")))[:, cols[1:3]]
    df_USA = DataFrame(CSV.File(in_datafolder("turbines_USA.csv")))[:, cols[1:3]]
    df_SE[:, :elec2018] .= missing    # MWh/year
    df_DE[:, :elec2018] .= missing    # MWh/year
    df_DE[:, :onshore] .= missing    # MWh/year
    df_USA[:, :elec2018] .= missing    # MWh/year
    df_USA[:, :onshore] .= missing    # MWh/year
    turbines = vcat(df_DK, df_SE, df_DE, df_USA)
    turbines = turbines[turbines[!,:capac].>=mincapac, :]

    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, downsample_masks, onshore_density, area_onshore, res = options
    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land,
        protected_mat, lonrange, latrange =
                read_datasets(options)
    mask_onshoreA, mask_onshoreB, mask_offshore,
        gridA, lowpopdens, goodland, protected_area =
            return_wind_masks(options, regions, offshoreregions, gridaccess, popdens,
                            topo, land, protected_mat, lonrange, latrange,
                            plotmasks=plotmasks, downsample=downsample_masks)
    mask_onshore = mask_onshoreA .| mask_onshoreB
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)

    windatlas = getwindatlas()[lonrange,latrange]
    onshoreclass, offshoreclass = makewindclasses(options, windatlas)

    nreg = length(regionlist)
    capac, offcapac, elec, pop = zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg)
    onshore, offshore = zeros(nreg), zeros(nreg)
    n_onshore, n_offshore = zeros(Int, nreg), zeros(Int, nreg)
    landcount, commonland, freq = zeros(Int, nreg, 16), fill("", nreg, 3), zeros(nreg, 3)
    maskedturbines, nogrid, highpop, badland, protected =
        zeros(Int, nreg), zeros(Int, nreg), zeros(Int, nreg), zeros(Int, nreg), zeros(Int, nreg)
    capac_ok, area, area_ok, maxcapac, maxcapac_ok, class, windspeed =
        zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg)

    for row in eachrow(turbines)
        reg, flon, flat = getregion_and_index(row[:lon], row[:lat], regions, lonrange, latrange)
        if (reg == 0 || reg == NOREGION)
            reg, flon, flat = getregion_and_index(row[:lon], row[:lat], offshoreregions, lonrange, latrange)
            (reg == 0 || reg == NOREGION) && continue
        end
        if ismissing(row[:onshore]) || row[:onshore]
            onshoreclass[flon, flat] < minclass && continue
            capac[reg] += row[:capac] / 1000
            capac_ok[reg] += row[:capac] / 1000 * mask_onshore[flon, flat]
            class[reg] += onshoreclass[flon, flat]
            windspeed[reg] += windatlas[flon, flat]
            onshore[reg] += mask_onshore[flon, flat]
            n_onshore[reg] += 1
            if land[flon, flat] != 0
                landcount[reg, land[flon, flat]] += 1
            end
            pop[reg] += capac[reg]*popdens[flon, flat]  # capacity-weighted
            if !mask_onshore[flon, flat]
                maskedturbines[reg] += 1
                nogrid[reg] += !gridA[flon, flat]
                highpop[reg] += !lowpopdens[flon, flat]
                badland[reg] += !goodland[flon, flat]
                protected[reg] += protected_area[flon, flat]
            end
        else
            offshoreclass[flon, flat] < minclass && continue
            offcapac[reg] += row[:capac] / 1000
            offshore[reg] += mask_offshore[flon, flat]
            n_offshore[reg] += 1
        end
        elec[reg] += coalesce(row[:elec2018] / 1000, 0.0)
    end
    freg = (regions.>0) .& (regions.!=NOREGION)
    all_onshore = sum(mask_onshore[freg]) / sum(freg)
    fregoff = (offshoreregions.>0) .& (offshoreregions.!=NOREGION)
    all_offshore = sum(mask_offshore[fregoff]) / sum(fregoff)
    println("\nTotal onshore OK: ", percentstring(sum(onshore)/sum(n_onshore)),
            " ($(percentstring(all_onshore)))")
    println("Total offshore OK: ", percentstring(sum(offshore)/sum(n_offshore)),
            " ($(percentstring(all_offshore)))")
    regionalpopdens = sum(popdens[freg]) / sum(freg)
    println("Masked turbines: ", percentstring(sum(maskedturbines)/sum(n_onshore)))
    println("\tno grid: ", percentstring(sum(nogrid)/sum(maskedturbines)))
    println("\thigh pop: ", percentstring(sum(highpop)/sum(maskedturbines)))
    println("\tbad land: ", percentstring(sum(badland)/sum(maskedturbines)))
    println("\tprotected: ", percentstring(sum(protected)/sum(maskedturbines)))
    println("Total population density (persons/km2): ", round(sum(pop)/sum(n_onshore)/sum(capac), digits=1),
            " ($(round(regionalpopdens, digits=1)))")
    onshore = round.(onshore ./ n_onshore * 100, digits=1)
    class = class ./ n_onshore
    windspeed = windspeed ./ n_onshore
    offshore = round.(offshore ./ n_offshore * 100, digits=1)
    pop = round.(pop ./ n_onshore ./ capac, digits=1)
    masked = round.(maskedturbines ./ n_onshore * 100, digits=1)
    nogrid = round.(nogrid ./ maskedturbines * 100, digits=1)
    highpop = round.(highpop ./ maskedturbines * 100, digits=1)
    badland = round.(badland ./ maskedturbines * 100, digits=1)
    protected = round.(protected ./ maskedturbines * 100, digits=1)
    println("Calculating results per municipality/county...")
    updateprogress = Progress(nreg, 1)
    for reg = 1:nreg
        area[reg] += sum((regions .== reg) .* cellarea')
        maxcapac[reg] = area[reg] * onshore_density * area_onshore
        area_ok[reg] += sum(((regions .== reg) .& mask_onshore) .* cellarea')
        maxcapac_ok[reg] = area_ok[reg] * onshore_density * area_onshore
        ss = sortslices([landcount[reg, :] collect(1:16)], dims=1, rev=true)
        commonland[reg,:] = getindex.(Ref(landtypes), ss[1:3, 2])
        freq[reg,:] = round.(ss[1:3, 1]/n_onshore[reg] * 100, digits=1)
        next!(updateprogress)
    end
    println("Turbine density of unmasked land (MW/km2): ",
            round(sum(capac_ok)/sum(area_ok), digits=3),
            " ($(round(sum(capac)/sum(area), digits=3)))")
    println("Exploited share of potential capacity (%): ",
            round(sum(capac_ok)/sum(maxcapac_ok)*100, digits=1),
            " ($(round(sum(capac)/sum(maxcapac)*100, digits=1)))")
    ss = sortslices([sum(landcount, dims=1)' collect(1:16)], dims=1, rev=true)
    for i = 1:3
        println("Landtype $i: $(landtypes[ss[i,2]]) ($(percentstring(ss[i,1]/sum(n_onshore))))")
    end
    df = DataFrame(region=regionlist, capac=round.(capac, digits=2),
            offcapac=round.(offcapac, digits=2), elec2018=round.(elec, digits=2),
            onshore=replace(onshore, NaN=>missing), offshore=replace(offshore, NaN=>missing),
            masked=masked, nogrid=nogrid, highpop=highpop, badland=badland, protected=protected,
            dens_tot=round.(capac./area, digits=3), dens_ok=round.(capac_ok./area_ok, digits=3),
            exploit_tot=round.(area_onshore*capac./maxcapac*100, digits=1),
            exploit_ok=round.(area_onshore*capac_ok./maxcapac_ok*100, digits=1),
            popdens=pop, land1=commonland[:,1], freq1=freq[:,1], land2=commonland[:,2], freq2=freq[:,2],
            land3=commonland[:,3], freq3=freq[:,3], area=round.(Int, area), class=class,
            windspeed=windspeed)
    outfile = "regionalwindGIS_$(gisregion)_minturbine=$(mincapac)_maxpopdens=$(Int(options.persons_per_km2)).csv"
    CSV.write(in_datafolder("output", outfile), df)
    df
end

percentstring(x) = "$(round(100*x, digits=1))%"

function turbinecharts(gisregion)
    td, pd = GISturbines_density(gisregion=gisregion);
    nz = (pd.>1) .| (td.>0);
    scene = plot(histogram(nbins = 50), td[td.>0])
    xlabel!("turbine density (MW/km2)")
    # scene = plot(histogram(nbins = 50), log10.(pd[nz]))
    # xlabel!("pop. density (log10 persons/km2)")

    # scene = StatsMakie.scatter(td[nz], log10.(pd[nz]), markersize=0.2)
    # xlabel!("turbine density (MW/km2)")
    # ylabel!("pop. density (log10 persons/km2)")

    scene

    # shadedmap("Denmark5", log.(1 .+tdm), downsample=2, colorscheme=:magma)
    # shadedmap("Denmark5", log10.(1 .+pd), downsample=2, colorscheme=:magma)
end

function myhist(data; normalize=false, args...)
    h = fit(Histogram, data; args...)
    weights = normalize ? h.weights/sum(h.weights) : h.weights
    barplot(h.edges[1][1:end-1], weights)
end

function GISturbines_density(; agg=1, optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    @unpack res, gisregion, scenarioyear = options

    println("\nReading regions and population dataset...")
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)

    pop = JLD.load(in_datafolder("population_$scenarioyear.jld"), "population")[lonrange,latrange]
    popdens = aggregate_array(pop, agg) ./ aggregate_array(cellarea, agg)'      # persons/km2

    df_DK = DataFrame(CSV.File(in_datafolder("turbines_DK.csv")))[:, [:lon, :lat, :capac]]
    df_USA = DataFrame(CSV.File(in_datafolder("turbines_USA.csv")))[:, [:lon, :lat, :capac]]
    df_SE = DataFrame(CSV.File(in_datafolder("turbines_SE.csv")))[:, [:lon, :lat, :capac]]
    df_UK = DataFrame(CSV.File(in_datafolder("turbines_UK.csv")))[:, [:lon, :lat, :capac]]
    turbines = vcat(df_DK, df_USA, df_SE, df_UK)

    turbinecapacity = aggregate_turbine_capacity(gisregion, turbines, regions, regionlist, lonrange, latrange)
    turbinedensity = aggregate_array(turbinecapacity, agg) ./ aggregate_array(cellarea, agg)'    # MW/km2

    # println("\nSaving...")

    # mkpath(in_datafolder("output"))
    # matopen(in_datafolder("output", "GISdata_turbines_$gisregion.mat"), "w") do file
    #     write(file, "turbinecapacity", turbinecapacity)
    # end
    return turbinedensity, popdens
end

# function aggregate_array(a::AbstractArray, n::Int)
#     sz = size(a) .÷ n
#     out = similar(a, sz)
#     R = CartesianIndices(A)
#     Ifirst, Ilast = first(R), last(R)
#     I1 = oneunit(Ifirst)
#     for I in R
#         n, s = 0, zero(eltype(out))
#         for J in max(Ifirst, I-I1):min(Ilast, I+I1)
#             s += A[J]
#             n += 1
#         end
#         out[I] = s/n
#     end
#     out
# end

function aggregate_array(a::AbstractArray, n::Int)
    n == 1 && return a
    sz = size(a)
    dims = length(sz)
    out = similar(a, size(a) .÷ n)
    R = CartesianIndices(a)[ntuple(i -> 1:n:sz[i]-n+1, dims)...]   # index to every nth element in all dims
    Rsub = CartesianIndices(ntuple(i -> 0:n-1, dims)) #
    for (i, I) in enumerate(R)
        s = zero(eltype(out))
        for J in I .+ Rsub
            s += a[J]
        end
        out[i] = s
    end
    out
end

function aggregate_turbine_capacity(gisregion, turbines, regions, regionlist, lonrange, latrange)
    println("Calculate turbine capacity per pixel...")

    numreg = length(regionlist)
    turbinecapacity = zeros(size(regions))       # MW

    for i = 1:size(turbines,1)
        reg, flon, flat = getregion_and_index(turbines[i,:lon], turbines[i,:lat], regions, lonrange, latrange)
        (reg == 0 || reg == NOREGION) && continue

        turbinecapacity[flon, flat] += turbines[i,:capac]/1000
    end

    return turbinecapacity
end

function getregion_and_index(lon, lat, regions, lonrange=1:36000, latrange=1:18000)
    res = 0.01
    res2 = res/2
    lons = (-180+res2:res:180-res2)[lonrange]         # longitude values (pixel center)
    lats = (90-res2:-res:-90+res2)[latrange]          # latitude values (pixel center)
    flon = findfirst(round(lon-res2, digits=5) .< lons .<= round(lon+res2, digits=5))
    flat = findfirst(round(lat-res2, digits=5) .< lats .<= round(lat+res2, digits=5))
    (flon == nothing || flat == nothing) && return 0, flon, flat
    # println(lons[flon],", ",lats[flat])
    # println(flon, " ", flat)
    if regions[flon, flat] > 0
        return regions[flon, flat], flon, flat
    else
        return NOREGION, flon, flat
    end
end

function shadedmap(gisregion, plotdata, regionlist, lons, lats, colors, source, dest, xs, ys;
                    resolutionscale=1)
    nreg = length(regionlist)
    scale = maximum(size(plotdata))/6500

    # data[regions.==NOREGION] .= nreg + 1

    xmin, xmax = extrema(xs)
    ymin, ymax = extrema(ys)
    aspect_ratio = (ymax - ymin) / (xmax - xmin)
    pngwidth = round(Int, resolutionscale*1.02*size(plotdata,1))     # allow for margins (+1% on both sides)
    pngsize = pngwidth, round(Int, pngwidth * aspect_ratio)     # use aspect ratio after projection transformation

    println("...constructing map...")
    scene = surface(xs, ys; color=float.(plotdata), colormap=colors,
                            shading=false, show_axis=false, scale_plot=false, interpolate=false)
    ga = geoaxis!(scene, lons[1], lons[end], lats[end], lats[1]; crs=(src=source, dest=dest,))[end]
    ga.x.tick.color = RGBA(colorant"black", 0.4)
    ga.y.tick.color = RGBA(colorant"black", 0.4)
    ga.x.tick.width = scale
    ga.y.tick.width = scale

    println("...saving...")
    mkpath(in_datafolder("output"))
    filename = in_datafolder("output", "$(gisregion)_shadedmap.png")
    isfile(filename) && rm(filename)
    Makie.save(filename, scene, resolution=pngsize)
end

function shadedmap(gisregion, plotdata; resolutionscale=1, downsample=1, colorscheme=:viridis)
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    plotdata = plotdata[1:downsample:end, 1:downsample:end]
    lonrange = lonrange[1:downsample:end]
    latrange = latrange[1:downsample:end]
    nreg = length(regionlist)

    colors = colorschemes[colorscheme].colors  #[RGB(0.2,0.3,0.4); colorschemes[:viridis].colors; RGB(0.4,0.4,0.4)]

    println("\nProjecting coordinates (Mollweide)...")
    res = 0.01
    res2 = res/2
    lons = (-180+res2:res:180-res2)[lonrange]         # longitude values (pixel center)
    lats = (90-res2:-res:-90+res2)[latrange]          # latitude values (pixel center)
    source = Projection("+proj=longlat +datum=WGS84")
    dest = Projection("+proj=moll +lon_0=$(mean(lons))")
    xs, ys = xygrid(lons, lats)
    Proj4.transform!(source, dest, vec(xs), vec(ys))

    println("\nOnshore map...")
    shadedmap(gisregion, plotdata, regionlist, lons, lats, colors, source, dest, xs, ys,
        resolutionscale=resolutionscale)
    # exit()
    return nothing
end

function return_wind_masks(options, regions, offshoreregions, gridaccess, popdens, topo, land, protected, lonrange, latrange; plotmasks=false, downsample=1)
    @unpack res, gisregion, exclude_landtypes, protected_codes, distance_elec_access, persons_per_km2,
                min_shore_distance, max_depth, classB_threshold = options

    println("Creating masks...")

    goodland = (regions .> 0)
    for i in exclude_landtypes
        goodland[land .== i] .= false
    end
    protected_area = zeros(Bool, size(protected))
    for i in protected_codes
        protected_area[protected .== i] .= true
    end

    # Pixels with electricity access for onshore wind A
    gridA = (gridaccess .> 0)

    # Pixels with electricity access for onshore wind B and offshore wind
    km_per_degree = π*2*6371/360
    disk = diskfilterkernel(distance_elec_access/km_per_degree/res)
    gridB = (imfilter(gridaccess, disk) .> max(1e-9, classB_threshold)) # avoid artifacts if classB_threshold == 0

    # println("MAKE SURE MASKS DON'T OVERLAP! (regions & offshoreregions, mask_*)")

    # all mask conditions
    mask_onshoreA = gridA .& (popdens .< persons_per_km2) .& goodland .& .!protected_area
    mask_onshoreB = (gridB .& .!gridA) .& (popdens .< persons_per_km2) .& goodland .& .!protected_area

    # shoreline mask for offshore wind
    disk = diskfilterkernel(min_shore_distance/km_per_degree/res)
    shore = (imfilter(regions .> 0, disk) .> 1e-6)

    # all mask conditions
    mask_offshore = gridB .& .!shore .& (topo .> -max_depth) .& (offshoreregions .> 0) .& .!protected_area

    if plotmasks != false   # can == :onlymasks as well
        # drawmap(land)
        isregion = (regions .> 0) .& (regions .!= NOREGION)

        # mask values refer to colors in ColorBrewer Set2_7:
        # https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/#colorbrewer-1
        masks = zeros(Int16, size(regions))
        masks[(masks .== 0) .& (popdens .> persons_per_km2)] .= 2
        masks[(masks .== 0) .& protected_area] .= 3
        masks[(masks .== 0) .& .!gridA .& .!gridB] .= 4
        masks[(masks .== 0) .& .!goodland] .= 1
        masks[(masks .== 0) .& .!gridA .& gridB] .= 8
        masks[(masks .== 0) .& isregion] .= 7
        masks[regions .== 0] .= 0
        masks[regions .== NOREGION] .= NOREGION
        legendtext = ["bad land type", "high population", "protected area", "no grid", "", "", "wind plant A", "wind plant B"]
        maskmap("$(gisregion)_masks_wind", masks, legendtext, lonrange, latrange; legend=true, downsample=downsample)

        # drawmap(.!shore .& (offshoreregions .> 0))
        # drawmap((topo .> -max_depth) .& (offshoreregions .> 0))
        # drawmap(.!protected_area .& (offshoreregions .> 0))
        # drawmap(mask_offshore)
    end

    return mask_onshoreA, mask_onshoreB, mask_offshore,
            gridA, (popdens .< persons_per_km2), goodland, protected_area
end

function make_turbineCSV_for_postdocs(; mincapac=0, minclass=0, firstyear=1978, plotmasks=false, optionlist...)
    cols = [:lon, :lat, :capac, :year, :onshore, :elec2018]
    turbines = DataFrame(CSV.File(in_datafolder("turbines_DK.csv")))[:, cols]
    turbines = turbines[turbines[:capac].>=mincapac, :]
    nrows = size(turbines, 1)

    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, downsample_masks, onshore_density, area_onshore, res = options
    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land,
        protected_mat, lonrange, latrange =
                read_datasets(options)
    mask_onshoreA, mask_onshoreB, mask_offshore,
        gridA, lowpopdens, goodland, protected_area =
            return_wind_masks(options, regions, offshoreregions, gridaccess, popdens,
                            topo, land, protected_mat, lonrange, latrange,
                            plotmasks=plotmasks, downsample=downsample_masks)
    mask_onshore = mask_onshoreA .| mask_onshoreB
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)

    windatlas = getwindatlas()[lonrange,latrange]

    nreg = length(regionlist)

    pop, area, area_ok, maxcapac, maxcapac_ok =
        zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg), zeros(nreg)

    turbines[!, :municipality] = fill("", nrows)
    turbines[!, :avg_windspeed] = fill(0.0, nrows)
    turbines[!, :landtype] = fill("", nrows)
    turbines[!, :popdensity] = fill(0.0, nrows)
    turbines[!, :masks_ok] = fill(false, nrows)

    for row in eachrow(turbines)
        reg, flon, flat = getregion_and_index(row[:lon], row[:lat], regions, lonrange, latrange)
        if (reg == 0 || reg == NOREGION)
            reg, flon, flat = getregion_and_index(row[:lon], row[:lat], offshoreregions, lonrange, latrange)
            (reg == 0 || reg == NOREGION) && continue
        end
        row[:municipality] = string(regionlist[reg])
        row[:avg_windspeed] = windatlas[flon,flat]
        row[:landtype] = land[flon,flat] > 0 ? landtypes[land[flon,flat]] : "Water"
        row[:popdensity] = popdens[flon, flat]
        row[:masks_ok] = row[:onshore] ? mask_onshore[flon, flat] : mask_offshore[flon, flat]
    end

    for reg = 1:nreg
        area[reg] += sum((regions .== reg) .* cellarea')
        area_ok[reg] += sum(((regions .== reg) .& mask_onshore) .* cellarea')
        pop[reg] = sum((regions .== reg) .* popdens) / area[reg]
    end

    munic = DataFrame(region=regionlist, avg_popdens=pop, area=area, area_masks_ok=area_ok)
    CSV.write("turbines_DK.csv", turbines)
    CSV.write("municipalities_DK.csv", munic)
    turbines, munic
end

ext: jl

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
# Autogenerated wrapper script for Kissat_jll for powerpc64le-linux-gnu
export kissat, libkissat

JLLWrappers.@generate_wrapper_header("Kissat")
JLLWrappers.@declare_executable_product(kissat)
JLLWrappers.@declare_library_product(libkissat, "libkissat.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        kissat,
        "bin/kissat",
    )

    JLLWrappers.@init_library_product(
        libkissat,
        "lib/libkissat.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function uhlig06()
#
#  test polynomial used by F. Uhlig
#
    p = [1.0;0;0;0;0;0;0;0;-1];
    z = [-1.00000000000000
         -0.70710678118655 + 0.70710678118655im
         -0.70710678118655 - 0.70710678118655im
          0 + 1.00000000000000im
          0 - 1.00000000000000im
          1.00000000000000
          0.70710678118655 + 0.70710678118655im
          0.70710678118655 - 0.70710678118655im];
    z = [z ones(8)]

    p, PolyZeros(z)
end


text:
module EDSUtil

using Base.Cartesian: @nloops, @nexprs

using Docile
@document

export normalize, PCA, un_transform, eds_ϵ, eds_M, eds_cheap, get_eds_grid,
       eval_density, qnwgh, n_complete, complete_inds, complete_polynomial,
       complete_polynomial!, qnwmonomial

function selectperm(v::AbstractVector,
                    k::Union{Int, OrdinalRange};
                    lt::Function=isless,
                    by::Function=identity,
                    rev::Bool=false,
                    order::Base.Order.Ordering=Base.Order.Forward)
    select!(collect(1:length(v)), k,
            Base.Order.Perm(Base.Order.ord(lt,by,rev,order),v))
end

# TODO: look into using the Distances.jl package
"""
Compute pairwise distances between a single point `x` (represented as
a d-element vector) and `n` points `y` (represented as an d × n
matrix).

Distances are stored in third argument `out`, which should be an
`n`-element vector
"""
function _pairwise!(x::Vector, y::Matrix, out::Vector)
    n_out = size(y, 2)
    nd = length(x)

    for i=1:n_out
        tmp = 0.0
        for j = 1:nd
            @inbounds tmp += (x[j] - y[j, i])^2
        end
        out[i] = tmp
    end
end

function _dist_lt_eps(D, eps)
    n = 0
    for d in D
        if d <= eps
            n += 1
        end
    end
    n
end

"""
Demean and convert to unit variance a matrix along a specified dimension

Dimension argument is optional, with default argument value `1`, which
collapses rows (i.e. the mean of each column is subtracted from the
array)
"""
normalize(x::Matrix, dim::Int=1) = (x .- mean(x, dim)) ./ std(x, dim)

"""
Normalize `x` by subtracting mean of `other` and dividing by the
standard deviation of `other`. The mean and standard deviation are computed
along dimension `dim`, which has a default value of `1`.
"""
normalize(x::Matrix, other::Matrix, dim::Int=1) =
    (x .- mean(other, dim)) ./ std(other, dim)
# NOTE: I wish I had a monad to normalize and de-normalize...
# TODO: I should instead do something like
#       `with_normalize(f::Function, x::Matrix, other::Matrix, dim::Int=1)`
#       that will normalize, call f, then denormalize.
# TODO: I could also do this with a
#       `macro with_normalization(x, other, dim, Expr)` where `Expr` is a
#       begin/end block containing the routine to be done under normalization.
#       the advantage here would be that I'm not actually using higher order
#       functions, just automating the normalization boiler-plate.

"""
The `PCA` type is used to work with the principle components of matrix of
data. The fields are

- `data`: original data
- `V`: the `V` component of the svd of the normalized `data`
- `PC`: The principal components of normalized `data`
- `PCn`: Normalized principal components
"""
immutable PCA
    data::Matrix{Float64}
    V::Matrix{Float64}
    PC::Matrix{Float64}  #  principal components
    PCn::Matrix{Float64}  #normalized principal components
end

"Construct a `PCA` object, given data"
function PCA(data::Matrix)
    datan = normalize(data)   # normalize data
    U, S, V = svd(datan)      # compute svd
    PC = datan*V              # compute principal components (PCs)
    PCn = (PC ./ std(PC, 1))  # normalize PCs
    PCA(data, V, PC, PCn)     # return object
end

"""
Given a matrix `data` and a `PCA` representation `other`, compute a
principal components representation of `data` using the `svd` of other
"""
function PCA(data::Matrix, other::PCA)
    datan = normalize(data, other.data)
    PC = datan*other.V
    PCn = (PC ./ std(other.PC, 1))
    return PCA(data, other.V, PC, PCn)
end

"""
Given `data` and `other`, construct a principal components
representation of `data`, using the svd of `other`.
"""
PCA(data::Matrix, other::Matrix) = PCA(data, PCA(other))

"""
Given `x_pcn` and a PCA representation `pca`, undo the `pca` and
normalization transformation to map `x_pcn` onto the domain of `pca.data`

Note that the `n` stands for normalized

This is typically used in this type of example:

```julia
pca = PCA(data)              # compute principal components
x_pc = do_stuff(pca)         # operate on data in PC space
x = un_transform(x_pc, pca)  # map back into domain of data
```
"""
function un_transform(x_pcn::Matrix, pca::PCA)
    x_pc = x_pcn .* std(pca.PC, 1)                     # de-normalize PCs
    x_n = x_pc * pca.V'                                 # unto PCA
    x = (x_n .* std(pca.data, 1)) .+ mean(pca.data, 1)  # de-normalize z
    x
end

"""
Constructs an epsilon distinguishable set (EDS) for a given `data` set
and minimum distance `ϵ` using algorithm P^ϵ from page 10 of MM
"""
eds_ϵ(data::Matrix{Float64}, ϵ::Float64) = eds_ϵ(PCA(data), ϵ)

"""
Constructs an epsilon distinguishable set (EDS), given the normalized
`PCA` for the `data` (obtained via `PCA(data)`) and a scalar `ϵ`
"""
function eds_ϵ(data_PCA::PCA, ϵ::Float64)
    n, d = size(data_PCA.data)

    # normalize data mean zero, unit std
    PCn = data_PCA.PCn'  # do `'` so I can operate on columns at a time

    # set up iteration objects
    EDS_PCn = Array(Float64, d, n)
    eps2 = ϵ^2
    M = 0
    N1 = n

    D_i2 = Array(Float64, N1)

    inds = Array(Int, N1)
    while N1 > 0
        M += 1
        EDS_PCn[:, M] = PCn[:, 1]  # step 3: add x_u to P^ϵ
        _pairwise!(PCn[:, 1], PCn, D_i2)  # step 1: compute D(x_i,x_j) ∀ j
        ndrop = _dist_lt_eps(D_i2, eps2)  # compute number of points to drop
        sortperm!(inds, D_i2)  # get indices that sort distances
        PCnsort = PCn[:, inds]  # sort principal components by inds
        PCn = PCnsort[:, ndrop+1:end]  # drop ones that are too close
        N1 = size(PCn, 2)  # recompte size of PCs
        D_i2 = D_i2[1:N1]  # trim distance vector
        inds = inds[1:N1]  # drop index vector
    end

    # now transpose back and undo PCA+normalization transformation
    EDS_PCn = EDS_PCn[:, 1:M]'
    EDS = un_transform(EDS_PCn, data_PCA)
    return EDS
end

"""
Given `data`, construct an EDS with a target number of points `M`.
Implements algorithm Mbar from page 13 of MM. `ϵ1` and `ϵ2` should satisfy
`ϵ1 < ϵ2` and

```julia
size(eds_ϵ(data, ϵ1), 1) > M > size(eds_ϵ(data, ϵ2), 1)
```

If no values are given, then suitable initial conditions are chosen such that:

```julia
ϵ1 = rmin / (M^(1/d) - 1)
ϵ2 = 0.5rmax * M^(1/d)
```

where `rmin` and `rmax` are the minimal and maximal norm of the normalized
principal components of `data` and `d = size(data, 2)`
"""
function eds_M(data::Matrix{Float64}, M::Int, ϵ1::Float64=NaN, ϵ2::Float64=NaN)
    d = size(data, 2)
    if any(map(isnan, [ϵ1, ϵ2]))  # find initial ϵ bounds as described in paper
        data_PCA = PCA(data)
        rmin, rmax = extrema(mapslices(norm, data_PCA.PCn, [2]))
        ϵ1 = rmin / (M^(1/d) - 1)
        ϵ2 = rmax * M^(1/d)
    else
        data_PCA = PCA(data)
    end

    ϵ1 > ϵ2  && error("ϵ1 must be smaller than ϵ2")
    # because ϵ1 < ϵ2 we will have that eds_1 has more points than eds_2
    eds_1 = eds_ϵ(data_PCA, ϵ1)
    eds_2 = eds_ϵ(data_PCA, ϵ2)

    # make sure initial guesses form a bound
    size(eds_1, 1) < M && error("ϵ1 not small enough for requested M")
    size(eds_2, 1) > M && error("ϵ2 not large enough for requested M")

    # now proceed with algorithm
    old_M_ϵ = typemax(Int)
    while true
        ϵ = (ϵ1 + ϵ2) / 2
        eds = eds_ϵ(data_PCA, ϵ)
        M_ϵ = size(eds, 1)

        # if we didn't change, stop iterating
        if M_ϵ == old_M_ϵ
            eds_1 = eds
            break
        end

        # do bisection update
        if M_ϵ > M
            ϵ1 = ϵ
        else
            ϵ2 = ϵ
        end

        # update old_M_ϵ
        old_M_ϵ = M_ϵ

    end

    return eds_1
end

# Small helper functions to handle intermediate steps of creating the eds
get_eds_grid(vectors::Vector...; δ::Float64=0.01, Mbar::Int=50) =
    get_eds_grid(hcat(vectors...); δ=δ, Mbar=Mbar)

function get_eds_grid(sim_data::Matrix; δ::Float64=0.01, Mbar::Int=50)
    density = eval_density(sim_data, sim_data)[1]
    capT = size(sim_data, 1)
    inds = sortperm(density)
    n_drop = ceil(Int, capT*δ)
    sort_data = sim_data[inds, :]
    sort_data = sort_data[n_drop:end, :]

    # now construct eds on this "thinned" out data
    eds = eds_M(sort_data, Mbar)
end

"""
Construct an EDS of `data` using the cheap algorithm for computing the
essentially ergodic set. Given a set of simulated data `P`, the algorithm
proceeds as follows:

1. "Thin" the data by selecting every `κ`th point
2. Select and EDS P^ϵ of M points using `eds_ϵ` from above
3. Estimate the density at each point in P^ϵ using `eval_density`
4. Remove a fraction of points `δ` from `P` that has lowest density

Note, there is a subtlety here. We are removing points from `P`, not
`P^ϵ`. This means that each time we throw a point out of `P^ϵ`, the mass
of discarded points is not `1/size(P, 1)`, but rater the density of that
point (from step 3), over cumulative density of all `M` points in `P^ϵ`.
This means that we can toss points out of `P^ϵ` one by one until the
total mass of discarded points is equal to `δ`.

**NOTE**: In this case `data` is assumed to be the raw data obtained via
simulation, not the essentially ergodic set (as is assumed in the other
`eds_*` methods). This routine will optionally select every `κ`th point,
but this argument has a default value of `1`
"""
function eds_cheap(data::Matrix{Float64}, ϵ::Float64, δ::Float64, κ::Int=1)
    # TODO: implement this
    nothing
end

"""
TODO: write docstring
"""
function eds_locally_adaptive(data::Matrix{Float64}, ϵ::Vector{Float64})
    n, d = size(data)
    @assert n == length(ϵ) "length(ϵ) must be same as size(data, 1)"

    # TODO: implement this algorithm. It might take some API work because
    #       we need information on function approximation errors at each point
    nothing
end

# Somehow the apple libm is faster than the one julia calls on OSX
@osx? (
         begin
             myexp(x::Float64) = ccall((:exp, :libm), Float64, (Float64,), x)
         end
       : begin
             myexp(x::Float64) = exp(x)
         end
       )

"""
Estimate the density function in a given set of `points` by applying
a multivariate kernel algorithm to the `data`. See formula (2) in MM.

First argument is the data points, with `size(⋅) = (n, d)`

Second argument is the points at which to evaluate the density, with
`size(⋅) = (npoints, d)`

Output is three items:

1. `density`: Estimated density function with `size(⋅) = (npoints,)`
2. `points_PCn`: normalized principal components of the matrix of
   data for which the density was estimated with `size(⋅)=(npoints, d)`.
   This is needed by the EDS algorithm so it is returned here to reduce
   duplicate computation
3. `di_min`: The distance from each point in `data` to the closest
   neighbor with `size(⋅)=(npoints,)`. Also needed by the EDS
   algorithm so it is returned here to reduce duplicate computation
"""
function eval_density(data::Matrix, points::Matrix)

    n, d = size(data)

    n_points = size(points, 1)

    # transform data
    data_PCA = PCA(data)
    PC = data_PCA.PC
    PCn = data_PCA.PCn'  # do `'` so I can operate on columns at a time

    # transform points
    points_PCA = PCA(points, data_PCA)
    pts_PCn = points_PCA.PCn'  # same comment about `'`

    # Define constants
    hbar = n^(-1.0/(d+4.0))
    constant = 1/(n*hbar^(d) * (2π)^(d/2))
    constant2 = -0.5/hbar^2

    # allocate space
    density = Array(Float64, n_points)
    Di_min = Array(Float64, n_points)

    # apply formula (2)
    for i=1:n_points  # loop over all points
        dens_i = 0.0
        min_di2 = Inf
        for j=1:n_points  # loop over all other points
            d_i2_j = 0.0
            for k=1:d  # loop over d
                @inbounds d_i2_j += (pts_PCn[k, i] - PCn[k, j])^2
            end
            dens_i += myexp(d_i2_j*constant2)
            if i != j && d_i2_j < min_di2
                min_di2 = d_i2_j
            end
        end
        density[i] = dens_i
        Di_min[i] = sqrt(min_di2)
    end

    density .*= constant

    return density, points_PCA.PCn, Di_min
end

"""
WARNING: this function doesn't do anything yet...

Constructs `n` clusters and computes clusters' centers for a given
`data` set.

See Judd, Maliar and Maliar, (2010),
"A Cluster-Grid Projection Method: Solving Problems with High Dimensionality",
NBER Working Paper 15965 (henceforth, JMM, 2010).
"""
function find_clusters(data::Matrix, n::Int)
    T, M = size(data)
    data_PCA = PCA(data)
    PCn = data_PCA.PCn'  # `'` so I can operate on columns at a time

    # TODO: don't have `clusterdata` function
end

# --------------------------------------------------------- #
# Stuff to construct basis matrices of complete polynomials #
# --------------------------------------------------------- #

"""
Construct basis matrix for complete polynomial of degree `d`, given
input data `z`. `z` is assumed to be the degree 1 realization of each
variable. For example, if variables are `q`, `r`, and `s`, then `z`
should be `z = [q r s]`

Output is a basis matrix. In our example, with `d` set to 2 we would have

TODO: update docstring to properly give order of terms

```julia
out = [ones(size(z,1)) q r s q.*r q.*s r.*s q.^2 q.*r q.*s r.^2 r.*s s.^2]
```
"""
:complete_polynomial

immutable Degree{N} end

function n_complete(n::Int, D::Int)
    out = 1
    for d=1:D
        tmp = 1
        for j=0:d-1
            tmp *= (n+j)
        end
        out += div(tmp, factorial(d))
    end
    out
end

@generated function complete_polynomial!{T,N}(z::Matrix, d::Degree{N}, out::Matrix{T})
    complete_polynomial_impl!(z, d, out)
end

function complete_polynomial_impl!{T,N}(z::Type{Matrix{T}},
                                        deg::Type{Degree{N}},
                                        out::Type{Matrix{T}})
    quote
        nobs, nvar = size(z)
        if size(out) != (nobs, n_complete(nvar, $N))
            error("z, out not compatible")
        end

        # reset first column to ones
        @inbounds for i=1:nobs
            out[i, 1] = one($T)
        end

        # set next nvar columns to input matrix
        @inbounds for n=2:nvar+1, i=1:nobs
            out[i, n] = z[i, n-1]
        end

        ix = nvar+1
        @nloops($N, # number of loops
                i,  # counter
                d->((d == $N ? 1 : i_{d+1}) : nvar),  # ranges
                d->(d == $N ? nothing :
                    (begin
                        ix += 1
                        for r=1:nobs
                            tmp = one($T)
                            @nexprs $N-d+1 j->(tmp *= z[r, i_{$N-j+1}])
                            out[r, ix]=tmp
                        end
                    end)),  # preexpr
                Expr(:block, :nothing)  # bodyexpr
                )
        out
    end
end

function complete_polynomial{T}(z::Matrix{T}, d::Int)
    nobs, nvar = size(z)
    out = Array(T, nobs, n_complete(nvar, d))
    complete_polynomial!(z, Degree{d}(), out)
    out
end

function complete_polynomial!{T}(z::Matrix{T}, d::Int,
                                 out::Matrix{T})
    complete_polynomial!(z, Degree{d}(), out)
    out
end

# ------------------------------------------------- #
# Routines to generate quadrature weights and nodes #
# ------------------------------------------------- #

"""
Gauss Hermite quadrature notes and weights in N dimensions. Limited to
no more than 10 nodes in each dimension.

TODO: I really don't like this. I'll probably swap it out for one of my
      CompEcon routines once I have verified that my code gives the
      same answer as the Matlab
"""
function qnwgh(n::Int=10, d::Int=1, vcv::Matrix=eye(d))
    if n == 1
        ϵ = [0.0]
        ω = [sqrt(pi)]
    elseif n == 2
        ϵ = [0.7071067811865475, -0.7071067811865475]
        ω = [0.8862269254527580, 0.8862269254527580]
    elseif n == 3
        ϵ = [1.224744871391589, 0, -1.224744871391589]
        ω = [0.2954089751509193, 1.181635900603677, 0.2954089751509193]
    elseif n == 4
        ϵ = [1.650680123885785, 0.5246476232752903,  -0.5246476232752903,
             -1.650680123885785]
        ω = [0.08131283544724518, 0.8049140900055128, 0.8049140900055128,
             0.08131283544724518]
    elseif n == 5
        ϵ = [2.020182870456086, 0.9585724646138185, 0, -0.9585724646138185,
             -2.020182870456086]
        ω = [0.01995324205904591, 0.3936193231522412, 0.9453087204829419,
             0.3936193231522412, 0.01995324205904591]
    elseif n == 6
        ϵ = [2.350604973674492, 1.335849074013697, 0.4360774119276165,
             -0.4360774119276165, -1.335849074013697, -2.350604973674492]
        ω = [0.004530009905508846, 0.1570673203228566, 0.7246295952243925,
             0.7246295952243925, 0.1570673203228566, 0.004530009905508846]
    elseif n == 7
        ϵ = [2.651961356835233, 1.673551628767471, 0.8162878828589647, 0,
             -0.8162878828589647, -1.673551628767471, -2.651961356835233]
        ω = [0.0009717812450995192, 0.05451558281912703, 0.4256072526101278,
             0.8102646175568073, 0.4256072526101278, 0.05451558281912703,
             0.0009717812450995192]
    elseif n == 8
        ϵ = [2.930637420257244, 1.981656756695843, 1.157193712446780,
             0.3811869902073221, -0.3811869902073221, -1.157193712446780,
             -1.981656756695843,-2.930637420257244]
        ω = [0.0001996040722113676, 0.01707798300741348, 0.2078023258148919,
             0.6611470125582413, 0.6611470125582413, 0.2078023258148919,
             0.01707798300741348, 0.0001996040722113676]
    elseif n == 9
        ϵ = [3.190993201781528, 2.266580584531843, 1.468553289216668,
             0.7235510187528376, 0,-0.7235510187528376, -1.468553289216668,
             -2.266580584531843, -3.190993201781528]
        ω = [0.00003960697726326438, 0.004943624275536947,
             0.08847452739437657, 0.4326515590025558, 0.7202352156060510,
             0.4326515590025558, 0.08847452739437657, 0.004943624275536947,
             0.00003960697726326438]
    elseif n == 10
        ϵ = [3.436159118837738, 2.532731674232790, 1.756683649299882,
             1.036610829789514, 0.3429013272237046, -0.3429013272237046,
             -1.036610829789514, -1.756683649299882, -2.532731674232790,
             -3.436159118837738]
        ω = [7.640432855232621e-06, 0.001343645746781233,
             0.03387439445548106, 0.2401386110823147, 0.6108626337353258,
             0.6108626337353258, 0.2401386110823147, 0.03387439445548106,
             0.001343645746781233, 7.640432855232621e-06]
    else
        error("n must be between 1 and 10")
    end

    n_nodes = n^d
    z1 = zeros(n_nodes, d)
    ω1 = ones(n_nodes)

    for i=1:d
        ix = 1
        for j=1:n^(d-i)
            for u=1:n
                n_new = n^(i-1)
                z1[ix:ix+n_new-1, i] = ϵ[u]
                ω1[ix:ix+n_new-1] .*= ω[u]
                ix += n_new
            end
        end
    end

    z = sqrt(2) .* z1
    weights = ω1 ./ (sqrt(π)^d)
    nodes = ifelse(length(vcv) == 1, z*chol(vcv)[1], z*chol(vcv))::Matrix{Float64}
    return nodes, weights
end

"""
Computes integration nodes and weights under an N-dimensional monomial
(non-product) integration rule. If `kind` is equal to `:first` (the
default), then `2n` nodes will be computed, otherwise an algorithm
producing `2n^2+1` nodes is used.
"""
:qnwmonomial

qnwmonomial(n::Int, vcv::Matrix{Float64}, kind::Symbol=:first) =
    kind == :first ? _qnwmonomial1(n, vcv) : _qnwmonomial2(n, vcv)

function _qnwmonomial1(n::Int, vcv::Matrix{Float64})
    n_nodes = 2n

    z1 = zeros(n_nodes, n)

    # In each node, random variable i takes value either 1 or -1, and
    # all other variables take value 0. For example, for N = 2,
    # z1 = [1 0; -1 0; 0 1; 0 -1]
    for i=1:n
        z1[2*(i-1)+1:2*i, i] = [1, -1]
    end

    sqrt_vcv = chol(vcv)
    R = sqrt(n)*sqrt_vcv
    ϵj = z1*R
    ωj = ones(n_nodes) ./ n_nodes
    ϵj, ωj
end

function _qnwmonomial2(n::Int, vcv::Matrix{Float64})
    n_nodes = 2n^2 + 1
    z0 = zeros(1, n)

    z1 = zeros(2n, n)
    # In each node, random variable i takes value either 1 or -1, and
    # all other variables take value 0. For example, for N = 2,
    # z1 = [1 0; -1 0; 0 1; 0 -1]
    for i=1:n
        z1[2*(i-1)+1:2*i, i] = [1, -1]
    end

    z2 = zeros(2n*(n-1), n)
    i = 0

    # In each node, a pair of random variables (p,q) takes either values
    # (1,1) or (1,-1) or (-1,1) or (-1,-1), and all other variables take
    # value 0. For example, for N = 2, `z2 = [1 1; 1 -1; -1 1; -1 1]`
    for p=1:n-1
        for q=p+1:n
            i += 1
            z2[4*(i-1)+1:4*i, p] = [1, -1, 1, -1]
            z2[4*(i-1)+1:4*i, q] = [1, 1, -1, -1]
        end
    end

    sqrt_vcv = chol(vcv)
    R = sqrt(n+2)*sqrt_vcv
    S = sqrt((n+2)/2)*sqrt_vcv
    ϵj = [z0; z1*R; z2*S]
    ωj = vcat(2/(n+2) * ones(size(z0, 1)),
              (4-n)/(2*(n+2)^2) * ones(size(z1, 1)),
               1/(n+2)^2 * ones(size(z2, 1)))
    return ϵj, ωj
end

end  # module

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
#=
Here we solve the equation:
```math
 q + dot(∇, uq) = 0
 p - dot(∇, up) = 0
```
on a sphere to test the conservation of the numerics

The boundary conditions are `p = q` when `dot(n, u) > 0` and
`q = p` when `dot(n, u) < 0` (i.e., `p` flows into `q` and vice-sersa).
=#

using MPI
using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using Random

using ClimateMachine.VariableTemplates
import ClimateMachine.DGMethods:
    BalanceLaw,
    vars_state_auxiliary,
    vars_state_conservative,
    vars_state_gradient,
    vars_state_gradient_flux,
    flux_first_order!,
    flux_second_order!,
    source!,
    boundary_state!,
    init_state_auxiliary!,
    init_state_conservative!,
    init_ode_state,
    LocalGeometry
import ClimateMachine.DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder,
    numerical_flux_first_order!,
    numerical_boundary_flux_first_order!

struct ConservationTestModel <: BalanceLaw end

vars_state_auxiliary(::ConservationTestModel, T) = @vars(vel::SVector{3, T})
vars_state_conservative(::ConservationTestModel, T) = @vars(q::T, p::T)

vars_state_gradient(::ConservationTestModel, T) = @vars()
vars_state_gradient_flux(::ConservationTestModel, T) = @vars()

function init_state_auxiliary!(
    ::ConservationTestModel,
    aux::Vars,
    g::LocalGeometry,
)
    x, y, z = g.coord
    r = x^2 + y^2 + z^2
    aux.vel = SVector(
        cos(10 * π * x) * sin(10 * π * y) + cos(20 * π * z),
        exp(sin(π * r)),
        sin(π * (x + y + z)),
    )
end

function init_state_conservative!(
    ::ConservationTestModel,
    state::Vars,
    aux::Vars,
    coord,
    t,
)
    state.q = rand()
    state.p = rand()
end

function flux_first_order!(
    ::ConservationTestModel,
    flux::Grad,
    state::Vars,
    auxstate::Vars,
    t::Real,
)
    vel = auxstate.vel
    flux.q = state.q .* vel
    flux.p = -state.p .* vel
end

flux_second_order!(::ConservationTestModel, _...) = nothing

source!(::ConservationTestModel, _...) = nothing

struct ConservationTestModelNumFlux <: NumericalFluxFirstOrder end

boundary_state!(
    ::CentralNumericalFluxSecondOrder,
    ::ConservationTestModel,
    _...,
) = nothing

function numerical_flux_first_order!(
    ::ConservationTestModelNumFlux,
    bl::BalanceLaw,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
) where {S, A}
    un⁻ = dot(n, aux⁻.vel)
    un⁺ = dot(n, aux⁺.vel)
    un = (un⁺ + un⁻) / 2

    if un > 0
        fluxᵀn.q = un * state⁻.q
        fluxᵀn.p = -un * state⁺.p
    else
        fluxᵀn.q = un * state⁺.q
        fluxᵀn.p = -un * state⁻.p
    end
end

function numerical_boundary_flux_first_order!(
    ::ConservationTestModelNumFlux,
    bl::BalanceLaw,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    bctype,
    t,
    state1⁻::Vars{S},
    aux1⁻::Vars{A},
) where {S, A}
    un = dot(n, aux⁻.vel)

    if un > 0
        fluxᵀn.q = un * state⁻.q
        fluxᵀn.p = -un * state⁻.q
    else
        fluxᵀn.q = un * state⁻.p
        fluxᵀn.p = -un * state⁻.p
    end
end

function run(mpicomm, ArrayType, N, Nhorz, Rrange, timeend, FT, dt)

    topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
    )
    dg = DGModel(
        ConservationTestModel(),
        grid,
        ConservationTestModelNumFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0); init_on_cpu = true)

    lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

    eng0 = norm(Q)
    sum0 = weightedsum(Q)
    @info @sprintf """Starting
    norm(Q₀) = %.16e
    sum(Q₀) = %.16e""" eng0 sum0

    max_mass_loss = FT(0)
    max_mass_gain = FT(0)
    cbmass = GenericCallbacks.EveryXSimulationSteps(1) do
        cbsum = weightedsum(Q)
        max_mass_loss = max(max_mass_loss, sum0 - cbsum)
        max_mass_gain = max(max_mass_gain, cbsum - sum0)
    end
    solve!(Q, lsrk; timeend = timeend, callbacks = (cbmass,))

    # Print some end of the simulation information
    engf = norm(Q)
    sumf = weightedsum(Q)
    @info @sprintf """Finished
    norm(Q)            = %.16e
    norm(Q) / norm(Q₀) = %.16e
    norm(Q) - norm(Q₀) = %.16e
    max mass loss      = %.16e
    max mass gain      = %.16e
    initial mass       = %.16e
    """ engf engf / eng0 engf - eng0 max_mass_loss max_mass_gain sum0
    max(max_mass_loss, max_mass_gain) / sum0
end

using Test
let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4

    Nhorz = 4

    tolerance = Dict(Float64 => 1e-15, Float32 => 1e-7)

    @testset "$(@__FILE__)" for FT in (Float64, Float32)
        dt = FT(1e-4)
        timeend = 100 * dt
        Rrange = range(FT(1), stop = FT(2), step = FT(1 // 4))

        @info (ArrayType, FT)
        delta_mass = run(
            mpicomm,
            ArrayType,
            polynomialorder,
            Nhorz,
            Rrange,
            timeend,
            FT,
            dt,
        )
        @test abs(delta_mass) < tolerance[FT]
    end
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
data = parse.(Int,permutedims(reduce(hcat, collect.(readlines(joinpath(@__DIR__, "./3.txt"))))))

function bitarr_to_int(arr)
  return sum(arr .* (2 .^ collect(length(arr)-1:-1:0)))
end

function gamma(arr)
  return BitArray(map((s) -> s >= size(arr)[1]/2,  sum(arr,dims=1)))
end

function epsilon(arr)
  return .!gamma(arr)
end

g = gamma(data)
e = epsilon(data)

# answer to part 1
bitarr_to_int(g') * bitarr_to_int(e')

function rating(data,rate_func,index=1)
  if size(data,1) == 1
    return bitarr_to_int(BitArray(data)')
  end
  rate = rate_func(data)
  filtered_data  = data[data[:,index].==rate[index],:]
  return rating(filtered_data, rate_func, index += 1)
end

o2 = rating(data, gamma)
co2 = rating(data, epsilon)

#answer to part 2
o2*co2


text:
# TODO: abstract to objects (spheres) with Transmission and Reflectivity
# TODO: rasterize rays

using Plots
using LinearAlgebra
using Test

pyplot()

function Base.isapprox(n1::Nothing, n2::Nothing)
    return true
end

abstract type Object end

mutable struct Ray
    # Velocity vector
    v::Vector{Float64}

    # Position vector
    p::Vector{Float64}

end

mutable struct Lens <: Object
    # Lens position
    p::Vector{Float64}

    # Lens radius
    r::Float64

    # Index of Refraction
    ior::Float64
end

function lens_normal_at(ray, lens)
    n = normalize(ray.p .- lens.p)
    if dot(-n, ray.v) < 0
        n *= -1
    end

    return n
end

function inside_of(ray::Ray, lens)
    return inside_of(ray.p, lens)
end

function inside_of(pos, lens)
    x = lens.p[1] - pos[1]
    y = lens.p[2] - pos[2]
    if (x^2 + y^2 <= lens.r^2)
        return true
    else
        return false
    end
end

# note: light moves as a particular speed with respect to the medium it is
#       moving through, so...
#       n_2*v = n_1*l + (n_1*cos(theta_1) - n_2*cos(theta_2))*n
#       Other approximations: ior = n_1/n_2, c = -n*l
function refract!(ray, lens, ior)
    n = lens_normal_at(ray, lens)
    c = dot(-n, ray.v);
    d = 1.0 - ior^2 * (1.0 - c^2);

    if (d < 0.0)
        return zeros(2);
    end

    ray.v = ior * ray.v + (ior * c - sqrt(d)) * n;
end

abstract type Wall <: Object end;

mutable struct Mirror <: Wall
    # Normal vector
    n::Vector{Float64}

    # Position of mirror
    p::Vector{Float64}

    # Mirror size
    scale::Float64

    Mirror(in_n, in_p) = new(in_n, in_p, 2.5)
end

function is_behind(ray, mirror)
    if dot(ray.p.-mirror.p, mirror.n) >= 0
        return true
    else
        return false
    end
end

# note: for reflection, l_x -> l_x, but l_y -> -l_y
#       In this case, the y component of l = cos(theta)*n
#       so new vector: v = l + 2cos(theta)*n
function reflect!(ray, n)
    ray.v .-= 2*dot(ray.v, n).*n
end

function draw_circle(p, r, res)
    return [x .+ (p[1], p[2]) for x in Plots.partialcircle(0, 2pi, res, r)]
end

function plot_rays(positions, objects::Vector{O},
                   filename) where {O <: Object}
    plt = plot(background_color=:black, aspect_ratio=:equal)

    for i = 1:size(positions)[1]
        plot!(plt, (positions[i,:,1], positions[i,:,2]); label = "ray",
              linecolor=:white)
    end

    for object in objects
        if typeof(object) == Mirror
            dir = [-object.n[2], object.n[1]]
            extents = zeros(2,2)
            extents[1,:] .= object.p .- object.scale * dir
            extents[2,:] .= object.p .+ object.scale * dir
            plot!(plt, (extents[:,1], extents[:,2]), label = "mirror")
        elseif typeof(object) == Lens
            circle = draw_circle(object.p, object.r, 100)
            plot!(circle; label="lens", linecolor=:lightblue)
        end
    end

    savefig(plt, filename)
end

function step!(ray::Ray, dt)
    ray.p .+= .+ ray.v.*dt
end

function intersection(ray::Ray, lens::Lens)
    # threshold to avoid floating precision errors
    thresh = 0.00001

    # distance from the center of the lens
    relative_dist = lens.p - ray.p

    # distance projected onto the direction vector
    projected_relative_dist = dot(relative_dist, ray.v)

    # distance of closest approach
    # Might make more sense to do a dot product instead of a norm
    closest_approach = abs(sqrt(norm(relative_dist)^2 -
                            projected_relative_dist^2))

    # TODO: arbitrary threshold, should be fixed.
    if abs(closest_approach - lens.r) < thresh
        closest_approach = lens.r
    end

    # Note: returns nothing if tangential
    if closest_approach >= lens.r
        return nothing
    end

    # half chord length
    half_chord = sqrt(lens.r^2 - closest_approach^2)

    # checks if ray is moving away from sphere by seeing if we are
    # outside of sphere (abs(projected_relative_dist) > half_chord), and if
    # the ray is moving away from the boundary
    if (projected_relative_dist < 0 &&
        abs(projected_relative_dist) > half_chord)
        return nothing
    end

    in_lens = half_chord < projected_relative_dist
    intersect = zeros(2)
    # Adding a threshold to output so we are definitely on the other side of
    # the boundary after intersecting
    if in_lens
        intersect = (projected_relative_dist - half_chord+thresh)*ray.v
    else
        intersect = (projected_relative_dist + half_chord+thresh)*ray.v
    end

    return intersect
end

function intersection_quadratic(ray::Ray, lens::Lens)
    relative_dist = ray.p-lens.p
    a = dot(ray.v, ray.v)
    b = 2.0 * dot(relative_dist, ray.v)
    c = dot(relative_dist, relative_dist) - lens.r*lens.r
    discriminant = b*b - 4*a*c

    if discriminant < 0
        return nothing
    elseif discriminant > 0
        roots = [(-b + sqrt(discriminant)) / (2*a),
                 (-b - sqrt(discriminant)) / (2*a)]
        min = minimum(roots)
        max = maximum(roots)

        if min > 0
            return (min)*ray.v
        elseif max > 0
            return (max)*ray.v
        else
            return nothing
        end
    else
        # Returns nothing if tangential
        return nothing
        #return (-b/(2*a))*ray.v
    end
end

function intersection(ray::Ray, wall::W) where {W <: Wall}
    intersection = -dot((ray.p .- wall.p),wall.n)/dot(ray.v, wall.n)

    if isfinite(intersection) && intersection > 0 && intersection != NaN
        return intersection*ray.v
    else
        return nothing
    end
end

function intersect_test()

    lens = Lens([3.0, 1.0], 1.0, 1.5)

    @testset "intersect lens tests" begin
        rays = [Ray([1.0, 0.0], [0.0, 1.0], zeros(2,2)),
                Ray([1.0, 0.0], [3.0, 1.0], zeros(2,2)),
                Ray([sqrt(0.5), sqrt(0.5)], [2.0, 0.0], zeros(2,2)),
                Ray([-1.0, 0.0], [0.0, 1.0], zeros(2,2)),
                Ray([0.0, 1.0], [0.0, 1.0], zeros(2,2)),
                Ray([1.0, 0.0], [0.0, 0.0], zeros(2,2))]
        answers = [[2.0, 0.0], [1.0, 0.0], [1-sqrt(0.5), 1-sqrt(0.5)],
                   nothing, nothing, [3.0, 0.0]]

        for i = 1:length(rays)
            pt = intersection(rays[i], lens)
            @test isapprox(pt, answers[i])
            pt = intersection_geometric(rays[i], lens)
            @test isapprox(pt, answers[i])
        end
    end

    wall = Mirror([-1.0, 0.0], [2.0, 1.0])

    @testset "intersect wall tests" begin
        rays = [Ray([1.0, 0.0], [0.0, 1.0], zeros(2,2)),
                Ray([sqrt(0.5), sqrt(0.5)], [0.0, 0.0], zeros(2,2)),
                Ray([-sqrt(0.5), sqrt(0.5)], [0.0, 0.0], zeros(2,2)),
                Ray([0.0, 1.0], [0.0, 1.0], zeros(2,2))]

        answers = [[2.0, 0.0], [2.0, 2.0], nothing, nothing]

        for i = 1:length(rays)
            pt = intersection(rays[i], wall)
            @test isapprox(pt, answers[i])
        end
    end
end

function propagate!(rays::Vector{Ray}, objects::Vector{O},
                    num_intersections, positions) where {O <: Object}
    for i = 2:num_intersections
        for j = 1:length(rays)

            intersect_final = [Inf, Inf]
            intersected_object = nothing
            for object in objects
                intersect = intersection(rays[j], object)
                if intersect != nothing &&
                   sum(intersect[:].^2) < sum(intersect_final[:].^2)
                    intersect_final = intersect
                    intersected_object = object
                end
            end

            if intersect_final != nothing
                rays[j].p .+= intersect_final
                positions[j, i, :] .= rays[j].p
                if typeof(intersected_object) == Lens
                    ior = intersected_object.ior
                    if inside_of(rays[j], intersected_object)
                        ior = 1/intersected_object.ior
                    end

                    refract!(rays[j], intersected_object, ior)
                    if rays[j].v == zeros(2)
                        reflect!(rays[j], lens_normal_at(rays[j], lens))
                    end

                elseif typeof(intersected_object) == Mirror
                    reflect!(rays[j], intersected_object.n)
                end

                println(rays[j].v, '\t', rays[j].p)
            end
        end
        println('\n')
    end
end

function parallel_propagate(ray_num, num_intersections, box_size;
                            filename="check.png")

    rays = [Ray([1, 0],
            [0.1, float(i)]) for i = 1:ray_num]

    lenses = [Lens([10, 5], 5, 1.5)]
    mirrors = [Mirror([-1.0, 0.0], [25.0, 5.0]),
               Mirror([1.0, 0.0], [0, 5.0]),
               Mirror([0.0, 1.0], [12.5, 0.0]),
               Mirror([0.0, -1.0], [12.5, 10.0])]

    objects = vcat(lenses, mirrors)
    positions = zeros(ray_num, num_intersections, 2)

    for i = 1:length(rays)
        positions[i, 1, :] .= rays[i].p
    end

    propagate!(rays, objects, num_intersections, positions)

    plot_rays(positions, objects, filename)

    return positions

end


text:
print("Please Tell Me Your Name: ")
a = readline()
println("hello $(a)")


text:
"""
    @check(expr, [message]) -> true

Throws an error with an optionally specified `message` if the given code
expression evaluates to `false`, otherwise returns `true`.
"""
macro check(expr, msgs...)
    msg_str = isempty(msgs) ? string(expr) : msgs[1]
    esc(:($expr || throw(AssertionError($msg_str))))
end

"""
    @checkeq(left, right, [message]) -> true

Throws an error with an optionally specified `message` if the given `left` code
expression does not evaluates to the same value as the given `right`
expression, otherwise returns `true`.
"""
macro checkeq(left, right, msgs...)
    msg_str = isempty(msgs) ? "$(left) == $(right)" : msgs[1]
    esc(:($left == $right || throw(AssertionError($msg_str))))
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using Hecke, DelimitedFiles

boundcond = ARGS[1]
gtype = ARGS[2]

_gtype = replace(replace(replace(replace(replace(gtype, '[' => '_'), ']' => '_'), ' ' => '_'), ',' => '_'), "__" => "_")

file = "conductors_" * boundcond * "_" * _gtype;

bound = fmpz(Meta.eval(Meta.parse(boundcond)))
gtype = convert(Vector{Int}, Meta.eval(Meta.parse(gtype)))

sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

gtype_right = snf(abelian_group(gtype))[1].snf
g_exp = gtype_right[end]
wild_primes = factor(g_exp).fac
for (p, v) in wild_primes
  if p == 2
    wild_primes[p] = v+2
  else
    wild_primes[p] = v+1
  end
end
conductors = Vector{Int}()
for i = 3:bound
  lf = factor(i)
  admissible = true
  for (p, v) in lf
    if p in keys(wild_primes)
      if v > 1 && v <= wild_primes[p]
        continue
      elseif v > wild_primes[p]
        admissible = false
        break
      else
        #we need to check that p-1 is not coprime to the other factors
        if length(wild_primes) == 1
          admissible = false
          break
        else
          t = p-1
          found = false
          for q in keys(wild_primes)
            if !is_coprime(t, q)
              found = true
              break
            end
          end
          if !found
            admissible = false
            break
          end
        end
      end
    else
      if v > 1
        admissible = false
        break
      end
      #p is tamely ramified...
      t = p-1
      found = false
      for q in keys(wild_primes)
        if !is_coprime(t, q)
          found = true
          break
        end
      end
      if !found
        admissible = false
        break
      end
    end
  end
  if admissible
    push!(conductors, i)
  end
end

println("DONE")

if length(ARGS) == 3
  file = ARGS[3]
else
  file = file * "_" * string(length(conductors))
end

if isfile(file)
  error("File $file already exists")
end

open(file, "w") do f
  print("Writing conductors to $file ... ")
  writedlm(f, conductors)
  println("DONE")
end

@assert countlines(file) == length(conductors)


text:
module BilevelJuMP

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
using JuMP
using Dualization
using LinearAlgebra
using IntervalArithmetic

export
BilevelModel, UpperToLower, LowerToUpper,
Upper, Lower, UpperOnly, LowerOnly,
DualOf

include("moi.jl")
include("jump.jl")

end


text:
using StatsPlots, Random
gr(size=(800,600))

ProjDir = @__DIR__

function chart(N::Int)

  shapes = range(0.5, stop=2.0, length=N)
  scales = range(0.2, stop=2.5, length=N)

  p = Array{Plots.Plot{Plots.GRBackend},2}(undef, N, N);
  for (shi,sh) in enumerate(shapes)
    for (sci, sc) in enumerate(scales)
      dgamma = Gamma(sh, sc)
      p[sci, shi] = density(rand(dgamma, 10000),
        title="Shape $(round(sh, digits=1)),
        scale $(round(sc, digits=1))",
        titlefontsize=6)
    end
  end

  plot(p..., layout=Plots.GridLayout(N, N), legend=false, yaxis=true, tickfontsize=6)
  savefig("$(ProjDir)/chart.png")

  plot(density(rand(Gamma(3, 15), 10000)), lab="Gamma(3, 15)")
  savefig("$(ProjDir)/chart_a.png")

end

chart(6)


text:
module Advent2020

using Dates: day, today
using Reexport: @reexport

include("FileIO.jl")
include("Utils.jl")
include("Machines.jl")
include(joinpath("EventUtils", "EventUtils.jl"))

@reexport using .EventUtils

for i in 1:25
    dir = joinpath(@__DIR__, "..", "scripts", "Day"*string(i))
    if isdir(dir)
        include(joinpath(dir, "code.jl"))
    else
        make_my_day(i)
    end
end

end