text:
using IICGP
using Cambrian
using CartesianGeneticProgramming
using Test

import Cambrian.mutate # function extension

function find_chr_index(chr::Vector{Float64}, v::Vector{Vector{Float64}})
    for i in eachindex(v)
        if v[i] == chr
            return i
        end
    end
end

function get_elites(e::DualCGPGAEvo)
    elite_echr = Vector{Vector{Float64}}()
    elite_cchr = Vector{Vector{Float64}}()
    elite_pairs = Vector{Tuple{Int64, Int64}}()
    ci = CartesianIndices(e.elites_matrix)
    for i in eachindex(e.elites_matrix)
        if e.elites_matrix[i]
            row, col = ci[i][1], ci[i][2]
            echr = deepcopy(e.encoder_sympop[row].chromosome)
            cchr = deepcopy(e.controller_sympop[col].chromosome)
            if !(echr in elite_echr)
                push!(elite_echr, echr)
            end
            if !(cchr in elite_cchr)
                push!(elite_cchr, cchr)
            end
            eindex = find_chr_index(echr, elite_echr)
            cindex = find_chr_index(cchr, elite_cchr)
            push!(elite_pairs, (eindex, cindex))
        end
    end
    elite_echr, elite_cchr, elite_pairs
end

function test_indexes(
    e::DualCGPGAEvo,
    elite_echr::Vector{Vector{Float64}},
    elite_cchr::Vector{Vector{Float64}},
    elite_pairs::Vector{Tuple{Int64, Int64}}
)
    indexes = select_indexes(e)
    nrows, ncols = size(e.fitness_matrix)
    # Test: elites are selected for next evaluation
    for pair in elite_pairs
        i, j = pair[1], pair[2]
        echr = elite_echr[i]
        cchr = elite_cchr[j]
        r = find_chr_index(echr, [ind.chromosome for ind in e.encoder_sympop])
        c = find_chr_index(cchr, [ind.chromosome for ind in e.controller_sympop])
        @test (r, c) ∈ indexes
    end
    # Test: each row and each col are to be evaluated
    isrowselected, iscolselected = falses(nrows), falses(ncols)
    for i in indexes
        isrowselected[i[1]] = true
        iscolselected[i[2]] = true
    end
    @test sum(isrowselected) == nrows
    @test sum(iscolselected) == ncols
end

function test_newpop(
    e::DualCGPGAEvo,
    elite_echr::Vector{Vector{Float64}},
    elite_cchr::Vector{Vector{Float64}},
    elite_pairs::Vector{Tuple{Int64, Int64}}
)
    @test length(e.encoder_sympop) == e.encoder_config.n_population
    @test length(e.controller_sympop) == e.controller_config.n_population
    @test sum(e.elites_matrix) == e.n_elite
    # Test: elites are kept between generations (same chromosomes)
    for i in eachindex(elite_echr)
        @test elite_echr[i] == e.encoder_sympop[i].chromosome
    end
    for j in eachindex(elite_cchr)
        @test elite_cchr[j] == e.controller_sympop[j].chromosome
    end
    # Test: indexes (elites / rows / cols)
    test_indexes(e, elite_echr, elite_cchr, elite_pairs)
    # Test: fitness matrix is reset to -Inf
    matsize = size(e.fitness_matrix)
    @test e.fitness_matrix == -Inf * ones(matsize...)
end

function test_newpop_custom(
    e::DualCGPGAEvo,
    elite_echr::Vector{Vector{Float64}},
    elite_cchr::Vector{Vector{Float64}}
)
    @test sum(e.elites_matrix) == length(custom_elites_indexes)
    expected_elites_matrix = falses(5, 5)
    for index in [(1,1), (1,2), (2,3)]
        expected_elites_matrix[index...] = true
    end
    @test e.elites_matrix == expected_elites_matrix
    @test e.encoder_sympop[1].chromosome == elite_echr[1]
    @test e.encoder_sympop[2].chromosome == elite_echr[3]
    @test e.controller_sympop[1].chromosome == elite_cchr[1]
    @test e.controller_sympop[2].chromosome == elite_cchr[2]
    @test e.controller_sympop[3].chromosome == elite_cchr[3]
end

function test_tournaments_ind(
    e::DualCGPGAEvo,
    elite_echr::Vector{Vector{Float64}},
    elite_cchr::Vector{Vector{Float64}},
    prev_echr::Vector{Vector{Float64}},
    prev_cchr::Vector{Vector{Float64}},
    comptype::String
)
    for esym in e.encoder_sympop
        if comptype == "in"
            # current chromosome was in previous pop
            @test esym.chromosome in prev_echr
        elseif comptype == "notin"
            if esym.chromosome in elite_echr
                # current chromosome was in previous pop
                @test esym.chromosome in prev_echr
            else
                # current chromosome was NOT in previous pop
                @test esym.chromosome ∉ prev_echr
            end
        end
    end
    for csym in e.controller_sympop
        if comptype == "in"
            @test csym.chromosome in prev_cchr
        elseif comptype == "notin"
            if csym.chromosome in elite_cchr
                @test csym.chromosome in prev_cchr
            else
                @test csym.chromosome ∉ prev_cchr
            end
        end
    end
end

cfg_filename = string(@__DIR__, "/dualcgpga_test.yaml")
game = "gravitar"
mcfg, ecfg, ccfg, reducer, bootstrap = IICGP.dualcgp_config(cfg_filename, game)
logid = "2021-test-logid"
resdir = dirname(@__DIR__)
function mutate(ind::CGPInd, ind_type::String)
    if ind_type == "encoder"
        return goldman_mutate(ecfg, ind, init_function=IPCGPInd)
    elseif ind_type == "controller"
        return goldman_mutate(ccfg, ind)
    end
end
fit(e::CGPInd, c::CGPInd, seed::Int64) = [sum(e.chromosome) * sum(c.chromosome)]
n_iter = 1000
n_seed = 5

@testset "CGP GA Populate" begin
    for t in 1:n_seed
        Random.seed!(t)
        evo = IICGP.DualCGPGAEvo(mcfg, ecfg, ccfg, fit, logid, resdir)
        prevfitnesses = [-Inf]
        for i in 1:n_iter
            evaluate(evo)
            maxfit = maximum(evo.fitness_matrix)
            @test maxfit > -Inf
            @test maxfit >= maximum(prevfitnesses)
            push!(prevfitnesses, maxfit)
            elite_echr, elite_cchr, elite_pairs = get_elites(evo)
            prev_echr = [esym.chromosome for esym in evo.encoder_sympop]
            prev_cchr = [csym.chromosome for csym in evo.controller_sympop]
            populate(evo)
            test_newpop(evo, elite_echr, elite_cchr, elite_pairs)
            test_tournaments_ind(evo, elite_echr, elite_cchr, prev_echr, prev_cchr, "notin")
        end
    end
end

custom_elites_indexes = [(2,2), (2,3), (4,4)]

@testset "CGP GA Populate with custom evaluation" begin
    evo = IICGP.DualCGPGAEvo(mcfg, ecfg, ccfg, fit, logid, resdir)
    for i in 1:n_iter
        evaluate(evo)
        # Custom evaluation
        ne, nc = ecfg.n_population, ccfg.n_population
        nmax = max(ne, nc)
        matsize = (ne, nc)
        evo.elites_matrix = falses(matsize...)
        evo.fitness_matrix = -Inf * ones(matsize...)
    	for i in 1:nmax
            i_modrow = i-ne*divrem(i-1,ne)[1]
            i_modcol = i-nc*divrem(i-1,nc)[1]
    		evo.fitness_matrix[i_modrow, i_modcol] = 1.0
    	end
        elite_echr = Vector{Vector{Float64}}()
        elite_cchr = Vector{Vector{Float64}}()
        if ne == 5 && nc == 5
            for index in custom_elites_indexes
                evo.elites_matrix[index...] = true
                evo.fitness_matrix[index...] = 3.0
                row, col = index
                push!(elite_echr, evo.encoder_sympop[row].chromosome)
                push!(elite_cchr, evo.controller_sympop[col].chromosome)
            end
        end
        populate(evo)
        test_newpop_custom(evo, elite_echr, elite_cchr)
    end
end

#=
function mutate(ind::CGPInd, ind_type::String)
    ind
end

@testset "CGP GA Populate without mutation" begin
    evo = IICGP.DualCGPGAEvo(mcfg, ecfg, ccfg, fit, logid, resdir)
    for i in 1:3
        evaluate(evo)
        elite_echr, elite_cchr, elite_pairs = get_elites(evo)
        prev_echr = [esym.chromosome for esym in evo.encoder_sympop]
        prev_cchr = [csym.chromosome for csym in evo.controller_sympop]
        populate(evo)
        test_newpop(evo, elite_echr, elite_cchr, elite_pairs)
        test_tournaments_ind(evo, elite_echr, elite_cchr, prev_echr, prev_cchr, "in")
    end
end
=#


text:
include("../src/conical_radon.jl")
C = conical_radon(reverse(rand(100,100)^0, dims=1), π/4, 1)

ext: jl


text:
include("mis1.jl")
include("mis2.jl")


text:
path_valid = "C:\\Users\\wilhe\\Dropbox\\My PC (DESKTOP-P6322LG)\\Desktop\\Final Training Work\\Good ANNs\\Subsea Classifier\\SILU_2l_2n_1sig.csv"
nn_data_raw_valid = CSV.read(path_valid, DataFrame)
size_df_valid = size(nn_data_raw_valid)

input_len = 4
hidden_len = 2
output_len = 1
depth = 2
g = swish1

xmin = [8.23365; 0.35171; 0.35; 0.10012]
xmax = [19.5046; 0.8; 0.499985; 0.25]

ymin1 = 0.0
ymax1 = 1.0

# extracts vector of W and b that
W = Matrix{Float64}[]
b = Vector{Float64}[]

# extract first hidden layer info
H1frame = filter(row -> row.Layer === "H1", nn_data_raw_valid)
H1frame_W = H1frame[:, [Symbol("W"*string(i)) for i=1:input_len]]
H1frame_b = H1frame[:, [:b]]
Welement = Matrix{Float64}(H1frame_W)
belement = Matrix{Float64}(H1frame_b)
push!(W, Welement)
push!(b, Float64[belement[i] for i = 1:size(belement,1)])

# extract info for 2 to last hidden layer
for i = 2:depth
    Hiframe = filter(row -> row.Layer === "H$i", nn_data_raw_valid)
    Hiframe_W = Hiframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
    Hiframe_b = Hiframe[:, [:b]]
    Welement2 = Matrix{Float64}(Hiframe_W)
    belement2 = Matrix{Float64}(Hiframe_b)
    push!(W, Welement2)
    push!(b, Float64[belement2[i] for i = 1:size(belement2,1)])
end

# extract info for output layer
Hoframe = filter(row -> row.Layer === "Output", nn_data_raw_valid)
Hoframe_W = Hoframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
Hoframe_b = Hoframe[:, [:b]]
Welement = Matrix{Float64}(Hoframe_W)
belement = Matrix{Float64}(Hoframe_b)
push!(W, Welement)
push!(b, Float64[belement[i] for i = 1:size(belement,1)])

function ANN_valid(x)
    xs = (x .- xmin)./(xmax .- xmin)

    yi1 = W[1]*xs .+ b[1]
    y1 = swish1.(yi1)

    yi2 = W[2]*y1 .+ b[2]
    y2 = swish1.(yi2)

    yi3 = W[3]*y2 .+ b[3]
    y3 = sigmoid.(yi3)

    return y3[1]*(ymax1 - ymin1) + ymin1
    #return y4[1]*(ymax1 - ymin1) + ymin1
end

path_perform = "C:\\Users\\wilhe\\Dropbox\\My PC (DESKTOP-P6322LG)\\Desktop\\Final Training Work\\Good ANNs\\Subsea Surrogate\\SILU_2l_12n_4sig.csv"
nn_data_raw_perform = CSV.read(path_perform, DataFrame)
size_df_perform = size(nn_data_raw_perform)

input_len = 4
hidden_len = 12
output_len = 4
depth = 2
g = swish1

# extracts vector of W and b that
W2 = Matrix{Float64}[]
b2 = Vector{Float64}[]

# extract first hidden layer info
H1frame = filter(row -> row.Layer === "H1", nn_data_raw_perform)
H1frame_W = H1frame[:, [Symbol("W"*string(i)) for i=1:input_len]]
H1frame_b = H1frame[:, [:b]]
Welement = Matrix{Float64}(H1frame_W)
belement = Matrix{Float64}(H1frame_b)
push!(W2, Welement)
push!(b2, Float64[belement[i] for i = 1:size(belement,1)])

# extract info for 2 to last hidden layer
for i = 2:depth
    Hiframe = filter(row -> row.Layer === "H$i", nn_data_raw_perform)
    Hiframe_W = Hiframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
    Hiframe_b = Hiframe[:, [:b]]
    Welement2 = Matrix{Float64}(Hiframe_W)
    belement2 = Matrix{Float64}(Hiframe_b)
    push!(W2, Welement2)
    push!(b2, Float64[belement2[i] for i = 1:size(belement2,1)])
end

# extract info for output layer
Hoframe = filter(row -> row.Layer === "Output", nn_data_raw_perform)
Hoframe_W = Hoframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
Hoframe_b = Hoframe[:, [:b]]
Welement = Matrix{Float64}(Hoframe_W)
belement = Matrix{Float64}(Hoframe_b)
push!(W2, Welement)
push!(b2, Float64[belement[i] for i = 1:size(belement,1)])

ymin2 = [0.462387; 0.00946356; 4e+06; 825.87]
ymax2 = [0.799078; 0.138537; 4.01e+06; 911.39]

function ANN_perform(x)
    xs = (x .- xmin)./(xmax .- xmin)

    yi1 = W2[1]*xs .+ b2[1]
    y1 = swish1.(yi1)

    yi2 = W2[2]*y1 .+ b2[2]
    y2 = swish1.(yi2)

    yi3 = W2[3]*y2 .+ b2[3]
    y3 = sigmoid.(yi3)

    ys5 = (ymax2 .- ymin2).*y3 .+ ymin2

    #return y5[1], y5[2], y5[3], y5[4]
    return ys5[1], ys5[2], ys5[3], ys5[4]
end

max_stars_repo_name: PSORLab/RobustHybridModels

max_issues_repo_name: PSORLab/RobustHybridModels

max_forks_repo_name: PSORLab/RobustHybridModels


text:
import Test

import CloudMicrophysics
import CLIMAParameters
import Thermodynamics

const TT = Test
const AM = CloudMicrophysics.AerosolModel
const AA = CloudMicrophysics.AerosolActivation
const CP = CLIMAParameters
const TD = Thermodynamics

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# Atmospheric conditions
T = 294.0       # air temperature K
p = 100000.0    # air pressure Pa
w = 0.5         # vertical velocity m/s

# We need the phase partition here only so that we can compute the
# moist air R_m and cp_m in aerosol activation module.
# We are assuming here saturated conditions and no liquid water or ice.
# This is consistent with the assumptions of the aerosol activation scheme.
p_vs = TD.saturation_vapor_pressure(param_set, T, TD.Liquid())
q_vs = 1 / (1 - CP.Planet.molmass_ratio(param_set) * (p_vs - p) / p_vs)
q = TD.PhasePartition(q_vs, 0.0, 0.0)

# Accumulation mode
r_dry_accum = 0.243 * 1e-6 # m
stdev_accum = 1.4          # -
N_accum = 100.0 * 1e6      # 1/m3

# Coarse Mode
r_dry_coarse = 1.5 * 1e-6  # m
stdev_coarse = 2.1         # -
N_coarse = 1.0 * 1e6       # 1/m3

# Abdul-Razzak and Ghan 2000 mode
r_dry_paper = 0.05 * 1e-6  # m
stdev_paper = 2.0          # -
N_1_paper = 100.0 * 1e6    # 1/m3

# TODO - move areosol properties to CLIMAParameters
# Sea Salt - universal parameters
M_seasalt = 0.058443
ρ_seasalt = 2170.0
ϕ_seasalt = 0.9
ν_seasalt = 2.0
ϵ_seasalt = 1.0
κ_seasalt = 1.12   # TODO - which value to take?
# Sulfate - universal parameters
M_sulfate = 0.132
ρ_sulfate = 1770.0
ϕ_sulfate = 1.0
ν_sulfate = 3.0
ϵ_sulfate = 1.0
κ_sulfate = 0.53   # TODO - which value to take?

# Aerosol size distribution modes
accum_seasalt_B = AM.Mode_B(
    r_dry_accum,
    stdev_accum,
    N_accum,
    (1.0,),
    (ϵ_seasalt,),
    (ϕ_seasalt,),
    (M_seasalt,),
    (ν_seasalt,),
    (ρ_seasalt,),
    1,
)
accum_seasalt_κ = AM.Mode_κ(
    r_dry_accum,
    stdev_accum,
    N_accum,
    (1.0,),
    (1.0,),
    (M_seasalt,),
    (κ_seasalt,),
    1,
)

coarse_seasalt_B = AM.Mode_B(
    r_dry_coarse,
    stdev_coarse,
    N_coarse,
    (1.0,),
    (ϵ_seasalt,),
    (ϕ_seasalt,),
    (M_seasalt,),
    (ν_seasalt,),
    (ρ_seasalt,),
    1,
)
coarse_seasalt_κ = AM.Mode_κ(
    r_dry_coarse,
    stdev_coarse,
    N_coarse,
    (1.0,),
    (1.0,),
    (M_seasalt,),
    (κ_seasalt,),
    1,
)

paper_mode_1_B = AM.Mode_B(
    r_dry_paper,
    stdev_paper,
    N_1_paper,
    (1.0,),
    (ϵ_sulfate,),
    (ϕ_sulfate,),
    (M_sulfate,),
    (ν_sulfate,),
    (ρ_sulfate,),
    1,
)
paper_mode_1_κ = AM.Mode_κ(
    r_dry_paper,
    stdev_paper,
    N_1_paper,
    (1.0,),
    (1.0,),
    (M_sulfate,),
    (κ_sulfate,),
    1,
)

# Aerosol size distributions
AM_1_B = AM.AerosolDistribution((accum_seasalt_B,))
AM_2_B = AM.AerosolDistribution((coarse_seasalt_B, accum_seasalt_B))
AM_3_B = AM.AerosolDistribution((accum_seasalt_B, coarse_seasalt_B))

AM_1_κ = AM.AerosolDistribution((accum_seasalt_κ,))
AM_2_κ = AM.AerosolDistribution((coarse_seasalt_κ, accum_seasalt_κ))
AM_3_κ = AM.AerosolDistribution((accum_seasalt_κ, coarse_seasalt_κ))

TT.@testset "callable and positive" begin

    for AM_t in (AM_3_B, AM_3_κ)
        TT.@test all(AA.mean_hygroscopicity_parameter(param_set, AM_t) .> 0.0)
        TT.@test AA.max_supersaturation(param_set, AM_t, T, p, w, q) > 0.0
        TT.@test all(
            AA.N_activated_per_mode(param_set, AM_t, T, p, w, q) .> 0.0,
        )
        TT.@test all(
            AA.M_activated_per_mode(param_set, AM_t, T, p, w, q) .> 0.0,
        )
        TT.@test AA.total_N_activated(param_set, AM_t, T, p, w, q) > 0.0
        TT.@test AA.total_M_activated(param_set, AM_t, T, p, w, q) > 0.0
    end
end

TT.@testset "same mean hygroscopicity for the same aerosol" begin

    TT.@test AA.mean_hygroscopicity_parameter(param_set, AM_3_B)[1] ==
             AA.mean_hygroscopicity_parameter(param_set, AM_1_B)[1]

    TT.@test AA.mean_hygroscopicity_parameter(param_set, AM_3_B)[2] ==
             AA.mean_hygroscopicity_parameter(param_set, AM_1_B)[1]

    TT.@test AA.mean_hygroscopicity_parameter(param_set, AM_3_κ)[1] ==
             AA.mean_hygroscopicity_parameter(param_set, AM_1_κ)[1]

    TT.@test AA.mean_hygroscopicity_parameter(param_set, AM_3_κ)[2] ==
             AA.mean_hygroscopicity_parameter(param_set, AM_1_κ)[1]

end

TT.@testset "B and kappa hygroscopicities are equivalent" begin

    TT.@test all(isapprox(
        AA.mean_hygroscopicity_parameter(param_set, AM_3_κ)[2],
        AA.mean_hygroscopicity_parameter(param_set, AM_3_B)[2],
        rtol = 0.1,
    ))
end

TT.@testset "order of modes does not matter" begin

    TT.@test AA.total_N_activated(param_set, AM_3_B, T, p, w, q) ==
             AA.total_N_activated(param_set, AM_2_B, T, p, w, q)
    TT.@test AA.total_M_activated(param_set, AM_3_B, T, p, w, q) ==
             AA.total_M_activated(param_set, AM_2_B, T, p, w, q)

    TT.@test AA.total_N_activated(param_set, AM_3_κ, T, p, w, q) ==
             AA.total_N_activated(param_set, AM_2_κ, T, p, w, q)
    TT.@test AA.total_M_activated(param_set, AM_3_κ, T, p, w, q) ==
             AA.total_M_activated(param_set, AM_2_κ, T, p, w, q)
end

TT.@testset "Abdul-Razzak and Ghan 2000 Fig 1" begin

    # data read from Fig 1 in Abdul-Razzak and Ghan 2000
    # using https://automeris.io/WebPlotDigitizer/
    N_2_obs = [
        18.74716810149539,
        110.41572270049846,
        416.00589034889026,
        918.1014952424102,
        1914.816492976891,
        4919.913910285455,
    ]
    N_act_obs = [
        0.7926937018577255,
        0.7161078386950611,
        0.5953670140462167,
        0.4850589034888989,
        0.34446080652469424,
        0.162630267331219,
    ]

    N_act_frac_B = Vector{Float64}(undef, 6)
    N_act_frac_κ = Vector{Float64}(undef, 6)

    it = 1
    for N_2_paper in N_2_obs
        paper_mode_2_B = AM.Mode_B(
            r_dry_paper,
            stdev_paper,
            N_2_paper * 1e6,
            (1.0,),
            (ϵ_sulfate,),
            (ϕ_sulfate,),
            (M_sulfate,),
            (ν_sulfate,),
            (ρ_sulfate,),
            1,
        )
        paper_mode_2_κ = AM.Mode_κ(
            r_dry_paper,
            stdev_paper,
            N_2_paper * 1e6,
            (1.0,),
            (1.0,),
            (M_sulfate,),
            (κ_sulfate,),
            1,
        )

        AD_B = AM.AerosolDistribution((paper_mode_1_B, paper_mode_2_B))
        AD_κ = AM.AerosolDistribution((paper_mode_1_κ, paper_mode_2_κ))

        N_act_frac_B[it] =
            AA.N_activated_per_mode(param_set, AD_B, T, p, w, q)[1] / N_1_paper
        N_act_frac_κ[it] =
            AA.N_activated_per_mode(param_set, AD_κ, T, p, w, q)[1] / N_1_paper
        it += 1
    end

    TT.@test all(isapprox(N_act_frac_B, N_act_obs, rtol = 0.05))
    TT.@test all(isapprox(N_act_frac_κ, N_act_obs, rtol = 0.1))

end


text:
import SymbolicUtils: symtype, term, hasmetadata
struct MTKParameterCtx end

isparameter(x::Num) = isparameter(value(x))
isparameter(x::Symbolic) = getmetadata(x, MTKParameterCtx, false)
isparameter(x) = false

"""
    toparam(s::Sym)

Maps the variable to a paramter.
"""
function toparam(s)
    if s isa Symbolics.Arr
        Symbolics.wrap(toparam(Symbolics.unwrap(s)))
    elseif s isa AbstractArray
        map(toparam, s)
    elseif symtype(s) <: AbstractArray
        Symbolics.recurse_and_apply(toparam, s)
    else
        setmetadata(s, MTKParameterCtx, true)
    end
end
toparam(s::Num) = Num(toparam(value(s)))

"""
    tovar(s::Sym)

Maps the variable to a state.
"""
tovar(s::Symbolic) = setmetadata(s, MTKParameterCtx, false)
tovar(s::Num) = Num(tovar(value(s)))

"""
$(SIGNATURES)

Define one or more known variables.
"""
macro parameters(xs...)
    Symbolics._parse_vars(:parameters,
                          Real,
                          xs,
                          toparam) |> esc
end


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

export is_unit, is_torsion_unit, is_independent, unit_group

add_verbose_scope(:UnitGroup)

################################################################################
#
#  Initialization
#
################################################################################

function _unit_group_init(O::NfMaxOrd)
  u = UnitGrpCtx{FacElem{nf_elem}}(O)
  return u
end

################################################################################
#
#  Field access
#
################################################################################

order(u::UnitGrpCtx) = u.order

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, U::UnitGrpCtx)
  print(io, "Unit group context of\n$(order(U))\n")
end

################################################################################
#
#  Unit rank
#
################################################################################

doc"""
***
    unit_rank(O::NfOrd) -> Int

> Returns the unit rank of $\mathcal O$, that is, the rank of the unit group
> $\mathcal O^\times$.
"""
function unit_rank(O::NfOrd)
  r1, r2 = signature(nf(O))
  return r1 + r2 - 1
end

################################################################################
#
#  Testing for invertibility
#
################################################################################

doc"""
***
    is_unit(x::NfOrdElem) -> Bool

> Returns whether $x$ is invertible or not.
"""
function is_unit(x::NfOrdElem)
  return abs(norm(x)) == 1
end

_is_unit(x::NfOrdElem) = is_unit(x)

function _is_unit{T <: Union{nf_elem, FacElem{nf_elem}}}(x::T)
  return abs(norm(x)) == 1
end

################################################################################
#
#  Test if units are independent
#
################################################################################

doc"""
***
    is_independent{T}(x::Array{T, 1})

> Given an array of non-zero elements in a number field, returns whether they
> are multiplicatively independent.
"""
function is_independent{T}(x::Array{T, 1}, p::Int = 32)
  return _is_independent(x, p)
end

function _is_independent{T}(x::Array{T, 1}, p::Int = 32)
  K = _base_ring(x[1])

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  # This can be made more memory friendly
  while true
    @assert p != 0

    conlog = conjugates_arb_log(x[1], p)

    A = MatrixSpace(parent(conlog[1]), length(x), rr)()

    for i in 1:rr
      A[1, i] = conlog[i]
    end

    Ar = base_ring(A)

    for k in 2:length(x)
      conlog = conjugates_arb_log(x[k], p)
      for i in 1:rr
        A[k, i] = conlog[i]
      end
    end

    B = A*transpose(A)
    @vprint :UnitGroup 1 "Computing det of $(rows(B))x$(cols(B)) matrix with precision $(p) ... \n"
    d = det(B)

    y = (Ar(1)//Ar(r))^r * (Ar(21)//Ar(128) * log(Ar(deg))//(Ar(deg)^2))^(2*r)
    if isfinite(d) && ispositive(y - d)
      return false, p
    elseif isfinite(d) && ispositive(d)
      return true, p
    end
    p = 2*p
  end
end

function _is_independent{T}(u::UnitGrpCtx{T}, y::FacElem{T})
  K = _base_ring(x[1])
  p = u.indep_prec

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  # This can be made more memory friendly
  while true
    @assert p != 0

#    if length(u.units) == 0
#      A = _conj_log_mat([y], p)
#      B = MatrixSpace(base_ring(A), rows(A), rows(A))
#      z.conj_log_mat_tranpose = (transpose(A), p)
#      mul!(B, A, z.conj_log_mat_transpose)
#      z.conj_log_mat = (A, p)
#      z.conj_log_mat_times_transpose = (B, p)
#    end
#
#    if z.conj_log_mat[2] == p # the old matrix has the same precision as our working precision
#      A = z.conj_log_mat[1]
#      conjlog = conjugates_arb_log(y, p)
#      newrow = MatrixSpace(base_ring(A), 1, cols(A))()
#      newcol = MatrixSpace(base_ring(A), cols(A), 1)()
#      z.conj_log_mat = (vcat(A, newrow), p)
#      z.conj_log_mat_transpose = (hcat(z.conj_log_mat_transpose[1], newcol), p)
#
#      for i in 1:cols(A)
#        z.conj_log_mat[1][rows(A) + 1, i] = conjlog[i]
#        z.conj_log_mat_transpose[1][i, rows(A) + 1]
#      end
#      A = z.conj_log_mat[1]
#    end
    A = _conj_log_mat(u.units, p)

    B = A*transpose(A)
    @vprint :UnitGroup 1 "Computing det of $(rows(B))x$(cols(B)) matrix with precision $(p) ... \n"
    d = det(B)

    y = (Ar(1)//Ar(r))^r * (Ar(21)//Ar(128) * log(Ar(deg))//(Ar(deg)^2))^(2*r)
    if isfinite(d) && ispositive(y - d)
      return false, p
    elseif isfinite(d) && ispositive(d)
      return true, p
    end
    p = 2*p
  end
end

# Checks whether x[1]^z[1] * ... x[n]^z[n]*y^[n+1] is a torsion unit
# This can be improved
function _check_relation_mod_torsion{T}(x::Array{T, 1}, y::T, z::Array{fmpz, 1})
  (length(x) + 1 != length(z)) && error("Lengths of arrays does not fit")
  r = x[1]^z[1]

  for i in 2:length(x)
    r = r*x[i]^z[i]
  end

  w = r*y^z[length(z)]

  b, _ = is_torsion_unit(w)
  return b
end

function _find_rational_relation!(rel::Array{fmpz, 1}, v::arb_mat, bound::fmpz)
  @vprint :UnitGroup 2 "Finding rational approximation in $v\n"
  r = length(rel) - 1

  z = Array(fmpq, r)

  # Compute an upper bound in the denominator of an entry in the relation
  # using Cramer's rule and lower regulator bounds

  # No comes the rational approximation phase

  # First a trivial check:
  # If the relation contains only integers, it does not yield any information

  i = 0

  is_integer = true

  while is_integer && i < r
    i = i + 1
    b, o = unique_integer(v[1, i])
    if b
      rel[i] = o
    end
    is_integer = is_integer && b
  end

  if is_integer
    rel[r + 1] = -1
    @vprint :UnitGroup 2 "Found rational relation.\n"
    return true
  end

  for i in 1:r
    if radius(v[1, i]) > 1
      # This is a strange case I cannot handle at the moment
      return false
      #no_change_matrix = MatrixSpace(ZZ, r, r)(1)
      #no_change_matrix = vcat(no_change_matrix, MatrixSpace(ZZ, 1, r)(0))
      #return x, no_change_matrix
    end

    app =  _frac_bounded_2(v[1, i], bound)
    if app[1]
      z[i] = app[2]
    else
      @vprint :UnitGroup 2 "Something went wrong with the approximation.\n"
      return false
    end
  end

  dlcm = den(z[1])

  for i in 2:length(z)
    dlcm = lcm(dlcm, den(z[i]))
  end

  for i in 1:r
    rel[i] = num(z[i]*dlcm)
  end

  rel[r + 1] = -dlcm

  # Check that relation is primitive
  g = rel[1]

  for i in 1:length(rel)
    g = gcd(g, rel[i])
    if g == 1
      break
    end
  end

  @assert g == 1

  @vprint :UnitGroup 2 "Found rational relation.\n"
  return true
end

# Given r elements x_1,...,x_r, where r is the unit rank, and y an additional unit,
# compute a basis z_1,...z_r such that <x_1,...x_r,y,T> = <z_1,...,z_r,T>,
# where T are the torsion units
function _find_relation{S, T}(x::Array{S, 1}, y::T, p::Int = 64)

  K = _base_ring(x[1])

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  R = ArbField(p)

  zz = Array(fmpz, r + 1)

  @vprint :UnitGroup 1 "Computing conjugates log matrix ... \n"
  A = _conj_log_mat_cutoff(x, p)

  Ar = base_ring(A)

  b = ArbMatSpace(Ar, 1, r)()

  conlog = conjugates_arb_log(y, p)

  for i in 1:r
    b[1, i] = conlog[i]
  end

  B = parent(A)()

  # I should do this using lu decomposition and caching
  # The inversion could go wrong,
  # Then we again have to increase the precision

  inv_succesful = false

  try
    @vprint :UnitGroup 1 "Inverting matrix ... \n"
    B = inv(A)
    inv_succesful = true
  catch e
    @vprint :UnitGroup 1 "Cannot invert matrix ... \n"
    rethrow(e)
  end

  v = b*B

  z = Array(fmpq, r)

  rreg = det(A)

  bound = _denominator_bound_in_relation(rreg, K)

  # Compute an upper bound in the denominator of an entry in the relation
  # using Cramer's rule and lower regulator bounds

  rel = Array(fmpz, r + 1)
  for i in 1:r+1
    rel[i] = zero(FlintZZ)
  end

  while !inv_succesful || !_find_rational_relation!(rel, v, bound)
    p =  2*p

    inv_succesful = false

    A = _conj_log_mat_cutoff(x, p)

    Ar = base_ring(A)

    b = ArbMatSpace(Ar, 1, r)()

    conlog = conjugates_arb_log(y, p)

    for i in 1:r
      b[1, i] = conlog[i]
    end

    if !inv_succesful
      try
        @vprint :UnitGroup 1 "Inverting matrix ... \n"
        B = inv(A)
        inv_succesful = true
      catch
        @vprint :UnitGroup 1 "Cannot invert matrix. Increasing precision to $(2*p)\n"
      end
    end

    v = b*B
  end

  # Check if it is a relation modulo torsion units!
  @vprint :UnitGroup 1 "Checking relation $rel \n"

  if !_check_relation_mod_torsion(x, y, rel)
    #error("Dirty approximation did not work")
    return _find_relation(x, y, 2*p)
    #rel[r + 1 ] = 0
    #return rel
  end

  @vprint :UnitGroup 1 "Found a valid relation!\n"
  return rel
end

function _denominator_bound_in_relation(rreg::arb, K::AnticNumberField)
  # Compute an upper bound in the denominator of an entry in the relation
  # using Cramer's rule and lower regulator bounds

  arb_bound = rreg * inv(lower_regulator_bound(K))

  # I want to get an upper bound as an fmpz
  tm = arf_struct(0, 0, 0, 0)
  ccall((:arf_init, :libarb), Void, (Ptr{arf_struct}, ), &tm)

  ccall((:arb_get_abs_ubound_arf, :libarb), Void, (Ptr{arf_struct}, Ptr{arb}, Int), &tm, &arb_bound, 64)

  bound = fmpz()

  # round towards +oo
  ccall((:arf_get_fmpz, :libarb), Void, (Ptr{fmpz}, Ptr{arf_struct}, Cint), &bound, &tm, 3)

  ccall((:arf_clear, :libarb), Void, (Ptr{arf_struct}, ), &tm)

  return bound
end

function _transform(x::Array{nf_elem}, y::fmpz_mat)
  n = length(x)
  @assert n == rows(y)
  m = cols(y)
  z = Array(nf_elem, m)
  for i in 1:m
    z[i] = x[1]^y[1, i]
    for j in 2:n
      z[i] = z[i]*x[j]^y[j, i]
    end
  end
  return z
end

function _frac_bounded_2(y::arb, bound::fmpz)
  p = prec(parent(y))
  x = _arb_get_fmpq(y)
  n = 1
  c = cfrac(x, n)[1]
  q = fmpq(c)

  new_q = q

  while nbits(num(new_q)) < div(p, 2) && nbits(den(new_q)) < div(p, 2) && den(new_q) < bound

    if contains(y, new_q)
      return true, new_q
    end

    n = n + 1
    c = cfrac(x, n)[1]
    new_q = fmpq(c)

  end
  return false, zero(FlintQQ)
end

function _max_frac_bounded(x::fmpq, b::fmpz)
  n = 2
  c = cfrac(x, n)[1]
  q = fmpq(c)

  while abs(den(q)) < b && q != x
    n = 2*n
    c = cfrac(x, n)[1]
    q = fmpq(c)
  end

  while abs(den(q)) > b
    n = n - 1
    c = cfrac(x, n)[1]
    q = fmpq(c)
  end

  return n
end

function _rel_add_prec(U)
  return U.rel_add_prec
end

function _add_dependent_unit{S, T}(U::UnitGrpCtx{S}, y::T)
  K = nf(order(U))
  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  p = _rel_add_prec(U)

  #p = 64

  zz = Array(fmpz, r + 1)

  @v_do :UnitGroup 1 pushindent()
  p, B = _conj_log_mat_cutoff_inv(U, p)
  @v_do :UnitGroup 1 popindent()
  @vprint :UnitGroup 2 "Precision is now $p\n"

  Ar = base_ring(B)

  b = ArbMatSpace(Ar, 1, r)()

  conlog = conjugates_arb_log(y, p)

  for i in 1:r
    b[1, i] = conlog[i]
  end

  # I should do this using lu decomposition and caching
  # The inversion could go wrong,
  # Then we again have to increase the precision

  inv_succesful = true

  @vprint :UnitGroup 3 "For $p element b: $b\n"
  v = b*B
  @vprint :UnitGroup 3 "For $p the vector v: $v\n"

  z = Array(fmpq, r)

  rreg = arb()

  if isdefined(U, :tentative_regulator)
    rreg = U.tentative_regulator
  else
    rreg = regulator(U.units, 64)
  end

  bound = _denominator_bound_in_relation(rreg, K)

    rel = Array(fmpz, r + 1)
  for i in 1:r+1
    rel[i] = zero(FlintZZ)
  end

  @vprint :UnitGroup 2 "First iteration to find a rational relation ... \n"
  while !_find_rational_relation!(rel, v, bound)
    @vprint :UnitGroup 2 "Precision not high enough, increasing from $p to $(2*p)\n"
    p =  2*p

    p, B = _conj_log_mat_cutoff_inv(U, p)

    conlog = conjugates_arb_log(y, p)

    b = ArbMatSpace(parent(conlog[1]), 1, r)()

    for i in 1:r
      b[1, i] = conlog[i]
    end

    @vprint :UnitGroup 3 "For $p element b: $b\n"

    v = b*B
    @vprint :UnitGroup 3 "For $p the vector v: $v\n"
  end

  @vprint :UnitGroup 3 "For $p rel: $rel\n"

  @vprint :UnitGroup 2 "Second iteration to check relation ... \n"
  while !_check_relation_mod_torsion(U.units, y, rel)
    @vprint :UnitGroup 2 "Precision not high enough, increasing from $p to $(2*p)\n"
    p = 2*p
    p, B = _conj_log_mat_cutoff_inv(U, p)

    conlog = conjugates_arb_log(y, p)

    b = ArbMatSpace(parent(conlog[1]), 1, r)()

    for i in 1:r
      b[1, i] = conlog[i]
    end

    @vprint :UnitGroup 3 "For $p element b: $b\n"
    v = b*B
    @vprint :UnitGroup 3 "For $p the vector v: $v\n"
    _find_rational_relation!(rel, v, bound)
    @vprint :UnitGroup 3 "For $p rel: $rel\n"
  end

  if abs(rel[r + 1]) == 1 || rel[r + 1] == 0
    U.rel_add_prec = p
    return false
  end

  m = MatrixSpace(FlintZZ, r + 1, 1)(reshape(rel, r + 1, 1))

  h, u = hnf_with_transform(m)

  @assert h[1,1] == 1

  u = inv(u)

  m = submat(u, 1:r+1, 2:r+1)

  U.units =  _transform(vcat(U.units, y), m)

  U.conj_log_mat_cutoff = Dict{Int, arb_mat}()
  U.conj_log_mat_cutoff_inv = Dict{Int, arb_mat}()
  U.tentative_regulator = regulator(U.units, 64)
  U.rel_add_prec = p
  return true
end

function  _conj_log_mat{T}(x::Array{T, 1}, p::Int)
  conlog = conjugates_arb_log(x[1], p)

  A = MatrixSpace(parent(conlog[1]), length(x), rr)()

  for i in 1:rr
    A[1, i] = conlog[i]
  end

  Ar = base_ring(A)

  for k in 2:length(x)
    conlog = conjugates_arb_log(x[k], p)
    for i in 1:rr
      A[k, i] = conlog[i]
    end
  end
  return
end

function _conj_log_mat_cutoff{T}(x::Array{T, 1}, p::Int)
  r = length(x)
  conlog = conjugates_arb_log(x[1], p)
  A = ArbMatSpace(parent(conlog[1]), r, r)()

  for i in 1:r
    A[1, i] = conlog[i]
  end

  for k in 2:r
    conlog = conjugates_arb_log(x[k], p)
    for i in 1:r
      A[k, i] = conlog[i]
    end
  end
  return A
end

function _conj_log_mat_cutoff(x::UnitGrpCtx, p::Int)
  if haskey(x.conj_log_mat_cutoff,  p)
    @vprint :UnitGroup 2 "Conj matrix for $p cached\n"
    return x.conj_log_mat_cutoff[p]
  else
    @vprint :UnitGroup 2 "Conj matrix for $p not cached\n"
    x.conj_log_mat_cutoff[p] = _conj_log_mat_cutoff(x.units, p)
    return x.conj_log_mat_cutoff[p]
  end
end

function _conj_log_mat_cutoff_inv(x::UnitGrpCtx, p::Int)
  @vprint :UnitGroup 2 "Computing inverse of conjugates log matrix (starting with prec $p) ... \n"
  if haskey(x.conj_log_mat_cutoff_inv,  p)
    @vprint :UnitGroup 2 "Inverse matrix cached for $p\n"
    return p, x.conj_log_mat_cutoff_inv[p]
  else
    @vprint :UnitGroup 2 "Inverse matrix not cached for $p\n"
    try
      @vprint :UnitGroup 2 "Trying to invert conj matrix with prec $p \n"
      @vprint :UnitGroup 3 "Matrix to invert is $(_conj_log_mat_cutoff(x, p))"
      x.conj_log_mat_cutoff_inv[p] = inv(_conj_log_mat_cutoff(x, p))
      @vprint :UnitGroup 2 "Successful. Returning with prec $p \n"
      @vprint :UnitGroup 3 "$(x.conj_log_mat_cutoff_inv[p])\n"
      return p, x.conj_log_mat_cutoff_inv[p]
    catch e
      # I should check that it really is that error
      @vprint :UnitGroup 2 "Increasing precision .."
      @v_do :UnitGroup 2 pushindent()
      r = _conj_log_mat_cutoff_inv(x, 2*p)
      @v_do :UnitGroup 2 popindent()
      return r
    end
  end
end

# Powering function for fmpz exponents
function _pow{T}(x::Array{T, 1}, y::Array{fmpz, 1})
  K = _base_ring(x[1])

  zz = deepcopy(y)

  z = Array(fmpz, length(x))

  for i in 1:length(x)
    z[i] = mod(zz[i], 2)
    zz[i] = zz[i] - z[i]
  end

  r = K(1)

  return zz
end

################################################################################
#
#  Free part of the unit group
#
################################################################################

doc"""
***
    regulator(x::Array{T, 1}, abs_tol::Int) -> arb

> Compute the regulator $r$ of the elements in $x$, such that the radius of $r$
> is less then `-2^abs_tol`.
"""
function regulator{T}(x::Array{T, 1}, abs_tol::Int)
  K = _base_ring(x[1])
  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  p = 32

  while true
    conlog = conjugates_arb_log(x[1], p)

    A = ArbMatSpace(parent(conlog[1]), r, r)()

    for j in 1:r
      A[1, j] = conlog[j]
    end

    for i in 2:r
      conlog = conjugates_arb_log(x[i], p)
      for j in 1:r
        A[i, j] = conlog[j]
      end
    end

    z = abs(det(A))

    if isfinite(z) && radiuslttwopower(z, -abs_tol)
      return z
    end

    p = 2*p
  end
end

function _make_row_primitive(x::fmpz_mat, j::Int)
  y = x[j, 1]
  for i in 1:cols(x)
    y = gcd(y, x[j, i])
  end
  if y > 1
    for i in 1:cols(x)
      x[j, i] = div(x[j, i], y)
    end
  end
end

function _make_row_primitive!(x::Array{fmpz, 1})
  y = x[1]
  for i in 2:length(x)
    y = gcd(y, x[i])
    if y == 1
      return x
    end
  end
  if y > 1
    for i in 1:cols(x)
      x[i] = div(x[i], y)
    end
  end
end

################################################################################
#
#  Compute unit group from class group context
#
################################################################################

function _unit_group(O::NfMaxOrd, c::ClassGrpCtx)
  u = UnitGrpCtx{FacElem{nf_elem}}(O)
  _unit_group_find_units(u, c)
  return u
end

# TH:
# Current strategy
# ================
# Compute some "random" kernel elements using the transformation matrix to first
# find r independent units, where r is the unit rank.
# In the second round, try to enlarge the unit group with some random kernel
# elements.
function _unit_group_find_units(u::UnitGrpCtx, x::ClassGrpCtx)
  @vprint :UnitGroup 1 "Processing ClassGrpCtx to find units ... \n"

  @vprint :UnitGroup 1 "Relation matrix has size $(rows(x.M)) x $(cols(x.M))\n"

  O = order(u)

  #ker = transpose(ker)

  K = nf(order(x.FB.ideals[1]))
  r = unit_rank(O)
  r1, r2 = signature(O)

  A = u.units

  j = 0

  dim_ker = rows(x.M) - rows(x.H)

  kelem = fmpz[ 0 for i in 1:rows(x.M) ]

  MAX_RND_RD_1 = 2*r

  time_indep = 0.0
  time_add_dep_unit = 0.0
  time_kernel = 0.0
  time_torsion = 0.0

  while(length(A) < r)

    for i in 1:length(kelem)
      kelem[i] = 0
    end

    if j > MAX_RND_RD_1
      return 0
    end
    @vprint :UnitGroup 1 "Found $(length(A)) independent units so far ($(r - length(A)) left to find)\n"
    j = j + 1

    for i in 1:dim_ker
      kelem[end - i + 1] = rand(0:1)
    end

    #time_kernel += @elapsed
    for i in length(x.H_trafo):-1:1
      apply_right!(kelem, x.H_trafo[i])
    end

    _make_row_primitive!(kelem)

    y = FacElem(x, kelem)

    time_torsion += @elapsed is_tors, p = is_torsion_unit(y, false, u.tors_prec)
    u.tors_prec = max(p, u.tors_prec)
    if is_tors
      continue
    end

    @vprint :UnitGroup 1 "Exponents are of bit size $(maximum([ nbits(o) for o in kelem]))\n"

    time_indep += @elapsed _add_unit(u, y)

  end

  @vprint :UnitGroup 1 "Found $r linear independent units \n"

  @vprint :UnitGroup 1 "Independent unit time: $time_indep\n"
  @vprint :UnitGroup 1 "Adding dependent unit time: $time_add_dep_unit\n"
  @vprint :UnitGroup 1 "Torsion test time: $time_torsion\n"
  @vprint :UnitGroup 1 "Kernel time: $time_kernel\n"

  u.full_rank = true

  j = 0

  not_larger = 0

  @vprint :UnitGroup 1 "Enlarging unit group by adding remaining kernel basis elements ...\n"
  while not_larger < 5

    for i in 1:length(kelem)
      kelem[i] = 0
    end

    for i in 1:dim_ker
      kelem[end - i + 1] = rand(-2:2)
    end

    time_kernel += @elapsed for i in length(x.H_trafo):-1:1
      apply_right!(kelem, x.H_trafo[i])
    end

    y = FacElem(x, kelem)

    #!is_unit(y) && throw(BlaError(x, kelem))

    @vprint :UnitGroup 2 "Test if kernel element yields torsion unit ... \n"
    @v_do :UnitGroup 2 pushindent()
    time_torsion += @elapsed is_tors, p = is_torsion_unit(y, false, u.tors_prec)
    u.tors_prec = max(p, u.tors_prec)
    if is_tors
      @v_do :UnitGroup 2 popindent()
      #println("torsion unit: $y")
      @vprint :UnitGroup 2 "Element is torsion unit\n"
      continue
    end
    @v_do :UnitGroup 2 popindent()

    @v_do :UnitGroup 2 pushindent()
    time_add_dep_unit += @elapsed m = _add_dependent_unit(u, y)
    @v_do :UnitGroup 2 popindent()

    if !m
      not_larger = not_larger + 1
    else
      not_larger = 0
    end
  end

  u.tentative_regulator = regulator(u.units, 64)

  @vprint :UnitGroup 1 "Finished processing\n"
  @vprint :UnitGroup 1 "Regulator of current unit group is $(u.tentative_regulator)\n"
  @vprint :UnitGroup 1 "-"^80 * "\n"
  @vprint :UnitGroup 1 "Independent unit time: $time_indep\n"
  @vprint :UnitGroup 1 "Adding dependent unit time: $time_add_dep_unit\n"
  @vprint :UnitGroup 1 "Torsion test time: $time_torsion\n"
  @vprint :UnitGroup 1 "Kernel time: $time_kernel\n"
  return 1
end

function _add_unit(u::UnitGrpCtx, x::FacElem{nf_elem})
  isindep, p = is_independent(vcat(u.units, [x]), u.indep_prec)
  u.indep_prec = max(p, u.indep_prec)
  if isindep
    push!(u.units, x)
  end
  nothing
end

################################################################################
#
#  Size reduction
#
################################################################################

function _reduce_size{T}(x::Array{T, 1}, prec::Int = 64)
  K = _base_ring(x[1])

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  conlog = conjugates_arb_log(x[1], prec)

  A = MatrixSpace(parent(conlog[1]), length(x), rr)()

  B = MatrixSpace(FlintZZ, rows(A), cols(A))()

  for i in 1:rr
    A[1, i] = conlog[i]
  end

  Ar = base_ring(A)

  for i in 1:rows(A)
    for j in 1:cols(A)
      b, y = unique_integer(ceil(ldexp(A[i, j], 64)))
      @assert b
      B[i, j] = y
    end
  end

  L, U = lll_with_transform(B)
end

################################################################################
#
#  Saturation
#
################################################################################

# TH:
# Let U = <u_1,...,u_n,z> with z a generator for Tor(U)
# For a prime p the group U/U^p is F_p-vector space of dimension
# rank(U) or rank(U) + 1 (depending on the order of z).
# if p divides N(P) - 1 = #F_P*, then F_P*/F_P*^p is a one-dimensional
# F_p-vector space. Thus the canonical map F_p-linear
#               U/U^p -> F_P*/F_P*^p
# can be described by a 1 x (rank(U)) or 1 x (rank(U) + 1) matrix over F_p,
# and can be computed by solving discrete logarithms in F_P
#
function _is_saturated(U::UnitGrpCtx, p::Int, B::Int = 2^30 - 1, proof::Bool = false)
  if proof
    error("Not yet implemented")
  end

  N = 3*unit_rank(order(U))

  @vprint :UnitGroup 1 "Computing $N prime ideals for saturation ...\n"

  primes =  _find_primes_for_saturation(order(U), p, N, B)

  m = _matrix_for_saturation(U, primes[1], p)

  for i in 2:N
    m = vcat(m, _matrix_for_saturation(U, primes[i], p))
  end

  @vprint :UnitGroup 1 "Computing kernel of p-th power map ...\n"
  (K, k) = _right_kernel(m)

  K = transpose(K)
  L = lift(K)
  T = typeof(L[1,1])

  nonzerorows = Array{Int, 1}()

  for j in 1:rows(L)
    if !is_zero_row(L, j)
      push!(nonzerorows, j)
    end
  end

  if k == 0
    return (true, zero(nf(order(U))))
  elseif k == 1 && sum(T[ L[nonzerorows[1], i]::T for i in 1:cols(L)-1]) == 0
    # Only one root, which is torsion.
    # We assume that the torsion group is the full torsion group
    return (true, zero(nf(order(U))))
  else
    for j in nonzerorows
      a = U.units[1]^(L[j, 1])
      for i in 2:length(U.units)
        a = a*U.units[i]^L[j, i]
      end

      if gcd(p, U.torsion_units_order) != 1
        a = a*elem_in_nf(U.torsion_units_gen)^L[j, length(U.units) + 1]
      end

      b = evaluate(a)

      @vprint :UnitGroup 1 "Testing/computing root ... \n"

      has_root, roota = root(b, p)

      if !has_root
        continue
      end

      return (false, roota)
    end
  end

  # try some random linear combination of kernel vectors

  MAX = 10

  for i in 1:MAX

    ra = rand(0:p-1, rows(K))
    v = MatrixSpace(base_ring(K), 1, cols(K))(0)
    for j in 1:cols(K)
      for l in 1:rows(K)
        v[1, j] = v[1, j] + ra[l]*K[l,j]
      end
    end

    if v == parent(v)(0)# || sum([v[1, j] for j in 1:rows(K)-1]) == 0
      continue
    end

    v = lift(v)

    a = U.units[1]^(v[1, 1])
    for j in 2:length(U.units)
      a = a*U.units[j]^v[1, j]
    end

    if gcd(p, U.torsion_units_order) != 1
      a = a*elem_in_nf(U.torsion_units_gen)^v[1, length(U.units) + 1]
    end

    b = evaluate(a)

    has_root, roota = root(b, p)

    if has_root
      return (false, roota)
    end
  end

  return (true, zero(nf(order(U))))
end

# The output will be of type
# elem_type(MatrixSpace(ResidueRing(ZZ, p), 1, rank(U) ( + 1))), so
# nmod_mat or fmpz_mod_mat
# THIS FUNCTION IS NOT TYPE STABLE
function _matrix_for_saturation(U::UnitGrpCtx, P::NfMaxOrdIdl, p::Int)
  O = order(U)
  K = nf(O)
  F, mF = ResidueField(O, P)
  mK = extend(mF, K)
  g = _primitive_element(F)

  # We have to add the generator of the torsion group
  if gcd(p, U.torsion_units_order) != 1
    res = MatrixSpace(ResidueRing(FlintZZ, p), 1, unit_rank(O) + 1)()
  else
    res = MatrixSpace(ResidueRing(FlintZZ, p), 1, unit_rank(O))()
  end

  t = K()

  for i in 1:length(U.units)
    u = U.units[i]
    y = one(F)

    # P.gen_two should be P-unformizer
    #println("$(P.gen_one), $b, $(P.gen_two)")

    for b in base(u)
      t = b*K(P.gen_two)^(-valuation(b, P))

      if mod(den(t), minimum(P)) == 0
        l = valuation(den(t), P)
        y = y*(mK(t*elem_in_nf(P.anti_uniformizer)^l)*mF(P.anti_uniformizer)^(-l))^u.fac[b]
      else
        y = y*mK(t)^u.fac[b]
      end
    end

    res[1, i] = disc_log(y, g, p)
  end

  if gcd(p, U.torsion_units_order) != 1
    res[1, unit_rank(O) + 1] = disc_log(mF(U.torsion_units_gen), g, p)
  end

  return res
end

# TH:
# This function finds n prime ideals P of O such that p divides N(P) - 1
# Moreover the prime ideals are unramified and min(P) does not divide
# the index of O in the equation order.
#
# The function loops through all prime ideals ordered by the minimum,
# starting at next_prime(st)
function _find_primes_for_saturation(O::NfMaxOrd, p::Int, n::Int,
                                     st::Int = 0)
  res = Array(NfMaxOrdIdl, n)
  i = 0

  q = st
  while i < n
    q = next_prime(q)

    if mod(index(O), q) == 0 || isramified(O, q)
      continue
    end

    lp = prime_decomposition(O, q)

    j = 1

    while j <= length(lp) && i < n
      Q = lp[j]
      if mod(norm(Q[1]) - 1, p) == 0
        i = i + 1
        res[i] = Q[1]
      end
      j = j + 1
    end
  end

  return res
end

function _primitive_element(F::FqNmodFiniteField)
  #println("Computing primitive element of $F")
  #println("Have to factor $(order(F) - 1)")
  fac = factor(order(F) - 1)
  while true
    a = rand(F)
    if iszero(a)
      continue
    end
    is_primitive = true
    for l in keys(fac)
      if isone(a^(div(order(F) - 1, l)))
        is_primitive = false
        break
      end
    end
    if is_primitive
      return a
    end
  end
end

function _refine_with_saturation(c::ClassGrpCtx, u::UnitGrpCtx)
  @vprint :UnitGroup "Enlarging unit group using saturation ... \n"

  b = _validate_class_unit_group(c, u)

  p = 2

  while b > 1
    @vprint :UnitGroup 1 "Saturating at $p ... \n"

    @v_do :UnitGroup 1 pushindent()
    issat, new_unit = _is_saturated(u, p)
    @v_do :UnitGroup 1 popindent()

    while !issat
      #println("I have found a new unit: $new_unit")
      _add_dependent_unit(u, FacElem(new_unit))
      #println("$(u.tentative_regulator)")

      @v_do :UnitGroup 1 pushindent()
      b = _validate_class_unit_group(c, u)
      @v_do :UnitGroup 1 popindent()

      if b == 1
        break
      end

      @v_do :UnitGroup 1 pushindent()
      issat, new_unit = _is_saturated(u, p)
      @v_do :UnitGroup 1 popindent()
    end

    @v_do :UnitGroup 1 pushindent()
    b = _validate_class_unit_group(c, u)
    @v_do :UnitGroup 1 popindent()

    p = next_prime(p)
    if p > b
      break
    end
  end
  return b
end

################################################################################
#
#  High level interface
#
################################################################################

doc"""
***
    unit_group(O::NfMaxOrd) -> Map

> Returns an isomorphism map $f \colon A \to \mathcal O^\times$. Let
> `A = domain(f)`. Then a set of fundamental units of $\mathcal O$ can be
> obtained via `[ f(A[i]) for i in 1:unit_rank(O) ]`.
"""
function unit_group(O::NfMaxOrd)
  if isdefined(O, :unit_group)
    return O.unit_group::AbToNfOrdUnitGrp{Nemo.nf_elem,NfOrdElem{NfMaxOrd}}
  else
    c, U, b = _class_unit_group(O)
    _refine_with_saturation(c, U)
    f = AbToNfOrdUnitGrp(O, U.units, U.torsion_units_gen, U.torsion_units_order)
    O.unit_group = f
    return f
  end
end

function lower_regulator_bound(K::AnticNumberField)
  return ArbField(64)("0.054")
end


text:
# This is a comment.

function fizzbuzz(i)
    if (i % 3 == 0 && i % 5 == 0)
        return "FizzBuzz"
    elseif (i % 3 == 0)
        return "Fizz"
    elseif (i % 5 == 0)
        return "Buzz"
    else
        return i
    end
end

for i = 1:100
    println(fizzbuzz(i))
end


text:
include("gkMCSetup.jl")

function makeSimpleLTD(ys::Vector{Float64}, N::Int64)
  rpio::RPIO{Float64Particle} = RPIO{Float64Particle}(N, length(ys), 1)
  function ltd(in::SVector{d, Float64})
    θ::GKθ = toGKθ(in)
    lp::Float64 = lglogprior(θ)
    if lp == -Inf return -Inf end
    model::Model = makeGKModel(θ, ys)
    return simpleEstimate(model, rpio)
  end
end

logtargetSimple = makeSimpleLTD(ys, N)
PSimple = makeAMKernel(logtargetSimple, d)

Random.seed!(12345)
setRNGs(54321)

simpleElapsed = @elapsed (chainSimple = simulateChainProgress(PSimple, t0, chainLength))
sar = PSimple(:acceptanceRate)
simpleCovEstimate = PSimple(:covEstimate)

@save "data/gkMCSimple.jld" chainSimple sar simpleCovEstimate simpleElapsed


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use MathOptInterface and
# Knitro to solve the following # simple nonlinear optimization problem.
# This model is test problem # HS15 from the Hock & Schittkowski collection.
#
# min   100(x1 - x0^2)^2 +(1 - x0)^2
# s.t.  x0 x1 >= 1
#       x0 + x1^2 >= 0
#       x0 <= 0.5
#
# The standard start point(-2, 1) usually converges to the standard
# minimum at(0.5, 2.0), with final objective = 306.5.
# Sometimes the solver converges to another local minimum
# at(-0.79212, -1.26243), with final objective = 360.4.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO, MathOptInterface
using Test

const MOI = MathOptInterface

struct HS15 <: MOI.AbstractNLPEvaluator
    enable_hessian::Bool
end

MOI.features_available(d::HS15) = [:Grad, :Hess]
MOI.initialize(d::HS15, features) = nothing

MOI.jacobian_structure(d::HS15) = []
MOI.hessian_lagrangian_structure(d::HS15) = Tuple{Int64,Int64}[(1, 1), (1, 2), (2, 2)]

function MOI.eval_objective(d::HS15, x)
    return 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
end

function MOI.eval_constraint(::HS15, g, x)
    return nothing
end
function MOI.eval_constraint_jacobian(::HS15, jac_g, x)
    return nothing
end

function MOI.eval_objective_gradient(d::HS15, grad_f, x)
    # Evaluate gradient of nonlinear objective
    dTmp = x[2] - x[1] * x[1]
    grad_f[1] = (-400.0 * dTmp * x[1]) - (2.0 * (1.0 - x[1]))
    return grad_f[2] = 200.0 * dTmp
end

function MOI.eval_hessian_lagrangian(d::HS15, H, x, σ, μ)
    H[1] = σ * ((-400.0 * x[2]) + (1200.0 * x[1] * x[1]) + 2.0) #(0,0)
    H[2] = σ * (-400.0 * x[1])  #(0,1)
    H[3] = σ * 200.0            #(1,1)

    return 0
end

function example_moi_nlp1(; verbose=true)
    options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
    # Define a Knitro solver instance through MOI.
    solver = KNITRO.Optimizer(outlev=3, opttol=1e-8)
    MOI.set(solver, MOI.Silent(), !verbose)

    lb = []
    ub = []
    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), HS15(false), true)

    n = 2
    v = MOI.add_variables(solver, n)

    u = [0.5, Inf]
    start = [-2.0, 1.0]

    for i in 1:2
        MOI.add_constraint(solver, v[i], MOI.LessThan(u[i]))
        MOI.set(solver, MOI.VariablePrimalStart(), v[i], start[i])
    end

    # Add the constraints and set their lower bounds.
    m = 2
    # First constraint: x0 * x1 >= 1.
    cf1 = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(1.0, v[1], v[2])],
        MOI.ScalarAffineTerm.(0.0, v),
        0.0,
    )
    c1 = MOI.add_constraint(solver, cf1, MOI.GreaterThan{Float64}(1.0))

    # Second constraint: x0 + x1^2 >= 0.
    cf2 = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(10.0, v[2], v[2])],
        [MOI.ScalarAffineTerm(1.0, v[1])],
        0.0,
    )
    c2 = MOI.add_constraint(solver, cf2, MOI.GreaterThan(0.0))

    MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    # Define NLP structure.
    MOI.set(solver, MOI.NLPBlock(), block_data)

    # Solve the problem.
    MOI.optimize!(solver)

    # Free solver environment properly
    return MOI.empty!(solver)
end

example_moi_nlp1(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
using Test, Flux, NeuralNetDiffEq , StochasticDiffEq
using DiffEqDevTools
#Initial case for u(0 , x)
function phi(xi)
    y = Float64[]
    for x in xi
        y = push!(y , 1.77*x -0.015*x^3)
    end
    y = reshape(y , size(xi)[1] , size(xi)[2] )
    return y
end

xspan = (-6.0 , 6.0)
tspan = (0.0 , 1.0)
#f = mu and g  = sigma
g(u , p , t) = 0.5*u
f(u , p , t) = 0.5*0.25*u
d = 1
sdealg = EM()
prob = KolmogorovPDEProblem(f , g, phi , xspan , tspan, d)
opt = Flux.ADAM(0.01)
m = Chain(Dense(1, 512, elu) , Dense(512 , 1024 , elu),Dense(1024 , 512 , elu), Dense(512 , 1))
sol = solve(prob, NeuralNetDiffEq.NNKolmogorov(m,opt , sdealg), verbose = true, dt = 0.001,
            abstol=1e-6, maxiters = 250)

function analytical(xi)
    y = Float64[]
    a = 1.77*exp(0.5*(0.5)^2*1.0)
    b = -0.015*exp(0.5*(0.5*3)^2*1.0)
    for x in xi
        y = push!(y , a*x + b*x^3)
    end
    y = reshape(y , size(xi)[1] , size(xi)[2] )
    return y
end
xs = xspan[1]:0.001:xspan[2]
x_val = rand(xs , d , 50)
errorl2 = Flux.mse(analytical(x_val) , m(x_val))
println("error_l2 = ", errorl2, "\n")
@test errorl2 < 0.3


text:
module Murrell

using TDMSReader, XLSX, NoiseRobustDifferentiation, Dates
using DelimitedFiles, OldTools

export M_file, M_read, M_reduce!, Youngs_mod!, YS!, SFP_read, JR!, M_interp!

include("base.jl")

end

max_stars_repo_name: chrisharb/Murrell

max_issues_repo_name: chrisharb/Murrell

max_forks_repo_name: chrisharb/Murrell


text:
using OffsetArrays

@inline function _bernstein(j::Int, p::Int, x::T) where {T}
    if j < 0 || j > p
        return zero(T)
    else
        if p == 0
            return one(T)
        else
            return _bernstein(j, p-1, x) * (1-x) + _bernstein(j-1, p-1, x) * x
        end
    end
end

"""
Bernstein basis on the interval [0..1].
"""
struct Bernstein{T, BT} <: Basis{T}
    b::BT
    n::Int

    function Bernstein{T}(n::Integer) where {T}
        p = n-1
        b = OffsetArray([y -> _bernstein(i, p, y) for i in 0:p], 0:p)
        new{T, typeof(b)}(b, n)
    end

end

Bernstein(::Type{T}, n::Integer) where {T} = Bernstein{T}(n)
Bernstein(n::Integer) = Bernstein(Float64, n)

(B::Bernstein)(x::Number, j::Integer) = B.b[j](x)

basis(B::Bernstein) = B.b
nbasis(B::Bernstein) = B.n
order(B::Bernstein) = nbasis(B)
degree(B::Bernstein) = nbasis(B) - 1

Base.eltype(::Bernstein{T}) where {T} = T
Base.eachindex(B::Bernstein) = eachindex(B.b)
Base.axes(B::Bernstein) = (Inclusion(0..1), eachindex(B))

Base.hash(B::Bernstein, h::UInt) = hash(B.n, h)
Base.:(==)(B1::Bernstein, B2::Bernstein) = (B1.n == B2.n)
Base.isequal(B1::Bernstein{T1}, B2::Bernstein{T2}) where {T1,T2} = (T1 == T2 && B1 == B2)

Base.getindex(B::Bernstein, x::Number, j::Integer) = B(x,j)
Base.getindex(B::Bernstein, x::Number,  ::Colon) = [b(x) for b in B.b]
Base.getindex(B::Bernstein, X::AbstractVector, j::Integer) = B.(X,j)
Base.getindex(B::Bernstein, X::AbstractVector,  ::Colon) = [b(x) for x in X, b in B.b]

## Derivative

function _eval_derivative(b::Bernstein{BT}, x::DT, i::Int) where {BT,DT}
    @assert i ≥ 0 && i < b.n
    (b.n-1) * ( _bernstein(i-1, b.n-2, x) - _bernstein(i, b.n-2, x) )
end

@simplify *(D::Derivative, B::Bernstein) = Mul(D,B)

const BernsteinDerivative = QMul2{<:Derivative,<:Bernstein}

Base.getindex(D::BernsteinDerivative, x::Number, j::Integer) = _eval_derivative(D.B, x, j)
Base.getindex(D::BernsteinDerivative, x::Number,  ::Colon) = [_eval_derivative(D.B, x, j) for j in eachindex(D.B)]
Base.getindex(D::BernsteinDerivative, X::AbstractVector, j::Integer) = [_eval_derivative(D.B, x, j) for x in X]
Base.getindex(D::BernsteinDerivative, X::AbstractVector,  ::Colon) = [_eval_derivative(D.B, x, j) for x in X, j in eachindex(D.B)]

Base.adjoint(B::Bernstein) = Derivative(axes(B,1)) * B


text:
using FastParzenWindows
using Test

@testset "FastParzenWindows.jl" begin
    # Write your tests here.
end


text:
using CSV, Distributions, HypothesisTests

data = CSV.read("../data/machine1.csv", header=false)[:,1]
xBar, n = mean(data), length(data)
sig = 1.2
alpha = 0.1
z = quantile(Normal(),1-alpha/2)

println("Calculating formula: ", (xBar - z*sig/sqrt(n), xBar + z*sig/sqrt(n)))
println("Using confint() function: ", confint(OneSampleZTest(xBar,sig,n),alpha))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
function create_grid(Np, Npoin_cg, Ne, x, DFloat)

    coord_cg=zeros(DFloat,Npoin_cg)
    # intma_cg=zeros(Int64,Np,Ne)
    # periodicity_cg=zeros(Int64,Npoin_cg)

    #Set some constants
    #Always assume that the domain is divided in equispaced intervals
    xmin=-1
    xmax=+1
    dx=(xmax-xmin)/Ne

    #Generate Grid Points
    # I=1
    # coord_cg[1]=xmin
    # for e=1:Ne
    #     x0=xmin + (e-1)*dx
    #     for i=2:Np
    #         I+=1
    #         coord_cg[I]= x0 + 0.5 * ( x[i] + 1 ) * dx
    #     end
    # end

    grid_points = copy(coord_cg)
    nelements = Ne
    npoints = Np
    ξ = copy(x)
    x_min = xmin

    grid_points[1] = x_min
    for i in 1:nelements
      # calculate left element edge location
      # and index offset
      i_left = (i - 1) * (npoints - 1)
      x_left = (i - 1) * dx + x_min

      # loop of number of interpolation points
      # within each element starting at the second
      # point to avoid duplication
      for j in 2:npoints
        println((i, i_left, j))
        # calculate interpolation point location
        # relative to left element edge
        x_relative = 0.5 * (ξ[j] + 1) * dx
        grid_points[i_left + j] = x_left + x_relative
      end
    end

    coord_cg = copy(grid_points)

    # #Generate INTMA
    # I=1
    # for e=1:Ne
    #     intma_cg[1,e]=I
    #     for i=2:Np
    #         I+=1
    #         intma_cg[i,e]=I
    #     end
    # end

    # #Generate Periodicity_CG array
    # for I=1:Npoin_cg
    #     periodicity_cg[I]=I
    # end
    # periodicity_cg[Npoin_cg]=1

    return coord_cg#, intma_cg, periodicity_cg)
end

# Np = poly_order + 1 (dof)
# Npoin_cg = num_elements * poly_order + 1
# ne = num_elements
# x = interpolation points on reference element (-1,1) there should be P+1

# |. . . , . . . , . . .| = Ne * (P + 1) - (Ne-1) = Ne * P + Ne - Ne + 1 = Ne * P + 1

"""
  QuadratureStyle
"""
abstract type AbstractQuadratureStyle end

Base.@kwdef struct LegendreGaussLobatto{𝒯} <: AbstractQuadratureStyle
    polynomial_order::𝒯
end

function Base.getproperty(quadrature_style::AbstractQuadratureStyle, name::Symbol)
    if name == :npoints
        return quadrature_style.polynomial_order + 1
    elseif name == :interpolation_points
        return create_interpolation_points(quadrature_style)
    elseif name == :quadrature_weights
        return create_quadrature_weights(quadrature_style)
    else
        return getfield(quadrature_style, name)
    end
end

function Base.propertynames(::AbstractQuadratureStyle)
    return (:polynomial_order, :npoints, :interpolation_points, :quadrature_weights)
end

function create_interpolation_points(quadrature_style::LegendreGaussLobatto)
    polynomial_order = quadrature_style.polynomial_order
    npoints          = quadrature_style.npoints

    xgl=zeros(DFloat,npoints)
    p= polynomial_order #Order of the Polynomials
    ph=floor(typeof(npoints), (p+1)/2 )

    for i=1:ph
        x=cos( (2*i-1)*pi/(2*p+1) )
        for k=1:20
            (L0,L0_1,L0_2)=legendre_poly(p,x); #Compute Nth order Derivatives of Legendre Polys
            #Get new Newton Iteration
            dx=-(1-x^2)*L0_1/(-2*x*L0_1 + (1-x^2)*L0_2)
            x=x+dx
            if (abs(dx) < 1.0e-20)
                break
            end
        end
        xgl[p+2-i]=x
    end

    #Check for Zero Root
    if (p+1 != 2*ph)
        x=0
        (L0,L0_1,L0_2)=legendre_poly(p,x)
        xgl[ph+1]=x
    end

    #Find remainder of roots via symmetry
    for i=1:ph
        xgl[i]=-xgl[p+2-i]
    end

    return xgl
end

function create_quadrature_weights(quadrature_style::LegendreGaussLobatto)
    polynomial_order = quadrature_style.polynomial_order
    npoints          = quadrature_style.npoints

    wgl=zeros(DFloat,npoints)
    p=polynomial_order; #Order of the Polynomials
    ph=floor(typeof(npoints), (p+1)/2 )

    for i=1:ph
        x=cos( (2*i-1)*pi/(2*p+1) )
        for k=1:20
            (L0,L0_1,L0_2)=legendre_poly(p,x); #Compute Nth order Derivatives of Legendre Polys
            #Get new Newton Iteration
            dx=-(1-x^2)*L0_1/(-2*x*L0_1 + (1-x^2)*L0_2)
            x=x+dx
            if (abs(dx) < 1.0e-20)
                break
            end
        end
        wgl[p+2-i]=2/(p*(p+1)*L0^2)
    end

    #Check for Zero Root
    if (p+1 != 2*ph)
        x=0
        (L0,L0_1,L0_2)=legendre_poly(p,x)
        wgl[ph+1]=2/(p*(p+1)*L0^2)
    end

    #Find remainder of roots via symmetry
    for i=1:ph
        wgl[i]=+wgl[p+2-i]
    end

    return wgl
end

function legendre_poly(p,x)
    L1, L1_1, L1_2 = 0,0,0
    L0, L0_1, L0_2 = 1,0,0

    for i=1:p
        L2=L1;L2_1=L1_1;L2_2=L1_2
        L1=L0;L1_1=L0_1;L1_2=L0_2
        a=(2*i-1)/i
        b=(i-1)/i
        L0=a*x*L1 - b*L2
        L0_1=a*(L1 + x*L1_1) - b*L2_1
        L0_2=a*(2*L1_1 + x*L1_2) - b*L2_2
    end
    return (L0,L0_1,L0_2)
end

"""
  Domain
"""
abstract type AbstractDomain end

Base.@kwdef struct IntervalDomain{𝒯} <: AbstractDomain
    min::𝒯
    max::𝒯

    function IntervalDomain{𝒯}(min, max)
        @assert min < max
        new{𝒯}(min, max)
    end
end

Base.length(domain::IntervalDomain) = domain.max - domain.min

Base.size(domain::IntervalDomain) = length(domain)

"""
  Topology
  "Informally, a topology tells how elements of a set of points relate spatially to each other." (wikipedia.org)
"""
abstract type AbstractTopology end

Base.@kwdef struct ElementTopology{𝒯,𝒰} <: AbstractTopology
    domain::𝒯
    nelements::𝒰
end

function Base.getproperty(topology::ElementTopology, name::Symbol)
    if name == :nfaces
        return topology.nelements + 1
    elseif name == :interface_locations
        return create_interface_locations(topology.domain, topology.nfaces)
    else
        return getfield(topology, name)
    end
end

function Base.propertynames(::ElementTopology)
    return (:domain, :nelements, :nfaces, :interface_locations)
end

function create_interface_locations(domain::IntervalDomain, nfaces::Int)
    return Array(range(domain.min, domain.max, length = nfaces))
end

"""
  ContinuousStyle
"""
abstract type AbstractContinuityStyle end

struct ContinuousStyle <: AbstractContinuityStyle end

"""
  Grid
"""
abstract type AbstractGrid end

Base.@kwdef struct SpectralElementGrid{𝒯,𝒰,𝒱} <: AbstractGrid
    topology::𝒯
    quadrature_style::𝒰
    continuity_style::𝒱
end

function Base.getproperty(grid::SpectralElementGrid, name::Symbol)
    if name == :points
        return create_grid_points(grid.topology, grid.quadrature_style, grid.continuity_style)
    else
        return getfield(grid, name)
    end
end

function Base.propertynames(::SpectralElementGrid)
    return (:topology, :quadrature_style, :continuity_style, :points)
end

function create_grid_points(topology::𝒯, quadrature_style::𝒰, ::ContinuousStyle) where {𝒯,𝒰}
    # Unpack relevant parameters
    nelements = topology.nelements
    x_faces   = topology.interface_locations

    npoints   = quadrature_style.npoints
    ξ_points  = quadrature_style.interpolation_points

    grid_points[1] = x_min
    for i in 1:nelements
      # calculate left element edge location
      # and index offset
      dx     = x_faces[i+1] - x_faces[i]
      i_left = (i - 1) * (npoints - 1)
      x_left = (i - 1) * dx + x_min

      # loop of number of interpolation points
      # within each element starting at the second
      # point to avoid duplication
      for j in 2:npoints
        println((i, i_left, j))
        # calculate interpolation point location
        # relative to left element edge
        x_relative = 0.5 * (ξ_points[j] + 1) * dx
        grid_points[i_left + j] = x_left + x_relative
      end
    end

    return grid_points
end