text:
### Generate measurements
@inbounds function rungenerator(system, measurements, num, numsys, settings, info; flow = 0)
    if settings.runflow == 1
        branch = flow.flow; bus = flow.main
        busi, Vi, Ti, Pinj, Qinj, branchi, fromi, toi, Pij,
            Qij, Pji, Qji, Iij, Dij, Iji, Dji = read_generator(bus, branch)

        for i = 1:numsys.Nbranch
            measurements.legacyFlow[i, 1] = branchi[i]
            measurements.legacyFlow[i + numsys.Nbranch, 1] = branchi[i]
            measurements.legacyFlow[i, 2] = fromi[i]
            measurements.legacyFlow[i + numsys.Nbranch, 2] = toi[i]
            measurements.legacyFlow[i, 3] = toi[i]
            measurements.legacyFlow[i + numsys.Nbranch, 3] = fromi[i]
            measurements.legacyFlow[i, 10] = Pij[i] / system.basePower
            measurements.legacyFlow[i + numsys.Nbranch, 10] = Pji[i] / system.basePower
            measurements.legacyFlow[i, 11] = Qij[i] / system.basePower
            measurements.legacyFlow[i + numsys.Nbranch, 11] = Qji[i] / system.basePower

            measurements.legacyCurrent[i, 1] = branchi[i]
            measurements.legacyCurrent[i + numsys.Nbranch, 1] = branchi[i]
            measurements.legacyCurrent[i, 2] = fromi[i]
            measurements.legacyCurrent[i + numsys.Nbranch, 2] = toi[i]
            measurements.legacyCurrent[i, 3] = toi[i]
            measurements.legacyCurrent[i + numsys.Nbranch, 3] = fromi[i]
            measurements.legacyCurrent[i, 7] = Iij[i]
            measurements.legacyCurrent[i + numsys.Nbranch, 7] = Iji[i]

            measurements.pmuCurrent[i, 1] = branchi[i]
            measurements.pmuCurrent[i + numsys.Nbranch, 1] = branchi[i]
            measurements.pmuCurrent[i, 2] = fromi[i]
            measurements.pmuCurrent[i + numsys.Nbranch, 2] = toi[i]
            measurements.pmuCurrent[i, 3] = toi[i]
            measurements.pmuCurrent[i + numsys.Nbranch, 3] = fromi[i]
            measurements.pmuCurrent[i, 10] = Iij[i]
            measurements.pmuCurrent[i + numsys.Nbranch, 10] = Iji[i]
            measurements.pmuCurrent[i, 11] = (pi / 180) * Dij[i]
            measurements.pmuCurrent[i + numsys.Nbranch, 11] = (pi / 180) * Dji[i]
        end

        for i = 1:numsys.Nbus
            measurements.legacyInjection[i, 1] = busi[i]
            measurements.legacyInjection[i, 8] = Pinj[i] / system.basePower
            measurements.legacyInjection[i, 9] = Qinj[i] / system.basePower

            measurements.legacyVoltage[i, 1] = busi[i]
            measurements.legacyVoltage[i, 5] = Vi[i]

            measurements.pmuVoltage[i, 1] = busi[i]
            measurements.pmuVoltage[i, 8] = Vi[i]
            measurements.pmuVoltage[i, 9] = (pi / 180) * Ti[i]
        end
    end

    names = Dict(:pmu => [:pmuVoltage, :pmuCurrent], :legacy => [:legacyFlow, :legacyCurrent, :legacyInjection, :legacyVoltage])
    measure = [:pmuVoltage, :pmuCurrent, :legacyFlow, :legacyCurrent, :legacyInjection, :legacyVoltage]
    pmuv = 1; pmuc = 1; legf = 1; legc = 1; legi = 1; legv = 1
    if num.pmuNv == 0 measure[1] = :no; pmuv = 0  end
    if num.pmuNc == 0 measure[2] = :no; pmuc = 0 end
    if num.legacyNf == 0 measure[3] = :no; legf = 0 end
    if num.legacyNc == 0 measure[4] = :no; legc = 0 end
    if num.legacyNi == 0 measure[5] = :no; legi = 0 end
    if num.legacyNv == 0 measure[6] = :no; legv = 0 end
    names[:legacy] = intersect(names[:legacy], measure)
    names[:pmu] = intersect(names[:pmu], measure)

    measurements = runset(system, measurements, numsys, settings, names)
    measurements = runvariance(system, measurements, settings, names)
    info = infogenerator(system, measurements, settings, info, names)

    gen = 2
    if numsys.Ngen == 0
        gen = 0
    end
    cost = 2
    if numsys.Ncost == 0
        cost = 0
    end
    group = (pmuVoltage = pmuv, pmuCurrent = pmuc, legacyFlow = legf, legacyCurrent = legc, legacyInjection = legi, legacyVoltage = legv, bus = 2, branch = 2, generator = gen, generatorcost = cost, basePower = 2)

    println("Measurment data is successfully generated.")
    if !isempty(settings.save)
        mheader = measureheader(); pheader = psheader(system); header = merge(mheader, pheader)
        savedata(measurements, system; group = group, header = header, path = settings.save, info = info)
    end

    return info
end

### Produce  measurement set
@inbounds function runset(system, measurements, numsys, settings, names)
    write = (pmuVoltage = [4, 7], pmuCurrent = [6, 9], legacyFlow = [6, 9], legacyCurrent = [6], legacyInjection = [4, 7], legacyVoltage = [4])
    type = (pmucomplete = :pmu,  legacycomplete = :legacy, pmuredundancy = :pmu,  legacyredundancy = :legacy)

    for set in keys(settings.set)
        if set in [:pmucomplete :legacycomplete]
            for i in names[type[set]]
                T::Array{Float64,2} = getfield(measurements, i)
                T[:, write[i]] = ones(size(T, 1), length(write[i]))
            end
        end
        if set in [:pmuVoltage, :pmuCurrent, :legacyFlow, :legacyCurrent, :legacyInjection, :legacyVoltage]
            for (pos, howMany) in enumerate(settings.set[set])
                success = false
                T::Array{Float64,2} = getfield(measurements, set)
                N = size(T, 1)
                if 0 <= howMany <= N
                    success = true
                    T[:, write[set][pos]] = zeros(N)
                    idxn = randperm(N)
                    idx = idxn[1:Int64(howMany)]
                    T[idx, write[set][pos]] = ones(N)[idx]
                elseif howMany > N || howMany == -200.0
                    success = true
                    T[:, write[set][pos]] = zeros(N)
                    T[:, write[set][pos]] = ones(N)
                elseif !success && howMany != -100
                    throw(ErrorException("the name-value pair setting is incorect, deployment measurements according to TYPE failed"))
                end
            end
        end
        if set in [:pmuredundancy :legacyredundancy]
            Nmax = 0
            for i in names[type[set]]
                Nmax += length(write[i]) * size(getfield(measurements, i), 1)
            end
            success = false
            howMany = settings.set[set][1]
            if 0 <= howMany <= Nmax / (2 * numsys.Nbus - 1)
                success = true
                idxn = randperm(Nmax)
                idx = idxn[1:trunc(Int, round((2 * numsys.Nbus - 1) * howMany))]
                total_set = fill(0.0, Nmax)
                total_set[idx] = ones(Nmax)[idx]
            elseif howMany > Nmax / (2 * numsys.Nbus - 1) || howMany == -200.0
                success = true
                total_set = fill(1.0, Nmax)
            end
            if success
                start = 0
                last = 0
                for i in names[type[set]]
                    for k in write[i]
                        T::Array{Float64,2} = getfield(measurements, i)
                        last += size(T, 1)
                        T[:, k] = total_set[(start + 1):last]
                        start += size(T, 1)
                    end
                end
            else
                throw(ErrorException("the name-value pair setting is incorect, deployment measurements according to REDUNDANCY failed"))
            end
        end
        if set == :pmudevice
            if :pmuVoltage in names[:pmu] && :pmuCurrent in names[:pmu]
                success = false
                howMany = settings.set[set][1]
                N = size(measurements.pmuVoltage, 1)
                if 0 <= howMany <= N
                    success = true
                    idxn = randperm(N)
                    idx = idxn[1:Int64(howMany)]
                elseif howMany > N || howMany == -200.0
                    success = true
                    idx = collect(1:N)
                end
                if success
                    Nc = size(measurements.pmuCurrent, 1)
                    measurements.pmuVoltage[:, [4, 7]] = zeros(N, 2)
                    measurements.pmuCurrent[:, [6, 9]] = zeros(Nc, 2)
                    labelVol = measurements.pmuVoltage[:, 1]
                    labelCurr = measurements.pmuCurrent[:, 2]
                    try
                        labelVol = measurements.pmuVoltage[:, 10]
                        labelCurr = measurements.pmuCurrent[:, 12]
                    catch
                    end
                    for i in idx
                        measurements.pmuVoltage[i, 4] = 1.0
                        measurements.pmuVoltage[i, 7] = 1.0
                        for j = 1:size(measurements.pmuCurrent, 1)
                            if labelVol[i] == labelCurr[j]
                                measurements.pmuCurrent[j, 6] = 1.0
                                measurements.pmuCurrent[j, 9] = 1.0
                            end
                        end
                    end
                else
                    throw(ErrorException("the name-value pair setting is incorect, deployment measurements according to DEVICE failed"))
                end
            else
                throw(ErrorException("the complete PMU measurement DATA is not found"))
           end
        end
        if set == :pmuoptimal
            if :pmuVoltage in names[:pmu] && :pmuCurrent in names[:pmu]
                Nvol = size(measurements.pmuVoltage, 1)
                Ncurr = size(measurements.pmuCurrent, 1)
                if Nvol == numsys.Nbus && Ncurr == 2 * numsys.Nbranch
                    bus = collect(1:numsys.Nbus)
                    numbering = false
                    fromi = convert(Array{Int64,1}, system.branch[:, 2])
                    toi = convert(Array{Int64,1}, system.branch[:, 3])
                    busi = @view(system.bus[:, 1])
                    for i = 1:numsys.Nbus
                        if bus[i] != busi[i]
                            numbering = true
                        end
                    end
                    from = renumber(fromi, busi, bus, numbering)
                    to = renumber(toi, busi, bus, numbering)

                    A = sparse([bus; from; to], [bus; to; from], fill(1, numsys.Nbus + 2 * numsys.Nbranch), numsys.Nbus, numsys.Nbus)
                    f  = fill(1, numsys.Nbus);
                    lb = fill(0, numsys.Nbus);
                    ub = fill(1, numsys.Nbus);
                    model = Model(GLPK.Optimizer)

                    @variable(model, x[i = 1:numsys.Nbus], Int)
                    @constraint(model, lb .<= x .<= ub)
                    @constraint(model, -A * x .<= -f)
                    @objective(model, Min, sum(x))
                    optimize!(model)
                    x = JuMP.value.(x)
                    idx = findall(x->x == 1, x)

                    measurements.pmuVoltage[:, [4, 7]] = zeros(Nvol, 2)
                    measurements.pmuCurrent[:, [6, 9]] = zeros(Ncurr, 2)
                    for i in idx
                        measurements.pmuVoltage[i, 4] = 1.0
                        measurements.pmuVoltage[i, 7] = 1.0
                        for j = 1:numsys.Nbranch
                            if busi[i] == fromi[j]
                                measurements.pmuCurrent[j, 6] = 1.0
                                measurements.pmuCurrent[j, 9] = 1.0
                            end
                            if busi[i] == toi[j]
                                measurements.pmuCurrent[j + numsys.Nbranch, 6] = 1.0
                                measurements.pmuCurrent[j + numsys.Nbranch, 9] = 1.0
                            end
                        end
                    end
                else
                    throw(ErrorException("the optimal algorithm requires the concordant number of buses and branches, and corresponding measurements"))
                end
            else
                throw(ErrorException("the complete PMU measurement DATA is not found"))
            end
        end
    end

    return measurements
end

### Produce measurement variances
@inbounds function runvariance(system, measurements, settings, names)
    write = (pmuVoltage = [2 8 3; 5 9 6], pmuCurrent = [4 10 5; 7 11 8], legacyFlow = [4 10 5; 7 11 8], legacyCurrent = [4 7 5], legacyInjection = [2 8 3; 5 9 6], legacyVoltage = [2 5 3])
    type = (pmucomplete = :pmu, legacycomplete = :legacy, pmurandom = :pmu, legacyrandom = :legacy)

    for var in keys(settings.variance)
        if var in [:pmucomplete :legacycomplete]
            howMany = settings.variance[var][1]
            if howMany > 0
                for i in names[type[var]]
                    T::Array{Float64,2} = getfield(measurements, i)
                    for row in eachrow(write[i])
                        T[:, row[3]] .= howMany
                        T[:, row[1]] .= T[:, row[2]] .+ howMany.^(1/2) .* randn!(zeros(size(T, 1)))
                    end
                end
            else
                error("the variance must be positive number")
            end
        end
        if var in [:pmurandom :legacyrandom]
            min = minimum(settings.variance[var])
            max = maximum(settings.variance[var])
            if min > 0 && max > 0
                for i in names[type[var]]
                    T::Array{Float64,2} = getfield(measurements, i)
                    for row in eachrow(write[i])
                        T[:, row[3]] .= min .+ (max - min) .* rand(size(T, 1))
                        T[:, row[1]] .= T[:, row[2]] .+ T[:, row[3]].^(1/2) .* randn!(zeros(size(T, 1)))

                    end
                end
            else
                throw(ErrorException("the variance value must be positive number"))
            end
        end
        if var in [:pmuVoltage, :pmuCurrent, :legacyFlow, :legacyCurrent, :legacyInjection, :legacyVoltage]
            for (pos, howMany) in enumerate(settings.variance[var])
                if howMany != -100.0
                    if howMany > 0
                        T::Array{Float64,2} = getfield(measurements, var)
                        T[:, write[var][pos, 3]] .= howMany
                        T[:, write[var][pos, 1]] .= T[:, write[var][pos, 2]] .+ howMany.^(1/2) .* randn!(zeros(size(T, 1)))
                    else
                        throw(ErrorException("the variance value must be positive number"))
                    end
                end
            end
        end
    end

    return measurements
end

### Read data
function read_generator(bus, branch)
    busi = @view(bus[:, 1])
    Vi = @view(bus[:, 2])
    Ti = @view(bus[:, 3])
    Pinj = @view(bus[:, 4])
    Qinj = @view(bus[:, 5])

    branchi = @view(branch[:, 1])
    fromi = @view(branch[:, 2])
    toi = @view(branch[:, 3])
    Pij = @view(branch[:, 4])
    Qij = @view(branch[:, 5])
    Pji = @view(branch[:, 6])
    Qji = @view(branch[:, 7])
    Iij = @view(branch[:, 11])
    Dij = @view(branch[:, 12])
    Iji = @view(branch[:, 13])
    Dji = @view(branch[:, 14])

    return busi, Vi, Ti, Pinj, Qinj,
        branchi, fromi, toi, Pij, Qij, Pji, Qji, Iij, Dij, Iji, Dji
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
function fft_t_to_ω(arr::AbstractArray, dt, dim = 1)
	rfft(arr, dim) .* dt
end
function fft_t_to_ω(arr::OffsetArray, dt, dim = 1)
	rfft(parent(arr), dim) .* dt
end

function fft_ω_to_t(arr::AbstractArray, dν, dim = 1, Nt = 2*(size(arr, dim)-1))
	brfft(arr, Nt, dim) .* dν
end
function fft_ω_to_t(arr::OffsetArray, dν, dim = 1, Nt = 2*(size(arr, dim)-1))
	brfft(parent(arr),Nt, dim) .* dν
end

ext: jl


text:
## Fit counties with significant sero data

## Do inference in parallel

using Distributed

## add number of processors desired
# addprocs(2)
params = (exename = `nice -19 /Applications/Julia-1.7.app/Contents/Resources/julia/bin/julia`,
      dir = "./Github/KenyaCoVaccinesPrivate")
addprocs([("samandfi@10.0.0.6", 7)]; params...)
addprocs(4)

@everywhere begin
      using Pkg
      Pkg.activate(".")
end

## Load relevant code into scope for each worker

@everywhere using JLD2, LinearAlgebra, TransformVariables, LogDensityProblems, MCMCChains, Dates
@everywhere import KenyaCoVaccines

## Load data
# @everywhere begin
using DataFrames
@load("data/linelist_data_with_pos_neg_20feb_to_27sept_c_age__pcr.jld2")
@load("data/cleaned_serology_data_c_age__combined.jld2")
@load("data/deaths_20210919.jld2")
@load("data/p_IH.jld2")
@load("data/p_ID.jld2")

##Load contact structure and population size data
M_Kenya_ho, M_Kenya_other, M_Kenya_school, M_Kenya_work = KenyaCoVaccines.get_rescaled_contact_matrices("data/agemixingmatrix_Kenya_all_types.jld2")
N_kenya = KenyaCoVaccines.get_population_size_matrix("data/2019_census_age_pyramids_counties.csv")

# Parameter transformations and initial KE guess for HMC

trans = as((
      β₀ = as(Real, 0.0, 5.0),  #β₀= as(Real, 0.0, 7.0),
      β_home = as(Real, 0.0, 5.0),
      β_school = as(Real, 0.0, 5.0),
      β_other = as(Real, 0.0, 5.0),
      β_work = as(Real, 0.0, 5.0),
      ϵ = as(Real, 0.0, 0.999),
      χ = as(Real, 0.1, 30.0),
      p_test = as(Real, 0.01, 5),
      E₀ = as(Real, 0.0, 10e3),
      inc_R_αvar = as(Real, 0.0, 1.0),
      time_scale_αvar = as(Real, 0.01, 0.5),
      mid_point_αvar = as(Real, 0.0, 4.0),
      inc_R_δvar = as(Real, 0.001, 1.5),
      time_scale_δvar = as(Real, 0.01, 0.5),
      mid_point_δvar = as(Real, 0.0, 4.0),
      init_scale = as(Real, 0.0, 2.0)
))

D = Diagonal(0.1 * ones(TransformVariables.dimension(trans)))

# end

## Name for each county with substantial serological data, and matching priors

counties_with_serodata = ["Nairobi", "Mombasa",
      "Nakuru", "Uasin Gishu", "Embu", "Kisumu", "Siaya", "Kisii", "Nyeri",
      "Kilifi", "Kwale"]

priors = [fill(KenyaCoVaccines.priors_onegroup_alpha_delta_variant_cities, 2)
      fill(KenyaCoVaccines.priors_onegroup_alpha_delta_variant_outside_cities, 9)]

#Baseline seropos for low density coverage counties at December 1st
#Based on the unweighted seroprevalence from Ngere et al in Nairobi mid-Sept

init_sero_baseline = [(34 + 74) / (179 + 244), (100 + 95 + 51) / (265 + 241 + 134), 21 / 61, 9 / 40, 9 / 40, 9 / 40]

## parallelised loop for fitting each county
@distributed for i = 1:11
      name = counties_with_serodata[i]
      prior = priors[i]
      model = KenyaCoVaccines.CoVAreaModel(name, KenyaCoVaccines.ll_onegroup_newvariant_infboost, prior;
            linelist_data_with_pos_neg = linelist_data_with_pos_neg,
            serology_data = serology_data,
            average_sero_init = init_sero_baseline,
            deaths = deaths_data,
            pop_data = N_kenya,
            M_county_ho = M_Kenya_ho,
            M_county_other = M_Kenya_other,
            M_county_school = M_Kenya_school,
            M_county_work = M_Kenya_work,
            rapid = false,
            scenario = 0,
            startdate = Date(2020, 12, 1),
            enddate = Date(2021, 9, 24))

      KenyaCoVaccines.inferparameters!(model, 2000, trans, 0.05, D; serowaningrate = 0.0)
      @save("modelfits/$(name)_model.jld2", model)

end
##
rmprocs(workers())

##Tests
@everywhere counties_with_serodata = ["Nairobi", "Mombasa",
      "Nakuru", "Uasin Gishu", "Embu", "Kisumu", "Siaya", "Kisii", "Nyeri",
      "Kilifi", "Kwale"]

@everywhere priors = [fill(KenyaCoVaccines.priors_onegroup_alpha_delta_variant_cities, 2)
      fill(KenyaCoVaccines.priors_onegroup_alpha_delta_variant_outside_cities, 9)]

r = @spawnat 2 begin
      name = counties_with_serodata[1]
      prior = priors[1]
      model = KenyaCoVaccines.CoVAreaModel(name, KenyaCoVaccines.ll_onegroup_newvariant_infboost, prior;
            linelist_data_with_pos_neg = linelist_data_with_pos_neg,
            serology_data = serology_data,
            average_sero_init = init_sero_baseline,
            deaths = deaths_data,
            pop_data = N_kenya,
            M_county_ho = M_Kenya_ho,
            M_county_other = M_Kenya_other,
            M_county_school = M_Kenya_school,
            M_county_work = M_Kenya_work,
            rapid = false,
            scenario = 0,
            startdate = Date(2020, 12, 1),
            enddate = Date(2021, 9, 24))

      # KenyaCoVaccines.inferparameters!(model, 2000, trans, 0.05, D; serowaningrate = 0.0)
      # @save("modelfits/$(name)_model.jld2", model)
end

model = fetch(r)

name = counties_with_serodata[1]
prior = priors[1]
model = KenyaCoVaccines.CoVAreaModel(name, KenyaCoVaccines.ll_onegroup_newvariant_infboost, prior;
      PCR_cases = linelist_data_with_pos_neg,
      sero_cases = serology_data,
      average_sero_init = init_sero_baseline,
      deaths = deaths_data,
      pop_data = N_kenya,
      M_county_ho = M_Kenya_ho,
      M_county_other = M_Kenya_other,
      M_county_school = M_Kenya_school,
      M_county_work = M_Kenya_work,
      rapid = false,
      scenario = 0,
      startdate = Date(2020, 12, 1),
      enddate = Date(2021, 9, 24))
KenyaCoVaccines.inferparameters!(model, 2000, trans, 0.05, D; serowaningrate = 0.0)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta2
#     language: julia
#     name: julia-1.7
# ---

# %%
@show VERSION
using BenchmarkTools
pimc_sum(n) = 4sum(rand()^2 + rand()^2 ≤ 1 for _ in 1:n)/n
@btime pimc_sum(10^7)

# %%
@show VERSION
using BenchmarkTools, Random
pimc_sum(n, rng) = 4sum(rand(rng)^2 + rand(rng)^2 ≤ 1 for _ in 1:n)/n
@btime pimc_sum(10^7, $(Random.default_rng()))

# %%
@show VERSION
using BenchmarkTools, Random
function pimc_for(n, rng)
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end
@btime pimc_for(10^7, $(Random.default_rng()))

# %%
@show VERSION
@btime pimc_sum(10^7, $(MersenneTwister()))
@btime pimc_for(10^7, $(MersenneTwister()))

# %%
@show VERSION
f3(n) = sum(k^3 for k in 1:n)
f4(n) = sum(k^4 for k in 1:n)
@btime f3(10^6)
@btime f4(10^6)

# %%
@show VERSION
f3(n) = sum(k^3 for k in 1:n)
f4(n) = sum(k^4 for k in 1:n)
@btime f3($(Int128(10^6)))
@btime f4($(Int128(10^6)))

# %%
@code_typed f3(10^6)

# %%
@code_typed f4(10^6)

# %%

# %%
@show VERSION
using BenchmarkTools
pimc_sum(n) = 4sum(rand()^2 + rand()^2 ≤ 1 for _ in 1:n)/n
@btime pimc_sum(10^7)

# %%
@show VERSION
using BenchmarkTools, Random
pimc_sum(n, rng) = 4sum(rand(rng)^2 + rand(rng)^2 ≤ 1 for _ in 1:n)/n
@btime pimc_sum(10^7, $(Random.default_rng()))

# %%
@show VERSION
using BenchmarkTools, Random
function pimc_for(n, rng)
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end
@btime pimc_for(10^7, $(Random.default_rng()))

# %%
@show VERSION
@btime pimc_sum(10^7, $(MersenneTwister()))
@btime pimc_for(10^7, $(MersenneTwister()))

# %%
@show VERSION
f3(n) = sum(k^3 for k in 1:n)
f4(n) = sum(k^4 for k in 1:n)
@btime f3(10^6)
@btime f4(10^6)

# %%
@show VERSION
f3(n) = sum(k^3 for k in 1:n)
f4(n) = sum(k^4 for k in 1:n)
@btime f3($(Int128(10^6)))
@btime f4($(Int128(10^6)))

# %%
@code_typed f3(10^6)

# %%
@code_typed f4(10^6)

# %%

max_stars_repo_path: 0007/sum of generator.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0007/sum of generator.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0007/sum of generator.jl

max_forks_repo_name: genkuroki/public


text:
#### Tested moist thermodynamic profiles

#=
This file contains functions to compute all of the
thermodynamic _states_ that Thermodynamics is
tested with in runtests.jl
=#

using Random

"""
    ProfileSet

A set of profiles used to test Thermodynamics.
"""
struct ProfileSet{AFT, QPT, PT}
    z::AFT          # Altitude
    T::AFT          # Temperature
    p::AFT          # Pressure
    RS::AFT         # Relative saturation
    e_int::AFT      # Internal energy
    ρ::AFT          # Density
    θ_liq_ice::AFT  # Liquid Ice Potential temperature
    q_tot::AFT      # Total specific humidity
    q_liq::AFT      # Liquid specific humidity
    q_ice::AFT      # Ice specific humidity
    q_pt::QPT       # Phase partition
    RH::AFT         # Relative humidity
    e_pot::AFT      # gravitational potential
    u::AFT          # velocity (x component)
    v::AFT          # velocity (y component)
    w::AFT          # velocity (z component)
    e_kin::AFT      # kinetic energy
    phase_type::PT  # Phase type (e.g., `PhaseDry`, `PhaseEquil`)
end

function Base.iterate(ps::ProfileSet, state = 1)
    state > length(ps.z) && return nothing
    nt = (
        z = ps.z[state],
        T = ps.T[state],
        p = ps.p[state],
        RS = ps.RS[state],
        e_int = ps.e_int[state],
        ρ = ps.ρ[state],
        θ_liq_ice = ps.θ_liq_ice[state],
        q_tot = ps.q_tot[state],
        q_liq = ps.q_liq[state],
        q_ice = ps.q_ice[state],
        q_pt = ps.q_pt[state],
        RH = ps.RH[state],
        e_pot = ps.e_pot[state],
        u = ps.u[state],
        v = ps.v[state],
        w = ps.w[state],
        e_kin = ps.e_kin[state],
        phase_type = ps.phase_type,
    )
    return (nt, state + 1)
end
Base.IteratorSize(::ProfileSet) = Base.HasLength()
Base.length(ps::ProfileSet) = length(ps.z)

# Since we use `rand` to generate the ProfileSet,
# just initialize on the CPU, and provide convert
# function to move arrays to the GPU.
convert_profile_set(ps::ProfileSet, ArrayType, slice) = ProfileSet(
    ArrayType(ps.z[slice]),
    ArrayType(ps.T[slice]),
    ArrayType(ps.p[slice]),
    ArrayType(ps.RS[slice]),
    ArrayType(ps.e_int[slice]),
    ArrayType(ps.ρ[slice]),
    ArrayType(ps.θ_liq_ice[slice]),
    ArrayType(ps.q_tot[slice]),
    ArrayType(ps.q_liq[slice]),
    ArrayType(ps.q_ice[slice]),
    PhasePartition.(ps.q_tot[slice], ps.q_liq[slice], ps.q_ice[slice]),
    ArrayType(ps.RH[slice]),
    ArrayType(ps.e_pot[slice]),
    ArrayType(ps.u[slice]),
    ArrayType(ps.v[slice]),
    ArrayType(ps.w[slice]),
    ArrayType(ps.e_kin[slice]),
    ps.phase_type,
)

"""
    input_config(
        ArrayType;
        n=50,
        n_RS1=10,
        n_RS2=20,
        T_min=150,
        T_surface=340
    )

Return input arguments to construct profiles
"""
function input_config(
    ArrayType;
    n = 50,
    n_RS1 = 10,
    n_RS2 = 20,
    T_surface = 340,
    T_min = 150,
)
    n_RS = n_RS1 + n_RS2
    z_range = ArrayType(range(0, stop = 2.5e4, length = n))
    relative_sat1 = ArrayType(range(0, stop = 1, length = n_RS1))
    relative_sat2 = ArrayType(range(1, stop = 1.02, length = n_RS2))
    relative_sat = vcat(relative_sat1, relative_sat2)
    return z_range, relative_sat, T_surface, T_min
end

"""
    shared_profiles(
        param_set::AbstractParameterSet,
        z_range::AbstractArray,
        relative_sat::AbstractArray,
        T_surface,
        T_min,
    )

Compute profiles shared across `PhaseDry`,
`PhaseEquil` and `PhaseNonEquil` thermodynamic
states, including:
 - `z` altitude
 - `T` temperature
 - `p` pressure
 - `RS` relative saturation
"""
function shared_profiles(
    param_set::AbstractParameterSet,
    z_range::AbstractArray,
    relative_sat::AbstractArray,
    T_surface,
    T_min,
)
    FT = eltype(z_range)
    n_RS = length(relative_sat)
    n = length(z_range)
    T = similar(z_range, n * n_RS)
    p = similar(z_range, n * n_RS)
    RS = similar(z_range, n * n_RS)
    z = similar(z_range, n * n_RS)
    linear_indices = LinearIndices((1:n, 1:n_RS))
    # We take the virtual temperature, returned here,
    # as the temperature, and then compute a thermodynamic
    # state consistent with that temperature. This profile
    # will not be in hydrostatic balance, but this does not
    # matter for the thermodynamic test profiles.
    profile =
        DecayingTemperatureProfile{FT}(param_set, FT(T_surface), FT(T_min))
    for i in linear_indices.indices[1]
        for j in linear_indices.indices[2]
            k = linear_indices[i, j]
            z[k] = z_range[i]
            T[k], p[k] = profile(param_set, z[k])
            RS[k] = relative_sat[j]
        end
    end
    return z, T, p, RS
end

####
#### PhaseDry
####

"""
    PhaseDryProfiles(param_set, ::Type{ArrayType})

Returns a `ProfileSet` used to test dry thermodynamic states.
"""
function PhaseDryProfiles(
    param_set::AbstractParameterSet,
    ::Type{ArrayType},
) where {ArrayType}
    phase_type = PhaseDry

    z_range, relative_sat, T_surface, T_min = input_config(ArrayType)
    z, T_virt, p, RS =
        shared_profiles(param_set, z_range, relative_sat, T_surface, T_min)
    T = T_virt
    FT = eltype(T)
    _R_d::FT = R_d(param_set)
    _grav::FT = grav(param_set)
    ρ = p ./ (_R_d .* T)

    # Additional variables
    q_tot = similar(T)
    fill!(q_tot, 0)
    q_pt = PhasePartition_equil.(param_set, T, ρ, q_tot, phase_type)
    e_int = internal_energy.(param_set, T, q_pt)
    θ_liq_ice = liquid_ice_pottemp.(param_set, T, ρ, q_pt)
    q_liq = getproperty.(q_pt, :liq)
    q_ice = getproperty.(q_pt, :ice)
    RH = relative_humidity.(param_set, T, p, phase_type, q_pt)
    e_pot = _grav * z
    Random.seed!(15)
    u = rand(FT, size(T)) * 50
    v = rand(FT, size(T)) * 50
    w = rand(FT, size(T)) * 50
    e_kin = (u .^ 2 .+ v .^ 2 .+ w .^ 2) / 2

    return ProfileSet{typeof(T), typeof(q_pt), typeof(phase_type)}(
        z,
        T,
        p,
        RS,
        e_int,
        ρ,
        θ_liq_ice,
        q_tot,
        q_liq,
        q_ice,
        q_pt,
        RH,
        e_pot,
        u,
        v,
        w,
        e_kin,
        phase_type,
    )
end

####
#### PhaseEquil
####

"""
    PhaseEquilProfiles(param_set, ::Type{ArrayType})

Returns a `ProfileSet` used to test moist states in thermodynamic equilibrium.
"""
function PhaseEquilProfiles(
    param_set::AbstractParameterSet,
    ::Type{ArrayType},
) where {ArrayType}
    phase_type = PhaseEquil

    # Prescribe z_range, relative_sat, T_surface, T_min
    z_range, relative_sat, T_surface, T_min = input_config(ArrayType)

    # Compute T, p, from DecayingTemperatureProfile, (reshape RS)
    z, T_virt, p, RS =
        shared_profiles(param_set, z_range, relative_sat, T_surface, T_min)
    T = T_virt

    FT = eltype(T)
    _R_d = FT(R_d(param_set))
    _grav::FT = grav(param_set)
    # Compute total specific humidity from temperature, pressure
    # and relative saturation, and partition the saturation excess
    # according to temperature.
    ρ = p ./ (_R_d .* T)
    q_tot = RS .* q_vap_saturation.(Ref(param_set), T, ρ, Ref(phase_type))
    q_pt = PhasePartition_equil.(Ref(param_set), T, ρ, q_tot, Ref(phase_type))

    # Extract phase partitioning and update pressure
    # to be thermodynamically consistent with T, ρ, q_pt
    q_liq = getproperty.(q_pt, :liq)
    q_ice = getproperty.(q_pt, :ice)
    p = air_pressure.(Ref(param_set), T, ρ, q_pt)

    e_int = internal_energy.(Ref(param_set), T, q_pt)
    θ_liq_ice = liquid_ice_pottemp.(Ref(param_set), T, ρ, q_pt)
    RH = relative_humidity.(Ref(param_set), T, p, Ref(phase_type), q_pt)
    e_pot = _grav * z
    Random.seed!(15)
    u = rand(FT, size(T)) * 50
    v = rand(FT, size(T)) * 50
    w = rand(FT, size(T)) * 50
    e_kin = (u .^ 2 .+ v .^ 2 .+ w .^ 2) / 2

    return ProfileSet{typeof(T), typeof(q_pt), typeof(phase_type)}(
        z,
        T,
        p,
        RS,
        e_int,
        ρ,
        θ_liq_ice,
        q_tot,
        q_liq,
        q_ice,
        q_pt,
        RH,
        e_pot,
        u,
        v,
        w,
        e_kin,
        phase_type,
    )
end


text:
using EpiJ
using Documenter

DocMeta.setdocmeta!(EpiJ, :DocTestSetup, :(using EpiJ); recursive=true)

makedocs(;
    modules=[EpiJ],
    authors="Jacob Zelko <jacobszelko@gmail.com>, Georgia Tech Research Institute",
    repo="https://github.com/JuliaHealth/EpiJ/blob/{commit}{path}#{line}",
    sitename="EpiJ",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliahealth.org/EpiJ",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaHealth/EpiJ.git",
    devbranch = "main",
    push_preview = false,
    target = "build"
)

max_stars_repo_name: JuliaHealth/EpiJ

max_issues_repo_name: JuliaHealth/EpiJ

max_forks_repo_name: JuliaHealth/EpiJ


text:
@testset "test utilities" begin
    mutable struct NoThrowTestSet <: Test.AbstractTestSet
        results::Vector
        NoThrowTestSet(desc) = new([])
    end
    Test.record(ts::NoThrowTestSet, t::Test.Result) = (push!(ts.results, t); t)
    Test.finish(ts::NoThrowTestSet) = ts.results
    fails = @testset NoThrowTestSet begin
        # OK
        @test_throws_cuerror CUDA.ERROR_UNKNOWN throw(CuError(CUDA.ERROR_UNKNOWN))
        # Fail, wrong CuError
        @test_throws_cuerror CUDA.ERROR_UNKNOWN throw(CuError(CUDA.ERROR_INVALID_VALUE))
        # Fail, wrong Exception
        @test_throws_cuerror CUDA.ERROR_UNKNOWN error()
    end
    @test isa(fails[1], Test.Pass)
    @test isa(fails[2], Test.Fail)
    @test isa(fails[3], Test.Fail)
end

@testset "@sync" begin
  t = Base.@elapsed ret = CUDA.@sync begin
    # TODO: do something that takes a while on the GPU
    #       (need to wrap clock64 for that)
    42
  end
  @test t >= 0
  @test ret == 42

  CUDA.@sync blocking=false begin end
end

@testset "versioninfo" begin
    CUDA.versioninfo(devnull)
end


text:
#=
Base julia sizehint! + push! garbage compared to
resize! + setindex!

Remove when sizehinting in base julia isn't garbage
=#

mutable struct NestedPushVector{T, K<:PushVector{T}, V<:AbstractVector{K}} <: AbstractVector{K}
    "Vector used for storage."
    parent::V
    "Number of elements held by `parent`."
    len::Int
    "Inner sizehint"
    inner_sh::Int
    "Outer sizehint"
    outer_sh::Int
end

# Simplifications
const PV{T} = PushVector{T, Vector{T}} where {T}
const NPV{T} = NestedPushVector{T, PV{T}, Vector{PV{T}}} where {T}

function NestedPushVector{T}(inner_sh::Integer = 4, outer_sh::Integer = 4) where {T}
    inner_sh ≥ 0 && outer_sh ≥ 0 || throw(DomainError((inner_sh, outer_sh), "Invalid initial size."))
    NestedPushVector(
        [PushVector{T}(inner_sh) for _ in 1:outer_sh],
        0,
        inner_sh,
        outer_sh
    )
end

@inline Base.length(v::NestedPushVector) = v.len
@inline Base.size(v::NestedPushVector) = (v.len, )

function extend!(v::NestedPushVector{T}, n::Int) where {T}
    ish = v.inner_sh
    prev_l = length(v.parent)
    resize!(v.parent, prev_l + n)
    for i in 1:n
        v.parent[prev_l + i] = PushVector{T}(ish)
    end
    nothing
end

@inline function Base.getindex(v::NestedPushVector, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.parent[i]
end

@inline function Base.setindex!(v::NestedPushVector, x::PushVector, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.parent[i] = x
end

function freenext!(v::NestedPushVector)
    v.len += 1
    if v.len > length(v.parent)
        extend!(v, v.len)
    end
    nothing
end

function Base.push!(v::NestedPushVector, x::PushVector)
    v.len += 1
    if v.len > length(v.parent)
        extend!(v, v.len)
    end
    v.parent[v.len] = x
    v
end

function Base.empty!(v::NestedPushVector)
    resize!(v.parent, v.outer_sh)
    v.len = 0
    for pv in v.parent
        resize!(pv.parent, v.inner_sh)
        pv.len = 0
    end
    nothing
end


text:
#Demo 2: Chaining ABCD matrices
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_MAGDB = "Magnitude (dB)"
LBLAX_FREQ = "Frequency (Hz)"
color1 = cons(:a, line = set(color=:red, width=2))
color2 = cons(:a, line = set(color=:blue, width=2))

#==Input data
===============================================================================#
_logspace(start, stop, n) = 10 .^ range(start, stop=stop, length=n)
C = capacitance(2e-12)
R = 20.0
#f = (1:10)*1e9
f = _logspace(log10(1e3), log10(20e9), 100)
xrange = cons(:a, xyaxes = set(xmin=1e6, xmax=100e9))

#==Computations
===============================================================================#
#Use DataF1 (function of 1 argument) to augment expressivity:
f = DataF1(f, f)
_1 = ones(f)

#Compute ABCD matrix of series R:
sng = true #Compute using single ABCD matrix?
if sng
	#Compute single ABCD matrix:
	@show TR = series(:ABCD, impedance(R))
else
	#Compute ABCD{DataF1}:
	TR = series(:ABCD, impedance(R*_1))
end

#Compute ABCD matrix of shunt capacitors:
TC1 = shunt(:ABCD, admittance(C, f=f))
#TC2 = shunt(:ABCD, impedance(2C, f=f)) #Add asymmetry #TODO: fix multiplication of C
TC2 = shunt(:ABCD, impedance(capacitance(2*C.v), f=f)) #Add asymmetry

#Compute Pi network:
T = TC1*TR*TC2

S = Network(:S, T)
(s11, s12, s21, s22) = mx2elem(S)

#Intermediate network parameters (verify conversion routines):
nplist = [:Y, :Z, :ABCD, :H, :G, :T]

#==Generate plot
===============================================================================#
is11 = 1; is22 = 2; i𝛵 = 3; #Strip indices

plot = cons(:plot, nstrips=3,
	ystrip1 = set(axislabel=LBLAX_MAGDB, striplabel="Input Reflection (s11)"),
	ystrip2 = set(axislabel=LBLAX_MAGDB, striplabel="Output Reflection (s22)"),
	ystrip3 = set(axislabel=LBLAX_MAGDB, striplabel="Transmission Coefficients (𝛵)"),
	xaxis = set(scale=:log, label=LBLAX_FREQ),
)
push!(plot,
	cons(:wfrm, dB20(s11), color1, label="s11", strip=is11),
	cons(:wfrm, dB20(s22), color1, label="s22", strip=is22),
	cons(:wfrm, dB20(s12), color1, label="s12", strip=i𝛵),
	cons(:wfrm, dB20(s21), color2, label="s21", strip=i𝛵),
)

Sref = S

let S, s11, s12, s21, s22 #HIDEWARN_0.7
for np in nplist
	X = Network(np, Sref)
	S = Network(:S, X)
	(s11, s12, s21, s22) = mx2elem(S)
	push!(plot,
		cons(:wfrm, dB20(s11), color1, label="s11 ($np)", strip=is11),
		cons(:wfrm, dB20(s22), color1, label="s22 ($np)", strip=is22),
		cons(:wfrm, dB20(s12), color1, label="s12 ($np)", strip=i𝛵),
		cons(:wfrm, dB20(s21), color2, label="s21 ($np)", strip=i𝛵),
	)
end
end

pcoll = push!(cons(:plotcoll, title="Pi Network Cascade + Param Conversion"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Behavior"))
$(vskip(1em))

# Functions
$(vskip(1em))

$(
   showdocs([hasstate, keypress, clickable, selectable, send, recv, wire])
)

$(vskip(4em))

"""

function main(window)
    docpage(intro)
end


text:
using Test
using RenderingGeometry

include("definitions.jl")
include("../test_logging.jl")
include("pair_constructor_test.jl")
include("pair_conversion_test.jl")
include("pair_operator_test.jl")

function run_pair_tests()
    @testset "Pair tests" begin
        test_all_pair_constructors()
        test_all_pair_conversions()
        test_all_pair_operators()
    end
end


text:
using Brillouin, PlotlyJS, Test
using Brillouin.KPaths: Bravais
using StaticArrays

# NB: For the plot functionality, there's not that many meaningful things to test against
#     in the output since it's a graphical thing. So, instead, we just check that things
#     run without erroring by checking that the return type is as expected

@testset "PlotlyJS" begin
    # ------------------------------------------------------------------------------------ #
    # Wigner-Seitz cells and k-space path visualization
    # 3D
    sgnum = 147
    Rs = convert(SVector{3,SVector{3,Float64}}, [[1.0,0,0], [-.5,√3/2,0], [0,0,1.25]])
    Gs = Bravais.reciprocalbasis(Rs)

    cell = wignerseitz(Gs)
    kp = irrfbz_path(sgnum, Rs)
    @test plot(cell) isa PlotlyJS.SyncPlot
    @test plot(kp) isa PlotlyJS.SyncPlot
    @test plot(cell, kp) isa PlotlyJS.SyncPlot

    # 2D
    sgnum = 5
    cntr = Bravais.centering(sgnum, 2) # c-type centered
    Rs = convert(SVector{2,SVector{2,Float64}}, [[1.0,0], [0,.8]])
    pGs = Bravais.primitivize(Bravais.reciprocalbasis(Rs), cntr)

    cell = wignerseitz(pGs)
    kp = irrfbz_path(sgnum, Rs)
    @test plot(cell) isa PlotlyJS.SyncPlot
    @test plot(kp) isa PlotlyJS.SyncPlot
    @test plot(cell, kp) isa PlotlyJS.SyncPlot

    # ------------------------------------------------------------------------------------ #
    # Band structure plots (example from docs)
    sgnum = 227
    Rs = [[1,0,0], [0,1,0], [0,0,1]]
    kp = irrfbz_path(sgnum, Rs)
    pGs = basis(kp)
    kpi = interpolate(kp, 100)

    function ϵ(k; γ::Real=1.0)
        kx = 2π*(-k[1]+k[2]+k[3])
        ky = 2π*(+k[1]-k[2]+k[3])
        kz = 2π*(+k[1]+k[2]-k[3])
        return 4γ * (cos(kx/2)*cos(ky/2) + cos(ky/2)*cos(kz/2) + cos(kz/2)*cos(kx/2))
    end

    band1 = ϵ.(kpi)
    band2 = 20 .- (1/2).*band1

    # basic functionality
    @test plot(kpi, [band1]) isa PlotlyJS.SyncPlot
    @test plot(kpi, reduce(hcat, [band1, band2])) isa PlotlyJS.SyncPlot # matrix input

    # annotations
    offset = 2maximum(band2)
    band3 = band1 .+ offset
    annotations = Dict(:X => [1:2 => "X, bands 1:2", 30.0 => "X, some frequency"],
                       :Γ => [1 => "Γ, band 1", 3 => "Γ, band 3"])
    @test plot(kpi, [band1, band2, band3]; annotations) isa PlotlyJS.SyncPlot

    # band highlights
    band_highlights = Dict(2:3 => attr(color=:black, width=3))
    @test plot(kpi, [band1, band2, band3]; band_highlights) isa PlotlyJS.SyncPlot

    # title/y-label
    @test plot(kpi, [band1]; title="hello world", ylabel="ωa/2πc") isa PlotlyJS.SyncPlot
end


text:
"""
    module FELinearAlgebra
Provides Finite Element linear algebra matrices
"""
module FELinearAlgebra

import LinearAlgebra
import Base: *, +, -, ==, \
import Base: getindex, setindex!, size

export
    # Types
    LocalMatrix,
    JacobianMatrix,
    SolutionVector,
    # Functions
    L2_error!,
    L2_norm,
    lp_error!,
    lp_norm,
    true_size

include("SolutionVector.jl")
include("LocalMatrix.jl")

end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
"""
This is Tutorial_Flegrei.jl. It contains all the necessary commands to import the data files from
geophysical models at the volcano and represent them in cartesian and UTM coordinates.
The following steps are performed:
1. Export earthquake data with a time label in both formats.
2. Export a velocity models in both formats.
3. Scale and interpolate a shear-wave model and export in both formats.

You will need to download the zipped folder containing all files from:
[https://seafile.rlp.net/f/ff2c8424274c4d56b1f7/](https://ngdc.noaa.gov/mgg/global/global.html)
Remember to work in the downloaded directory.

"""
using DelimitedFiles, GeophysicalModelGenerator, Glob, GeoStats

# 1. define where the two files are located on your computer
data_80s            = readdlm("SeismicLocations/Seismicity_UTM_1983_1984.txt", '\t', skipstart=0, header=false);
data_00s            = readdlm("SeismicLocations/Seismicity_UTM_2005_2016.txt", ' ', skipstart=0, header=false);

# 1.1 create a single file and detine the timing and locations
data                = vcat(data_80s, data_00s)
time                = data[:,1];
WE                  = data[:,2];
SN                  = data[:,3];
depth               = data[:,4];

# 1.2 save in both formats
EQ_Data_Cart        = CartData(WE, SN, depth, (Depth = depth * m, Time = time * yr,));
Write_Paraview(EQ_Data_Cart, "CF_Earthquakes_Cartesian", PointsData=true)
EQ_Data_UTM         = UTMData(WE, SN, depth, 33, true, (Depth = depth * m, Time = time * yr,));
Data_set_UTM        =   convert(GeophysicalModelGenerator.GeoData, EQ_Data_UTM)
Write_Paraview(Data_set_UTM, "CF_Earthquakes_UTM", PointsData=true)

# 2 load the velocity model and assign the locations of nodes and velocities + Vp/Vs
data            =   readdlm("TravelTimeTomography/modvPS.dat", '\t', Float64, skipstart=0, header=false);
WE              =   data[:,1];
SN              =   data[:,2];
depth           =   data[:,3];
Vp              =   data[:,4];
Vs              =   data[:,5];
VpVs            =   data[:,6];

# 2 Reshape everything and save in both formats.
resolution      =   (length(unique(depth)),  length(unique(SN)), length(unique(WE)))
dim_perm        =   [3 2 1]
we              =   permutedims(reshape(WE, resolution), dim_perm);
sn              =   permutedims(reshape(SN, resolution), dim_perm);
depth           =   permutedims(reshape(depth, resolution), dim_perm);
Vp3d            =   permutedims(reshape(Vp, resolution), dim_perm);
Vs3d            =   permutedims(reshape(Vs, resolution), dim_perm);
Vp_Vs3d         =   permutedims(reshape(VpVs, resolution), dim_perm);
Data_set_Cartesian  =   CartData(we, sn, depth, (vp = Vp3d * (km / s), vs = Vs3d * (km / s), vpvs = Vp_Vs3d,))
Write_Paraview(Data_set_Cartesian, "CF_Velocity_Cartesian")
Data_set        =   UTMData(we, sn, depth, 33, true, (vp = Vp3d * (km / s), vs = Vs3d * (km / s), vpvs = Vp_Vs3d,))
Data_set_UTM    =   convert(GeophysicalModelGenerator.GeoData, Data_set)
Write_Paraview(Data_set_UTM, "CF_Velocity_UTM")

# 3 Create a list of text files containing sections and loop through them
list_files        = glob("AmbientNoiseTomography/*.txt");
li                = size(list_files, 1);
for i = 1:li
    nameFile        = list_files[i];
    name_vts        = nameFile[24:26];
    data            = readdlm(nameFile, '\t', Float64);
    WE              = data[:,1];
    SN              = data[:,2];
    depth           = data[:,3];
    Vs              = data[:,4];

# 3.1 Constrain them on a smaller grid:
    ind             = findall((WE .>= 419000) .& (WE .<= 435000) .& (SN .>= 4514000) .& (SN .<= 4528000));
    WE              = WE[ind];
    SN              = SN[ind];
    depth           = depth[ind];
    Vs              = Vs[ind];

  # Then create a 3D lat,long regular grid from UTM - measurements are all in the same time zone:
    l               = length(WE);

    n_WE            = minimum(WE):100:maximum(WE);
    n_SN            = minimum(SN):100:maximum(SN);
    we, sn, Depth   = XYZGrid(n_WE, n_SN, depth[1]);
    Vs_3D           = zeros(size(Depth));
    Cgrid           = GeoStats.CartesianGrid((size(we, 1), size(we, 2)), (minimum(we), minimum(sn)), (we[2,2,1] - we[1,1,1], sn[2,2,1] - sn[1,1,1]))
    coord           = PointSet([WE[:]'; SN[:]']);
    Geo             = georef((Vs = Vs[:],), coord);
    P               = EstimationProblem(Geo, Cgrid, :Vs);
    S               = IDW(:Vs => (;neighbors=2));
    sol             = solve(P, S);
    sol_Vs          = values(sol).Vs;
    Vs_2D           = reshape(sol_Vs, size(domain(sol)));
    Vs_3D[:,:,1]    = Vs_2D;
    Data_set_Cart   = CartData(we, sn, Depth, (Vs = Vs_3D  * (km / s),))
    Write_Paraview(Data_set_Cart, "CF_Noise" * name_vts * "_Cartesian")
    Data_set        = UTMData(we, sn, Depth, 33, true, (Vs = Vs_3D * (km / s),));
    Data_set_UTM    = convert(GeophysicalModelGenerator.GeoData, Data_set);
    Write_Paraview(Data_set_UTM, "CF_Noise" * name_vts * "_UTM")
end


text:
# ---
# title: 1556. Thousand Separator
# id: problem1556
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/thousand-separator/description/>
# hidden: true
# ---
#
# Given an integer `n`, add a dot (".") as the thousands separator and return it
# in string format.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 987
#     Output: "987"
#
#
# **Example 2:**
#
#
#
#     Input: n = 1234
#     Output: "1.234"
#
#
# **Example 3:**
#
#
#
#     Input: n = 123456789
#     Output: "123.456.789"
#
#
# **Example 4:**
#
#
#
#     Input: n = 0
#     Output: "0"
#
#
#
#
# **Constraints:**
#
#   * `0 <= n < 2^31`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ITensors,
      Test

@test LatticeBond(1, 2) == LatticeBond(1, 2, 0.0, 0.0, 0.0, 0.0, "")
@testset "Square lattice" begin
    sL = squareLattice(3, 4)
    @test length(sL) == 17
end

@testset "Triangular lattice" begin
    tL = triangularLattice(3, 4)
    @test length(tL) == 23
end


text:
using MeshBasisFunctions
using Meshes
using Test

points = Point2[(0,0), (1,0), (0,1), (1,1),
                (0.25,0.5), (0.75,0.5), (0.5, 0), (0.5, 1.0)]
tris = connect.([(1,5,3), (1, 7, 5), (7, 6, 5), (7, 2, 6),
                 (4, 6, 2), (4, 8, 6), (8, 5, 6), (8, 3, 5)],
                Triangle)
mesh = SimpleMesh(points, tris)
md = meshdata(mesh, Dict(0 => (z = randn(8), )))

basis = PiecewiseLinearBasis{2}()

interp_pts = [Point(0.0, 0.0),
              Point(0.5, 0.5),
              Point(0.34, 0.98),
              Point(0.12, 0.84)]
interp_op = interpolation_operator(interp_pts,
                                   mesh,
                                   basis)
interp1 = interp_op * values(md, 0)[:z]
interp2 = [interpolate(pt, md, :z, basis) for pt in interp_pts]

@testset "MeshBasisFunctions.jl" begin
    @test all(interp1 .== interp2)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

using .Convex
using .ECOS

function _kliep_coeffs(K_nu, K_de, dre::KLIEP, optlib::Type{ConvexLib})
  # retrieve parameters
  σ, b = dre.σ, size(K_de, 2)

  # number of numerator and denominator samples
  n_nu, n_de = size(K_nu, 1), size(K_de, 1)

  # constants for objective and constraints
  K = K_nu
  k = sum(K_de, dims=1)

  # objective function and constraints
  α = Convex.Variable(b); w = K*α
  objective   = sum(log(w[i]) for i in 1:n_nu)
  constraints = [α ≥ 0, dot(α,k) == n_de]

  # optimization problem
  problem = Convex.maximize(objective, constraints)

  # solve problem with ECOS solver
  Convex.solve!(problem, ECOS.Optimizer(verbose=false))

  # optimal coefficients
  vec(α.value)
end


text:
using FractionalDiffEq
using Plots, LaTeXStrings

s="\$ODE\\ Example\$"

T = 30
h = 0.05
tspan = collect(0.05:h:T)

f(x)=1/2*(-exp(-x)-sin(x)-cos(x)+2)
target=f.(tspan)

rightfun(x) = sin(x)
prob = MultiTermsFODEProblem([1, 1], [2, 1], rightfun, [0, 0], T)
sol = solve(prob, h, FODEMatrixDiscrete())

plot(sol, title=s, legend=:bottomright, label="ODE Numerical Solution!")

plot!(tspan, target, lw=3,ls=:dash,label="ODE Analytical Solution!")


text:
geneCols = Dict{String,Symbol}(
   "aliases" => :GeneSymbolAlias,
   "symbol" => :OfficialGeneSymbol,
   "region" => :GeneRegion,
   "regionStart" => :GeneStart,
   "regionEnd" => :GeneEnd,
   "id" => :EnsemblGeneId,
   "version" => :EnsemblGeneVersion,
   "biotype" => :EnsemblGeneBiotype,
   "name" => :GeneName
);
transcriptCols = Dict{String,Symbol}(
   "region" => :TranscriptRegion,
   "regionStart" => :TranscriptStart,
   "regionEnd" => :TranscriptEnd,
   "id" => :EnsemblTranscriptId,
   "version" => :EnsemblTranscriptVersion,
   "biotype" => :EnsemblTranscriptBiotype
);
proteinCols = Dict{String,Symbol}(
   "id" => :EnsemblProteinId,
   "version" => :EnsemblProteinVersion
);
elTypeDict = Dict{String, Union}(
   # "aliases" => Vector{String},
   "aliases" => Union{String, Missing},
   "symbol" => Union{String, Missing},
   "region" => Union{String, Missing},
   "regionStart" => Union{Int, Missing},
   "regionEnd" => Union{Int, Missing},
   "id" => Union{String, Missing},
   "version" => Union{Int, Missing},
   "biotype" => Union{String, Missing},
   "name" => Union{String, Missing}
);

function parseGngJson(json::Array{Any,1})

   #########################
   # Count number or rows and create columns
   #  In the case of the columns only look at the first entry
   #  First four columns are inputDB, inputID, targetID and targetDB
   nCols = 4;
   colNames = Symbol["InputId", "InputSourceDb", "TargetId", "TargetDb"];
   colNamesUnstack = Symbol["InputId", "InputSourceDb"];
   colElTypes = Union[Union{String, Missing}, Union{String, Missing}, Union{String, Missing}, Union{String, Missing}];
   # Gene column names and types
   map(i->haskey(geneCols,i) && (push!(colNames, geneCols[i]),
                                 push!(colNamesUnstack, geneCols[i]),
                                 push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]))
   # Transcript:
   # Dict{String,Any} with 4 entries:
   # "peptide"   => Dict{String,Any}(Pair{String,Any}("id", "ENSP00000381185"),Pair{String,Any}("targetIds", Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]),Pair{String,A…
   # "id"        => "ENST00000398117"
   # "targetIds" => Any[]
   # "region"    => "18"
   map(i->haskey(transcriptCols,i) && (push!(colNames, transcriptCols[i]),
                                       push!(colNamesUnstack, transcriptCols[i]),
                                       push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]["transcript"][1]))
   # number of cols minus 1 because of targetIds, e.g.:
   # Dict{String,Any} with 3 entries:
   # "id"        => "ENSP00000381185"
   # "targetIds" => Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]
   # "version"   => 1
   map(i->haskey(proteinCols,i) && (push!(colNames, proteinCols[i]),
                                    push!(colNamesUnstack, proteinCols[i]),
                                    push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]["transcript"][1]["peptide"]))
   ##################
   # Rows
   nRows = 0;
   for entry in json
      # For each gene entry
      # nRows += 1;
      # Number of external gene identifiers
      gene = entry["row"][1]["gene"];
      nRows += (length(gene["targetIds"]) * length(gene["transcript"]));
      for transcript in gene["transcript"]
         # Transcript:
         # Dict{String,Any} with 4 entries:
         # "peptide"   => Dict{String,Any}(Pair{String,Any}("id", "ENSP00000381185"),Pair{String,Any}("targetIds", Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]),Pair{String,A…
         # "id"        => "ENST00000398117"
         # "targetIds" => Any[]
         # "region"    => "18"
         nRows += length(transcript["targetIds"]);
         # Peptide:
         # Dict{String,Any} with 3 entries:
         # "id"        => "ENSP00000381185"
         # "targetIds" => Any[Dict{String,Any}(Pair{String,Any}("dbName", "UniProtKB/TrEMBL"),Pair{String,Any}("id", "A0A024R2B3"))]
         # "version"   => 1
         nRows += length(transcript["peptide"]["targetIds"]);
      end
      # println(entry["row"][1])
   end
   nRows

   x = DataFrames.DataFrame(
      colElTypes,
      colNames,
      nRows
      # 100
   );

   lowG = 1;
   highG = 1;
   lowT = 1;
   highT = 1;
   for entry in json
      # For each gene entry
      lowT = lowG;
      highT = lowT;
      highG = lowG;
      # Number of external gene identifiers
      gene = entry["row"][1]["gene"];
      externalIds = String[];
      externalDbs = String[];
      # In the case of external IDs add entries to arrays
      !isempty(gene["targetIds"]) && map(i->(push!(externalIds, i["id"]), push!(externalDbs, i["dbName"])), gene["targetIds"]);
      externalGeneIds = copy(externalIds);
      externalGeneDbs = copy(externalDbs);
      # nRows += length(gene["targetIds"])
      for transcript in gene["transcript"]
         # In the case of external IDs add entries to arrays
         !isempty(transcript["targetIds"]) && map(i->(push!(externalIds, i["id"]), push!(externalDbs, i["dbName"])), transcript["targetIds"]);
         peptide = transcript["peptide"];
         !isempty(peptide["targetIds"]) && map(i->(push!(externalIds, i["id"]), push!(externalDbs, i["dbName"])), peptide["targetIds"]);
         highT += (length(externalIds) - 1);
         # for
         # map(i->haskey(geneCols,i) && (push!(colNames, geneCols[i]), push!(colElTypes, elTypeDict[i])), keys(json[1]["row"][1]["gene"]))
         # for epEntry in transcript["peptide"]["targetIds"]

         map(key->haskey(transcriptCols,key) && (x[transcriptCols[key]][lowT:highT] = transcript[key]), keys(transcript))
         map(key->haskey(proteinCols,key) && (x[proteinCols[key]][lowT:highT] = peptide[key]), keys(peptide))

         # map(key->haskey(geneCols,key) && (println("col: $(geneCols[key]) = $(gene[key])")), keys(gene))
         # map(key->haskey(geneCols,key) && (println("$(collect(lowT:highT)) -> $(fill(gene[key], high))")), keys(gene))

         x[:TargetId][lowT:highT] = externalIds;
         x[:TargetDb][lowT:highT] = externalDbs;
         externalIds = copy(externalGeneIds);
         externalDbs = copy(externalGeneDbs);
         # x[low,:EnsemblGeneId] = gene["id"];
         # x[low,:EnsemblTranscriptId] = transcript["id"];
         # x[low,:EnsemblProteinId] = peptide["id"];
         # println("[$i] : $externalId, $externalDb")
         # println.("enspId: $(transcript["peptide"]["id"]) -> $(transcript["peptide"]["targetIds"])")
         lowT = highT + 1;
         highT = lowT;
      end
      highG = highT-1;
      # Fill in gene information
      for (key, value) in gene
         if haskey(geneCols, key) && !isempty(value)
            val = value;
            if key == "aliases"
               val = join(value, ",");
            end
            x[geneCols[key]][lowG:highG] = val;
         end
      end
      # map(key->haskey(geneCols,key) && (x[geneCols[key]][lowG:highG] = fill(join(gene[key], ","), (highG - lowG + 1))), keys(gene));
      x[:InputId][lowG:highG] = entry["row"][1]["inputId"];
      x[:InputSourceDb][lowG:highG] = entry["row"][1]["inputDb"];
      lowG = highG + 1;
   end
   x

   return x;

end

function unstackDf(x::D, keyCols::S = Vector{Symbol}(), keyColsElType::T = Vector{DataType}()) where {
   D<:DataFrames.DataFrame,
   S<:Vector{Symbol},
   T<:Vector{DataType}}

   if isempty(keyCols)
      keyCols = names(x)[map(i->!in(i,["TargetId", "TargetDb"]), collect(String,names(x)))];
      keyColsElType = DataFrames.eltypes(x)[map(i->!in(i,["TargetId", "TargetDb"]), collect(String,names(x)))];
   elseif length(keyCols != length(keyColsElType))
      error("keyCols and keyColsElType have to be of same length.!");
   end
   colNames = copy(keyCols);
   colElTypes = copy(keyColsElType);

   append!(colNames, DataFrames.identifier.(collect(unique(x[:TargetDb]))));
   append!(colElTypes, fill(Union{String, Missing}, (length(colNames) - length(keyCols))));
   nRows = DataFrames.nrow(unique(x[keyCols]));

   xDf = DataFrames.DataFrame(
      colElTypes,
      colNames,
      nRows
   );
   sortByCols = copy(keyCols);
   push!(sortByCols, :TargetDb);
   sort!(x, cols = sortByCols);
   idx = 0;
   key = "";
   for row in DataFrames.eachrow(x)
      keyCurr = join(map(i->i[2], row[keyCols]));
      if key != keyCurr
         idx += 1;
         key = keyCurr;
         map(i->xDf[idx,i[1]] = i[2], row[keyCols]);
      end
      colId = DataFrames.identifier(row[:TargetDb])
      xDf[idx, colId] = DataFrames.ismissing(xDf[idx, colId]) ?
                           row[:TargetId] :
                           xDf[idx, colId] * "," * row[:TargetId];
   end
   return xDf;
end

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
@testset "NonOverlappingOrthotopes2DisjunctionLP: $(fct_type), orthotopes $(n_ortho), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], n_ortho in [2, 3], dim in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(DisjunctionModel{T}())
    model = COIB.NonOverlappingOrthotopes2DisjunctionLP{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    if fct_type == "vector of variables"
        @test MOI.supports_constraint(
            model,
            MOI.VectorOfVariables,
            CP.Disjunction{NTuple{4, MOI.LessThan{T}}},
        )
        @test MOI.supports_constraint(
            model,
            MOI.VectorOfVariables,
            CP.Disjunction{NTuple{6, MOI.LessThan{T}}},
        )
        @test MOIB.supports_bridging_constraint(
            model,
            MOI.VectorOfVariables,
            CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES},
        )
    elseif fct_type == "vector affine function"
        @test MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            CP.Disjunction{NTuple{4, MOI.LessThan{T}}},
        )
        @test MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            CP.Disjunction{NTuple{6, MOI.LessThan{T}}},
        )
        @test MOIB.supports_bridging_constraint(
            model,
            MOI.VectorAffineFunction{T},
            CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES},
        )
    else
        @assert false
    end

    if T == Int
        x_pos, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(dim * n_ortho)])
        x_sze, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(dim * n_ortho)])
        x_end, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:(dim * n_ortho)])
    elseif T == Float64
        x_pos = MOI.add_variables(model, dim * n_ortho)
        x_sze = MOI.add_variables(model, dim * n_ortho)
        x_end = MOI.add_variables(model, dim * n_ortho)
    end
    x_ortho = [
        [
            x_pos[(1 + (i - 1) * dim) : (i * dim)],
            x_sze[(1 + (i - 1) * dim) : (i * dim)],
            x_end[(1 + (i - 1) * dim) : (i * dim)],
        ]
        for i in 1:n_ortho
    ]
    x = vcat([vcat(x_ortho[i]...) for i in 1:n_ortho]...)

    # Test the construction of x before going on (no loops to ease checking
    # for correctness).
    @test length(x_pos) == dim * n_ortho
    @test length(x_sze) == dim * n_ortho
    @test length(x_end) == dim * n_ortho
    @test length(x) == 3 * dim * n_ortho
    if n_ortho >= 1
        @test x[(0 * dim + 1):(1 * dim)] == x_pos[(0 * dim + 1):(1 * dim)]
        @test x[(1 * dim + 1):(2 * dim)] == x_sze[(0 * dim + 1):(1 * dim)]
        @test x[(2 * dim + 1):(3 * dim)] == x_end[(0 * dim + 1):(1 * dim)]
    end
    if n_ortho >= 2
        @test x[(3 * dim + 1):(4 * dim)] == x_pos[(1 * dim + 1):(2 * dim)]
        @test x[(4 * dim + 1):(5 * dim)] == x_sze[(1 * dim + 1):(2 * dim)]
        @test x[(5 * dim + 1):(6 * dim)] == x_end[(1 * dim + 1):(2 * dim)]
    end
    if n_ortho >= 3
        @test x[(6 * dim + 1):(7 * dim)] == x_pos[(2 * dim + 1):(3 * dim)]
        @test x[(7 * dim + 1):(8 * dim)] == x_sze[(2 * dim + 1):(3 * dim)]
        @test x[(8 * dim + 1):(9 * dim)] == x_end[(2 * dim + 1):(3 * dim)]
    end
    if n_ortho >= 4
        @assert false
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end

    c = MOI.add_constraint(model, fct, CP.NonOverlappingOrthotopes(n_ortho, dim))

    for i in 1:(3 * dim *  n_ortho)
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.NonOverlappingOrthotopes{CP.UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
            (MOI.VectorAffineFunction{T}, CP.Disjunction{NTuple{n, MOI.LessThan{T}} where n}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == n_ortho * dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Disjunction{NTuple{n, MOI.LessThan{T}} where n}}()) == Int(n_ortho * (n_ortho - 1) / 2)

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == collect(values(bridge.cons_ends))
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Disjunction{NTuple{n, MOI.LessThan{T}} where n}}()) == collect(values(bridge.cons_disjunction))
    end

    @testset "End-point constraints" begin
        @test length(bridge.cons_ends) == n_ortho * dim
        for i in 1:n_ortho
            for d in 1:dim
                @test MOI.is_valid(model, bridge.cons_ends[i, d])
                @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_ends[i, d]) == MOI.EqualTo(zero(T))
                f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_ends[i, d])

                @test length(f.terms) == 3

                t1 = f.terms[1]
                @test t1.coefficient === one(T)
                @test t1.variable === x_pos[(i - 1) * dim + d]

                t2 = f.terms[2]
                @test t2.coefficient === one(T)
                @test t2.variable === x_sze[(i - 1) * dim + d]

                t3 = f.terms[3]
                @test t3.coefficient === -one(T)
                @test t3.variable === x_end[(i - 1) * dim + d]
            end
        end
    end

    @testset "Disjunction constraints" begin
        @test length(bridge.cons_disjunction) == Int(n_ortho * (n_ortho - 1) / 2)
        for i in 1:n_ortho
            for j in 1:n_ortho
                if i < j
                    @test MOI.is_valid(model, bridge.cons_disjunction[i, j])
                    @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_disjunction[i, j]) isa CP.Disjunction

                    f_vector = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_disjunction[i, j])
                    @test length(f_vector.terms) == 3 * 2 * dim

                    f_scalars = MOIU.scalarize(f_vector)
                    @test length(f_scalars) == 2 * dim

                    for d in 1:dim
                        f1 = f_scalars[(d - 1) * 2 + 1]
                        f2 = f_scalars[(d - 1) * 2 + 2]

                        # Doesn't work, because == requires that the order
                        # of the terms is the same.
                        # @test f1 == one(T) * x_pos[(i - 1) * dim + d] + one(T) * x_sze[(i - 1) * dim + d] - one(T) * x_pos[(j - 1) * dim + d]

                        t = f1.terms[1]
                        @test t.coefficient === one(T)
                        @test t.variable === x_pos[(i - 1) * dim + d]

                        t = f1.terms[2]
                        @test t.coefficient === -one(T)
                        @test t.variable === x_pos[(j - 1) * dim + d]

                        t = f1.terms[3]
                        @test t.coefficient === one(T)
                        @test t.variable === x_sze[(i - 1) * dim + d]

                        t = f2.terms[1]
                        @test t.coefficient === -one(T)
                        @test t.variable === x_pos[(i - 1) * dim + d]

                        t = f2.terms[2]
                        @test t.coefficient === one(T)
                        @test t.variable === x_pos[(j - 1) * dim + d]

                        t = f2.terms[3]
                        @test t.coefficient === one(T)
                        @test t.variable === x_sze[(j - 1) * dim + d]
                    end
                end
            end
        end
    end
end


text:
using IcanProj
using PyPlot
using DataFrames
using PyCall
using JFSM2
using CSV
using Statistics

import IcanProj.project_results

function project_results(df::DataFrame, on::Symbol)

    map_values = fill(NaN, 1550, 1195)

    map_values[df[:ind_julia]] = df[on]

    return map_values

end

function plot_map(df, variable, cb_label, cb_unit, title_label, figpath)

    resmap = project_results(df, variable)

    file = joinpath(figpath, String(variable) * ".png")

    cb_text =  cb_label * " " * cb_unit

    py"""

    import matplotlib.pyplot as plt

    fig, ax = plt.subplots()

    #fig.figsize = (6, 12)

    ax.set_title($(title_label))

    img_plot = ax.imshow($(resmap), cmap = 'jet', vmin = -0.1, vmax = 1.1)

    cbar = fig.colorbar(img_plot)

    cbar.set_label($(cb_text))

    #cbar.set_clim(0, 1)

    ax.grid(linestyle='dotted')

    plt.tick_params(axis='both', left=False, top=False, right=False, bottom=False,
                    labelleft=False, labeltop=False, labelright=False, labelbottom=False)

    plt.tight_layout()

    fig.savefig($(file), dpi = 200, bbox_inches='tight')

    plt.close('all')

    """

    @info "finished plotting"

end

# Load results

path = dirname(pathof(IcanProj))

df_res = CSV.File(joinpath(path, "..", "data", "table_errors_50km.txt")) |> DataFrame

df_links = CSV.read(joinpath(path, "..", "data", "df_links.csv")) |> DataFrame

df_all = join(df_links, df_res, on = :ind_50km)

cfgs = cfg_table()

# Plot average error for all cfgs with exchng = 0

for variable in ["swe", "latmo", "hatmo", "melt"]

    cfg_subset = cfgs[convert(Array{Bool}, cfgs[:exchng] .== 0), :]

    colnames = Symbol.("nse_$(variable)_cfg" .* string.(cfg_subset[:cfg]))

    df_all[Symbol("nse_$(variable)_exchng0")] =  mean(convert(Array, df_all[colnames]), dims = 2)[:]

end

# Plot average error for all cfgs with exchng = 1

for variable in ["swe", "latmo", "hatmo", "melt"]

    cfg_subset = cfgs[convert(Array{Bool}, cfgs[:exchng] .== 1), :]

    colnames = Symbol.("nse_$(variable)_cfg" .* string.(cfg_subset[:cfg]))

    df_all[Symbol("nse_$(variable)_exchng1")] =  mean(convert(Array, df_all[colnames]), dims = 2)[:]

end

# Plot averaged results

info = [(:nse_swe_exchng0,   "NSE", "(-)", "SWE exchnge = 0"),
        (:nse_latmo_exchng0, "NSE", "(-)", "Latmo exchnge = 0"),
        (:nse_hatmo_exchng0, "NSE", "(-)", "Hatmo exchnge = 0"),
        (:nse_melt_exchng0,  "NSE", "(-)", "Melt exchnge = 0"),
        (:nse_swe_exchng1,   "NSE", "(-)", "SWE exchnge = 1"),
        (:nse_latmo_exchng1, "NSE", "(-)", "Latmo exchnge = 1"),
        (:nse_hatmo_exchng1, "NSE", "(-)", "Hatmo exchnge = 1"),
        (:nse_melt_exchng1,  "NSE", "(-)", "Melt exchnge = 1")]

figpath = joinpath(path, "..", "plots", "error_maps_1")

for (v, cl, cu, t) in info

    plot_map(df_all, v, cl, cu, t, figpath)

end

#variable = :nse_rnet_cfg1

#cb_label = "Correlation coeffcient"

#cb_unit = "(-)"

#title_label = "SWE"

#figpath = "dummy"

#plot_map(df_all, variable, cb_label, cb_unit, title_label, figpath)

#mapres = project_results(df, :r_swe_cfg1)

#figure()
#imshow(mapres)


text:
using Statistics
x = rand(10^7);
data = 1 ./ (x .^ 2);
a,b,c,d = Statistics.quantile!(data, [0, .001, .999, 1]);
#used_bits(d) # Naive required bit-width
println("min=",a,", 0.1%=", b, ", 99.9%=",c,", max=",d)


text:
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)
@testset "test ac polar opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5907; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 17551; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 194.14; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 397.36; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

    @test isapprox(result["objective"], 150228.15; atol = 1e0)
    end
end

@testset "test dc opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5638.97; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 17479; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 171.69; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 371.01; atol = 1e0)
    end

    @testset "24-bus rts ac dc case with three zones" begin
       result = run_acdcopf("../test/data/case24_3zones_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

       @test result["termination_status"] == LOCALLY_SOLVED
       @test isapprox(result["objective"], 143714.75; atol = 1e0)
    end
end

@testset "test soc (BIM) opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 14999; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf("../test/data/case5_2grids.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 363.50; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", SOCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.24; atol = 1e0)
    end
end

@testset "test soc distflow opf_bf" begin
    @testset "3-bus case" begin
        result = run_acdcopf_bf("../test/data/case3.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf_bf("../test/data/case5_asym.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 14999; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf_bf("../test/data/case5_acdc.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.91; atol = 1e0)
    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcopf_bf("../test/data/case5_2grids.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 363.50; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf_bf("../test/data/case24_3zones_acdc.m", SOCBFPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.26; atol = 1e0)
    end
end

@testset "test qc opf" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5780; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_acdcopf("../test/data/case5_asym.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 15921; atol = 1e0)
    end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", QCWRPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 150156.25; atol = 1e0)
    end
end

@testset "test qc opf with trilinear convexhull relaxation" begin
    @testset "3-bus case" begin
        result = run_acdcopf("../test/data/case3.m", QCWRTriPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 5817.58; atol = 1e0)
    end
    # @testset "5-bus asymmetric case" begin
    #     result = run_acdcopf("../test/data/case5_asym.m", QCWRTriPowerModel, ipopt_solver)
    #
    #     @test result["termination_status"] == LOCALLY_SOLVED
    #     @test isapprox(result["objective"], 15816.9; atol = 1e0)
    # end
    @testset "5-bus ac dc case" begin
        result = run_acdcopf("../test/data/case5_acdc.m", QCWRTriPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 183.76; atol = 1e0)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcopf("../test/data/case24_3zones_acdc.m", QCWRTriPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"],  150156.25; atol = 1e0)
    end
end

# @testset "test sdp opf" begin
#     @testset "3-bus case" begin
#         result = run_acdcopf("../test/data/case3.m", SDPWRMPowerModel, scs_solver; setting = s)
#
#         @test result["termination_status"] == :Optimal
#         @test isapprox(result["objective"], 5851.3; atol = 2e0)
#     end
#     # @testset "5-bus ac dc case" begin
#     #     result = run_acdcopf("../test/data/case5_acdc.m", SDPWRMPowerModel, mosek)
#     #
#     #     @test result["termination_status"] == :Optimal
#     #     @test isapprox(result["objective"], 194.12; atol = 1e0)
#     # end
#     # TODO replace this with smaller case, way too slow for unit testing
#     # @testset "24-bus rts ac dc case with three zones" begin
#     #    result = run_acdcopf("../test/data/case24_3zones_acdc.m", SDPWRMPowerModel, scs_solver)
#     #
#     #    @test result["termination_status"] == :Optimal
#     #    @test isapprox(result["objective"], 150179.96; atol = 1e0)
#     # end
# end


text:
"""
    info(registry)

Show information about a registry, specifically its name, description
and fields.

## Examples

{cell}
```julia
using FeatureRegistries: exampleregistry, info
registry = exampleregistry()
info(registry)
```
"""
function info end

struct RegistryInfo
    registry::Registry
end

function Base.show(io::IO, registryinfo::RegistryInfo)
    registry = registryinfo.registry
    println(io,
        "Information on registry ",
        crayon"underline", getfield(registry, :name), crayon"reset", "\n\n",
        "$(crayon"bold")Fields$(crayon"reset"): ")

    fs = registry.fields
    data = reduce(vcat, [
            reshape([
                f.name,
                AnsiTextCell(k in getfield(registry, :id) ? "$(crayon"red bold"):$k (ID)" : ":$k"),
                AnsiTextCell("$(crayon"dark_gray")$(_fieldtype(f))"),
                f.description],
            1, :)
            for (k, f) in pairs(fs)
        ])
    pretty_table(
        io,
        data,
        alignment=:l,
        vlines = :all,
        hlines = :none,
        vcrop_mode=:middle,
        tf=PrettyTables.tf_borderless,
        header = ["Name", "Column", "Type", "Description"]
    )

    if !isempty(registry.description)
        println(io, "\n$(crayon"bold")Description$(crayon"reset"): ")
        show(io, RichCell(Markdown.parse(registry.description)))
    end
end

_fieldtype(::Field{T}) where T = T
info(registry::Registry) = RegistryInfo(registry)


text:
module DocTestPkg

push!(LOAD_PATH, joinpath(Pkg.dir("DocTestPkg"), "src"))
using SupportMod

# DocTestPkg
export Equation, func, calcIt

"""
  Type holding all data for an instance of an equation
"""
type Equation{Tsol, Tres}
  q::Array{Tsol, 3}
  res::Array{Tres, 3}
  t::Float64
end

"""
  Outer constructor for [`Equation`](@ref)
"""
function Equation(numDofPerNode::Integer, numNodesPerElement::Integer,
                 numEl::Integer)
  q = zeros(numDofPerNode, numNodesPerElement, numEl)
  res = zeros(numDofPerNode, numNodesPerElement, numEl)
  t = 0.0

  return Equation{Float64, Float64}(q, res, t)
end

"""
Returns double the number `x` plus `1`
"""
func(y) = 2y + 1

"""
  This function calculates something
"""
function calcIt(x::Integer)
  return x + 2
end

end # module


text:
using QuantumStateBase
using Test
using LinearAlgebra
using ClassicalOrthogonalPolynomials

const QSB = QuantumStateBase
const DIM = 100

@testset "QuantumStateBase.jl" begin
    # polynomial
    include("polynomial.jl")

    # state
    include("operator.jl")
    include("state.jl")

    # wigner
    include("wigner_util.jl")
    include("wigner.jl")
end


text:
using ExaTron
using LinearAlgebra
using SparseArrays
using Random

using Profile
using PProf

Random.seed!(0)
n = 8
A = ExaTron.TronDenseMatrix{Array}(n)
L = ExaTron.TronDenseMatrix{Array}(n)
for j=1:n, i=j:n
    v = rand(1)[1]
    L.vals[i,j] = v
end
A.vals .= (L.vals * transpose(L.vals))

g = ones(n)
x = zeros(n)
xc = zeros(n)
xl = -Inf*ones(n)
xu = Inf*ones(n)
eval_f_cb(x) = 0.5*(transpose(x)*A.vals*x) + transpose(g)*x
function eval_grad_f_cb(x, grad_f)
    grad_f .= A.vals*x .+ g
end
function eval_h_cb(x, mode, rows, cols, scale, lambda, values)
    if mode == :Structure
        nz = 1
        for j=1:n,i=j:n
            rows[nz] = i
            cols[nz] = j
            nz += 1
        end
    else
        nz = 1
        for j=1:n,i=j:n
            values[nz] = A.vals[i,j]
            nz += 1
        end
    end
end

prob = ExaTron.createProblem(n, xl, xu, Int(floor((n*(n+1))/2)),
                             eval_f_cb, eval_grad_f_cb, eval_h_cb)#; :matrix_type=>:Dense)

howmany = 10
for i=1:howmany
    fill!(prob.x, 0.0)
    @profile ExaTron.solveProblem(prob)
end
pprof()


text:
using NPZ, Plots, NipponColors
cd(@__DIR__)

D = dict_color()

# sampler
gmv = npzread("generated_macroscopic_variables_normalized.npy")
∇gmv = npzread("generated_macroscopic_variables_normalized_gradients.npy")

# solver generator
smv = npzread("sod_macroscopic_variables_normalized.npy")
∇smv = npzread("sod_macroscopic_variables_gradients.npy")

# plot
p1 = scatter(gmv[:, 1], gmv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-1, 1), ylims=(0, 1), xlabel="U", ylabel="T")
savefig(p1, "illu_generator_u.pdf")

p2 = scatter(smv[:, 1], smv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-1, 1), ylims=(0, 1), xlabel="U", ylabel="T")
savefig(p2, "illu_solver_u.pdf")

p3 = scatter(∇gmv[:, 1], ∇gmv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-10, 10), ylims=(-10, 10), xlabel="∇U", ylabel="∇T")
savefig(p3, "illu_generator_du.pdf")

p4 = scatter(∇smv[:, 1], ∇smv[:, 2], marker = (:circle, 3, 0.6, D["ukon"], stroke(1, 0.1, :gray27, :dot)), legend=:none, xlims=(-10, 10), ylims=(-10, 10), xlabel="∇U", ylabel="∇T")
savefig(p4, "illu_solver_du.pdf")

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
module Surrogates
using LinearAlgebra
using Distributions

abstract type AbstractSurrogate <: Function end
include("utils.jl")
include("Radials.jl")
include("Kriging.jl")
include("Sampling.jl")
include("Optimization.jl")
include("Lobachesky.jl")
include("LinearSurrogate.jl")
include("InverseDistanceSurrogate.jl")
include("SecondOrderPolynomialSurrogate.jl")
include("SthenoKriging.jl")
include("RandomForestSurrogate.jl")
include("NeuralSurrogate.jl")
include("Wendland.jl")
include("MOE.jl")
include("VariableFidelity.jl")
include("PolynomialChaos.jl")
include("Earth.jl")

current_surrogates = ["Kriging","LinearSurrogate","LobacheskySurrogate","NeuralSurrogate",
                          "RadialBasis","RandomForestSurrogate","SecondOrderPolynomialSurrogate","Wendland"]
export current_surrogates_MOE
export AbstractSurrogate, SamplingAlgorithm
export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial
export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, LowDiscrepancySample
export RandomSample, KroneckerSample, GoldenSample
export SRBF,LCBS,EI,DYCORS,SOP,surrogate_optimize
export LobacheskySurrogate, lobachesky_integral, lobachesky_integrate_dimension
export LinearSurrogate
export RandomForestSurrogate
export SVMSurrogate
export NeuralSurrogate
export InverseDistanceSurrogate
export SecondOrderPolynomialSurrogate
export SthenoKriging
export Wendland
export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure
export LobacheskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure
export WendlandStructure
export MOE
export VariableFidelitySurrogate
export PolynomialChaosSurrogate
export EarthSurrogate

end


text:
module ScalarWrapper

using JavaCall

const BaseMultiStartMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultiStartMultivariateOptimizer
const BaseMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultivariateOptimizer
const BaseOptimizer = @jimport org.hipparchus.optim.BaseOptimizer
const Class = @jimport java.lang.Class
const ConvergenceChecker = @jimport org.hipparchus.optim.ConvergenceChecker
const Enum = @jimport java.lang.Enum
const GoalType = @jimport org.hipparchus.optim.nonlinear.scalar.GoalType
const GradientMultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.GradientMultivariateOptimizer
const JString = @jimport java.lang.String
const LeastSquaresConverter = @jimport org.hipparchus.optim.nonlinear.scalar.LeastSquaresConverter
const LineSearch = @jimport org.hipparchus.optim.nonlinear.scalar.LineSearch
const MultiStartMultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.MultiStartMultivariateOptimizer
const MultivariateFunction = @jimport org.hipparchus.analysis.MultivariateFunction
const MultivariateFunctionMappingAdapter = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
const MultivariateFunctionPenaltyAdapter = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
const MultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateOptimizer
const MultivariateVectorFunction = @jimport org.hipparchus.analysis.MultivariateVectorFunction
const Object = @jimport java.lang.Object
const ObjectiveFunction = @jimport org.hipparchus.optim.nonlinear.scalar.ObjectiveFunction
const ObjectiveFunctionGradient = @jimport org.hipparchus.optim.nonlinear.scalar.ObjectiveFunctionGradient
const OptimizationData = @jimport org.hipparchus.optim.OptimizationData
const Optional = @jimport java.util.Optional
const PointValuePair = @jimport org.hipparchus.optim.PointValuePair
const RandomVectorGenerator = @jimport org.hipparchus.random.RandomVectorGenerator
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const UnivariatePointValuePair = @jimport org.hipparchus.optim.univariate.UnivariatePointValuePair

export GoalType
export GradientMultivariateOptimizer
export LeastSquaresConverter
export LineSearch
export MultiStartMultivariateOptimizer
export MultivariateFunctionMappingAdapter
export MultivariateFunctionPenaltyAdapter
export MultivariateOptimizer
export ObjectiveFunction
export ObjectiveFunctionGradient
export bounded_to_unbounded
export compare_to
export compute_objective_value
export describe_constable
export equals
export get_class
export get_convergence_checker
export get_declaring_class
export get_evaluations
export get_goal_type
export get_iterations
export get_lower_bound
export get_max_evaluations
export get_max_iterations
export get_objective_function
export get_objective_function_gradient
export get_optima
export get_start_point
export get_upper_bound
export hash_code
export name
export notify
export notify_all
export optimize
export ordinal
export search
export to_string
export unbounded_to_bounded
export value
export value_of
export values
export wait

include("GradientWrapper/GradientWrapper.jl")
include("NoderivWrapper/NoderivWrapper.jl")
include("goal_type.jl")
include("gradient_multivariate_optimizer.jl")
include("least_squares_converter.jl")
include("line_search.jl")
include("multi_start_multivariate_optimizer.jl")
include("multivariate_function_mapping_adapter.jl")
include("multivariate_function_penalty_adapter.jl")
include("multivariate_optimizer.jl")
include("objective_function.jl")
include("objective_function_gradient.jl")

end


text:
module NOVAS

using DataDeps, LinearAlgebra, Memoize, StaticArrays, DelimitedFiles, FortranFiles,
      Polynomials
import Base.show

include("data.jl")
include("constants.jl")
include("utils.jl")
include("nutation.jl")
include("daf.jl")
include("naif_codes.jl")
include("spk.jl")
include("solarsystem.jl")
include("main.jl")
# Thin C wrapper
include("LibNOVAS.jl")

end