text:
using Documenter
using PrettyTables

makedocs(
    modules = [PrettyTables],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://ronisbr.github.io/PrettyTables.jl/stable/",
    ),
    sitename = "Pretty Tables",
    authors = "Ronan Arraes Jardim Chagas",
    pages = [
        "Home"               => "index.md",
        "Usage"              => "man/usage.md",
        "Back-ends"          => Any[
            "Text"           => "man/text_backend.md",
            "HTML"           => "man/html_backend.md",
        ],
        "Alignment"          => "man/alignment.md",
        "Filters"            => "man/filters.md",
        "Formatters"         => "man/formatter.md",
        "Examples"           => Any[
            "Text back-end"  => "man/text_examples.md",
            "HTML back-end"  => "man/html_examples.md",
        ],
    ],
)

deploydocs(
    repo = "github.com/ronisbr/PrettyTables.jl.git",
    target = "build",
)


text:
using ODE, DiffEqBase, DiffEqProblemLibrary, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
using LinearAlgebra

dt=1/2^(4)

algs = [ode23(),ode45(),ode78(),ode4(),ode4ms(),ode4s()] # no ode23s

# Check for errors

prob = prob_ode_linear

for alg in algs
    sol =solve(prob,alg;dt=dt,abstol=1e-6,reltol=1e-3)
    @test typeof(sol[2]) <: Number
end

prob = prob_ode_2Dlinear

for alg in algs
    sol =solve(prob,alg;dt=dt,dtmin=eps(),abstol=1e-6,reltol=1e-3)
    @test size(sol[2]) == (4,2)
end

# Check that warnings are displayed
for alg in algs
    sol = solve(prob,alg;dt=dt,abstol=1e-6,reltol=1e-3, dense=true, verbose=true,
                save_idxs = true, progress = true, beta1 = 1.23, gamma = 0.5)
end

ext: jl


text:
using MambaModels
using Test

println("Tests are executed in docs job.")

@test 1 == 1


text:
using NPZ
using JuMP
using Ipopt
using LinearAlgebra
using .Threads
using WoodburyMatrices
using Statistics
using UnicodePlots

function main(B)
    X1 = npzread("Kernels/G1.npy")
    X2 = npzread("Kernels/G2.npy")
    X3 = npzread("Kernels/G3.npy")
    nfactor = npzread("Testing/nfactor.npy")

    n = size(X3, 1)
    ndas = size(X1, 1)

    T = Float64
    ϵ = 1.0
    λ = 100.0
    α = 100.0

    X = [X3 zeros(T, size(X3));
        zeros(T, size(X3)) X3;
        nfactor*X1 nfactor*X2]
    Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    #Xn = sqrt.(sum(X.^2, dims=2)) #normalize all rows to be 1?
    X ./= Xn

    v, V = eigen(X'*X + ϵ*I)

    #eigtol = 5#quantile(v, 0.99)

    n_eigs = 5#sum(v.>=eigtol)

    Γ = V[:, (end-n_eigs+1):end]

    μ_init = ones(n) ./ n .* B
    μn = sum(μ_init)/B
    μ_init ./= μn
    U_init = copy(Γ)

    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        return (M, Mc)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function subproblem_sμ(μ_init, u, X, λ, B, maxiter=1000)

        function f(μ, u, λ)
            M, Mc = Ml(μ)
            return -logdet(M) + λ*tr(u'*(Mc\u))
        end

        f(μ) = f(μ, u, λ)

        jump_wrapper_f(x...) = f(vcat([x[i] for i in 1:length(x)], [x[i] for i in 1:length(x)], ones(T, ndas)))

        function gi(i, Mc, μ, u, λ)
            ld_term = @views -X[i,:]'*(Mc\X[i,:])
            tmp_A = SymWoodbury(Mc, X[i,:], -μ[i])
            tmp_b = @views tmp_A \ X[i,:]
            #tmp_1 = (u'*tmp_b)^2
            tmp_1 = sum((u'*tmp_b).^2) # think that this should be correct for tr(U'MU) for matrix U but need to check
            tmp_d = @views X[i,:]'*tmp_b
            tmp_den = (1+μ[i]*tmp_d)^2
            eig_term = ((1+μ[i]*tmp_d)*tmp_1 - μ[i]*tmp_1*tmp_d) / tmp_den
            return ld_term - λ*eig_term #apparently need a -ve sign here, work out why
        end

        function g!(gvec, μ, X, u, λ, n)
            M, Mc = Ml(μ)
            @threads for i = 1:n
                g1 = gi(i, Mc, μ, u, λ) # in respect to x component
                g2 = gi(i+n, Mc, μ, u, λ) # in respect to y component
                gvec[i] = g1 + g2
            end
        end

        g!(gvec, μ) = g!(gvec, μ, X, u, λ, n)

        jump_wrapper_g!(gvec, x...) = g!(gvec, vcat([x[i] for i in 1:length(x)], [x[i] for i in 1:length(x)], ones(T, ndas)))

        model = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter"=>maxiter, "tol"=>1e-4/B))
        JuMP.register(model, :f, n, jump_wrapper_f, jump_wrapper_g!)

        @variable(model, 0.0<=x[1:n]<= 1.0)

        set_start_value.(x, μ_init)

        @NLconstraint(model, +(x...) <= B)

        @NLobjective(model, Min, f(x...))

        optimize!(model)

        return (value.(x), vcat(value.(x), value.(x), ones(T, ndas)), objective_value(model))

    end

    function subproblem_sO(u_init, μ, γ, X, λ, α, maxiter=1000)
        _, Mc = Ml(μ)

        f(u) = α*sum((u.-γ).^2) + λ*tr(u'*(Mc\u))
        function g!(gvec, u)
            gvec .= 2*α*(u-γ) + 2*λ*(Mc\u)
        end
        jump_wrapper_f(x...) = f([x[i] for i in 1:length(x)])
        jump_wrapper_g!(gvec, x...) = g!(gvec, [x[i] for i in 1:length(x)])

        sumsquares(x...) = sum(x[i]^2 for i in 1:length(x))
        function sumsquares_grad!(gvec, x...)
            for i in 1:length(x)
                gvec[i] = 2*x[i]
            end
        end

        model = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter"=>maxiter))
        JuMP.register(model, :f, length(γ), jump_wrapper_f, jump_wrapper_g!)
        JuMP.register(model, :sumsquares, length(γ), sumsquares, sumsquares_grad!)

        @variable(model, x[1:length(γ)])

        set_start_value.(x, u_init)

        @NLconstraint(model, sumsquares(x...) == 1)

        @NLobjective(model, Min, f(x...))

        optimize!(model)

        return (value.(x), objective_value(model))

    end

    function subproblem_sO(U_init, μ, Γ, X, λ, α, m, n, maxiter=1000)
            _, Mc = Ml(μ)

            f(U) = α*sum((U.-Γ).^2) + λ*tr(U'*(Mc\U))

            function g!(gvec, U, m, n)
                @threads for i = 1:n
                    @views gvec[((i-1)*m+1):i*m] .= 2*α*(U[1:m, i].-Γ[1:m, i]) .+ 2*λ*(Mc\U[1:m, i])
                end
            end

            L = LinearIndices((m, n))
            jump_wrapper_f(x...) = f([x[L[i,j]] for i in 1:m, j in 1:n])
            jump_wrapper_g!(gvec, x...) = g!(gvec, [x[L[i,j]] for i in 1:m, j in 1:n], m, n)

            sumprod(x...) = sum(x[i]*x[i+m] for i in 1:m)
            function sumprod_grad!(gvec, x...)
                for i in 1:m
                    gvec[i] = x[i+m]
                    gvec[i+m] = x[i]
                end
            end

            model = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter"=>maxiter))
            JuMP.register(model, :f, m*n, jump_wrapper_f, jump_wrapper_g!)
            JuMP.register(model, :sumprod, 2*m, sumprod, sumprod_grad!)

            @variable(model, x[1:m, 1:n])

            set_start_value.(x, U_init)

            for i = 1:n
                for j = i:n
                    x_i = x[1:m, i]
                    x_j = x[1:m, j]
                    @NLconstraint(model, sumprod(x_i..., x_j...) == ((i==j) ? 1 : 0))            end
            end

            @NLobjective(model, Min, f(x...))

            optimize!(model)

            return (value.(x), objective_value(model)) # want to have all of the x here as a matrix
    end

    function optimize_spectral(μ_init, u_init, iters, γ, X, λ, α, B)
        μ = copy(μ_init)
        u = copy(u_init)
        ob1l = Float64[]
        ob2l = Float64[]
        for iter in 1:iters
            μ, μl, ob1 = subproblem_sμ(μ, u, X, λ, B)
            u, ob2 = subproblem_sO(u, μl, γ, X, λ, α)
            push!(ob1l, ob1)
            push!(ob2l, ob2)
        end
        μ, μl, ob1 = subproblem_sμ(μ, u, X, λ, B)
        push!(ob1l, ob1)

        return (μ, u, ob1l, ob2l)
    end

    function optimize_spectral(μ_init, U_init, iters, Γ, X, λ, α, B, m, n)
        μ = copy(μ_init)
        U = copy(U_init)
        ob1l = Float64[]
        ob2l = Float64[]
        for iter in 1:iters
            μ, μl, ob1 = subproblem_sμ(μ, U, X, λ, B)
            U, ob2 = subproblem_sO(U, μl, Γ, X, λ, α, m, n)
            push!(ob1l, ob1)
            push!(ob2l, ob2)
        end
        μ, μl, ob1 = subproblem_sμ(μ, U, X, λ, B)
        push!(ob1l, ob1)

        return (μ, U, ob1l, ob2l)
    end

    μ, U_res, ob1l, ob2l = optimize_spectral(μ_init, U_init, 5, Γ, X, λ, α, B, size(Γ)...)

    M, _ = Ml(vcat(μ,μ, ones(ndas)))
    vμ, Vμ = eigen(M)

    npzwrite("SpectralResults/res_spectral_$B.npy", μ)

    lineplot(ob1l)

    lineplot(ob2l)
end

for i in [4,8,16,32,64,128,256]
    main(i)
end

#yeh

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
include("../examples/baker.jl")
using PyPlot
function plot_x_and_y()
    Npts = 10000
    s = [0.,0.,0.,0]
	eps = 1.e-4
	u = Array(LinRange(eps, π - eps, Npts))
	upπ = Array(LinRange(π + eps, 2π - eps, Npts))
	x = collect([u upπ]')
	x2 = collect([upπ upπ]')
	x3 = collect([u u]')
	x1 = collect([upπ u]')
    for i = 1:Npts
	    x[:,i] = next(x[:,i], s)
	    x1[:,i] = next(x1[:,i], s)
	    x2[:,i] = next(x2[:,i], s)
	    x3[:,i] = next(x3[:,i], s)
    end
	fig, ax = subplots()
	ax.plot(u, x[1,:],"r.",ms=1.9)
	ax.plot(upπ, x1[1,:],"r.",ms=1.9)
    ax.xaxis.set_tick_params(labelsize=28)
    ax.yaxis.set_tick_params(labelsize=28)
	ax.grid(true)
	ax.set_xlabel(L"$U_0$", fontsize=28, labelpad=-20)
	ax.set_ylabel(L"$U_1$", fontsize=28)

    ax.axis("scaled")

	fig, ax = subplots()
	ax.plot(upπ, x[2,:],"k--",lw=2.9)
	ax.plot(u, x1[2,:],"b",lw=2.9,label=L"$U_0 > \pi$")
	ax.plot(upπ, x2[2,:],"b",lw=2.9)
	ax.plot(u, x3[2,:],"k--",lw=2.9, label=L"$U_0 < \pi$")
    ax.xaxis.set_tick_params(labelsize=28)
    ax.yaxis.set_tick_params(labelsize=28)
	ax.grid(true)
    ax.legend(fontsize=28)
	ax.set_xlabel(L"$S_0$", fontsize=28, labelpad=-20)
	ax.set_ylabel(L"$S_1$", fontsize=28)

    ax.axis("scaled")
end
function plot_one_step()
    n_xq, n_yq = 20, 40
    eps = 1.e-14
    x_q = LinRange(0, pi - eps, n_xq)
    y_q = LinRange(0, 2*pi, n_yq)
    #x_q = x_q[1:end-1]
    #y_q = y_q[1:end-1]
    #n_xq = n_xq - 1
    #n_yq = n_yq - 1
    cm = plt.get_cmap("winter")
    x_q = repeat(x_q, 1, n_yq)
    y_q = repeat(y_q', n_xq, 1)
    fig, ax = subplots(1,1)
    fig1, ax1 = subplots(1,1)

    cplot = ax.contour(x_q, y_q, y_q, n_yq,
                       cmap=cm)
    cbar = fig.colorbar(cplot, ax=ax)
    x1_q = zeros(n_xq, n_yq)
    y1_q = zeros(n_xq, n_yq)

    for j = 1:n_yq, i = 1:n_xq
        x1_q[i,j], y1_q[i,j] =
            step([x_q[i,j], y_q[i,j]]
                 , zeros(4), 1)[:,end]
    end
    cplot1 = ax1.contour(x1_q, y1_q, y_q, n_yq,
                       cmap=cm)
    x_q = x_q'
    y_q = y_q'
    x1_q = x1_q'
    y1_q = y1_q'
    ax.contour(x_q, y_q, x_q, n_xq,
                       cmap=cm, vmin=0.,vmax=2*pi)
    ax1.contour(x1_q, y1_q, x_q, n_xq,
                       cmap=cm, vmin=0., vmax=2*pi)

    #n_xq = n_xq + 1
    #n_yq = n_yq + 1
    x_q = LinRange(pi + eps, 2*pi, n_xq)
    #x_q = x_q[1:end-1]
    y_q = LinRange(0., 2*pi, n_yq)
    #y_q = y_q[1:end-1]
    #n_xq = n_xq - 1
    #n_yq = n_yq - 1
    x_q = repeat(x_q, 1, n_yq)
    y_q = repeat(y_q', n_xq, 1)

    ax.contour(x_q, y_q, y_q, n_yq,
                cmap=cm)
    x1_q = zeros(n_xq, n_yq)
    y1_q = zeros(n_xq, n_yq)

    for j = 1:n_yq, i = 1:n_xq
        x1_q[i,j], y1_q[i,j] =
            step([x_q[i,j], y_q[i,j]]
                 , zeros(4), 1)[:,end]
    end
    cplot1 = ax1.contour(x1_q, y1_q, y_q, n_yq,
                       cmap=cm)
    cbar1 = fig1.colorbar(cplot, ax=ax1)
    x_q = x_q'
    y_q = y_q'
    x1_q = x1_q'
    y1_q = y1_q'
    ax.contour(x_q, y_q, x_q, n_xq,
                       cmap=cm,vmin=0.,vmax=2*pi)
    ax1.contour(x1_q, y1_q, x_q, n_xq,
                       cmap=cm, vmin=0., vmax=2*pi)

    ax1.xaxis.set_tick_params(labelsize=28)
    ax1.yaxis.set_tick_params(labelsize=28)

    ax.xaxis.set_tick_params(labelsize=28)
    ax.yaxis.set_tick_params(labelsize=28)
    cbar.ax.tick_params(labelsize=28)
    cbar1.ax.tick_params(labelsize=28)
    ax.axis("scaled")
    ax1.axis("scaled")
    ax.set_xlim([0.,2*pi])
    ax.set_ylim([0.,2*pi])
    ax1.set_xlim([0.,2*pi])
    ax1.set_ylim([0.,2*pi])
end
function plot_mapping_vs_params()
    n_xpts, n_ypts = 75, 75
    x_pts = LinRange(0, 2*pi, n_xpts)
    #x_pts = x_pts[1:end-1]
    #n_xpts = n_xpts - 1

    y_pts = LinRange(0, 2*pi, n_ypts)
    #y_pts = y_pts[1:end-1]
    #n_ypts = n_ypts - 1

    cm = plt.get_cmap("winter")
    x_gr = repeat(x_pts, 1, n_ypts)
    y_gr = repeat(y_pts', n_xpts, 1)
    x_gr = x_gr[:]
    y_gr = y_gr[:]
    n_gr = length(x_gr)

    #n_xpts = n_xpts + 1
    #n_ypts = n_ypts + 1
    eps = 1.e-1
    x_pts = LinRange(eps, 2*pi, n_xpts)
    #x_pts = x_pts[1:end-1]
    #n_xpts = n_xpts - 1

    y_pts = LinRange(eps, 2*pi, n_ypts)
    #y_pts = y_pts[1:end-1]
    #n_ypts = n_ypts - 1

    x_gr1 = repeat(x_pts, 1, n_ypts)
    y_gr1 = repeat(y_pts', n_xpts, 1)

    x_gr1 = x_gr1[:]
    y_gr1 = y_gr1[:]

    for n = 1:4
        s = zeros(4)
        s[n] = 1.0

        x1_gr = zeros(n_gr)
        y1_gr = zeros(n_gr)

        for i = 1:n_gr
            x1_gr[i], y1_gr[i] =
                step([x_gr[i], y_gr[i]]
                     , s, 1)[:,end]
        end
        # Horizontal lines
        fig1, ax1 = subplots(1,1)
        clrs = cm(y_gr/(2*pi))
        sp1 = ax1.scatter(x1_gr, y1_gr, s=40.0,
                          c=clrs, marker="_",
                             cmap=cm)
        #cbar = fig.colorbar(sp1, ax=ax1)

        # Vertical lines
        x1_gr1 = zeros(n_gr)
        y1_gr1 = zeros(n_gr)

        for i = 1:n_gr
            x1_gr1[i], y1_gr1[i] =
                step([x1_gr1[i], y1_gr1[i]]
                     , s, 1)[:,end]
        end

        clrs = cm(x_gr1/(2*pi))
        ax1.scatter(x_gr1, y_gr1, s=40.0,
                    marker="|",c=clrs, cmap=cm)
        ax1.xaxis.set_tick_params(labelsize=28)
        ax1.yaxis.set_tick_params(labelsize=28)
        ax1.axis("scaled")
        ax1.set_xlim([0.,2*pi])
        ax1.set_ylim([0.,2*pi])
        ax1.set_title("\$s_$n = 1\$", fontsize=28)

    end
end
function plot_annotated_map_rectangles(s, n)
    N = 500000
    n_gr_x = 2
    n_gr_y = 8
    dx = 2π/n_gr_x
    dy = 2π/n_gr_y
    fig = figure(figsize=(8,6))
    ax = fig.add_subplot(111)
    clrs = matplotlib.cm["tab20c"]
    fig1 = figure(figsize=(8,6))
    ax1 = fig1.add_subplot(111)
    ax.set_xlim([0.,2π])
    ax.set_ylim([0.,2π])
    ax1.set_xlim([0.,2π])
    ax1.set_ylim([0.,2π])

    for j = 1:n_gr_x
        x = (j-1)*dx
        for k = 1:n_gr_y
            y = (k-1)*dy
            u1 = [[x + rand()*dx, y + rand()*dy] for i=1:N]
            u1_next = copy(u1)
            for m = 1:n
                u1_next = next.(u1_next, Ref(s))
            end
            u1_next = hcat(u1_next...)
            u1 = hcat(u1...)

            ind = k + (j-1)*n_gr_y
            if (ind > 20)
                    ind = mod(ind, 20)
            end
            ax1.plot(u1_next[1,:], u1_next[2,:],".",color=clrs.colors[ind],ms=1.)
            ax.plot(u1[1,:], u1[2,:],".", color=clrs.colors[ind], ms=1.0)
        end
    end
    ax.xaxis.set_tick_params(labelsize=28)
    ax1.xaxis.set_tick_params(labelsize=28)
    ax.yaxis.set_tick_params(labelsize=28)
    ax1.yaxis.set_tick_params(labelsize=28)
    ax1.set_title("\$ \\varphi^{$n} \$",fontsize=32)
    ax.axis("scaled")
    ax1.axis("scaled")

end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
############################################
#
#    mshr
#
# https://bitbucket.org/dolfin-project/mshr/wiki/Home
#
############################################

# no need to export
# not exported in `mshr`
@pyclass mshr CSGGeometry

# 2D primitives
@pyclass mshr Rectangle CSGGeometry
@pyclass mshr Circle CSGGeometry
@pyclass mshr Polygon CSGGeometry
@pyclass mshr Ellipse CSGGeometry

# 3D primitives
@pyclass mshr Cylinder CSGGeometry
@pyclass mshr Box CSGGeometry
@pyclass mshr Surface3D CSGGeometry
@pyclass mshr Cone CSGGeometry
@pyclass mshr Ellipsoid CSGGeometry
@pyclass mshr Sphere CSGGeometry
@pyclass mshr Tetrahedron CSGGeometry

@pyfunc mshr generate_mesh

+(geo1::CSGGeometry, geo2::CSGGeometry) = CSGGeometry(geo1.pyobject + geo2.pyobject)
-(geo1::CSGGeometry, geo2::CSGGeometry) = CSGGeometry(geo1.pyobject - geo2.pyobject)

export Rectangle, Circle, Polygon, Ellipse,
       Cylinder, Box, Surface3D, Cone, Ellipsoid, Sphere, Tetrahedron,
       generate_mesh


text:
using Contractions
using Test

@testset "Contractions.jl" begin
    # Write your tests here.
    docx1 = "He'll be there, if not I can't"

    # Results
    results_for_fix_contractions = "He'll be there, if not I cannot"
    results_for_fix_contractions_lowercase = "he will be there, if not i cannot"
    @test Contractions.fix_contractions(docx1) == results_for_fix_contractions
    @test Contractions.fix_contractions(lowercase(docx1)) == results_for_fix_contractions_lowercase

end


text:
# ---
# title: 695. Max Area of Island
# id: problem695
# author: Indigo
# date: 2021-06-29
# difficulty: Medium
# categories: Array, Depth-first Search
# link: <https://leetcode.com/problems/max-area-of-island/description/>
# hidden: true
# ---
#
# Given a non-empty 2D array `grid` of 0's and 1's, an **island** is a group of
# `1`'s (representing land) connected 4-directionally (horizontal or vertical.)
# You may assume all four edges of the grid are surrounded by water.
#
# Find the maximum area of an island in the given 2D array. (If there is no
# island, the maximum area is 0.)
#
# **Example 1:**
#
#
#
#     [[0,0,1,0,0,0,0,1,0,0,0,0,0],
#      [0,0,0,0,0,0,0,1,1,1,0,0,0],
#      [0,1,1,0,1,0,0,0,0,0,0,0,0],
#      [0,1,0,0,1,1,0,0, **1** ,0, **1** ,0,0],
#      [0,1,0,0,1,1,0,0, **1** , **1** , **1** ,0,0],
#      [0,0,0,0,0,0,0,0,0,0, **1** ,0,0],
#      [0,0,0,0,0,0,0,1,1,1,0,0,0],
#      [0,0,0,0,0,0,0,1,1,0,0,0,0]]
#
#
# Given the above grid, return `6`. Note the answer is not 11, because the
# island must be connected 4-directionally.
#
# **Example 2:**
#
#
#
#     [[0,0,0,0,0,0,0,0]]
#
# Given the above grid, return `0`.
#
# **Note:** The length of each dimension in the given `grid` does not exceed 50.
#
#
## @lc code=start
using LeetCode

function max_area_of_island(grid::Matrix{Int})
    function dfs(grid::Matrix{Int}, I)
        grid[I] = 0
        res = 1
        idcs = CartesianIndices(grid)
        for neighb in ((1, 0), (0, 1), (-1, 0), (0, -1))
            new_I = I + CartesianIndex(neighb)
            (new_I ∈ idcs && grid[new_I] == 1) && (res += dfs(grid, new_I))
        end
        return res
    end
    res = 0
    for I in CartesianIndices(grid)
        (grid[I] == 1) && (res = max(res, dfs(grid, I)))
    end
    return res
end
## @lc code=end


text:
nn = S.NN([19,1], [S.Relu()])
ref = S.calc(nn, vec(vars[12345,:]))
ts = S.cycle(nn, vec(vars[12345,:]), [res[12456]])
nn.ws[1][1,1] += 1e-8
dref = ( (S.calc(nn, vec(vars[12345,:]))-[res[12456]]).^2 - (ref-[res[12456]]).^2)*1e8  # 1.536
[ ts.dws[1][1,1] dref ]

nn = S.NN([19,1], [S.SoftRelu()])
ref = S.calc(nn, vec(vars[12345,:]))
ts = S.cycle(nn, vec(vars[12345,:]), [res[12456]])
nn.ws[1][1,1] += 1e-8
dref = ( (S.calc(nn, vec(vars[12345,:]))-[res[12456]]).^2 - (ref-[res[12456]]).^2)*1e8  # 1.536
[ ts.dws[1][1,1] dref ]

nn = S.NN([19,1], [S.Sigmoid()])
ref = S.calc(nn, vec(vars[12345,:]))
ts = S.cycle(nn, vec(vars[12345,:]), [res[12456]])
nn.ws[1][1,1] += 1e-8
dref = ( (S.calc(nn, vec(vars[12345,:]))-[res[12456]]).^2 - (ref-[res[12456]]).^2)*1e8  # 1.536
[ ts.dws[1][1,1] dref ]  # faux

nn = S.NN([19,1], [S.Tanh()])
ref = S.calc(nn, vec(vars[12345,:]))
ts = S.cycle(nn, vec(vars[12345,:]), [res[12456]])
nn.ws[1][1,1] += 1e-8
dref = ( (S.calc(nn, vec(vars[12345,:]))-[res[12456]]).^2 - (ref-[res[12456]]).^2)*1e8  # 1.536
[ ts.dws[1][1,1] dref ]

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
module Diagnostics

    using GeometricBase

    export compute_one_form, compute_invariant, compute_invariant_error, compute_momentum_error, compute_error_drift

    """
    Takes a 1d DataSeries holding an invariant and computes the relative error `(inv(t)-inv(0))/inv(0)`.

    Returns a 1d DataSeries similar to the argument holding the time series of the relativ errors.
    """
    function compute_relative_error(invds::AbstractDataSeries{T,1}) where {T}
        errds = similar(invds)
        for i in eachindex(errds, invds)
            errds[i] = (invds[i] - invds[0]) / invds[0]
        end
        return errds
    end

    """
    Compute the one-form (symplectic potential) for the solution of a Lagrangian system.

    Arguments: `(t::TimeSeries, q::AbstractDataSeries, one_form::Function)`

    The `one_form` function needs to take three arguments `(t,q,k)` where `k` is the index
    of the one-form component.

    Returns a DataSeries similar to `q` holding the time series of the one-form.
    """
    function compute_one_form(t::TimeSeries, q::AbstractDataSeries, one_form::Function)
        ϑ = similar(q)
        try
            for i in axes(p,2)
                for k in axes(p,1)
                    ϑ[k,i] = one_form(t[i], q[:,i], k)
                end
            end
        catch ex
            if isa(ex, DomainError)
                @warn("DOMAIN ERROR: One-form diagnostics crashed.")
            else
                throw(ex)
            end
        end
        return ϑ
    end

    """
    Compute an invariant for the solution of an ODE or DAE system.

    Arguments: `(t::TimeSeries, q::AbstractDataSeries{T}, invariant::Function)`

    The `invariant` functions needs to take two arguments `(t,q)` and return the
    corresponding value of the invariant.

    Returns a 1d DataSeries holding the time series of the invariant.
    """
    function compute_invariant(t::TimeSeries, q::AbstractDataSeries{<:AbstractArray{T}}, invariant::Function) where {T}
        invds = DataSeries(T, ntime(q))
        try
            for i in eachindex(invds)
                invds[i] = invariant(t[i], q[i])
            end
        catch ex
            if isa(ex, DomainError)
                @warn("DOMAIN ERROR: Invariant diagnostics crashed.")
            else
                throw(ex)
            end
        end
        return invds
    end

    """
    Compute an invariant for the solution of a partitioned ODE or DAE system.

    Arguments: `(t::TimeSeries, q::AbstractDataSeries{T}, p::AbstractDataSeries{T}, invariant::Function)`

    The `invariant` functions needs to take three arguments `(t,q,p)` and return the
    corresponding value of the invariant.

    Returns a 1d DataSeries holding the time series of the invariant.
    """
    function compute_invariant(t::TimeSeries, q::AbstractDataSeries{<:AbstractArray{T}}, p::AbstractDataSeries{<:AbstractArray{T}}, invariant::Function) where {T}
        invds = DataSeries(T, ntime(q))
        try
            for i in eachindex(invds)
                invds[i] = invariant(t[i], q[i], p[i])
            end
        catch ex
            if isa(ex, DomainError)
                @warn("DOMAIN ERROR: Invariant diagnostics crashed.")
            else
                throw(ex)
            end
        end
        return invds
    end

    """
    Compute the relative error of an invariant for the solution of an ODE or DAE system.

    Arguments: `(t::TimeSeries, q::AbstractDataSeries{T}, invariant::Function)`

    The `invariant` functions needs to take two arguments `(t,q)` and return the
    corresponding value of the invariant.

    Returns a tuple of two 1d DataSeries holding the time series of the invariant and the relativ error, respectively.
    """
    function compute_invariant_error(t::TimeSeries, q::AbstractDataSeries{T}, invariant::Function) where {T}
        invds = compute_invariant(t, q, invariant)
        errds = compute_relative_error(invds)
        (invds, errds)
    end

    """
    Compute the relative error of an invariant for the solution of a partitioned ODE or DAE system.

    Arguments: `(t::TimeSeries, q::AbstractDataSeries{T}, p::AbstractDataSeries{T}, invariant::Function)`

    The `invariant` functions needs to take three arguments `(t,q,p)` and return the
    corresponding value of the invariant.

    Returns a tuple of two 1d DataSeries holding the time series of the invariant and the relativ error, respectively.
    """
    function compute_invariant_error(t::TimeSeries, q::AbstractDataSeries{T}, p::AbstractDataSeries{T}, invariant::Function) where {T}
        invds = compute_invariant(t, q, p, invariant)
        errds = compute_relative_error(invds)
        (invds, errds)
    end

    """
    Computes the difference of the momentum and the one-form of an implicit ODE or DAE system.

    Arguments: `(t::TimeSeries, q::AbstractDataSeries{T}, p::AbstractDataSeries{T}, one_form::Function)`

    The `one_form` function needs to take three arguments `(t,q,k)` where `k` is the index
    of the one-form component.

    Returns a DataSeries similar to `p` holding the time series of the difference between the momentum and the one-form.
    """
    function compute_momentum_error(t::TimeSeries, q::AbstractDataSeries{T}, p::AbstractDataSeries{T,2}, one_form::Function) where {T}
        e = similar(p)

        for i in axes(p,2)
            for n in axes(p,1)
                for k in eachindex(p[n,i])
                    e[n,i][k] = p[n,i][k] - one_form(t[n], q[n,i], k)
                end
            end
        end

        return e
    end

    function compute_momentum_error(t::TimeSeries, q::AbstractDataSeries{T}, p::AbstractDataSeries{T,1}, one_form::Function) where {T}
        e = similar(p)

        for n in axes(p,1)
            for k in eachindex(p[n])
                e[n][k] = p[n][k] - one_form(t[n], q[n], k)
            end
        end

        return e
    end

    """
    Computes the difference of the momentum and the one-form of an implicit ODE or DAE system.

    Arguments: `(p::AbstractDataSeries{DT}, ϑ::AbstractDataSeries{DT})`

    Returns a DataSeries similar to `p` holding the time series of the difference between `p` and `ϑ`.
    """
    function compute_momentum_error(p::AbstractDataSeries{DT}, ϑ::AbstractDataSeries{DT}) where {DT}
        @assert axes(p) == axes(ϑ)

        e = similar(p)
        parent(e) .= parent(p) .- parent(ϑ)

        return e
    end

    """
    Computes the drift in an invariant error.

    Arguments: `(t::TimeSeries, invariant_error::AbstractDataSeries{T,1}, interval_length=100)`

    The time series of the solution is split into intervals of `interval_length` time steps.
    In each interval, the maximum of the absolute value of the invariant error is determined.
    Returns a tuple of a TimeSeries that holds the centers of all intervals and a 1d DataSeries
    that holds the maxima.

    This is useful to detect drifts in invariants that are not preserved exactly but whose error
    is oscillating such as the energy error of Hamiltonian systems with symplectic integrators.
    """
    function compute_error_drift(t::TimeSeries, invariant_error::AbstractDataSeries{T,1}, interval_length=100) where {T}
        @assert ntime(t) == ntime(invariant_error)
        @assert mod(t.n, interval_length) == 0

        nint   = div(ntime(invariant_error), interval_length)
        Tdrift = TimeSeries(nint, (t[end] - t[begin]) / nint)
        Idrift = DataSeries(T, nint)

        Tdrift[0] = t[0]

        for i in 1:nint
            i1 = interval_length*(i-1)+1
            i2 = interval_length*i
            Idrift[i] = maximum(abs.(invariant_error[i1:i2]))
            Tdrift[i] = div(t[i1] + t[i2], 2)
        end

        (Tdrift, Idrift)
    end

end


text:
using ComoniconTestUtils
using Comonicon.ZSHCompletions
using Test
using Random

Random.seed!(42)

@testset "test completion" for _ in 1:5
    cmd = rand_command()
    script = ZSHCompletions.emit(cmd)
    @test occursin("#compdef _$(cmd.root.name) $(cmd.root.name)", script)
    @test occursin("function _$(cmd.root.name)() {", script)
end


text:
using Test
using MLDatasets
using MLDatasets.FastAIDatasets
using MLDatasets.FastAIDatasets: datasetpath, loaddataclassification,
    mapobs, filterobs, groupobs, joinobs, getobs, nobs

@testset "datasetpath" begin
    @test_nowarn FastAIDatasets.datasetpath("mnist_var_size_tiny")
end

@testset "loaddataset" begin
    @test_nowarn loaddataclassification(datasetpath("mnist_var_size_tiny"))
end

@testset "Data container transformations" begin
    @testset "mapobs" begin
        data = 1:10
        mdata = mapobs(-, data)
        @test getobs(mdata, 8) == -8

        mdata2 = mapobs((-, x -> 2x), data)
        @test getobs(mdata2, 8) == (-8, 16)

        nameddata = mapobs((x = sqrt, y = log), data)
        @test getobs(nameddata, 10) == (x = sqrt(10), y = log(10))
        @test getobs(nameddata.x, 10) == sqrt(10)
    end

    @testset "filterobs" begin
        data = 1:10
        fdata = filterobs(>(5), data)
        @test nobs(fdata) == 5
    end

    @testset "groupobs" begin
        data = -10:10
        datas = groupobs(>(0), data)
        @test length(datas) == 2
    end

    @testset "joinobs" begin
        data1, data2 = 1:10, 11:20
        jdata = joinobs(data1, data2)
        @test getobs(jdata, 15) == 15
    end
end


text:
# curand functions:

using Random
import Random: rand!, randn!
rand!(a::KnetArray{Float32})=(@curand(curandGenerateUniform,(Cptr,Ptr{Cfloat},Csize_t),rng(),a,length(a)); a)
rand!(a::KnetArray{Float64})=(@curand(curandGenerateUniformDouble,(Cptr,Ptr{Cdouble},Csize_t),rng(),a,length(a)); a)

function randn!(a::KnetArray{Float32}, mean = 0, stdev = 1)
    @curand(curandGenerateNormal,(Cptr,Ptr{Cfloat},Csize_t, Cfloat, Cfloat),
        rng(), a, length(a), Cfloat(mean), Cfloat(stdev))
    a
end

function randn!(a::KnetArray{Float64}, mean = 0, stdev = 1)
    @curand(curandGenerateNormalDouble,(Cptr,Ptr{Cdouble},Csize_t, Cdouble, Cdouble),
        rng(), a, length(a), Cdouble(mean), Cdouble(stdev))
    a
end

let RNG=0
global rng
function rng(init=false)
    if RNG==0 || init
        ptr = Cptr[0]
        # CURAND_RNG_PSEUDO_DEFAULT = 100, ///< Default pseudorandom generator
        @curand(curandCreateGenerator,(Cptr,Cint),ptr,100)
        RNG = ptr[1]
    end
    return RNG
end
end

"""
    Knet.seed!(n::Integer)

Run seed!(n) on both cpu and gpu.
"""
function seed!(n::Integer)
    # need to regenerate RNG for the seed to take effect for some reason
    if gpu() >= 0
        @curand(curandSetPseudoRandomGeneratorSeed,(Cptr,Culonglong),rng(true),n)
    end
    Random.seed!(n)
end

@deprecate setseed Knet.seed!


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "BlueZ"
version = v"5.54.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.kernel.org/pub/linux/bluetooth/bluez-$(version.major).$(version.minor).tar.xz", "68cdab9e63e8832b130d5979dc8c96fdb087b31278f342874d992af3e56656dc")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd bluez-*
# Hint to find libstc++, required to link against C++ libs when using C compiler
if [[ "${target}" == *-linux-* ]]; then
    if [[ "${nbits}" == 32 ]]; then
        export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib";
    else
        export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64";
    fi;
fi
# linux/if_alg doesn't seem to work; prevent configure from finding it
sed -i -e "s~ linux/if_alg.h~~" configure.ac
autoconf
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-systemd --enable-library
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("btmon", :btmon),
    ExecutableProduct("l2ping", :l2ping),
    ExecutableProduct("hid2hci", :hid2hci, "lib/udev"),
    ExecutableProduct("l2test", :l2test),
    ExecutableProduct("bccmd", :bccmd),
    ExecutableProduct("btattach", :btattach),
    ExecutableProduct("hex2hcd", :hex2hcd),
    ExecutableProduct("bluetoothd", :bluetoothd, "libexec/bluetooth"),
    ExecutableProduct("bluemoon", :bluemoon),
    ExecutableProduct("rctest", :rctest),
    ExecutableProduct("bluetooth", :bluetooth, "lib/cups/backend"),
    ExecutableProduct("mpris-proxy", :mpris_proxy),
    ExecutableProduct("obexd", :obexd, "libexec/bluetooth"),
    ExecutableProduct("bluetoothctl", :bluetoothctl),
    LibraryProduct("libbluetooth", :libbluetooth),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Dbus_jll", uuid="ee1fde0b-3d02-5ea6-8484-8dfef6360eab"))
    Dependency(PackageSpec(name="eudev_jll", uuid="35ca27e7-8b34-5b7f-bca9-bdc33f59eb06"))
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"), v"2.68.1"; compat="2.68.1")
    Dependency(PackageSpec(name="Libical_jll", uuid="bce108ef-3f60-5dd0-bcd6-e13a096cb796"))
    Dependency(PackageSpec(name="Readline_jll", uuid="05236dd9-4125-5232-aa7c-9ec0c9b2c25a"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module DocThemePC2

using DocumenterTools: Themes

function install(docs_path)
    # create the themes
    for w in ("light", "dark")
        header = read(joinpath(@__DIR__, "../assets/pc2-header.scss"), String)
        theme = read(joinpath(@__DIR__, "../assets/pc2-$(w)defs.scss"), String)
        write(joinpath(docs_path, "pc2-$(w).scss"), header * "\n" * theme)
    end
    # compile the themes
    Themes.compile(joinpath(docs_path, "pc2-light.scss"), joinpath(docs_path, "src/assets/themes/documenter-light.css"))
    Themes.compile(joinpath(docs_path, "pc2-dark.scss"), joinpath(docs_path, "src/assets/themes/documenter-dark.css"))
    # remove tmp files
    rm(joinpath(docs_path, "pc2-light.scss"))
    rm(joinpath(docs_path, "pc2-dark.scss"))
end

end


text:
"""
测试kagome的相关功能
"""

using MARY_fRG.Basics
using MARY_fRG.Fermi
using MARY_fRG.Fermi.Kagome
using MARY_fRG.Fermi.Patch
using MARY_fRG.Fermi.Surface
using MARY_fRG.Drawers
using Plots
using MARY_fRG.FlowEquation

"""
运行
"""
function run()
    nu1, nu2, nu3 = get_kagome_ν(0.1, 0.3)
    println(sum(nu1 .* nu1))
    println(sum(nu1 .* nu2))
    println(sum(nu1 .* nu3))
    println(nu1[1]*nu1[1] + nu2[1]*nu2[1] + nu3[1]*nu3[1])
    println(nu1[1]*nu1[2] + nu2[1]*nu2[2] + nu3[1]*nu3[2])
    #费米面的路径
    mpt = Point2D(0, 2*3.1415926 / 1.7320508)
    mpt2 = Point2D(-3.1415926, 3.1415926 / 1.7320508)
    mpt3 = Point2D(-3.1415926, -3.1415926 / 1.7320508)
    mpt4 = Point2D(0, -2*3.1415926 / 1.7320508)
    mpt5 = Point2D(3.1415926, -3.1415926 / 1.7320508)
    mpt6 = Point2D(3.1415926, 3.1415926 / 1.7320508)
    xvals = [0.0]
    pps = [mpt]
    for idx in 0:1:18
        omega = (idx+1) / 20
        push!(xvals, omega)
        push!(pps, middle_point(mpt, mpt2; sc1=1-omega, sc2=omega))
    end
    #println(pps)
    #nu1, nu2, nu3 = get_kagome_ν(mpt2.x, mpt2.y)
    #println(nu1)
    #println(pps[11])
    #nu1, nu2, nu3 = get_kagome_ν(pps[11].x, pps[11].y)
    #println(nu2[1])
    #return
    push!(xvals, 1.0)
    push!(pps, mpt2)
    for idx in 0:1:18
        omega = (idx+1) / 20
        push!(xvals, 1+omega)
        push!(pps, middle_point(mpt2, mpt3; sc1=1-omega, sc2=omega))
    end
    push!(xvals, 2.0)
    push!(pps, mpt3)
    for idx in 0:1:18
        omega = (idx+1) / 20
        push!(xvals, 2+omega)
        push!(pps, middle_point(mpt3, mpt4, sc1=1-omega, sc2=omega))
    end
    push!(xvals, 3.0)
    push!(pps, mpt4)
    for idx in 0:1:18
        omega = (idx+1) / 20
        push!(xvals, 3+omega)
        push!(pps, middle_point(mpt4, mpt5, sc1=1-omega, sc2=omega))
    end
    push!(xvals, 4.0)
    push!(pps, mpt5)
    for idx in 0:1:18
        omega = (idx+1) / 20
        push!(xvals, 4+omega)
        push!(pps, middle_point(mpt5, mpt6, sc1=1-omega, sc2=omega))
    end
    push!(xvals, 5.0)
    push!(pps, mpt6)
    for idx in 0:1:18
        omega = (idx+1) / 20
        push!(xvals, 5+omega)
        push!(pps, middle_point(mpt6, mpt, sc1=1-omega, sc2=omega))
    end
    push!(xvals, 6.0)
    push!(pps, mpt)
    println(pps)
    sitea = []
    siteb = []
    sitec = []
    for kkt in pps
        nu1, nu2, nu3 = get_kagome_ν(kkt.x, kkt.y)
        push!(sitea, abs(nu2[1]))
        push!(siteb, abs(nu2[2]))
        push!(sitec, abs(nu2[3]))
    end
    println(siteb)
    plt = plot(xvals, sitea)
    plot!(plt, xvals, siteb)
    plot!(plt, xvals, sitec)
    png(plt, "kag_band_contrib1")
    #绘制一下d能带的贡献
    sitea = []
    siteb = []
    sitec = []
    #nu1, nu2, nu3 = get_nu_numpy(pps[10].coord[0], pps[10].coord[1])
    #print(nu2[1])
    #print(check_patches_converge(pps[10]))
    for kkt in pps
        nu1, nu2, nu3 = get_kagome_ν(kkt.x, kkt.y)
        push!(sitea, abs(nu3[1]))
        push!(siteb, abs(nu3[2]))
        push!(sitec, abs(nu3[3]))
    end
    plt = plot(xvals, sitea)
    plot!(plt, xvals, siteb)
    plot!(plt, xvals, sitec)
    png(plt, "kag_band_contrib2")
    #print(pps[10])
    #print(siteb[10])
    #nu1, nu2, nu3 = get_nu_numpy(pps[10].coord[0], pps[10].coord[1])
    #print(nu2[1])
    #print(check_patches_converge(pps[10]))
end

#run()

"""
相互作用
"""
function interact()
    ks = upperband_kagome_lattice(0.)
    pats = patches_under_vonhove(ks.brillouin, ks.dispersion[1], 24)
    #pats = [
    #    Point2D(0.1, -0.1),
    #    Point2D(1., 1.),
    #    Point2D(-1., -1.)
    #]
    plt = draw_lines(ks.brillouin.edges)
    draw_points!(plt, pats)
    @savepng plt "kag_surf1"
    #pats = get_von_hove_patches(24)
    #vmat = band_v(6.0, pats)
    #draw_heatmap(vmat[:, :, 0])
    #println(vmat[1, 2, 1], vmat[2, 1, 2])
    #println(vmat[2, 0, 0], vmat[0, 2, 2])
    #println(vmat[2, 0, 1], vmat[0, 2, 1])
    println(pats)
    umat = get_kagome_U_mt(6.0, pats)
    println(all(isfinite.(umat)))
    maxv, maxi = findmin(isfinite.(umat))
    println(maxi, maxv)
    println(umat[maxi])
    #
    #k1v = pats[maxi[1]]
    #k2v = pats[maxi[2]]
    #k3v = pats[maxi[3]]
    #k4v = Fermi.hexagon_kadd(k1v, k2v)
    #k4v = Fermi.hexagon_kadd(k4v, -k3v)
    #_, k1nu, _ = get_kagome_ν(k1v.x, k1v.y)
    #_, k2nu, _ = get_kagome_ν(k2v.x, k2v.y)
    #_, k3nu, _ = get_kagome_ν(k3v.x, k3v.y)
    #_, k4nu, _ = get_kagome_ν(k4v.x, k4v.y)
    #println(k1v)
    #println(k2v)
    #println(k3v)
    #println(k4v)
    #println(k1nu)
    #println(k2nu)
    #println(k3nu)
    #println(k4nu)
    #
    plt = heatmap(umat[:, :, 1])
    png(plt, "kag_interact")
    #amat = umat + vmat
    #draw_heatmap(amat[:, :, 0])
    #smat = get_sublattice_u(umat, pats)
    #print("site")
    #draw_heatmap(smat[0, 0, 0, 0, :, :, 0])
    #draw_heatmap(smat[1, 1, 1, 1, :, :, 2])
    #draw_heatmap(smat[0, 1, 1, 0, :, :, 0])
    #draw_heatmap(smat[2, 2, 2, 2, :, :, 5])
end

#interact()

#get_kagome_ν(0., 0.)

"""
计算ult
"""
function get_ult_bubb(Γ4; maxnum=1000000)
    Γ4_ult = fliter_away_surface(Γ4, 2.0)
    Γ4_ult = refine_to_surface(Γ4_ult)
    while length(Γ4_ult.ltris) < maxnum
        Γ4_ult = refine_to_surface(Γ4_ult)
    end
    #
    maxi, mini = engpeak_to_surface(Γ4_ult)
    bubb_pp_u, bubb_fs_u, bubb_ex_u = all_bubble_tf_ult_mt(Γ4_ult, maxi)
    return bubb_pp_u, bubb_fs_u, bubb_ex_u
end

function kprofile()
    model = upperband_kagome_lattice(0.)
    Γ4 = TFGamma4(
        model, 6.0, 24, 200
    )
    @time bpu, bfu, beu = get_ult_bubb(Γ4)
    @time bpu, bfu, beu = get_ult_bubb(Γ4)
end

kprofile()


text:
cone_product(a, b) = [dot(a, b); a[1] * b[2:end] + b[1] * a[2:end]]

function cone_product!(a, z, s; reset=false)
    reset && fill!(a, 0.0)
    n = length(a)
    a[1] += z' * s
    for i = 2:n
        a[i] += z[1] * s[i]
        a[i] += s[1] * z[i]
    end
    return nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ResourceGroupExportResult <: SwaggerModel
    template::Any # spec type: Union{ Nothing, Any } # spec name: template
    error::Any # spec type: Union{ Nothing, ErrorResponse } # spec name: error

    function ResourceGroupExportResult(;template=nothing, error=nothing)
        o = new()
        validate_property(ResourceGroupExportResult, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(ResourceGroupExportResult, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type ResourceGroupExportResult

const _property_map_ResourceGroupExportResult = Dict{Symbol,Symbol}(Symbol("template")=>Symbol("template"), Symbol("error")=>Symbol("error"))
const _property_types_ResourceGroupExportResult = Dict{Symbol,String}(Symbol("template")=>"Any", Symbol("error")=>"ErrorResponse")
Base.propertynames(::Type{ ResourceGroupExportResult }) = collect(keys(_property_map_ResourceGroupExportResult))
Swagger.property_type(::Type{ ResourceGroupExportResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ResourceGroupExportResult[name]))}
Swagger.field_name(::Type{ ResourceGroupExportResult }, property_name::Symbol) =  _property_map_ResourceGroupExportResult[property_name]

function check_required(o::ResourceGroupExportResult)
    true
end

function validate_property(::Type{ ResourceGroupExportResult }, name::Symbol, val)
end


text:
using Test

using MarketData
using TimeSeries

using MarketTechnicals

@testset "Candlesticks" begin

@testset "doji" begin
    @test timestamp(ohlc[findall(doji(ohlc))])[1] == Date(2000, 5, 30)
    @test timestamp(ohlc[findall(doji(ohlc))])[2] == Date(2000, 12, 11)
    @test length(ohlc[findall(doji(ohlc))])      == 5
end

end  # @testset "Candlesticks"


text:
# using PyPlot

# using Interpolations

"""
Equations utilized are from Terman, 1943 Radio Engineers' Handbook

StrandDiam - Strand diameter (m) if one conductor is used, this is the conductor thickness

CoilTurns

CoilLength -Axial length of the coil (m)

CoilDiam - Diameter of the full solenoid (mean for elliptical ones)

NumStrands (default= 1) - Number of strands

BundleDiam (default= StrandDiam)

;
f = 1:1000:1e6, Test freuencies

PlotOn=true
"""
function LitzWire_ACRes(
        StrandDiam,
        CoilTurns,
        CoilLength,
        CoilDiam,
        NumStrands = 1,
        BundleDiam = StrandDiam;
        f = 1:1000:1e6,
        PlotOn=false
)

        if BundleDiam < (NumStrands * StrandDiam)
                error("Too many strands--will not fit if they are side-by-side")
        end

        ρ = 1.7e-8 #Ohm-meter
        WireL = pi * CoilDiam * CoilTurns
        CoilProps = Dict("WireL" => WireL)
        CoilProps["StrandDiam"] = StrandDiam
        CoilProps["BundleDiam"] = BundleDiam
        CoilProps["NumStrands"] = NumStrands
        StrandArea = π/4 *StrandDiam^ 2
        Rdc = ρ * WireL ./ (StrandArea) ./ NumStrands
        CoilProps["Rdc"] = Rdc

        u = Table20Terman(CoilTurns)

        μ₀ = 4*π*1e-7
        SkinDepth = sqrt.(ρ ./ (π * f * μ₀))
        Pitch = CoilLength / CoilTurns
        if Pitch < BundleDiam
                error("Pitch cannot be finer than a single bundle of wire")
        end
        Tmp = [[Table18Terman(StrandDiam,ρ,ff;units="meters")[1],
        Table18Terman(StrandDiam,ρ,ff;units="meters")[2],
        Table18Terman(StrandDiam,ρ,ff;units="meters")[3]] for ff in f]
        Tmp = hcat(Tmp...)
        X =Tmp[1,:]
        G = Tmp[2,:]
        H = Tmp[3,:]
        k = Table23Terman(NumStrands)
        Rac =[Rdc .* ( H[i] .+ (k + u .* (StrandDiam ./ Pitch) .^ 2) .*
                (StrandDiam/BundleDiam)^2 .* NumStrands^2 * G[i])
                for i in 1:length(f)]

        if PlotOn==true
                figure(1)
                PyPlot.plot(f, Rac,label="AC Resistance","k")

                PyPlot.plot([f[1], f[end]], [Rdc, Rdc] , "b--", label = "DC")
                xlabel("Freq. (Hz)")
                ylabel("Resistance (Ω)")
                legend()
        end

        return CoilProps, Rac, SkinDepth
end

function Table18Terman(X::AbstractFloat)
Table = [0.5   1       0.00097
         0.6   1.001   0.00202
         0.8   1.002   0.00632
         1.0   1.005   0.01519
         1.1   1.008   0.02196
         1.3   1.015   0.04127
         1.5   1.026   0.0691
         1.7   1.042   0.1055
         1.9   1.064   0.1489
         2.2   1.111   0.2214
         2.5   1.175   0.2949
         2.8   1.256   0.3412
         3.0   1.318   0.4049
         3.3   1.420   0.4626
         3.5   1.492   0.4987
         3.8   1.603   0.5503
         4.0   1.679   0.5842
         4.2   1.752   0.618
         4.5   1.863   0.669
         4.8   1.971   0.720
         5.0   2.043   0.755
         5.4   2.184   0.826
         6.0   2.394   0.932
         6.6   2.603   1.038
         7.0   2.743   1.109
         7.6   2.954   1.216
         8.0   3.094   1.287
         8.6   3.306   1.393
         9.0   3.446   1.464
         10.0  3.799   1.641
         13.0  4.856   2.171
         16.0  5.915   2.702
         20.0  7.328   3.409
         23.0  8.388   3.940
         25.0  9.094   4.294
         30.0  10.86   5.177
         40.0  14.40   6.946
         50.0  17.93   8.713
         60.0  21.46   10.48
         70.0  25.00   12.25
         80.0  28.54   14.02
         90.0  32.07   15.78
         100.  35.61   17.55]

         if X<0.5
                 G=X^4/64
                 H = 1.0
         elseif X>100
                 G = (sqrt(2)*X+1)/4
                 H = (sqrt(2)*X-1)/8
         else
                 interp_linearG = LinearInterpolation(Table[:,1], Table[:,3])
                 interp_linearH = LinearInterpolation(Table[:,1], Table[:,2])
                 G =interp_linearG(X)
                 H =interp_linearH(X)
         end

         return X,G,H
end

function Table18Terman(d,ρ,f;units="meters")
        if units=="meters"
                d = d*100 ##Units are cm
                ρ = ρ*1e11 ## units are abΩ⋅cm
        elseif (units=="centimeters")|(units=="cm")
        else
                error("Units not recognized")
        end
        X = π*d*sqrt(2*f/ρ)
        return Table18Terman(X)
end

"""
This function interpolates Table 20 in Terman Radio Engineers' Handbook
        It applies for solenoids with relatively low density (pitch > 2x wire diameter)
        It takes in number of turns, and returns `u`
"""
function Table20Terman(N)
Table = [2 1
         4 1.8
         6 2.16
         8 2.37
         10 2.51
         12 2.61
         16 2.74
         24 2.91
         32 3.0]
         if N>30
                 return 3.29
         else
                 interp_linear = LinearInterpolation(Table[:,1], Table[:,2])
                 return interp_linear(N)
         end
 end

function Table23Terman(n)
        Table = [3   1.55
                 9   1.84
                 27  1.92]
                 if n>27
                         return 2
                 elseif n<3
                         return 0
                 else
                         interp_linear = LinearInterpolation(Table[:,1], Table[:,2])
                         return interp_linear(n)
                 end

end


text:
###### problem parameters ####
k= # number of defence faciity
I= # number of customer
J= # number of facility services
###### initialize paramaters
m= # number of memplex
n= # number of frog in each memplex
q= # number of submemlex members
###### generate initial population

###### partition population into memplexs

###### memetic evolution withing each memplex

##### shuffling all frogs

#### check termination consdition

max_stars_repo_name: imanpanahi/MilProject

max_issues_repo_name: imanpanahi/MilProject

max_forks_repo_name: imanpanahi/MilProject


text:
module Renderer

export respond, json, redirect_to, html, flax, include_asset, has_requested, css_asset, js_asset
export respond_with_json, respond_with_html, respond_with, html!, json!, http_error, error!

using Nullables, JSON, HTTP
using Genie, Genie.Util, Genie.Configuration, Genie.Loggers, Genie.Macros

eval(:(include("$(Genie.config.html_template_engine).jl")))
Genie.config.html_template_engine != Genie.config.json_template_engine && Core.eval(:(include("$(Genie.config.json_template_engine).jl")))

eval(:(using .$(Genie.config.html_template_engine), .$(Genie.config.json_template_engine)))
eval(:(const HTMLTemplateEngine = $(Genie.config.html_template_engine)))
eval(:(const JSONTemplateEngine = $(Genie.config.json_template_engine)))

export HTMLTemplateEngine, JSONTemplateEngine

const DEFAULT_LAYOUT_FILE = Genie.config.renderer_default_layout_file

const CONTENT_TYPES = Dict{Symbol,String}(
  :html   => "text/html",
  :plain  => "text/plain",
  :json   => "application/json",
  :js     => "text/javascript",
  :xml    => "text/xml",
)
const DEFAULT_CONTENT_TYPE = :html

"""
"""
function respond_with(response_type::Symbol, args...; kargs...)
  if lowercase(string(response_type)) == "html"
    respond_with_html(args...; kargs...)
  elseif lowercase(string(response_type)) == "json"
    respond_with_json(args...; kargs...)
  end
end
function respond_with(response_type::Symbol, err::T) where {T<:Exception}
  if lowercase(string(response_type)) == "html"
    Genie.Router.error_404(err.msg)
  elseif lowercase(string(response_type)) == "json"
    respond(Dict(:json => JSON.json(err)), 404, Dict{String,String}("Content-Type" => "application/json"))
  end
end

"""
    html(resource::Symbol, action::Symbol, layout::Symbol = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}

Invokes the HTML renderer of the underlying configured templating library.
"""
function html(resource::Symbol, action::Symbol, layout::Symbol = DEFAULT_LAYOUT_FILE; vars...) :: Dict{Symbol,String}
  HTMLTemplateEngine.render_html(resource, action, layout; vars...)
end
function html(view::String, layout::String; vars...) :: Dict{Symbol,String}
  HTMLTemplateEngine.render_html(view, layout; vars...)
end

"""
    respond_with_html(resource::Symbol, action::Symbol, layout::Symbol = DEFAULT_LAYOUT_FILE, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Response

Invokes the HTML renderer of the underlying configured templating library and wraps it into a `HTTP.Response`.
"""
function respond_with_html(resource::Symbol, action::Symbol, layout::Symbol = DEFAULT_LAYOUT_FILE; vars...) :: HTTP.Response
  html(resource, action, layout; vars...) |> respond
end
function respond_with_html(view::String, layout::String; vars...) :: HTTP.Response
  html(view, layout; vars...) |> respond
end
const html! = respond_with_html

function html!(output::String) :: HTTP.Response
  respond(output)
end

function flax(resource::Union{Symbol,String}, action::Union{Symbol,String}, layout::Union{Symbol,String} = DEFAULT_LAYOUT_FILE; vars...) :: Dict{Symbol,String}
  HTMLTemplateEngine.render_flax(resource, action, layout; vars...)
end

"""
    json(resource::Symbol, action::Symbol, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Dict{Symbol,String}

Invokes the JSON renderer of the underlying configured templating library.
"""
function json(resource::Union{Symbol,String}, action::Union{Symbol,String}; vars...) :: Dict{Symbol,String}
  JSONTemplateEngine.render_json(resource, action; vars...)
end

"""
    respond_with_json(resource::Symbol, action::Symbol, check_nulls::Vector{Pair{Symbol,Nullable}} = Vector{Pair{Symbol,Nullable}}(); vars...) :: Response

Invokes the JSON renderer of the underlying configured templating library and wraps it into a `HTTP.Response`.
"""
function respond_with_json(resource::Union{Symbol,String}, action::Union{Symbol,String}; vars...) :: HTTP.Response
  json(resource, action; vars...) |> respond
end
const json! = respond_with_json

function json!(output::String) :: HTTP.Response
  respond(output, :json)
end
function json!(output::Any) :: HTTP.Response
  JSON.json(output) |> json!
end

"""
    redirect_to(location::String, code::Int = 302, headers = Dict{String,String}()) :: Response

Sets redirect headers and prepares the `Response`.
"""
function redirect_to(location::String, code = 302, headers = Dict{String,String}()) :: HTTP.Response
  headers["Location"] = location
  respond(Dict{Symbol,String}(:plain => "Redirecting you to $location"), code, headers)
end
function redirect_to(named_route::Symbol, code = 302, headers = Dict{String,String}()) :: HTTP.Response
  redirect_to(Genie.Router.link_to(named_route), code, headers)
end

"""
    has_requested(content_type::Symbol) :: Bool

Checks wheter or not the requested content type matches `content_type`.
"""
function has_requested(content_type::Symbol) :: Bool
  task_local_storage(:__params)[:response_type] == content_type
end

"""
    respond{T}(body::Dict{Symbol,T}, code::Int = 200, headers = Dict{String,String}()) :: Response

Constructs a `Response` corresponding to the content-type of the request.
"""
function respond(body::Dict{Symbol,T}, code::Int = 200, headers = Dict{String,String}())::HTTP.Response where {T}
  sbody::String =   if haskey(body, :json)
                      headers["Content-Type"] = CONTENT_TYPES[:json]
                      body[:json]
                    elseif haskey(body, :html)
                      headers["Content-Type"] = CONTENT_TYPES[:html]
                      body[:html]
                    elseif haskey(body, :js)
                      headers["Content-Type"] = CONTENT_TYPES[:js]
                      body[:js]
                    elseif haskey(body, :plain)
                      headers["Content-Type"] = CONTENT_TYPES[:plain]
                      body[:plain]
                    else
                      log("Unsupported Content-Type", :err)
                      log(body)
                      Loggers.@location

                      error("Unsupported Content-Type")
                    end

                    HTTP.Response(code, [h for h in headers], body = sbody)
end
function respond(body::String, content_type::Symbol) :: HTTP.Response
  HTTP.Response(200, ["Content-Type" => CONTENT_TYPES[content_type]], body = body)
end
function respond(body::String, content_type::String) :: HTTP.Response
  HTTP.Response(200, ["Content-Type" => content_type], body = body)
end
function respond(response::Tuple, headers = Dict{AbstractString,AbstractString}()) :: HTTP.Response
  respond(response[1], response[2], [h for h in headers])
end
function respond(response::HTTP.Response) :: HTTP.Response
  response
end
function respond(body::String, params::Dict{Symbol,T})::HTTP.Response where {T}
  r = params[:RESPONSE]
  r.data = body

  r |> respond
end
function respond(body::String) :: HTTP.Response
  respond(HTTP.Response(body))
end
function respond(args...; kargs...) :: HTTP.Response
  respond_with(Genie.Router.response_type(), args...; kargs...)
end
function respond(err::T)::HTTP.Response where {T<:Exception}
  respond_with(Genie.Router.response_type(), err)
end

"""
    http_error(status_code; id = "resource_not_found", code = "404-0001", title = "Not found", detail = "The requested resource was not found")

Constructs an error `Response`.
"""
function http_error(status_code; id = "", code = "", title = "", msg = "") :: HTTP.Response
  respond(msg, status_code, Dict{String,String}())
end
const error! = http_error

"""
    include_asset(asset_type::Symbol, file_name::String; fingerprinted = Genie.config.assets_fingerprinted) :: String

Returns the path to an asset. `asset_type` can be one of `:js`, `:css`. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicated wheter or not fingerprinted (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function include_asset(asset_type::Symbol, file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  suffix = fingerprinted ? "-" * Genie.ASSET_FINGERPRINT * ".$(asset_type)" : ".$(asset_type)"
  "/$asset_type/$(file_name)$(suffix)"
end
function include_asset(asset_type::Symbol, file_name::Symbol; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  include_asset(asset_type, string(file_name), fingerprinted = fingerprinted)
end

"""
    css_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String

Path to a css asset. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicated wheter or not fingerprinted (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function css_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  include_asset(:css, file_name, fingerprinted = fingerprinted)
end

"""
    js_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String

Path to a js asset. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicated wheter or not fingerprinted (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function js_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  include_asset(:js, file_name, fingerprinted = fingerprinted)
end

end


text:
mtg = read_mtg("files/simple_plant.mtg")

@testset "names" begin
    @test get_attributes(mtg) == [:scales, :description, :symbols, :Length, :Width, :XEuler]
    @test names(mtg) == get_attributes(mtg)
    @test names(DataFrame(mtg))[8:end] == string.(names(mtg))
end

@testset "get attribute value" begin
    @test mtg[1][1][1][:Width] == 0.02
    @test mtg[1][1][1].attributes[:Width] == 0.02
end

@testset "set attribute value" begin
    mtg[1][1][1][:Width] = 1.0
    @test mtg[1][1][1][:Width] == 1.0
end


text:
using Distributed
@everywhere using Dates
@everywhere using DelimitedFiles
@everywhere using CSV
#@everywhere using Plots
@everywhere using Statistics
@everywhere using DocStringExtensions
@everywhere using DataFrames
@everywhere using Random

@everywhere module_dir = "/home/sarah/HBVModel/src/"
@everywhere push!(LOAD_PATH, $module_dir)

# load list of structs
@everywhere include(module_dir*"Model_Core/structs.jl")
# load components of models represented by buckets
@everywhere include(module_dir*"Model_Core/processes_buckets.jl")
# load functions that combine all components of one HRU
@everywhere include(module_dir*"Model_Core/elevations.jl")
# load functions for combining all HRUs and for running the model
@everywhere include(module_dir*"Model_Core/allHRU.jl")
# load function for running model which just returns the necessary output for calibration
@everywhere include(module_dir*"Model_Core/run_model.jl")
# load functions for preprocessing temperature and precipitation data
@everywhere include(module_dir*"Model_Core/Preprocessing.jl")
# load functions for calculating the potential evaporation
@everywhere include(module_dir*"Model_Core/Potential_Evaporation.jl")
# load objective functionsM
@everywhere include(module_dir*"Model_Core/ObjectiveFunctions.jl")
# load parameterselection
@everywhere include(module_dir*"Model_Core/parameterselection.jl")
# load running model in several precipitation zones
@everywhere include(module_dir*"Model_Core/runmodel_Prec_Zones.jl")

@everywhere function run_MC(ID, nmax)

        local_path = "/home/sarah/"
        # ------------ CATCHMENT SPECIFIC INPUTS----------------
        ID_Prec_Zones = [17700, 114926]
        # size of the area of precipitation zones
        Area_Zones = [235811198.0, 31497403.0]
        Area_Catchment = sum(Area_Zones)
        Area_Zones_Percent = Area_Zones / Area_Catchment
        Snow_Threshold = 600
        Height_Threshold = 2700

        Mean_Elevation_Catchment =  2300 # in reality 2233.399986
        Elevations_Catchment = Elevations(200.0, 1000.0, 3600.0, 1385., 1385.) # take temp at 17700
        Sunhours_Vienna = [8.83, 10.26, 11.95, 13.75, 15.28, 16.11, 15.75, 14.36, 12.63, 10.9, 9.28, 8.43]
        # where to skip to in data file of precipitation measurements
        Skipto = [0, 24]
        # get the areal percentage of all elevation zones in the HRUs in the precipitation zones
        Areas_HRUs =  CSV.read(local_path*"HBVModel/Defreggental/HBV_Area_Elevation_round.csv", skipto=2, decimal='.', delim = ',')
        # get the percentage of each HRU of the precipitation zone
        Percentage_HRU = CSV.read(local_path*"HBVModel/Defreggental/HRU_Prec_Zones.csv", header=[1], decimal='.', delim = ',')
        Elevation_Catchment = convert(Vector, Areas_HRUs[2:end,1])
        startyear = 1983
        endyear = 2005
        scale_factor_Discharge = 0.65
        # timeperiod for which model should be run (look if timeseries of data has same length)
        Timeseries = collect(Date(startyear, 1, 1):Day(1):Date(endyear,12,31))
        #------------ TEMPERATURE AND POT. EVAPORATION CALCULATIONS ---------------------

        Temperature = CSV.read(local_path*"HBVModel/Defreggental/prenner_tag_17700.dat", header = true, skipto = 3, delim = ' ', ignorerepeated = true)

        # get data for 20 years: from 1987 to end of 2006
        # from 1986 to 2005 13669: 20973
        #hydrological year 13577:20881
        Temperature = dropmissing(Temperature)
        Temperature_Array = Temperature.t / 10
        Precipitation_17700 = Temperature.nied / 10
        Timeseries_Temp = Date.(Temperature.datum, Dates.DateFormat("yyyymmdd"))
        startindex = findfirst(isequal(Date(startyear, 1, 1)), Timeseries_Temp)
        endindex = findfirst(isequal(Date(endyear, 12, 31)), Timeseries_Temp)
        Temperature_Daily = Temperature_Array[startindex[1]:endindex[1]]
        Dates_Temperature_Daily = Timeseries_Temp[startindex[1]:endindex[1]]
        Precipitation_17700 = Precipitation_17700[startindex[1]:endindex[1]]
        Precipitation_17700[findall(x -> x == -0.1, Precipitation_17700)] .= 0.0

        #
        # Dates_Temperature_Daily = Timeseries_Temp[startindex[1]:endindex[1]]
        # Dates_missing_Temp = Dates_Temperature_Daily[findall(x-> x == 999.9, Temperature_Daily)]
        # # get the temperature data at each elevation
        Elevation_Zone_Catchment, Temperature_Elevation_Catchment, Total_Elevationbands_Catchment = gettemperatureatelevation(Elevations_Catchment, Temperature_Daily)
        # # get the temperature data at the mean elevation to calculate the mean potential evaporation
        Temperature_Mean_Elevation = Temperature_Elevation_Catchment[:,findfirst(x-> x==Mean_Elevation_Catchment, Elevation_Zone_Catchment)]
        Potential_Evaporation = getEpot_Daily_thornthwaite(Temperature_Mean_Elevation, Dates_Temperature_Daily, Sunhours_Vienna)
        # ------------ LOAD OBSERVED DISCHARGE DATA ----------------
        Discharge = CSV.read(local_path*"HBVModel/Defreggental/Q-Tagesmittel-212100.csv", header= false, skipto=26, decimal=',', delim = ';', types=[String, Float64])
        Discharge = convert(Matrix, Discharge)
        startindex = findfirst(isequal("01.01."*string(startyear)*" 00:00:00"), Discharge)
        endindex = findfirst(isequal("31.12."*string(endyear)*" 00:00:00"), Discharge)
        Observed_Discharge = Array{Float64,1}[]
        push!(Observed_Discharge, Discharge[startindex[1]:endindex[1],2])
        Observed_Discharge = Observed_Discharge[1]
        Observed_Discharge = Observed_Discharge * 1000 / Area_Catchment * (3600 * 24)

        # ------------ LOAD TIMESERIES DATA AS DATES ------------------
        #Timeseries = Date.(Discharge[startindex[1]:endindex[1],1], Dates.DateFormat("d.m.y H:M:S"))
        firstyear = Dates.year(Timeseries[1])
        lastyear = Dates.year(Timeseries[end])

        # ------------- LOAD OBSERVED SNOW COVER DATA PER PRECIPITATION ZONE ------------
        # # find day wehere 2000 starts for snow cover calculations
        start2000 = findfirst(x -> x == Date(2000, 01, 01), Timeseries)
        length_2000_end = length(Observed_Discharge) - start2000 + 1
        observed_snow_cover = Array{Float64,2}[]
        for ID in ID_Prec_Zones
                current_observed_snow = readdlm(local_path*"HBVModel/Defreggental/snow_cover_fixed_Zone"*string(ID)*".csv",',', Float64)
                current_observed_snow = current_observed_snow[1:length_2000_end,3: end]
                push!(observed_snow_cover, current_observed_snow)
        end
        # ------------- LOAD PRECIPITATION DATA OF EACH PRECIPITATION ZONE ----------------------
        # get elevations at which precipitation was measured in each precipitation zone
        Elevations_17700 = Elevations(200., 1200., 3600., 1385., 1140)
        Elevations_114926 = Elevations(200, 1000, 2800, 1110., 1140)
        Elevations_All_Zones = [Elevations_17700, Elevations_114926]

        #get the total discharge
        Total_Discharge = zeros(length(Temperature_Daily))
        Inputs_All_Zones = Array{HRU_Input, 1}[]
        Storages_All_Zones = Array{Storages, 1}[]
        Precipitation_All_Zones = Array{Float64, 2}[]
        Precipitation_Gradient = 0.0
        Elevation_Percentage = Array{Float64, 1}[]
        Nr_Elevationbands_All_Zones = Int64[]
        Elevations_Each_Precipitation_Zone = Array{Float64, 1}[]
        Glacier_All_Zones = Array{Float64, 2}[]

        for i in 1: length(ID_Prec_Zones)
                if ID_Prec_Zones[i] == 114926
                        #print(ID_Prec_Zones[i])
                        Precipitation = CSV.read(local_path*"HBVModel/Defreggental/N-Tagessummen-"*string(ID_Prec_Zones[i])*".csv", header= false, skipto=Skipto[i], missingstring = "L\xfccke", decimal=',', delim = ';')
                        Precipitation_Array = convert(Matrix, Precipitation)
                        startindex = findfirst(isequal("01.01."*string(startyear)*" 07:00:00   "), Precipitation_Array)
                        endindex = findfirst(isequal("31.12."*string(endyear)*" 07:00:00   "), Precipitation_Array)
                        Precipitation_Array = Precipitation_Array[startindex[1]:endindex[1],:]
                        Precipitation_Array[:,1] = Date.(Precipitation_Array[:,1], Dates.DateFormat("d.m.y H:M:S   "))
                        # find duplicates and remove them
                        df = DataFrame(Precipitation_Array)
                        df = unique!(df)
                        # drop missing values
                        df = dropmissing(df)
                        Precipitation_Array = convert(Matrix, df)
                        Elevation_HRUs, Precipitation, Nr_Elevationbands = getprecipitationatelevation(Elevations_All_Zones[i], Precipitation_Gradient, Precipitation_Array[:,2])
                        push!(Precipitation_All_Zones, Precipitation)
                        push!(Nr_Elevationbands_All_Zones, Nr_Elevationbands)
                        push!(Elevations_Each_Precipitation_Zone, Elevation_HRUs)
                elseif ID_Prec_Zones[i] == 17700
                        Precipitation_Array = Precipitation_17700
                        # for all non data values use values of other precipitation zone
                        Elevation_HRUs, Precipitation, Nr_Elevationbands = getprecipitationatelevation(Elevations_All_Zones[i], Precipitation_Gradient, Precipitation_Array)
                        push!(Precipitation_All_Zones, Precipitation)
                        push!(Nr_Elevationbands_All_Zones, Nr_Elevationbands)
                        push!(Elevations_Each_Precipitation_Zone, Elevation_HRUs)
                end

                #glacier area only for 17700, for 114926 file contains only zeros
                Glacier_Area = CSV.read(local_path*"HBVModel/Defreggental/Glaciers_Elevations_"*string(ID_Prec_Zones[i])*"_evolution_69_06.csv",  header= true, delim=',')
                Years = collect(startyear:endyear)
                glacier_daily = zeros(Total_Elevationbands_Catchment)
                for current_year in Years
                        glacier_current_year = Glacier_Area[!, string(current_year)]
                        current_glacier_daily = repeat(glacier_current_year, 1, Dates.daysinyear(current_year))
                        glacier_daily = hcat(glacier_daily, current_glacier_daily)
                end
                push!(Glacier_All_Zones, glacier_daily[:,2:end])

                index_HRU = (findall(x -> x==ID_Prec_Zones[i], Areas_HRUs[1,2:end]))
                # for each precipitation zone get the relevant areal extentd
                Current_Areas_HRUs = convert(Matrix, Areas_HRUs[2: end, index_HRU])
                # the elevations of each HRU have to be known in order to get the right temperature data for each elevation
                Area_Bare_Elevations, Bare_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,1], Elevation_Catchment, Elevation_HRUs)
                Area_Forest_Elevations, Forest_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,2], Elevation_Catchment, Elevation_HRUs)
                Area_Grass_Elevations, Grass_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,3], Elevation_Catchment, Elevation_HRUs)
                Area_Rip_Elevations, Rip_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,4], Elevation_Catchment, Elevation_HRUs)
                #print(Bare_Elevation_Count, Forest_Elevation_Count, Grass_Elevation_Count, Rip_Elevation_Count)
                @assert 1 - eps(Float64) <= sum(Area_Bare_Elevations) <= 1 + eps(Float64)
                @assert 1 - eps(Float64) <= sum(Area_Forest_Elevations) <= 1 + eps(Float64)
                @assert 1 - eps(Float64) <= sum(Area_Grass_Elevations) <= 1 + eps(Float64)
                @assert 1 - eps(Float64) <= sum(Area_Rip_Elevations) <= 1 + eps(Float64)

                Area = Area_Zones[i]
                Current_Percentage_HRU = Percentage_HRU[:,1 + i]/Area
                # calculate percentage of elevations
                Perc_Elevation = zeros(Total_Elevationbands_Catchment)
                for j in 1 : Total_Elevationbands_Catchment
                        for h in 1:4
                                Perc_Elevation[j] += Current_Areas_HRUs[j,h] * Current_Percentage_HRU[h]
                        end
                end
                Perc_Elevation = Perc_Elevation[(findall(x -> x!= 0, Perc_Elevation))]
                @assert 0.99 <= sum(Perc_Elevation) <= 1.01
                push!(Elevation_Percentage, Perc_Elevation)
                # calculate the inputs once for every precipitation zone because they will stay the same during the Monte Carlo Sampling
                bare_input = HRU_Input(Area_Bare_Elevations, Current_Percentage_HRU[1],zeros(length(Bare_Elevation_Count)) , Bare_Elevation_Count, length(Bare_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100), (Snow_Threshold, Height_Threshold), 0, [0], 0, [0], 0, 0)
                forest_input = HRU_Input(Area_Forest_Elevations, Current_Percentage_HRU[2], zeros(length(Forest_Elevation_Count)) , Forest_Elevation_Count, length(Forest_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100), (Snow_Threshold, Height_Threshold), 0, [0], 0, [0],  0, 0)
                grass_input = HRU_Input(Area_Grass_Elevations, Current_Percentage_HRU[3], zeros(length(Grass_Elevation_Count)) , Grass_Elevation_Count,length(Grass_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100), (Snow_Threshold, Height_Threshold), 0, [0], 0, [0],  0, 0)
                rip_input = HRU_Input(Area_Rip_Elevations, Current_Percentage_HRU[4], zeros(length(Rip_Elevation_Count)) , Rip_Elevation_Count, length(Rip_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100),(Snow_Threshold, Height_Threshold), 0, [0], 0, [0],  0, 0)
                all_inputs = [bare_input, forest_input, grass_input, rip_input]
                #print(typeof(all_inputs))
                push!(Inputs_All_Zones, all_inputs)

                bare_storage = Storages(0, zeros(length(Bare_Elevation_Count)), zeros(length(Bare_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)
                forest_storage = Storages(0, zeros(length(Forest_Elevation_Count)), zeros(length(Forest_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)
                grass_storage = Storages(0, zeros(length(Grass_Elevation_Count)), zeros(length(Grass_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)
                rip_storage = Storages(0, zeros(length(Rip_Elevation_Count)), zeros(length(Rip_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)

                all_storages = [bare_storage, forest_storage, grass_storage, rip_storage]
                push!(Storages_All_Zones, all_storages)
        end
        # ---------------- CALCULATE OBSERVED OBJECTIVE FUNCTIONS -------------------------------------
        # calculate the sum of precipitation of all precipitation zones to calculate objective functions
        Total_Precipitation = Precipitation_All_Zones[1][:,1]*Area_Zones_Percent[1] + Precipitation_All_Zones[2][:,1]*Area_Zones_Percent[2]
        # end of spin up time is 3 years after the start of the calibration and start in the month October
        index_spinup = findfirst(x -> Dates.year(x) == firstyear + 2 && Dates.month(x) == 10, Timeseries)
        # evaluations chouls alsways contain whole year
        index_lastdate = findfirst(x -> Dates.year(x) == lastyear && Dates.month(x) == 10, Timeseries) - 1

        delete_days = readdlm(local_path*"HBVModel/Defreggental/Delete_Days.csv", ',', Int)
        Timeseries_Obj = Timeseries[index_spinup: index_lastdate]
        deleteat!(Timeseries_Obj, delete_days)
        Observed_Discharge_Obj = Observed_Discharge[index_spinup: index_lastdate] .* scale_factor_Discharge
        observed_AC_1day = autocorrelation(Observed_Discharge_Obj, 1)
        observed_AC_90day = autocorrelationcurve(Observed_Discharge_Obj, 90)[1]

        deleteat!(Observed_Discharge_Obj, delete_days)
        Total_Precipitation_Obj = Total_Precipitation[index_spinup: index_lastdate]
        deleteat!(Total_Precipitation_Obj, delete_days)
        #calculating the observed FDC; AC; Runoff
        observed_FDC = flowdurationcurve(log.(Observed_Discharge_Obj))[1]
        # observed_AC_1day = autocorrelation(Observed_Discharge_Obj, 1)
        # observed_AC_90day = autocorrelationcurve(Observed_Discharge_Obj, 90)[1]
        observed_monthly_runoff = monthlyrunoff(Area_Catchment, Total_Precipitation_Obj, Observed_Discharge_Obj, Timeseries_Obj)[1]
        # ---------------- START MONTE CARLO SAMPLING ------------------------
        #All_Goodness_new = []
        All_Goodness = zeros(29)
        #All_Parameter_Sets = Array{Any, 1}[]

        GWStorage = 55.0
        print("worker ", ID, " preparation finished", "\n")
        count = 1
        number_Files = 0
        #best_calibrations = readdlm("/home/sarah/Master/Thesis/Calibrations/Defreggental_less_dates/Defreggental_Parameterfit_All_runs_less_dates_best_500010.csv", ',')
        # parameters_best_calibrations = nmax[1+71430*(ID-1):71430*ID,:]
        # print("startvalue ", 1+71430*(ID-1), "endvalue", 71430*ID)

        #All_discharge = Array{Any, 1}[]
        for n in 1 : 1:nmax#size(parameters_best_calibrations)[1]
                Current_Inputs_All_Zones = deepcopy(Inputs_All_Zones)
                Current_Storages_All_Zones = deepcopy(Storages_All_Zones)
                Current_GWStorage = deepcopy(GWStorage)
                parameters, slow_parameters, parameters_array = parameter_selection_defreggental()

                #beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh = parameters_best_calibrations[n, :]
                # bare_parameters = Parameters(beta_Bare, Ce, 0, 0.0, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
                # forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
                # grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
                # rip_parameters = Parameters(beta_Rip, Ce, 0.0, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
                # slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)
                #
                # parameters = [bare_parameters, forest_parameters, grass_parameters, rip_parameters]
                # parameters_array = parameters_best_calibrations[n, :]

                Discharge, Snow_Extend = runmodelprecipitationzones_glacier(Potential_Evaporation, Glacier_All_Zones, Precipitation_All_Zones, Temperature_Elevation_Catchment, Current_Inputs_All_Zones, Current_Storages_All_Zones, Current_GWStorage, parameters, slow_parameters, Area_Zones, Area_Zones_Percent, Elevation_Percentage, Elevation_Zone_Catchment, ID_Prec_Zones, Nr_Elevationbands_All_Zones, observed_snow_cover, start2000)
                Discharge = Discharge * 1000 / Area_Catchment * (3600 * 24)
                # don't calculate the goodness of fit for the spinup time!
                Discharge_Obj = Discharge[index_spinup:index_lastdate]
                deleteat!(Discharge_Obj, delete_days)
                Goodness_Fit, ObjFunctions = objectivefunctions_delete_days(Discharge[index_spinup:index_lastdate], Discharge_Obj, Snow_Extend, Observed_Discharge_Obj, observed_FDC, observed_AC_1day, observed_AC_90day, observed_monthly_runoff, Area_Catchment, Total_Precipitation_Obj, Timeseries_Obj)
                #if goodness higher than -9999 save it
                #if goodness higher than -9999 save it
                if Goodness_Fit != -9999
                        Goodness = [Goodness_Fit, ObjFunctions, parameters_array]
                        Goodness = collect(Iterators.flatten(Goodness))
                        All_Goodness = hcat(All_Goodness, Goodness)
                        #append!(All_Goodness_new, Goodness_new)
                        #push!(All_Goodness, Goodness)
                        if size(All_Goodness)[2]-1 == 100
                                All_Goodness = transpose(All_Goodness[:, 2:end])
                                if count != 100
                                        open(local_path*"HBVModel/Defreggental_Parameterfit_less_dates_snow_redistr_"*string(ID)*"_"*string(number_Files)*".csv", "a") do io
                                                writedlm(io, All_Goodness,",")
                                        end
                                        count+= 1
                                else
                                        open(local_path*"HBVModel/Defreggental_Parameterfit_less_dates_snow_redistr_"*string(ID)*"_"*string(number_Files)*".csv", "a") do io
                                                writedlm(io, All_Goodness,",")
                                        end
                                        count = 1
                                        number_Files += 1
                                end

                                #print("worker ", ID, " wrote 100 tested parameter sets to file.", "\n")
                                All_Goodness = zeros(29)
                        end
                end
                if mod(n, 1000) == 0
                        print("number of runs", n, "\n")
                end
        end
        All_Goodness = transpose(All_Goodness[:, 2:end])
        open(local_path*"HBVModel/Defreggental_Parameterfit_less_dates_snow_redistr_"*string(ID)*".csv", "a") do io
                writedlm(io, All_Goodness,",")
        end
end
#nmax = readdlm("/home/sarah/Master/Thesis/Calibrations/Defreggental_less_dates/Defreggental_Parameterfit_All_runs_less_dates_best_500010.csv", ',')[:,10:29]
@time begin
#run_MC(1,100)
nmax = 10
pmap(ID -> run_MC(ID, nmax) , [1])
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
module SimpleHttpIO

export readlinebare,
       readlinestr,
       readlinebarestr,
       AbstractIOSocket,
       writeln,
       CRLF,
       AbstractIOSocket,
       IOSocket,
       readbyte,
       BYTE_TYPE,
       stop,
       LINE_READERS

const CRLF = "\r\n"
const N = Void
const N_TYPE = Type{Void}
const STR_TYPE = isdefined(:AbstractString)? AbstractString : String
const BYTE_TYPE = isdefined(:UInt8)? UInt8 : Uint8

abstract AbstractIOSocket

type IOSocket <: AbstractIOSocket
   io::IO
   env::Associative
end

IOSocket(io::IO; kwargs...) = IOSocket(io, Dict{Any, Any}(kwargs))

function writeln(io::AbstractIOSocket, args...)
    s = Base.write(io.io, args...)
    s2 = Base.write(io.io, CRLF.data)
    (isa(s, Integer) ? s : 0) + s2
end

Base.read{T}(io::AbstractIOSocket, t::Type{T}, bs::Int32) = Base.read(io.io, t , bs)

Base.readbytes(io::AbstractIOSocket, size::Integer=1024) = Base.read(io.io, BYTE_TYPE, size)
readbyte(io::AbstractIOSocket) = readbytes(io, 1)[1]
readline_(io::AbstractIOSocket) = begin
    buf = BYTE_TYPE[]
    last = N_TYPE

    while true
        b = read(io.io, BYTE_TYPE, 1)[1]
        push!(buf, b)

        if last == '\r' && b == '\n'
           break
        end

        last = b
    end
    buf
end

readlinestr(io::AbstractIOSocket) = UTF8String(readline_(io))
readlinebare(io::AbstractIOSocket) = readline_(io)[1:end-2]
Base.readline(io::AbstractIOSocket) = readline_(io)

readlinebarestr(io::AbstractIOSocket) = UTF8String(readlinebare(io))

const LINE_READERS = [:s => readlinestr, :s! => readlinebarestr,
                      :b => readline_, :b! => readlinebare]

type EachLineStopIterator <: Exception end

type EachLine
    stream::AbstractIOSocket
    reader::Function
    count::Integer
end

stop() = throw(EachLineStopIterator())

Base.eachline(stream::AbstractIOSocket, reader::Function=Base.readline) = EachLine(stream, reader, 0)
Base.eachline(stream::AbstractIOSocket, reader::Symbol=:b) = EachLine(stream, LINE_READERS[reader], 0)

Base.eachline(cb::Function, stream::AbstractIOSocket, reader::Function=Base.readline) = begin
    itr = Base.eachline(stream, reader)

    try
        for line in itr
            cb(itr, line)
        end
    catch e
        if !isa(e, EachLineStopIterator)
            rethrow()
        end
    end

    itr.count
end

Base.eachline(cb::Function, stream::AbstractIOSocket, reader::Symbol) = begin
    Base.eachline(cb, stream, LINE_READERS[reader])
end

Base.start(itr::EachLine) = (itr.count = 0; N)
Base.done(itr::EachLine, nada) = false

Base.next(itr::EachLine, nada) = begin
    itr.count += 1
    itr.reader(itr.stream), N
end

end


text:
# Returns the occupied orbitals, the occupation array and optionally the eigenvalues without
# virtual states (or states with small occupation level for metals).
# threshold is a parameter to distinguish between states we want to keep and the
# others when using temperature. It is set to 0.0 by default, to treat with insulators.
function select_occupied_orbitals(basis, ψ, occupation, eigenvalues=nothing; threshold=0.0)
    N = [something(findlast(x -> x > threshold, occk), 0) for occk in occupation]
    selected_ψ   = [@view ψk[:, 1:N[ik]] for (ik, ψk)   in enumerate(ψ)]
    selected_occ = [      occk[1:N[ik]]  for (ik, occk) in enumerate(occupation)]

    # if we have an insulator, sanity check that the orbitals we kept are the
    # occupied ones
    if threshold == 0.0
        model   = basis.model
        n_spin  = model.n_spin_components
        n_bands = div(model.n_electrons, n_spin * filled_occupation(model), RoundUp)
        @assert n_bands == size(selected_ψ[1], 2)
    end

    if isnothing(eigenvalues)
        (ψ=selected_ψ, occupation=selected_occ)
    else
        selected_evals = [evalk[1:N[ik]] for (ik, evalk) in enumerate(eigenvalues)]
        (ψ=selected_ψ, occupation=selected_occ, eigenvalues=selected_evals)
    end
end

function select_occupied_orbitals(scfres::NamedTuple; threshold=0.0)
    truncated = select_occupied_orbitals(scfres.basis, scfres.ψ, scfres.occupation,
                                         scfres.eigenvalues; threshold)

    min_trunc_bands = minimum(length, scfres.occupation) - maximum(length, truncated.occupation)
    @assert min_trunc_bands ≥ scfres.n_ep_extra  # Ensure extra bands are truncated
    merge(scfres, truncated, (; n_ep_extra=0))
end

# Packing routines used in direct_minimization and newton algorithms.
# They pack / unpack sets of ψ's (or compatible arrays, such as hamiltonian
# applies and gradients) to make them compatible to be used in algorithms
# from IterativeSolvers.
# Some care is needed here : some operators (for instance K in newton.jl)
# are real-linear but not complex-linear. To overcome this difficulty, instead of
# seeing them as operators from C^N to C^N, we see them as
# operators from R^2N to R^2N. In practice, this is done with the
# reinterpret function from julia.
# /!\ pack_ψ does not share memory while unpack_ψ does

reinterpret_real(x) = reinterpret(real(eltype(x)), x)
reinterpret_complex(x) = reinterpret(Complex{eltype(x)}, x)

function pack_ψ(ψ)
    # TODO as an optimization, do that lazily? See LazyArrays
    vcat([vec(ψk) for ψk in ψ]...)
end

# Returns pointers into the unpacked ψ
# /!\ The resulting array is only valid as long as the original x is still in live memory.
function unsafe_unpack_ψ(x, sizes_ψ)
    lengths = prod.(sizes_ψ)
    ends = cumsum(lengths)
    # We unsafe_wrap the resulting array to avoid a complicated type for ψ.
    map(1:length(sizes_ψ)) do ik
        unsafe_wrap(Array{complex(eltype(x))},
                    pointer(@views x[ends[ik]-lengths[ik]+1:ends[ik]]),
                    sizes_ψ[ik])
    end
end
unpack_ψ(x, sizes_ψ) = deepcopy(unsafe_unpack_ψ(x, sizes_ψ))

function random_orbitals(basis::PlaneWaveBasis{T}, kpt::Kpoint, howmany) where {T}
    ortho_qr(randn(Complex{T}, length(G_vectors(basis, kpt)), howmany))
end


text:
#InspectDR: Defaults
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
#Identifies whether a plot should apply f1 acceleration to drop points:
const PDM_NEVER = PointDropMatrix([false false; false false])
const PDM_ALWAYS = PointDropMatrix([true true; true true])
const PDM_NOGLYPH = PointDropMatrix([false false; true false])
const PDM_HASLINE = PointDropMatrix([false false; true true])
const PDM_DEFAULTS = Dict(
	:never => PDM_NEVER,
	:always => PDM_ALWAYS,
	:noglyph => PDM_NOGLYPH,
	:hasline => PDM_HASLINE,
)

#==Types
===============================================================================#
mutable struct Defaults
	rendersvg::Bool #Might want to dissalow SVG renderings for performance reasons
	xaxiscontrol_visible::Bool
	pointdropmatrix::PointDropMatrix
	colorscale::ColorScale

	plotlayout::PlotLayout
	mplotlayout::MultiplotLayout
end
Defaults() = Defaults(
	false, false, PDM_NEVER, ColorScale(),
	PlotLayout(PREDEFAULTS), MultiplotLayout(PREDEFAULTS)
)

#==Data
===============================================================================#
const global defaults = Defaults()

#==Constructors
===============================================================================#

PlotLayout() = PlotLayout(defaults.plotlayout)
MultiplotLayout() = MultiplotLayout(defaults.mplotlayout)

#==Set style of Defaults
===============================================================================#
function setstyle!(d::Defaults, styleid::Symbol; refresh::Bool=true, kwargs...)
	s.plotlayout = getstyle(PlotLayout, StyleID(styleid); kwargs...)
	s.mplotlayout = getstyle(MultiplotLayout, StyleID(styleid); kwargs...)
end

#==Data Initialization
===============================================================================#
#Initialize InspectDR.defaults (To be called in __init__()):
function _initialize(dflt::Defaults)
	local userdefaults = Dict{Symbol, Any}()
	try
		userdefaults = copy(Main.DEFAULTS_INSPECTDR)
	catch
	end

	#Get dict value, or default if not exist.
	#Note: Pop key from dict to avoid double processing
	function condget(dict, key::Symbol, T::Type, default)
		if haskey(dict, key)
			return T(pop!(dict, key))
		else
			return T(default)
		end
	end
	function condget(dict, key::Symbol, T::Type, default::KeepPrevType)
		if haskey(dict, key)
			return T(pop!(dict, key))
		else
			return default
		end
	end

	dflt.rendersvg = condget(userdefaults, :rendersvg, Bool, true)
	dflt.xaxiscontrol_visible = condget(userdefaults, :xaxiscontrol_visible, Bool, true)
	notation_x = condget(userdefaults, :notation_x, Symbol, :ENG)
	notation_y = condget(userdefaults, :notation_y, Symbol, :ENG)
	notation_z = condget(userdefaults, :notation_z, Symbol, :ENG)
	droppoints = condget(userdefaults, :droppoints, Symbol, :noglyph)
	dflt.pointdropmatrix = PDM_DEFAULTS[droppoints]
	cmapid = condget(userdefaults, :colorscale, String, "Oranges")
	dflt.colorscale = ColorScale(Colors.colormap(cmapid))

	fontname = condget(userdefaults, :fontname, String, DEFAULT_FONTNAME)
	fontscale = condget(userdefaults, :fontscale, Float64, 1.0)

	dflt.mplotlayout = getstyle(MultiplotLayout, :screen)
	overwritefont!(dflt.mplotlayout, fontname=fontname, fontscale=fontscale)

	#Manually get SUPPORTED defaults for Multiplot layout (seed from :screen):
	condoverwrite(key::Symbol, T::Type) =
		(overwriteprop!(dflt.mplotlayout, key, condget(userdefaults, key, T, KEEP_PREV)))
	condoverwrite(:ncolumns, Int)
	condoverwrite(:valloc_plot, Float64)
	condoverwrite(:halloc_plot, Float64)

	#Automatically get defaults for Plot layout (seed from :screen):
	dflt.plotlayout = getstyle(PlotLayout, :screen,
		notation_x=notation_x, notation_y=notation_x, enable_legend=false
	)
	overwritefont!(dflt.plotlayout, fontname=fontname, fontscale=fontscale)
	for (k, v) in userdefaults
		try
			dflt.plotlayout[k] = v
		catch
			@warn("Cannot set default PlotLayout.$k = $v.")
		end
	end

	return
end

#Last line


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "msolve"
version = v"0.1.3"

# Collection of sources required to complete build
sources = [
    GitSource("https://gitlab.lip6.fr/safey/msolve.git", "0e6de1bab447381b7fdf42c3c414caa6ed2a8f88"),
    #= ArchiveSource("https://www.mathematik.uni-kl.de/~ederc/msolve-0.1.2.tar.gz", "ce6454b28477cb3b5670042faf7b3282e234fe1e0ee5a62c184d0512ef4126e1") =#
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd msolve/
./autogen.sh
ac_cv_func_malloc_0_nonnull=yes ac_cv_func_realloc_0_nonnull=yes ./configure --with-gnu-ld --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-openmp
make -j${nprocs}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc"),
    Platform("x86_64", "macos"; )
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libmsolve", :libmsolve),
    LibraryProduct("libneogb", :libneogb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.1.2"; compat="6.1.2")
    Dependency(PackageSpec(name="FLINT_jll"), compat = "~200.800")
    Dependency(PackageSpec(name="MPFR_jll", uuid="3a97d323-0669-5f0c-9066-3539efd106a3"))
]

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

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil