text:
"""
ALL_FT+ fileformat

Defines a list of flight trajectories (planned/actual/etc) for a given day

See EUROCONTROL NEST Manual Section 9.7.5 for ALL_FT+ fileformat description

The headers of the DataFrame follow the naming convention from the EUROCONTROL
NEST Manual + '_' + index_number from manual.
"""
module Allftplus

export read

# include("utility.jl")
using ..util
using Format
using CSV
using Dates
using DataFrames

const fileformat = Dict(1=>String, 2=>String, 3=>String, 4=>String, 5=>String,
6=>String, 7=>String, 8=>String, 9=>String, 10=>String, 11=>String, 12=>String,
13=>String, 14=>String, 15=>String, 16=>String, 17=>String, 18=>String,
19=>String, 20=>String, 21=>String, 22=>String, 23=>Int64, 24=>String,
25=>String, 26=>String, 27=>String, 28=>String, 29=>String, 30=>Int64,
31=>Int64, 32=>String, 33=>String, 34=>String, 35=>String, 36=>String,
37=>String, 38=>String, 39=>Int64, 40=>Int64, 41=>Int64, 42=>String,
43=>String, 44=>String, 45=>String, 46=>String, 47=>String, 48=>String,
49=>Int64, 50=>Int64, 51=>String, 52=>String, 53=>String, 54=>Int64,
55=>Int64, 56=>String, 57=>String, 58=>String, 59=>String, 60=>String,
61=>String, 62=>String, 63=>String, 64=>String, 65=>String, 66=>String,
67=>String, 68=>String, 69=>String, 70=>String, 71=>String, 72=>String,
73=>Int64, 74=>String, 75=>Int64, 76=>String, 77=>Int64, 78=>Int64,
79=>String, 80=>String, 81=>Int64, 82=>String, 83=>Float64, 84=>Float64,
85=>Int64, 86=>String, 87=>Int64, 88=>String, 89=>Int64, 90=>String,
91=>Int64, 92=>Int64, 93=>String, 94=>String, 95=>Int64, 96=>String,
97=>Float64, 98=>Float64, 99=>Int64, 100=>String, 101=>Int64, 102=>String,
103=>Int64, 104=>String, 105=>Int64, 106=>Int64, 107=>String, 108=>String,
109=>Int64, 110=>String, 111=>Float64, 112=>Float64, 113=>Int64, 114=>String,
115=>Int64, 116=>String, 117=>Int64, 118=>String, 119=>String, 120=>String,
121=>Float64, 122=>Float64, 123=>Int64, 124=>String, 125=>Int64, 126=>String,
127=>Int64, 128=>String, 129=>String, 130=>Float64, 131=>Float64, 132=>Int64,
133=>String, 134=>Int64, 135=>String, 136=>Int64, 137=>String, 138=>String,
139=>Float64, 140=>Float64, 141=>Int64, 142=>String, 143=>Int64, 144=>String,
145=>Int64, 146=>String, 147=>String, 148=>Float64, 149=>Float64, 150=>Int64,
151=>String, 152=>Int64, 153=>String, 154=>Int64, 155=>String, 156=>String,
157=>Float64, 158=>Float64, 159=>Int64, 160=>String, 161=>Int64, 162=>String,
163=>Int64, 164=>String, 165=>String, 166=>String, 167=>String, 168=>Int64,
169=>String, 170=>String, 171=>String, 172=>String, 173=>String, 174=>String,
175=>Int64, 176=>Int64, 177=>String, 178=>Int64, 179=>String, 180=>Int64,
181=>String)

const header_format = ["departureAerodromeIcaoId_0", "arrivalAerodromeIcaoId_1",
"aircraftId_2", "aircraftOperatorIcaoId_3", "aircraftTypeIcaoId_4", "aobt_5",
"ifpsId_6", "iobt_7", "originalFlightDataQuality_8", "flightDataQuality_9",
"source_10", "exemptionReasonType_11", "exemptionReasonDistance_12", "lateFiler_13",
"lateUpdater_14", "northAtlanticFlight_15", "cobt_16", "eobt_17", "lobt_18",
"flightState_19", "previousToActivationFlightState_20", "suspensionStatus_21",
"tactId_22", "samCtot_23", "samSent_24", "sipCtot_25", "sipSent_26", "slotForced_27",
"mostPenalizingRegulationId_28", "regulationsAffectedByNrOfInstances_29",
"excludedFromNrOfInstances_30", "lastReceivedAtfmMessageTitle_31",
"lastReceivedMessageTitle_32", "lastSentAtfmMessageTitle_33",
"manualExemptionReason_34", "sensitiveFlight_35", "readyForImprovement_36",
"readyToDepart_37", "revisedTaxiTime_38", "tis_39", "trs_40",
"toBeSentSlotMessageTitle_41", "toBeSentProposalMessageTitle_42",
"lastSentSlotMessageTitle_43", "lastSentProposalMessageTitle_44",
"lastSentSlotMessage_45", "lastSentProposalMessage_46", "flightCountOption_47",
"normalFlightTactId_48", "proposalFlightTactId_49",
"operatingAircraftOperatorIcaoId_50", "reroutingWhy_51", "reroutedFlightState_52",
"runwayVisualRange_53", "numberIgnoredErrors_54", "arcAddrSource_55", "arcAddr_56",
"ifpsRegistrationMark_57", "flightType_58", "aircraftEquipment_59", "cdmStatus_60",
"cdmEarlyTtot_61", "cdmAoTtot_62", "cdmAtcTtot_63", "cdmSequencedTtot_64",
"cdmTaxiTime_65", "cdmOffBlockTimeDiscrepancy_66", "cdmDepartureProcedureId_67",
"cdmAircraftTypeId_68", "cdmRegistrationMark_69", "cdmNoSlotBefore_70",
"cdmDepartureStatus_71", "ftfmEetFirNrOfInstances_72", "ftfmEetFirList_73",
"ftfmEetPtNrOfInstances_74", "ftfmEetPtList_75", "ftfmAiracCycleReleaseNumber_76",
"ftfmEnvBaselineNumber_77", "ftfmDepartureRunway_78", "ftfmArrivalRunway_79",
"ftfmReqFlightlevelSpeedNrOfInstances_80", "ftfmReqFlightlevelSpeedList_81",
"ftfmConsumedFuel_82", "ftfmRouteCharges_83", "ftfmAllFtPointNrOfInstances_84",
"ftfmAllFtPointProfile_85", "ftfmAllFtAirspaceNrOfInstances_86",
"ftfmAllFtAirspaceProfile_87", "ftfmAllFtCircleIntersectionsNrOfInstances_88",
"ftfmAllFtCircleIntersections_89", "rtfmAiracCycleReleaseNumber_90",
"rtfmEnvBaselineNumber_91", "rtfmDepartureRunway_92", "rtfmArrivalRunway_93",
"rtfmReqFlightlevelSpeedNrOfInstances_94", "rtfmReqFlightlevelSpeedList_95",
"rtfmConsumedFuel_96", "rtfmRouteCharges_97", "rtfmAllFtPointNrOfInstances_98",
"rtfmAllFtPointProfile_99", "rtfmAllFtAirspaceNrOfInstances_100",
"rtfmAllFtAirspaceProfile_101", "rtfmAllFtCircleIntersectionsNrOfInstances_102",
"rtfmAllFtCircleIntersections_103", "ctfmAiracCycleReleaseNumber_104",
"ctfmEnvBaselineNumber_105", "ctfmDepartureRunway_106", "ctfmArrivalRunway_107",
"ctfmReqFlightlevelSpeedNrOfInstances_108", "ctfmReqFlightlevelSpeedList_109",
"ctfmConsumedFuel_110", "ctfmRouteCharges_111", "ctfmAllFtPointNrOfInstances_112",
"ctfmAllFtPointProfile_113", "ctfmAllFtAirspaceNrOfInstances_114",
"ctfmAllFtAirspaceProfile_115", "ctfmAllFtCircleIntersectionsNrOfInstances_116",
"ctfmAllFtCircleIntersections_117", "noCPGCPFReason_118", "scrObt_119",
"scrConsumedFuel_120", "scrRouteCharges_121", "scrAllFtPointNrOfInstances_122",
"scrAllFtPointProfile_123", "scrAllFtAirspaceNrOfInstances_124",
"scrAllFtAirspaceProfile_125", "scrAllFtCircleIntersectionsNrOfInstances_126",
"scrAllFtCircleIntersections_127", "srrObt_128", "srrConsumedFuel_129",
"srrRouteCharges_130", "srrAllFtPointNrOfInstances_131", "srrAllFtPointProfile_132",
"srrAllFtAirspaceNrOfInstances_133", "srrAllFtAirspaceProfile_134",
"srrAllFtCircleIntersectionsNrOfInstances_135", "srrAllFtCircleIntersections_136",
"surObt_137", "surConsumedFuel_138", "surRouteCharges_139",
"surAllFtPointNrOfInstances_140", "surAllFtPointProfile_141",
"surAllFtAirspaceNrOfInstances_142", "surAllFtAirspaceProfile_143",
"surAllFtCircleIntersectionsNrOfInstances_144", "surAllFtCircleIntersections_145",
"dctObt_146", "dctConsumedFuel_147", "dctRouteCharges_148",
"dctAllFtPointNrOfInstances_149", "dctAllFtPointProfile_150",
"dctAllFtAirspaceNrOfInstances_151", "dctAllFtAirspaceProfile_152",
"dctAllFtCircleIntersectionsNrOfInstances_153", "dctAllFtCircleIntersections_154",
"cpfObt_155", "cpfConsumedFuel_156", "cpfRouteCharges_157",
"cpfAllFtPointNrOfInstances_158", "cpfAllFtPointProfile_159",
"cpfAllFtAirspaceNrOfInstances_160", "cpfAllFtAirspaceProfile_161",
"cpfAllFtCircleIntersectionsNrOfInstances_162", "cpfAllFtCircleIntersections_163",
"aircraftidIATA_164", "intentionFlight_165",
"intentionRelatedRouteAssignmentMethod_166", "intentionUID_167",
"intentionEditionDate_168", "intentionSource_169", "associatedIntentions_170",
"enrichmentOutput_171", "eventID_172", "eventTime_173", "flightVersionNr_174",
"ftfmNrTvProfiles_175", "ftfmTvProfile_176", "rtfmNrTvProfiles_177",
"rtfmTvProfile_178", "ctfmNrTvProfiles_179", "ctfmTvProfile_180"]

const yymmdd = DateFormat("YYmmdd")
const hhmmss = DateFormat("HHMMSS")
const mmmmss = DateFormat("MMMMSS")
const yyyymmddhhmmss = DateFormat("YYYYmmddHHMMSS")
const year2000 = Year(2000)

function read(file)
    df = CSV.read(file, delim=";", types=fileformat, header=header_format,
    copycols=true, datarow=2)
    # remove_unused!(df)
    reformat!(df)
    return df
end

struct FlightlevelSpeed
    FL::Int64
    Spd::String
    Value::Int64
    function FlightlevelSpeed(x::AbstractString)
        flspdvalue = split(x, ':')
        FL = parse(Int64, flspdvalue[1][2:end])
        Spd = flspdvalue[2]
        Value = parse(Int64, flspdvalue[3])
        new(FL, Spd, Value)
    end
end

function reqFlightlevelSpeedList(items::Union{AbstractString, Missing})
    if items === missing
        return missing
    else
        return [FlightlevelSpeed(item) for item in split(items)]
    end
end

struct AllFtPointProfile
    datetime::Union{DateTime, Missing}
    point::AbstractString
    route::AbstractString
    FL::Union{Int64, Missing}
    pointDistance::Union{Int64, Missing}
    pointType::AbstractString
    geoPointId::Union{Point_deg, AbstractString}
    ratio::Union{Int64, Missing}
    isVisible::Bool  # Y indicated IFR/GAT/IFPSTART, N indicates VFR/OAT/IFPSTOP/STAY
    function AllFtPointProfile(x::AbstractString)
        items = split(x, ':')
        datetime = items[1] == "" ? missing : format_datetime(items[1], yyyymmddhhmmss)
        point = items[2]
        route = items[3]
        FL = items[4] == "" ? missing : parse(Int64, items[4])
        pointDistance = items[5] == "" ? missing : parse(Int64, items[5])
        pointType = items[6]
        geoPointId = occursin(r"\d{6}[NS]\d{7}[EW]", items[7]) ? latlon(items[7]) : items[7]
        ratio = items[8] == "" ? missing : parse(Int64, items[8])
        isVisible = items[9] == "Y" ? true : false
        new(datetime, point, route, FL, pointDistance, pointType, geoPointId,
        ratio, isVisible)
    end
end

function AllFtPointProfileList(items::Union{AbstractString, Missing})
    if items === missing
        return missing
    else
        return [AllFtPointProfile(item) for item in split(items)]
    end
end

struct AllFtAirspaceProfile
    entry_datetime::Union{DateTime, Missing}
    sector::AbstractString
    exit_datetime::Union{DateTime, Missing}
    fir::AbstractString
    entry_geoPointId::Union{Point_deg, AbstractString}
    exit_geoPointId::Union{Point_deg, AbstractString}
    # entry_geoPointId::Union{Point, AbstractString}
    # exit_geoPointId::Union{Point, AbstractString}
    entry_FL::Union{Int64, Missing}
    exit_FL::Union{Int64, Missing}
    entry_pointDistance::Union{Int64, Missing}
    exit_pointDistance::Union{Int64, Missing}
    function AllFtAirspaceProfile(x::AbstractString)
        items = split(x, ':')
        entry_datetime = items[1] == "" ? missing : format_datetime(items[1], yyyymmddhhmmss)
        sector = items[2]
        exit_datetime = items[3] == "" ? missing : format_datetime(items[3], yyyymmddhhmmss)
        fir = items[4]
        entry_geoPointId = occursin(r"\d{6}[NS]\d{7}[EW]", items[5]) ? latlon(items[5]) : items[5]
        exit_geoPointId = occursin(r"\d{6}[NS]\d{7}[EW]", items[6]) ? latlon(items[6]) : items[6]
        entry_FL = items[7] == "" ? missing : parse(Int64, items[7])
        exit_FL = items[8] == "" ? missing : parse(Int64, items[8])
        entry_pointDistance = items[9] == "" ? missing : parse(Int64, items[9])
        exit_pointDistance = items[10] == "" ? missing : parse(Int64, items[10])
        new(entry_datetime, sector, exit_datetime, fir, entry_geoPointId,
        exit_geoPointId, entry_FL, exit_FL, entry_pointDistance,
        exit_pointDistance)
    end
end

function AllFtAirspaceProfileList(items::Union{AbstractString, Missing})
    if items === missing
        return missing
    else
        return [AllFtAirspaceProfile(item) for item in split(items)]
    end
end

function reformat!(df)
    #Date conversion
    df[:,:TEMP] = format_datetime.(df[:,:aobt_5], yyyymmddhhmmss)
    select!(df, Not(:aobt_5))
    df[:,:aobt_5] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:iobt_7], yyyymmddhhmmss)
    select!(df, Not(:iobt_7))
    df[:,:iobt_7] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:cobt_16], yyyymmddhhmmss)
    select!(df, Not(:cobt_16))
    df[:,:cobt_16] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:eobt_17], yyyymmddhhmmss)
    select!(df, Not(:eobt_17))
    df[:,:eobt_17] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:lobt_18], yyyymmddhhmmss)
    select!(df, Not(:lobt_18))
    df[:,:lobt_18] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:samCtot_23], yyyymmddhhmmss)
    select!(df, Not(:samCtot_23))
    df[:,:samCtot_23] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:sipCtot_25], yyyymmddhhmmss)
    select!(df, Not(:sipCtot_25))
    df[:,:sipCtot_25] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:lastSentSlotMessage_45], yyyymmddhhmmss)
    select!(df, Not(:lastSentSlotMessage_45))
    df[:,:lastSentSlotMessage_45] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:lastSentProposalMessage_46], yyyymmddhhmmss)
    select!(df, Not(:lastSentProposalMessage_46))
    df[:,:lastSentProposalMessage_46] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:cdmEarlyTtot_61], yyyymmddhhmmss)
    select!(df, Not(:cdmEarlyTtot_61))
    df[:,:cdmEarlyTtot_61] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:cdmAoTtot_62], yyyymmddhhmmss)
    select!(df, Not(:cdmAoTtot_62))
    df[:,:cdmAoTtot_62] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:cdmAtcTtot_63], yyyymmddhhmmss)
    select!(df, Not(:cdmAtcTtot_63))
    df[:,:cdmAtcTtot_63] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:cdmSequencedTtot_64], yyyymmddhhmmss)
    select!(df, Not(:cdmSequencedTtot_64))
    df[:,:cdmSequencedTtot_64] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_time.(df[:,:cdmTaxiTime_65], mmmmss)
    select!(df, Not(:cdmTaxiTime_65))
    df[:,:cdmTaxiTime_65] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:cdmNoSlotBefore_70], yyyymmddhhmmss)
    select!(df, Not(:cdmNoSlotBefore_70))
    df[:,:cdmNoSlotBefore_70] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:scrObt_119], yyyymmddhhmmss)
    select!(df, Not(:scrObt_119))
    df[:,:scrObt_119] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:srrObt_128], yyyymmddhhmmss)
    select!(df, Not(:srrObt_128))
    df[:,:srrObt_128] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:surObt_137], yyyymmddhhmmss)
    select!(df, Not(:surObt_137))
    df[:,:surObt_137] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:dctObt_146], yyyymmddhhmmss)
    select!(df, Not(:dctObt_146))
    df[:,:dctObt_146] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:eventTime_173], yyyymmddhhmmss)
    select!(df, Not(:eventTime_173))
    df[:,:eventTime_173] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

end

end # module


text:
using JuMP, IJulia, Gurobi, DataFrames, CSV, Random, LinearAlgebra, Distributions

function generate_synthetic_data(β0, β1; SNR = 30)
    Random.seed!(99)
    X1 = rand(40,4)
    Y1 = β0 .+ X1 * β1
    noise = randn(40)
    k = sqrt(var(Y1) / (SNR*var(noise)))
    Y_with_noise = Y1 + k.*noise
    df = hcat(Y_with_noise,X1)
    return df
end

function create_data(n, p, contam_prct, sig_noise_ratio)
    Random.seed!(99)
    ## Create indices to contaminate X and y
    contam = randperm(n*(p+1))[1:Int(round(n*contam_prct))]
#     X = []
#     y = []

    ## Create X vars where the mean and sd of each successive column increases (no real reason)
    ## If the count of the nested for loop iterator == anything in contam indices, add from U(80,200)
    count = 0
    X = randn(n,p)
#     ## Reshape to create matrix
#     X = convert(Matrix, reshape(X, n, p))

    ## Create true 𝛽 values
    if p == 4
        𝛽 = [0.4, -0.5, 0.8, 0.5]
    elseif p == 5 && contam_prct == .12
        𝛽 = [0.3, -0.6, 1, 0.7, -0.5]
    elseif p == 5 && contam_prct == .24
        𝛽 = [0.3, -0.6, 1.3, 0.7, -0.7]
    elseif p == 6
        𝛽 = [0.3, -0.4, 0.6, 0.4, -0.6, 0.6]
    end

    y_pre_noise = X*𝛽
#     @show y
    noise = randn(n)
    e = sqrt(var(y_pre_noise)/(sig_noise_ratio*var(noise)))
    @show e
    y = y_pre_noise + e.*noise
#     @show y

    for i in 1:n
        count += 1
        if count in contam
            y[i] += rand(Uniform(6, 10), 1)[1]
        else
            y[i] = y[i]
        end
    end

    ## Perturb some of the X values now that Y has been created uncontaminated
    X = vec(reshape(X, n*p, 1))
    count2 = 0
    for i in 1:(n*p)
        count2 += 1
        if count2 in contam
            X[i] += rand(Uniform(4,6), 1)[1]
        end
    end
    X = convert(Matrix, reshape(X, n, p))

    ## Store outlier location in row, col form:
    contam_ind = sort(contam)
    cols = []
    rows = []
    for i in 1:length(contam_ind)
        append!(cols, ceil(Int, contam_ind[i]/n))
        if contam_ind[i] % n == 0
            append!(rows, n)
        else
            append!(rows, contam_ind[i] % n)
        end
    end
    outliers = hcat(cols, rows)

    return X, y, 𝛽, contam, outliers
end

function get_outliers(contam, n, p)
    ## Be able to return the indices of outliers
    contam_ind = sort(contam)
    cols = []
    rows = []
    for i in 1:length(contam_ind)
        append!(cols, ceil(Int, contam_ind[i]/n))
        if contam_ind[i] % n == 0
            append!(rows, n)
        else
            append!(rows, contam_ind[i] % n)
        end
    end

    for i in 1:length(contam_ind)
        if cols[i] <= p
            println("Outlier ", i, " in X", cols[i], ", row ", rows[i])
        else
            println("Outlier ", i, " in Y, row ", rows[i])
        end
    end
end

X1, y1, 𝛽1, contam1, outliers1 = create_data(25, 4, 0, 40);
X2, y2, 𝛽2, contam2, outliers2 = create_data(40, 5, 0.12, 40);
X3, y3, 𝛽3, contam3, outliers3 = create_data(40, 5, 0.24, 40);
X4, y4, 𝛽4, contam4, outliers4 = create_data(80, 6, 0.36, 40);

#get_outliers(contam1, 25, 4)
#get_outliers(contam2, 40, 5)
get_outliers(contam3, 40, 5)
#get_outliers(contam4, 80, 6)

# Test:
mean(y3 - X3*𝛽3)

sd_X1 = map(std, (X1[:,j] for j=1:size(X1,2)))
sd_X2 = map(std, (X2[:,j] for j=1:size(X2,2)))
sd_X3 = map(std, (X3[:,j] for j=1:size(X3,2)))
sd_X4 = map(std, (X4[:,j] for j=1:size(X4,2)))

𝛽1s = 𝛽1.*(sd_X1/std(y1))
𝛽2s = 𝛽2.*(sd_X2/std(y2))
𝛽3s = 𝛽3.*(sd_X3/std(y3))
𝛽4s = 𝛽4.*(sd_X4/std(y4));

## Create dataframes with X and y
df1 = hcat(y1, X1);
df2 = hcat(y2, X2);
df3 = hcat(y3, X3);
df4 = hcat(y4, X4);
#df5 = hcat(y5, X5);

CSV.write("Data/df1.csv", DataFrame(df1); writeheader=false);
CSV.write("Data/df2.csv", DataFrame(df2); writeheader=false);
CSV.write("Data/df3.csv", DataFrame(df3); writeheader=false);
CSV.write("Data/df4.csv", DataFrame(df4); writeheader=false);
#CSV.write("Data/df5.csv", DataFrame(df5); writeheader=false);

ext: jl

max_stars_repo_path: Data Creation Function Final.jl

max_issues_repo_path: Data Creation Function Final.jl

max_forks_repo_path: Data Creation Function Final.jl


text:
is_question(elem::HTMLElement) = any(is_question, elem.children)
is_question(elem::HTMLText) = occursin(r"^[0-9]+\. ", elem.text)

struct Question
    content
    responses::Vector{Response}
    label
    id
end

function question(contents, responses, label; detect_label=false)
    matched = 0
    for elem in contents
        for node in PreOrderDFS(elem)
            if node isa HTMLText
                m = match(r"([0-9]+)\. ", node.text)
                if detect_label && m !== nothing
                    label = m.captures[1]
                    matched += 1
                end
                node.text = replace(node.text, r"[0-9]+\. " => "")
            end
        end
    end

    detect_label && matched == 0 && @warn "Failed to detect label, using default ($label)"
    matched > 1 && @warn "Multiple lables detected, using last one ($label)"

    Question(contents, responses, label, string(uuid4()))
end

function Base.show(io::IO, q::Question)
    print(io, q.label, ": ")

    for c in q.content
        println(io, c)
    end

    print(io, "\n")

    for r in q.responses
        println(io, "  ", r)
    end
end


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

using Test
using LibCURL_jll

@testset "LibCURL_jll" begin
    v = unsafe_string(ccall((:curl_version, libcurl), Cstring, ()))
    @test startswith(v, "libcurl/")
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
module DevUtils

import Pkg
import LibGit2

using ..Helpers
import ..OscarCI.find_branch, ..OscarCI.github_repo_exists

export oscar_develop, oscar_update, oscar_branch, oscar_add_remotes

function fetch_project(repo::LibGit2.GitRepo)
   for remote in LibGit2.remotes(repo)
      LibGit2.fetch(repo; remote=remote)
   end
end

function add_remote(repo::LibGit2.GitRepo; fork=nothing)
   pkg = basename(LibGit2.path(repo))
   remote = isnothing(fork) ? "origin" : fork
   isnothing(LibGit2.lookup_remote(repo, remote)) || return

   @info "$pkg: adding remote '$remote'"
   url = pkg_url(pkg, fork=fork)
   if !github_repo_exists(pkg,fork)
      @warn "$pkg: no such repository at org '$fork'"
      return
   end
   LibGit2.GitRemote(repo, remote, url,
                              "+refs/heads/*:refs/remotes/$remote/*")
   LibGit2.set_remote_push_url(repo, remote, pkg_giturl(pkg, fork=fork))
   LibGit2.fetch(repo; remote=remote)
end

add_remote(pkgdir::AbstractString; kwargs...) =
   add_remote(LibGit2.GitRepo(pkgdir); kwargs...)

function create_tracking_branch(repo::LibGit2.GitRepo, branch::AbstractString; remote::AbstractString="origin")
   pkg = basename(LibGit2.path(repo))
   @info "$pkg: creating branch '$branch' tracking from '$remote'"
   hash = LibGit2.GitHash(repo, "refs/remotes/$remote/$branch")
   LibGit2.branch!(repo, branch, string(hash); track="$remote")
   # for some reason LibGit2.jl does not set the correct upstream branch...
   LibGit2.with(LibGit2.GitConfig, repo) do cfg
      LibGit2.set!(cfg, "branch.$branch.remote", remote)
   end
end

function create_branch(repo::LibGit2.GitRepo, branch::AbstractString; start="origin/master")
   pkg = basename(LibGit2.path(repo))
   @info "$pkg: creating branch '$branch' (based on '$start')"
   fetch_project(repo)
   commit = start == "HEAD" ?
      LibGit2.head(repo) :
      LibGit2.lookup_branch(repo,"$start")
   isnothing(commit) && (commit = LibGit2.lookup_branch(repo,"$start", true))
   if !isnothing(commit)
      LibGit2.branch!(repo, branch, string(LibGit2.GitHash(commit)))
   else
      @warn "$pkg: start point '$start' not found"
   end
end

create_branch(pkgdir::AbstractString, branch; kwargs...) = create_branch(LibGit2.GitRepo(pkgdir), branch; kwargs...)

function update_project(repo::LibGit2.GitRepo)
   pkg = basename(LibGit2.path(repo))
   head = LibGit2.head(repo)
   upstream = LibGit2.upstream(head)
   if !isnothing(upstream)
      @info "$pkg: updating from '$(LibGit2.shortname(upstream))'"
      upstream = LibGit2.GitAnnotated(repo, upstream)
      fetch_project(repo)
      LibGit2.ffmerge!(repo, upstream)
   else
      @warn "$pkg: branch '$(LibGit2.shortname(head))' has no upstream"
   end
end

update_project(pkgdir::AbstractString) = update_project(LibGit2.GitRepo(pkgdir))

function checkout_branch(pkg::AbstractString, devdir::AbstractString; branch::AbstractString, fork=nothing, url=nothing)
   @info "$pkg: checkout branch '$branch'"
   # libgit2 is quite confusing...
   repo = LibGit2.GitRepo(devdir)
   LibGit2.fetch(repo; remote="origin")
   remote = LibGit2.lookup_remote(repo, isnothing(fork) ? "origin" : fork)
   # add new remote for fork if necessary
   if !isnothing(fork)
      add_remote(repo; fork=fork)
   end
   remote = LibGit2.name(remote)
   # check for branch
   commit = LibGit2.lookup_branch(repo,"$branch")
   if commit != nothing
      # switch branch and update
      LibGit2.branch!(repo, branch)
      update_project(devdir)
   else
      # create new tracking branch
      create_tracking_branch(repo, branch; remote=remote)
   end
end

function merge_branch(pkg::AbstractString, devdir::AbstractString, branch::Any)
   branch = branch == true ? "origin/master" : "$branch"
   @info "$pkg: trying to merge '$branch'"
   repo = LibGit2.GitRepo(devdir)
   fetch_project(repo)
   branch = LibGit2.GitAnnotated(repo, branch)
   if !LibGit2.merge!(repo, [branch])
      @error "$pkg: unable to merge '$branch'"
   end
end

function clone_project(pkg::AbstractString, devdir::AbstractString; branch::AbstractString, fork=nothing, url=nothing)
   @info "$pkg: cloning to '$devdir' (branch '$branch'):"
   if !isnothing(fork)
      # we create the main repo as well as the fork as upstream
      repo = LibGit2.clone(pkg_url(pkg; full=true), devdir)
      add_remote(repo; fork=fork)
      create_tracking_branch(repo, branch; remote=fork)
   else
      repo = LibGit2.clone(url, devdir, branch=branch)
   end
   LibGit2.set_remote_push_url(repo,"origin", pkg_giturl(pkg))
end

function checkout_project(pkg::AbstractString, dir::AbstractString; branch::AbstractString="master", fork=nothing, merge=false)
   (url, branch, fork) = find_branch(pkg, branch; fork=fork)
   devdir = abspath(joinpath(dir, pkg))
   if isdir(devdir)
      checkout_branch(pkg, devdir; branch=branch, fork=fork, url=url)
   else
      clone_project(pkg, devdir; branch=branch, fork=fork, url=url)
   end
   if merge !=false && (branch != "master" || !isnothing(fork))
      merge_branch(pkg, devdir, merge)
   end
   return devdir
end

function pkg_array_to_dict(pkgs::Array{String})
   pkgdict = Dict{String,Any}()
   for pkg in pkgs
      pkgre = match(r"(\w+)(?:#(.*))?", pkg)
      push!(pkgdict, pkgre[1] => pkgre[2])
   end
   return pkgdict
end

"""
    oscar_develop(pkgs::Array{String}; <keyword arguments>)
    oscar_develop(pkgs::Dict{String,Any}; <keyword arguments>)

For each of the Oscar packages given in `pkgs`, create a new checkout in `dir`,
and try to create new tracking branches for _existing_ upstream branches `branch`
(or fall back to `master`).

_Note:_ To create new branches for a new feature please use [`oscar_branch`](@ref).

If `fork` is given the branch will be looked up in `https://github.com/fork/pkg.jl`
and a second remote is created automatically; in addition to `origin` which will
always point to the main repository.

The push-urls are set to `git@github.com:org-name/PackageName.jl` to facilitate
pushing via ssh (both for `origin` and the optional `fork`).

These package checkouts are then added (dev'd) to a new julia project in `dir/project`
which you can then use by running julia with:
```
julia --project=dir/project
```

# Arguments
- `pkgs`: list of packages to operate on, please omit the `.jl`.
- `branch::String="master"`: branch for checkout.
- `dir="oscar-dev"`: development subdirectory.
- `fork=nothing`: github organisation/user for branch lookup for all packages.
- `active_repo=nothing`: used in CI to reuse the existing checkout of that package,
  corresponding to the github variable `\$GITHUB_REPOSITORY`.
- `merge=false`: used for CI: if `true` try to merge the latest `origin/master` into
  each checked out project; or any other branch if given a String specifying a `LibGit2.GitAnnotated`.
  Please note that this will not create the merge-commit (similar to `--no-commit`).

Each package name can optionally contain a branchname and a fork url:
- `PackageName#somebranch` will checkout `somebranch` from the default upstream.
- `PackageName#https://github.com/myfork/PackageName.jl#otherbranch` will use the
  `myfork` user for this package.
The package name and branch can also be given as dictionary mapping `PackageName.jl`
to `[forkurl#]branchname`.

# Examples
```julia-repl
julia> oscar_develop(["Oscar","Polymake"]; branch="some_feature")
```
```julia-repl
julia> oscar_develop(["Oscar","Singular#more_rings"]; dir="dev_more_rings")
```
"""
function oscar_develop(pkgs::Dict{String,Any}; dir=default_dev_dir, branch::AbstractString="master", fork=nothing, active_repo=nothing, merge=false)
   mkpath(dir)
   active_pkg = pkg_from_repo(active_repo)
   if isnothing(fork) && !isnothing(active_repo) && fork_from_repo(active_repo) != pkg_org(active_pkg)
      fork = fork_from_repo(active_repo)
   end
   withenv("JULIA_PKG_DEVDIR"=>"$dir") do
      Pkg.activate(joinpath(dir,"project")) do
         @info "populating development directory '$dir':"
         if !isnothing(active_pkg) && !in(active_pkg, Helpers.non_jl_repo)
            # during CI we always need to dev the active checkout
            @info "  reusing current dir for $active_pkg"
            Pkg.develop(Pkg.PackageSpec(path="."))
         end
         for (pkg, pkgbranch) in pkgs
            if pkg === active_pkg
               continue
            else
               if pkgbranch == "release"
                  # make sure we have that package added explicitly
                  Pkg.add(pkg)
               else
                  isnothing(pkgbranch) && (pkgbranch=branch)
                  devdir = checkout_project(pkg, dir; branch=pkgbranch, fork=fork, merge=merge)
                  Pkg.develop(Pkg.PackageSpec(path=devdir))
               end
            end
         end
      end
   end
   @info "Please start julia with:\njulia --project=$(abspath(joinpath(dir,"project")))"
end

oscar_develop(pkgs::Array{String}; kwargs...) =
   oscar_develop(pkg_array_to_dict(pkgs); kwargs...)

"""
    oscar_update(; <keyword arguments>)

For each Oscar package in `dir` fetch all remotes and do a fast forward merge
with the currently tracked upstream branch.
Similiar to doing `git pull --ff-only` in each directory.

# Arguments
- `dir="oscar-dev"`: development subdirectory.
"""
function oscar_update(; dir=default_dev_dir)
   for pkgdir in joinpath.(dir, pkg_names(dir))
      if isdir(joinpath(pkgdir,".git"))
         update_project(pkgdir)
      end
   end
   # update project after switching branches
   Pkg.activate(joinpath(dir,"project")) do
      Pkg.update()
   end
end

"""
    oscar_add_remotes(fork::String; <keyword arguments>)
    oscar_add_remotes(pkgs::Array{String}, fork::String; <keyword arguments>)

For each Oscar package in `pkgs` (or existing in `dir` if `pkgs` is not given)
add a new git remote for `https://github.com/fork/PackageName.jl`.
The push-url is set to `git@github.com:fork/PackageName.jl` to facilitate
pushing via ssh.

# Arguments
- `pkgs::Array{String}`: list of packages to operate on, please omit the `.jl`.
- `fork::String`: github organisation/user for new remote
- `dir="oscar-dev"`: development subdirectory.
"""
function oscar_add_remotes(pkgs::Array{String}, fork::AbstractString; dir=default_dev_dir)
   for pkgdir in joinpath.(dir, pkgs)
      if isdir(joinpath(pkgdir, ".git"))
         add_remote(pkgdir; fork=fork)
      end
   end
end

oscar_add_remotes(fork::AbstractString; dir=default_dev_dir) =
   oscar_add_remotes(pkg_names(dir), fork; dir=dir)

"""
    oscar_branch(branch::AbstractString; <keyword arguments>)
    oscar_branch(pkgs::Array{String}, branch::AbstractString; <keyword arguments>)

For each of the Oscar packages given in `pkgs` create a new local branch named
`branch` with start point `start`, in the working copies under the current
development directory given via `dir`.
If no `pkgs` array is given it will run on all subdirectories of `dir`.

If the directory `dir` does not exist and `pkgs is given it will checkout all
`pkgs` first by calling [`oscar_develop`](@ref).

# Arguments
- `pkgs::Array{String}`: list of packages to operate on, please omit the `.jl`.
- `branch::String`: new branch name.
- `dir="oscar-dev"`: development subdirectory.
- `start="origin/master"`: start point for the branches, use `HEAD` for the current
                           head in each directory.
"""
function oscar_branch(pkgs::Array{String}, branch::AbstractString; dir=default_dev_dir, start="origin/master")
   # if path is not there, develop it first
   if !ispath(dir)
     oscar_develop(pkgs; dir=dir)
   end
   for pkgdir in joinpath.(dir, pkgs)
      if isdir(joinpath(pkgdir, ".git"))
         create_branch(pkgdir, branch; start=start)
      end
   end
   # update project after switching branches
   Pkg.activate(joinpath(dir,"project")) do
      Pkg.update()
   end
end

oscar_branch(branch::AbstractString; dir=default_dev_dir, start="origin/master") =
   oscar_branch(pkg_names(dir), branch; dir=dir, start=start)

end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module MarketData

using TimeSeries

export AAPL, BA, CAT, DELL, EBAY, F, GE,
       cl, op, ohlc, ohlcv, datetime1, datetime2, mdata,
       o, h, l, c, v

###### include ##################

include("const.jl")
include("helpers.jl")

end


text:
abstract type JuliaConTalkType end
struct LightningTalk <: JuliaConTalkType end
struct SponsorTalk <: JuliaConTalkType end
struct Talk <: JuliaConTalkType end
struct Keynote <: JuliaConTalkType end
struct BoF <: JuliaConTalkType end
struct Minisymposia <: JuliaConTalkType end
struct Workshop <: JuliaConTalkType end

struct JuliaConTalk{T<:JuliaConTalkType}
    start::String
    duration::String
    url::String
    title::String
    type::T
    speaker::Vector{String}
end

struct JuliaConTrack
    name::String
    talks::Vector{JuliaConTalk}
end

struct JuliaConDay
    number::Int
    date::Date
    start::DateTime
    stop::DateTime
    tracks::Vector{JuliaConTrack}
end

struct JuliaConSchedule
    start::Date
    stop::Date
    ndays::Int
    days::Vector{JuliaConDay}
end

function json2struct(conf)
    jcon_days = Vector{JuliaConDay}(undef, length(conf["days"]))
    k = 1
    for day in conf["days"]
        jcon_tracks = Vector{JuliaConTrack}(undef, length(day["rooms"]))
        j = 1
        for (track, talks) in day["rooms"]
            jcon_talks = Vector{JuliaConTalk}(undef, length(talks))
            for (i, talk) in enumerate(talks)
                jcon_talks[i] = JuliaConTalk(
                    talk["start"],
                    talk["duration"],
                    talk["url"],
                    talk["title"],
                    talktype_from_str(talk["type"]),
                    String[p["public_name"] for p in talk["persons"]],
                )
            end
            jcon_tracks[j] = JuliaConTrack(string(track), jcon_talks)
            j += 1
        end
        jcon_days[k] = JuliaConDay(
            k,
            Date(day["date"]),
            DateTime(day["day_start"][1:(end - 6)]),
            DateTime(day["day_end"][1:(end - 6)]),
            jcon_tracks,
        )
        k += 1
    end

    jcon = JuliaConSchedule(
        Date(conf["start"]), Date(conf["end"]), length(conf["days"]), jcon_days
    )
    return jcon
end

function talktype_from_str(str)
    if str == "Talk"
        return Talk()
    elseif str == "Lightning Talk"
        return LightningTalk()
    elseif str == "Sponsor Talk"
        return SponsorTalk()
    elseif str == "Keynote"
        return Keynote()
    elseif str == "Birds of Feather"
        return BoF()
    elseif str == "Minisymposia"
        return Minisymposia()
    elseif contains(str, "Workshop")
        return Workshop()
    else
        error("Unknown JuliaCon talk type \"$str\".")
    end
end

function show(io::IO, x::JuliaConSchedule)
    return print(io, "JuliaCon ", Dates.Year(x.start).value, " Schedule")
end
show(io::IO, ::MIME"text/plain", x::JuliaConSchedule) = show(io, x)

function show(io::IO, x::JuliaConDay)
    return print(io, "JuliaCon ", Dates.Year(x.start).value, ": Day ", x.number)
end
show(io::IO, ::MIME"text/plain", x::JuliaConDay) = show(io, x)

show(io::IO, x::JuliaConTrack) = print(io, "JuliaCon Track: ", x.name)
function show(io::IO, ::MIME"text/plain", x::JuliaConTrack)
    show(io, x)
    if length(x.talks) > 1
        print(io, " (", length(x.talks), " talks)")
    else
        print(io, " (", length(x.talks), " talk)")
    end
end

show(io::IO, x::JuliaConTalk{T}) where {T<:JuliaConTalkType} = print(io, T)
function show(io::IO, ::MIME"text/plain", x::JuliaConTalk{T}) where {T<:JuliaConTalkType}
    show(io, x)
    println(io)
    println(io, "├ Title: ", x.title)
    println(io, "├ Time: ", x.start)
    return print(io, "└ URL: ", x.url)
end

show(io::IO, x::LightningTalk) = print(io, "Lightning Talk")
show(io::IO, x::SponsorTalk) = print(io, "Sponsor Talk")
show(io::IO, x::Talk) = print(io, "Talk")
show(io::IO, x::Keynote) = print(io, "Keynote")
show(io::IO, x::BoF) = print(io, "Birds of Feather")
show(io::IO, x::Minisymposia) = print(io, "Minisymposia")
show(io::IO, x::Workshop) = print(io, "Workshop")

abbrev(::Type{LightningTalk}) = "L"
abbrev(::Type{SponsorTalk}) = "S"
abbrev(::Type{Talk}) = "T"
abbrev(::Type{Keynote}) = "K"
abbrev(::Type{BoF}) = "BoF"
abbrev(::Type{Minisymposia}) = "M"
abbrev(::Type{Workshop}) = "W"

abbrev(x::JuliaConTalkType) = abbrev(typeof(x))


text:
const NodeInt = UInt32

struct SGNode
   gene::NodeInt
   node::NodeInt
end

const SGNodeSet = Vector{SGNode}

struct Edges{K}
   left::Vector{SGNodeSet}
   right::Vector{SGNodeSet}
end

Base.:<( a::SGNode, b::SGNode ) = <( a.gene, b.gene )
Base.:>( a::SGNode, b::SGNode ) = >( a.gene, b.gene )
Base.:<=( a::SGNode, b::SGNode ) = <=( a.gene, b.gene )
Base.:>=( a::SGNode, b::SGNode ) = >=( a.gene, b.gene )

sortlt( a::SGNode, b::SGNode ) = a.gene == b.gene ? <( a.node, b.node ) : <( a.gene, b.gene )

Base.convert( ::Type{Edges{K}}, graphs::Vector{SpliceGraph} ) where {K <: Integer} = build_edges( graphs, K )

# Is the EdgeType a connecting edge between two nodes?
# Check a donor splice site by:   is_edge( edge, true )
# and an acceptor splice site by: is_edge( edge, false )
# Return: Bool
function is_edge( edge::EdgeType, left::Bool )
   eint = convert(UInt8, edge)
   eint == 0x04 && return true
   left || (eint ⊻= 0b011)
   (eint == 0x02 || eint == 0x05) && return true
   false
end

Base.intersect( arrA::Vector{SGNode}, arrB::Vector{SGNode} ) = intersect_sorted( arrA, arrB )

# Intersect two sorted arrays of SGNode by gene entry
# Return the list of intersected arrB entries
function intersect_sorted( arrA::Vector{T}, arrB::Vector{T} ) where T
   res = Vector{T}()
   i,j = 1,1
   @inbounds while i <= length(arrA) && j <= length(arrB)
      if arrA[i] > arrB[j]
         j += 1
      elseif arrA[i] < arrB[j]
         i += 1
      else
         push!( res, arrB[j] )
         i += 1
         j += 1
      end
   end
   res
end

function add_kmer_edge!( kmers::Vector{SGNodeSet},
                         seq::S, l, r, left::Bool,
                         entry::SGNode ) where S <: SGSequence
   (l <= 0 || r > length(seq)) && return(zero(UInt64))
   s = copy(seq[l:r])
   ksize = r-l+1
   ind = 0 #default
   try
      curkmer = sgkmer( s )
      ind = kmer_index(curkmer)
   catch
      abstr = String(s)
      ismatch( r"S|N", abstr ) && return(zero(UInt64))
      sub = replace( abstr, r"D|R", "" )
      curl,curr = l-1,r+1
      while length(sub) < ksize && curl >= 1 && curr <= length(seq)
         sub = left ? String(copy(seq[curl:curl])) * sub :
                      sub * String(copy(seq[curr:curr]))
         sub = replace( sub, r"D|R", "" )
         curr += 1
         curl -= 1
      end
      try
        curkmer = sgkmer( s )
        ind = kmer_index(curkmer)
      catch
      end
   end
   if ind > 0
      if !isassigned(kmers, ind)
         kmers[ind] = SGNodeSet()
      end
      push!(kmers[ind], entry)
      kmers[ind] = sort( kmers[ind], lt=sortlt ) |> unique
   end
   UInt64(max(0, ind-1))
end

function build_edges( graphs::Vector{SpliceGraph}, k::Integer )
   @assert(1 <= k <= 32, "$k is out of bounds for Edges{K} set!")
   left  = Vector{SGNodeSet}(4^k)
   right = Vector{SGNodeSet}(4^k)

   for (i,g) in enumerate(graphs)
      for (j,n) in enumerate(g.nodeoffset)

         if is_edge( g.edgetype[j], true ) #&& isvalid( g.seq, (n-k-2):(n-3) ) # left edge
            lkmer = add_kmer_edge!( left, g.seq, n-k, n-1, true,  SGNode(i,j) )
            g.edgeleft[j] = convert(SGKmer{k}, lkmer)
         end
         if is_edge( g.edgetype[j], false ) #&& isvalid( g.seq, n:(n+k-1) )# right edge
            rkmer = add_kmer_edge!( right, g.seq, n,  n+k-1, false, SGNode(i,j) )
            g.edgeright[j] = convert(SGKmer{k}, rkmer)
         end

      end
   end
   Edges{k}(left, right)
end


text:
using LinearAlgebra
@testset "Newton–Schulz iteration" begin
    A = [2.7 4; 5.3 6.6] / 7

    k = 4
    (graph, cref) = graph_newton_schulz(k)
    (graph_2, cref_2) = graph_newton_schulz_degopt(k)
    Vk = A
    for i = 1:k
        Vkk = Vk * (2 * I - A * Vk)
        Vk = Vkk
    end
    @test eval_graph(graph, A) ≈ Vk
    @test eval_graph(graph_2, A) ≈ Vk

    (graph, cref) = graph_newton_schulz(15)
    @test eval_graph(graph, A) ≈ inv(A)

    (graph_2, cref_2) = graph_newton_schulz_degopt(15)
    @test eval_graph(graph, A) ≈ eval_graph(graph_2, A)
end


text:
using MonteCarloMarkovKernels
using Random
import Statistics: mean, var
using Test

seed!(12345)

function makeAR1Kernel(c::Float64, φ::Float64, σ::Float64)
  function P(x::Float64)
    return c + φ * x + σ * randn()
  end
  return P
end

const c = 2.32
const φ = 0.98
const σ = 2.3
P = makeAR1Kernel(c, φ, σ)

chain = simulateChain(P, 0.0, 1024*1024)

@test mean(chain) ≈ c/(1-φ) rtol=0.1
@test var(chain) ≈ σ^2/(1-φ^2) rtol=0.1

@test estimateBM(chain) ≈ σ^2/(1-φ^2) * (1+φ)/(1-φ) rtol=0.5
@test estimateSV(chain, :ModifiedBartlett) ≈
  σ^2/(1-φ^2) * (1+φ)/(1-φ) rtol=0.5


text:
module GeoRDD
    using Statistics
    import Statistics: mean
    using GaussianProcesses
    using PDMats
    using Optim
    import GeoInterface
    using GeoInterface: coordinates, xcoord, ycoord
    using LibGEOS: nearestPoints, interpolate, distance
    using LibGEOS: MultiPolygon, envelope
    import LibGEOS
    import LibGEOS: interpolate
    import Combinatorics
    using LinearAlgebra
    import Base: getindex, keys, values, start, done, next, iterate, eltype, length, size, convert
    using Distributions: Normal, MultivariateNormal, ccdf, cdf
    using PDMats: AbstractPDMat
    import GaussianProcesses: update_mll!, update_mll_and_dmll!,
                              get_params, set_params!, num_params,
                              optimize!, GPE
    using GaussianProcesses: grad_stack, grad_stack!, grad_slice!, get_ααinvcKI!,
                             Mean, Kernel, KernelData, LinIso, MeanZero,
                             cov!, cov, cov_ij, dmll_kern!,
                             mat, cholfactors, wrap_cK, make_posdef!
    import StatsModels
    using StatsModels: Formula
    import DataFrames
    using DataFrames: DataFrame, AbstractDataFrame
    import NLopt

    # utilities
    include("gp_utils.jl")
    include("geometry.jl")
    # Types for multiple realizations of GPs
    include("regiondata.jl")
    include("GPrealisations.jl")
    include("GPrealisationsCovars.jl")
    include("region_gp_fit.jl")
    # Cliff Face (treatment effect estimation)
    include("cliff_face.jl")
    # Local average treatment effect estimation
    include("border_projection.jl")
    include("average_treatment_effect.jl")
    include("weight_at_units.jl")
    # Hypothesis testing
    include("hypotest.jl")
    include("hypotest_chi2.jl")
    include("hypotest_invvar.jl")
    include("hypotest_mll.jl")
    include("placebo_geometry.jl")
    include("power.jl")
end


text:
# Test items in FWF.jl file

# Ensure options objects get created
@testset "Options Tests" begin
    x = FWF.Options()
    @test x.usemissings == true
    @test x.trimstrings == true
    @test x.errorlevel == :parse
    @test x.unitbytes == true
    @test x.skip == 0
    @test x.missingvals == Set{String}()
    @test x.dateformats == Dict{Int, DateFormat}()
    @test x.columnrange == Vector{UnitRange{Int}}()

    m = Set{String}(["***"])
    d = Dict{Int, DateFormat}(1 => DateFormat("mmddyy"))
    x = [1:10, 2:20]
    tmp = FWF.Options(usemissings=false, trimstrings=false, errorlevel=:skip,
                      unitbytes=false, skip=10, missingvals=m, dateformats=d, columnrange=x)
    @test tmp.usemissings == false
    @test tmp.trimstrings == false
    @test tmp.errorlevel == :skip
    @test tmp.unitbytes == false
    @test tmp.skip == 10
    @test tmp.missingvals == m
    @test tmp.dateformats == d
    @test tmp.columnrange == x

    @test_throws ArgumentError FWF.Options(usemissings=false, trimstrings=false,
                                           errorlevel=:x, unitbytes=false, skip=10,
                                           missingvals=m, dateformats=d, columnrange=x)
end


text:
@testset "warlpiri-addition" begin

  #commutativity test.
  for a in WP, b in WP, c in WP, d in WP
      @test ((a → b) + (c → d)) == ((c → d) + (a → b))
  end

  const warlpiri_sum_matrix = [
   #W0_0 W0_1 W0_2 W0_3 W0_4 W0_5 W0_6 W0_7 W0_8   W1_1 W1_2 W1_3 W1_4 W1_5 W1_6 W1_7 W1_8   W2_2 W2_3 W2_4 W2_5 W2_6 W2_7 W2_8   W3_3 W3_4 W3_5 W3_6 W3_7 W3_8   W4_4 W4_5 W4_6 W4_7 W4_8   W5_5 W5_6 W5_7 W5_8   W6_6 W6_7 W6_8   W7_7 W7_8   W8_8
    W0_0 W0_1 W0_2 W0_3 W0_4 W0_5 W0_6 W0_7 W0_8   W1_1 W1_2 W1_3 W1_4 W1_5 W1_6 W1_7 W1_8   W2_2 W2_3 W2_4 W2_5 W2_6 W2_7 W2_8   W3_3 W3_4 W3_5 W3_6 W3_7 W3_8   W4_4 W4_5 W4_6 W4_7 W4_8   W5_5 W5_6 W5_7 W5_8   W6_6 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_0
    W0_1 W0_3 W0_3 W0_5 W0_5 W0_7 W0_7 W0_7 W0_8   W1_3 W1_3 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W2_3 W2_5 W2_5 W2_7 W2_7 W2_7 W2_8   W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W4_5 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_1
    W0_2 W0_3 W0_4 W0_5 W0_5 W0_7 W0_7 W0_7 W0_8   W1_3 W1_4 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W2_4 W2_5 W2_5 W2_7 W2_7 W2_7 W2_8   W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W4_5 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_2
    W0_3 W0_5 W0_5 W0_5 W0_5 W0_7 W0_7 W0_7 W0_8   W1_5 W1_5 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W2_5 W2_5 W2_5 W2_7 W2_7 W2_7 W2_8   W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W4_5 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_3
    W0_4 W0_5 W0_5 W0_5 W0_6 W0_7 W0_7 W0_7 W0_8   W1_5 W1_5 W1_5 W1_6 W1_7 W1_7 W1_7 W1_8   W2_5 W2_5 W2_6 W2_7 W2_7 W2_7 W2_8   W3_5 W3_6 W3_7 W3_7 W3_7 W3_8   W4_6 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_4
    W0_5 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_5
    W0_6 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_6
    W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_8   W7_7 W7_8   W8_8 # + W0_7
    W0_8 W0_8 W0_8 W0_8 W0_8 W0_8 W0_8 W0_8 WALL   W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W2_8 W2_8 W2_8 W2_8 W2_8 W2_8 WALL   W3_8 W3_8 W3_8 W3_8 W3_8 WALL   W4_8 W4_8 W4_8 W4_8 WALL   W5_8 W5_8 W5_8 WALL   W6_8 W6_8 WALL   W7_8 WALL   WALL # + W0_8

    W1_1 W1_3 W1_3 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W1_3 W1_3 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W3_3 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W5_5 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_1
    W1_2 W1_3 W1_4 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W1_3 W1_4 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W3_4 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W5_5 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_2
    W1_3 W1_5 W1_5 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W1_5 W1_5 W1_5 W1_5 W1_7 W1_7 W1_7 W1_8   W3_5 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W5_5 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_3
    W1_4 W1_5 W1_5 W1_5 W1_6 W1_7 W1_7 W1_7 W1_8   W1_5 W1_5 W1_5 W1_6 W1_7 W1_7 W1_7 W1_8   W3_5 W3_5 W3_6 W3_7 W3_7 W3_7 W3_8   W3_5 W3_6 W3_7 W3_7 W3_7 W3_8   W5_6 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_4
    W1_5 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_5
    W1_6 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_6
    W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W1_7
    W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 WALL   W3_8 W3_8 W3_8 W3_8 W3_8 WALL   W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W1_8

    W2_2 W2_3 W2_4 W2_5 W2_5 W2_7 W2_7 W2_7 W2_8   W3_3 W3_4 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W4_4 W4_5 W4_5 W4_7 W4_7 W4_7 W4_8   W5_5 W5_5 W5_7 W5_7 W5_7 W5_8   W5_5 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W2_2
    W2_3 W2_5 W2_5 W2_5 W2_5 W2_7 W2_7 W2_7 W2_8   W3_5 W3_5 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W4_5 W4_5 W4_5 W4_7 W4_7 W4_7 W4_8   W5_5 W5_5 W5_7 W5_7 W5_7 W5_8   W5_5 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W2_3
    W2_4 W2_5 W2_5 W2_5 W2_6 W2_7 W2_7 W2_7 W2_8   W3_5 W3_5 W3_5 W3_6 W3_7 W3_7 W3_7 W3_8   W4_5 W4_5 W4_6 W4_7 W4_7 W4_7 W4_8   W5_5 W5_6 W5_7 W5_7 W5_7 W5_8   W5_6 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W2_4
    W2_5 W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W2_5
    W2_6 W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W2_6
    W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_7 W2_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W2_7
    W2_8 W2_8 W2_8 W2_8 W2_8 W2_8 W2_8 W2_8 WALL   W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 WALL   W4_8 W4_8 W4_8 W4_8 W4_8 W4_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W2_8

    W3_3 W3_5 W3_5 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W3_5 W3_5 W3_5 W3_5 W3_7 W3_7 W3_7 W3_8   W5_5 W5_5 W5_5 W5_7 W5_7 W5_7 W5_8   W5_5 W5_5 W5_7 W5_7 W5_7 W5_8   W5_5 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W3_3
    W3_4 W3_5 W3_5 W3_5 W3_6 W3_7 W3_7 W3_7 W3_8   W3_5 W3_5 W3_5 W3_6 W3_7 W3_7 W3_7 W3_8   W5_5 W5_5 W5_6 W5_7 W5_7 W5_7 W5_8   W5_5 W5_6 W5_7 W5_7 W5_7 W5_8   W5_6 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W3_4
    W3_5 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W3_5
    W3_6 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W3_6
    W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_7 W3_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W3_7
    W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 WALL   W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 W3_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W3_8

    W4_4 W4_5 W4_5 W4_5 W4_6 W4_7 W4_7 W4_7 W4_8   W5_5 W5_5 W5_5 W5_6 W5_7 W5_7 W5_7 W5_8   W5_5 W5_5 W5_6 W5_7 W5_7 W5_7 W5_8   W5_5 W5_6 W5_7 W5_7 W5_7 W5_8   W6_6 W6_7 W6_7 W6_7 W6_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W4_4
    W4_5 W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_7 W6_7 W6_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W4_5
    W4_6 W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_7 W6_7 W6_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W4_6
    W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_7 W4_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W6_7 W6_7 W6_7 W6_7 W6_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W4_7
    W4_8 W4_8 W4_8 W4_8 W4_8 W4_8 W4_8 W4_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W6_8 W6_8 W6_8 W6_8 WALL   W7_8 W7_8 W7_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W4_8

    W5_5 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W5_5
    W5_6 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W5_6
    W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W5_7 W5_7 W5_7 W5_7 W5_7 W5_8   W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W5_7
    W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W5_8 W5_8 W5_8 W5_8 W5_8 WALL   W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W5_8

    W6_6 W6_7 W6_7 W6_7 W6_7 W6_7 W6_7 W6_7 W6_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W6_6
    W6_7 W6_7 W6_7 W6_7 W6_7 W6_7 W6_7 W6_7 W6_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W6_7
    W6_8 W6_8 W6_8 W6_8 W6_8 W6_8 W6_8 W6_8 WALL   W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W6_8

    W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_7 W7_8   W7_7 W7_7 W7_8   W7_7 W7_8   W8_8 # + W7_7
    W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 W7_8 WALL   W7_8 W7_8 W7_8 WALL   W7_8 W7_8 WALL   W7_8 WALL   WALL # + W7_8

    W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 WALL   W8_8 W8_8 WALL   W8_8 WALL   WALL # + W8_8
    ]

    warlpiri_add_test_matrix = [a + b for a in pos_warlpiri, b in pos_warlpiri]
    warlpiri_label_matrix = [(a, b) for a in pos_warlpiri, b in pos_warlpiri]

    for (expected, evaluated, l) in zip(warlpiri_sum_matrix, warlpiri_add_test_matrix, warlpiri_label_matrix)
        @test (expected, l) == (evaluated, l)
    end

    #test some edge cases for warlpiri addition
    #Inf → 1 + 2 == allreals
    @test (WP1000 → WP0100) + W6_6 == WALL
    @test (WP0100 → WP1100) + (WP1100 → WP0100) == WALL

    for x in pos_warlpiri
        @test x + WALL == WALL
        @test x + WNUL == WNUL
    end
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
include("./imports.jl")

@testset ExtendedTestSet "apply!(buf, ::Map, ...)" begin
    newitem() = ArrayItem(randn(5, 5))
    tfm = MapElem(x -> x + 1)

    buf = apply(tfm, newitem())
    buf_copy = deepcopy(buf)

    @test buf.data ≈ buf_copy.data
    @test_nowarn apply!(buf, tfm, newitem())
    @test !(buf.data ≈ buf_copy.data)
end

@testset ExtendedTestSet "copyitemdata!" begin
    newitem() = ArrayItem(randn(5, 5))
    @testset ExtendedTestSet "" begin
        a = newitem()
        b = newitem()
        @test !(a.data ≈ b.data)
        DataAugmentation.copyitemdata!(a, b)
        @test a.data ≈ b.data
    end

    @testset ExtendedTestSet "" begin
        a = newitem()
        b = newitem()
        @test !(a.data ≈ b.data)
        DataAugmentation.copyitemdata!([a], [b])
        @test a.data ≈ b.data
    end

    @testset ExtendedTestSet "" begin
        a = newitem()
        b = newitem()
        @test !(a.data ≈ b.data)
        DataAugmentation.copyitemdata!((a,), (b,))
        @test a.data ≈ b.data
    end
end

@testset ExtendedTestSet "`Buffered`" begin
    newitem() = ArrayItem(randn(5, 5))
    T = typeof(newitem())
    # buffer should be created
    tb = Buffered(MapElem(x -> x + one(typeof(x))))
    @test !haskey(tb.buffers, T)
    @test_nowarn apply(tb, newitem())
    @test haskey(tb.buffers, T)

    # buffer should change
    buf = deepcopy(tb.buffers[T])
    @test_nowarn apply(tb, newitem())
    @test !(buf.data ≈ tb.buffers[T].data)

    # inplace version
    buf2 = deepcopy(buf)
    @test_nowarn apply!(buf, tb, newitem())
    @test !(buf2.data ≈ tb.buffers[T].data)
    testapply(Buffered(MapElem(x -> x + one(typeof(x)))), ArrayItem)
end

@testset ExtendedTestSet "`BufferedThreadsafe`" begin
    newitem() = ArrayItem(randn(5, 5))
    T = typeof(newitem())
    # buffer should be created
    tbt = BufferedThreadsafe(MapElem(x -> x + 1))
    tb = tbt.buffereds[1]
    @test !haskey(tb.buffers, T)
    @test_nowarn apply(tb, newitem())
    @test haskey(tb.buffers, T)

    # buffer should change
    buf = deepcopy(tb.buffers[T])
    @test_nowarn apply(tb, newitem())
    @test !(buf.data ≈ tb.buffers[T].data)

    # inplace version
    buf2 = deepcopy(buf)
    @test_nowarn apply!(buf, tb, newitem())
    @test !(buf2.data ≈ tb.buffers[T].data)

    testapply(BufferedThreadsafe(MapElem(x -> x + one(typeof(x)))), ArrayItem)
end

struct PlusRand <: Transform
end
DataAugmentation.getrandstate(::PlusRand) = rand()
function DataAugmentation.apply(tfm::PlusRand, item::DataAugmentation.AbstractItem; randstate = getrandstate(tfm))
    return DataAugmentation.setdata(item, map(x -> x + randstate, itemdata(item)))
end
function DataAugmentation.apply!(buf, tfm::PlusRand, item::DataAugmentation.AbstractItem; randstate = getrandstate(tfm))
    map!(x -> x + randstate, itemdata(buf), itemdata(item))
    return buf
end

@testset ExtendedTestSet "" begin
    item1 = ArrayItem(ones(10))
    item2 = ArrayItem(ones(10))
    tfm = PlusRand()
    @test_nowarn apply(tfm, item1)
    titem1, titem2 = apply(tfm, (item1, item2))
    a1, a2 = itemdata.((titem1, titem2))
    b1, b2 = copy.((a1, a2))
    @test a1 == a2
    apply!((titem1, titem2), tfm, (item1, item2))
    @test a1 == a2
    @test a1 != b1
end


text:
using Base.Test
import WaterData

module test_eos_resources
import WaterData

struct NoTempInfo <: WaterData.EOS; end
notempinfo = NoTempInfo()

end

@testset "EOS base tests" begin
    res = test_eos_resources
    @testset "Tagging whether an EOS is temperature dependent" begin
        @test_throws WaterData.NotImplementedError WaterData.istempdependent(res.notempinfo)
    end
end


text:
md"""
# Mutation Operations

Mutation is the way one candidate is transformed to a slightly different candidate. NaiveGAflux supports doing this while
preserving parameters and alignment between layers, thus reducing the impact of mutating an already trained candidate.

The following basic mutation operations are currently supported:
1. Change the output size of vertices using [`NoutMutation`](@ref).
2. Remove vertices using [`RemoveVertexMutation`](@ref).
3. Add vertices using [`AddVertexMutation`](@ref).
4. Remove edges between vertices using [`RemoveEdgeMutation`](@ref).
5. Add edges between vertices using [`AddEdgeMutation`](@ref).
6. Mutation of kernel size for conv layers using [`KernelSizeMutation`](@ref).
7. Change of activation function using [`ActivationFunctionMutation`](@ref).
8. Change the type of optimizer using [`OptimizerMutation`](@ref).
9. Add an optimizer using [`AddOptimizerMutation`](@ref).

Mutation operations are exported as structs rather than functions since they are designed to be composed with more generic utilities. Here are a few examples:
"""

@testset "Mutation examples" begin #src
Random.seed!(NaiveGAflux.rng_default, 3) #src

# Start with a simple model to mutate.
invertex = denseinputvertex("in", 3)
layer1 = fluxvertex(Dense(nout(invertex), 4), invertex)
layer2 = fluxvertex(Dense(nout(layer1), 5), layer1)
graph = CompGraph(invertex, layer2)

# Create an [`NoutMutation`](@ref) to mutate it.
mutation = NoutMutation(-0.5, 0.5)

@test nout(layer2) == 5

mutation(layer2)

@test nout(layer2) == 7

# [`VertexMutation`](@ref) applies the wrapped mutation to all vertices in a `CompGraph`
mutation = VertexMutation(mutation)

@test nout.(vertices(graph)) == [3,4,7]

mutation(graph)

@test nout.(vertices(graph)) == [3,6,10]

# Input vertex is never mutated, but the other two changed.
# Use the [`MutationShield`](@ref) trait to protect otherwise mutable vertices from mutation.
outlayer = fluxvertex(Dense(nout(layer2), 10), layer2, traitfun = MutationShield)
graph = CompGraph(invertex, outlayer)

mutation(graph)

@test nout.(vertices(graph)) == [3,9,6,10]

# In most cases it makes sense to mutate with a certain probability.
mutation = VertexMutation(MutationProbability(NoutMutation(-0.5, 0.5), 0.5))

mutation(graph)

@test nout.(vertices(graph)) == [3,9,3,10]

# Or just chose to either mutate the whole graph or don't do anything.
mutation = MutationProbability(VertexMutation(NoutMutation(-0.5, 0.5)), 0.98)

mutation(graph)

@test nout.(vertices(graph)) == [3,10,2,10]
@test size(graph(ones(3,1))) == (10, 1)

# Mutation can also be conditioned:
mutation = VertexMutation(MutationFilter(v -> nout(v) < 10, RemoveVertexMutation()))

mutation(graph)

@test nout.(vertices(graph)) == [3,10,10]

# When adding vertices it is probably a good idea to try to initialize them as identity mappings.
addmut = AddVertexMutation(VertexSpace(DenseSpace(5, identity)), IdentityWeightInit())

# Chaining mutations is also useful:
noutmut = NoutMutation(-0.8, 0.8)
mutation = VertexMutation(MutationChain(addmut, noutmut))

mutation(graph)

@test nout.(vertices(graph)) == [3,6,10,10]
end #src


text:
"""
    DynamicStall{TI,TV} <: AbstractModel

Dynamic stall model constructed by coupling an invisicid flow model and a
viscous flow model
"""
struct DynamicStall{TI,TV} <: AbstractModel
    inviscid::TI # inviscid (attached flow) model
    viscous::TV # viscous flow model
end

# Viscous Flow Models
struct Onera{TF}
    a0::TF = 0.3
    a2::TF = 0.2
    e2::TF = -2.86
    r0::TF = 0.2
    r2::TF = 0.2
end

# --- Constructors --- #

"""
    DynamicStall(inviscid, viscous)

Initialize an object of type [`DynamicStall`](@ref)
"""
DynamicStall()

# --- Traits --- #

function number_of_states(::Type{DynamicStall{TI,TV}}) where {TI,TV}
    return number_of_states(TI) + number_of_states(TV)
end
function number_of_inputs(::Type{DynamicStall{TI,TV}}) where {TI,TV}
    return number_of_inputs(TI) + number_of_inputs(TV)
end
function number_of_parameters(::Type{DynamicStall{TI,TV}}) where {TI,TV}
    return number_of_parameters(TI) + number_of_parameters(TV)
end
inplaceness(::Type{<:DynamicStall}) = OutOfPlace()
mass_matrix_type(::Type{<:DynamicStall}) = Identity()
state_jacobian_type(::Type{<:DynamicStall}) = Linear()
input_jacobian_type(::Type{<:DynamicStall}) = Nonlinear()

# --- Methods --- #

function get_rates(model::Wagner, λ, d, p, t) where {N,TF,SV,SA}
    # extract states
    λ1, λ2 = λ
    # extract inputs
    u, v, ω = d
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    C1 = model.C1
    C2 = model.C2
    ε1 = model.eps1
    ε2 = model.eps2
    # calculate rates
    return wagner_rates(a, b, α0, C1, C2, ε1, ε2, u, v, ω, λ1, λ2)
end

# --- Performance Overloads --- #

function get_state_jacobian(model::Wagner, λ, d, p, t)
    # extract inputs
    u, v, ω = d
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    ε1 = model.eps1
    ε2 = model.eps2
    # jacobian with respect to aerodynamic states
    return wagner_state_jacobian(a, b, ε1, ε2, u)
end

function get_input_jacobian(model::Wagner, λ, d, p, t)
    # extract states
    λ1, λ2 = λ
    # extract inputs
    u, v, ω = d
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    C1 = model.C1
    C2 = model.C2
    ε1 = model.eps1
    ε2 = model.eps2
    # return jacobian
    return wagner_input_jacobian(a, b, α0, C1, C2, ε1, ε2, u, v, ω, λ1, λ2)
end

# --- Unit Testing Methods --- #

get_lhs(::Wagner, dλ, λ, d, p, t) = dλ

# --- Convenience Methods --- #

set_states!(y, model::Wagner; lambda) = y .= lambda

function set_inputs!(y, model::Wagner; u, v, omega)

    y[1] = u
    y[2] = v
    y[3] = omega

    return y
end

function set_parameters!(p, model::Wagner; a, b, a0, alpha0)

    p[1] = a
    p[2] = b
    p[3] = a0
    p[4] = alpha0

    return p
end

function separate_states(model::Wagner, x)

    return (lambda = x,)
end

function separate_inputs(model::Wagner, y)

    return (u=y[1], v=y[2], omega=y[3])
end

function separate_parameters(model::Wagner, p)

    return (a=p[1], b=p[2], a0=p[3], alpha0=p[4])
end

# --- Internal Methods --- #


text:
using Documenter, Convex, Literate, Pkg

# Needed to run GR headless on Travis
previous_GKSwstype = get(ENV, "GKSwstype", "")
ENV["GKSwstype"] = "100"

@info "Building examples..."

filename(str) = first(splitext(last(splitdir(str))))
filename_to_name(str) = uppercasefirst(replace(replace(filename(str), "-" => " "), "_" => " "))

fix_math_md(content) = replace(content, r"\$\$(.*?)\$\$"s => s"```math\1```")

literate_path = joinpath(@__DIR__(), "examples_literate")
build_path =  joinpath(@__DIR__, "src", "examples")
rm(build_path; force=true, recursive=true)
isdir(build_path) || mkdir(build_path)

@info "[Examples] Preparing notebooks..."

notebooks_path = joinpath(@__DIR__, "notebooks")
rm(notebooks_path, recursive=true, force=true)
mkdir(notebooks_path)

for dir in readdir(literate_path)
    dir_path = joinpath(literate_path, dir)
    isdir(dir_path) || continue
    @info "Processing directory $dir"
    notebook_dir = joinpath(notebooks_path, dir)
    isdir(notebook_dir) || mkdir(notebook_dir)
    for file in readdir(dir_path)
        file_path = joinpath(dir_path, file)
            out_path = joinpath(notebooks_path, dir, file)
        if endswith(file, ".jl")
            Literate.notebook(file_path, notebook_dir, execute=false)
        else
            cp(file_path, out_path)
        end
    end
end

# Copy `Project.toml` to notebooks
cp(joinpath(@__DIR__, "Project.toml"), joinpath(notebooks_path, "Project.toml"))

# Add a README file to notebooks
open(joinpath(notebooks_path, "README.md"), "w") do io
    print(io, """
    # Convex.jl example notebooks

    Start Julia in this directory and set the project flag to point to this directory. E.g. run the command

    ```julia
    julia --project=.
    ```

    in this directory.

    Then add `IJulia` if it's not installed already in your global environment by

    ```julia
    pkg> add IJulia
    ```

    Also call `instantiate` to download the required packages:

    ```julia
    pkg> instantiate
    ```

    Then launch Jupyter:

    ```julia
    julia> using IJulia

    julia> notebook(dir=pwd(); detached=true)
    ```

    This should allow you to try any of the notebooks.
    """)
end

# zip up the notebooks directory
zip_path = joinpath(build_path, "notebooks.zip")
run(Cmd(`zip $zip_path -r notebooks`; dir = @__DIR__))

@info "[Examples] Preparing markdown files..."

for dir in readdir(literate_path)
    dir_path = joinpath(literate_path, dir)
    isdir(dir_path) || continue
    @info "Processing directory $dir"
    build_dir = joinpath(build_path, dir)
    isdir(build_dir) || mkdir(build_dir)
    for file in readdir(dir_path)
        file_path = joinpath(dir_path, file)
            out_path = joinpath(build_path, dir, file)
        if endswith(file, ".jl")
            postprocess = function(content)
                """
                All of the examples can be found in Jupyter notebook form [here](../$(filename(zip_path)).zip).

                ```@setup $(filename(file))
                __START_TIME = time_ns()
                @info "Starting example $(filename(file))"
                ```
                """ * content * """
                ```@setup $(filename(file))
                __END_TIME = time_ns()
                elapsed = string(round((__END_TIME - __START_TIME)*1e-9; sigdigits = 3), "s")
                @info "Finished example $(filename(file)) after " * elapsed
                ```
                """
            end
            Literate.markdown(file_path, build_dir; preprocess = fix_math_md, documenter = true, postprocess =  postprocess)
        else
            cp(file_path, out_path)
        end
    end
end

@info "Starting `makedocs`"

# Build nav tree for examples
function nav_dir(dir, path)
    sort([ joinpath("examples", dir, file) for file in readdir(path) if endswith(file, ".md") && file != "index.md" ])
end

examples_nav = [ filename_to_name(dir) => nav_dir(dir, joinpath(build_path, dir)) for dir in readdir(build_path) if isdir(joinpath(build_path, dir)) ]

makedocs(;
    modules = [Convex],
    format = Documenter.HTML(),
    pages = [
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Quick Tutorial" => "quick_tutorial.md",
        "Basic Types" => "types.md",
        "Supported Operations" => "operations.md",
        "Complex-domain Optimization" => "complex-domain_optimization.md",
        "Solvers" => "solvers.md",
        "FAQ" => "faq.md",
        "Advanced" => "advanced.md",
        "Problem Depot" => "problem_depot.md",
        "Contributing" => "contributing.md",
        "Credits" => "credits.md",
        "Examples" => examples_nav,
    ],
    repo = "https://github.com/JuliaOpt/Convex.jl/blob/{commit}{path}#L{line}",
    sitename = "Convex.jl")

deploydocs(repo = "github.com/JuliaOpt/Convex.jl.git", push_preview = true)

# restore the environmental variable `GKSwstype`.
ENV["GKSwstype"] = previous_GKSwstype;


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ e0e61a2a-5b42-11eb-3e21-a5af0c3ba8fb
begin
	import Pkg
	Pkg.add("InteractiveUtils")
	Pkg.add("Interact")
 using Interact
	Pkg.add("Images")
 using Images
	Pkg.add("WebIO")
	using WebIO
end

# ╔═╡ 5cd77430-5b9d-11eb-0661-eb3cb1b47151
Pkg.installed()

# ╔═╡ fc21549e-5683-11eb-07bf-51233c3f5e9b

md"""
## My Notebook $~~~~~~~$ Jan 13 2021
[^Page1]:[^Index]

For no brake spaces use tildes to

get spaces between two dollar signs.

Leave a blank line for a line feed.

This text has a footnote [^1]

This is a reference [^Barringer]

Math symbols are rendered between two $ signs.

## $ \alpha = \beta $

$\mathscr David~\mathscr Marc~\mathscr Tan~\mathscr~\mathscr Creti$

$\mathfrak David~\mathfrak Marc~\mathfrak Tan~\mathfrak Creti$

$\mathcal David~\mathcal Marc~\mathcal Tan~\mathcal Creti$

$$\mathbb Greetings. $$

$$This~display~has~~~~~spaces~~~~~made~with~tildes!$$

---

Hello from between horizontal rulers.

---

``\LaTeX\ \ \ \ A = \pi r^2   ``

```math
f(a) = \frac{1}{2\pi}\int_{0}^{2\pi} (\alpha+R\cos(\theta))d\theta
```

[^Page2]: [^Index]
# Level One
## Level Two
### Level Three
#### Level Four
##### Level Five
###### Level Six

Here's a quote:
> Julia is a high-level, high-performance dynamic programming language for
> technical computing, with syntax that is familiar to users of other
> technical computing environments.

+ Here are some lists.
  + item one   two spaces before
    + sub 1
    + sub 2
    + sub 3
  + item two
  + item three

  1)
  1)
  1)

---
!!! note "Note"
    This is the content of the note.

!!! warning "Beware!"
    This is your last warning!

!!! hint "MyHint"
    This is the content of the hint.

!!! help "Help!"
    How may we help you?

!!! advice "Advice"
    A word to the wise is sufficient!

!!! more "Additional Help"
    Try again

    and again!

!!! divider ""

```julia
 for s ∈ ["foo","bar","baz"]
           println(s)
       end
```

|----Column One---|----Column Two----|----Column Three----|
|:----------:      | :----------      |:------------:      |
| *Row* 2         | **Row** 2        |``C2``Column,Row     |
| *Row* 3         | **Row** 3        |``C3``Column,Row     |

### Footnotes & Index
[^Index]:**INDEX:** [^Page1] Top [^Page2] Mid
[^1]:Footnote One
[^Barringer]:Reference Article

"""

# ╔═╡ 58484c68-5bb9-11eb-2cea-6df30770eaaa
md"""
## Don't Pet the Kitty!
!["A pretty tiger"](https://upload.wikimedia.org/wikipedia/commons/5/56/Tiger.50.jpg)
"""

# ╔═╡ a634c51e-5bb9-11eb-185b-7777c0b5a778
md"""
## My bull Bovis.
!["A pretty tiger"](https://thumbs2.imgbox.com/c6/6f/T13ZjubW_t.jpg)
"""

# ╔═╡ b61675e2-5a4b-11eb-1983-b96c71b4e0d1
myVariable =66

# ╔═╡ 2ad843c6-5b00-11eb-3092-51f414ed4efe
myString = "My variable is"

# ╔═╡ 9cb89e92-59f6-11eb-14b4-0dc23f5c914d
s2 = "$myString: $myVariable"

# ╔═╡ 02b028e4-5b02-11eb-0fcf-9d8d07ff2817
s2

# ╔═╡ Cell order:
# ╠═e0e61a2a-5b42-11eb-3e21-a5af0c3ba8fb
# ╠═5cd77430-5b9d-11eb-0661-eb3cb1b47151
# ╠═fc21549e-5683-11eb-07bf-51233c3f5e9b
# ╠═58484c68-5bb9-11eb-2cea-6df30770eaaa
# ╠═a634c51e-5bb9-11eb-185b-7777c0b5a778
# ╠═b61675e2-5a4b-11eb-1983-b96c71b4e0d1
# ╠═2ad843c6-5b00-11eb-3092-51f414ed4efe
# ╠═9cb89e92-59f6-11eb-14b4-0dc23f5c914d
# ╠═02b028e4-5b02-11eb-0fcf-9d8d07ff2817

max_stars_repo_name: paradocs/my-first-binder

max_issues_repo_name: paradocs/my-first-binder

max_forks_repo_name: paradocs/my-first-binder


text:
@define HeaderCoefficientsRK begin
    name::Symbol
    o::Int
    s::Int
end

@define CoefficientsRK begin
    a::Matrix{T}
    b::Vector{T}
    c::Vector{T}
end


text:
include("../utils.jl")

function Ω(::SymplecticStiefel, p, X)
    Q = SymplecticMatrix(X, p)
    pT_p = lu(p' * p)

    inv_pTp_pT = pT_p \ p'
    X_inv_pTp_pT = X * inv_pTp_pT

    Ω = Q * (X_inv_pTp_pT') * (I - Q' * p * inv_pTp_pT * Q) * Q
    Ω .+= X_inv_pTp_pT
    return Ω
end

function exp_naiive!(M::SymplecticStiefel, q, p, X)
    Ω_X = Ω(M, p, X)
    q .= exp(Ω_X - Ω_X') * exp(Array(Ω_X')) * p
    return q
end

@testset "SymplecticStiefel" begin
    @testset "Real" begin
        SpSt_6_4 = SymplecticStiefel(2 * 3, 2 * 2)

        p_6_4 = Array{Float64}(
            [
                0 0 -5 -1
                0 0 9 -2
                0 0 -2 1
                -2 -9 -3 6
                -3 -13 -21 9
                -8 -36 18 -6
            ],
        )
        q_6_4 = Array{Float64}([
            0 0 -3 1
            0 0 8 -3
            0 0 -2 1
            -1 -4 -6 3
            -1 -3 -21 9
            -2 -6 18 -6
        ])
        X1 = [
            0.0 0.0 4.25 4.25
            0.0 0.0 0.125 0.125
            0.0 0.0 -1.125 -1.125
            3.625 18.125 -10.875 -10.875
            5.0 25.0 -9.0 -9.0
            13.5 67.5 4.5 4.5
        ]
        X2 = [
            -0.02648060 0.00416977 0.01130802 0.01015956
            0.01718954 -0.00680433 0.02364406 -0.00083272
            0.00050392 0.00191916 -0.01035902 -0.00079734
            0.01811917 -0.02307032 -0.04297277 -0.05409099
            -0.02529516 0.00959934 -0.08594555 -0.06117803
            -0.02823014 0.00029946 -0.04196034 -0.04145413
        ]

        points = [
            [
                0 0 2 6
                0 0 1 -4
                0 0 0 1
                -1 0 -7 1
                1 0 -4 -24
                10 -1 4 -4
            ],
            [
                0 0 5 8
                0 0 -2 -6
                0 0 0 1
                1 0 -4 3
                3 0 -19 -34
                10 -1 1 -6
            ],
            [
                5 -3 0 0
                11 -6 0 0
                -2 1 0 0
                8 -5 -3 -3
                0 0 2 1
                5 -3 3 -2
            ],
            [
                1 -1 0 0
                3 -2 0 0
                0 0 0 0
                2 -2 -2 -3
                0 0 1 1
                1 -1 1 -2
            ],
            [
                0 0 0 0
                0 0 -1 0
                0 0 4 -1
                -1 -10 -5 2
                1 4 -2 0
                0 1 -16 4
            ],
            [
                0 0 0 0
                0 0 -1 0
                0 0 4 -1
                0 -6 -5 2
                1 4 -2 0
                0 1 -16 4
            ],
        ]

        close_points = [
            [
                -2 1 -2 -14
                1 0 1 11
                1 0 4 9
                0 0 0 1
                0 0 0 -1
                0 0 1 3
            ],
            [
                -3 1 -2 -10
                1 0 1 7
                1 0 4 10
                0 0 0 1
                0 0 0 0
                0 0 1 3
            ],
            [
                -2 1 -6 -26
                1 0 5 23
                3 -1 2 3
                0 0 1 4
                0 0 0 -1
                0 0 1 3
            ],
        ]

        @testset "Basics" begin
            @test repr(SpSt_6_4) == "SymplecticStiefel{6, 4, ℝ}()"
            @test representation_size(SpSt_6_4) == (6, 4)
            @test base_manifold(SpSt_6_4) === SpSt_6_4

            @test is_point(SpSt_6_4, p_6_4)
            @test_throws DomainError is_point(SpSt_6_4, 2 * p_6_4, true)

            @test is_vector(SpSt_6_4, p_6_4, X1; atol=1.0e-12)
            @test is_vector(SpSt_6_4, p_6_4, X2; atol=1.0e-6)
            @test_throws DomainError is_vector(SpSt_6_4, p_6_4, X2, true; atol=1.0e-12)
            @test is_vector(SpSt_6_4, p_6_4, X1 + X2; atol=1.0e-6)
            @test_throws DomainError is_vector(SpSt_6_4, p_6_4, X1 + p_6_4, true)
        end
        @testset "Symplectic Inverse" begin
            I_2k = Array(I, 4, 4)
            @test Manifolds.symplectic_inverse_times(SpSt_6_4, p_6_4, p_6_4) == I_2k
            @test Manifolds.symplectic_inverse_times!(
                SpSt_6_4,
                zeros(4, 4),
                p_6_4,
                p_6_4,
            ) == I_2k
            @test inv(SpSt_6_4, p_6_4) * p_6_4 == I_2k
            @test inv!(SpSt_6_4, copy(p_6_4'), p_6_4) * p_6_4 == I_2k
        end
        @testset "Embedding" begin
            x = [
                1 1 -9 7
                -1 -1 -13 11
                7 8 -22 19
                0 0 7 -6
                0 0 -1 1
                0 0 -1 1
            ]
            y = similar(x)
            z = embed(SpSt_6_4, x)
            @test z == x

            Y = similar(X1)
            embed!(SpSt_6_4, Y, p_6_4, X1)
            @test Y == X1
        end
        @testset "Retractions and Exponential Mapping" begin
            @test isapprox(retract(SpSt_6_4, p_6_4, X1), q_6_4; atol=1.0e-12)
            @test isapprox(
                retract(SpSt_6_4, p_6_4, X1, CayleyRetraction()),
                q_6_4;
                atol=1.0e-12,
            )

            X_inv_cayley_retraction = inverse_retract(SpSt_6_4, p_6_4, q_6_4)
            X_inv_cayley_retraction_2 =
                inverse_retract(SpSt_6_4, p_6_4, q_6_4, CayleyInverseRetraction())
            @test isapprox(X_inv_cayley_retraction, X_inv_cayley_retraction_2; atol=1.0e-16)
            @test isapprox(X_inv_cayley_retraction, X1; atol=1.0e-12)
        end
        @testset "Riemannian Metric" begin
            X1_norm = 37.85466645
            @test isapprox(norm(SpSt_6_4, p_6_4, X1), X1_norm; atol=1.0e-8)
            @test isapprox(
                norm(SpSt_6_4, p_6_4, X1),
                √inner(SpSt_6_4, p_6_4, X1, X1);
                atol=1.0e-8,
            )

            X2_norm = 1.0
            @test isapprox(norm(SpSt_6_4, p_6_4, X2), X2_norm; atol=1.0e-6)
            @test isapprox(
                norm(SpSt_6_4, p_6_4, X2),
                √inner(SpSt_6_4, p_6_4, X2, X2);
                atol=1.0e-6,
            )

            # Project Project matrix A ∈ ℝ^{6 × 4} onto (T_pSpSt):
            A_6_4 = Array{Float64}([
                -7 2 12 0
                4 0 1 -2
                -1 -1 4 0
                -18 4 -1 5
                7 0 -2 11
                2 2 -2 9
            ])
            A_6_4_proj = similar(A_6_4)
            Manifolds.project!(SpSt_6_4, A_6_4_proj, p_6_4, A_6_4)
            @test is_vector(SpSt_6_4, p_6_4, A_6_4_proj, true; atol=1.0e-12)
        end
        @testset "Generate random points/tangent vectors" begin
            M_big = SymplecticStiefel(20, 10)
            p_big = rand(M_big)
            @test is_point(M_big, p_big, true; atol=1.0e-14)
            X_big = rand(M_big; vector_at=p_big, hamiltonian_norm=1.0)
            @test is_vector(M_big, p_big, X_big, true; atol=1.0e-14)
        end
        @testset "test_manifold(Symplectic(6), ...)" begin
            types = [Matrix{Float64}]
            TEST_FLOAT32 && push!(types, Matrix{Float32})
            TEST_STATIC_SIZED && push!(types, MMatrix{6,4,Float64,24})
            for type in types
                @testset "Type $(type)" begin
                    @testset "CayleyRetraction" begin
                        test_manifold(
                            SpSt_6_4,
                            convert.(type, points);
                            retraction_methods=[CayleyRetraction()],
                            default_retraction_method=CayleyRetraction(),
                            default_inverse_retraction_method=CayleyInverseRetraction(),
                            test_inplace=true,
                            is_point_atol_multiplier=1.0e4,
                            is_tangent_atol_multiplier=1.0e3,
                            retraction_atol_multiplier=1.0e1,
                            test_project_tangent=(type != MMatrix{6,4,Float64,24}),
                            test_injectivity_radius=false,
                            test_exp_log=false,
                            test_representation_size=true,
                        )
                    end

                    @testset "ExponentialRetraction" begin
                        test_manifold(
                            SpSt_6_4,
                            convert.(type, close_points);
                            retraction_methods=[ExponentialRetraction()],
                            default_retraction_method=ExponentialRetraction(),
                            default_inverse_retraction_method=CayleyInverseRetraction(),
                            test_inplace=true,
                            is_point_atol_multiplier=1.0e11,
                            is_tangent_atol_multiplier=1.0e2,
                            retraction_atol_multiplier=1.0e4,
                            test_project_tangent=(type != MMatrix{6,4,Float64,24}),
                            test_injectivity_radius=false,
                            test_exp_log=false,
                            test_representation_size=true,
                        )
                    end
                end
            end # for
        end

        @testset "Gradient Computations" begin
            Q_grad = SymplecticMatrix(points[1])
            function test_f(p)
                k = size(p)[2]
                return tr(p[1:k, 1:k])
            end
            function analytical_grad_f(p)
                n, k = size(p)
                euc_grad_f = [Array(I, k, k); zeros((n - k), k)]
                return Q_grad * p * (euc_grad_f') * Q_grad * p + euc_grad_f * p' * p
            end
            p_grad = convert(Array{Float64}, points[1])
            fd_diff = RiemannianProjectionBackend(Manifolds.FiniteDifferencesBackend())

            @test isapprox(
                Manifolds.gradient(SpSt_6_4, test_f, p_grad, fd_diff),
                analytical_grad_f(p_grad);
                atol=1.0e-9,
            )

            grad_f_p = similar(p_grad)
            Manifolds.gradient!(SpSt_6_4, test_f, grad_f_p, p_grad, fd_diff)
            @test isapprox(grad_f_p, analytical_grad_f(p_grad); atol=1.0e-9)
        end
    end
end


text:
include("../SimpleTasks.jl")

module RunHybridDaemon

using SimpleTasks.Types
using SimpleTasks.Services.AWSQueue
using SimpleTasks.Services.CLIBucket
using SimpleTasks.Services.GCSCLIProvider
using SimpleTasks.Services.FileSystemCache
using SimpleTasks.Services.BucketCacheDatasource
using SimpleTasks.Services.Daemon
using SimpleTasks.Examples.NoOpTask

import AWS

type RunConfig
    task_queue_name::ASCIIString
    error_queue_name::ASCIIString
    done_queue_name::ASCIIString
    bucket_name::ASCIIString
    cache_directory::ASCIIString
    poll_frequency_seconds::Int64
end

#=
 = Create the queue and bucket service and start the daemon
 =#
function run(task_queue_name, error_queue_name, done_queue_name,
        bucket_name, cache_directory, poll_frequency_seconds)
    # Load AWS credentials via AWS library (either through environment
    # variables or ~/.awssecret or query permissions server)
    env = AWS.AWSEnv()

    # Create a queue to read tasks from
    task_queue = AWSQueueService(env, task_queue_name)

    # Create a queue to write errors to
    error_queue = AWSQueueService(env, error_queue_name)

    # Create a queue to report task completion
    done_queue = AWSQueueService(env, done_queue_name)

    # Create a datasource to read and write data from
    bucket = CLIBucketService(GCSCLIProvider.Details(), bucket_name)
    cache = FileSystemCacheService(cache_directory)
    datasource = BucketCacheDatasourceService(bucket, cache)

    # create a daemon to run tasks
    daemon = DaemonService(task_queue, error_queue, done_queue,
        bucket, datasource, poll_frequency_seconds)

    # Register the NOOP task
    register!(daemon, NoOpTask.NAME, NoOpTaskDetails)

    Daemon.run(daemon)
end

#=
 =Parse ARG into run daemon configuration.
 =Return RunConfig
 =#
function parse_args()
    if length(ARGS) < 6
        error("Not enough arguments given, (given $ARGS) sample usage:
            -- julia daemon.jl task_queue_name error_queue_name " *
            "done_queue_name bucket_name cache_directory " *
            "poll_frequency_seconds")
    end

    run_config = RunConfig(
        ASCIIString(ARGS[1]),
        ASCIIString(ARGS[2]),
        ASCIIString(ARGS[3]),
        ASCIIString(ARGS[4]),
        ASCIIString(ARGS[5]),
        parse(Int64, ARGS[6])
    )
    return run_config
end

#Run main!
function __init__()
    run_config = parse_args()
    #=
     =run_config = RunConfig("task-queue-TEST", "error-queue-TEST",
     =    "done-queue-TEST", "seunglab-test", "/var/tmp/taskdaemon", 5)
     =#
    run(run_config.task_queue_name, run_config.error_queue_name,
        run_config.done_queue_name, run_config.bucket_name,
        run_config.cache_directory, run_config.poll_frequency_seconds)
end

end # end module RunDaemon

max_forks_repo_name: seung-lab/SimpleTasks


text:
function Sigmoid(x)
   g=1 ./ (1 .+exp.(-x))
    return g
end

function SigmoidGradient(x)

    g=Sigmoid(x).*(1 .- Sigmoid(x))
    return g

end


text:
module TestMightReturnReduced

using Test
using Transducers
using Transducers: _might_return_reduced

@testset "_might_return_reduced" begin
    @test !_might_return_reduced(reducingfunction(Map(identity), +), 0, 1:1)
    @test !_might_return_reduced(reducingfunction(Filter(isodd), +), 0, 1:1)
    @test _might_return_reduced(reducingfunction(ReduceIf(isodd), +), 0, 1:1)
    @test _might_return_reduced(reducingfunction(ReduceIf(ismissing), +), 0, [0, missing])
    @test _might_return_reduced(reducingfunction(ReduceIf(ismissing), +), 0, Any[])
    @test !_might_return_reduced(reducingfunction(ReduceIf(ismissing), +), 0, 1:1)
    @test _might_return_reduced(reducingfunction(ReduceIf(ismissing), +), 0, (0, missing))
    @test _might_return_reduced(reducingfunction(ReduceIf(ismissing), +), 0, (missing, 0))
    @test !_might_return_reduced(
        reducingfunction(ReduceIf(ismissing), +),
        0,
        ntuple(identity, 20),
    )
    @test _might_return_reduced(
        reducingfunction(ReduceIf(ismissing), +),
        0,
        (ntuple(identity, 20)..., missing),
    )
    # It doesn't matter if it's going to throw:
    @test _might_return_reduced(reducingfunction(Map(throw), +), 0, 1:1) isa Bool
    @test _might_return_reduced(reducingfunction(Map(_ -> _uninferable_), +), 0, 1:1)
end

end  # module


text:
"""
    DLV(ϵ::Float64)

DLV searches layer by layer for counter examples in hidden layers.

# Problem requirement
1. Network: any depth, any activation (currently only support ReLU)
2. Input: hyperrectangle
3. Output: abstractpolytope

# Return
`CounterExampleResult`

# Method
The following operations are performed layer by layer. for layer i
1. determine a reachable set from the reachable set in layer i-1
2. determine a search tree in the reachable set by refining the search tree in layer i-1
3. Verify
    - True -> continue to layer i+1
    - False -> counter example

The argument `ϵ` is the resolution of the initial search tree. Default `1.0`.

# Property
Sound but not complete.

# Reference
[X. Huang, M. Kwiatkowska, S. Wang, and M. Wu,
"Safety Verification of Deep Neural Networks,"
in *International Conference on Computer Aided Verification*, 2017.](https://arxiv.org/abs/1610.06940)

[https://github.com/VeriDeep/DLV](https://github.com/VeriDeep/DLV)
"""
@with_kw struct DLV
    ϵ::Float64 = 1.0
end
# TODO: create types for the two mapping cases, since they are now both unstable and boxed
# also check out how to get performance out of closures, since that can be an issue in julia

function solve(solver::DLV, problem::Problem)
    # The list of etas
    η = get_bounds(problem)
    # The list of sample intervals
    δ = Vector{Vector{Float64}}(undef,length(η))
    δ[1] = fill(solver.ϵ, dim(η[1]))

    if issubset(last(η), problem.output)
        return CounterExampleResult(:SAT)
    end

    output = compute_output(problem.network, problem.input.center)
    for (i, layer) in enumerate(problem.network.layers)
        δ[i+1] = get_manipulation(layer, δ[i], η[i+1])
        if i == length(problem.network.layers)
            mapping = x -> (x ∈ problem.output)
        else
            forward_nnet = Network(problem.network.layers[i+1:end])
            mapping = x -> (compute_output(forward_nnet, x) ∈ problem.output)
        end
        var, y = bounded_variation(η[i+1], mapping, δ[i+1])  # TODO rename "var"

        if var
            backward_nnet = Network(problem.network.layers[1:i])
            status, x = backward_map(y, backward_nnet, η[1:i+1])
            if status
                return CounterExampleResult(:UNSAT, x)
            end
        end
    end
    return ReachabilityResult(:UNSAT, [last(η)])
end

# For simplicity, we just cut the sample interval into half
function get_manipulation(layer::Layer, δ::Vector{Float64}, bound::Hyperrectangle)
    δ_new = abs.(layer.weights) * δ ./ 2
    return δ_new
end

# Try to find an input x that arrives at output y
function backward_map(y::Vector{Float64}, nnet::Network, bounds::Vector{Hyperrectangle})
    output = Hyperrectangle(y, zeros(size(y)))
    input = first(bounds)
    model = Model(solver = GLPKSolverMIP())
    neurons = init_neurons(model, nnet)
    deltas  = init_deltas(model, nnet)
    add_set_constraint!(model, input, first(neurons))
    add_set_constraint!(model, output, last(neurons))
    encode_mip_constraint!(model, nnet, bounds, neurons, deltas)
    o = max_disturbance!(model, first(neurons) - input.center)
    status = solve(model, suppress_warnings = true)
    if status == :Optimal
        return (true, getvalue(first(neurons)))
    else
        return (false, [])
    end
end

function bounded_variation(bound::Hyperrectangle, mapping::Function, δ::Vector{Float64})
    # step 1: check whether the boundary points have the same class
    var, y = uniform_boundary_class(bound, mapping)
    var && return (var, y)
    # step 2: check the 0-variation in all dimension
    var, y = zero_variation(bound, mapping, δ)
    var && return (var, y)
    return (false, similar(y, 0))
end

function uniform_boundary_class(bound::Hyperrectangle, mapping::Function)
    y = bound.center
    for i = 1:dim(bound)
        y[i] += bound.radius[i]
        mapping(y) || return (true, y)
        y[i] -= 2 * bound.radius[i]
        mapping(y) || return (true, y)
        y[i] += bound.radius[i]
    end
    return (false, similar(y, 0))
end

function zero_variation(bound::Hyperrectangle, mapping::Function, δ::Vector{Float64})
    y = bound.center
    for i = 1:dim(bound)
        z = deepcopy(y)
        while maximum(z - high(bound)) < 0
            z[i] += δ[i]
            mapping(z) || return (true, z)
        end
        z = deepcopy(y)
        while minimum(z - high(bound)) > 0
            z[i] -= δ[i]
            mapping(z) || return (true, z)
        end
    end
    return (false, similar(y, 0))
end


text:
module RocketBufferActorTest

using Test
using Rocket

@testset "BufferActor" begin

    println("Testing: actor BufferActor")

    @testset begin
        source = of([ 1, 2, 3 ])
        actor  = BufferActor{Int}(3)

        subscribe!(source, actor)

        @test actor.values == [ 1, 2, 3 ]
        @test getvalues(actor) == [ 1, 2, 3 ]
    end

    @testset begin
        source = from([ [ 1, 2, 3 ], [ 4, 5, 6 ] ])
        actor  = BufferActor{Int}(3)

        subscribe!(source, actor)

        @test actor.values == [ 4, 5, 6 ]
        @test getvalues(actor) == [ 4, 5, 6 ]
    end

    @testset begin
        source = faulted(Vector{Int}, "Error")
        actor  = BufferActor{Int}(3)

        @test_throws ErrorException subscribe!(source, actor)
    end

    @testset begin
        @test buffer(Int, 3) isa BufferActor{Int}
    end
end

end