text:
using Distributions
using Random
using LinearAlgebra
using CSV

mutable struct myPoint
    x::Float64
    y::Float64
    close::Int16
    far::Int16
    row_num::Int16
    col_num::Int16
    key::Int64
end

mutable struct myGrid
    points::Array{myPoint}
end

mutable struct Centroid
    x::Float64
    y::Float64
    close_num::Int16
end

function assignPoint(i::Int, x::Float64, y::Float64, close::Int, far::Int, grid_size::Float64)
    row = floor(Int8, x/grid_size)
    col = floor(Int8, y/grid_size)
    max_num =  ceil(107 / grid_size)
    key = Int(row*max_num + col)
    myPoint(x, y, close, far, row, col, key)
end

function myGrid()
    myGrid([])
end

groupby(f, list::Array ) = begin
    foldl(list; init = Dict()) do dict, v
        push!(get!(dict, f(v), []), v)
        dict
    end
end

function get_point_key(p::myPoint)
    p.key
end

function get_centroid(points)
    x_sum = 0
    y_sum = 0
    count = 0
    close_num = 0
    for p in points
        x_sum += p.x *p.close
        y_sum += p.y *p.close
        close_num += p.close
        count += 1
    end
    Centroid(x_sum/close_num, y_sum/close_num, floor(close_num/count))
end

function sample_grid(one_centroid::Centroid, grid_size::Float64)
    r_tmp = rand(Uniform(0, grid_size/sqrt(pi)), one_centroid.close_num)
    theta_tmp = rand(Uniform(0, 2pi), one_centroid.close_num)
    X_tmp = one_centroid.x .+ cos.(theta_tmp) .* r_tmp
    Y_tmp = one_centroid.y .+ sin.(theta_tmp) .* r_tmp
    [X_tmp Y_tmp]
end

function generate_samples(project_path, filename)
    grid_size = 1.25
    arr = CSV.File(project_path * "/" * filename)
    grid = myGrid()
    for (i, row) in enumerate(arr)
        if row.Close == 0
            continue
        end
        point = assignPoint(i, row.X, row.Y, row.Close, row.Far, grid_size)
        push!(grid.points, point)
    end

    dict_1 = groupby(get_point_key, grid.points)
    centroids = map(get_centroid, collect(values(dict_1)))
    array_list = map(sample_grid, centroids, ones(length(centroids)) .+ grid_size)
    mapconcat2(v) = vcat((v)...)
    samples = mapconcat2(array_list)
    return samples;
end

function mlt_pdf(data, mu, cov)
    try
        distribution = MvNormal(mu, Matrix(Hermitian(cov)));
        res = pdf(distribution, data);
        if isnan(res)
            return 1e6;
        end
        return res;
    catch
        return 1e-6;
    end
end

function k_means(k, tol, max_iter, data)
    centroids = data[rand(1:size(data)[1], k), :];
    for i = 1:max_iter
        classifications = Dict();
        for j = 1:k
            classifications[j] = [0 0];
        end
        for j = 1:length(data[:,1])
            distances = [];
            for jj = 1:k
                push!(distances, norm(data[j,:] .- centroids[jj, :]));
            end
            classifications[mapslices(argmin, distances, dims=1)[1]] = [classifications[mapslices(argmin, distances, dims=1)[1]]; reshape(data[j,:],(1,2))];
        end

        for j = 1:k
            classifications[j] = classifications[j][2:length(classifications[j][:,1]), :]
        end
        prev_centroids = centroids;

        for j = 1:k
            if (size(classifications[j])[1] != 0)
                centroids[j,:] = mean(classifications[j], dims = 1);
            end
        end

        optimized = true;
        for j in 1:length(centroids[:,1])
            original_centroid = prev_centroids[j];
            current_centroid = centroids[j];
            if abs(sum((current_centroid .- original_centroid) ./ original_centroid .*100.0)) > tol
                    optimized = false;
            end
        end
        if optimized
            break
        end
    end
    return centroids
end

function gmm_em(max_iterations, num_clusters, data)
    pi_ = ones(num_clusters)/num_clusters;
    mu = sort(k_means(num_clusters, 0.001, 200, data), dims = 1);
    for i = 1:2
        mu = mu .+ sort(k_means(num_clusters, 0.001, 200, data), dims = 1);
    end
    mu = mu ./ 3;
    cov = []
    for i = 1:num_clusters
        push!(cov, [[5. 0.]; [0. 5.]]);
    end
    ric = zeros((length(data[:,1]), num_clusters));
    reg_cov = 1e-6 .* [[5. 0.];[0. 5.]];
    log_likelihoods = [];
    for i = 1:max_iterations
        ric = zeros((length(data[:,1]), num_clusters));

        for j = 1:num_clusters
            covc = cov[j] .+ reg_cov;
            mn = [0];
            for num = 1:length(data[:,1])
                mn = [mn; [mlt_pdf(data[num, :], mu[j,:], covc)]];
            end

            mn = mn[2:length(mn)];
            ric[:, j] = pi_[j] .* mn;
        end

        for j = 1:length(data[:,1])
            ric[j, :] = ric[j, :] ./ sum(ric[j, :]);
        end

        mc = sum(ric, dims=1);
        pi_ = mc ./ sum(mc);
        mu = transpose((transpose(data) * ric) ./ [mc; mc]);

        cov = []
        for j = 1:num_clusters
            var1 = 1/mc[j] * dot(ric[:, j], (data[:,1] .- mu[j, 1]).^2);
            var2 = 1/mc[j] * dot(ric[:, j], (data[:,2] .- mu[j, 2]).^2);
            push!(cov, [[var1 0.]; [0. var2]]);
        end

        likelihood_sum = 0;
        for j = 1:num_clusters
            mn = [0];
            for num = 1:length(data[:,1])
                mn = [mn; [mlt_pdf(data[num,:], mu[j,:], cov[j] + reg_cov)]];
            end
            mn = mn[2:length(mn)];
            likelihood_sum = likelihood_sum .+ ric[:,j].*mn;
        end
        push!(log_likelihoods, sum(log.(likelihood_sum)));
        if i >= 20 && abs((log_likelihoods[i] - log_likelihoods[i-1])/log_likelihoods[i]) < 0.01
            break;
        end
    end

    likelihood = log_likelihoods[length(log_likelihoods)];
    bic = - 2*likelihood + log(size(data)[1])*(15*num_clusters);
    return pi_, mu, cov, ric, log_likelihoods, bic;
end

function main()
    # path = "/Users/clara/Desktop"
    # f_name = "data.csv"
    path = "D:/JI/2020 summer/VE414 Bayesian Analysis/Project/zzp";
    f_name = "data_proj_414.csv";
    data = generate_samples(path, f_name);
    println("----Samples are generated successfully----");
    max_iterations = 50;

    println("-----------Start training model-----------");
    bics = [0];
    for k=1:30
        pi_, mu, cov, ric, log_likelihoods, bic = gmm_em(max_iterations , k, data);
        pi_, mu, cov, ric, log_likelihoods, bic_ = gmm_em(max_iterations , k, data);
        bics = [bics;minimum([bic, bic_])];
        println("--Clusters: $(k)  BIC: $(mean([bic, bic_]))--");
    end
    bics = bics[2:length(bics[:,1]),:];
    println("-----------Stop training model------------");

    original_data = CSV.read(path*"/"*f_name);
    positions = [original_data["X"] original_data["Y"]];
    area_grid = zeros((107, 107));
    for i = 1:size(positions)[1]
        area_grid[Int32(positions[i, 1] - positions[i, 1] % 1) + 1, Int32(positions[i, 2] - positions[i, 2] % 1) + 1] = 1
    end
    area = sum(area_grid);
    lamb =  (mapslices(argmin,bics[:,1],dims=1)[1] * 107 * 107 / area);

    total = 0;
    mode_ = pdf(Poisson(lamb), 0);

    for i = 1:2*round(lamb)
        if pdf(Poisson(lamb), i) > mode_
            mode_ = pdf(Poisson(lamb), i);
            total = i;
        end
    end
    println("-------Total number of Jiulings $(total)------");
    return total;
end


text:
using LiquidsStructure
using Test

const AHS = AttractiveHardSpheres
const k∞ = inv(eps()) # large wavevector

@testset "Hard disks, Rosenfeld FMT" begin
    η = 0.1
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.6604201250395436
    @test f(π/2) ≈ 0.7269796280941714
    @test f(1π ) ≈ 0.9054985794257596
    @test f(2π ) ≈ 1.0331423908135913
    @test f(k∞ ) ≈ 1.0

    η = 0.2
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.42121335772860696
    @test f(π/2) ≈ 0.4957408483427397
    @test f(1π ) ≈ 0.7616284799236219
    @test f(2π ) ≈ 1.0838625683290253
    @test f(k∞ ) ≈ 1.0

    η = 0.3
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.256818185902224
    @test f(π/2) ≈ 0.3150257470860998
    @test f(1π ) ≈ 0.5746880320600511
    @test f(2π ) ≈ 1.1663158306646833
    @test f(k∞ ) ≈ 1.0

    η = 0.4
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.14747909074224583
    @test f(π/2) ≈ 0.18467032538432618
    @test f(1π ) ≈ 0.37610672082612917
    @test f(2π ) ≈ 1.312642950246751
    @test f(k∞ ) ≈ 1.0

    η = 0.5
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.0779666302822392
    @test f(π/2) ≈ 0.09810277945694483
    @test f(1π ) ≈ 0.2083203632199568
    @test f(2π ) ≈ 1.608912503571731
    @test f(k∞ ) ≈ 1.0

    η = 0.6
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.036533750792325716
    @test f(π/2) ≈ 0.04568348215003821
    @test f(1π ) ≈ 0.09561925118936167
    @test f(2π ) ≈ 2.325299703948229
    @test f(k∞ ) ≈ 1.0

    η = 0.7
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.01413529677840881
    @test f(π/2) ≈ 0.01744107861016397
    @test f(1π ) ≈ 0.03474517028685856
    @test f(2π ) ≈ 3.298291540535111
    @test f(k∞ ) ≈ 1.0

    η = 0.8
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.003851154422049552
    @test f(π/2) ≈ 0.004670418158317173
    @test f(1π ) ≈ 0.008718525753638387
    @test f(2π ) ≈ 0.572989315611843
    @test f(k∞ ) ≈ 1.0

    η = 0.9
    f = StructureFactor(HardDisks(η), RosenfeldFMT)
    @test f(0.0) ≈ 0.0004441340638213542
    @test f(π/2) ≈ 0.0005286284359062185
    @test f(1π ) ≈ 0.0009226519112744395
    @test f(2π ) ≈ 0.025405875690152686
    @test f(k∞ ) ≈ 1.0
end

@testset "Hard spheres, Percus–Yevick" begin
    η = 0.1
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.455625
    @test f(π/2) ≈ 0.5162698774690883
    @test f(1π ) ≈ 0.7166968096587476
    @test f(2π ) ≈ 1.0831625291823719
    @test f(k∞ ) ≈ 1.0

    η = 0.2
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.20897959183673484
    @test f(π/2) ≈ 0.2492889697353935
    @test f(1π ) ≈ 0.42259178223764265
    @test f(2π ) ≈ 1.2365374806798426
    @test f(k∞ ) ≈ 1.0

    η = 0.3
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.0937890624999999
    @test f(π/2) ≈ 0.113580452287522
    @test f(1π ) ≈ 0.2083407060258776
    @test f(2π ) ≈ 1.513704692974072
    @test f(k∞ ) ≈ 1.0

    η = 0.4
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.04
    @test f(π/2) ≈ 0.04834624242434464
    @test f(1π ) ≈ 0.0889840427994623
    @test f(2π ) ≈ 1.7888899102094844
    @test f(k∞ ) ≈ 1.0

    η = 0.5
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.015625000000000007
    @test f(π/2) ≈ 0.01870914640113364
    @test f(1π ) ≈ 0.0333665751636566
    @test f(2π ) ≈ 1.0433145711209884
    @test f(k∞ ) ≈ 1.0

    η = 0.6
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.005289256198347106
    @test f(π/2) ≈ 0.006257471729168288
    @test f(1π ) ≈ 0.010689789920403656
    @test f(2π ) ≈ 0.23121862372838764
    @test f(k∞ ) ≈ 1.0

    η = 0.7
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.0014062500000000021
    @test f(π/2) ≈ 0.00164315820030793
    @test f(1π ) ≈ 0.002687286985730639
    @test f(2π ) ≈ 0.03729828280050181
    @test f(k∞ ) ≈ 1.0

    η = 0.8
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 0.00023668639053254386
    @test f(π/2) ≈ 0.0002733505741836147
    @test f(1π ) ≈ 0.00042959576011198774
    @test f(2π ) ≈ 0.004276959009112015
    @test f(k∞ ) ≈ 1.0

    η = 0.9
    f = StructureFactor(HardSpheres(η), PercusYevick)
    @test f(0.0) ≈ 1.2755102040816327e-5
    @test f(π/2) ≈ 1.4577626809289744e-5
    @test f(1π ) ≈ 2.2128300716842136e-5
    @test f(2π ) ≈ 0.0001727400947607701
    @test f(k∞ ) ≈ 1.0
end

@testset "Hard spheres, Verlet–Weis" begin
    η = 0.1
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.4578452915332574
    @test f(π/2) ≈ 0.5182637197480197
    @test f(1π ) ≈ 0.7174855200551944
    @test f(2π ) ≈ 1.0829101226236082
    @test f(k∞ ) ≈ 1.0

    η = 0.2
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.2131237008975026
    @test f(π/2) ≈ 0.25367291404685227
    @test f(1π ) ≈ 0.4267618525871056
    @test f(2π ) ≈ 1.233011655529535
    @test f(k∞ ) ≈ 1.0

    η = 0.3
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.09821655794190104
    @test f(π/2) ≈ 0.11862134915114293
    @test f(1π ) ≈ 0.21541648542267264
    @test f(2π ) ≈ 1.4916754221537918
    @test f(k∞ ) ≈ 1.0

    η = 0.4
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.04369978853680091
    @test f(π/2) ≈ 0.05268107258287156
    @test f(1π ) ≈ 0.09605449462293068
    @test f(2π ) ≈ 1.7231617060108197
    @test f(k∞ ) ≈ 1.0

    η = 0.5
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.018237055564413282
    @test f(π/2) ≈ 0.021789977028257794
    @test f(1π ) ≈ 0.03857987865992043
    @test f(2π ) ≈ 1.0818137287598686
    @test f(k∞ ) ≈ 1.0

    η = 0.6
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.006861393142599895
    @test f(π/2) ≈ 0.008104719915885968
    @test f(1π ) ≈ 0.01377843599859209
    @test f(2π ) ≈ 0.2882778877829879
    @test f(k∞ ) ≈ 1.0

    η = 0.7
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.0021846162335262737
    @test f(π/2) ≈ 0.0025501353769391737
    @test f(1π ) ≈ 0.004158905666715577
    @test f(2π ) ≈ 0.05749168504321367
    @test f(k∞ ) ≈ 1.0

    η = 0.8
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 0.0005224489795918363
    @test f(π/2) ≈ 0.0006030737594212403
    @test f(1π ) ≈ 0.0009465346587026628
    @test f(2π ) ≈ 0.00944076344580056
    @test f(k∞ ) ≈ 1.0

    η = 0.9
    f = StructureFactor(HardSpheres(η), VerletWeis)
    @test f(0.0) ≈ 7.067452999457804e-5
    @test f(π/2) ≈ 8.075895958572265e-5
    @test f(1π ) ≈ 0.00012253818964376312
    @test f(2π ) ≈ 0.0009578116657578688
    @test f(k∞ ) ≈ 1.0
end

@testset "Yukawa Attractive hard spheres, Sharma–Sharma" begin
    Z  = 2.0
    λ  = 1.5
    U  = Yukawa{3}(Z)

    η  = 0.1
    T′ = 1.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 2.603191157463335
    @test f(π/2) ≈ 0.6863273179451208
    @test f(1π ) ≈ 0.6382449330348262
    @test f(2π ) ≈ 1.1517581024741745

    η  = 0.2
    T′ = 1.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.9156580819599524
    @test f(π/2) ≈ 0.33365513773093114
    @test f(1π ) ≈ 0.37184289969239975
    @test f(2π ) ≈ 1.4273016612357994
    @test f(k∞ ) ≈ 1.0

    η  = 0.3
    T′ = 0.75
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.3353923486908502
    @test f(π/2) ≈ 0.15289981479991796
    @test f(1π ) ≈ 0.1874649216578271
    @test f(2π ) ≈ 2.2242590379730283
    @test f(k∞ ) ≈ 1.0

    η  = 0.4
    T′ = 0.75
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.07528189724199456
    @test f(π/2) ≈ 0.06074518898941603
    @test f(1π ) ≈ 0.08823291923872414
    @test f(2π ) ≈ 3.497363781631624
    @test f(k∞ ) ≈ 1.0

    η  = 0.5
    T′ = 1.5
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.0204777822845714
    @test f(π/2) ≈ 0.022564349488829176
    @test f(1π ) ≈ 0.037740052742288845
    @test f(2π ) ≈ 1.3506683217180429
    @test f(k∞ ) ≈ 1.0

    η  = 0.6
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.007125400066847413
    @test f(π/2) ≈ 0.008198910387096473
    @test f(1π ) ≈ 0.013680583265695689
    @test f(2π ) ≈ 0.3027297985049469
    @test f(k∞ ) ≈ 1.0

    η  = 0.7
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.002215102876656909
    @test f(π/2) ≈ 0.002560935319716773
    @test f(1π ) ≈ 0.004148456520065383
    @test f(2π ) ≈ 0.05813743812055229
    @test f(k∞ ) ≈ 1.0

    η  = 0.8
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.0005244216612367169
    @test f(π/2) ≈ 0.0006037619166927905
    @test f(1π ) ≈ 0.0009459149383432805
    @test f(2π ) ≈ 0.009460483948091025
    @test f(k∞ ) ≈ 1.0

    η  = 0.9
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 7.071501177136412e-5
    @test f(π/2) ≈ 8.077282906281272e-5
    @test f(1π ) ≈ 0.0001225264983879417
    @test f(2π ) ≈ 0.0009580396021185613
    @test f(k∞ ) ≈ 1.0

end

@testset "Square Well Attractive hard spheres, Sharma–Sharma" begin
    λ  = 1.5
    U  = SquareWell{3}(λ)

    η  = 0.1
    T′ = 1.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 3.5193435613655701
    @test f(π/2) ≈ 0.92661127755109707
    @test f(1π ) ≈ 0.58331796691647120
    @test f(2π ) ≈ 1.2962503853125780

    η  = 0.2
    T′ = 1.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 1.1209371060679321
    @test f(π/2) ≈ 0.44614027958622848
    @test f(1π ) ≈ 0.33507838902406228
    @test f(2π ) ≈ 1.9721579006903056
    @test f(k∞ ) ≈ 1.0

    η  = 0.3
    T′ = 0.75
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 0.38735928482005272
    @test f(π/2) ≈ 0.19885016413408196
    @test f(1π ) ≈ 0.16879151714668816
    @test f(2π ) ≈ 16.010240228035098
    @test f(k∞ ) ≈ 1.0

    η  = 0.4
    T′ = 0.75
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 7.8430925994702708e-2
    @test f(π/2) ≈ 6.9217956774394074e-2
    @test f(1π ) ≈ 8.2504930426001266e-2
    @test f(2π ) ≈ -4.3432773350119831
    @test f(k∞ ) ≈ 1.0

    η  = 0.5
    T′ = 1.5
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 2.0618522825051694e-2
    @test f(π/2) ≈ 2.3224344950697608e-2
    @test f(1π ) ≈ 3.7052366895499976e-2
    @test f(2π ) ≈ 2.3936751366867295
    @test f(k∞ ) ≈ 1.0

    η  = 0.6
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 7.1406640933698879e-3
    @test f(π/2) ≈ 8.2758206645957878e-3
    @test f(1π ) ≈ 1.3598245019352462e-2
    @test f(2π ) ≈ 0.33190285202602476
    @test f(k∞ ) ≈ 1.0

    η  = 0.7
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 2.2168215473822354e-3
    @test f(π/2) ≈ 2.5696376277179267e-3
    @test f(1π ) ≈ 4.1395889487837005e-3
    @test f(2π ) ≈ 5.9305356062453886e-2
    @test f(k∞ ) ≈ 1.0

    η  = 0.8
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 5.2453169154913833e-4
    @test f(π/2) ≈ 6.0431333832571899e-4
    @test f(1π ) ≈ 9.4538720482588478e-4
    @test f(2π ) ≈ 9.4952595505868632e-3
    @test f(k∞ ) ≈ 1.0

    η  = 0.9
    T′ = 2.0
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ 7.0717262118774555e-5
    @test f(π/2) ≈ 8.0783923331053070e-5
    @test f(1π ) ≈ 1.2251653217610584e-4
    @test f(2π ) ≈ 9.5843950557686168e-4
    @test f(k∞ ) ≈ 1.0

end

@testset "SharmaSharma{VerletWeis} with no attraction versus VerletWeis" begin
    η = 0.1

    f₀ = StructureFactor(HardSpheres(η), VerletWeis)

    T′ = 1.0
    Z  = 2.0
    U  = Yukawa{3}(Z, ϵ = 0.0)
    f  = StructureFactor(AHS(η, T′, U), SharmaSharma{VerletWeis})
    @test f(0.0) ≈ f₀(0.0)
    @test f(π/2) ≈ f₀(π/2)
    @test f(1π ) ≈ f₀(1π )
    @test f(2π ) ≈ f₀(2π )
    @test f(k∞ ) ≈ f₀(k∞ )
end

@testset "Dipolar hard spheres, MSA" begin
    η = 0.1
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.455625, 0.6291407983765894, 1.2666736611367042) )
    @test all( f(π/2) .≈ (0.5162698774690883, 0.6850767960631556, 1.1938720996192682) )
    @test all( f(1π ) .≈ (0.7166968096587476, 0.8416086006775434, 1.0666596582612127) )
    @test all( f(2π ) .≈ (1.0831625291823719, 1.0426349627796052, 0.9836365170156721) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.2
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.20897959183673484, 0.4431390440180322, 1.5229321819837567) )
    @test all( f(π/2) .≈ (0.2492889697353935, 0.5035359373018747, 1.3590195670624214) )
    @test all( f(1π ) .≈ (0.42259178223764265, 0.7056963342919332, 1.1112828002460329) )
    @test all( f(2π ) .≈ (1.2365374806798426, 1.0871649764504834, 0.9728720629632218) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.3
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.0937890624999999, 0.33702806816959363, 1.7633203545450626) )
    @test all( f(π/2) .≈ (0.113580452287522, 0.391990692356982, 1.4979049530982467) )
    @test all( f(1π ) .≈ (0.2083407060258776, 0.5976735425934862, 1.142269051998251) )
    @test all( f(2π ) .≈ (1.513704692974072, 1.1316981741146153, 0.9653800070336875) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.4
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.04, 0.2700545515466438, 1.9882472781755554) )
    @test all( f(π/2) .≈ (0.04834624242434464, 0.3184635771562895, 1.6157348484418494) )
    @test all( f(1π ) .≈ (0.0889840427994623, 0.5132862854019397, 1.164806497625769) )
    @test all( f(2π ) .≈ (1.7888899102094844, 1.1754992038280996, 0.959870575710502) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.5
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.015625000000000007, 0.2244848150807989, 2.19952173374085) )
    @test all( f(π/2) .≈ (0.01870914640113364, 0.26704088038841756, 1.717078039617591) )
    @test all( f(1π ) .≈ (0.0333665751636566, 0.44699097699243284, 1.1818864428745708) )
    @test all( f(2π ) .≈ (1.0433145711209884, 1.218279674681617, 0.9556335505378297) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.6
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.005289256198347106, 0.19168706560349014, 2.399024573085883) )
    @test all( f(π/2) .≈ (0.006257471729168288, 0.22934004449014636, 1.805425696890689) )
    @test all( f(1π ) .≈ (0.010689789920403656, 0.3942033193842698, 1.1952745418340207) )
    @test all( f(2π ) .≈ (0.23121862372838764, 1.259913401542428, 0.9522584797097164) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.7
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.0014062500000000021, 0.1670469612108733, 2.5883898638161775) )
    @test all( f(π/2) .≈ (0.00164315820030793, 0.20064493447924217, 1.8833761824596493) )
    @test all( f(1π ) .≈ (0.002687286985730639, 0.35151915456843674, 1.2060595083159515) )
    @test all( f(2π ) .≈ (0.03729828280050181, 1.3003308303131165, 0.9494945625210339) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.8
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (0.00023668639053254386, 0.14790394033902815, 2.7689702226837616) )
    @test all( f(π/2) .≈ (0.0002733505741836147, 0.17813835606095949, 1.952872588160253) )
    @test all( f(1π ) .≈ (0.00042959576011198774, 0.31648036971131643, 1.2149430519305693) )
    @test all( f(2π ) .≈ (0.004276959009112015, 1.3394797925053805, 0.9471805606574033) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )

    η = 0.9
    f = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    @test all( f(0.0) .≈ (1.2755102040816327e-5, 0.13262741072118495, 2.9418730238949813) )
    @test all( f(π/2) .≈ (1.4577626809289744e-5, 0.1600487295395634, 2.015388817420257) )
    @test all( f(1π ) .≈ (2.2128300716842136e-5, 0.28731317162020314, 1.2223958555431893) )
    @test all( f(2π ) .≈ (0.0001727400947607701, 1.3773111025661495, 0.9452081841208045) )
    @test all( f(k∞ ) .≈ (1.0, 1.0, 1.0) )
end

@testset "DipolarHardSpheres, MSA–PercusYevick–VerletWeis" begin
    η = 0.1
    f₁ = StructureFactor(HardSpheres(η), PercusYevick)
    g₁ = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    f₂ = StructureFactor(HardSpheres(η), VerletWeis)
    g₂ = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{VerletWeis})
    k = 0.0
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = π / 2
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = 1π
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = 2π
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]

    η = 0.45
    f₁ = StructureFactor(HardSpheres(η), PercusYevick)
    g₁ = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    f₂ = StructureFactor(HardSpheres(η), VerletWeis)
    g₂ = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{VerletWeis})
    k = 0.0
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = π / 2
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = 1π
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = 2π
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]

    η = 0.9
    f₁ = StructureFactor(HardSpheres(η), PercusYevick)
    g₁ = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{PercusYevick})
    f₂ = StructureFactor(HardSpheres(η), VerletWeis)
    g₂ = StructureFactor(DipolarHardSpheres(η, 1.0), MSA{VerletWeis})
    k = 0.0
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = π / 2
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = 1π
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
    k = 2π
    v₁, v₂ = g₁(k), g₂(k)
    @test v₁[1] == f₁(k)
    @test v₂[1] == f₂(k)
    @test v₁[2] == v₂[2]
    @test v₁[3] == v₂[3]
end

@testset "Utils" begin
    S = StructureFactor(HardDisks(0.1), RosenfeldFMT)
    @test dimensionality(S) == dimensionality(S.liquid) == 2
    @test repr(S.liquid) == "HardDisks(0.1, 1.0, HardCore{2}(1))"
    @test repr(S.scheme) == "RosenfeldFMT"
    @test repr(S) == "StructureFactor(HardDisks(0.1, 1.0, HardCore{2}(1)), RosenfeldFMT)"

    S = StructureFactor(HardSpheres(0.2), PercusYevick)
    @test dimensionality(S) == dimensionality(S.liquid) == 3
    @test repr(S.liquid) == "HardSpheres(0.2, 1.0, HardCore{3}(1))"
    @test repr(S.scheme) == "PercusYevick"
    @test repr(S) == "StructureFactor(HardSpheres(0.2, 1.0, HardCore{3}(1)), PercusYevick)"

    S = StructureFactor(HardSpheres(0.3), VerletWeis)
    @test repr(S.liquid) == "HardSpheres(0.3, 1.0, HardCore{3}(1))"
    @test repr(S.scheme) == "VerletWeis"
    @test repr(S) == "StructureFactor(HardSpheres(0.3, 1.0, HardCore{3}(1)), VerletWeis)"

    S = StructureFactor(DipolarHardSpheres(0.4, 1.0), MSA{PercusYevick})
    @test dimensionality(S) == dimensionality(S.liquid) == 3
    potential = "CompositePotential(HardCore{3}(1), DipoleDipole(1))"
    @test repr(S.liquid) == "DipolarHardSpheres(0.4, 1.0, $potential)"
    @test repr(S.scheme) == "MSA{PercusYevick}"
    @test repr(S) == "StructureFactor(DipolarHardSpheres(0.4, 1.0, $potential), MSA{PercusYevick})"
end

ext: jl


text:
using AlphaStableDistributions
using Test, Random, Distributions

@testset "AlphaStableDistributions.jl" begin

    stabletypes = [AlphaStable,SymmetricAlphaStable]
    αs = [0.6:0.1:2,1:0.1:2]
    for (i, stabletype) in enumerate(stabletypes)
        for α in αs[i]
            d1 = AlphaStable(α=α)
            s = rand(d1, 100000)

            d2 = fit(stabletype, s)

            @test d1.α ≈ d2.α rtol=0.1
            stabletype != SymmetricAlphaStable && @test d1.β ≈ d2.β atol=0.2
            @test d1.scale ≈ d2.scale rtol=0.1
            @test d1.location ≈ d2.location atol=0.1
        end

        xnormal = rand(Normal(3.0, 4.0), 96000)
        d = fit(stabletype, xnormal)
        @test d.α ≈ 2 rtol=0.2
        stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
        @test d.scale ≈ 4/√2 rtol=0.2
        @test d.location ≈ 3 rtol=0.1

        xcauchy = rand(Cauchy(3.0, 4.0), 96000)
        d = fit(stabletype, xcauchy)
        @test d.α ≈ 1 rtol=0.2
        stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
        @test d.scale ≈ 4 rtol=0.2
        @test d.location ≈ 3 rtol=0.1
    end

    for α in 1.1:0.1:1.9
        d = AlphaSubGaussian(n=96000, α=α)
        x = rand(d)
        x2 = copy(x)
        rand!(d, x2)
        @test x != x2

        d3 = fit(AlphaStable, x)
        @test d3.α ≈ α rtol=0.2
        @test d3.β ≈ 0 atol=0.2
        @test d3.scale ≈ 1 rtol=0.2
        @test d3.location ≈ 0 atol=0.1
    end

    d4 = AlphaSubGaussian(n=96000)
    m = size(d4.R, 1)-1
    x = rand(d4)
    d5 = fit(AlphaSubGaussian, x, m, p=1.0)
    @test d4.α ≈ d5.α rtol=0.1
    @test d4.R ≈ d5.R rtol=0.1

end
# 362.499 ms (4620903 allocations: 227.64 MiB)
# 346.520 ms (4621052 allocations: 209.62 MiB) # StaticArrays in outer fun
# 345.925 ms (4225524 allocations: 167.66 MiB) # tempind to tuple
# 395.606 ms (3637770 allocations: 164.76 MiB) # x1 to SVector
# 336.877 ms (10125987 allocations: 236.71 MiB) # typeassert on subgprt
# 328.315 ms (3636312 allocations: 164.69 MiB) # revert x1 svector
# 320.845 ms (3440006 allocations: 161.71 MiB)
# 210.449 ms (3438629 allocations: 86.64 MiB) # full typeinfo in x creation
#
#
# @code_warntype rand(Random.GLOBAL_RNG, AlphaSubGaussian(n=96000))
#
#
#
#
# d = AlphaSubGaussian(n=96)
# using SpecialFunctions
# rand(d)
# rng = Random.GLOBAL_RNG
#
# α=d.α; R=d.R; n=d.n
# α ∈ 1.10:0.01:1.98 || throw(DomainError(α, "α must lie within `1.10:0.01:1.98`"))
# m = size(R, 1)-1
# funk1 = x -> (2^α)*sin(π*α/2)*gamma((α+2)/2)*gamma((α+x)/2)/(gamma(x/2)*π*α/2)
# funk2 = x -> 4*gamma(x/α)/((α*2^2)*gamma(x/2)^2)
# funkmarg = x -> gamma(x/2)/(gamma((x-1)/2)*sqrt(π))
# c = 1.2
# k1 = (funk1(m), funk1(m+1))
# k2 = (funk2(m), funk2(m+1))
#
# kmarg = funkmarg(m+1)
# onetoendm1 = StaticArrays.SOneTo(size(R,1)-1)
# kappa = det(R)/det(R[onetoendm1, onetoendm1])
# invR = inv(R)
# invRx1 = inv(R[onetoendm1, onetoendm1])
# sigrootx1 = cholesky(R[onetoendm1, onetoendm1]).L
# modefactor = R[end, onetoendm1]'*inv(R[onetoendm1, onetoendm1])
#
# @code_warntype AlphaStableDistributions.subgausscondprobtabulate(α, SVector(1.,2,3,4), 3., invRx1, invR, randn(2,2), 0.1, 0.1, 0.1, [1,2], 0.1, k1, k2, kmarg)
#
#
# @MVector zeros(n)
#
#
#
# using Cthulhu
# d = AlphaSubGaussian(n=96)
# @descend_code_warntype rand(Random.GLOBAL_RNG, d)
#
# using AlphaStableDistributions
# d1 = AlphaStable(α=1.5)
# s = rand(d1, 100000)
# using ThreadsX
# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort)


text:
module SimplePosetRepresentations

# Nothing here yet. We'll do IntervalOrder, SemiOrder, CircleOrder

using SimplePosets, ClosedIntervals

import Base.show
export IntervalOrder, SemiOrder
export Circle, inside, CircleOrder

export poset_builder

### GENERAL MACHINERY

function poset_builder(Xdict::Dict{S,T}, comparator::Function) where {S,T}
  P = SimplePoset{S}()
  for v in keys(Xdict)
    add!(P,v)
  end
  elts = collect(keys(Xdict))
  n = length(elts)
  for i=1:n
    a = elts[i]
    u = Xdict[a]
    for j=1:n
      b = elts[j]
      v = Xdict[b]
      if i!=j && comparator(u,v)
        add!(P,a,b)
      end
    end
  end
  return P
end

function poset_builder(Xlist::Vector{S}, comparator::Function) where S
  d = Dict{Int,S}()
  n = length(Xlist)
  for i=1:n
    d[i] = Xlist[i]
  end
  return poset_builder(d,comparator)
end

function poset_builder(Xset::Set{S}, comparator::Function) where S
  d = Dict{S,S}()
  for s in Xset
    d[s] = s
  end
  return poset_builder(d,comparator)
end

### INTERVAL ORDERS

_interval_compare(I::ClosedInterval{S}, J::ClosedInterval{S}) where S = I << J

"""
`IntervalOrder(Jmap)` creates an interval order from a dictionary
mapping elements to closed intervals.
"""
function IntervalOrder(Jmap::Dict{S,ClosedInterval{T}}) where {S,T}
  return poset_builder(Jmap,_interval_compare)
end

"""
`IntervalOrder(Jlist)` creates an interval order from a list of
intervals. The elements of the poset are named `1:n` where `n`
is the length of the `Jlist`.
"""
function IntervalOrder(Jlist::Vector{ClosedInterval{T}}) where T
  return poset_builder(Jlist,_interval_compare)
end

"""
`IntervalOrder(Jset)` creates an interval order from a set of
closed intervals. The elements of the poset are the intervals.
"""
function IntervalOrder(Jset::Set{ClosedInterval{T}}) where T
  return poset_builder(Jset,_interval_compare)
end

### SEMIORDERS

_unit_compare(x::Real,y::Real) = x+1 < y

"""
`SemiOrder(Xmap)` creates a semiorder whose elements are the
keys in the dictionary `Xmap`. We have `a<b` provided
`Xmap[a]+1 < Xmap[b]`.
"""
function SemiOrder(Xmap::Dict{S,T}) where {S,T<:Real}
  return poset_builder(Xmap, _unit_compare)
end

"""
`SemiOrder(Xlist)` creates a semiorder from a list of real numbers.
The elements are `1:n` and we have `a<b` provided `Xlist[a]+1 < Xlist[b]`.
"""
function SemiOrder(Xlist::Vector{T}) where {T<:Real}
  return poset_builder(Xlist, _unit_compare)
end

"""
`SemiOrder(Xset)` creates a semiorder from a set of real numbers.
The elements of the poset are the elements of `Xset`. For two elements
`a` and `b` we have `a<b` in the poset iff `a+1 < b` as real numbers.
"""
function SemiOrder(Xset::Set{T}) where {T<:Real}
  return poset_builder(Xset,_unit_compare)
end

SemiOrder(Xset::BitSet) = SemiOrder(Set(Xset))

### CIRCLE ORDERS
"""
The `Circle` datatype represents a circle in the
plane specified by a center and a radius like this:
`Circle(x,y,r)`.
"""
struct Circle{X<:Real, Y<:Real, R<:Real}
  x::X     # x-coordinate of center
  y::Y     # y-coordinate of center
  r::R     # radius
end

function show(io::IO, C::Circle)
  print(io, "Circle($(C.x), $(C.y), $(C.r))")
end

"""
`inside(C,D)` returns `true` provided circle `C`
is entirely contained inside circle `D`.
"""
function inside(C::Circle, D::Circle)
  dx = C.x - D.x
  dy = C.y - D.y
  dr = C.r - D.r
  return (dr<=0)&&(dx*dx + dy*dy <= dr*dr)
end

function CircleOrder(Cmap::Dict{S,Circle}) where S
  return poset_builder(Cmap,inside)
end

function CircleOrder(Clist::Vector{Circle})
  return poset_builder(Clist,inside)
end

function CircleOrder(Cset::Set{Circle})
  return poset_builder(Cset,inside)
end

end  # end of module


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function quarterly2monthly(a,n)

    #=
    Filename: monthly2quarterly.jl
    a: quarterly data
    n: length of resulting monthly array
    =#

    if ndims(a) == 1
        nQ = 1;
    else
        nQ = size(a,2);
    end

    a[ismissing.(a)] .= NaN;
    b = hcat([kron(a[:,i], [NaN; NaN; 1]) for i=1:nQ]...) |> Array{Union{Float64, Missing}};
    b[isnan.(b)] .= missing;

    if size(b,1) < n
        b = [b; missing.*ones(n-size(b,1), nQ)];
    end

    return b;
end


text:
using InflationFunctions
using Documenter

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

makedocs(;
    modules=[InflationFunctions],
    authors="Rodrigo Chang <rrcp777@gmail.com> and contributors",
    repo="https://github.com/DIE-BG/InflationFunctions.jl/blob/{commit}{path}#{line}",
    sitename="InflationFunctions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://die-bg.github.io/InflationFunctions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DIE-BG/InflationFunctions.jl",
    devbranch="main",
)


text:
using PowerModels
using .PP2PM

function run_powermodels(json_path)
    pm = PP2PM.load_pm_from_json(json_path)

    model = PP2PM.get_model(pm["pm_model"])

    solver = PP2PM.get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    result = PowerModels.run_pf(pm, model, solver,
                                    setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

max_stars_repo_name: bergkvist/pandapower

max_issues_repo_name: bergkvist/pandapower

max_forks_repo_name: bergkvist/pandapower


text:
function compare_curves(args...; tmin = 1e-8, tmax = 1, npoints = 50)
    ts = 10 .^ range(log(10, tmin), stop = log(10, tmax), length = npoints)

    curves = Dict(Symbol(r) => [] for r in args if r isa Function)
    curves[:slope2] = []
    curves[:slope3] = []

    for (i, t) in enumerate(ts)
        for curve in args
            if curve isa Function
                val = curve(t)
                @assert val >= 0 "curve should be positive, value is $(curve(t))."
                push!(curves[Symbol(curve)], (t, val))
            end
        end
        push!(curves[:slope2], (t, t^2))
        push!(curves[:slope3], (t, t^3))
    end
    return curves
end

function display_curvescomparison(curves)
    ps = []

    for (k, v) in curves
        legendentry = string(k)
        opts = occursin("slope", legendentry) ? Dict("no marks" => nothing) : Dict()
        push!(ps, PlotInc(PGFPlotsX.Options(opts...), Coordinates(v)))
        push!(ps, LegendEntry(replace(string(k), "_" => " ")))
    end

    ## Plot errors:
    return fig_time_subopt = @pgf Axis(
        {
            xmode = "log",
            ymode = "log",
            height = "12cm",
            width = "12cm",
            xlabel = "t",
            # ylabel = L"f(R_x(t\eta))-f(x)-t\langle gradf(x),\eta\rangle",
            legend_pos = "outer north east",
        },
        ps...,
    )
end

function check_regularizer_gradient_hessian(M, reg; x=nothing, ξ=nothing, η=nothing)

    isnothing(x) && (x = randomMPoint(M))
    isnothing(ξ) && (ξ = randomTVector(M, x))
    isnothing(η) && (η = randomTVector(M, x))

    f(x) = g(reg, x)
    gradf_x = ∇M_g(reg, M, x)
    hessf_x_ξ = ∇²M_g_ξ(reg, M, x, ξ)
    hessf_x_η = ∇²M_g_ξ(reg, M, x, η)

    # 1. Gradient is tangent.
    println("- gradf(x) ∈ T_x M:\t\t\t\t", is_tangent_vector(M, x, gradf_x; atol = 1e-10))

    # 2. Hessian is symetric
    println(
        "- Hess f(x)[η] ∈ T_x M:\t\t\t\t",
        is_tangent_vector(M, x, hessf_x_η, atol = 1e-10),
    )
    println(
        "- Hess f(x)[ξ] ∈ T_x M:\t\t\t\t",
        is_tangent_vector(M, x, hessf_x_ξ, atol = 1e-10),
    )
    println(
        "- ⟨Hess f(x)[ξ], η⟩ - ⟨Hess f(x)[η], ξ⟩:\t",
        inner(M, x, hessf_x_η, ξ) - inner(M, x, hessf_x_ξ, η),
    )

    f_x = f(embed(M, x))
    ηgradfx = inner(M, x, η, gradf_x)
    ηHessf_xη = inner(M, x, η, hessf_x_η)

    @show ηgradfx, ηHessf_xη

    frgrad(t) = abs(f(embed(M, retract(M, x, t * η))) - (f(embed(M, x)) + t * ηgradfx))
    function frhess(t)
        return abs(
            f(embed(M, retract(M, x, t * η))) - (f_x + t * ηgradfx + 0.5 * t^2 * ηHessf_xη),
        )
    end

    return compare_curves(frgrad, frhess)
end

function check_e2r_gradient_hessian(M)
    x = randomMPoint(M)
    p = embed(M, retract(M, x, randomTVector(M, x)))

    ## function, rgrad, rhess
    f(x) = exp(sum(x .* (x .- p) .^ 2))

    function get_rgrad(M::Manifold, x)
        egrad = ForwardDiff.gradient(f, embed(M, x))
        return egrad_to_rgrad(M, x, egrad)
    end

    function get_rhess(M::Manifold, x, η)
        x_emb = embed(M, x)
        η_emb = embed(M, x, η)
        ∇f_x = ForwardDiff.gradient(f, x_emb)
        ∇²f_x_η =
            reshape(ForwardDiff.hessian(f, x_emb) * vec(η_emb), representation_size(M))
        return ehess_to_rhess(M, x, ∇f_x, ∇²f_x_η, η)
    end

    ## Gradient check
    η = randomTVector(M, x)

    # f1(t) =
    # f_x = f(x)
    gradf_x = get_rgrad(M, x)

    # 1. gradf ∈ T_x M
    println("- gradf(x) ∈ T_x M:\t\t\t\t", is_tangent_vector(M, x, gradf_x; atol = 1e-10))
    println("- η ∈ T_x M:\t\t\t\t\t", is_tangent_vector(M, x, η; atol = 1e-10))

    # 2. Hessian is symetric
    η = randomTVector(M, x)
    ξ = randomTVector(M, x)

    hessf_x_η = get_rhess(M, x, η)
    hessf_x_ξ = get_rhess(M, x, ξ)

    println(
        "- Hess f(x)[η] ∈ T_x M:\t\t\t\t",
        is_tangent_vector(M, x, hessf_x_η, atol = 1e-10),
    )
    println(
        "- Hess f(x)[ξ] ∈ T_x M:\t\t\t\t",
        is_tangent_vector(M, x, hessf_x_ξ, atol = 1e-10),
    )
    println(
        "- ⟨Hess f(x)[ξ], η⟩ - ⟨Hess f(x)[η], ξ⟩:\t",
        inner(M, x, hessf_x_η, ξ) - inner(M, x, hessf_x_ξ, η),
    )

    f_x = f(embed(M, x))
    ηgradfx = inner(M, x, η, gradf_x)
    ηHessf_xη = inner(M, x, get_rhess(M, x, η), η)

    @show ηgradfx, ηHessf_xη

    frgrad(t) = abs(f(embed(M, retract(M, x, t * η))) - (f(embed(M, x)) + t * ηgradfx))
    function frhess(t)
        return abs(
            f(embed(M, retract(M, x, t * η))) - (f_x + t * ηgradfx + 0.5 * t^2 * ηHessf_xη),
        )
    end

    return compare_curves(frgrad, frhess)
end

function check_retraction(M)
    Random.seed!(1234)
    x = randomMPoint(M)
    Random.seed!(4312)
    η = randomTVector(M, x)

    x_emb = embed(M, x)
    η_emb = embed(M, x, η)

    function retractiontη(t)
        retracted = retract(M, x, t * η)
        return embed(M, retracted)
    end
    xtη(t) = x_emb + t * η_emb

    function retractionerror(t; x = x)
        return norm(project(M, x, retractiontη(t) - xtη(t)))
    end

    return compare_curves(retractionerror)
end

function remove_small_functionvals(values; threshold = 1e-12)
    return filter(x -> x[2] >= threshold, values)
end

function get_sloperesidual(data; threshold = 1e-12)
    n = length(data)
    x = [data[i][1] for i in 1:n]
    logy = [log(data[i][2]) for i in 1:n]

    ## Fit linear regressor
    # model : y = βx + c
    # explain Y = [log.(y) 1] by X = [log.(x)]
    # Y = X β
    # β ∈ ℝ 2x1 - slope and absciss at origin
    # X ∈ ℝ nx2 - absciss and intercept
    # Y ∈ ℝ nx1 - ordonate
    X = ones(n, 2)
    X[:, 1] = log.(x)
    Y = logy

    F = factorize(X' * X)
    β = F \ X' * Y
    predictedslope = β[1]
    residual = 0.5 * norm((Y - X * β)[:, 1])^2
    return predictedslope, residual
end

export get_sloperesidual, remove_small_functionvals


text:
using Sobol
using Roots

using Random
using Distributions
using SparseArrays
using LinearAlgebra

"""
Describe meta particules, represented by a Dirac disbtibution in (x, v), with a weight (~a high) of wei
"""
struct Particles
    x ::Vector{Float64}     # list of the positions
    v ::Vector{Float64}     # list of the velocities
    wei ::Vector{Float64}   # list of the weights of the particules
    nbpart ::Int64          # nmber of particules
end

"""
    Defines Runge-Kutta-Nystrom time integrator via its Butcher tableau,
    and holds some pre-allocated arrays used for the time integration only.
"""
struct rkn_order4
    a :: Array{Float64, 2}
    b̄ :: Vector{Float64}
    c :: Vector{Float64}
    b :: Vector{Float64}
    dt :: Float64
    fg ::       Array{Float64, 2}
    G ::        Array{Float64, 1}

    function rkn_order4(X, dt)
        # a, b̄, c, b correspond to the Butcher tableau of Runge-Kutta-Nystrom 3steps order4.
        a = [0.0        0.0       0.0;
            (2-√3)/12   0.0           0.0;
            0.0         √(3)/6      0.0]
        b̄ = [(5 - 3*√3)/24,     (3+√3)/12,  (1+√3)/24]
        c = [(3+√3)/6,          (3-√3)/6,   (3+√3)/6]
        b = [(3-2*√3)/12,       1/2,        (3+2*√3)/12]

        new(a .* dt^2, b̄ .* dt^2, c .* dt, b .* dt, dt,
            zeros(Float64, length(X), 3),  # fg
            similar(X) # G
        )
    end
end

"""
Holds pre-allocated arrays
"""
struct ParticleMover
    Φ ::Vector{Float64}
    ∂Φ :: Vector{Float64}
    meshx
    kx
    K
    C ::Vector{Float64}
    S ::Vector{Float64}
    tmpcosk ::Vector{Float64}
    tmpsink ::Vector{Float64}
    tmpcoskimplicit ::Vector{Float64}
    tmpsinkimplicit ::Vector{Float64}
    poisson_matrix
    ρ ::Vector{Float64}
    Φ_grid ::Vector{Float64}
    idxonmesh :: Vector{Int64}
    idxonmeshp1 :: Vector{Int64}
    rkn :: rkn_order4
    dt :: Float64

    function ParticleMover(particles::Particles, meshx, K, dt; kx=1)
        Φ = similar(particles.x)
        ∂Φ = similar(particles.x)
        tmpcosk = similar(particles.x)
        tmpsink = similar(particles.x)
        nx = meshx.len
        np = particles.nbpart

        matrix_poisson = spdiagm(  -1 => .- ones(Float64,nx-1),
                                    0 => 2 .* ones(Float64,nx),
                                    1 => .- ones(Float64,nx-1))
        matrix_poisson[1, nx] = -1
        matrix_poisson[nx, 1] = -1

        matrix_poisson ./= meshx.step^2

        new(Φ, ∂Φ, meshx, kx, K, Vector{Float64}(undef, K), Vector{Float64}(undef, K), tmpcosk, tmpsink, similar(tmpcosk), similar(tmpsink), matrix_poisson, Vector{Float64}(undef, nx), Vector{Float64}(undef, nx), Vector{Int64}(undef, np), Vector{Int64}(undef, np), rkn_order4(particles.x, dt), dt)
    end
end

#==== Time steppers ====#
"""symplectic_RKN_order4!(X, V, F, rkn, kx)

    Advect (X, V) on a time step dt using symplectic Runge-Kutta-Nystrom method of order4 [Feng, Qin (2010), sect.7.3, p.327, scheme1].

    The equation satisfied by X is
    ```math
    \\frac{d^2 X(t)}{dt^2} = C(t)\\cos(X(t)) - S(t)\\sin(X(t))
    ```

    RKN method considers Ẋ = V as a variable, and updates both X and V.

    Args:
    - X: matrix of positions at time t_n
    - V: matrix of velocities at time t_n
    - F: values of initial condition at time t_0
    - rkn: rkn_order_4 struct, storing butcher tableau and pre-allocated arrays (holds the value of dt)
    - kx: 2π/L

    Updates X, V in place, and returns coefficients C, S at current time.

"""
function symplectic_RKN_order4!(p, pmover, kernel)
    @views begin
        for s=1:3
            @. pmover.rkn.G = p.x + p.v * pmover.rkn.c[s] + pmover.rkn.a[s, 1] * pmover.rkn.fg[:, 1] +
                                                             pmover.rkn.a[s, 2] * pmover.rkn.fg[:, 2] +
                                                             pmover.rkn.a[s, 3] * pmover.rkn.fg[:, 3]

            kernel(pmover.rkn.fg[:, s], pmover.rkn.G, p, pmover; coeffdt=0)
        end

        @. p.x += pmover.dt * p.v + pmover.rkn.b̄[1] * pmover.rkn.fg[:, 1] +
                                    pmover.rkn.b̄[2] * pmover.rkn.fg[:, 2] +
                                    pmover.rkn.b̄[3] * pmover.rkn.fg[:, 3]
        @. p.v += pmover.rkn.b[1] * pmover.rkn.fg[:, 1] +
                    pmover.rkn.b[2] * pmover.rkn.fg[:, 2] +
                    pmover.rkn.b[3] * pmover.rkn.fg[:, 3]

        kernel(pmover.∂Φ, p.x, p, pmover)
    end
end

"""strang_splitting!(particles, pmover, kernel)

    Other method for advecting (X, V) on a time step.

    Uses Verlet scheme (of order 2).

    Args:
    - particles: Particle struct
    - pmover: ParticleMover struct
    - kernel: how to compute the flow field.

    Updates X, V in place, and returns coefficients C, S at current time.
"""
function strang_splitting!(particles, pmover, kernel)
    pmover.Φ .= 0
    pmover.∂Φ .= 0
    @. particles.x += particles.v * pmover.dt / 2
    kernel(pmover.∂Φ, particles.x, particles, pmover)
    @. particles.v += pmover.dt * pmover.∂Φ
    @. particles.x += particles.v * pmover.dt / 2
    kernel(pmover.∂Φ, particles.x, particles, pmover)
end

function strang_splitting_implicit!(particles, pmover, kernel)
    pmover.Φ .= 0
    pmover.∂Φ .= 0
    κ = 1/2
    kernel(pmover.∂Φ, particles.x, particles, pmover; coeffdt=κ)
    @. particles.v += pmover.dt / 2 * pmover.∂Φ
    @. particles.x += particles.v * pmover.dt
    kernel(pmover.∂Φ, particles.x, particles, pmover; coeffdt=κ)
    @. particles.v += pmover.dt / 2 * pmover.∂Φ
    kernel(pmover.∂Φ, particles.x, particles, pmover)
end

# ===== Kernel computations ==== #
"""
Compute -∂_x Φ[f](`x` + `pmover.dt*coeffdt*p.v`) and stores it in `dst`. Also updates `pmover.Φ`
"""
function kernel_poisson!(dst, x, p, pmover; coeffdt=0)
    dst .= 0

    if coeffdt > 0
        pmover.tmpcoskimplicit .= similar(pmover.tmpcosk)
        pmover.tmpsinkimplicit .= similar(pmover.tmpsink)
    end

    pmover.Φ .= 0

    for k=1:pmover.K
        @. pmover.tmpcosk = cos(x * (k * pmover.kx))
        @. pmover.tmpsink = sin(x * (k * pmover.kx))

        if coeffdt > 0
            @. pmover.tmpcoskimplicit = cos((x + pmover.dt * coeffdt * p.v) * (k * pmover.kx))
            @. pmover.tmpsinkimplicit = sin((x + pmover.dt * coeffdt * p.v) * (k * pmover.kx))
        end

        pmover.C[k] = sum(pmover.tmpcosk .* p.wei)
        pmover.S[k] = sum(pmover.tmpsink .* p.wei)

        if coeffdt > 0
            pmover.C[k] += -sum(p.v .* pmover.tmpsink .* p.wei) * k * 2π / pmover.meshx.stop * pmover.dt * coeffdt
            pmover.S[k] +=  sum(p.v .* pmover.tmpcosk .* p.wei) * k * 2π / pmover.meshx.stop * pmover.dt * coeffdt
        end

        if coeffdt > 0
            pmover.tmpcosk .= pmover.tmpcoskimplicit
            pmover.tmpsink .= pmover.tmpsinkimplicit
        end

        @. pmover.Φ += (pmover.C[k] * pmover.tmpcosk + pmover.S[k] * pmover.tmpsink) / k^2
        # The line below computes -∂Φ[f](`x`) and stores it to `dst`.
        # Changing dynamics :
        #   "+=": repulsive potential (plasmas dynamics)
        #   "-=": attractive potential (galaxies dynamics)
        @. dst += (pmover.C[k] * pmover.tmpsink - pmover.S[k] .* pmover.tmpcosk) / k
    end

    dst ./= π
    pmover.Φ .*= pmover.meshx.stop / (2*π^2)
end

function kernel_gyrokinetic!(dst, x, p, pmover)
    dst .= .-x
end

# ===== Some quantities we can compute at each step ==== #
"""compute_electricalenergy²(p, pmover)

    Returns the square of the electrical energy
"""
function compute_electricalenergy²(p, pmover)
    return sum(p.wei .* pmover.Φ)
end

function compute_momentum(particles)
    return sum(particles.wei .* particles.v)
end

function compute_totalenergy²(particles, Eelec²)
    return 1/2 * (Eelec² .+ sum(particles.v.^2 .* particles.wei))
end

""" periodic_boundary_conditions!(p, pmover)

    Impose periodic boundary conditions in space.
"""
function periodic_boundary_conditions!(p, pmover)
    p.x[findall(>(pmover.meshx.stop), p.x)] .-= pmover.meshx.stop
    p.x[findall(<(0)                , p.x)] .+= pmover.meshx.stop
end

function PF_step!(p::Particles, pmover::ParticleMover; kernel=kernel_poisson!)
    symplectic_RKN_order4!(p, pmover, kernel)
    # strang_splitting!(p, pmover, kernel)
    # strang_splitting_implicit!(p, pmover, kernel)

    periodic_boundary_conditions!(p, pmover)

    E² = compute_electricalenergy²(p, pmover)
    return E², compute_momentum(p), compute_totalenergy²(p, E²)
end

function PIC_step!(p, meshx, meshv, dt, dΦdx)
    #=
    PIC step when supplying the derivative of the electrostatic potential
    =#
    # Use a 3-step splitting, to be of order 2.
    # We interpolate the derivative of the potential, defined on a grid, to obtain an approximation
    # of its value at each particle position.

    idxonmesh = Int64.(fld.(p.x, meshx.step)) .+ 1
    t = (p.x .- (idxonmesh.-1) .* meshx.step) ./ meshx.step
    idxonmesh[findall(i -> i > meshx.len,  idxonmesh)] .-= meshx.len
    idxonmesh[findall(i -> i < 1               ,  idxonmesh)] .+= meshx.len
    idxonmeshp1 = idxonmesh .+ 1
    idxonmeshp1[findall(i -> i > meshx.len,  idxonmeshp1)] .-= meshx.len
    p.v .-= dt ./ 2 .* (dΦdx[idxonmesh] .* (1 .- t) .+ dΦdx[idxonmeshp1] .* t)

    # for ipart = 1:p.nbpart
    #     idxgridx = Int64(fld(p.x[ipart], meshx.step)) + 1
    #     t = (p.x[ipart] - (idxgridx-1) * meshx.step) / meshx.step
    #     p.v[ipart] -= dt/2 * (dΦdx[idxgridx] * (1-t) + dΦdx[idxgridx < meshx.len ? idxgridx+1 : 1] * t)
    #     # # Periodic boundary conditions in velocity
    #     # if p.v[ipart] > meshv.stop
    #     #     p.v[ipart] -= meshv.stop - meshv.start
    #     # elseif p.v[ipart] < meshv.start
    #     #     p.v[ipart] += meshv.stop - meshv.start
    #     # end
    # end
    update_positions!(p, meshx, dt)
    # for ipart = 1:p.nbpart
    #     idxgridx = Int64(fld(p.x[ipart], meshx.step)) + 1
    #     t = (p.x[ipart] - (idxgridx-1) * meshx.step) / meshx.step
    #     p.v[ipart] -= dt/2 * (dΦdx[idxgridx] * (1-t) + dΦdx[idxgridx < meshx.len ? idxgridx+1 : 1] * t)
    #     # Periodic boundary conditions in velocity
    #     # if p.v[ipart] > meshv.stop
    #     #     p.v[ipart] -= meshv.stop - meshv.start
    #     # elseif p.v[ipart] < meshv.start
    #     #     p.v[ipart] += meshv.stop - meshv.start
    #     # end
    # end

    idxonmesh .= Int64.(fld.(p.x, meshx.step)) .+ 1
    t = (p.x .- (idxonmesh.-1) .* meshx.step) ./ meshx.step
    idxonmesh[findall(i -> i > meshx.len,  idxonmesh)] .-= meshx.len
    idxonmesh[findall(i -> i < 1               ,  idxonmesh)] .+= meshx.len
    idxonmeshp1 .= idxonmesh .+ 1
    idxonmeshp1[findall(i -> i > meshx.len,  idxonmeshp1)] .-= meshx.len
    p.v .-= dt ./ 2 .* (dΦdx[idxonmesh] .* (1 .- t) .+ dΦdx[idxonmeshp1] .* t)

    return sqrt.(sum(dΦdx.^2) * meshx.stop / p.nbpart)
end

#= ================= #
    CLASSICAL PIC
#  ================ =#

function samples(nsamples, kx, α::Float64, μ::Float64, β::Float64)
    #=
    Sample distribution defined by 1+α\cos(kx * x) in space and a gaussian of variance beta and mean μ in velocity.
    =#
    x0 = Array{Float64}(undef, nsamples)
    seq = SobolSeq(1)
    for i = 1:nsamples
        r = Sobol.next!(seq)
        resx = find_zero(x-> x + α/kx * sin(kx*x) - r[1]*2π/kx, 0.0)
        x0[i] = resx
    end
    v0 = rand(Normal(μ, 1/√β), nsamples)
    wei = (2π/kx)/nsamples .* ones(Float64, nsamples)
    return x0, v0, wei
end

function update_positions!(p, mesh, dt)
    p.x .+= p.v .* dt
    p.x[findall(x -> x >= mesh.stop,  p.x)] .-= mesh.stop - mesh.start
    p.x[findall(x -> x <  mesh.start, p.x)] .+= mesh.stop - mesh.start
end

"""
update particle velocities vp (Φ_v)
"""
function update_velocities!(p, pmover, dt)
    compute_S_C!(p, pmover)
    pmover.Φ .= 0
    pmover.∂Φ .= 0

    for k = 1:pmover.K
        pmover.tmpcosk .= cos.(k .* pmover.kx .* p.x)
        pmover.tmpsink .= sin.(k .* pmover.kx .* p.x)
        pmover.Φ   .+= (pmover.tmpcosk .* pmover.C[k] .+ pmover.tmpsink .* pmover.S[k]) ./ k^2
        pmover.∂Φ  .+= (.-pmover.tmpsink .* pmover.C[k] .+ pmover.tmpcosk .* pmover.S[k]) ./ k
    end
    pmover.Φ .*= pmover.meshx.stop / 2*π^2
    pmover.∂Φ ./= π
    p.v .-= dt .* pmover.∂Φ
end

""" S[k] = \\sum_l=1^n {\\beta_l * sin(k kx x_l)} et C[k] = \\sum_l=1^n {\\beta_l * cos(k kx x_l)}
utile pour le calcul de la vitesse et du potentiel electrique Φ"""
function compute_S_C!(p, pmover)
    pmover.S .= 0
    pmover.C .= 0

    for k = 1:pmover.K
        pmover.S[k] = sum(p.wei .* sin.(k .* 2π ./ pmover.meshx.stop .* p.x))
        pmover.C[k] = sum(p.wei .* cos.(k .* 2π ./ pmover.meshx.stop .* p.x))
    end
end

# ##### functions used to perform projection on grid. ##### #

"""
compute rho, charge density (ie int f dv)

Projection on mesh is done here.
"""
function compute_rho!(p, pmover)
    dx = pmover.meshx.step
    pmover.ρ .= 0.0

    pmover.idxonmesh .= Int64.(fld.(p.x, dx)) .+ 1
    t = (p.x .- (pmover.idxonmesh.-1) .* dx) ./ dx

    pmover.idxonmesh[findall(i -> i > pmover.meshx.len,  pmover.idxonmesh)] .-= pmover.meshx.len
    pmover.idxonmesh[findall(i -> i < 1               ,  pmover.idxonmesh)] .+= pmover.meshx.len

    pmover.idxonmeshp1 .= pmover.idxonmesh .+ 1
    pmover.idxonmeshp1[findall(i -> i > pmover.meshx.len,  pmover.idxonmeshp1)] .-= pmover.meshx.len

    pmover.ρ[pmover.idxonmesh]   .+= p.wei .* (1 .- t)
    pmover.ρ[pmover.idxonmeshp1] .+= p.wei .* t

    println((sum(pmover.ρ.^2), sum((p.wei .* t).^2), sum((p.wei .* (1 .-t)).^2)))

    pmover.ρ .-= sum(pmover.ρ) / pmover.meshx.stop
    return pmover.ρ
end

function compute_Φ!(pmover) # solving poisson equation with FD solver
    dx = pmover.meshx.step
    L = pmover.meshx.stop

    pmover.Φ_grid .= pmover.poisson_matrix \ pmover.ρ
    pmover.Φ_grid .-= sum(pmover.Φ_grid) * dx / L
end

function compute_int_E(p, pmover) # energy computed with a projection on the grid
    compute_rho!(p, pmover)
    compute_Φ!(pmover)
    E = compute_E(pmover)
    return sum(E.^2) * pmover.meshx.step
end

function compute_E(pmover)
    return -(circshift(pmover.Φ_grid, 1) .- circshift(pmover.Φ_grid, -1)) ./ (2*pmover.meshx.step)
end

# Carybe function
"""
compute rho, charge density (ie int f dv)
"""
function compute_rho( p, m )

   L = m.stop
   nx = m.len
   dx = m.step
   rho = zeros(nx)

   for ipart=1:p.nbpart
        xp = p.x[ipart]
        i = Int(floor(Real(xp) / dx) + 1) #number of the cell with the particule
        dum = p.wei[ipart] / dx
        mx_i = (i-1) * dx
        mx_j = i * dx
        a1 = (mx_j-xp) * dum
        a2 = (xp-mx_i) * dum
        if i < 1
            i += nx
        elseif i > nx
            i -= nx
        end
        rho[i]   += a1
        rho[i < nx ? i+1 : i+1 - nx] += a2
   end

#    rho[nx] += rho[1]
#    rho[1] = rho[nx]

   rho ./= dx

   rho_total  = sum(rho) * dx / L
   return rho, rho_total
end


text:
using DataStructures

unscrambled = Dict{String,Int}(
    Set(collect("abcefg")) => 0,
    Set(collect("cf")) => 1,
    Set(collect("acdeg")) => 2,
    Set(collect("acdfg")) => 3,
    Set(collect("bcdf")) => 4,
    Set(collect("abdfg")) => 5,
    Set(collect("abdefg")) => 6,
    Set(collect("acf")) => 7,
    Set(collect("abcdefg")) => 8,
    Set(collect("abcdfg")) => 9
)

by_count = DefaultDict(Set)
for (k,v) in unscrambled
    push!(by_count[length(k)], k)
end

by_count

parse_file(filename) = [map(split,split(line, " | ")) for line in readlines(filename)]

parse_file("8eg.txt")[1]

function simple(s)
    l = length(s)
    l == 2 || l == 4 || l == 3 || l == 7
end

function simple2(s)
    length(s) in [2,4,3,7]
end

function lookup_simple(s)
    l = length(s)
    result = Dict{Char,Char}()
    @assert length(by_count[l]) == 1
    target = collect(by_count[l])[1]
    for (idx, char) in enumerate(s)
        result[char] = target[idx]
    end
    result
end

lookup_simple("ae")

solve1(filename) = reduce(+, [length(filter(simple, output)) for (egs, output) in parse_file(filename)])

solve1("8eg.txt")
solve1("8.txt")

function possible_mapping(eg)
    result = DefaultDict(Set{Char})
    l = length(eg)
    targets = by_count[l]
    for target in targets
        for char in eg
            for target_char in target
                push!(result[char], target_char)
            end
        end
    end
    result
end

function get_result(mapping, output)
    display(mapping)
    to_lookup = map(x->mapping[x], collect(output))
    unscrambled[to_lookup]
end

all_chars = collect("abcdefg")
function solveline(egs, outputs)
    mapping = Dict()
    for eg in egs
        if simple(eg)
            for (k,v) in lookup_simple(eg)
                mapping[k] = v
            end
        end
    end
    @assert length(mapping) == 7
    [(mapping, output) for output in outputs]
end

[solveline(line...) for line in parse_file("8eg.txt")]

mapping = Dict('f' => 'e', 'g' => 'f', 'c' => 'c', 'd' => 'c', 'e' => 'd', 'a' => 'g', 'b' => 'b')
eg = "gbdfcae"

map(x->mapping[x], collect(eg))

get_result(mapping, eg)


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

export FinGenGrpAb, FinGenGrpAbElem, parent, is_finite, is_infinite, rank,
       getindex, show, +, *, call, ngens, snf_with_transform, nrels,
       -, ==, is_trivial, order, exponent, AbelianGroup, DiagonalGroup,
       quo, sub

import Base.+, Nemo.snf, Nemo.parent, Base.rand

################################################################################
#
#  Parent and element types IO, types are in HeckeTypes
#
################################################################################

function show(io::IO, A::FinGenGrpAbGen)
  print(io, "(General) abelian group with relation matrix $(A.rels)")
  if isdefined(A, :snf_map)
    println(io, "\nwith structure of ", codomain(A.snf_map))
  end
end

function show(io::IO, A::FinGenGrpAbSnf)
  println(io, "Abelian group")
  nz = 0
  for i=1:length(A.snf)
    if A.snf[i] == 0
      nz += 1
    end
  end

  if A.snf[1] != 0
    print(io, "Z/$(A.snf[1])")
  end

  i = 2
  while i<= length(A.snf) && A.snf[i] != 0
    print(io, " x Z/$(A.snf[i])")
    i += 1
  end
  if i== length(A.snf)
    print(io, " x Z")
  elseif i < length(A.snf)
    print(io, " x Z^$(length(A.snf)-i+1)")
  end
end

function show(io::IO, a::FinGenGrpAbElem)
  print(io, "Element of\n$(a.parent)\n with components\n$(a.coeff)")
#  print(io, a.coeff)
end

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

doc"""
***
    parent(x::FinGenGrpAbElem ) -> FinGenGrpAb

> Returns the parent of $x$.
"""
function parent(x::FinGenGrpAbElem)
  return x.parent
end

doc"""
***
    ngens(G::FinGenGrpAbSnf ) -> Int

> Returns the number of generators of G in the current representation.
"""
function ngens(A::FinGenGrpAbGen)
  return cols(A.rels)
end

function ngens(A::FinGenGrpAbSnf)
  return length(A.snf)
end

doc"""
***
    nrels(G::FinGenGrpAbSnf ) -> Int

> Returns the number of relations of G in the current representation.
"""
function nrels(A::FinGenGrpAbGen)
  return rows(A.rels)
end

function nrels(A::FinGenGrpAbSnf)
  return length(A.snf)
end

doc"""
***
    getindex(x::FinGenGrpAbElem, i::Int) -> fmpz

> Returns the $i$-th component of the element $x$.
"""
function getindex(x::FinGenGrpAbElem, i::Int)
  return x.coeff[1,i]
end

################################################################################
################################################################################
################################################################################
function assert_hnf(A::FinGenGrpAbGen)
  isdefined(A, :hnf) && return
  A.hnf = hnf(A.rels)
end

function reduce_mod_hnf!(a::fmpz_mat, H::fmpz_mat)
  j = 1
  for i=1:min(rows(H), cols(H))
    while H[i,j] == 0 && j <= cols(H)
      j+=1
    end
    if j > cols(H)
      return
    end
    q = fdiv(a[1, j], H[i,j])
    for k=j:cols(a)
      a[1,k] =  a[1, k] -q*H[i,k]
    end
  end
end

function FinGenGrpAbElemCreate(A::FinGenGrpAbGen, a::fmpz_mat)
  assert_hnf(A::FinGenGrpAbGen)
  reduce_mod_hnf!(a, A.hnf)
  return FinGenGrpAbElem{FinGenGrpAbGen}(A, a)
  # reduce mod hnf
end

function FinGenGrpAbElemCreate(A::FinGenGrpAbSnf, a::fmpz_mat)
  for i=1:ngens(A)
    if A.snf[i] != 0
      a[1,i] = a[1,i] % A.snf[i]
    end
  end
  return FinGenGrpAbElem{FinGenGrpAbSnf}(A, a)
end
################################################################################
#
# Arithmetic
#
################################################################################

function ==(a::FinGenGrpAbElem, b::FinGenGrpAbElem)
  a.parent == b.parent || error("elements must belong to the same group")
  return a.coeff == b.coeff
end

doc"""
***
    +(x::FinGenGrpAbElem, y::FinGenGrpAbElem) -> FinGenGrpAbElem

> Computes $x + y$.
"""
function +(x::FinGenGrpAbElem, y::FinGenGrpAbElem)
  x.parent == y.parent || error("elements must belong to the same group")
  n = FinGenGrpAbElemCreate(x.parent, x.coeff + y.coeff)
  return n
end

doc"""
***
    -(x::FinGenGrpAbElem, y::FinGenGrpAbElem) -> FinGenGrpAbElem

> Computes $x - y$.
"""
function -(x::FinGenGrpAbElem, y::FinGenGrpAbElem)
  x.parent == y.parent || error("elements must belong to the same group")
  n = FinGenGrpAbElemCreate(x.parent, x.coeff - y.coeff)
  return n
end

doc"""
***
    -(x::FinGenGrpAbElem) -> FinGenGrpAbElem

> Computes $-x$.
"""
function -(x::FinGenGrpAbElem)
  n = FinGenGrpAbElemCreate(x.parent, -x.coeff)
  return n
end

doc"""
***
    *(x::fmpz, y::FinGenGrpAbElem) -> FinGenGrpAbElem

> Computes $x \cdot y$.
"""
function *(x::fmpz, y::FinGenGrpAbElem)
  n = x*y.coeff
  return FinGenGrpAbElemCreate(y.parent, n)
end

doc"""
***
    *(x::Integer, y::FinGenGrpAbElem) -> FinGenGrpAbElem

> Computes $x \cdot y$.
"""
function *(x::Integer, y::FinGenGrpAbElem)
  n = x*y.coeff
  return FinGenGrpAbElemCreate(y.parent, n)
end

################################################################################
#
#  Parent object overloading
#
################################################################################

doc"""
***
    call(A::FinGenGrpAbSnf, x::Array{fmpz, 1}) -> FinGenGrpAbElem

> Given an array of elements of type `fmpz` of the same length as ngens($A$),
> this function returns the element of $A$ with components `x`.
"""
function call(A::FinGenGrpAbSnf, x::Array{fmpz, 1})
  ngens(A) != length(x) && error("Lengths do not coincide")
  y = MatrixSpace(ZZ, 1, ngens(A))(x)
  z = FinGenGrpAbElemCreate(A, y)
  return z
end

doc"""
***
    call{T <: Integer}(A::FinGenGrpAbSnf, x::Array{T, 1}) -> FinGenGrpAbElem

> Given an array of elements of type `Integer` of the same length as ngens($A$),
> this function returns the element of $A$ with components `x`.
"""
function call{T <: Integer}(A::FinGenGrpAbSnf, x::Array{T, 1})
  ngens(A) != length(x) && error("Lengths do not coincide")
  z = A(map(fmpz, x))
  return z
end

doc"""
***
    call(A::FinGenGrpAbSnf, x::fmpz_mat) -> FinGenGrpAbElem

> Given a 1 x ngens(A) fmpz_mat,
> this function returns the element of $A$ with components `x`.
"""
function call(A::FinGenGrpAbSnf, x::fmpz_mat)
  ngens(A) != cols(x) && error("Lengths do not coincide")
  rows(x) != 1 && error("Matrix should have only one row")
  z = FinGenGrpAbElemCreate(A, x)
  return z
end

doc"""
***
    call(A::FinGenGrpAbGen, x::Array{fmpz, 1}) -> FinGenGrpAbElem

> Given an array of elements of type `fmpz` of the same length as ngens($A$),
> this function returns the element of $A$ with components `x`.
"""
function call(A::FinGenGrpAbGen, x::Array{fmpz, 1})
  ngens(A) != length(x) && error("Lengths do not coincide")
  y = MatrixSpace(ZZ, 1, ngens(A))(x)
  z = FinGenGrpAbElemCreate(A, y)
  return z
end

doc"""
***
    call{T <: Integer}(A::FinGenGrpAbGen, x::Array{T, 1}) -> FinGenGrpAbElem

> Given an array of elements of type `Integer` of the same length as ngens($A$),
> this function returns the element of $A$ with components `x`.
"""
function call{T <: Integer}(A::FinGenGrpAbGen, x::Array{T, 1})
  ngens(A) != length(x) && error("Lengths do not coincide")
  z = A(map(fmpz, x))
  return z
end

doc"""
***
    call(A::FinGenGrpAbGen, x::fmpz_mat) -> FinGenGrpAbElem

> Given a 1 x ngens(A) fmpz_mat,
> this function returns the element of $A$ with components `x`.
"""
function call(A::FinGenGrpAbGen, x::fmpz_mat)
  ngens(A) != cols(x) && error("Lengths do not coincide")
  rows(x) != 1 && error("Matrix should have only one row")
  z = FinGenGrpAbElemCreate(A, x)
  return z
end

doc"""
***
    getindex(A::FinGenGrpAbSnf, i::Int) -> FinGenGrpAbElem

> Returns the element of $A$ with components $(0,\dotsc,0,1,0,\dotsc,0)$,
> where the $1$ is at the $i$-th position.
"""
function getindex(A::FinGenGrpAb, i::Int)
  (i<1 || i > ngens(A)) && error("Index out of range")
  z = MatrixSpace(ZZ, 1, ngens(A))()
  for j in 1:ngens(A)
    z[1,j] = fmpz()
  end
  z[1,i] = fmpz(1)
  return A(z)
end

elem_type(G::FinGenGrpAbSnf) = FinGenGrpAbElem{FinGenGrpAbSnf}
elem_type(G::FinGenGrpAbGen) = FinGenGrpAbElem{FinGenGrpAbGen}
elem_type(::Type{FinGenGrpAbGen}) = FinGenGrpAbElem{FinGenGrpAbGen}
elem_type(::Type{FinGenGrpAbSnf}) = FinGenGrpAbElem{FinGenGrpAbSnf}
elem_type(::Type{FinGenGrpAb}) = FinGenGrpAbElem # I have no idea what this
                                    # does, but it appears to be important
                                    # U, m = UnitGroup(ResidueRing(ZZ, 2^10))
                                    # m(U[2])
                                    # preimage(m, ans)
                                    # fails without...

doc"""
***
  snf(A::FinGenGrpAbSnf) -> FinGenGrpAb, Map

> Given some abelian group, return the "same" group in canonical form,
> i.e. as the direct product of cyclic groups with dividing orders.
> The second return value is the map to translate between the new and old
> groups.
"""
function snf(G::FinGenGrpAbSnf)
  return G, IdentityMap(G)
end

doc"""
***
  is_diag(A::fmpz_mat)

> Tests if A is diagonal.
"""
function is_diag(A::fmpz_mat)
  for i = 1:cols(A)
    for j= 1:rows(A)
      if i != j && A[j,i]!=0
        return false
      end
    end
  end
  return true
end

#=
g, e,f = gcdx(a, b)
U = [1 0 ; -divexact(b, g)*f 1]*[1 1; 0 1];
V = [e -divexact(b, g) ; f divexact(a, g)];

then U*[ a 0; 0 b] * V = [g 0 ; 0 l]
=#
doc"""
***
  snf_with_transform(A::fmpz_mat; l::Bool = true, r::Bool = true) -> fmpz_mat [, fmpz_mat [, fmpz_mat]]

> Given some integer matrix A, compute the Smith normal form (elementary
> divisor normal form) of A. If l and/ or r are true, then the corresponding
> left and/ or right transformation matrices are computed as well.
"""
function snf_with_transform(A::fmpz_mat; l::Bool = true, r::Bool = true)
  if r
    R = MatrixSpace(ZZ, cols(A), cols(A))(1)
  end

  if l
    L = MatrixSpace(ZZ, rows(A), rows(A))(1)
  end
  # TODO: if only one trafo is required, start with the HNF that does not
  #       compute the trafo
  #       Rationale: most of the work is on the 1st HNF..
  #
  S = deepcopy(A)
  while !is_diag(S)
    if l
      S, T = hnf_with_transform(S)
      L = T*L
    else
      S = hnf(S)
    end

    if is_diag(S)
      break
    end
    if r
      S, T = hnf_with_transform(S')
      R = T*R
    else
      S = hnf(S')
    end
    S = S'
  end
  #this is probably not really optimal...
  for i=1:min(rows(S), cols(S))
    if S[i,i] == 1
      continue
    end
    for j=i+1:min(rows(S), cols(S))
      if S[j,j] == 0
        continue
      end
      if S[i,i] != 0 && S[j,j] % S[i,i] == 0
        continue
      end
      g, e,f = gcdx(S[i,i], S[j,j])
      a = divexact(S[i,i], g)
      S[i,i] = g
      b = divexact(S[j,j], g)
      S[j,j] *= a
      if l
        # U = [1 0; -b*f 1] * [ 1 1; 0 1] = [1 1; -b*f -b*f+1]
        # so row i and j of L will be transformed. We do it naively
        # those 2x2 transformations of 2 rows should be a c-primitive
        # or at least a Nemo/Hecke primitive
        for k=1:cols(L)
          x = -b*f
#          L[i,k], L[j,k] = L[i,k]+L[j,k], x*L[i,k]+(x+1)*L[j,k]
          L[i,k], L[j,k] = L[i,k]+L[j,k], x*(L[i,k]+L[j,k])+L[j,k]
        end
      end
      if r
        # V = [e -b ; f a];
        # so col i and j of R will be transformed. We do it naively
        # careful: at this point, R is still transposed
        for k=1:rows(R)
          R[i, k], R[j, k] = e*R[i,k]+f*R[j,k], -b*R[i,k]+a*R[j,k]
        end
      end
    end
  end

  if l
    if r
      return S, L, R'
    else
      return S, L
    end
  elseif r
    return S, R'
  else
    return S
  end
end

function snf(G::FinGenGrpAbGen)
  if isdefined(G, :snf_map)
    return codomain(G.snf_map), G.snf_map
  end
  S, T = snf_with_transform(G.rels, l=false, r=true)
  d = fmpz[S[i,i] for i=1:min(rows(S), cols(S))]
  while length(d) < ngens(G)
    push!(d, 0)
  end
  s = Array{fmpz, 1}()
  for i=1:length(d)
    if d[i] != 1
      push!(s, d[i])
    end
  end
  TT = MatrixSpace(ZZ, rows(T), length(s))()
  j = 1
  for i=1:length(d)
    if d[i] != 1
      for k=1:rows(T)
        TT[k, j] = T[k, i]
      end
      j += 1
    end
  end
  TTi = MatrixSpace(ZZ, length(s), rows(T))()
  Ti = inv(T)

  j = 1
  for i=1:length(d)
    if d[i] != 1
      for k=1:rows(T)
        TTi[j, k] = Ti[i, k]
      end
      j += 1
    end
  end
  H = FinGenGrpAbSnf(s)
  mp = FinGenGrpAbMap(G, H, TT, TTi)
  G.snf_map = mp
  return H, mp
end

################################################################################
#
#  Predicates and basic invariants
#
################################################################################

doc"""
***
    is_finite(A::FinGenGrpAbSnf) -> Bool

> Returns whether $A$ is finite.
"""
is_finite(A::FinGenGrpAbSnf) = A.snf[end] != 0
is_finite(A::FinGenGrpAbGen) = is_finite(snf(A)[1])

doc"""
***
    is_infinite(A::FinGenGrpAbSnf) -> Bool
    is_infinite(A::FinGenGrpAbGen) -> Bool

> Returns whether $A$ is infinite.
"""
is_infinite(A::FinGenGrpAbSnf) = !is_finite(A)
is_infinite(A::FinGenGrpAbGen) = !is_finite(A)

doc"""
***
    rank(A::FinGenGrpAbSnf) -> Int
    rank(A::FinGenGrpAbGen) -> Int

> Returns the rank of $A$, that is, the dimension of the
> $\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
"""
rank(A::FinGenGrpAbSnf) = length(find(x->x==0, A.snf))
rank(A::FinGenGrpAbGen) = rank(snf(A)[1])

doc"""
***
    order(A::FinGenGrpAbElem, i::Int) -> fmpz

> Returns the order of $A$. It is assumed that the order is finite.
"""
function order(a::FinGenGrpAbElem)
  G, m = snf(a.parent)
  b = m(a)
  o = fmpz(1)
  for i=1:ngens(G)
    if G.snf[i] == 0 && b[i] != 0
      error("element has inifinite order")
    end
    o = lcm(o, divexact(G.snf[i], gcd(G.snf[i], b[i])))
  end
  return o
end

doc"""
***
    order(A::FinGenGrpAbSnf) -> fmpz
    order(A::FinGenGrpAbGen) -> fmpz

> Returns the order of $A$. It is assumed that $A$ is not infinite.
"""
function order(A::FinGenGrpAbSnf)
  is_infinite(A) && error("Group must be finite")
  return prod(A.snf)
end
order(A::FinGenGrpAbGen) = order(snf(A)[1])

doc"""
***
    exponent(A::FinGenGrpAbSnf) -> fmpz
    exponent(A::FinGenGrpAbGen) -> fmpz

> Returns the exponent of $A$. It is assumed that $A$ is not infinite.
"""
function exponent(A::FinGenGrpAbSnf)
  is_infinite(A) && error("Group must be finite")
  return A.snf[end]
end
exponent(A::FinGenGrpAbGen) = exponent(snf(A)[1])

doc"""
***
    is_trivial(A::FinGenGrpAbSnf, i::Int) -> bool
    is_trivial(A::FinGenGrpAbGen, i::Int) -> bool

> Checks if A is the trivial group.
"""
is_trivial(A::FinGenGrpAbSnf) = order(A)==1
is_trivial(A::FinGenGrpAbGen) = order(A)==1

##############################################################################
#
#  Creation
#
##############################################################################

doc"""
***
  AbelianGroup(M::fmpz_mat) -> FinGenGrpAb
  AbelianGroup(M::Array{fmpz, 2}) -> FinGenGrpAb
  AbelianGroup(M::Array{Integer, 2}) -> FinGenGrpAb
  AbelianGroup(M::Array{fmpz, 1}) -> FinGenGrpAb
  AbelianGroup(M::Array{Integer, 1}) -> FinGenGrpAb

> Creates the abelian group with M as relation matrix. That is, the group will
> have cols(M) generators and each row describes one relation.
"""
function AbelianGroup(M::fmpz_mat)
  return FinGenGrpAbGen(M)
end

function AbelianGroup(M::Array{fmpz, 2})
  return AbelianGroup(MatrixSpace(FlintZZ, size(M)[1], size(M)[2])(M))
end

function AbelianGroup{T <: Integer}(M::Array{T, 2})
  return AbelianGroup(MatrixSpace(FlintZZ, size(M)[1], size(M)[2])(M))
end

function AbelianGroup(M::Array{fmpz, 1})
  return AbelianGroup(MatrixSpace(FlintZZ, 1, length(M))(M'))
end

function AbelianGroup{T <: Integer}(M::Array{T, 1})
  return AbelianGroup(MatrixSpace(FlintZZ, 1, length(M))(M'))
end

doc"""
***
  DiagonalGroup(M::fmpz_mat) -> FinGenGrpAb
  DiagonalGroup(M::Array{fmpz, 1}) -> FinGenGrpAb
  DiagonalGroup(M::Array{Integer, 1}) -> FinGenGrpAb

> Creates the abelian group as the product of cyclic groups Z/M[i].
"""
function DiagonalGroup(M::fmpz_mat)
  if rows(M) != 1
    error("The argument must have only one row")
  end

  N = MatrixSpace(FlintZZ, cols(M), cols(M))()
  for i=1:cols(M)
    N[i,i] = M[1,i]
  end
  if is_snf(N)
    return FinGenGrpAbSnf(fmpz[M[1,i] for i=1:cols(M)])
  else
    return FinGenGrpAbGen(N)
  end
end

function DiagonalGroup{T <: Union{Integer, fmpz}}(M::Array{T, 1})
  N = MatrixSpace(FlintZZ, length(M), length(M))()
  for i=1:length(M)
    N[i,i] = M[i]
  end
  if is_snf(N)
    return FinGenGrpAbSnf(M)
  else
    return FinGenGrpAbGen(N)
  end
end

# missing:
# is_torsion, torsion_subgroup, subgroups, ...

##############################################################################
#
##############################################################################

doc"""
***
  rand(G::FinGenGrpAbSnf) -> FinGenGrpAbElem
  rand(G::FinGenGrpAbGen) -> FinGenGrpAbElem

> For a finite abelian group G, return an element chosen uniformly at random.
"""
function rand(G::FinGenGrpAbSnf)
  if !is_finite(G)
    error("Group is not finite")
  end
  return G([rand(1:G.snf[i]) for i in 1:ngens(G)])
end

function rand(G::FinGenGrpAbGen)
  S, mS = snf(G)
  return preimage(mS, rand(S))
end

doc"""
***
  rand(G::FinGenGrpAbSnf, B::fmpz) -> FinGenGrpAbElem
  rand(G::FinGenGrpAbSnf, B::Integer) -> FinGenGrpAbElem
  rand(G::FinGenGrpAbGen, B::fmpz) -> FinGenGrpAbElem
  rand(G::FinGenGrpAbGen, B::Integer) -> FinGenGrpAbElem

> For a (potentially infinite) abelian group G, return an element
> chosen uniformly at random with coefficients bounded by B.
"""
function rand(G::FinGenGrpAbSnf, B::fmpz)
  return G([rand(1:(G.snf[i]==0 ? B : min(B, G.snf[i]))) for i in 1:ngens(G)])
end

function rand(G::FinGenGrpAbSnf, B::Integer)
  return rand(G, fmpz(B))
end

function rand(G::FinGenGrpAbGen, B::fmpz)
  S, mS = snf(G)
  return preimage(mS, rand(S, fmpz(B)))
end

function rand(G::FinGenGrpAbGen, B::Integer)
  return rand(G, fmpz(B))
end

##############################################################################
#
##############################################################################

doc"""
***
  sub(G::FinGenGrpAb, s::Array{FinGenGrpAbElem, 1}) -> FinGenGrpAb, Map

> Create the subgroup of G generated by the elements in s.
"""
function sub{T <: FinGenGrpAb}(G::T, s::Array{FinGenGrpAbElem{T}, 1})
  p = s[1].parent
  @assert G == p
  m = MatrixSpace(FlintZZ, length(s)+nrels(p), ngens(p)+length(s))()
  for i=1:length(s)
    for j=1:ngens(p)
      m[i + nrels(p),j] = s[i][j]
    end
    m[i + nrels(p), i+ngens(p)] = 1
  end
  if typeof(p) == FinGenGrpAbSnf
    for i=1:ngens(p)
      m[i, i] = p.snf[i]
    end
  else
    for i=1:nrels(p)
      for j=1:ngens(p)
        m[i, j] = p.rels[i,j]
      end
    end
  end
  h = hnf(m)
  r = sub(h, nrels(p)+1:rows(h), ngens(p)+1:cols(h))
  S = AbelianGroup(r)
  return S, FinGenGrpAbMap(S, p, sub(m, nrels(p)+1:rows(h), 1:ngens(p)))
end

doc"""
***
  quo(G::FinGenGrpAb, s::Array{FinGenGrpAbElem, 1}) -> FinGenGrpAb, Map

> Create the quotient of G by the elements in s.
"""
function quo{T <: FinGenGrpAb}(G::FinGenGrpAb, s::Array{FinGenGrpAbElem{T}, 1})
  p = s[1].parent
  @assert G == p
  m = MatrixSpace(FlintZZ, length(s)+nrels(p), ngens(p))()
  for i=1:length(s)
    for j=1:ngens(p)
      m[i + nrels(p),j] = s[i][j]
    end
  end
  if typeof(p) == FinGenGrpAbSnf
    for i=1:ngens(p)
      m[i, i] = p.snf[i]
    end
  else
    for i=1:nrels(p)
      for j=1:ngens(p)
        m[i, j] = p.rels[i,j]
      end
    end
  end

  Q = AbelianGroup(m)
  I = MatrixSpace(FlintZZ, ngens(p), ngens(p))(1)
  m = FinGenGrpAbMap(p, Q, I, I)
  return Q, m
end

doc"""
***
  quo(G::FinGenGrpAb, n::Integer}) -> FinGenGrpAb, Map
  quo(G::FinGenGrpAb, n::fmpz}) -> FinGenGrpAb, Map

> Create the quotient of G by n*G
"""
function quo(G::FinGenGrpAbSnf, n::Union{fmpz, Integer})
  r = [gcd(x, n) for x = G.snf]
  I = MatrixSpace(FlintZZ, ngens(G), ngens(G))(1)
  Q = DiagonalGroup(r)
  m = FinGenGrpAbMap(G, Q, I, I)
  return Q, m
end

function quo(G::FinGenGrpAbGen, n::Union{fmpz, Integer})
  m = vcat(G.rels, MatrixSpace(FlintZZ, ngens(G), ngens(G))(n))
  Q = AbelianGroup(m)
  I = MatrixSpace(FlintZZ, ngens(G), ngens(G))(1)
  m = FinGenGrpAbMap(G, Q, I, I)
  return Q, m
end

#=  <example>
M = MatrixSpace(ZZ, 2, 3)([1 2 3; 4 5 6])
G = Hecke.FinGenGrpAbGen(M)
a = Hecke.FinGenGrpAbElemCreate(G, MatrixSpace(ZZ, 1, 3)([1 2 3]))
100*a
200*a
300*a
H, mp = snf(G)
mp(a)

G = Hecke.FinGenGrpAbSnf(fmpz[2,3,4])
=#


text:
using Printf
using LinearAlgebra
using CSV
using StatsBase
import Base.*, Base.+

dn_source = dirname(Base.source_path())
push!(LOAD_PATH, dn_source)
using DL

# 関数の合成
∘(f::Function, g::Function) = x::Tensor-> f(g(x::Tensor))

# Arrayをシャッフルする関数
function shuffle(vx)
    n = length(vx)
    return sample(vx, n, replace=false)
end

function main()
    filename = "iris.csv"
    data = CSV.read(filename)
    n_data, n_type = size(data)
    # n_train = 150 * 0.8 = 120
    n_train::Int64 = n_data * 0.8
    # n_test = 150 - 120  = 30
    n_test::Int64 = n_data - n_train

    # 花の名前 -> Onehot
    # irisデータセットでは3種類の花の名前がある
    花の名前s = unique(data[end])
    # n_hana = 3
    n_hana = length(花の名前s)
    vid2onehot = function(x::Int64)
        onehot = zeros(n_hana, 1)
        onehot[x, 1] = 1.0
        return onehot
    end
    # 花の名前をonehot形式に変換する
    花toOnehot = x::String -> vid2onehot(findfirst(花の名前s .== x))
    # argmax(onehot)を花の名前に変える
    id2name = x -> 花の名前s[x]

    # 入力の次元
    in_features = n_type - 1
    # 出力の次元
    out_features = n_hana

    data = map(i -> [data[i, :]...], 1:n_data)
    # データをshuffle
    data = shuffle(data)
    # 学習データ
    vtrain_x = map(i -> Float64.(data[i][1:in_features]), 1:n_train)
    vtrain_y = map(i -> data[i][n_type], 1:n_train)
    vtrain_y = map(花toOnehot, vtrain_y)
    traindataset = collect(zip(vtrain_x, vtrain_y))
    # テストデータ
    vtest_x = map(i -> Float64.(data[i][1:in_features]), n_train+1:n_data)
    vtest_y = map(i -> data[i][n_type], n_train+1:n_data)
    vtest_y = map(花toOnehot, vtest_y)
    testdataset = collect(zip(vtest_x, vtest_y))

    # 層
    W1 = VariableTensor(randn(32, in_features)/sqrt(32*in_features))
    b1 = VariableTensor(zeros(32, 1))
    W2 = VariableTensor(randn(out_features, 32)/sqrt(out_features*32))
    b2 = VariableTensor(zeros(out_features, 1))
    l1 = x -> W1 * x + b1
    l2 = x -> W2 * x + b2
    # モデル
    model = l2 ∘ ReLU ∘ l1
    # 損失関数
    loss_func = (predict::Tensor, data::Tensor) -> MSE(predict, data)
    # 最適化関数
    optimizer = SGD(1e-4)

    f = open("loss.dat", "w")
    for epoch in 1:1000
        # train
        @show epoch
        train_loss::Float64 = 0
        # traindatasetをシャッフルする
        traindataset = shuffle(traindataset)
        for (batch_idx, (x, y)) in enumerate(traindataset)
            x = ConstantTensor(x)
            y = ConstantTensor(y)
            predict = model(x)
            loss = loss_func(predict, y)
            train_loss += loss.x
            # W1.grad, W2.grad, b1.grad, b2.gradに微分を足しこんでいく
            backward(loss)

            # .gradを使ってW1.x, W2.x, b1.x, b2.xを更新する
            # batch_size = 30で学習する.
            if batch_idx % 30 == 0 || batch_idx == n_train
                optimizer(loss)
                # .gradを0にする
                zero_grad(loss)
            end
        end

        @show train_loss

        # test
        正解数::Int64 = 0
        for (x, y) in testdataset
            # 入力
            x = ConstantTensor(x)
            # 教師データ
            y = ConstantTensor(y)
            p = model(x)
            pname = id2name(argmax(p.x))
            name = id2name(argmax(y.x))
            正解数 += pname == name
        end
        @printf(f, "%d, % .3e, %.3f\n", epoch, train_loss, 正解数 / n_test)

    end
    close(f)
end

main()

max_stars_repo_name: Soyukke/julia-ayame-machine-learning

max_issues_repo_name: Soyukke/julia-ayame-machine-learning

max_forks_repo_name: Soyukke/julia-ayame-machine-learning


text:
using Documenter, MethodOfLines

makedocs(
    sitename="MethodOfLines.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[MethodOfLines],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://methodoflines.sciml.ai/stable/"),

    pages=[
        "MethodOfLines.jl: Automated Finite Difference for Phyiscs-Informed Learning" => "index.md",
     ]
)

deploydocs(
    repo="github.com/SciML/MethodOfLines.jl";
    push_preview=true
)


text:
module TestDual

using Compat
using Compat.Test
using IntervalArithmetic
using EAGODomainReduction

Y = Interval{Float64}[Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0)]
ymult_lo = [50, 1.0, 2.0, 3.0]
ymult_hi = [0, 1.0, 2.0, 3.0]
yLBD = 1.0
yUBD = 3.0
Variable_DR!(Y,ymult_lo,ymult_hi,yLBD,yUBD)
println("Y post full duality-based bound tighteing: ", Y)
@test 3.95999-1E-4 <= Y[1].lo <= 3.95999+1E-4
@test 4.0-1E-4 <= Y[1].hi <= 4.0+1E-4
@test 2.0-1E-4 <= Y[2].lo <= 2.0+1E-4
@test 4.0-1E-4 <= Y[2].hi <= 4.0+1E-4
@test 3.0-1E-4 <= Y[3].lo <= 3.0+1E-4
@test 4.0-1E-4 <= Y[3].hi <= 4.0+1E-4
@test 3.33333-1E-4 <= Y[4].lo <= 3.33333+1E-4
@test 4.0-1E-4 <= Y[4].hi <= 4.0+1E-4

X = Interval{Float64}[Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0),Interval(1.0,4.0)]
mult_lo = [50, 1.0]
mult_hi = [-1.0, 0]
LBD = 1.0
UBD = 2.0
nx = 2
Variable_DR_Imp!(X,mult_lo,mult_hi,LBD,UBD,nx)
println("X post implicit duality-based bound tighteing: ", X)
@test 1.0-1E-4 <= X[1].lo <= 1.0+1E-4
@test 4.0-1E-4 <= X[1].hi <= 4.0+1E-4
@test 1.0-1E-4 <= X[2].lo <= 1.0+1E-4
@test 4.0-1E-4 <= X[2].hi <= 4.0+1E-4
@test 3.97999-1E-4 <= X[3].lo <= 3.97999+1E-4
@test 4.0-1E-4 <= X[3].hi <= 4.0+1E-4
@test 3.0-1E-4 <= X[4].lo <= 3.0+1E-4
@test 4.0-1E-4 <= X[4].hi <= 4.0+1E-4

end


text:
using YeSQL
using Documenter

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

makedocs(;
    modules=[YeSQL],
    authors="Todd Hay",
    repo="https://github.com/tah-alv/YeSQL.jl/blob/{commit}{path}#{line}",
    sitename="YeSQL.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tah-alv.github.io/YeSQL.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "pages/api.md",
    ],
)

deploydocs(;
    repo="github.com/tah-alv/YeSQL.jl",
)