text:
function cvd_plot(res::Array{Int64,2}; sird=false)
    l, _ = size(res)
    if sird
        plot(0:l - 1, res,
            label = ["Susceptible" "Infected" "Recovered" "Deceased"],
            xlabel = "Time",
            ylabel = "Count",
            linewidth = 2)
    else
        rec = res[:, 3] .+ res[:, 4]
        dat = hcat(res[:, 2], rec)
        plot(0:l - 1, dat,
            label = ["Infected" "Recovered"],
            xlabel = "Time",
            ylabel = "Count",
            linewidth = 2)
    end
end

function cvd_plot(res::Array{Float64,2}; sird=false)
    l, _ = size(res)
    p1 = if sird
        plot(0:l - 1, res[:, 1:4],
            label = ["Susceptible" "Infected" "Recovered" "Deceased"],
            ylabel = "Count",
            linewidth = 2)
    else
        rec = res[:, 3] .+ res[:, 4]
        dat = hcat(res[:, 2], rec)
        plot(0:l - 1, dat,
            label = ["Infected" "Recovered"],
            xlabel = "Time",
            ylabel = "Count",
            linewidth = 2)
    end
    p2 = plot(0:l - 1, res[:, 5], legend = false, xlabel = "Time", ylabel = "Weight", linewidth = 2)
    return plot(p1, p2, layout = (2, 1))
end

ext: jl


text:
# constants for use with Float64 values

const PRECISION         =    53
const BIAS              =  1023
const EXPO_MAX          =  1023
const EXPO_MIN          = -1074
const EXPO_LEAST_NORMAL = -1022
const SIGN_MASK         =  0x7ff0000000000000

const SPLITTER          =  134217729.0  # 2.0^27 + 1.0

# constant vector used int scaledSqureRoot
const sqrtGammaToM = [
    5.26013590154837350724098988288012866555033980282317e210, # gamma^-2/2 = (2^-350)^-2
    2.29349861599007151161082089530208694079656498916828e105, # gamma^-1/2 = (2^-350)^-1
    1.0,                                                      # gamma^0/2  = (2^-350)^0
    4.3601508761683463371878950054372442606045150559341e-106, # gamma^1/2  = (2^-350)^+1
    1.9010915662951598235150724058351031092648712063735e-211  # gamma^2/2  = (2^-350)^+2
];

# constants used in faithfulFinal()
const bound1    = 3.08469742733169167070816004443201143863168634177267e-179  # beta_lo^2 / eps^3 = 2^-593
const bound2    = 1.4396524142538228424993723224595141948383030778566e+163   # beta_hi^2 * eps^2 = 2^+542
const gamma     = 1.9010915662951598235150724058351031092648712063735e-211   # gamma             = 2^-700
const gammaToM1 = 5.26013590154837350724098988288012866555033980282317e+210  # gamma^-1          = 2^+700

# constants used in faithfulNorm()
const absMask = 0x7fffffffffffffff # mask that strips the sign of a Float64
const boundA  = 0x5450000000000000 # 2^+326 bound
const boundB  = 0x2890000000000000 # 2^-374 bound
const twoM700 = 0x1430000000000000 # 2^-700 scaling factor
const two700  = 0x6bb0000000000000 # 2^700 scaling factor
const splits  = 0x41a0000002000000 # 2^27 + 1


text:
# Julia wrapper for header: /usr/include/libavutil/hwcontext.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_hwdevice_find_type_by_name,
    av_hwdevice_get_type_name,
    av_hwdevice_iterate_types,
    av_hwdevice_ctx_alloc,
    av_hwdevice_ctx_init,
    av_hwdevice_ctx_create,
    av_hwdevice_ctx_create_derived,
    av_hwframe_ctx_alloc,
    av_hwframe_ctx_init,
    av_hwframe_get_buffer,
    av_hwframe_transfer_data,
    av_hwframe_transfer_get_formats,
    av_hwdevice_hwconfig_alloc,
    av_hwdevice_get_hwframe_constraints,
    av_hwframe_constraints_free,
    av_hwframe_map,
    av_hwframe_ctx_create_derived

function av_hwdevice_find_type_by_name(name)
    ccall((:av_hwdevice_find_type_by_name, libavutil), AVHWDeviceType, (Cstring,), name)
end

function av_hwdevice_get_type_name(_type::AVHWDeviceType)
    ccall((:av_hwdevice_get_type_name, libavutil), Cstring, (AVHWDeviceType,), _type)
end

function av_hwdevice_iterate_types(prev::AVHWDeviceType)
    ccall((:av_hwdevice_iterate_types, libavutil), AVHWDeviceType, (AVHWDeviceType,), prev)
end

function av_hwdevice_ctx_alloc(_type::AVHWDeviceType)
    ccall((:av_hwdevice_ctx_alloc, libavutil), Ptr{AVBufferRef}, (AVHWDeviceType,), _type)
end

function av_hwdevice_ctx_init(ref)
    ccall((:av_hwdevice_ctx_init, libavutil), Cint, (Ptr{AVBufferRef},), ref)
end

function av_hwdevice_ctx_create(device_ctx, _type::AVHWDeviceType, device, opts, flags::Integer)
    ccall((:av_hwdevice_ctx_create, libavutil), Cint, (Ptr{Ptr{AVBufferRef}}, AVHWDeviceType, Cstring, Ptr{AVDictionary}, Cint), device_ctx, _type, device, opts, flags)
end

function av_hwdevice_ctx_create_derived(dst_ctx, _type::AVHWDeviceType, src_ctx, flags::Integer)
    ccall((:av_hwdevice_ctx_create_derived, libavutil), Cint, (Ptr{Ptr{AVBufferRef}}, AVHWDeviceType, Ptr{AVBufferRef}, Cint), dst_ctx, _type, src_ctx, flags)
end

function av_hwframe_ctx_alloc(device_ctx)
    ccall((:av_hwframe_ctx_alloc, libavutil), Ptr{AVBufferRef}, (Ptr{AVBufferRef},), device_ctx)
end

function av_hwframe_ctx_init(ref)
    ccall((:av_hwframe_ctx_init, libavutil), Cint, (Ptr{AVBufferRef},), ref)
end

function av_hwframe_get_buffer(hwframe_ctx, frame, flags::Integer)
    ccall((:av_hwframe_get_buffer, libavutil), Cint, (Ptr{AVBufferRef}, Ptr{AVFrame}, Cint), hwframe_ctx, frame, flags)
end

function av_hwframe_transfer_data(dst, src, flags::Integer)
    ccall((:av_hwframe_transfer_data, libavutil), Cint, (Ptr{AVFrame}, Ptr{AVFrame}, Cint), dst, src, flags)
end

function av_hwframe_transfer_get_formats(hwframe_ctx, dir::AVHWFrameTransferDirection, formats, flags::Integer)
    ccall((:av_hwframe_transfer_get_formats, libavutil), Cint, (Ptr{AVBufferRef}, AVHWFrameTransferDirection, Ptr{Ptr{AVPixelFormat}}, Cint), hwframe_ctx, dir, formats, flags)
end

function av_hwdevice_hwconfig_alloc(device_ctx)
    ccall((:av_hwdevice_hwconfig_alloc, libavutil), Ptr{Cvoid}, (Ptr{AVBufferRef},), device_ctx)
end

function av_hwdevice_get_hwframe_constraints(ref, hwconfig)
    ccall((:av_hwdevice_get_hwframe_constraints, libavutil), Ptr{AVHWFramesConstraints}, (Ptr{AVBufferRef}, Ptr{Cvoid}), ref, hwconfig)
end

function av_hwframe_constraints_free(constraints)
    ccall((:av_hwframe_constraints_free, libavutil), Cvoid, (Ptr{Ptr{AVHWFramesConstraints}},), constraints)
end

function av_hwframe_map(dst, src, flags::Integer)
    ccall((:av_hwframe_map, libavutil), Cint, (Ptr{AVFrame}, Ptr{AVFrame}, Cint), dst, src, flags)
end

function av_hwframe_ctx_create_derived(derived_frame_ctx, format::AVPixelFormat, derived_device_ctx, source_frame_ctx, flags::Integer)
    ccall((:av_hwframe_ctx_create_derived, libavutil), Cint, (Ptr{Ptr{AVBufferRef}}, AVPixelFormat, Ptr{AVBufferRef}, Ptr{AVBufferRef}, Cint), derived_frame_ctx, format, derived_device_ctx, source_frame_ctx, flags)
end


text:
# ----------------------------------------------------
# --- Import modules
# ----------------------------------------------------
using DigitalComm
using Test
# ----------------------------------------------------
# --- Tests
# ----------------------------------------------------
println("Tests for multicarrier waveforms");

# ----------------------------------------------------
# --- Overall parameters
# ----------------------------------------------------
# --- Overall PHY parameters
nbSymb 			= 14;			  # --- Number of symbols (one frame)
nFFT 			= 1024;			  # --- Base FFT size
samplingFreq	= 15.36;		  # --- Frequency value (MHz)
snrVect			= (-10:30);
# --- Frequency allocation
allocatedSubcarriers= getLTEAlloc(nFFT);
nbSubcarriers	      = length(allocatedSubcarriers);

# ----------------------------------------------------
# --- Waveform contender
# ----------------------------------------------------
# --- Init OFDM structure
ofdm  = initOFDM(
				 nFFT,						        # --- nFFT                 : FFT size
				 72,						        # --- nCP                  : CP size
				 allocatedSubcarriers		        # --- allocatedSubcarriers : Subcarrier allocation
				 );
# --- Init SCFDMA structure
scfdma  = initSCFDMA(
					 nFFT,						        # --- nFFT                 : FFT size
					 72,						        # --- nCP                  : CP size
					 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
					 12;								# --- sizeDFT			   : DFT preprocessing size
					 );
# --- Init UF-OFDM structure
ufofdm  = initUFOFDM(
					 nFFT,					        # --- nFFT                 : FFT size
					 73,					        # --- L                    : Filter length (same size +1 due to conv)
					 allocatedSubcarriers,	        # --- allocatedSubcarriers : Subcarrier allocation
					 applyPD=1,				        # --- applyPD              : Do predistortion at Tx stage
					 attenuation=40,		        # --- attenuation          : Filter attenuation in dB
					 );
# --- Init BF-OFDM structure
bfofdm	= initBFOFDM(
					 32,				            # --- nFBMC                : PPN size (max number of carriers)
					 64,				            # --- nOFDM                : Precoder size (OFDM sizer)
					 3,				            	# --- K                    : Overlapping factor
					 9,					            # --- GI                   : CP size of precoder
					 0.5,				            # --- δ                    : compression factor
					 allocatedSubcarriers,          # --- allocatedSubcarriers : Subcarrier allocation
					 "gaussian",	            # --- filterName           : Pulse shape name
					 BT=0.36,				        # --- BT                   : Potential BT value for Gaussian
					 filterStopBand = 110,			# --- filterStopBand       : DC stopband value
					 fS=[],				            # --- fS                   : Potential frequency coefficient for FS filter
					 nFFT= 1024,		            # --- nFFT                 : associated FFT value in Rx
					 nCP= 72,			            # --- nCP                  : extended CP size
					 );
# --- Init WOLA-OFDM structure
wola  = initWOLA(
				 nFFT,						        # --- nFFT                 : FFT size
				 72,						        # --- nCP                  : CP size
				 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
				 "triangle",						# --- Window type @Tx side
				 20,								# --- Window size @Tx side
				 "triangle",						# --- Window type @Rx side
				 20,								# --- Window size @Rx side
				 );
fbmc  = initFBMC(
				 nFFT,						        # --- nFFT                 : FFT size
				 4,									# --- K					   : Overlapping factor
				 allocatedSubcarriers		        # --- allocatedSubcarriers : Subcarrier allocation
				 );

# ----------------------------------------------------
# --- Merging structures
# ----------------------------------------------------
# Create  a dictionnary to rule them all
waveforms 	= initWaveforms(ofdm,
							scfdma,
							ufofdm,
							bfofdm,
							wola,
							fbmc,
							);

for  (name,struc) in waveforms
	@testset  "$name" begin
		for mcs=[4,16,64,256]
			nbBits			      = nbSymb * nbSubcarriers * Int(log2(mcs));
			# --- Binary sequence
			bitSeq	      = genBitSequence(nbBits);
			# Mapping
			qamSeq		  = bitMappingQAM(mcs,bitSeq);
			# --- T/F matrix
			qamMat		  = reshape(qamSeq,nbSubcarriers,nbSymb);
			# --- Signal
			sigId		  = genSig(qamMat,struc);
			# --- Waveform demodulator
			qamDec	  = decodeSig(sigId,struc);
			# --- Binary demapper
			bitDec	  = bitDemappingQAM(mcs,qamDec[:]);
			# --- BER measure
			nbE	 = sum(xor.(bitDec,bitSeq));
			@test nbE == 0
		end
	end
end


text:
# Bogumił Kamiński, 2022

# Codes for chapter 11

# Code for section 11.1

# deserialization of source data frame

using DataFrames
using Serialization
walk = deserialize("walk.bin")

# Code for a note on conversion

x = [1.5]
x[1] = 1
x

# Code from section 11.1.1

Matrix(walk)
Matrix{Any}(walk)
Matrix{String}(walk)

using Plots
plot(walk)

plot(Matrix(walk); labels=["x" "y"] , legend=:topleft)

# Code from section 11.1.2

Tables.columntable(walk)

using BenchmarkTools
function mysum(table)
           s = 0 #A
           for v in table.x #B
               s += v
           end
           return s
       end
df = DataFrame(x=1:1_000_000);
@btime mysum($df)

tab = Tables.columntable(df);
@btime mysum($tab)

@code_warntype mysum(df)

@code_warntype mysum(tab)

typeof(tab)

function barrier_mysum2(x)
    s = 0
    for v in x
        s += v
    end
    return s
end
mysum2(table) = barrier_mysum2(table.x)
@btime mysum2($df)

df = DataFrame(a=[1, 1, 2], b=[1, 1, 2])
unique(df)

tab = Tables.columntable(df)
unique(tab)

# Code from section 11.1.3

Tables.rowtable(walk)

nti = Tables.namedtupleiterator(walk)
for v in nti
    println(v)
end

er = eachrow(walk)
er[1]
er[end]
ec = eachcol(walk)
ec[1]
ec[end]

identity.(eachcol(walk))

df = DataFrame(x=1:2, b=["a", "b"])
identity.(eachcol(df))

# Code from section 11.2

using CSV
raw_data = """
city,date,rainfall
Olecko,2020-11-16,2.9
Olecko,2020-11-17,4.1
Olecko,2020-11-19,4.3
Olecko,2020-11-20,2.0
Olecko,2020-11-21,0.6
Olecko,2020-11-22,1.0
Ełk,2020-11-16,3.9
Ełk,2020-11-19,1.2
Ełk,2020-11-20,2.0
Ełk,2020-11-22,2.0
""";
rainfall_df = CSV.read(IOBuffer(raw_data), DataFrame)

gdf_city = groupby(rainfall_df, "city")

gdf_city_date = groupby(rainfall_df, Not("rainfall"))

keys(gdf_city_date)

gk1 = keys(gdf_city_date)[1]
g1_t = Tuple(gk1)
g1_nt = NamedTuple(gk1)
g1_dict = Dict(gk1)

gdf_city_date[1]
gdf_city_date[gk1]
gdf_city_date[g1_t]
gdf_city_date[g1_nt]
gdf_city_date[g1_dict]

gdf_city[("Olecko",)]
gdf_city[(city="Olecko",)]

using BenchmarkTools
bench_df = DataFrame(id=1:10^8);
bench_gdf = groupby(bench_df, :id);
@btime groupby($bench_df, :id);
bench_i = 1_000_000;
bench_gk = keys(bench_gdf)[bench_i];
bench_t = Tuple(bench_gk);
bench_nt = NamedTuple(bench_gk);
bench_dict = Dict(bench_gk);
@btime $bench_gdf[$bench_i];
@btime $bench_gdf[$bench_gk];
@btime $bench_gdf[$bench_t];
@btime $bench_gdf[$bench_nt];
@btime $bench_gdf[$bench_dict];

gdf_city[[2, 1]]
gdf_city[[1]]

[nrow(df) for df in gdf_city]

for p in pairs(gdf_city)
    println(p)
end

Dict(key.city => nrow(df) for (key, df) in pairs(gdf_city))

combine(gdf_city, nrow)

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
function parse_block(block; failure=true, test_warn=true)

    bname = match(r"^\s*testcase\s+(\S+)\s+\{\s*$", block[1])[1]

    testset = """@testset "$bname" begin
            """
    ind = "    "
    for i in 2:length(block)-1
        line = strip(block[i])
        (isempty(line) || startswith(line, "//")) && continue
        command = parse_command(line; failure=failure, test_warn=test_warn)
        testset = """$testset
        $ind$command
        """
    end
    testset = """$testset
            end

            """
end

"""

This function parses a line into julia code, e.g.

```
add [1, 2] [1, 2] = [2, 4]
```

is parsed into
```
@test +(Interval(1, 2), Interval(1, 2)) === Interval(2, 4)
```
"""
function parse_command(line; failure=true, test_warn=true)
    # extract parts in line
    m = match(r"^(.+)=(.+);$", line)
    lhs = m[1]
    rhs = m[2]
    rhs = split(rhs, "signal")
    warn = length(rhs) > 1 ? rhs[2] : ""
    rhs = rhs[1]

    lhs = parse_lhs(lhs)
    rhs = parse_rhs(rhs)

    expr = build_expression(lhs, rhs)
    if failure
        try
            res = eval(Meta.parse(expr))
            command = res ? "@test $expr" : "@test_skip $expr"
        catch
            command = "@test_skip $expr"
        end
    else
        command = "@test $expr"
    end
    if test_warn
        command = isempty(warn) ? command : "@test_logs (:warn, ) $command"
    end
    return command
end

function parse_lhs(lhs)
    lhs = strip(lhs)
    m =  match(r"^(\S+) (.+)$", lhs)
    fname = m[1]
    args = m[2]

    #special case, input text
    fname == "b-textToInterval" && return "@interval($args)"
    fname == "d-textToInterval" && return "@decorated($args)"

    # input numbers
    args = replace(args, "infinity" => "Inf")
    args = replace(args, "X" => "x")
    if fname == "b-numsToInterval"
        args = join(split(args), ',')
        return "interval($args)"
    end

    if fname == "d-numsToInterval"
        args = join(split(args), ',')
        return "DecoratedInterval($args)"
    end

    # input intervals
    rx = r"\[([^\]]+)\](?:_(\w+))?"
    for m in eachmatch(rx, args)
        args = replace(args, m.match => parse_interval(m[1], m[2]))
    end
    args = replace(args, " " => ", ")
    args = replace(args, ",," => ",")
    args = replace(args, "{" => "[")
    args = replace(args, "}" => "]")
    return functions[fname](args)

end

function parse_rhs(rhs)
    rhs = strip(rhs)
    rhs = replace(rhs, "infinity" => "Inf")
    rhs = replace(rhs, "X" => "x")
    if '[' ∉ rhs # one or more scalar/bolean values separated by space
        return split(rhs)
    else # one or more intervals
        rx = r"\[([^\]]+)\](?:_(\w+))?"
        ivals = [parse_interval(m[1], m[2]; check=false) for m in eachmatch(rx, rhs)]
        return ivals
    end
end

function parse_interval(ival, dec; check=true)

    ival == "nai" && return "nai()"
    if ival == "entire"
        ival =  "entireinterval()"
    elseif ival == "empty"
        ival = "emptyinterval()"
    else
        ival = check ? "interval($ival)" : "Interval($ival)"
    end
    isnothing(dec) || (ival = "DecoratedInterval($ival, $dec)")
    return ival
end

function build_expression(lhs, rhs::AbstractString)
    rhs == "nai()" && return "isnai($lhs)"
    rhs == "NaN" && return "isnan($lhs)"
    return "$lhs === $rhs"
end

function build_expression(lhs, rhs::Vector)
    length(rhs) == 1 && return build_expression(lhs, rhs[1])
    expr = [build_expression(lhs*"[$i]", r) for (i, r) in enumerate(rhs)]
    return join(expr, " && ")
end


text:
# # Nonlocal diffusion on $\mathbb{S}^2$
# This example calculates the spectrum of the nonlocal diffusion operator:
# ```math
# \mathcal{L}_\delta u = \int_{\mathbb{S}^2} \rho_\delta(|\mathbf{x}-\mathbf{y}|)\left[u(\mathbf{x}) - u(\mathbf{y})\right] \,\mathrm{d}\Omega(\mathbf{y}),
# ```
# defined in Eq. (2) of
#
# R. M. Slevinsky, H. Montanelli, and Q. Du, [A spectral method for nonlocal diffusion operators on the sphere](https://doi.org/10.1016/j.jcp.2018.06.024), *J. Comp. Phys.*, **372**:893--911, 2018.
#
# In the above, $0<\delta<2$, $-1<\alpha<1$, and the kernel:
# ```math
# \rho_\delta(|\mathbf{x}-\mathbf{y}|) = \frac{4(1+\alpha)}{\pi \delta^{2+2\alpha}} \frac{\chi_{[0,\delta]}(|\mathbf{x}-\mathbf{y}|)}{|\mathbf{x}-\mathbf{y}|^{2-2\alpha}},
# ```
# where $\chi_I(\cdot)$ is the indicator function on the set $I$.
#
# This nonlocal operator is diagonalized by spherical harmonics:
# ```math
# \mathcal{L}_\delta Y_\ell^m(\mathbf{x}) = \lambda_\ell(\alpha, \delta) Y_\ell^m(\mathbf{x}),
# ```
# and its eigenfunctions are given by the generalized Funk--Hecke formula:
# ```math
# \lambda_\ell(\alpha, \delta) = \frac{(1+\alpha) 2^{2+\alpha}}{\delta^{2+2\alpha}}\int_{1-\delta^2/2}^1 \left[P_\ell(t)-1\right] (1-t)^{\alpha-1} \,\mathrm{d} t.
# ```
# In the paper, the authors use Clenshaw--Curtis quadrature and asymptotic evaluation of Legendre polynomials to achieve $\mathcal{O}(n^2\log n)$ complexity for the evaluation of the first $n$ eigenvalues. With a change of basis, this complexity can be reduced to $\mathcal{O}(n\log n)$.
#
# First, we represent:
# ```math
# P_n(t) - 1 = \sum_{j=0}^{n-1} \left[P_{j+1}(t) - P_j(t)\right] = -\sum_{j=0}^{n-1} (1-t) P_j^{(1,0)}(t).
# ```
# Then, we represent $P_j^{(1,0)}(t)$ with Jacobi polynomials $P_i^{(\alpha,0)}(t)$ and we integrate using [DLMF 18.9.16](https://dlmf.nist.gov/18.9.16):
# ```math
# \int_x^1 P_i^{(\alpha,0)}(t)(1-t)^\alpha\,\mathrm{d}t = \left\{ \begin{array}{cc} \frac{(1-x)^{\alpha+1}}{\alpha+1} & \mathrm{for~}i=0,\\ \frac{1}{2i}(1-x)^{\alpha+1}(1+x)P_{i-1}^{(\alpha+1,1)}(x), & \mathrm{for~}i>0.\end{array}\right.
# ```
# The code below implements this algorithm, making use of the Jacobi--Jacobi transform `plan_jac2jac`.
# For numerical stability, the conversion from Jacobi polynomials $P_j^{(1,0)}(t)$ to $P_i^{(\alpha,0)}(t)$ is divided into conversion from $P_j^{(1,0)}(t)$ to $P_k^{(0,0)}(t)$, before conversion from $P_k^{(0,0)}(t)$ to $P_i^{(\alpha,0)}(t)$.

using FastTransforms, LinearAlgebra

function oprec!(n::Integer, v::AbstractVector, alpha::Real, delta2::Real)
    if n > 0
        v[1] = 1
    end
    if n > 1
        v[2] = (4*alpha+8-(alpha+4)*delta2)/4
    end
    for i = 1:n-2
        v[i+2] = (((2*i+alpha+2)*(2*i+alpha+4)+alpha*(alpha+2))/(2*(i+1)*(2*i+alpha+2))*(2*i+alpha+3)/(i+alpha+3) - delta2/4*(2*i+alpha+3)/(i+1)*(2*i+alpha+4)/(i+alpha+3))*v[i+1] - (i+alpha+1)/(i+alpha+3)*(2*i+alpha+4)/(2*i+alpha+2)*v[i]
    end
    return v
end

function evaluate_lambda(n::Integer, alpha::T, delta::T) where T
    delta2 = delta*delta
    scl = (1+alpha)*(2-delta2/2)

    lambda = Vector{T}(undef, n)

    if n > 0
        lambda[1] = 0
    end
    if n > 1
        lambda[2] = -2
    end

    oprec!(n-2, view(lambda, 3:n), alpha, delta2)

    for i = 2:n-1
        lambda[i+1] *= -scl/(i-1)
    end

    p = plan_jac2jac(T, n-1, zero(T), zero(T), alpha, zero(T))

    lambda[2:end] .= p'lambda[2:end]

    for i = 2:n-1
        lambda[i+1] = ((2i-1)*lambda[i+1] + (i-1)*lambda[i])/i
    end

    for i = 2:n-1
        lambda[i+1] += lambda[i]
    end

    return lambda
end

# The spectrum in `Float64`:
lambda = evaluate_lambda(10, -0.5, 1.0)

# The spectrum in `BigFloat`:
lambdabf = evaluate_lambda(10, parse(BigFloat, "-0.5"), parse(BigFloat, "1.0"))

# The $\infty$-norm relative error:
norm(lambda-lambdabf, Inf)/norm(lambda, Inf)


text:
# Defining resolution
dx, dy = 0.5, 0.5

include.(["Geometry.jl", "Math.jl"])
using Plots

# Problem definition
k(_, y) = y >= 0.5 ? 2 : 1
g(x, _) = cbrt(x^2)

shape = [(-1, -1), (1, -1), (1, 0), (0, 0), (0, 1), (-1, 1)]
Γ = [g, 1-g, 0, 0, g, -1-g]
Γn, Γd = prepareEdge(shape, Γ)

# Base functions
e(x₀, y₀) = (x, y) -> (x₀ - dx, y₀ - dy) < (x, y) < (x₀ + dx, y₀ + dy) ? (1 - abs(x - x₀) / dx) * (1 - abs(y - y₀) / dy) : 0

L(v) = ∑([∫dl(k * gᵢ * v, Γᵢ) for (gᵢ, Γᵢ) in Γn])
B(u, v) = ∫∫dS(k * (∂x(u) * ∂x(v) + ∂y(u) * ∂y(v)), (-1, -1), (1, 1))

eᵢ = [e(x, y) for x ∈ -1:dx:1, y ∈ -1:dy:1 if inpolygon((x, y), shape) != 0 && !onEdge((x, y), Γd)]
wᵢ = [B(u, v) for u ∈ eᵢ, v ∈ eᵢ] \ [L(v) for v ∈ eᵢ]

f(x, y) = ∑(wᵢ .* [e(x, y) for e ∈ eᵢ])

# Drawing plot
X, Y = -1:0.01:1, -1:0.01:1
Z = [f(x, y) for y ∈ Y, x ∈ X]
heatmap(X, Y, Z, aspect_ratio = 1)

max_stars_repo_name: Pan-Maciek/FEM

max_issues_repo_name: Pan-Maciek/FEM

max_forks_repo_name: Pan-Maciek/FEM


text:
module MIDI

include("trackevent.jl")
include("midievent.jl")
include("metaevent.jl")
include("sysexevent.jl")
include("note.jl")
include("miditrack.jl")
include("midifile.jl")
include("constants.jl")
include("variablelength.jl")

end


text:
using BinaryBuilder

name = "GDAL"
version = v"3.0.4"

# Collection of sources required to build GDAL
sources = [
    ArchiveSource("https://github.com/OSGeo/gdal/releases/download/v$version/gdal-$version.tar.gz",
        "fc15d2b9107b250305a1e0bd8421dd9ec1ba7ac73421e4509267052995af5e83"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gdal-*/

if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
    # Apply patch to customise PROJ library
    atomic_patch -p1 "$WORKSPACE/srcdir/patches/configure_ac_proj_libs.patch"
    autoreconf -vi
    export PROJ_LIBS="proj_7_0"
elif [[ "${target}" == *-linux-* ]]; then
    # Make sure GEOS is linked against libstdc++
     atomic_patch -p1 "$WORKSPACE/srcdir/patches/geos-m4-extra-libs.patch"
     atomic_patch -p1 "$WORKSPACE/srcdir/patches/configure_ac_curl_libs.patch"
    export EXTRA_GEOS_LIBS="-lstdc++"
    export EXTRA_CURL_LIBS="-lstdc++"
    if [[ "${target}" == powerpc64le-* ]]; then
        atomic_patch -p1 "$WORKSPACE/srcdir/patches/sqlite3-m4-extra-libs.patch"
        export EXTRA_GEOS_LIBS="${EXTRA_GEOS_LIBS} -lm"
        export EXTRA_SQLITE3_LIBS="-lm"
        # libpthread and libldl are needed for libgdal, so let's always use them
        export LDFLAGS="-lpthread -ldl"
    fi
    autoreconf -vi
fi

# Clear out `.la` files since they're often wrong and screw us up
rm -f ${prefix}/lib/*.la

./configure --prefix=$prefix --host=$target \
    --with-geos=${bindir}/geos-config \
    --with-proj=$prefix \
    --with-libz=$prefix \
    --with-sqlite3=$prefix \
    --with-curl=${bindir}/curl-config \
    --with-openjpeg \
    --with-python=no \
    --enable-shared \
    --disable-static
# Make sure that some important libraries are found
grep "HAVE_GEOS='yes'" config.log
grep "HAVE_SQLITE='yes'" config.log
grep "CURL_SETTING='yes'" config.log

make -j${nproc}
make install
"""

platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libgdal", :libgdal),
    ExecutableProduct("gdal_contour", :gdal_contour_path),
    ExecutableProduct("gdal_grid", :gdal_grid_path),
    ExecutableProduct("gdal_rasterize", :gdal_rasterize_path),
    ExecutableProduct("gdal_translate", :gdal_translate_path),
    ExecutableProduct("gdaladdo", :gdaladdo_path),
    ExecutableProduct("gdalbuildvrt", :gdalbuildvrt_path),
    ExecutableProduct("gdaldem", :gdaldem_path),
    ExecutableProduct("gdalinfo", :gdalinfo_path),
    ExecutableProduct("gdallocationinfo", :gdallocationinfo_path),
    ExecutableProduct("gdalmanage", :gdalmanage_path),
    ExecutableProduct("gdalsrsinfo", :gdalsrsinfo_path),
    ExecutableProduct("gdaltindex", :gdaltindex_path),
    ExecutableProduct("gdaltransform", :gdaltransform_path),
    ExecutableProduct("gdalwarp", :gdalwarp_path),
    ExecutableProduct("nearblack", :nearblack_path),
    ExecutableProduct("ogr2ogr", :ogr2ogr_path),
    ExecutableProduct("ogrinfo", :ogrinfo_path),
    ExecutableProduct("ogrlineref", :ogrlineref_path),
    ExecutableProduct("ogrtindex", :ogrtindex_path),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GEOS_jll"),
    Dependency("PROJ_jll"),
    Dependency("Zlib_jll"),
    Dependency("SQLite_jll"),
    Dependency("LibCURL_jll"),
    Dependency("OpenJpeg_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")


text:
# Something it seems I need - from https://github.com/JuliaLang/julia/issues/32567
import Base: redirect_stdout
function redirect_both(f::Function, stdout_buf, stderr_buf)
    old_stdout = stdout
    old_stderr = stderr

    stdout_rd,stdout_wr = redirect_stdout()
    stderr_rd,stderr_wr = redirect_stderr()

    # Have to jump through hoops to clear the pipes and to not run into blocking
    # issues.
    transfer_stdout = @async begin
        while isopen(stdout_wr)
            write(stdout_buf, read(stdout_rd))
        end
    end
    transfer_stderr = @async begin
        while isopen(stderr_wr)
            write(stderr_buf, read(stderr_rd))
        end
    end

    try
        ret = f()

        flush(stdout_wr)
        flush(stderr_wr)
        Libc.flush_cstdio()
        flush(stdout)
        flush(stderr)
        close(stdout_wr)
        close(stderr_wr)
        close(stdout_rd)
        close(stderr_rd)

        wait(transfer_stdout)
        wait(transfer_stderr)

        return ret
    finally
        redirect_stdout(old_stdout)
        redirect_stderr(old_stderr)
    end
end

function PrettyTime(secs, sigfigs=2)
    time = secs*1e6

    factors = [1000, 1000, 60, 60, 24, Inf]
    cumfactors = cumprod(factors)
    names = ["μs", "ms", "secs", "mins", "hrs", "days"]
    ind = findfirst(time .< cumfactors)

    out = time
    if ind > 1
        out /= cumfactors[ind-1]
    end
    out = round(out, sigdigits=sigfigs)

    return string(out) * " " * names[ind]
end


text:
    function xgrid2c(ψ::Array{Complex{Float64},1},cinfo::ProjectedGPE.Cinfo)
        error("implement")
    end

    function xgrid2c(ψ::Array{Complex{Float64},2},x,y,cinfo::ProjectedGPE.Cinfo)
        Mx,My = cinfo.M
        dx = x[2]-x[1];dy = y[2]-y[1]
        Tx = eigmat(Mx,x)
        Ty = eigmat(My,y)
        return Tx'*(ψ*Ty)*dx*dy
    end

    function x2c(ψ::Array{Complex{Float64},3},cinfo::ProjectedGPE.Cinfo)
        error("implement")
    end


text:
using Core.Compiler: AbstractInterpreter, InferenceResult, InferenceParams, InferenceState, OptimizationParams, MethodInstance
using GPUCompiler: CodeCache, WorldView

struct EnzymeInterpeter <: AbstractInterpreter
    global_cache::CodeCache
    method_table::Union{Nothing,Core.MethodTable}

    # Cache of inference results for this particular interpreter
    local_cache::Vector{InferenceResult}
    # The world age we're working inside of
    world::UInt

    # Parameters for inference and optimization
    inf_params::InferenceParams
    opt_params::OptimizationParams

    function EnzymeInterpeter(cache::CodeCache, mt::Union{Nothing,Core.MethodTable}, world::UInt)
        @assert world <= Base.get_world_counter()

        return new(
            cache,
            mt,

            # Initially empty cache
            Vector{InferenceResult}(),

            # world age counter
            world,

            # parameters for inference and optimization
            InferenceParams(unoptimize_throw_blocks=false),
            VERSION >= v"1.8.0-DEV.486" ? OptimizationParams() :
                                          OptimizationParams(unoptimize_throw_blocks=false),
        )
    end
end

Core.Compiler.InferenceParams(interp::EnzymeInterpeter) = interp.inf_params
Core.Compiler.OptimizationParams(interp::EnzymeInterpeter) = interp.opt_params
Core.Compiler.get_world_counter(interp::EnzymeInterpeter) = interp.world
Core.Compiler.get_inference_cache(interp::EnzymeInterpeter) = interp.local_cache
Core.Compiler.code_cache(interp::EnzymeInterpeter) = WorldView(interp.global_cache, interp.world)

# No need to do any locking since we're not putting our results into the runtime cache
Core.Compiler.lock_mi_inference(interp::EnzymeInterpeter, mi::MethodInstance) = nothing
Core.Compiler.unlock_mi_inference(interp::EnzymeInterpeter, mi::MethodInstance) = nothing

function Core.Compiler.add_remark!(interp::EnzymeInterpeter, sv::InferenceState, msg)
end

Core.Compiler.may_optimize(interp::EnzymeInterpeter) = true
Core.Compiler.may_compress(interp::EnzymeInterpeter) = true
# From @aviatesk:
#     `may_discard_trees = true`` means a complicated (in terms of inlineability) source will be discarded,
#      but as far as I understand Enzyme wants "always inlining, except special cased functions",
#      so I guess we really don't want to discard sources?
Core.Compiler.may_discard_trees(interp::EnzymeInterpeter) = false
if VERSION >= v"1.7.0-DEV.577"
Core.Compiler.verbose_stmt_info(interp::EnzymeInterpeter) = false
end

if isdefined(Base.Experimental, Symbol("@overlay"))
Core.Compiler.method_table(interp::EnzymeInterpeter, sv::InferenceState) =
    Core.Compiler.OverlayMethodTable(interp.world, interp.method_table)
else
Core.Compiler.method_table(interp::EnzymeInterpeter, sv::InferenceState) =
    GPUCompiler.WorldOverlayMethodTable(interp.world)
end

function is_primitive_func(@nospecialize(TT))
    isa(TT, DataType) || return false
    ft = TT.parameters[1]
    if ft === typeof(Base.cbrt) || ft === typeof(Base.sin) || ft === typeof(Base.cos) ||
       ft === typeof(Base.tan) || ft === typeof(Base.exp) || ft === typeof(Base.log) ||
       ft === typeof(Base.asin) || ft === typeof(Base.tanh) || ft === typeof(Base.FastMath.tanh_fast) ||
       ft === typeof(Base.string) || ft === typeof(Base.sqrt)

        if TT <: Tuple{ft, Float32} || TT <: Tuple{ft, Float64} || TT <: Tuple{ft, Float16}
            return true
        end
    end
    if ft === typeof(Base.:^)
        if TT <: Tuple{ft, Float32, Float32} || TT <: Tuple{ft, Float64, Float64}
            return true
        end
        if TT <: Tuple{ft, Float32, <:Integer} || TT <: Tuple{ft, Float64, <:Integer}
            return true
        end
    end
    # FIXME(@wsmoses): For which types should we not inline?
    if ft === typeof(Base.wait) || ft === typeof(Base._wait) || ft === typeof(Base.enq_work) ||
       ft === typeof(Base.Threads.threadid) || ft == typeof(Base.Threads.nthreads) ||
       ft === typeof(Base.Threads.threading_run)
        return true
    end
    return false
end

# branch on https://github.com/JuliaLang/julia/pull/41328
@static if isdefined(Core.Compiler, :is_stmt_inline)

function Core.Compiler.inlining_policy(
    interp::EnzymeInterpeter, @nospecialize(src), stmt_flag::UInt8,
    mi::MethodInstance, argtypes::Vector{Any})

    if is_primitive_func(mi.specTypes)
        return nothing
    end

    return Base.@invoke Core.Compiler.inlining_policy(
        interp::AbstractInterpreter, src::Any, stmt_flag::UInt8,
        mi::MethodInstance, argtypes::Vector{Any})
end

elseif isdefined(Core.Compiler, :inlining_policy)

import Core.Compiler: InliningTodo, InliningState
enzyme_inlining_policy(@nospecialize(src)) = Core.Compiler.default_inlining_policy(src)
Core.Compiler.inlining_policy(::EnzymeInterpeter) = enzyme_inlining_policy
function Core.Compiler.resolve_todo(todo::InliningTodo, state::InliningState{S, T, <:typeof(enzyme_inlining_policy)}) where {S<:Union{Nothing, Core.Compiler.EdgeTracker}, T}
    mi = todo.mi
    if is_primitive_func(mi.specTypes)
        return Core.Compiler.compileable_specialization(state.et, todo.spec.match)
    end

    return Base.@invoke Core.Compiler.resolve_todo(
        todo::InliningTodo, state::InliningState)
end

end # @static if isdefined(Core.Compiler, :is_stmt_inline)


text:
module FewShotAnomalyDetection

using Flux
using NNlib
using Distributions
using SpecialFunctions
using Adapt
using Random
using LinearAlgebra
using FluxExtensions
using EvalCurves
using Pkg
using Statistics

include("utils.jl")
include("KNNmemory.jl")
include("bessel.jl")
include("svae.jl")
include("svaebase.jl")
include("svaetwocaps.jl")
include("svae_vamp.jl")
include("svae_vamp_means.jl")
include("vae.jl")
include("vae_vamp_wass.jl")
include("vae_ard.jl")
include("jacodeco.jl")

# if in("CuArrays",keys(Pkg.installed())) && in("CUDAnative",keys(Pkg.installed()))
#     using CuArrays
#     using CUDAnative
#     include("svae_gpu.jl")
# end

export SVAE, SVAEbase, SVAEtwocaps, SVAEvamp, SVAEvampmeans, VAE, VAE_ard, x_from_z, elbo_loss, elbo_loss_ard, decomposed_elbo_loss, wloss, decomposed_wloss, rloss, log_pxexpectedz, pz, log_pz, log_px, log_pz_jacobian_encoder, log_pz_jacobian_decoder, samplez, sampleVamp, zparams, zfromx, printing_wloss, mem_wloss, log_det_jacobian_encoder, log_det_jacobian_decoder, hsplit1softp, VAEvamp, add_labeled_centroid!, correctedjacodeco, manifoldz, jacodeco

end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
# This file was generated, do not modify it. # hide
# In this tutorial, we are exploring the application of Ridge and Lasso


text:
#=

    EstimateBinding.jl

    Perform inference using the deterministic model where the antibody binding
    rate, γ, is also estimated (this model is non-identifiable).

=#

using Model

using CSV, DataFrames, DataFramesMeta
using Distributions
using ForwardDiff
using JLD2
using LinearAlgebra
using Optim
using Roots
using StatsBase
using StatsPlots

##############################################################
## Run DeterministicModel.jl
##############################################################

include("DeterministicModel.jl")

##############################################################
## Fit ODE model to data
##############################################################

# Likelihood
function loglike_alt(θ)
    loglike = 0.0
    for t in T
        # Solve ODE model
        A,I = solve_alt_ode_model(t,θ[1:4])
        # Predicted unquenched signal
        m_unquenched = A + e / θ[5]
        # Predicted quenched signal
        m_quenched = I + (1 - η) * (A - I) + e / θ[5]
        # Data ~ N(model,σ²)
        loglike += loglikelihood(Normal(θ[5]*m_unquenched,σ),@subset(data_summary,:Time .== t,(!).(:Quenched)).Signal_Cy5)
        loglike += loglikelihood(Normal(θ[5]*m_quenched,σ),@subset(data_summary,:Time .== t,:Quenched).Signal_Cy5)
    end
    return loglike
end

# Obtain maximum likelihood estimate (just about meaningless, since model is non-identifiable)
θalt = optimize(p -> -loglike_alt(p),ones(5)).minimizer

##############################################################
## Profile (log) binding rate γ
##############################################################

log_γ = range(-1.0,5.0,length=50)
LLopt = [-optimize(p -> -loglike_alt([γ;p]),collect(θdet)).minimum for γ in exp.(log_γ)]

# Get 95% CI
pll = γ -> -optimize(p -> -loglike_alt([γ;p]),collect(θdet)).minimum - maximum(LLopt)
find_zero(γ -> pll(γ) + quantile(Chisq(1),0.95) / 2, [0.0,1.0])

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using Distances

using mjWrap
using MuJoCo

################################################## Experiment Centric Functions
## accepts model and trajectory data
## returns sum of trajectory rewards for each traj; saves per timestep reward
function rewardfunction(x::mjWrap.mjSet, tid::Int,
                        s0::AbstractVector{Float64}, # pass in states
                        s1::AbstractVector{Float64},
                        o0::AbstractVector{Float64}, # observations
                        o1::AbstractVector{Float64},
                        ctrl::AbstractVector{Float64})::Float64
   d = x.datas[tid]

   reward = 2
   reward -= 1*s1[1]         #(10*s1[1])^2
   reward -= 1*s1[2]         #(10*s1[2])^2
   reward -= 4*abs(1.1 - s1[3]) #(10*(1.1 - s1[3]))^2

   reward -= 1e-7*mjWrap.scaledcontrol(x.m, d, d.qfrc_actuator)
   #reward -= 1e-2*norm(ctrl)^2

   #reward -= 1e-7*norm(o1[4:end]-o0[4:end])^2 # min jerk
   #reward -= 1e-6*norm(o1-o0)^2 # min jerk

   return reward
end

function obsfunction(x::mjWrap.mjSet, tid::Int,
                     s::AbstractVector{Float64},
                     o::AbstractVector{Float64})
   #o[1:x.nq] .=
   #o .= s
   d = x.datas[tid]
   o .= d.qacc
end

## accepts model,
## default state,
## preallocated matrix of states (nqnv x numT)
function initfunction!(x::mjWrap.mjSet,
                       s::mjWrap.TrajSamples)
   mag = 0.2
   init_state = view( s.state, :, 1, : )
   for t=1:s.numT
      @. init_state[:,t] = s.s0 + rand() * mag - mag/2.0
      init_state[3,t] -= 0.1
   end
end

function ctrlfunction!(x::mjWrap.mjSet, ctrl::Array{Float64})
   randn!(ctrl)
end

# modelfunc should takes in time index (iteration index) to know change curric
function modelfunction!(x::mjWrap.mjSet, iter::Int=0)
end

## evaluation function: reward function shapes exploration, evaluation
## function defines success or failure
function evalfunction(x::mjWrap.mjSet, d::mjWrap.TrajSamples)
   return mean(d.state[3,end,:]) # x position at time T, for all numT
end

using ExpFunctions
myfuncs = FunctionSet(modelfunction!,
                      initfunction!,
                      ctrlfunction!,
                      obsfunction,
                      rewardfunction,
                      evalfunction)


text:
using Distributions
using JSON

mu = 0
s = 0.19
d = Logistic( mu,s )

x = [ -5, -2.5, 0, 2.5, 5 ]

dmgf(t) = mgf(d, t )
y = map( dmgf, x )
println( y )

data = Dict([
	("mu", mu),
	("s", s),
	("data", x),
	("expected", y)
])

outfile = open("./test/fixtures/number.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )

max_stars_repo_name: distributions-io/logistic-mgf

max_issues_repo_name: distributions-io/logistic-mgf

max_forks_repo_name: distributions-io/logistic-mgf


text:
using Pkg
pkg"precompile"
using DoubleFloats

using MPI
using Test

import MPI: mpiexec

# Code coverage command line options; must correspond to src/julia.h
# and src/ui/repl.c
const JL_LOG_NONE = 0
const JL_LOG_USER = 1
const JL_LOG_ALL = 2
const coverage_opts =
    Dict{Int, String}(JL_LOG_NONE => "none",
                      JL_LOG_USER => "user",
                      JL_LOG_ALL => "all")

# Files to run with mpiexec -n 1
singlefiles = ["test_spawn.jl"]

function runtests()
    nprocs = clamp(Sys.CPU_THREADS, 2, 4)
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir = dirname(@__FILE__)
    istest(f) = endswith(f, ".jl") && startswith(f, "test_")
    testfiles = sort(filter(istest, readdir(testdir)))
    extra_args = []
    @static if !Sys.iswindows()
        if occursin( "OpenRTE", read(`$mpiexec --version`, String))
            push!(extra_args,"--oversubscribe")
        end
    end

    nfail = 0
    printstyled("Running MPI.jl tests\n"; color=:white)

    for f in testfiles
        coverage_opt = coverage_opts[Base.JLOptions().code_coverage]
        if f ∈ singlefiles
            run(`$mpiexec $extra_args -n 1 $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
        else
            run(`$mpiexec $extra_args -n $nprocs $exename --code-coverage=$coverage_opt $(joinpath(testdir, f))`)
        end
        Base.with_output_color(:green,stdout) do io
            println(io,"\tSUCCESS: $f")
        end
    end
    return nfail
end

exit(runtests())


text:
"""
    kPointRecombinationCross(parent1::GEPGenotype, parent2::GEPGenotype,
                            gpExperimentInfo::GEPInfo, rng::Random.AbstractRNG,
                            k::Integer=1)

documentation
"""
function kPointRecombinationCross(parent1::GEPGenotype, parent2::GEPGenotype,
                                  gpExperimentInfo::GEPInfo, rng::Random.AbstractRNG,
                                  k::Integer=1)

    if k < 1
        error("k_points ($k) should be greater than 0")
    end

    parent1Rep = parent1._representation
    parent2Rep = parent2._representation
    genotypeLen = length(parent1Rep)
    if k > genotypeLen
        error("k_points ($k) should be minor than the size of the genoytpe ($genotypeLen)")
    end

    child1 = Array{Node}(undef, genotypeLen)
    child2 = Array{Node}(undef, genotypeLen)

    points = Array{Integer}(undef, k+2)
    points[1] = 1
    points[end] = genotypeLen

    for i=1:k
        points[i+1] = (rand(rng, UInt16) % genotypeLen) + 1
    end

    points = sort(points)

    for i=1:k
        indexes = points[i]:points[i+1]-1
        if isodd(i)
            child1[indexes] = parent1Rep[indexes]
            child2[indexes] = parent2Rep[indexes]
        else
            child1[indexes] = parent2Rep[indexes]
            child2[indexes] = parent1Rep[indexes]
        end
    end

    indexes = points[k+1]:points[k+2]
    if isodd(k)
        child1[indexes] = parent2Rep[indexes]
        child2[indexes] = parent1Rep[indexes]
    else
        child1[indexes] = parent1Rep[indexes]
        child2[indexes] = parent2Rep[indexes]
    end

    return [GEPGenotype(child1), GEPGenotype(child2)]
end

"""
    onePointRecombinationCross(parent1::GEPGenotype, parent2::GEPGenotype,
                               rng::Random.AbstractRNG)

documentation
"""
function onePointRecombinationCross(parent1::GEPGenotype, parent2::GEPGenotype,
                                    rng::Random.AbstractRNG)

    return kPointRecombinationCross(parent1, parent2, rng)
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
const ≅ = isequal

"""
    recode!(dest::AbstractArray, src::AbstractArray[, default::Any], pairs::Pair...)

Fill `dest` with elements from `src`, replacing those matching a key of `pairs`
with the corresponding value.

For each `Pair` in `pairs`, if the element is equal to (according to `isequal`))
the key (first item of the pair) or to one of its entries if it is a collection,
then the corresponding value (second item) is copied to `dest`.
If the element matches no key and `default` is not provided or `nothing`, it is copied as-is;
if `default` is specified, it is used in place of the original element.
`dest` and `src` must be of the same length, but not necessarily of the same type.
Elements of `src` as well as values from `pairs` will be `convert`ed when possible
on assignment.
If an element matches more than one key, the first match is used.

    recode!(dest::CategoricalArray, src::AbstractArray[, default::Any], pairs::Pair...)

If `dest` is a `CategoricalArray` then the ordering of resulting levels is determined
by the order of passed `pairs` and `default` will be the last level if provided.

    recode!(dest::AbstractArray, src::AbstractArray{>:Missing}[, default::Any], pairs::Pair...)

If `src` contains missing values, they are never replaced with `default`:
use `missing` in a pair to recode them.
"""
function recode! end

recode!(dest::AbstractArray, src::AbstractArray, pairs::Pair...) =
    recode!(dest, src, nothing, pairs...)
# To fix ambiguity
recode!(dest::CategoricalArray, src::AbstractArray, pairs::Pair...) =
    recode!(dest, src, nothing, pairs...)
recode!(dest::CategoricalArray, src::CategoricalArray, pairs::Pair...) =
    recode!(dest, src, nothing, pairs...)

function recode!(dest::AbstractArray{T}, src::AbstractArray, default::Any, pairs::Pair...) where {T}
    if length(dest) != length(src)
        throw(DimensionMismatch("dest and src must be of the same length (got $(length(dest)) and $(length(src)))"))
    end

    @inbounds for i in eachindex(dest, src)
        x = src[i]

        for j in 1:length(pairs)
            p = pairs[j]
            if ((isa(p.first, Union{AbstractArray, Tuple}) && any(x ≅ y for y in p.first)) ||
                x ≅ p.first)
                dest[i] = p.second
                @goto nextitem
            end
        end

        # Value not in any of the pairs
        if ismissing(x)
            eltype(dest) >: Missing ||
                throw(MissingException("missing value found, but dest does not support them: " *
                                       "recode them to a supported value"))
            dest[i] = missing
        elseif default isa Nothing
            try
                dest[i] = x
            catch err
                isa(err, MethodError) || rethrow(err)
                throw(ArgumentError("cannot `convert` value $(repr(x)) (of type $(typeof(x))) to type of recoded levels ($T). " *
                                    "This will happen with recode() when not all original levels are recoded " *
                                    "(i.e. some are preserved) and their type is incompatible with that of recoded levels."))
            end
        else
            dest[i] = default
        end

        @label nextitem
    end

    dest
end

function recode!(dest::CategoricalArray{T}, src::AbstractArray, default::Any, pairs::Pair...) where {T}
    if length(dest) != length(src)
        throw(DimensionMismatch("dest and src must be of the same length (got $(length(dest)) and $(length(src)))"))
    end

    vals = T[p.second for p in pairs]
    default !== nothing && push!(vals, default)

    levels!(dest.pool, filter!(!ismissing, unique(vals)))
    # In the absence of duplicated recoded values, we do not need to lookup the reference
    # for each pair in the loop, which is more efficient (with loop unswitching)
    dupvals = length(vals) != length(levels(dest.pool))

    drefs = dest.refs
    pairmap = [ismissing(v) ? 0 : get(dest.pool, v) for v in vals]
    defaultref = default === nothing || ismissing(default) ? 0 : get(dest.pool, default)
    @inbounds for i in eachindex(drefs, src)
        x = src[i]

        for j in 1:length(pairs)
            p = pairs[j]
            if ((isa(p.first, Union{AbstractArray, Tuple}) && any(x ≅ y for y in p.first)) ||
                x ≅ p.first)
                drefs[i] = dupvals ? pairmap[j] : j
                @goto nextitem
            end
        end

        # Value not in any of the pairs
        if ismissing(x)
            eltype(dest) >: Missing ||
                throw(MissingException("missing value found, but dest does not support them: " *
                                       "recode them to a supported value"))
            drefs[i] = 0
        elseif default === nothing
            try
                dest[i] = x # Need a dictionary lookup, and potentially adding a new level
            catch err
                isa(err, MethodError) || rethrow(err)
                throw(ArgumentError("cannot `convert` value $(repr(x)) (of type $(typeof(x))) to type of recoded levels ($T). " *
                                    "This will happen with recode() when not all original levels are recoded "*
                                    "(i.e. some are preserved) and their type is incompatible with that of recoded levels."))
            end
        else
            drefs[i] = defaultref
        end

        @label nextitem
    end

    # Put existing levels first, and sort them if possible
    # for consistency with CategoricalArray
    oldlevels = setdiff(levels(dest), vals)
    filter!(!ismissing, oldlevels)
    L = eltype(oldlevels)
    if Base.OrderStyle(L) isa Base.Ordered
        sort!(oldlevels)
    elseif hasmethod(isless, (L, L))
        # isless may throw an error, e.g. for AbstractArray{T} of unordered T
        try
            sort!(oldlevels)
        catch e
            e isa MethodError || rethrow(e)
        end
    end
    levels!(dest, union(oldlevels, levels(dest)))

    dest
end

function recode!(dest::CategoricalArray{T, N, R}, src::CategoricalArray,
                 default::Any, pairs::Pair...) where {T, N, R<:Integer}
    if length(dest) != length(src)
        throw(DimensionMismatch("dest and src must be of the same length " *
                                "(got $(length(dest)) and $(length(src)))"))
    end

    vals = T[p.second for p in pairs]
    if default === nothing
        srclevels = levels(src)

        # Remove recoded levels as they won't appear in result
        firsts = (p.first for p in pairs)
        keptlevels = Vector{T}(undef, 0)
        sizehint!(keptlevels, length(srclevels))

        for l in srclevels
            if !(any(x -> x ≅ l, firsts) ||
                 any(f -> isa(f, Union{AbstractArray, Tuple}) && any(l ≅ y for y in f), firsts))
                try
                    push!(keptlevels, l)
                catch err
                    isa(err, MethodError) || rethrow(err)
                    throw(ArgumentError("cannot `convert` value $(repr(l)) (of type $(typeof(l))) to type of recoded levels ($T). " *
                                        "This will happen with recode() when not all original levels are recoded " *
                                        "(i.e. some are preserved) and their type is incompatible with that of recoded levels."))
                end
            end
        end
        levs, ordered = mergelevels(isordered(src), keptlevels, filter!(!ismissing, unique(vals)))
    else
        push!(vals, default)
        levs = filter!(!ismissing, unique(vals))
        # The order of default cannot be determined
        ordered = false
    end

    srclevels = src.pool === dest.pool ? copy(levels(src.pool)) : levels(src.pool)
    if length(levs) > length(srclevels) && view(levs, 1:length(srclevels)) == srclevels
        levels!(dest.pool, levs)
    else
        dest.pool = CategoricalPool{nonmissingtype(T), R}(levs, isordered(dest))
    end

    drefs = dest.refs
    srefs = src.refs

    origmap = [get(dest.pool, v, 0) for v in srclevels]
    levelsmap = Vector{DefaultRefType}(undef, length(srclevels)+1)
    # For missing values (0 if no missing in pairs' keys)
    levelsmap[1] = 0
    for p in pairs
        if ((isa(p.first, Union{AbstractArray, Tuple}) && any(ismissing, p.first)) ||
            ismissing(p.first))
            levelsmap[1] = get(dest.pool, p.second)
            break
        end
    end
    pairmap = [ismissing(p.second) ? 0 : get(dest.pool, p.second) for p in pairs]
    # Preserving ordered property only makes sense if new order is consistent with previous one
    ordered && (ordered = issorted(pairmap))
    ordered!(dest, ordered)
    defaultref = default === nothing || ismissing(default) ? 0 : get(dest.pool, default)
    @inbounds for (i, l) in enumerate(srclevels)
        for j in 1:length(pairs)
            p = pairs[j]
            if ((isa(p.first, Union{AbstractArray, Tuple}) && any(l ≅ y for y in p.first)) ||
                l ≅ p.first)
                levelsmap[i+1] = pairmap[j]
                @goto nextitem
            end
        end

        # Value not in any of the pairs
        if default === nothing
            levelsmap[i+1] = origmap[i]
        else
            levelsmap[i+1] = defaultref
        end

        @label nextitem
    end

    @inbounds for i in eachindex(drefs)
        v = levelsmap[srefs[i]+1]
        if !(eltype(dest) >: Missing)
            v > 0 || throw(MissingException("missing value found, but dest does not support them: " *
                                            "recode them to a supported value"))
        end
        drefs[i] = v
    end

    dest
end

"""
    recode!(a::AbstractArray[, default::Any], pairs::Pair...)

Convenience function for in-place recoding, equivalent to `recode!(a, a, ...)`.

# Examples
```jldoctest
julia> using CategoricalArrays

julia> x = collect(1:10);

julia> recode!(x, 1=>100, 2:4=>0, [5; 9:10]=>-1);

julia> x
10-element Array{Int64,1}:
 100
   0
   0
   0
  -1
   6
   7
   8
  -1
  -1
```
"""
recode!(a::AbstractArray, default::Any, pairs::Pair...) =
    recode!(a, a, default, pairs...)
recode!(a::AbstractArray, pairs::Pair...) = recode!(a, a, nothing, pairs...)

cat_promote_valuetype(x::Pair{K, V}) where {K, V} = V
cat_promote_valuetype(x::Pair{K, V}, y::Pair...) where {K, V} =
    cat_promote_type(V, cat_promote_valuetype(y...))

keytype_hasmissing(x::Pair{K}) where {K} = K === Missing
keytype_hasmissing(x::Pair{K}, y::Pair...) where {K} = K === Missing || keytype_hasmissing(y...)

"""
    recode(a::AbstractArray[, default::Any], pairs::Pair...)

Return a copy of `a`, replacing elements matching a key of `pairs` with the corresponding value.
The type of the array is chosen so that it can
hold all recoded elements (but not necessarily original elements from `a`).

For each `Pair` in `pairs`, if the element is equal to (according to `isequal`)
or `in` the key (first item of the pair), then the corresponding value
(second item) is used.
If the element matches no key and `default` is not provided or `nothing`, it is copied as-is;
if `default` is specified, it is used in place of the original element.
If an element matches more than one key, the first match is used.

    recode(a::CategoricalArray[, default::Any], pairs::Pair...)

If `a` is a `CategoricalArray` then the ordering of resulting levels is determined
by the order of passed `pairs` and `default` will be the last level if provided.

# Examples
```jldoctest
julia> using CategoricalArrays

julia> recode(1:10, 1=>100, 2:4=>0, [5; 9:10]=>-1)
10-element Array{Int64,1}:
 100
   0
   0
   0
  -1
   6
   7
   8
  -1
  -1

```

     recode(a::AbstractArray{>:Missing}[, default::Any], pairs::Pair...)

If `a` contains missing values, they are never replaced with `default`:
use `missing` in a pair to recode them. If that's not the case, the returned array
will accept missing values.

# Examples
```jldoctest
julia> using CategoricalArrays

julia> recode(1:10, 1=>100, 2:4=>0, [5; 9:10]=>-1, 6=>missing)
10-element Array{Union{Missing, Int64},1}:
 100
   0
   0
   0
  -1
    missing
   7
   8
  -1
  -1

```
"""
function recode end

recode(a::AbstractArray, pairs::Pair...) = recode(a, nothing, pairs...)
# To fix ambiguity
recode(a::CategoricalArray, pairs::Pair...) = recode(a, nothing, pairs...)

function recode(a::AbstractArray, default::Any, pairs::Pair...)
    V = cat_promote_valuetype(pairs...)
    # T cannot take into account eltype(src), since we can't know
    # whether it matters at compile time (all levels recoded or not)
    # and using a wider type than necessary would be annoying
    T = default isa Nothing ? V : cat_promote_type(typeof(default), V)
    # Exception 1: if T === Missing and default not missing,
    # assume the caller wants to recode only some values to missing,
    # but accept original values
    if T === Missing && !isa(default, Missing)
        dest = Array{Union{eltype(a), Missing}}(undef, size(a))
    # Exception 2: if original array accepted missing values and missing does not appear
    # in one of the pairs' LHS, result must accept missing values
    elseif T >: Missing || default isa Missing || (eltype(a) >: Missing && !keytype_hasmissing(pairs...))
        dest = Array{Union{T, Missing}}(undef, size(a))
    else
        dest = Array{nonmissingtype(T)}(undef, size(a))
    end
    recode!(dest, a, default, pairs...)
end

function recode(a::CategoricalArray{S, N, R}, default::Any, pairs::Pair...) where {S, N, R}
    V = cat_promote_valuetype(pairs...)
    # T cannot take into account eltype(src), since we can't know
    # whether it matters at compile time (all levels recoded or not)
    # and using a wider type than necessary would be annoying
    T = default isa Nothing ? V : cat_promote_type(typeof(default), V)
    # Exception 1: if T === Missing and default not missing,
    # assume the caller wants to recode only some values to missing,
    # but accept original values
    # Example: recode(categorical([0,1]), 0=>missing)
    if T === Missing && !isa(default, Missing)
        dest = CategoricalArray{Union{S, Missing}, N, R}(undef, size(a))
    # Exception 2: if original array accepted missing values and missing does
    # not appear in one of the pairs' LHS, result must accept missing values
    # Example: recode(categorical([missing,1]), 1=>0)
    elseif T >: Missing || default isa Missing || (eltype(a) >: Missing && !keytype_hasmissing(pairs...))
        dest = CategoricalArray{Union{T, Missing}, N, R}(undef, size(a))
    else
        dest = CategoricalArray{nonmissingtype(T), N, R}(undef, size(a))
    end
    recode!(dest, a, default, pairs...)
end

function Base.replace(a::CategoricalArray{S, N, R}, pairs::Pair...) where {S, N, R}
    # Base.replace(a::Array, pairs::Pair...) uses a wider type promotion than
    # recode. It promotes the source type S with the replaced types T.
    T = cat_promote_valuetype(pairs...)
    # Exception: replacing missings
    # Example: replace(categorical([missing,1.5]), missing=>0)
    if keytype_hasmissing(pairs...)
        dest = CategoricalArray{cat_promote_type(nonmissingtype(S), T), N, R}(undef, size(a))
    else
        dest = CategoricalArray{cat_promote_type(S, T), N, R}(undef, size(a))
    end
    recode!(dest, a, nothing, pairs...)
end

if VERSION >= v"0.7.0-"
    Base.replace!(a::CategoricalArray, pairs::Pair...) = recode!(a, pairs...)
else
    replace!(a::CategoricalArray, pairs::Pair...) = recode!(a, pairs...)
end


text:
"Find out where different methods are."
function where(k)
    f = eval(k)
    a = (rand(), rand(2), rand(2,2), rand(2,2,2))
    for x in a
        try println(@which f(x)); end
        for y in a
            try println(@which f(x,y)); end
        end
    end
end

"See which scalar/array args f accepts."
function testargs(k)
    print("$k:")
    f = eval(k)
    try f(rand()); print(" (N,)"); end
    try f(rand(2)); print(" (V,)"); end
    try f(rand(2,2)); print(" (M,)"); end
    try f(rand(2,2,2)); print(" (T,)"); end
    try f(rand(),rand()); print(" (N,N)"); end
    try f(rand(),rand(2)); print(" (N,V)"); end
    try f(rand(),rand(2,2)); print(" (N,M)"); end
    try f(rand(),rand(2,2,2)); print(" (N,T)"); end
    try f(rand(2),rand()); print(" (V,N)"); end
    try f(rand(2),rand(2)); print(" (V,V)"); end
    try f(rand(2),rand(2,2)); print(" (V,M)"); end
    try f(rand(2),rand(2,2,2)); print(" (V,T)"); end
    try f(rand(2,2),rand()); print(" (M,N)"); end
    try f(rand(2,2),rand(2)); print(" (M,V)"); end
    try f(rand(2,2),rand(2,2)); print(" (M,M)"); end
    try f(rand(2,2),rand(2,2,2)); print(" (M,T)"); end
    try f(rand(2,2,2),rand()); print(" (T,N)"); end
    try f(rand(2,2,2),rand(2)); print(" (T,V)"); end
    try f(rand(2,2,2),rand(2,2)); print(" (T,M)"); end
    try f(rand(2,2,2),rand(2,2,2)); print(" (T,T)"); end
    println()
end

function test1arg(x; o...)
    for k in keys(math1arg)
        fx = eval(k)
        f = isa(x,Number) ? fx : a->sum(fx(a))
        (f===fx) || name(f,(:sum,k))
        @test check_grads(f, x; o...)
    end
end

function test2arg(x1,x2,flist=keys(math2arg); o...)
    for k in flist
        fx = eval(k)
        f = isa(x1,Number) && isa(x2,Number) ? fx : (a,b)->sum(fx(a,b))
        (f===fx) || name(f,(:sum,k))
        @test check_grads(f, x1, x2; o...)
    end
end

function test2arg_old(x1,x2;eps=1e-5)
    # First, second, or both arguments can be arrays
    # Grad only works for scalar output, so we'll sum the results
    psum(a)=(isa(a,Number) ? a : (s=0;for i=1:length(a);s+=a[i];end;s))
    for k in keys(math2arg)
        fx = eval(k)
        f = (x1,x2)->psum(fx(x1,x2))
        y = f(x1,x2)
        dbg((k,:x,x1,x2,:y,y))
        g1,g2 = grad(f,1),grad(f,2)
        dx = [g1(x1,x2),g2(x1,x2)]
        # Need to separate here based on array or not, maybe write a gcheck?
        Dx = [(f(x1+eps,x2)-f(x1-eps,x2))/2eps, Dx2 = (f(x1,x2+eps)-f(x1,x2-eps))/2eps]
        println((k,:diff,abs(dx-Dx),:dx,dx,:Dx,Dx))
        isapprox2(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function test2arg_scalar(eps=1e-5)
    for k in keys(math2arg)
        f = eval(k)
        g1,g2 = grad(f,1),grad(f,2)
        x1,x2 = rand(2)
        y = f(x1,x2)
        dbg((k,:x,x1,x2,:y,y))
        dx = [g1(x1,x2),g2(x1,x2)]
        Dx = [(f(x1+eps,x2)-f(x1-eps,x2))/2eps, Dx2 = (f(x1,x2+eps)-f(x1,x2-eps))/2eps]
        println((k,:diff,abs(dx-Dx),:dx,dx,:Dx,Dx))
        # @test_approx_eq_eps dx Dx 1e-6
        isapprox2(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function test1arg_scalar(eps=1e-5)
    for k in keys(math1arg)
        f = eval(k)
        g = grad(f)
        x = rand()
        y = f(x)
        dbg((k,:x,x,:y,y))
        dx = g(x)
        Dx = (f(x+eps)-f(x-eps))/2eps
        println((k,:diff,abs(dx-Dx),:dx,dx,:Dx,Dx))
        # @test_approx_eq_eps dx Dx 1e-6
        isapprox(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function test1arg_array(eps=1e-5)
    for k in keys(math1arg)
        fx = eval(k)
        f = x->(a=fx(x);a[1]+a[2])
        g = grad(f)
        x = rand(2,3)
        y = f(x)
        dbg((k,:x,x,:y,y))
        dx = g(x)
        Dx = similar(x)
        for i=1:length(x)
            xi=x[i]
            x[i]=xi+eps; f1=f(x)
            x[i]=xi-eps; f2=f(x)
            x[i]=xi
            Dx[i]=(f1-f2)/2eps
        end
        println((k,:maxdiff,maximum(abs(dx-Dx)),:avgdx,mean(abs(dx)))) # :dx,dx,:Dx,Dx))
        isapprox2(dx, Dx) || error("$k derivative is wrong: $((dx,Dx))")
    end
    return true
end

function isapprox2(x, y;
                   maxeps::Real = max(eps(eltype(x)), eps(eltype(y))),
                   rtol::Real=maxeps^(1/3), atol::Real=maxeps^(1/2))
    size(x) == size(y) || (warn("isapprox: $(size(x))!=$(size(y))"); return false)
    x = convert(Array, x)
    y = convert(Array, y)
    d = abs(x-y)
    s = abs(x)+abs(y)
    maximum(d - rtol * s) <= atol
end

# include("math1arg.jl")
# @test test1arg_scalar()
# @test test1arg_array()

function test1(x=(1.,2.))
    # foo(x)=sin(x[1])+cos(x[2])
    foo(x) = x[1]+x[2]
    goo = grad(foo)
    @show goo(x)
    # This does not work because goo is not scalar valued!
    # hoo = grad(goo)
    # @show hoo(x)
end

function test2()
    gsin = grad(sin)
    hsin = grad(gsin)
    #@show sin(1.0)
    #@show gsin(1.0)
    @show hsin(1.0)
end

function test3()
    foo2(x,y)=sin(x)+cos(y)
    goo2 = grad(foo2)
    goo22 = grad(foo2, 2)
    @show goo2(1,2)
    @show goo22(1,2)
end


text:
function norms(R, xroot, lroot, ϵroot)

    @info "Computing normalization constants..."

    # Coefficient depending on ϵ which appears in the normalization constant
    function coefnorm(l, ϵ)
        if ϵ==1
            return -(l+1)
        else
            return l
        end
    end

    norm = @. 1/√(R^3*2*sphericalbessel(xroot, lroot)^2*(1+coefnorm(lroot, ϵroot)/xroot))

    @info "Done!"

    return norm

end

max_stars_repo_name: avramescudana/QuarkDroplets

max_issues_repo_name: avramescudana/QuarkDroplets

max_forks_repo_name: avramescudana/QuarkDroplets


text:
using CSV, Primes, DataFrames
a=[1 3 7 9];

P=Int64.(zeros(1,11));
X=[1 2 3 4 5 6 7 8 9 10 11];

for i=1:4
    for k=0:30000
        for j=1:30000
            for l=1:11
                P[l]=a[i]+10*k+10*j*(l-1);
            end
            if [isprime(P[1]) isprime(P[2]) isprime(P[3]) isprime(P[4]) isprime(P[5]) isprime(P[6]) isprime(P[7]) isprime(P[8]) isprime(P[9]) isprime(P[10]) isprime(P[11])]==[Bool(1) Bool(1)  Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1) Bool(1)]
               X=[X;P];
            end
        end
    end
end

Y=DataFrame(X,:auto);
CSV.write("gt11.csv", Y);


text:
# # Helper functions
to_point2(a::Vector{<: Number}) = Point2f(a[1], a[2])

"""
    geo2basic(obj::GeoInterface)

Converts any GeoInterface object to the equivalent GeometryBasics, which can be plotted by Makie.
"""
function geo2basic(poly::GeoInterface.AbstractPolygon)
    polygon_coordinates = GeoInterface.coordinates(poly)
    linestrings = map(x-> to_point2.(x), polygon_coordinates)
    return GeometryBasics.Polygon(linestrings[1], linestrings[2:end])
end

function geo2basic(poly::Vector{Vector{Vector{Float64}}})
    linestrings = map(x-> to_point2.(x), poly)
    return GeometryBasics.Polygon(linestrings[1], linestrings[2:end])
end

function geo2basic(linestring::Vector{Vector{Float64}})
    return to_point2.(linestring)
end

function geo2basic(point::GeoInterface.AbstractPoint)
    return to_point2(GeoInterface.coordinates(point))
end

function geo2basic(multi_poly::GeoInterface.AbstractMultiPolygon)
    polygons = GeoInterface.coordinates(multi_poly)
    return geo2basic.(polygons)
end

geo2basic(feature::GeoInterface.Feature) = geo2basic(GeoInterface.geometry(feature))
function geo2basic(feature::GeoInterface.AbstractLineString)
    return GeometryBasics.LineString(geo2basic(GeoInterface.coordinates(feature)))
end

to_multipoly(poly::Polygon) = GeometryBasics.MultiPolygon([poly])
to_multipoly(mp::MultiPolygon) = mp
to_multipoly(any) = GeometryBasics.MultiPolygon(any)

# Only converts polygons and multipolygons
function geo2basic(fc::GeoInterface.AbstractFeatureCollection)
    return map(geo2basic, GeoInterface.features(fc))
end

function convert_arguments(P::Type{<: Union{Poly, Mesh}}, geom::GeoInterface.AbstractGeometry)
    return convert_arguments(P, geo2basic(geom))
end

function convert_arguments(P::Type{<:Poly}, geom::GeoInterface.AbstractFeatureCollection)
    return convert_arguments(P, to_multipoly.(geo2basic(geom)))
end

# set the default plot type for Vectors of polygons,
# so that they are plotted using the most efficient method!
plottype(::Vector{<: GeoInterface.AbstractMultiPolygon}) = Mesh
plottype(::Vector{<: GeoInterface.AbstractPolygon}) = Mesh


text:
struct Animal
    name::String
end

sayname(a::Animal) = print("My name is $(a.name).")

mydog = Animal("Tom")

sayname(mydog)


text:
################################################################################
#
#  Interface
#
################################################################################
@doc Markdown.doc"""
    absolute_automorphism_group(C::ClassField)

 Computes a generating set for the automorphisms of the
   number field corresponding to $C$. It assumes that the base field is normal.
   If "check" is true, the function checks if the extension is normal.
"""
function absolute_automorphism_group(C::ClassField, check::Bool = false)
  L = number_field(C)
  K = base_field(C)
  autK = automorphisms(K)
  @assert length(autK) == degree(K)
  if check
    @assert isnormal(C)
  end
  autK_gen = small_generating_set(autK)
  return absolute_automorphism_group(C, autK_gen)
end

function absolute_automorphism_group(C::ClassField, aut_gen_of_base_field::Array{NfToNfMor, 1})
  L = number_field(C)
  @vprint :ClassField 1 "Computing rel_auto "
  @vtime :ClassField 1 aut_L_rel = rel_auto(C)::Vector{NfRelNSToNfRelNSMor_nf_elem}
  if iscyclic(C) && length(aut_L_rel) > 1
    aut = aut_L_rel[1]
    for i = 2:length(aut_L_rel)
      aut *= aut_L_rel[i]
    end
    aut_L_rel = NfRelNSToNfRelNSMor_nf_elem[aut]
  end
  @vprint :ClassField 1 "Extending automorphisms"
  @vtime :ClassField 1 rel_extend = Hecke.new_extend_aut(C, aut_gen_of_base_field)
  rel_gens = vcat(aut_L_rel, rel_extend)::Vector{NfRelNSToNfRelNSMor_nf_elem}
  C.AbsAutGrpA = rel_gens
  return rel_gens::Vector{NfRelNSToNfRelNSMor_nf_elem}
end

function automorphism_groupQQ(C::ClassField)
  return rel_auto(C)
end

###############################################################################
#
#  Automorphisms of abelian extension
#
###############################################################################

function rel_auto_easy(A::ClassField_pp)

  # sqrt[n](a) -> zeta sqrt[n](a) on A.A
  #on A.K, the Kummer: sqrt[n](a) = gen(A.K) -> zeta gen(A.K)
  #we have the embedding A.A -> A.K : gen(A.A) -> A.pe
  M = sparse_matrix(base_field(A.K))
  b = A.K(1)
  push!(M, SRow(b))
  for i=2:degree(A)
    b *= A.pe
    push!(M, SRow(b))
  end
  tau = hom(A.K, A.K, A.bigK.zeta*gen(A.K), check = true)
  N = SRow(tau(A.pe))
  C = cyclotomic_extension(base_field(A), degree(A))
  Mk = _expand(M, pseudo_inv(C.mp[1]))
  Nk = _expand(N, pseudo_inv(C.mp[1]))
  s = solve(Mk, Nk) # will not work, matrix non-square...
  im = A.A()
  r = degree(C.Kr)
  for (i, c) = s
    setcoeff!(im, i-1, c)
  end
  return hom(A.A, A.A, im, check = false)

end

function rel_auto_intersect(A::ClassField_pp)

  # In the computation of the class field, I saved the
  # automorphisms of A.K over k.
  # Now, I have to search for the one that generates the Galois
  # group of the target field over k
  C = cyclotomic_extension(base_field(A), degree(A))
  if !isdefined(A, :AutG)
    _aut_A_over_k(C, A)
  end
  G, mG = snf(abelian_group(A.AutR))
  #Now, I restrict them to A.A
  M = sparse_matrix(base_field(A.K))
  b = A.K(1)
  push!(M, SRow(b))
  for i = 2:degree(A)
    b *= A.pe
    push!(M, SRow(b))
  end
  Mk = _expand(M, pseudo_inv(C.mp[1]))
  # One of the automorphisms must generate the group, so I check the order.
  for j = 1:ngens(G)
    if !divisible(G.snf[j], fmpz(degree(A)))
      continue
    end
    #Construct the automorphism
    gener = mG(G[j])
    elem = A.pe
    for i = 1:ncols(A.AutR)
      if !iszero(gener[i])
        for s = 1:Int(gener[i])
          elem = A.AutG[i](elem)
        end
      end
    end
    N = SRow(elem)
    Nk = _expand(N, pseudo_inv(C.mp[1]))
    s = solve(Mk, Nk) # will not work, matrix non-square...
    im = A.A()
    for (i, c) = s
      setcoeff!(im, i-1, c)
    end
    return hom(A.A, A.A, im, check = false)
  end
  error("I can't find the automorphism!")

end

function rel_auto_generic(A::ClassField_pp)
  K = A.A
  imgs = roots(K.pol, K)
  homs = morphism_type(K)[hom(K, K, x, check = false) for x in imgs]
  return small_generating_set(homs, *)[1]
end

function rel_auto(A::ClassField_pp)
  @assert isdefined(A, :A)
  if !isdefined(A, :K)
    return rel_auto_generic(A)
  elseif degree(A) == degree(A.K)
    #If the cyclotomic extension and the target field are linearly disjoint, it is easy.
    return rel_auto_easy(A)
  else
    #Tricky case
    return rel_auto_intersect(A)
  end
end

function rel_auto(A::ClassField)
  aut = Vector{morphism_type(NfRel{nf_elem})}(undef, length(A.cyc))
  for i = 1:length(aut)
    aut[i] = rel_auto(A.cyc[i])
  end
  K = number_field(A)
  g = gens(K)
  Aut = Vector{morphism_type(K)}(undef, length(aut))
  for i = 1:length(aut)
    imgs = Vector{NfRelNSElem{nf_elem}}(undef, length(aut))
    for j = 1:length(imgs)
      if j == i
        imgs[j] = image_primitive_element(aut[i]).data(g[j])
      else
        imgs[j] = g[j]
      end
    end
    Aut[i] = hom(K, K, imgs)
  end
  return Aut
end

###############################################################################
#
#  Extension of automorphisms from the base field to the class field
#
###############################################################################

@doc Markdown.doc"""
    extend_to_cyclotomic(C::CyclotomicExt, tau::NfToNfMor) -> NfRelToNfRelMor

Given a cyclotomic extension $C$ of a number field $K$ and an automorphism $\tau$ of $K$,
  computes an extension of $\tau$ to $C$.

"""
function extend_to_cyclotomic(C::CyclotomicExt, tau::NfToNfMor)
  K = domain(tau)
  @assert K == base_field(C.Kr)
  gKr = gen(C.Kr)
  if euler_phi(C.n) == degree(C.Kr)
    #The extension with the roots of unity is disjoint from K
    #Therefore, the minimal polynomial has coefficient over QQ.
    return hom(C.Kr, C.Kr, tau, gKr)
  end
  g = C.Kr.pol
  tau_g = parent(g)([tau(coeff(g, i)) for i=0:degree(g)])
  i = 1
  z = deepcopy(gKr)
  while gcd(i, C.n) != 1 || !iszero(tau_g(z))
    i += 1
    mul!(z, z, gKr)
  end
  return hom(C.Kr, C.Kr, tau, z)

end

function new_extend_aut(A::ClassField, auto::T) where T <: Map
  return new_extend_aut(A, T[auto])[1]
end

function new_extend_aut(A::ClassField, autos::Array{T, 1}) where T <: Map
  # tau: k -> k
  k = domain(autos[1])
  @assert k == codomain(autos[1])
  @assert k == base_field(A)
  lp = factor(fmpz(degree(A)))
  L = number_field(A)
  # I call number field because to extend the automorphism I need the defining polynomials
  all_imgs = Array{Array{NfRelNSElem{nf_elem}, 1}, 1}(undef, length(autos))
  #Initialize the array
  for i=1:length(autos)
    all_imgs[i] = Vector{NfRelNSElem{nf_elem}}(undef, length(A.cyc))#[L() for i=1:length(A.cyc)]
  end
  lG = gens(L)
  #P-Sylow subgroups are invariant, I can reduce to the prime power case.
  res = Array{NfRelNSToNfRelNSMor_nf_elem, 1}(undef, length(autos))
  for (p, v) = lp.fac
    @vprint :ClassField 1 "Extending auto pp ..."
    @vtime :ClassField 1 imgs = extend_aut_pp(A, autos, p)
    # The output are the images of the cyclic components in A.A
    indices = Array{Int, 1}(undef, length(imgs[1]))
    j = 1
    for i = 1:length(imgs[1])
      while degree(A.cyc[j]) % Int(p) != 0
        j += 1
      end
      indices[i] =  j
      j += 1
    end
    #I need to embed Ap in L
    Ap = parent(imgs[1][1])
    emb = hom(Ap, L, NfRelNSElem{nf_elem}[lG[indices[i]] for i = 1:length(indices)])
    for j = 1:length(autos)
      for i = 1:length(imgs[j])
        all_imgs[j][indices[i]] = emb(imgs[j][i])
      end
    end
  end
  for i = 1:length(res)
    res[i] = hom(L, L, autos[i], all_imgs[i])
    #@hassert :NfOrd 1 isconsistent(res[i])
  end
  return res

end

################################################################################
#
#  Frobenius generating the automorphisms of the Kummer extension
#
################################################################################

#Find a prime ideal P such that the Frobenius generates the Galois group of the extension.
function find_frob(A::ClassField_pp, K::KummerExt, emb::NfToNfMor)

  m = defining_modulus(A)[1]
  d = A.o
  K1 = kummer_extension(d, [A.a])
  k1 = base_field(K1)
  O1 = maximal_order(k1)
  k = base_field(K)
  O = maximal_order(k)
  Sp = Hecke.PrimesSet(200, -1)
  cp = lcm([discriminant(O), minimum(m, copy = false), index(O), index(O1)])
  P = ideal(O, 1)
  for p in Sp
    if cp % p == 0
      continue
    end
    lp = prime_decomposition(O1, p)
    for i = 1:length(lp)
      try
        z = canonical_frobenius(lp[i][1], K1)
        if order(z) != d
          continue
        end
        lP = prime_decomposition_nonindex(emb, lp[i][1])
        P = lP[1][1]
        zK = canonical_frobenius(P, K)
      catch e
        if typeof(e) != BadPrime
          rethrow(e)
        end
        continue
      end
      return P
    end
  end
  error("Something strange is happening")
end

function find_frob(A::ClassField_pp)

  m = defining_modulus(A)[1]
  d = A.o
  K1 = kummer_extension(d, [A.a])
  k1 = base_field(K1)
  O = maximal_order(k1)
  Sp = Hecke.PrimesSet(200, -1)
  cp = lcm([minimum(m), index(O), discriminant(O)])
  for p in Sp
    if cp % p == 0
      continue
    end
    lp = prime_decomposition(O, p)
    for i = 1:length(lp)
      try
        zK1 = canonical_frobenius(lp[i][1], K1)
        if order(zK1) != d
          continue
        end
      catch e
        if typeof(e) != BadPrime
          rethrow(e)
        end
        continue
      end
      return lp[i][1]
    end
  end
  error("Something strange is happening")
end

#Finds prime such that the Frobenius automorphisms generate the automorphism group of the kummer extension
function find_gens(KK::KummerExt, gens_imgs::Array{Array{FacElem{nf_elem, AnticNumberField}, 1}, 1}, coprime_to::fmpz, idx::Int = 1)
  K = base_field(KK)
  O = maximal_order(K)
  els = GrpAbFinGenElem[]
  Q, mQ = quo(KK.AutG, els, false)
  s, ms = snf(Q)
  Sp = Hecke.PrimesSet(1000, -1)
  cp = lcm(discriminant(O), coprime_to)
  frob_gens = NfOrdIdl[]
  for q in Sp
    if cp % q == 0
      continue
    end
    lp = prime_decomposition(O, q)
    for i = 1:length(lp)
      try
        z = canonical_frobenius(lp[i][1], KK)
        el_in_quo = ms\(mQ(z))
        if iszero(el_in_quo)
          continue
        end
        found = false
        for i = 1:ngens(s)
          if iscoprime(s.snf[i], el_in_quo[i])
            found = true
            break
          end
        end
        if !found
          continue
        end
        for x in gens_imgs
          for y in x
            canonical_frobenius(lp[i][1], KK, y)
          end
        end
        push!(frob_gens, lp[i][1])
        push!(els, z)
      catch e
        if typeof(e) != BadPrime
          rethrow(e)
        end
        continue
      end
      Q, mQ = quo(KK.AutG, els, false)
      s, ms = snf(Q)
      if order(s) == idx
        break
      end
    end
    if order(s) == idx
      break
    end
  end
  return frob_gens
end

#extension of automorphisms in the case of extensions of exponent 2
function extend_aut2(A::ClassField, autos::Array{NfToNfMor, 1})
  Cp = [x for x in A.cyc if degree(x) % 2 == 0]
  autos_extended = Vector{Vector{NfRelNSElem{nf_elem}}}(undef, length(autos))
  AA, gAA = number_field([c.A.pol for c in Cp], check = false)
  if length(Cp) == 1
    for i = 1:length(autos)
      fl, el = ispower(autos[i](Cp[1].a)*inv(Cp[1].a), 2)
      @assert fl
      autos_extended[i] = NfRelNSElem{nf_elem}[AA(evaluate(el))*gAA[1]]
    end
    return autos_extended
  end
  KK = kummer_extension([2 for i = 1:length(Cp)], [x.a for x in Cp])
  act_on_gens = Array{Array{FacElem{nf_elem, AnticNumberField}, 1}, 1}(undef, length(KK.gen))
  for i = 1:length(KK.gen)
    act_on_gen_i = Array{FacElem{nf_elem, AnticNumberField}, 1}(undef, length(autos))
    for j = 1:length(autos)
      act_on_gen_i[j] = FacElem(Dict(autos[j](ke) => v for (ke, v) in KK.gen[i]))
    end
    act_on_gens[i] = act_on_gen_i
  end
  frob_gens = find_gens(KK, act_on_gens, minimum(defining_modulus(A)[1]))
  #I will compute a possible image cyclic component by cyclic component
  for w = 1:length(autos)
    images_KK = Array{Tuple{GrpAbFinGenElem, FacElem{nf_elem, AnticNumberField}}, 1}(undef, length(Cp))
    for i = 1:length(Cp)
      fl, coord, rt = _find_embedding(KK, act_on_gens[i][w], 2, frob_gens)
      @assert fl
      images_KK[i] = (coord, rt)
    end

    #Now, I can define the automorphism on AA
    images_K = Array{NfRelNSElem{nf_elem}, 1}(undef, length(images_KK))
    for i = 1:length(images_K)
      s = AA(evaluate(images_KK[i][2]))
      for j = 1:length(Cp)
        mul!(s, s, gAA[j]^Int(images_KK[i][1][j]))
      end
      images_K[i] = s
    end
    autos_extended[w] = images_K
  end
  return autos_extended

end

function extend_generic(A::ClassField, autos::Vector{NfToNfMor}, p::fmpz)
  Cp = [x1 for x1 in A.cyc if degree(x1) % Int(p) == 0]
  A, gA = number_field([c.A.pol for c in Cp], check = false)
  rts = Vector{Vector{NfRelNSElem{nf_elem}}}(undef, length(autos))
  for i = 1:length(autos)
    imgs = Vector{NfRelNSElem{nf_elem}}(undef, length(Cp))
    for j = 1:length(gA)
      pol = map_coefficients(autos[i], Cp[j].A.pol)
      imgs[j] = roots(pol, A)[1]
    end
    rts[i] = imgs
  end
  return rts
end

function check_disjoint_cyclotomic(A::ClassField, p::fmpz)
  e = ppio(fmpz(exponent(A)), p)[1]
  K = base_field(A)
  mr = A.rayclassgroupmap
  mq = A.quotientmap
  x = PolynomialRing(FlintZZ, "x")[2]
  f = cyclotomic(Int(e), x)
  fK = map_coefficients(K, f)
  s, ms = norm_group(fK, mr, false, cached = false)
  mp = ms*mq
  i, mi = image(mp)
  return Int(divexact(order(codomain(mq)), order(i)))
end

function extend_aut_pp(A::ClassField, autos::Array{NfToNfMor, 1}, p::fmpz)
  Cp = [x1 for x1 in A.cyc if degree(x1) % Int(p) == 0]
  if !all(x -> isdefined(x, :a), Cp)
    return extend_generic(A, autos, p)
  end
  d = maximum(degree(x) for x in Cp)
  if d == 2
    return extend_aut2(A, autos)
  end

  m = minimum(defining_modulus(A)[1])
  ind_image = 1
  if !isone(gcd(d, m)) && d != minimum(degree(x) for x in Cp)
    #Difficult case. First, we check that the extension and
    #the cyclotomic extension are disjoint
    ind_image = check_disjoint_cyclotomic(A, p)
    if !isone(ind_image)
      return extend_autos_hard_case(A, autos, p, ind_image)
    end
  end

  AA, gAA = number_field([c.A.pol for c in Cp], check = false)
  #Main Idea: I extend tau to the big kummer extension KK and then I restrict it to AA.
  k = base_field(A)
  C = cyclotomic_extension(k, d)
  KC = absolute_simple_field(C)
  # C is the base field of the kummer extension generated
  # by all the cyclic components.
  # I extend the automorphisms to C

  Autos_abs = Array{NfToNfMor, 1}(undef, length(autos))
  for i = 1:length(autos)
    aut = extend_to_cyclotomic(C, autos[i])
    Autos_abs[i] = hom(KC, KC, C.mp[1]\(aut(C.mp[1](gen(KC)))), check = false)
  end

  #I compute the embeddings of the small cyclotomic extensions into the others
  abs_emb = Array{NfToNfMor, 1}(undef, length(Cp))
  for i = 1:length(Cp)
    dCp = degree(Cp[i])
    if dCp == d
      abs_emb[i] = id_hom(KC)
    else
      Cs = cyclotomic_extension(k, dCp)
      emb = hom(Cs.Kr, C.Kr, gen(C.Kr)^div(d, dCp), check = false)
      img = C.mp[1]\(emb(Cs.mp[1](gen(Cs.Ka))))
      abs_emb[i] = hom(Cs.Ka, KC, img, check = false)
    end
  end

  #Now, I can compute the corresponding Kummer extension over the big cyclotomic field.
  exps = Array{Int, 1}(undef, length(Cp))
  gens = Array{FacElem{nf_elem, AnticNumberField}, 1}(undef, length(Cp))
  for i = 1:length(Cp)
    if degree(Cp[i]) == d
      gens[i] = Cp[i].a
      exps[i] = Cp[i].o
    else
      D = Dict{nf_elem, fmpz}()
      for (ke,v) in Cp[i].a
        D[abs_emb[i](ke)] = v
      end
      a = FacElem(D)
      exps[i] = Cp[i].o
      gens[i] = a
    end
  end
  KK = kummer_extension(exps, gens)
  K, gK = number_field(KK)
  #I need the inclusions of the single extensions Cp[i].K in K
  incs = Array{NfRelToNfRelNSMor_nf_elem, 1}(undef, length(Cp))
  for i = 1:length(Cp)
    incs[i] = hom(Cp[i].K, K, abs_emb[i], gK[i])
  end

  # I want extend the automorphisms to KK
  # First, I find a set of primes such that their Frobenius generates the Galois group of KK
  act_on_gens = Array{Array{FacElem{nf_elem, AnticNumberField}, 1}, 1}(undef, length(KK.gen))
  for i = 1:length(KK.gen)
    act_on_gen_i = Array{FacElem{nf_elem, AnticNumberField}, 1}(undef, length(autos))
    for j = 1:length(autos)
      D1 = Dict{nf_elem, fmpz}()
      for (ke, v) in KK.gen[i]
        D1[Autos_abs[j](ke)] = v
      end
      act_on_gen_i[j] = FacElem(D1)
    end
    act_on_gens[i] = act_on_gen_i
  end
  frob_gens = find_gens(KK, act_on_gens, minimum(defining_modulus(A)[1]))

  autos_extended = Array{morphism_type(K, K), 1}(undef, length(autos))
  #I will compute a possible image cyclic component by cyclic component
  for w = 1:length(autos)
    images_KK = Array{Tuple{GrpAbFinGenElem, FacElem{nf_elem, AnticNumberField}}, 1}(undef, length(Cp))
    for i = 1:length(Cp)
      fl, coord_img_emb, rt_img_emb = _find_embedding(KK, act_on_gens[i][w], Int(order(KK.AutG[i])), frob_gens)
      @assert fl
      images_KK[i] = (coord_img_emb, rt_img_emb)
    end

    #Now, I can define the automorphism on K
    images_K = Array{NfRelNSElem{nf_elem}, 1}(undef, length(images_KK))
    for i = 1:length(images_K)
      s = K(evaluate(images_KK[i][2]))
      for j = 1:length(Cp)
        mul!(s, s, gK[j]^Int(images_KK[i][1][j]))
      end
      images_K[i] = s
    end
    autos_extended[w] = hom(K, K, Autos_abs[w], images_K)
    #@hassert :NfOrd 1 isconsistent(autos_extended[w])
  end
  res = restriction(K, Cp, autos_extended, incs)
  return res

end

###############################################################################
#
#  Restriction of automorphisms
#
###############################################################################

#This function restricts the automorphisms in autos to the number field generated by the class fields in Cp
# incs are the inclusions of the class fields in K
function restriction(K::NfRelNS{nf_elem}, Cp::Vector{ClassField_pp{S, T}}, autos::Vector{NfRelNSToNfRelNSMor_nf_elem}, incs::Vector{NfRelToNfRelNSMor_nf_elem}) where {S, T}

  C = cyclotomic_extension(base_field(Cp[1]), maximum(degree(x) for x in Cp))
  #First, I compute the images in K of the generators of the class fields
  # and their images under the automorphisms
  gK = gens(K)
  all_pe = Array{Tuple{NfRelNSElem{nf_elem}, Array{NfRelNSElem{nf_elem}, 1}}, 1}(undef, length(Cp))
  for j = 1:length(Cp)
    pe = incs[j](Cp[j].pe)
    tau_pe = Array{typeof(pe), 1}(undef, length(autos))
    for i = 1:length(tau_pe)
      tau_pe[i] = autos[i](pe)
    end
    all_pe[j] = (pe, tau_pe)
  end
  #AA is the target field
  AA, gAA = number_field([c.A.pol for c = Cp], cached = false, check = false)
  #And now, linear algebra to compute the restriction
  #I need the product basis fo all the primitive elements of Cp
  B = Array{NfRelNSElem, 1}(undef, degree(AA))
  B[1] = K(1)
  for i = 2:degree(Cp[1])
    B[i] = all_pe[1][1]*B[i-1]
  end
  ind = degree(Cp[1])
  for jj = 2:length(Cp)
    el = all_pe[jj][1]
    for i = 2:degree(Cp[jj])
      for j = 1:ind
        B[(i-1)* ind + j] = B[j]* el
      end
      el *= all_pe[jj][1]
    end
    ind *= degree(Cp[jj])
  end

  #Now, I construct the corresponding sparse matrix
  M = sparse_matrix(base_field(K))
  for i = 1:length(B)
    push!(M, SRow(B[i]))
  end

  b_AA = basis(AA)
  Mk = _expand(M, pseudo_inv(C.mp[1]))
  all_im = Array{Array{NfRelNSElem{nf_elem}, 1}, 1}(undef, length(autos))
  for i = 1:length(autos)
    all_imCp = Array{NfRelNSElem{nf_elem}, 1}(undef, length(Cp))
    for jj=1:length(Cp)
      N = SRow(all_pe[jj][2][i])
      Nk = _expand(N, pseudo_inv(C.mp[1]))
      n = solve(Mk, Nk)
      im = sum(v*b_AA[l] for (l, v) = n)
      all_imCp[jj] = im
    end
    all_im[i] = all_imCp
  end
  return all_im

end

function _find_embedding(KK::KummerExt, el::FacElem{nf_elem, AnticNumberField}, ord_el::Int, frob_gens::Array{NfOrdIdl, 1})
  #Compute the action of the Frobenius on the generators and on tau(a)
  imgs_rhs = Array{Int, 1}(undef, length(frob_gens))
  imgs_lhs = Array{GrpAbFinGenElem, 1}(undef, length(frob_gens))
  i = 0
  for P in frob_gens
    i += 1
    imgs_lhs[i] = canonical_frobenius(P, KK)
    imgs_rhs[i] = canonical_frobenius(P, KK, el)*divexact(KK.n, ord_el)
  end
  # Now, I have to solve the system.
  # Careful! I have to multiply the components with their difference with the exponent :(
  G = KK.AutG
  #In H, I need a copy for every relation I have
  H = abelian_group(fmpz[KK.n for i = 1:length(imgs_rhs)])
  imgs = Array{GrpAbFinGenElem, 1}(undef, ngens(G))
  for i = 1:length(KK.gen)
    m = Array{Int, 1}(undef, length(imgs_lhs))
    d = divexact(KK.n, order(G[i]))
    for j = 1:length(imgs_lhs)
      m[j] = imgs_lhs[j][i]*d
    end
    imgs[i] = H(m)
  end
  mp = hom(gens(G), imgs, check = true)
  b = H(imgs_rhs)
  fl, coord = haspreimage(mp, b)
  if !fl
    return false, coord, el
  end

  #Now, I need the element of the base field
  prod_gens = KK.gen[1]^(div(-coord[1]*ord_el, Int(order(KK.AutG[1]))))
  for i = 2:length(KK.gen)
    mul!(prod_gens, prod_gens, KK.gen[i]^(div(-coord[i]*ord_el, Int(order(KK.AutG[i])))))
  end
  mul!(prod_gens, prod_gens, el)
  fl2, rt = ispower(prod_gens, ord_el, with_roots_unity = true)
  return fl2, coord, rt
end

################################################################################
#
#  Embeddings
#
################################################################################

function extend_hom(A::ClassField, B::ClassField, tau::T) where T <: Map
  # tau: k1       -> k2
  #global last_extend = (A, tau)
  k1 = domain(tau)
  k2 = codomain(tau)
  @assert k1 == base_field(A)
  @assert k2 == base_field(B)
  @assert degree(B) % degree(A) == 0 #actually, this should hold for the exponent
  lp = factor(fmpz(degree(B)))
  all_h = [A.A() for x in A.cyc]
  for (p, v) = lp.fac
    Cp = [Ap for Ap = A.cyc if degree(Ap) % Int(p) == 0]
    Dp = [Bp for Bp = B.cyc if degree(Bp) % Int(p) == 0]
    h = [extend_hom(X, Cp, tau) for x = Dp]
  end
end

function extend_hom(C::ClassField_pp, D::Array{ClassField_pp, 1}, tau)
    #if it works, then Cp -> Dp should also work
    k2 = codomain(tau)
    k1 = domain(tau)
    i = 1
    om = 0
    im = 0
    while i <= length(D)
      if degree(D[i]) > om
        om = degree(D[i])
        im = i
      end
      i += 1
    end
    # now Dp[im] is of maximal exponent - hence, it should have the maximal
    # big Kummer extension. By construction (above), the set of s-units
    # SHOULD guarantee this....
    # om defintely has the maximal base field, ie. the most roots of 1.
    # Now I want (the images) for all generators in terms of this large Kummer field.
    #
    Dy = cyclotomic_extension(k2, Int(om))
    Cy = cyclotomic_extension(k1, C.degree)
    g = Cy.Kr.pol
    tau_g = k2["x"][1]([tau(coeff(g, i)) for i=0:degree(g)])
    println("g: $g")
    println("tau(g): $tau_g")
    i = 1
    z = gen(Dy.Kr)
    while gcd(i, om) != 1 || !iszero(tau_g(z))
      i *= 1
      z *= gen(Dy.Kr)
    end
    z_i = i

    z_i_inv = invmod(z_i, om)

    Tau = NfRelToNfRelMor(Cy.Kr, Dy.Kr, tau, z)
    @show tau_Ka = hom(Cy.Ka, Dy.Ka, Dy.mp[1]\(Tau(Cy.mp[1](gen(Cy.Ka)))), check = false)

    lp = collect(keys(D[im].bigK.frob_cache))
    pp = maximum(minimum(x) for x = lp)
    S = Base.Iterators.flatten((lp, PrimeIdealsSet(order(lp[1]), pp, fmpz(-1), indexdivisors=false, ramified=false, degreebound = 1)))

    @assert Dy.Ka == base_field(D[im].K)

    all_s = []
    all_tau_s = []
    all_emb = []
    for c in D
#      println("om: $om -> ", degree(c), " vs ", c.o)
      Cs = cyclotomic_extension(k2, Int(degree(c)))
      Emb = hom(Cs.Kr, Dy.Kr, gen(Dy.Kr)^div(om, degree(c)), check = false)
      emb = Cs.mp[1] * Emb * pseudo_inv(Dy.mp[1])
      a = FacElem(Dict(emb(k) => v for (k,v) = c.a.fac))
      push!(all_emb, (a, emb, divexact(om, c.o)))
    end
    b = FacElem(Dict(tau_Ka(k) => v for (k,v) = C.a.fac))

    G = abelian_group([om for i=1:length(D[im].bigK.gen)])
    Q, mQ = quo(G, elem_type(G)[])
    U = abelian_group([om for i = D])
    s_gen = elem_type(U)[]
    tau_b = fmpz[]

    for p = S
      local f
      local fa
      local tfa
      try
        f = canonical_frobenius(p, D[im].bigK).coeff
        fa = [canonical_frobenius(p, D[im].bigK, a[1]) for a = all_emb]
        tfa = canonical_frobenius(p, D[im].bigK, b)
      catch e
        if typeof(e) != BadPrime
          rethrow(e)
        end
        continue
      end
      el = mQ(G(f))
      if iszero(el)
        continue
      end
      Q, mmQ = quo(Q, [el])
      mQ = mQ*mmQ
      push!(s_gen, U(fa))
      push!(tau_b, (tfa))
      if order(Q) == 1
        break
      end
    end

    T_grp = abelian_group([om for i= s_gen])
    @show t_gen = [T_grp([x[i] for x = s_gen]) for i=1:length(D)]
    @show t_tau_g = T_grp(tau_b)
    @show t_corr = [gcd(content(x.coeff), om) for x = t_gen]
    @show t_corr_b = gcd(gcd(tau_b), om)
    @assert t_corr_b == 1
    #if any entry in t_corr is != 1, then the degree of the kummer
    #extension has to be lower:
    #support C2 x C8, the generator for the C2 is in the Cylo(8),
    #thus over the larger base field, the extension is trivial

    q, mq = quo(T_grp, divexact(C.o, Int(t_corr_b)))
    @assert domain(mq) == T_grp
    _, ms = sub(q, [mq(x) for x = t_gen])
    @show fl, lf = haspreimage(ms, mq(t_tau_g))
    @assert fl
    mu = prod(all_emb[j][1]^lf[j] for j=1:length(D)) * inv(b)
    fl, rt = ispower(mu, divexact(C.o, Int(t_corr_b)))
    @assert fl
    all_b = (evaluate(rt), lf)

    Ka = Dy.Ka
    KaT, X = PolynomialRing(Ka, "T", cached = false)
    KK, gKK = number_field([X^Int(divexact(D[j].o, t_corr[j])) - root(evaluate(all_emb[j][1]), Int(t_corr[j])) for j=1:length(D)], check = false)
    s = gKK[1]
    s = s^Int(divexact(D[1].o, C.o)*all_b[2][1])
    for j in 2:length(D)
      s = s * gKK[j]^Int(divexact(D[j].o, C.o)*all_b[2][j])
    end
    h = NfRelToNfRelNSMor(C.K, KK, tau_Ka, inv(all_b[1]) * s)

    # now "all" that remains is to restrict h to the subfield, using lin. alg..

    all_pe = []
    for jj=1:length(D)
      emb = hom(D[jj].K, KK, tau_Ka, gens(KK)[jj], check = false)
      pe = emb(D[jj].pe)
      push!(all_pe, pe)
    end

    B = [KK(1), all_pe[1]]
    d = degree(D[1])
    while length(B) < degree(D[1])
      push!(B, B[end]*all_pe[1])
    end

    for jj=2:length(D)
      d *= degree(D[jj])
      _D = copy(B)
      while length(B) < d
        _D = [x*all_pe[jj] for x = _D]
        append!(B, _D)
      end
    end
    M = sparse_matrix(Ka)
    for i=1:d
      push!(M, SRow(B[i]))
    end
    AA, gAA = number_field([c.A.pol for c = D], check = false)
    @assert d == degree(AA)
    @assert d == length(B)
    b_AA = basis(AA)
    Mk = _expand(M, pseudo_inv(Dy.mp[1]))
    #@hassert :ClassField 2 nullspace(Mk')[1] == 0
    N = SRow(h(C.pe))
    Nk = _expand(N, pseudo_inv(Dy.mp[1]))
    n = solve(Mk, Nk)
    all_im = sum(v*b_AA[l] for (l, v) = n)

      return all_im

      #=

    im = NfRelNSElem{nf_elem}[]
    i = 1
    j = 1
    while j<=length(A.cyc)
      if i<= length(Cp) && degree(A.cyc[j]) == degree(Cp[i])
        push!(im, gens(A.A)[j])
        i += 1
        j += 1
      else
        j += 1
      end
    end
    emb = NfRelNSToNfRelNSMor(KK, A.A, im)
    i = 1
    j = 1
    while j<=length(A.cyc)
      if i<= length(Cp) && degree(A.cyc[j]) == degree(Cp[i])
        all_h[j] = emb(all_im[i])
        i += 1
        j += 1
      else
        j += 1
      end
    end
  end
  return NfRelNSToNfRelNSMor(A.A, A.A, tau, all_h)
  =#
end

#M is over K, mp: K -> K/k, expand M into a matrix over k
function _expand(M::Generic.Mat{nf_elem}, mp::Map)
  Kr = domain(mp)
  Ka = codomain(mp)
  k = base_field(Kr)
  d = degree(Kr)
  N = zero_matrix(k, nrows(M), ncols(M) * d)
  for i=1:nrows(M)
    for j = 1:ncols(M)
      a = mp\M[i,j]
      for l=0:d-1
        N[i, (j-1)*d+l+1] = coeff(a, l)
      end
    end
  end
  return N
end

function _expand(M::SRow{nf_elem}, mp::Map)
  Kr = domain(mp)
  k = base_field(Kr)
  d = degree(Kr)
  sr = SRow(k)
  for (j, v) = M
    a = mp\v
    for l=0:d-1
      c = coeff(a, l)
      if !iszero(c)
        push!(sr.pos, (j-1)*d+1+l)
        push!(sr.values, c)
      end
    end
  end
  return sr
end

function _expand(M::SMat{nf_elem}, mp::Map)
  Kr = domain(mp)
  k = base_field(Kr)
  N = sparse_matrix(k)
  for i=1:nrows(M)
    sr = _expand(M[i], mp)
    push!(N, sr)
  end
  return N
end

################################################################################
#
#  Extend auto - intersection with cyclotomic extension not trivial
#
################################################################################

function extend_autos_hard_case(A::ClassField, autos::Vector{NfToNfMor}, p::fmpz, deg_intersection::Int)
  error("Not yet implemented")
  Cp = [x1 for x1 in A.cyc if degree(x1) % Int(p) == 0]
  k = base_field(A)
  C = cyclotomic_extension(k, d)
  KC = absolute_simple_field(C)
  # C is the base field of the kummer extension generated
  # by all the cyclic components.
  # I extend the automorphisms to C
  Autos_abs = Array{NfToNfMor, 1}(undef, length(autos))
  for i = 1:length(autos)
    aut = extend_to_cyclotomic(C, autos[i])
    Autos_abs[i] = hom(KC, KC, C.mp[1]\(aut(C.mp[1](gen(KC)))), check = false)
  end

  #I compute the embeddings of the small cyclotomic extensions into the others
  abs_emb = Array{NfToNfMor, 1}(undef, length(Cp))
  for i = 1:length(Cp)
    dCp = degree(Cp[i])
    if dCp == d
      abs_emb[i] = id_hom(KC)
    else
      Cs = cyclotomic_extension(k, dCp)
      emb = hom(Cs.Kr, C.Kr, gen(C.Kr)^div(d, dCp), check = false)
      img = C.mp[1]\(emb(Cs.mp[1](gen(Cs.Ka))))
      abs_emb[i] = hom(Cs.Ka, KC, img, check = false)
    end
  end

  #Now, I can compute the corresponding Kummer extension over the big cyclotomic field.
  exps = Array{Int, 1}(undef, length(Cp))
  gens = Array{FacElem{nf_elem, AnticNumberField}, 1}(undef, length(Cp))
  for i = 1:length(Cp)
    if degree(Cp[i]) == d
      gens[i] = Cp[i].a
      exps[i] = Cp[i].o
    else
      D = Dict{nf_elem, fmpz}()
      for (ke,v) in Cp[i].a
        D[abs_emb[i](ke)] = v
      end
      a = FacElem(D)
      exps[i] = Cp[i].o
      gens[i] = a
    end
  end
  KK = kummer_extension(exps, gens)
  #Now, KK is not a real kummer extension. There are relations. We need to find them.
  #We know that the right degree of the translation of Cp to C is degree(Cp)/idx
  #So we compute Frobenius automorphisms until we generate a subgroup with the right subgroup.

  act_on_gens = Array{Array{FacElem{nf_elem, AnticNumberField}, 1}, 1}(undef, length(KK.gen))
  for i = 1:length(KK.gen)
    act_on_gen_i = Array{FacElem{nf_elem, AnticNumberField}, 1}(undef, length(autos))
    for j = 1:length(autos)
      D1 = Dict{nf_elem, fmpz}()
      for (ke, v) in KK.gen[i]
        D1[Autos_abs[j](ke)] = v
      end
      act_on_gen_i[j] = FacElem(D1)
    end
    act_on_gens[i] = act_on_gen_i
  end
  frob_gens = find_gens(KK, act_on_gens, minimum(defining_modulus(A)[1]), deg_intersection)
  #Now, I can create a proper kummer extension.
  s, ms = sub(KK.AutG, GrpAbFinGenElem[canonical_frobenius(x, KK) for x in frob_gens])
  @assert order(s) == divexact(order(KK.AutG), deg_intersection)
  S, mS = snf(s)
  gens_real_KK = Vector{FacElem{nf_elem, AnticNumberField}}(undef, ngens(S))
  for i = 1:ngens(S)
    imgSi = mS(S[i])
    gens_real_KK[i] = prod(KK.gens[j]^imgSi[j] for j = 1:ngens(KK))
  end
  KK_real = kummer_extension(gens_real_KK, S.snf)
  new_act_on_gens = Array{Array{FacElem{nf_elem, AnticNumberField}, 1}, 1}(undef, length(KK.gen))
  for i = 1:length(KK.gen)
    act_on_gen_i = Array{FacElem{nf_elem, AnticNumberField}, 1}(undef, length(autos))
    for j = 1:length(autos)
      D1 = Dict{nf_elem, fmpz}()
      for (ke, v) in KK.gen[i]
        D1[Autos_abs[j](ke)] = v
      end
      act_on_gen_i[j] = FacElem(D1)
    end
    new_act_on_gens[i] = act_on_gen_i
  end
  @show find_frob = find_gens(KK_real, new_act_on_gens, minimum(defining_modulus(A)[1]))
  error("stop")
  AA, gAA = number_field(KK_real)
  #I now need the embedding of Cp into KK_real

end