text:
###############################################################################
####################         JLD2Saver Callback            ####################
###############################################################################

"""
A callback to save the TMLE estimation results to disk instead of keeping them in memory
"""
mutable struct JLD2Saver <: TMLE.Callback
    file::String
    save_machines::Bool
end

function TMLE.after_tmle(callback::JLD2Saver, report::TMLEReport, target_id::Int, query_id::Int)
    jldopen(callback.file, "a") do io
        group = haskey(io, "TMLEREPORTS") ? io["TMLEREPORTS"] : JLD2.Group(io, "TMLEREPORTS")
        group[string(target_id, "_", query_id)] = report
    end
end

function TMLE.after_fit(callback::JLD2Saver, mach::Machine, id::Symbol)
    if callback.save_machines
        jldopen(callback.file, "a") do io
            group = haskey(io, "MACHINES") ? io["MACHINES"] : JLD2.Group(io, "MACHINES")
            group[string(id)] = mach
        end
    end
end

function TMLE.finalize(callback::JLD2Saver, estimation_report::NamedTuple)
    jldopen(callback.file, "a") do io
        io["low_propensity_scores"] = estimation_report.low_propensity_scores
    end
    return estimation_report
end

###############################################################################
####################             Utilities                 ####################
###############################################################################

function sample_ids_per_phenotype(data, phenotypes_columns)
    sample_ids = Dict()
    for colname in phenotypes_columns
        sample_ids[colname] =
            dropmissing(data[!, ["SAMPLE_ID", colname]]).SAMPLE_ID

    end
    return sample_ids
end

convert_target(x, ::Type{Bool}) = categorical(x)
convert_target(x, ::Type{Real}) = x

function genotypes_combinations(query)
    snps = keys(query.case)
    snp_levels = NamedTuple{snps}(
        [[getfield(query.control, snp), getfield(query.case, snp)] for snp in snps]
    )
    combinations = DataFrame(NamedTuple{snps}([[] for _ in snps]))
    for comb in reduce(vcat, Iterators.product(snp_levels...))
        push!(combinations, comb)
    end
    return combinations
end

function actualqueries(genotypes, queries)
    unique_genotypes = unique(dropmissing(genotypes))
    actual_queries = []
    for query in queries
        required_genotypes = genotypes_combinations(query)
        joined = innerjoin(unique_genotypes, required_genotypes, on=names(genotypes))
        if size(joined, 1) == size(required_genotypes, 1)
            push!(actual_queries, query)
        else
            @warn "Query: ", query.name, " will not be processed due to missing genotypic data."
        end
    end
    return actual_queries
end

function preprocess(genotypes, confounders, phenotypes, target_type, queries)
    # Make sure data SAMPLE_ID types coincide
    genotypes.SAMPLE_ID = string.(genotypes.SAMPLE_ID)
    confounders.SAMPLE_ID = string.(confounders.SAMPLE_ID)
    phenotypes.SAMPLE_ID = string.(phenotypes.SAMPLE_ID)

    # columns names
    genotypes_columns = filter(!=("SAMPLE_ID"), names(genotypes))
    confounders_columns = filter(!=("SAMPLE_ID"), names(confounders))
    phenotypes_columns = filter(!=("SAMPLE_ID"), names(phenotypes))

    # Join all elements together by SAMPLE_ID
    data = innerjoin(
            innerjoin(genotypes, confounders, on=:SAMPLE_ID),
            phenotypes,
            on=:SAMPLE_ID
            )

    # Drop missing values based on genotypes and covariates
    dropmissing!(data, vcat(genotypes_columns, confounders_columns))

    # Retrieve T and convert to categorical data
    # The use of the query he is so that the order of the columns in both
    # T and queries are matching which is currently a technical requirement of the TMLE package
    T = DataFrame()
    for name in keys(first(queries).control)
        T[:, name] = categorical(data[!, name])
    end

    # Retrive sample_ids per phenotype
    sample_ids = sample_ids_per_phenotype(data, phenotypes_columns)

    # Retrieve W
    W = DataFrame()
    for col in confounders_columns
        W[:, col] = convert(Vector{Float64}, data[!, col])
    end

    # Retrieve Y and convert to categorical if needed
    Y = DataFrame()
    for name in phenotypes_columns
        Y[:, name] = convert_target(data[!, name], target_type)
    end

    return T, W, Y, sample_ids
end

###############################################################################
####################           Main Function               ####################
###############################################################################

function UKBBVariantRun(parsed_args)
    v = parsed_args["verbosity"]
    target_type = parsed_args["target-type"] == "Real" ? Real : Bool

    queries = parse_queries(parsed_args["queries"])

    v >= 1 && @info "Loading data."
    # Load Genotypes
    genotypes = CSV.read(
        parsed_args["genotypes"],
        DataFrame,
        select=[:SAMPLE_ID, keys(first(queries).control)...]
    )

    # Read Confounders
    confounders = CSV.File(parsed_args["confounders"]) |> DataFrame

    # Load phenotypes
    phenotypes = load_phenotypes(parsed_args["phenotypes"], parsed_args["phenotypes-list"])

    # Preprocessing
    v >= 1 && @info "Preprocessing data."
    genotypes, confounders, phenotypes, sample_ids =
        preprocess(genotypes, confounders, phenotypes, target_type, queries)

    # Each genotype should have probability > 0
    queries = actualqueries(genotypes, queries)
    if size(queries, 1) == 0
        jldopen(parsed_args["out"], "a") do io
            JLD2.Group(io, "TMLEREPORTS")
        end
        return 0
    end

    # Build estimator
    tmle_config = TOML.parsefile(parsed_args["estimator"])
    tmle = estimator_from_toml(tmle_config, queries, target_type, adaptive_cv=parsed_args["adaptive-cv"])

    # Run TMLE
    v >= 1 && @info "TMLE Estimation."
    TMLE.fit(tmle, genotypes, confounders, phenotypes;
             verbosity=v-1,
             cache=false,
             callbacks=[JLD2Saver(parsed_args["out"], parsed_args["save-full"])]
    )

    # Save sample_ids
    jldopen(parsed_args["out"], "a") do io
        io["SAMPLE_IDS"] = sample_ids
    end

    v >= 1 && @info "Done."
    return 0
end


text:
function npart_plot(sim; species="electron")
    ts = get_time.(sim) .|> unit_t
    npart = get_npart.(sim, (species,))

    Plots.plot(ts, npart, xlabel="t", ylabel="N")
end

ext: jl


text:
cd("/home/albyback/git/dotlabelproject/OuterProject/SparseInverseProblems.jl/")

# this uses a forked version
# git clone https://github.com/marius311/SparseInverseProblems.jl.git
# which includes a toml file and some
# modifications for this to run in Julia properly

using Pkg
Pkg.activate(".")
using SparseInverseProblems
using Compose
# cd("/home/albyback/git/dotlabelproject/OuterProject/SparseInverseProblems.jl/examples/simple_example/")
include("/home/albyback/git/dotlabelproject/OuterProject/SparseInverseProblems.jl/examples/simple_example/simple_example.jl")

#Feel free to override
function localUpdate(sim :: ForwardModel,lossFn :: Loss,
    thetas :: Matrix{Float64}, y :: Vector{Float64}, tau :: Float64, max_iters)
  for cd_iter = 1:max_iters
    weights = solveFiniteDimProblem(sim, lossFn, thetas, y, tau)
    #remove points with zero weight
    if any(weights.==0.0)
      println("Removing ",sum(weights.==0.0), " zero-weight points.")
      thetas = thetas[:,weights.!= 0.0]
      weights = weights[weights.!= 0.0]
    end
    #local minimization over the support
    new_thetas = localDescent(sim, lossFn, thetas,weights, y)
    #break if termination condition is met
    if length(thetas) == length(new_thetas) && maximum(abs.(vec(thetas)-vec(new_thetas))) <= 1E-7
        break
    end
    thetas = new_thetas
  end
  #final prune
  weights = solveFiniteDimProblem(sim, lossFn, thetas, y, tau)
  if any(weights.==0.0)
    println("Removing ",sum(weights.==0.0), " zero-weight points.")
    thetas = thetas[:,weights.!= 0.0]
    weights = weights[weights.!= 0.0]
  end
  return thetas, weights
end

function run_sample()
      evaluation_points = -5.0:0.25:5.0
      model = SimpleExample(evaluation_points, -10.0:0.5:10.0)
      means = [-2.0,1.0,3.5]
      k = length(means)
      weights = ones(k)+randn(k)*0.3
      weights[1] = 0.2
      target = max.(phi(model,Matrix(means'),weights) + randn(length(evaluation_points))*0.1,0.0);
      # (means_est,weights_est) = ADCG(model, LSLoss(), target, 100.0;  callback=callback)

      sim = model
      lossFn =LSLoss()
      y = target
      tau = 100.0;
      max_iters = 50
      min_optimality_gap = 1E-5
      max_cd_iters = 200
      fully_corrective = false
      @assert(tau > 0.0)

      bound = -Inf
      thetas = zeros(0,0) #hack
      weights = zeros(0)
      #cache the forward model applied to the current measure.
      output = zeros(length(y))

      iter = 1

      #compute the current residual
      residual = output - y
      #evalute the objective value and gradient of the loss
      objective_value, grad = loss(lossFn, residual)
      #compute the next parameter value to add to the support
      theta,score = lmo(sim,grad)
      print(theta)
      print(score)

      old_thetas = thetas
      thetas = iter == 1 ? reshape(theta, length(theta),1) : [thetas theta]
      #run local optimization over the support.
      old_weights = copy(weights)
      thetas,weights = localUpdate(sim,lossFn,thetas,y,tau,max_cd_iters)

      #score is - |<\psi(theta), gradient>|
      #update the lower bound on the optimal value
      # bound = max(bound, objective_value+score*tau-dot(output,grad))
      #check if the bound is met.
      # if (objective_value < bound + min_optimality_gap || score >= 0.0)
      # return thetas,weights
      # end
end

run_sample()

using Plots

plot(evaluation_points,y)

# # want to run a 1D example
# abstract BoxConstrainedDifferentiableModel <: ForwardModel
#
# #compute the measurement model
# psi(model :: BoxConstrainedDifferentiableModel, theta :: Vector{Float64})
#
# #compute the jacobian of the forward model
# dpsi(model :: BoxConstrainedDifferentiableModel, theta :: Vector{Float64})
#
# # Initial starting point for continuous optimization for the FW step.
# # Should return a good guess for $\arg\min_\theta \langle \psi(theta), v \rangle.$
# # Often computed using a grid.
# getStartingPoint(model :: BoxConstrainedDifferentiableModel, v :: Vector{Float64}) =
#   error("getStartingPoint not implemented for model $(typeof(model)).")
#
# # Box constraints on the parameters.
# # Returns a tuple of two vectors : lower bounds and upper bounds.
# parameterBounds(model :: BoxConstrainedDifferentiableModel)


text:
## big.jl  A test with a larger set of data
## (c) 2014--2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

using Test
import Random

## Note that due to the random sample, some of these test may fail every once in a while,
## So let's set a seed
Random.seed!(1)

## test some statistics on large data sets
## calibrated scores (i.e., these scores can be interpreted as log-likelihood-ratios)
tar =  2 .+ 2randn(1_000_000)
non = -2 .+ 2randn(1_000_000)

## eer
@test eer(tar, non) ≈ pnorm(-1) atol=0.01

## eerch
@time r = roc(tar, non)
@test eerch(r) ≈ pnorm(-1) atol=0.01

## AUC
pauc = 0.0786496
@test auc(r) ≈ pauc atol=0.01
auc(r, pfa=0.1)
AUC(r, pfa=0.1)

## Cllr
cllr2 = 0.5140558 ## Cllr for distributions with μ = ±2 and σ = 2, i.e., d'=2
@test cllr(tar, non) ≈ cllr2 atol=0.01
@test mincllr(tar, non) ≈ cllr2 atol=0.01


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule jlqml_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("jlqml")
JLLWrappers.@generate_main_file("jlqml", UUID("6b5019fb-a83d-5b4e-a9f7-678a36c28df7"))
end  # module jlqml_jll


text:
# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::CurvedMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype<:Real}
  elements = init_elements(mesh, equations, dg.basis, uEltype)

  cache = (; elements)

  return cache
end

# Extract contravariant vector Ja^i (i = index) as SVector
@inline function get_contravariant_vector(index, contravariant_vectors, indices...)

  SVector(ntuple(dim -> contravariant_vectors[index, dim, indices...], ndims(contravariant_vectors) - 3))
end

@inline ndofs(mesh::CurvedMesh, dg::DG, cache) = nelements(cache.elements) * nnodes(dg)^ndims(mesh)

include("containers.jl")
include("dg_1d.jl")
include("dg_2d.jl")
include("dg_3d.jl")


text:
# This file contains functions for the Cost Efficiency Reverse DDF DEA model
"""
CostReverseDDFDEAModel
An data structure representing a cost RDDF DEA model.
"""
struct CostReverseDDFDEAModel <: AbstractCostDEAModel
    n::Int64
    m::Int64
    s::Int64
    measure::Symbol
    rts::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
    Gxrddf::Union{Vector,Matrix}
    Gyrddf::Union{Vector,Matrix}
end

"""
    deacostrddf(X, Y, W, measure)
Compute profit efficiency using data envelopment analysis Reverse DDF model for
inputs `X`, outputs `Y`, price of inputs `W`, and efficiency `measure`.

# Measure specification:
- `:ERG`: Enhanced Russell Graph Slack Based Measure.

# Optional Arguments
- `rts=:VRS`: choose between constant returns to scale `:CRS` or variable returns to scale `:VRS`.
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];

julia> Y = [1; 1; 1; 1; 1; 1; 1; 1];

julia> W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

julia> deacostrddf(X, Y, W, :ERG)
Cost Reverse DDF DEA Model
DMUs = 8; Inputs = 2; Outputs = 1
Returns to Scale: VRS
Associated efficiency measure = ERG
──────────────────────────────────
       Cost  Technical  Allocative
──────────────────────────────────
1  0.0        0.0        0.0
2  0.5        0.0        0.5
3  0.5        0.0        0.5
4  0.416667   0.416667   0.0
5  0.6        0.6        0.0
6  0.25       0.166667   0.0833333
7  0.525      0.35       0.175
8  0.532609   0.4375     0.0951087
──────────────────────────────────
```
"""
function deacostrddf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector},
    measure::Symbol;
    rts::Symbol = :VRS, monetary::Bool = false, atol::Float64 = 1e-6,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::CostReverseDDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end

    n = nx

    # Get directions based on efficiency measure
    Gxrddf, Gyrddf = DataEnvelopmentAnalysis.dearddfdirections(X, Y, measure, orient = :Input, rts = rts, optimizer = optimizer)

    # Get normalization of efficiency measure model
    if measure == :ERG
        measurecostmodel = deacostrussell(X, Y, W, rts = rts, optimizer = optimizer)
    end

    normmeasurecost = normfactor(measurecostmodel)

    # Cost model using calculated directions
    rddfcostmodel = deacostddf(X, Y, W, Gx = Gxrddf, rts = rts, optimizer = optimizer)

    cefficiency = efficiency(rddfcostmodel)
    clambdaeff = peersmatrix(rddfcostmodel)
    techefficiency = efficiency(rddfcostmodel, :Technical)
    allocefficiency = efficiency(rddfcostmodel, :Allocative)
    normalization = normfactor(rddfcostmodel)
    Xtarget = targets(rddfcostmodel, :X)
    Ytarget = targets(rddfcostmodel, :Y)

    # Change normalization factor of efficient DMUs
    for i=1:n
        if abs(techefficiency[i] ) <= atol
            cefficiency[i] = cefficiency[i] ./ normmeasurecost[i] .* normalization[i]
            normalization[i] = normmeasurecost[i]
        end
    end

    if monetary
        cefficiency = cefficiency .* normalization
        techefficiency = techefficiency .* normalization
    end

    allocefficiency = cefficiency - techefficiency

    return CostReverseDDFDEAModel(n, m, s, measure, rts, monetary, names, cefficiency, clambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget, Gxrddf, Gyrddf)

end

function Base.show(io::IO, x::CostReverseDDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Cost Reverse DDF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale: ", string(x.rts))
        print(io, "\n")
        print(io, "Associated efficiency measure = ", string(x.measure))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Cost", "Technical", "Allocative"], dmunames))
    end

end


text:
## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-audio-analysis

"""
    tracks_get_audio_analysis(track_id)

**Summary**: Get a detailed audio analysis for a single track identified by it's unique Spotify ID.

# Arguments
- `track_id`: The Spotify ID for the track

# Example
```julia-repl
julia> Spotify.tracks_get_audio_analysis("6rqhFgbbKwnb9MLmUQDhG6")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :meta     => {…
  :track    => {…
  :bars     => JSON3.Object[{…
  :beats    => JSON3.Object[{…
```
"""
function tracks_get_audio_analysis(track_id)
    return spotify_request("audio-analysis/$track_id")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-audio-features

"""
    tracks_get_audio_features(track_id)

**Summary**: Get audio feature information for a single track identified by it's unique Spotify ID.

# Arguments
- `track_id`: The Spotify ID for the track

# Example
```julia-repl
julia> Spotify.tracks_get_audio_features("6rqhFgbbKwnb9MLmUQDhG6")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 18 entries:
  :danceability     => 0.592
  :energy           => 0.0196
  :key              => 1
  :loudness         => -33.35
  :mode             => 1
  :speechiness      => 0.0358
  :acousticness     => 0.362
  :instrumentalness => 0.854
```
"""
function tracks_get_audio_features(track_id)
    return spotify_request("audio-features/$track_id")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-track

"""
    tracks_get(track_id; market="US")

**Summary**: Get a spotify catalog information for a single track identified by it's unique Spotify ID.

# Arguments
- `track_id` : The Spotify ID for the track.

# Optional keywords
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".

# Example
```julia-repl
julia> Spotify.tracks_get("6rqhFgbbKwnb9MLmUQDhG6")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 18 entries:
  :album         => {…
  :artists       => JSON3.Object[{…
  :disc_number   => 1
  :duration_ms   => 65314
  :explicit      => true
  :external_ids  => {…
  :external_urls => {…
```
"""
function tracks_get_single(track_id; market::String="US")
    return spotify_request("tracks/$track_id?market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-several-tracks

"""
    tracks_get_multiple(ids::String; market::String="US")

**Summary**: Get Spotify catalog information for multiple tracks based on their Spotify IDs.

# Arguments
- `ids` : A comma-separated list of the Spotify IDs.

# Optional keywords
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".

# Example
```julia-repl
julia> Spotify.tracks_get_multiple("4iV5W9uYEdYUVa79Axb7Rh,1301WleyT98MSxVHPZCA6")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 1 entry:
  :tracks => Union{Nothing, JSON3.Object}[{…
```
"""
function tracks_get_multiple(ids; market::String="US")
    return spotify_request("tracks?ids=$ids&market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-users-saved-tracks

"""
    tracks_get_saved(;limit::Int64=20, market::String="US", offset::Int64=0)

**Summary**: Get a list of the songs saved in the current Spotify user's 'Your Music' library.

# Optional keywords
- `limit::Int64` : Maximum number of items to return, default is set to 20
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".
- `offset::Int64` : Index of the first item to return, default is set to 0

# Example
```julia-repl
julia> Spotify.tracks_get_saved()[1]
[ Info: We try the request without checking if current grant includes scope user-library-read.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :href     => "https://api.spotify.com/v1/me/tracks?offset=0&limit=20&market=US"
  :items    => JSON3.Object[{…
  :limit    => 20
  :next     => "https://api.spotify.com/v1/me/tracks?offset=20&limit=20&market=US"
```
"""
function tracks_get_saved(;limit::Int64=20, market::String="US", offset::Int64=0)
    return spotify_request("me/tracks?limit=$limit&market=$market&offset=$offset";
                           scope = "user-library-read")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-users-saved-tracks

"""
    tracks_get_contains(ids)

**Summary**: Check if one or more tracks is already saved in the current Spotify user's 'Your Music' library.

# Arguments
- `ids` : A comma-separated list of the Spotify track IDs.

# Example
```julia-repl
julia> Spotify.tracks_get_contains("4iV5W9uYEdYUVa79Axb7Rh, 4VqPOruhp5EdPBeR92t6lQ")[1]
2-element JSON3.Array{Bool, Base.CodeUnits{UInt8, String}, Vector{UInt64}}:
 0
 0
```
"""
function tracks_get_contains(ids)
    return spotify_request("me/tracks/contains?ids=$ids"; scope = "user-library-read")
end


text:
using Bridge
 using Random
 using DataFrames
 using CSV

 SRC_DIR = joinpath(Base.source_dir(), "..", "src")
 include(joinpath(SRC_DIR, "auxiliary", "data_simulation_fns.jl"))
 OUT_DIR = joinpath(Base.source_dir(), "..", "output")
 mkpath(OUT_DIR)

 include(joinpath(SRC_DIR, "examples", "fitzHughNagumo.jl"))
 param = :simpleConjug
 FILENAME_OUT = joinpath(OUT_DIR,
                         "test_path_part_obs_"*String(param)*".csv")

 P = FitzhughDiffusion(param, 10.0, -8.0, 15.0, 0.0, 3.0)
 # starting point under :regular parametrisation
 x0 = ℝ{2}(-0.5, 0.6)
 # translate to conjugate parametrisation
 x0 = regularToConjug(x0, P.ϵ, 0.0)

 dt = 1/50000
 T = 10.0
 tt = 0.0:dt:T

 Random.seed!(4)
 XX, _ = simulateSegment(0.0, x0, P, tt)

 num_obs = 100
 skip = div(length(tt), num_obs)
 Time = collect(tt)[1:skip:end]
 df = DataFrame(time=Time, x1=[x[1] for x in XX.yy[1:skip:end]],
                x2=[(i==1 ? x0[2] : NaN) for (i,t) in enumerate(Time)])
 CSV.write(FILENAME_OUT, df)


text:
using Test
import DataFrames: DataFrame, DataFrameRow

using JMcDM

const testUtilityFunctions = true
const testMCDMFunctions = true
const testSCDMFunctions = true
const testLPBasedFunctions = true

if testUtilityFunctions
    include("./testutility.jl")
end

if testSCDMFunctions
    include("./testscdm.jl")
end # Test SCDM Tools

if testMCDMFunctions
    include("./testmcdm.jl")
end  # Test MCDM Tools

if testLPBasedFunctions
    include("./testlp.jl")
end # Test LP Based Tools

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
# It'st most appropriate to define a Julia iterable object for this task
# Julia doesn't have Python'st yield, the closest to it is produce/consume calls with Julia tasks
# but for various reasons they don't work out for this task
# This solution works with two integers, a Julia rational or a real

mutable struct ContinuedFraction{T<:Integer}
    n1::T # numerator or real
    n2::T # denominator or 1 if real
    t1::T # generated coefficient
end

# Constructors for all possible input types
ContinuedFraction{T<:Integer}(n1::T, n2::T) = ContinuedFraction(n1, n2, 0)
ContinuedFraction(n::Rational) = ContinuedFraction(numerator(n), denominator(n))
ContinuedFraction(n::AbstractFloat) = ContinuedFraction(Rational(n))

# Methods to make our object iterable
Base.start(::ContinuedFraction) = nothing
# Returns true if we've prepared the continued fraction
Base.done(cf::ContinuedFraction, st) = cf.n2 == 0
# Generates the next coefficient
function Base.next(cf::ContinuedFraction, st)
    cf.n1, (cf.t1, cf.n2) = cf.n2, divrem(cf.n1, cf.n2)
    return cf.t1, nothing
end

# Tell Julia that this object always returns ints (all coeffs are integers)
Base.eltype{T}(::Type{ContinuedFraction{T}}) = T

# Overload the default collect function so that we can collect the first maxiter coeffs of infinite continued fractions
# array slicing doesn't work as Julia crashes before the slicing due to our infinitely long array
function Base.collect(itr::ContinuedFraction, maxiter::Integer = 100)
    r = Array{eltype(itr)}(maxiter)
    i = 1
    for v in itr
        r[i] = v
        i += 1
        if i > maxiter break end
    end
    return r[1:i-1]
end

# Test cases according to task description with outputs in comments
println(collect(ContinuedFraction(1, 2)))       # => [0, 2]
println(collect(ContinuedFraction(3, 1)))       # => [3]
println(collect(ContinuedFraction(23, 8)))      # => [2, 1, 7]
println(collect(ContinuedFraction(13, 11)))     # => [1, 5, 2]
println(collect(ContinuedFraction(22, 7)))      # => [3, 7]
println(collect(ContinuedFraction(14142, 10000)))       # => [1, 2, 2, 2, 2, 2, 1, 1, 29]
println(collect(ContinuedFraction(141421, 100000)))     # => [1, 2, 2, 2, 2, 2, 2, 3, 1, 1, 3, 1, 7, 2]
println(collect(ContinuedFraction(1414214, 1000000)))   # => [1, 2, 2, 2, 2, 2, 2, 2, 3, 6, 1, 2, 1, 12]
println(collect(ContinuedFraction(14142136, 10000000))) # => [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 6, 1, 2, 4, 1, 1, 2]

println(collect(ContinuedFraction(13 // 11)))   # => [1, 5, 2]
println(collect(ContinuedFraction(√2), 20))     # => [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# using Plots: plot, plot!, heatmap, heatmap!, @layout, cgrad
# using Plots.PlotMeasures
# # using Plots.grid
# using Plots: plot, plot!, heatmap, @layout, cgrad, grid, heatmap!
# using AbstractPlotting, # GLMakie
# using UnicodePlots

export plot_ε, compare_fields, plot_field, plot_d⃗, ylims

################################################################################
#                              Utility Functions                               #
################################################################################

# ε plotting functions

function plot_ε(ε1,x::AbstractVector{Float64},y::AbstractVector{Float64}; outlines=false,cmap_diag=cgrad(:viridis),cmap_offdiag=cgrad(:bluesreds),scale_each_component=true)
	tmi,bmi,lmi,rmi = -4mm, -4mm, -4mm, -4mm
	tmo,bmo,lmo,rmo = 2mm, 0mm, 0mm, 0mm

	labels = [	"ε₁₁"	"ε₁₂"	"ε₁₃"
				"ε₂₁"	"ε₂₂"	"ε₂₃"
				"ε₃₁"	"ε₃₂"	"ε₃₃"	]

	showaxes = [
			:y 	:no :no
			:y 	:no :no
		 	:yes :x :x
	]
	lm = [
			lmo lmi lmi
			lmo lmi lmi
			lmo lmi lmi
	]
	rm = [
			rmi rmi rmo
			rmi rmi rmo
			rmi rmi rmo
	]
	tm = [
			tmo tmo tmo
			tmi tmi tmi
			tmi tmi tmi
	]
	bm = [
			bmi bmi bmi
			bmi bmi bmi
			bmo bmo bmo
	]

	plot_grid =[
		heatmap(
			x,
			y,
			# transpose([real(ε1[a,b,1][i,j]) for a=1:length(x),b=1:length(y)]),
			transpose([real(ε1[i,j,a,b,1]) for a=1:length(x),b=1:length(y)]),
			c = (i==j ? cmap_diag : cmap_offdiag), #cgrad(:cherry),
			aspect_ratio=:equal,
			colorbar = false,
			showaxis = showaxes[j,i],
			top_margin= tm[j,i],
			bottom_margin=bm[j,i],
			left_margin = lm[j,i],
			right_margin = rm[j,i],
			# clim = ( min(ε1...), max(ε1[:,:,i,j]) ),
			clim = ( i==j ? (scale_each_component ? ( 0., maximum(real(ε1[i,j,:,:,1]))) : ( 0., maximum(real(ε1))) ) : (scale_each_component ? ( min( minimum(real(ε1[i,j,:,:,1])), -maximum(real(ε1[i,j,:,:,1]))), max( maximum(real(ε1[i,j,:,:,1])), -minimum(real(ε1[i,j,:,:,1])) ) ) : ( minimum(real(ε1)), maximum(real(ε1))))),
		) for i=1:3,j=1:3 ]

	zoom_plot = heatmap(
			x,
			y,
			# transpose([real(ε1[a,b,1][1,1]) for a=1:length(x),b=1:length(y)]),
			transpose([real(ε1[1,1,a,b,1]) for a=1:length(x),b=1:length(y)]),
			# transpose(ε1[:,:][1,1]),
			c=cmap_diag, #cgrad(:cherry),
			aspect_ratio=:equal,
			legend=false,
			colorbar = true,
			clim = (scale_each_component ? ( 0., maximum(real(ε1[1,1,:,:,1]))) : ( 0., maximum(real(ε1)))),
		)

	if outlines
		plot!(zoom_plot,shape(b),color=nothing,linecolor=:yellow)
		plot!(zoom_plot,shape(s),color=nothing,linecolor=:black)
		plot!(zoom_plot,shape(t),color=nothing,linecolor=:black)
	end
	l = @layout [  a
					grid(3,3)
				]
	plot(
		zoom_plot,
		plot_grid...,
		layout=l,
		size = (600,850),
	)
end

function plot_ε(ε1,g::MaxwellGrid; outlines=false, cmap=cgrad(:cherry))
	plot_ε(ε1,g.x,g.y,outlines=outlines,cmap=cmap)
end

function plot_ε(shapes::AbstractVector{GeometryPrimitives.Shape{2,4,D}},Δx=6.,Δy=4.,Nx=64,Ny=64; outlines=false, cmap=cgrad(:cherry)) where D
	g = MaxwellGrid(Δx,Δy,Nx,Ny)
	plot_ε(εₛ(shapes,g),g.x,g.y,outlines=outlines,cmap=cmap)
end

function plot_ε(shapes::AbstractVector{GeometryPrimitives.Shape{2,4,D}},g::MaxwellGrid; outlines=false, cmap=cgrad(:cherry)) where D
	plot_ε(εₛ(shapes,g),g.x,g.y,outlines=outlines,cmap=cmap)
end

################################################################################
#                        Modesolver plotting functions                         #
################################################################################

function compare_fields(f_mpb,f,ds)
	xlim = (minimum(ds.x),maximum(ds.x))
	ylim = (minimum(ds.y),maximum(ds.y))

    hm_f_mpb_real = [ heatmap(ds.x,ds.y,[real(f_mpb[i,j][ix]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
    hm_f_mpb_imag = [ heatmap(ds.x,ds.y,[imag(f_mpb[i,j][ix]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

    hm_f_real = [ heatmap(ds.x,ds.y,[real(f[ix,i,j]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
    hm_f_imag = [ heatmap(ds.x,ds.y,[imag(f[ix,i,j]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

    hm_f_ratio_real = [ heatmap(ds.x,ds.y,[real(f[ix,i,j])/real(f_mpb[i,j][ix]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
    hm_f_ratio_imag = [ heatmap(ds.x,ds.y,[imag(f[ix,i,j])/imag(f_mpb[i,j][ix]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

    l = @layout [   a   b   c
                    d   e   f
                    g   h   i
                    k   l   m
                    n   o   p
                    q   r   s    ]
    plot(hm_f_mpb_real...,
        hm_f_mpb_imag...,
        hm_f_real...,
        hm_f_imag...,
        hm_f_ratio_real...,
        hm_f_ratio_imag...,
        layout=l,
        size = (1300,1300),
    )
end

function plot_field(f::Array{ComplexF64,4},ds::MaxwellGrid;zind=1)
	xlim = (minimum(ds.x),maximum(ds.x))
	ylim = (minimum(ds.y),maximum(ds.y))

    hm_f_real = [ heatmap(ds.x,ds.y,[real(f[ix,i,j,zind]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
    hm_f_imag = [ heatmap(ds.x,ds.y,[imag(f[ix,i,j,zind]) for i=1:ds.Nx,j=1:ds.Ny]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

    l = @layout [   a   b
                    c   d
					e   f    ]

	# l = @layout [   a   b   c
    #                 d   e	f 	]

    plot(
        hm_f_real[1],
        hm_f_imag[1],
		hm_f_real[2],
        hm_f_imag[2],
		hm_f_real[3],
        hm_f_imag[3],
        layout=l,
        size = (1000,800),
    )
end

function plot_d⃗(H::Array{ComplexF64,2},k::Float64,g::MaxwellGrid;zind=1)
	mn,kpg_mag = calc_kpg(k,g.Δx,g.Δy,g.Δz,g.Nx,g.Ny,g.Nz)
	d = fft(kx_t2c(reshape(H,(2,g.Nx,g.Ny,g.Nz)),kpg_mag,mn),(2:4))
	plot_field(d,g;zind)
end


text:
@testset "Mesh, timestep, and Triangulate utils" begin
    tol = 5e2*eps()
    @testset "Timestep utils" begin
        rk4a,rk4b,rk4c = ck45()
        @test rk4c[1] ≈ 0.0 && rk4c[end] ≈ 1.0
        @test rk4a[1] ≈ 0.0
        @test all(rk4b .> 0)
    end

    @testset "Gmsh reading" begin
        VXY,EToV = readGmsh2D("squareCylinder2D.msh")
        @test size(EToV)==(3031,3)
    end

    @testset "Triangulate utils/example meshes" begin
        # test triangulate
        meshIO = triangulate_domain(SquareDomain())
        (VX,VY),EToV = triangulateIO_to_VXYEToV(meshIO)
        rd = RefElemData(Tri(),2)
        md = MeshData((VX,VY),EToV,rd)
        @test size(EToV,1)==md.num_elements==620
        @test length(VX)==length(VY)==338
        @test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2,3,4]

        h = .1
        meshIO = triangulate_domain(RectangularDomainWithHole()) # square_hole_domain(h)
        (VX,VY),EToV = triangulateIO_to_VXYEToV(meshIO)
        rd = RefElemData(Tri(),2)
        md = MeshData((VX,VY),EToV,rd)
        @test size(EToV,1)==md.num_elements==598
        @test length(VX)==length(VY)==327
        @test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2]
        meshIO2 = refine(meshIO,h)
        @test sort(unique(get_node_boundary_tags(meshIO2,rd,MeshData(triangulateIO_to_VXYEToV(meshIO2)...,rd))))==[0,1,2]

        meshIO = triangulate_domain(Scramjet())
        (VX,VY),EToV = triangulateIO_to_VXYEToV(meshIO)
        rd = RefElemData(Tri(),2)
        md = MeshData((VX,VY),EToV,rd)
        @test size(EToV,1)==md.num_elements==1550
        @test length(VX)==length(VY)==871
        @test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2,3]

        boundary_faces = tag_boundary_faces(meshIO, rd, md, Dict(:wall=>1, :inflow=>2, :outflow=>3))
        @test length(boundary_faces[:wall])==173
        @test length(boundary_faces[:inflow])==14
        @test length(boundary_faces[:outflow])==5
    end
end


text:
using COSMO, Test, LinearAlgebra, SparseArrays, Random, QDLDL, Pkg
rng = Random.MersenneTwister(2401)

function make_test_kkt(P, A, sigma, rho)

    R = length(rho)   == 1 ? ((1.) ./ rho[1]) * I : Diagonal((1.) ./ rho)
    S = length(sigma) == 1 ? (sigma[1]) * I : Diagonal(sigma)

    #compute the full KKT matrix
    K    = [P+S A'; A -R]
    return K
end

solver_types = [COSMO.QdldlKKTSolver
                COSMO.CholmodKKTSolver]

# optional dependencies
if in("Pardiso",keys(Pkg.installed()))
    using Pardiso
    Pardiso.MKL_PARDISO_LOADED[] && push!(solver_types, COSMO.MKLPardisoKKTSolver)
    Pardiso.PARDISO_LOADED[]     && push!(solver_types, COSMO.PardisoDirectKKTSolver)
    Pardiso.PARDISO_LOADED[]     && push!(solver_types, COSMO.PardisoIndirectKKTSolver)
end

solver_tols   = [1e-10, 1e-10, 1e-10, 1e-10, 5e-5]

 @testset "$(solver_types[i]) : KKT solver tests" for i = 1:length(solver_types)

    m = 10
    n = 20

    for rho1 in [rand(1), rand(1)[], rand(m)],
        rho2 in [rand(1), rand(m)],
        sigma in [rand(1), rand(1)[], rand(n)]

        P  = sparse(generate_pos_def_matrix(n, rng))
        A  = sprandn(m, n, 0.2)
        b = randn(m + n)

        F = solver_types[i](P, A, sigma, rho1)
        J = make_test_kkt(P, A, sigma, rho1)
        x = copy(b)

        #test a single solve
        COSMO.solve!(F, x, b)
        @test norm(x - J \ b) <= solver_tols[i]

        #test a rho update and solve
        J = make_test_kkt(P, A, sigma, rho2)
        x = copy(b)
        COSMO.update_rho!(F, rho2)
        COSMO.solve!(F, x, b)
        @test norm(x - J \ b) <= solver_tols[i]

     end

 end

nothing


text:
# Compass search as described on page 18 (402) in Kolda2003:
#  Kolda, Tamara G., Robert Michael Lewis, and Virginia Torczon. "Optimization
#  by direct search: New perspectives on some classical and modern methods.",
#  SIAM review 45.3 (2003): 385-482.
function compass_search(f, n; max_fevals = 1e6, delta_tol = 1e-20,
  step_size = 1.0, x = false, known_fmin = false)

  @assert delta_tol > 0
  @assert step_size > delta_tol

  if x == false
    x = randn(n, 1)
  end

  directions = [Matrix{Float64}(I, n,n) -Matrix{Float64}(I, n, n)]

  num_fevals = 0
  fc = fbest = Inf
  candidate = zeros(n, 1)

  while(step_size > delta_tol && num_fevals <= max_fevals)
    # Check all directions to find a better point
    found_better = false
    for direction in 1:(2*n)
      candidate = x + step_size .* directions[:, direction]

      fc = f(candidate)
      num_fevals += 1

      if fc < fbest
        found_better = true
        break
      end

    end

    if found_better
      fbest = fc
      x = candidate

      if fbest == known_fmin
        break
      end
    else
      step_size = step_size / 2
    end
  end

  return x, fbest, num_fevals
end


text:
export step!, run!

"""
update e.elites with the best individuals from the current population or
existing elites, based on fitness
"""
function elites_generation(e::AbstractEvolution)
    pop = sort([e.population; e.elites])
    n_elites = length(e.elites)
    for i in 1:n_elites
        e.elites[i] = deepcopy(pop[length(e.population)+i])
    end
end

"""
generation should be used for any additional record keeping necessary in an
algorithm. It is called after populate and evaluate, so all individuals should
have genes and fitness reflecting their current state. The default function
calls elites_generation if e.elites is defined, otherwise does nothing.
"""
function generation(e::AbstractEvolution)
    if hasfield(typeof(e), :elites)
        elites_generation(e)
    end
end

"""
The generic iteration of an evolution. Calls populate, evaluate, and generation.
Also calls log_gen and save_gen based on the provided config values. Subclasses
of AbstractEvolution should override the populate, evaluate, or generation
functions rather than overriding this function.
"""
function step!(e::AbstractEvolution)
    e.gen += 1
    if e.gen > 1
        populate(e)
    end
    evaluate(e)
    generation(e)
    if ((e.config.log_gen > 0) && mod(e.gen, e.config.log_gen) == 0)
        log_gen(e)
    end
    if ((e.config.save_gen > 0) && mod(e.gen, e.config.save_gen) == 0)
        save_gen(e)
    end
end

"Call step!(e) e.config.n_gen times consecutively"
function run!(e::AbstractEvolution)
    for i in (e.gen+1):e.config.n_gen
        step!(e)
    end
end


text:
module PkgVerifierServerPrototype

import HttpCommon
import Mux
import PkgVerifierPrototype
import JSON

export serve

repo = ""
certificates = Dict{ASCIIString,AbstractString}()

function get_certificate(req)
    package = req[:params][:package]
    cert = ""
    try
        cert = certificates[package]
    catch error
        return "\"FAILURE\""
    end

    return cert
end

function put_certificate(req)
    package = req[:params][:package]

    cert = ASCIIString(req[:data])

    try
        PkgVerifierPrototype.open_data_certificate(repo, cert)
    catch error
        return string(error)
    end

    certificates[package] = cert

    return "\"SUCCESS\""
end

@Mux.app verifier_server_app = (
    Mux.defaults,
    Mux.page("/cert/:package", get_certificate),
    Mux.page("/push/:package", put_certificate),
    Mux.notfound())

function serve(local_repo::AbstractString)
    global repo = local_repo
    Mux.serve(verifier_server_app)
end

end # module

max_stars_repo_name: LachlanGunn/PkgVerifierServerPrototype

max_issues_repo_name: LachlanGunn/PkgVerifierServerPrototype

max_forks_repo_name: LachlanGunn/PkgVerifierServerPrototype


text:
using Roots
using Test
using JSON
using Printf

## testing/benchmarking derivative-free methods

mutable struct Func1
    name :: Symbol
    val :: Function
    x0 :: Function
    alpha # a value or values
    params :: Vector{Any}
end
function Base.show(io::IO, f::Func1)
    @printf io "Func(%s)" f.name
end

#Construct a function object,
macro Func1(name)
    @gensym f p x
    esc(quote
        $f = Func1($name, val, x0, alpha, params)
        push!(known_functions, $f)
        $f
    end)
end

known_functions = Func1[]

## This set of tests is useful for benchmarking the number of function
## calls, failures, and max errors for the various derivative_free methods.

# easy ones
func1 = let
    val = (_, x) -> cos(x) - x/2
    a,b, N = 0.0, pi/2, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 1.0298665293222589
    @Func1 :func1
end

func2 = let
    val = (_, x) -> exp(x) - x^4
    a, b, N = 5.0, 20.0, 11
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = find_zeros(x->val(0,x), -5.0, 20.0)
    @Func1 :func2
end

func3 = let
    val = (_, x) -> x^5 - x - 1
    a, b, N = 0.5, 2.0, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = .1673039782614187
    @Func1 :func3
end

# wider range
func4 = let
    val = (_,x) -> (1 + cos(x)) * (exp(x) - 2)
    a, b, N = 0.0, 1.7, 15
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.6931471805599453
    @Func1 :func4
end

func5 = let
    val = (_,x) -> 2*x - exp(-x)
    a, b, N = -5.0, 5.0, 15
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.35173371124919584
    @Func1 :func5
end

# once over hump, hard to get back to 0
func6 = let
    val = (_,x) -> x * exp(-x)
    a, b, N = -5.0, 5.0, 15
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.0
    @Func1 :func6
end

# once outside of extrema, points away from zero
func7 = let
    val = (_,x) ->  20.0 * x / ( 100.0 * x^2 + 1.0)
    a, b, N = -0.2, 0.2, 15
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.0
    @Func1 :func7
end

# should have some failures here

func8 = let
 val = (_,x) -> cos(x) - x/2
    a, b, N = 0.0, 3pi/4, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 1.0298665293222589
    @Func1 :func8
end

func9 = let
    val = (_,x) -> tanh(x) - tan(x)
    a, b, N = 5.0, 10.0, 8
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = find_zeros(x -> val(0,x), 5, 12)
    @Func1 :func9
end

func10 = let
    val = (_,x) -> asin(x^2 - 1) - x/2 + 1
    a, b, N = -0.5, 0.9, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = find_zeros(x -> val(0,x), -0.5, 1.0)
    @Func1 :func10
end

func11 = let
    val = (_,x) -> exp(-x) - cos(x)
    a, b, N = -2.0, 2.0, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = find_zeros(x -> exp(-x) - cos(x), -2,2)
    @Func1 :func11
end

func12 = let
    val = (_,x) -> sqrt(x) - 1/x - 3
    a, b, N = 5, 20, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 9.633595562832696
    @Func1 :func12
end

## multiplicity: http://ir.igsnrr.ac.cn/bitstream/311030/8840/1/%E4%BE%AF%E9%BA%9F%E7%A7%91(SCI)2.pdf
func13 = let
 val = (_,x) -> (x - sqrt(5))^4 / ((x-1)^2 + 2)
    a, b, N = 3.0, 4.5, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha =  2.236067977499790
    @Func1 :func13
end

func14 = let
 val = (_,x) -> (sin(x)^2 - 2x + 1)^5
    a, b, N = 3.0, 5.0, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.71483582544138924
    @Func1 :func14
end

func15 = let
 val = (_,x) -> (8x*exp(-x^2) -2x - 3)^8
    a, b, N = -2.0,-1.5, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = -1.7903531791589544
    @Func1 :func15
end

# this is 1e6 * (x-1)^7, has many "zeros" near 1
func16 = let
    val = (_,x) -> 1e6 * (x^7 -7x^6 +21x^5 -35x^4 +35x^3-21x^2+7x-1)
    a, b, N = -0.0, 2.0, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 1.0
    @Func1 :func16
end

func17 = let
 val = (_,x) -> (exp(-x^2 + x + 3) - x + 2)^9
    a, b, N = 2.4, 2.5, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha =  2.4905398276083051
    @Func1 :func17
end

func18 = let
 val = (_,x) -> (exp(-x) + 2sin(x))^4
    a, b, N = 3.0, 3.4, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 3.1627488709263654
    @Func1 :func18
end

## hard w/o aggressive bracketing
func19 = let
 val = (_,x) -> cbrt(x)
    a, b, N = -3.0, 3.0, 10
    x0(n) = range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.0
    @Func1 :func19
end

func20 = let
    val = (_,x) -> x < -1/4 ? 3/4*x-0.3125 : x < 1/4 ? 2x : 3/4*x+0.3125
    a,b,N = -1.0, 1.0, 10
    x0(n) =  range(a, stop=b, length=N)[n]
    params = 1:N
    alpha = 0.0
    @Func1 :func20
end

mutable struct MethodResults1
    name
    problems::Int
    evalcount :: Int
maxresidual :: Float64
maxdifference :: Float64
    failures :: Vector{Tuple{Func1, Int}}
end
MethodResults1() = MethodResults1(nothing, 0, 0, 0.0, Inf,  Tuple{Func1, Int}[])
Base.show(io::IO, results::MethodResults1) =
    print(io, "MethodResults($(results.name), evalcount=$(results.evalcount), numfailures=$(length(results.failures)), maxresidual=$(results.maxresidual))")

## Run a method on all known functions.
mindiff(a, alpha) =  minimum([a - i for i in alpha])
function run_tests(method; verbose=false, trace=false, name=nothing, abandon=false)
    results = MethodResults1()
    results.name = name
    results.problems = 0
    for f in known_functions
        for i in 1:length(f.params)
            p = f.params[i]
            evalcount = 0
            function feval(x)
                evalcount += 1
                result = f.val(p, x)
                trace && @printf "%s[%d]: %s ⇒ %s\n" f i x result
                result
            end
            result, residual = nothing, nothing
            try
                result = method(feval, f.x0(p))
                md = abs(mindiff(result, f.alpha) )
                isnan(result) && error("NaN")
                (md > 1) && error("ran away")

                results.problems += 1
                results.evalcount += evalcount
                residual = f.val(p, result)

                verbose && @printf "%s[%d] ⇒ %d / %s, residual %.5e\n" f.name i evalcount result residual

                ## Some functions might return non-real values on failures
                if isa(result, AbstractFloat) && isa(residual, AbstractFloat) && isfinite(residual)
                    results.maxresidual = max(results.maxresidual, abs(residual))
                    results.maxdifference = max(results.maxdifference, mindiff(result, f.alpha))
                end

                verbose && abs(residual) >1/10 && @printf "Large Residual [%s] %d/%d\n" f.name i p

            catch ex
                verbose && @printf "%s[%d] ⇒ FAILED: %s\n" f.name i ex
                push!(results.failures, (f, i))
                abandon && rethrow(ex)
            end
        end
    end
    results
end

avg(x) = sum(x)/length(x)
D(f, h=1e-4) = x -> (f(x+h) - f(x-h))/(2h)
D2(f, h=1e-4) = x -> (f(x+h) - 2f(x) + f(x-h))/h^2

if !isinteractive()

    @testset "derivative free methods" begin

        ## Test for failures, ideally all of these would be 0
        ## test for residual, ideally small
        ## test for evaluation counts, ideally not so low for these problems

        ## basic methods
        Ms = [Roots.Order1(), Roots.Order1B(), Roots.Order2(), Roots.Order2B(),
              Roots.Order5(), Roots.Order8(), Roots.Order16()]
        results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]

        failures = [length(result.failures) for result in results]
        residuals = [result.maxresidual for result in results]
        cnts = [result.evalcount/result.problems for result in results]

        @test maximum(failures) <= 60
        @test maximum(residuals) <= 5e-5
        @test avg(cnts) <= 40

        ## methods which fall back to bisection when bracket found
        Ms = [Roots.Order0()]
        results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]

        failures = [length(result.failures) for result in results]
        residuals = [result.maxresidual for result in results]
        cnts = [result.evalcount/result.problems  for result in results]

        @test maximum(failures) <= 40
        @test maximum(residuals) <= 1e-5
        @test avg(cnts) <= 30

        ## Newton and Halley
        Fs = [(f,b) -> find_zero((f, D(f)), b, Roots.Newton()),
              (f,b) -> find_zero((f, D(f), D2(f)), b, Roots.Halley()),
              (f, b) -> find_zero((f, D(f), D2(f)), b, Roots.Schroder())
              ]
        results = [run_tests(F) for F in Fs]

        failures = [length(result.failures) for result in results]
        residuals = [result.maxresidual for result in results]
        cnts = [result.evalcount/result.problems  for result in results]

        @test maximum(failures) <= 70
        @test maximum(residuals) <= 5e-4
        @test avg(cnts) <= 125

    end

    @testset "derivative free, non Float64" begin

        Ms = [Roots.Order0(), Roots.Order1(), Roots.Order1B(), Roots.Order2(), Roots.Order2B(),
              Roots.Order5(), Roots.Order8(), Roots.Order16()]
        Ts = [Float16, Float32, BigFloat]

        fn = x -> x^3 - 2x - 5
        alpha = 2.094551481542326591482386540579302963857306105628239180304128529045312189983499
        x0 = 2.0

        for T in Ts
            for M in Ms
                xstar = find_zero(fn, T(x0), M)
                @test xstar ≈ T(alpha) atol=max(sqrt(eps(T)), eps())
                @test isa(xstar, T)
            end
        end

        for T in Ts
            xstar = Roots.find_zero((fn,D(fn,sqrt(eps(T)))), T(x0), Roots.Newton())
            @test xstar ≈ T(alpha)  atol=max(sqrt(eps(T)), eps())
            @test isa(xstar, T)

            xstar = Roots.find_zero((fn,D(fn,sqrt(eps(T))),D2(fn,sqrt(eps(T)))), T(x0), Roots.Halley())
            @test xstar ≈ T(alpha) atol=max(sqrt(eps(T)), eps())
            @test isa(xstar, T)

            xstar = Roots.find_zero((fn,D(fn,sqrt(eps(T))),D2(fn,sqrt(eps(T)))), T(x0), Roots.Schroder())
            @test xstar ≈ T(alpha) atol=max(sqrt(eps(T)), eps())
            @test isa(xstar, T)
        end
    end
end


text:
# Autogenerated wrapper script for GAP_pkg_normalizinterface_jll for x86_64-linux-musl-cxx03
export NormalizInterface

using GAP_jll
using GAP_lib_jll
using normaliz_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_normalizinterface")
JLLWrappers.@declare_file_product(NormalizInterface)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll, normaliz_jll)
    JLLWrappers.@init_file_product(
        NormalizInterface,
        "lib/gap/NormalizInterface.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DashHtmlComponents

function Section(title::String, links::String, description::Union{Nothing, String}, headerStyle = ())
    return(
        html_div(
            className = "toc--section",
            children = (
                html_h2(
                    title,
                    style = merge((borderBottom = "thin lightgrey solid", marginTop = 50,), headerStyle, )
                ),
                html_div(description),
                html_ul(
                    links,
                    className = "toc--chapters"
                )
            )
        )
    )
end


text:
export stack, unstack

"""
`stack(t, by = pkeynames(t); select = excludecols(t, by), variable = :variable, value = :value)`

Reshape a table from the wide to the long format. Columns in `by` are kept as indexing columns.
Columns in `select` are stacked. In addition to the id columns, two additional columns labeled `variable` and `value`
are added, containg the column identifier and the stacked columns.

## Examples

```jldoctest stack
julia> t = table(1:4, [1, 4, 9, 16], [1, 8, 27, 64], names = [:x, :xsquare, :xcube], pkey = :x);

julia> stack(t)
Table with 8 rows, 3 columns:
x  variable  value
──────────────────
1  :xsquare  1
1  :xcube    1
2  :xsquare  4
2  :xcube    8
3  :xsquare  9
3  :xcube    27
4  :xsquare  16
4  :xcube    64
```
"""
function stack(t::D, by = pkeynames(t); select = isa(t, NDSparse) ? valuenames(t) : excludecols(t, by),
    variable = :variable, value = :value) where {D<:Dataset}

    (by != pkeynames(t)) && return stack(reindex(t, by, select); variable = :variable, value = :value)

    valuecols = columns(t, select)
    valuecol = [valuecol[i] for i in 1:length(t) for valuecol in valuecols]

    labels = fieldnames(valuecols)
    labelcol = [label for i in 1:length(t) for label in labels]

    bycols = map(arg -> repeat(arg, inner = length(valuecols)), columns(t, by))
    convert(collectiontype(D), Columns(bycols), Columns(labelcol, valuecol, names = [variable, value]))
end

function unstack(::Type{D}, ::Type{T}, key, val, cols::AbstractVector{S}) where {D <:Dataset, T, S}
    dest_val = Columns((DataValues.DataValueArray{T}(length(val)) for i in cols)...; names = cols)
    for (i, el) in enumerate(val)
        for j in el
            k, v = j
            isnull(columns(dest_val, S(k))[i]) || error("Repeated values with same label are not allowed")
            columns(dest_val, S(k))[i] = v
        end
    end
    convert(collectiontype(D), key, dest_val)
end

"""
`unstack(t, by = pkeynames(t); variable = :variable, value = :value)`

Reshape a table from the long to the wide format. Columns in `by` are kept as indexing columns.
Keyword arguments `variable` and `value` denote which column contains the column identifier and
which the corresponding values.

## Examples

```jldoctest unstack
julia> t = table(1:4, [1, 4, 9, 16], [1, 8, 27, 64], names = [:x, :xsquare, :xcube], pkey = :x);

julia> long = stack(t)
Table with 8 rows, 3 columns:
x  variable  value
──────────────────
1  :xsquare  1
1  :xcube    1
2  :xsquare  4
2  :xcube    8
3  :xsquare  9
3  :xcube    27
4  :xsquare  16
4  :xcube    64

julia> unstack(long)
Table with 4 rows, 3 columns:
x  xsquare  xcube
─────────────────
1  1        1
2  4        8
3  9        27
4  16       64
```
"""
function unstack(t::D, by = pkeynames(t); variable = :variable, value = :value) where {D<:Dataset}
    tgrp = groupby((value => identity,), t, by, select = (variable, value))
    S = eltype(colnames(t))
    cols = S.(union(columns(t, variable)))
    T = eltype(columns(t, value))
    unstack(D, T isa Type{<:DataValue} ? eltype(T) : T, pkeys(tgrp), columns(tgrp, value), cols)
end


text:
abstract AbstractReferenceElement

abstract AbstractReferenceElement_1D <: AbstractReferenceElement
abstract AbstractReferenceElement_2D <: AbstractReferenceElement
abstract AbstractReferenceElement_3D <: AbstractReferenceElement

vertices(refele::AbstractReferenceElement) = refele.vertices
edges(refele::AbstractReferenceElement) = refele.edges
faces(refele::Union(AbstractReferenceElement_2D,
                    AbstractReferenceElement_3D)) = refele.faces

cell(refele::AbstractReferenceElement_3D) = refele.cell
facet(refele::AbstractReferenceElement_2D) = refele.edges
facet(refele::AbstractReferenceElement_3D) = refele.faces

dim(::AbstractReferenceElement_1D) = 1
dim(::AbstractReferenceElement_2D) = 2
dim(::AbstractReferenceElement_3D) = 3

# TODO this can fail for entity(refline, 3) for example
function entity(refele::AbstractReferenceElement, dim::Int)
    if dim == 1
        return refele.edges
    elseif dim == 2
        return refele.faces
    elseif dim == 3
        return refele.cell
    end
end

immutable ReferenceLine <: AbstractReferenceElement_1D
    vertices::Vector{Vector{Float64}}
    edges::Vector{Vector{Int}}
end

function ReferenceLine()
    vertices = Array(Vector{Float64}, 2)
    vertices[1] = [-1]
    vertices[2] = [1]
    edges = Array(Vector{Int}, 1)
    edges[1] = [1, 2]

    ReferenceLine(vertices, edges)
end
volume(::ReferenceLine) = 2.0

immutable ReferenceTriangle <: AbstractReferenceElement_2D
    vertices::Vector{Vector{Float64}}
    edges::Vector{Vector{Int}}
    faces::Vector{Vector{Int}}
end

function ReferenceTriangle()
    vertices = Array(Vector{Float64}, 3)
    vertices[1] = [0, 0]
    vertices[2] = [1, 0]
    vertices[3] = [0, 1]

    edges = Array(Vector{Int}, 3)
    edges[1] = [1, 2]
    edges[2] = [2, 3]
    edges[3] = [3, 1]

    faces = Array(Vector{Int}, 1)
    faces[1] = [1, 2, 3]

    ReferenceTriangle(vertices, edges, faces)
end
volume(::ReferenceTriangle) = 0.5

immutable ReferenceSquare <: AbstractReferenceElement_2D
    vertices::Vector{Vector{Float64}}
    edges::Vector{Vector{Int}}
    faces::Vector{Vector{Int}}
end

function ReferenceSquare()
    vertices = Array(Vector{Float64}, 4)
    vertices[1] = [-1, -1]
    vertices[2] = [ 1, -1]
    vertices[3] = [ 1,  1]
    vertices[4] = [-1,  1]

    edges = Array(Vector{Int}, 4)
    edges[1] = [1, 2]
    edges[2] = [2, 3]
    edges[3] = [3, 4]
    edges[4] = [4, 1]

    faces = Array(Vector{Int}, 1)
    faces[1] = [1, 2, 3, 4]

    ReferenceSquare(vertices, edges, faces)
end
volume(::ReferenceSquare) = 4.0

immutable ReferenceTetrahedron <: AbstractReferenceElement_3D
    vertices::Vector{Vector{Float64}}
    edges::Vector{Vector{Int}}
    faces::Vector{Vector{Int}}
    cell::Vector{Int}
end

function ReferenceTetrahedron()
    vertices = Array(Vector{Float64}, 4)
    vertices[1] = [0, 0 ,0]
    vertices[2] = [1, 0, 0]
    vertices[3] = [0, 1, 0]
    vertices[4] = [0, 0, 1]

    edges = Array(Vector{Int}, 6)
    edges[1] = [1, 2]
    edges[2] = [1, 3]
    edges[3] = [1, 4]
    edges[4] = [2, 3]
    edges[5] = [2, 4]
    edges[6] = [3, 4]

    faces = Array(Vector{Int}, 4)
    faces[1] = [1, 3, 2]
    faces[2] = [1, 2, 4]
    faces[3] = [1, 4, 3]
    faces[4] = [2, 3, 4]

    cell = [1, 2, 3, 4]

    ReferenceTetrahedron(vertices, edges, faces, cell)
end
volume(::ReferenceTetrahedron) = 1/6

immutable ReferenceCube <: AbstractReferenceElement_3D
    vertices::Vector{Vector{Float64}}
    edges::Vector{Vector{Int}}
    faces::Vector{Vector{Int}}
    cell::Vector{Int}
end

function ReferenceCube()
    vertices = Array(Vector{Float64}, 8)
    vertices[1] = [-1, -1, -1]
    vertices[2] = [ 1, -1, -1]
    vertices[3] = [ 1,  1, -1]
    vertices[4] = [-1,  1, -1]
    vertices[5] = [-1, -1,  1]
    vertices[6] = [ 1, -1,  1]
    vertices[7] = [ 1,  1,  1]
    vertices[8] = [-1,  1,  1]

    edges = Array(Vector{Int}, 12)
    edges[1] = [1, 2]
    edges[2] = [1, 3]
    edges[3] = [1, 5]
    edges[4] = [2, 4]
    edges[5] = [2, 6]
    edges[6] = [3, 4]
    edges[7] = [3, 7]
    edges[8] = [4, 8]
    edges[9] = [5, 6]
    edges[10] = [5, 7]
    edges[11] = [6, 8]
    edges[12] = [7, 8]

    faces = Array(Vector{Int}, 6)
    faces[1] = [1, 2, 6, 5]
    faces[2] = [1, 3, 4, 2]
    faces[3] = [1, 5, 7, 3]
    faces[4] = [2, 4, 8, 6]
    faces[5] = [3, 7, 8, 4]
    faces[6] = [5, 6, 8, 7]

    cell = [1, 2, 3, 4, 5, 6, 7, 8]

    ReferenceCube(vertices, edges, faces, cell)
end
volume(::ReferenceCube) = 8

max_stars_repo_name: JuliaTagBot/JFET

max_issues_repo_name: JuliaTagBot/JFET

max_forks_repo_name: JuliaTagBot/JFET


text:
"""
    newton(x0::T,xL::T,xU::T,f::Function,df::Function,envp1::T,envp2::T)

Defines a local 1D newton method to solve for the root of `f` between the bounds
`xL` and `xU` using `x0` as a starting point. The derivative of `f` is `df`. The
inputs `envp1` and `envp2` are the envelope calculation parameters.
"""
@inline function newton(x0::T,xL::T,xU::T,f::Function,df::Function,envp1::T,envp2::T) where {T<:Real}
  dfk::T = zero(x0)

  xk::T = max(xL,min(x0,xU))
  fk::T = f(xk,envp1,envp2)

  for i=1:MC_param.env_max_int
    dfk = df(xk,envp1,envp2)
    if (abs(fk)<MC_param.env_tol)
      return xk
    end
    if (dfk == zero(x0))
      error("NEWTON EXCEPTION")
    elseif (xk==xL && fk/dfk>zero(x0))
      return xk
    elseif (xk==xU && fk/dfk<zero(x0))
      return xk
    end
    xk = max(xL,min(xU,xk-fk/dfk))
    fk = f(xk,envp1,envp2)
  end
  error("NEWTON EXCEPTION")
end

"""
    secant(x0::T,x1::T,xL::T,xU::T,f::Function,envp1::T,envp2::T)  where {T<:Real}

Defines a local 1D secant method to solve for the root of `f` between
the bounds `xL` and `xU` using `x0` and `x1` as a starting points. The inputs
`envp1` and `envp2` are the envelope calculation parameters.
"""
@inline function secant(x0::T,x1::T,xL::T,xU::T,f::Function,envp1::T,envp2::T) where {T<:Real}
  xkm::T = max(xL,min(xU,x0))
  xk::T = max(xL,min(xU,x1))
  fkm::T = f(xkm,envp1,envp2)

  for i=1:MC_param.env_max_int
    fk = f(xk,envp1,envp2)
    Bk::T = (fk-fkm)/(xk-xkm)
    if (abs(fk)<MC_param.env_tol)
      return xk
    end
    if (Bk == zero(x0))
      error("SECANT EXCEPTION")
    elseif ((xk==xL)&(fk/Bk>zero(x0)))
      return xk
    elseif ((xk==xU)&(fk/Bkzero(x0)))
      return xk
    end
    xkm = xk
    fkm = fk
    xk = max(xL,min(xU,xk-fk/Bk))
  end
  error("SECANT EXCEPTION")
end

"""
    golden_section(xL::T,xU::T,f::Function,envp1::T,envp2::T) where {T<:Real}

Defines a local 1D golden section method to solve for the root of `f` between
the bounds `xL` and `xU` using `x0` as a starting point. Define iteration used
in golden section method. The inputs `envp1` and `envp2` are the envelope
calculation parameters.
"""
@inline function golden_section(xL::T,xU::T,f::Function,envp1::T,envp2::T) where {T<:Real}
  fL::T = f(xL,envp1,envp2)
  fU::T = f(xU,envp1,envp2)

  if (fL*fU > zero(fL))
    error("GOLDEN EXCEPTION")
  end
  xm::T = xU-(2.0-golden)*(xU-xL)
  fm::T = f(xm,envp1,envp2)
  return golden_section_it(1,xL,fL,xm,fm,xU,fU,f,envp1,envp2)
end
"""
    golden_section_it(init::Int64,a::T,fa::T,b::T,fb::T,c::T,
                      fc::T,f::Function,envp1::T,envp2::T) where {T<:Real}

Define iteration used in golden section method. The inputs `fa`,`fb`, and `fc`,
are the function `f` evaluated at `a`,`b`, and `c` respectively. The inputs
`envp1` and `envp2` are the envelope calculation parameters. The value `init` is
the iteration number of the golden section method.
"""
@inline function golden_section_it(init::Int64,a::T,fa::T,b::T,fb::T,c::T,
                                   fc::T,f::Function,envp1::T,envp2::T) where {T<:Real}
  b_t_x::Bool = (c-b > b-a)
  if (b_t_x)
    x::T = b + (2.0-golden)*(c-b)
  else
    x = b - (2.0-golden)*(b-a)
  end
  itr::Int64 = init
  if (abs(c-a)<MC_param.env_tol*(abs(b)+abs(x))||(itr>MC_param.env_max_int))
    return (c+a)/2.0
  end
  itr += 1
  fx::T = f(x,envp1,envp2)
  if (b_t_x)
    if (fa*fx<zero(fa))
      golden_section_it(itr,a,fa,b,fb,x,fx,f,envp1,envp2)
    else
      golden_section_it(itr,b,fb,x,fx,c,fc,f,envp1,envp2)
    end
  else
    if (fa*fb<(fa))
      golden_section_it(itr,a,fa,x,fx,b,fb,f,envp1,envp2)
    else
      golden_section_it(itr,x,fx,b,fb,c,fc,f,envp1,envp2)
    end
  end
end


text:
struct LatPhysLattice{LT <: LatPhysBase.AbstractLattice} <: AbstractLattice
    lattice::LT
    neighs::Matrix{Int}
end

function LatPhysLattice(lattice::LatPhysBase.AbstractLattice)
    # Build lookup table for neighbors
    # neighs[:, site] = list of neighoring site indices
    bs = sort(LatPhysBase.bonds(lattice), by = LatPhysBase.from)
    bonds_per_site = findfirst(b -> LatPhysBase.from(b) > 1, bs) - 1
    neighs = Matrix{Int}(undef, bonds_per_site, LatPhysBase.numSites(lattice))
    for site_idx in 1:LatPhysBase.numSites(lattice)
        for bond_idx in 1:bonds_per_site
            bond = bs[bonds_per_site * (site_idx - 1) + bond_idx]
            @assert from(bond) == site_idx
            neighs[bond_idx, site_idx] = to(bond)
        end
    end

    LatPhysLattice(lattice, neighs)
end

@inline Base.length(l::LatPhysLattice) = LatPhysBase.numSites(l.lattice)
@inline function neighbors(l::LatPhysLattice, directed::Val{true})
    ((LatPhysBase.from(b), LatPhysBase.to(b)) for b in LatPhysBase.bonds(l))
end
@inline function neighbors(l::LatPhysLattice, directed::Val{false})
    (
        (LatPhysBase.from(b), LatPhysBase.to(b))
        for b in LatPhysBase.bonds(l)
        if LatPhysBase.from(b) < LatPhysBase.to(b)
    )
end
@inline function neighbors(l::LatPhysLattice, site_index::Integer)
    l.neighs[:, site_index]
end

@inline neighbors_lookup_table(l::LatPhysLattice) = copy(l.neighs)


text:
#! /usr/bin/julia

# Rosetta Code, Standard deviation

function makerunningstd()
    a = zero(Float64)
    b = zero(Float64)
    n = zero(Int64)
    function runningstd(x)
        a += x
        b += x^2
        n += 1
        std = sqrt(n*b - a^2)/n
        return std
    end
    return runningstd
end

test = [2, 4, 4, 4, 5, 5, 7, 9]

rstd = makerunningstd()

println("Perform a running standard deviation of ", test)
for i in test
    println(i, " => ", rstd(i))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# abstract types
abstract type VLAbstractModelObject end

# concrete types -
struct VLResult{T}
    value::T
end

struct VLJuliaModelObject <: VLAbstractModelObject

    # data -
    path_to_model_file::String
    path_to_output_dir::String
    defaults_file_name::String
    model_type::Symbol

    # constructor -
    function VLJuliaModelObject(path_to_model_file, path_to_output_dir;
        defaults_file_path::String="Defaults.toml", model_type::Symbol=:static)

        # build new model object -
        this = new(path_to_model_file, path_to_output_dir, defaults_file_path, model_type)
    end
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
using ParallelAccelerator

@acc function bicg(m::Int, n::Int, A::Matrix{Float64},
     s::Vector{Float64}, q::Vector{Float64}, p::Vector{Float64}, r::Vector{Float64})

     for i = 1:m
          s[i] = 0.0
     end
     for i = 1:n
          q[i] = 0.0
          for j = 1:m
               s[j] += r[i] * A[i, j]
               q[i] += A[i, j] * p[j]
          end
     end
end

@inbounds function bicg_o(m::Int, n::Int, A::Matrix{Float64},
     s::Vector{Float64}, q::Vector{Float64}, p::Vector{Float64}, r::Vector{Float64})

     tic()
     for i = 1:m
          @fuse 1 s[i] = 0.0
     end
     for i = 1:n
          q[i] = 0.0
          for j = 1:m
               s[j] += r[i] * A[i, j]
               q[i] += A[i, j] * p[j]
          end
     end
     println("Julia $(toq())")
end

m = 20000
n = m

A = rand(n, m)
st = Vector{Float64}(m)
qt = Vector{Float64}(n)
p = rand(m)
r = rand(n)

bicg(m, n, A, st, qt, p, r)

s = Vector{Float64}(m)
q = Vector{Float64}(n)

bicg_o(m, n, A, s, q, p, r)

function compare(A, B)
    for i in 1:n
        if A[i] != B[i]
            println("different $((i)) : $(A[i]) != $(B[i])")
            return false
        end
    end
    return true
end

println(compare(s, st))
println(compare(q, qt))

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
using .Flux: Flux, fmap, functor

# like gpu() from flux
tojax(m) = fmap(x -> adapt(JaxArray, x), m)

# Jax will recursively call this on any jl object (like jl_wrap<identity>)
# until he gets only fundamental types out. THankfully, that's mostly
# what functor does.
_flatten(x) = begin
    d, f = functor(x)
    tup = tuple(d...)
    #@info "Called Chain_flatten: $x" tup

    return tup, f
end

_unflatten(f, x) = begin
    #@info("called chain_unf $f - $x")
    return f(x)
end


text:
# ==============
# API extensions
# ==============

function MathematicalSystems.statedim(H::HybridSystem)
    ns = [statedim(H, i) for i in 1:nmodes(H)]
    if !all(==(ns[1]), ns)
        @warn("the number of state dimensions differs across locations")
    end
    return maximum(ns)
end

function MathematicalSystems.inputdim(H::HybridSystem)
    ms = [inputdim(H, i) for i in 1:nmodes(H)]
    if !all(==(ms[1]), ms)
        @warn("the number of input dimensions differs across locations")
    end
    return maximum(ms)
end

function MathematicalSystems.isconstrained(X::LazySet)
    return !isuniversal(X)
end

function MathematicalSystems.isconstrained(v::AbstractVector{<:LazySet})
    return !all(isuniversal, v)
end

# ===========
# Indentation
# ===========

mutable struct Indentation
    ind::String
    step::String

    # default indentation: two spaces
    function Indentation(ind="", step="  ")
        return new(ind, step)
    end
end

function indent!(ind::Indentation)
    ind.ind = ind.ind * ind.step
end

function dedent!(ind::Indentation)
    ind.ind = ind.ind[1:end-(length(ind.step))]
end

function _write_indented(io, string, indentation)
    write(io, "$(indentation.ind)$string")
end

function _dedent(indentation)
    return
end

# ==============
# Variable names
# ==============

function _variable_name(id, dictionary)
    get(dictionary, id, "x" * string(id))
end

function _input_name(id, dictionary)
    get(dictionary, id, "u" * string(id))
end

# ======================
# Creation of model file
# ======================

function writesxmodel(filename, system; dictionary=Dict())
    open(filename, "w") do io
        write(io,
"<?xml version=\"1.0\" encoding=\"iso-8859-1\"?>\n" *
"<sspaceex xmlns=\"http://www-verimag.imag.fr/xml-namespaces/sspaceex\" version=\"0.2\" math=\"SpaceEx\">\n"
        )
        indentation = Indentation()
        indent!(indentation)
        _write_system(io, system, dictionary, indentation)
        dedent!(indentation)
        write(io, "</sspaceex>\n")
    end
    nothing
end

function _write_system(io, system, dictionary, indentation)
    _write_component(io, system, dictionary, indentation)
end

function _write_component(io, system, dictionary, indentation)
    _write_indented(io, "<component id=\"system\">\n", indentation)
    indent!(indentation)
    _write_parameters(io, system, dictionary, indentation)
    _write_locations(io, system, dictionary, indentation)
    _write_transitions(io, system, dictionary, indentation)
    dedent!(indentation)
    _write_indented(io, "</component>\n", indentation)
end

function _write_parameters(io, system, dictionary, indentation)
    _write_state_variables(io, system, dictionary, indentation)
    _write_input_variables(io, system, dictionary, indentation)
end

function _write_input_variables(io, system, dictionary, indentation)
    m = inputdim(system)
    controlled = false
    for i in 1:m
        name = _input_name(i, dictionary)
        _write_parameter(io, name, controlled, dictionary, indentation)
    end
end

function _write_state_variables(io, system, dictionary, indentation)
    n = statedim(system)
    controlled = true
    for i in 1:n
        name = _variable_name(i, dictionary)
        _write_parameter(io, name, controlled, dictionary, indentation)
    end
end

function _write_parameter(io, name, controlled, dictionary, indentation)
    _write_indented(io, "<param name=\"$name\" type=\"real\" d1=\"1\" d2=\"1\" " *
        "local=\"false\" dynamics=\"any\" controlled=\"$controlled\" />\n",
        indentation)
end

function _write_locations(io, system::AbstractContinuousSystem, dictionary,
                          indentation)
    _write_location(io, system, 1, dictionary, indentation)
end

function _write_locations(io, system::HybridSystem, dictionary,
                          indentation)
    for (i, loc) in enumerate(system.modes)
        _write_location(io, loc, i, dictionary, indentation)
    end
end

function _write_location(io, system, id, dictionary, indentation)
    name = "loc" * string(id)
    _write_indented(io, "<location id=\"$id\" name=\"$name\">\n", indentation)
    indent!(indentation)
    _write_invariant(io, system, dictionary, indentation)
    _write_flow(io, system, dictionary, indentation)
    dedent!(indentation)
    _write_indented(io, "</location>\n", indentation)
end

function _write_invariant(io, system, dictionary, indentation)
    if !isconstrained(system)
        return  # nothing to write
    end

    _write_indented(io, "<invariant>", indentation)
    indent!(indentation)

    X = stateset(system)
    if !isnothing(X)
        _write_state_constraints_specific(io, system, X, dictionary, indentation)
    end

    U = inputset(system)
    if !isnothing(U)
        _write_input_constraints_specific(io, system, U, dictionary, indentation)
    end

    dedent!(indentation)
    write(io, "</invariant>\n")
end

function _write_state_constraints_specific(io, system, X, dictionary, indentation)
    println("WARNING: state constraints of type $(typeof(X)) are currenctly " *
        "not supported and will be ignored")
end

function _write_state_constraints_specific(io, system, X::Universe, dictionary,
                                           indentation)
    # nothing to write
end

function _write_state_constraints_specific(io, system, X::AbstractHyperrectangle,
                                           dictionary, indentation)
    for i in 1:dim(X)
        xi = _variable_name(i, dictionary)
        l = low(X, i)
        u = high(X, i)
        if i > 1
            write(io, " & ")
        end
        write(io, "$l <= $xi & $xi <= $u")
    end
end

function _write_state_constraints_specific(io, system,
                                           H::Union{HalfSpace, Hyperplane},
                                           dictionary, indentation)
    first = true
    for (i, ai) in enumerate(H.a)
        if ai == 0
            continue
        end
        xi = _variable_name(i, dictionary)
        if ai < 0
            if first
                prefix = "-"
            else
                prefix = " - "
            end
        else
            if first
                prefix = ""
            else
                prefix = " + "
            end
        end
        write(io, "$prefix$(abs(ai)) * $xi")
        first = false
    end
    if first
        # a = 0
        write(io, "0")
    end
    if H isa HalfSpace
        operator = "<="
    elseif H isa Hyperplane
        operator = "=="
    else
        error("unsupported set type $(typeof(H))")
    end
    b = H.b
    if b == 0
        b = 0  # turns -0.0 into 0 (for aesthetic reasons only)
    end
    write(io, " $operator $b")
end

function _write_state_constraints_specific(io, system, X::AbstractVector{<:LazySet},
                                           dictionary, indentation)
    first = true
    for (i, Xi) in enumerate(X)
        if first
            first = false
        else
            write(io, " & ")
        end
        _write_state_constraints_specific(io, system, Xi, dictionary, indentation)
    end
end

function _write_input_constraints_specific(io, system, U, dictionary, indentation)
    println("WARNING: input constraints of type $(typeof(U)) are currenctly " *
        "not supported and will be ignored")
end

function _write_input_constraints_specific(io, system, U::Universe, dictionary,
                                           indentation)
    # nothing to write
end

function _write_input_constraints_specific(io, system, U::AbstractHyperrectangle,
                                           dictionary, indentation)
    for i in 1:dim(U)
        ui = _input_name(i, dictionary)
        l = low(U, i)
        u = high(U, i)
        if i > 1
            write(io, " & ")
        end
        write(io, "$l <= $ui & $ui <= $u")
    end
end

function _write_flow(io, system, dictionary, indentation)
    _write_indented(io, "<flow>", indentation)
    _write_flow_specific(io, system, dictionary)
    write(io, "</flow>\n")
end

function _write_flow_specific(io,
                              system::Union{LinearContinuousSystem,
                                            LinearControlContinuousSystem,
                                            ConstrainedLinearContinuousSystem,
                                            ConstrainedLinearControlContinuousSystem,
                                            AffineContinuousSystem,
                                            AffineControlContinuousSystem,
                                            ConstrainedAffineContinuousSystem,
                                            ConstrainedAffineControlContinuousSystem},
                              dictionary)
    A = state_matrix(system)
    n = statedim(system)

    B = input_matrix(system)  # returns `nothing` for systems with no inputs
    m = inputdim(system)

    c = affine_term(system)
    if isnothing(c)
        c = zeros(eltype(A), n)
    end

    for i in 1:n
        xi = _variable_name(i, dictionary)
        write(io, "$xi' == ")
        first = true
        for j in 1:n
            Aij = A[i, j]
            if Aij == 0
                continue
            end
            xj = _variable_name(j, dictionary)
            if Aij < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(Aij)) * $xj")
            first = false
        end
        for j in 1:m
            Bij = B[i, j]
            if Bij == 0
                continue
            end
            uj = _input_name(j, dictionary)
            if Bij < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(Bij)) * $uj")
            first = false
        end
        if c[i] != 0
            if c[i] < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(c[i]))")
            first = false
        end
        if first
            # rhs was zero everywhere
            write(io, "0")
        end
        if i < n
            write(io, " & ")
        end
    end
end

function _write_transitions(io, system::AbstractContinuousSystem, dictionary,
                            indentation)
    # nothing to write
end

function _write_transitions(io, system::HybridSystem, dictionary, indentation)
    for transition in transitions(system)
        _write_transition(io, system, transition, dictionary, indentation)
    end
end

function _write_transition(io, H::HybridSystem, transition, dictionary, indentation)
    s = string(source(H, transition))
    t = string(target(H, transition))
    _write_indented(io, "<transition source=\"$s\" target=\"$t\">\n", indentation)
    indent!(indentation)
    _write_transition_label(io, H, transition, dictionary, indentation)
    _write_guard(io, H, transition, dictionary, indentation)
    _write_assignment(io, H, transition, dictionary, indentation)
    dedent!(indentation)
    _write_indented(io, "</transition>\n", indentation)
end

function _write_transition_label(io, H, transition, dictionary, indentation)
    # labels are ignored
end

function _write_guard(io, H, transition, dictionary, indentation)
    G = guard(H, transition)
    if !isconstrained(G)
        return  # nothing to write
    end
    _write_indented(io, "<guard>", indentation)
    _write_state_constraints_specific(io, system, G, dictionary, indentation)
    write(io, "</guard>\n")
end

function _write_assignment(io, H, transition, dictionary, indentation)
    asgn = resetmap(H, transition)  # workaround because transitions are parsed as system and not as map
    if isnothing(asgn)
        return  # nothing to write
    end

    n = statedim(H)
    if asgn isa AbstractMap
        if (asgn isa IdentityMap) || (asgn isa ConstrainedIdentityMap)
            return  # nothing to write
        elseif (asgn isa MathematicalSystems.LinearMap) || (asgn isa ConstrainedLinearMap)
            A = asgn.A
            c = zeros(n)
        elseif (asgn isa MathematicalSystems.AffineMap) || (asgn isa ConstrainedAffineMap)
            A = asgn.A
            c = asgn.c
        elseif (asgn isa MathematicalSystems.ResetMap) || (asgn isa ConstrainedResetMap)
            A = zeros(n, n)
            c = zeros(n)
            @inbounds for i in 1:n
                if i in keys(asgn.dict)
                    c[i] = asgn.dict[i]
                else
                    A[i, i] = 1.0
                end
            end
        else
            @warn("only linear assignments are supported at the moment")
        end
    else
        A = state_matrix(asgn)
        if !iszero(input_matrix(asgn))
            @warn("only linear assignments are supported at the moment")
        end
    end

    _write_indented(io, "<assignment>", indentation)
    for i in 1:n
        xi = _variable_name(i, dictionary)
        write(io, "$xi' == ")
        first = true
        for j in 1:n
            Aij = A[i, j]
            if Aij == 0
                continue
            end
            xj = _variable_name(j, dictionary)
            if Aij < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(Aij)) * $xj")
            first = false
        end
        if c[i] != 0
            if c[i] < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(c[i]))")
            first = false
        end
        if first
            # rhs was zero everywhere
            write(io, "0")
        end
        if i < n
            write(io, " & ")
        end
    end
    write(io, "</assignment>\n")
end


text:
using Bukdu # ApplicationController Conn pipeline routes get
using HTTP.Messages: setheader

struct RESTController <: ApplicationController
    conn::Conn
end

function index(c::RESTController)
    @info c.conn.request.response.headers
    render(JSON, 42)
end

pipeline(:cors) do conn::Conn
    setheader(conn.request.response, "Access-Control-Allow-Origin" => "*")
end

routes(:cors) do
    get("/", RESTController, index)
end

Bukdu.start(8080)


text:
function crossover(p_int1::Vector{Int64},p_float1::Matrix{Float64},
    p_int2::Vector{Int64},p_float2::Matrix{Float64})

    p_int3 = zeros(Int,size(p_int1,1)*size(p_int2,1))
    p_float3 = zeros((size(p_float1,1)*size(p_float2,1),size(p_float1,2)))

    add_on = 0

    min_index = min(size(p_int1,1),size(p_int2,1))
    max_index = max(size(p_int1,1),size(p_int2,1))

    for i in 1:min_index
        if rand() < 0.5
            p_int3[i] = p_int1[i]
            p_float3[i:i,:] = p_float1[i:i,:]
        else
            p_int3[i] = p_int2[i]
            p_float3[i:i,:] = p_float2[i:i,:]
        end
    end

    if min_index < size(p_int2,1)

        p_int3[min_index+1:max_index] = p_int2[min_index+1:max_index]
        p_float3[min_index+1:max_index,:] = p_float2[min_index+1:max_index,:]

        p_int3 = p_int3[1:max_index + add_on]
        p_float3 = p_float3[1:max_index + add_on,:]

    elseif min_index < size(p_int1,1)

        p_int3[min_index+1:max_index] = p_int1[min_index+1:max_index]
        p_float3[min_index+1:max_index,:] = p_float1[min_index+1:max_index,:]

        p_int3 = p_int3[1:max_index + add_on]
        p_float3 = p_float3[1:max_index + add_on,:]

    else
        p_int3 = p_int3[1:max_index + add_on]
        p_float3 = p_float3[1:max_index + add_on,:]
    end

    return p_int3,p_float3
end

function duplication_mutation(p_float::Matrix{Float64},p_int::Vector{Int64},p::Float64)
    if rand() < p
        idx = rand(1:size(p_int,1))
        c_int = [p_int;p_int[idx]]
        c_float = [p_float;p_float[idx:idx,:]]

        return c_int,c_float
    else
        return p_int,p_float
    end
end

function deletion_mutation(p_float::Matrix{Float64},p_int::Vector{Int64},save::Int,p::Float64)
    if rand() < p
        idx = rand(save+1:size(p_int,1))
        c_int = p_int[1:size(p_int,1) .!= idx]
        c_float = p_float[1:size(p_float,1) .!= idx,: ]

        return c_int,c_float
    else
        return p_int,p_float
    end
end

function creep_mutation!(p_float::Matrix{Float64},p::Float64)
    if rand() < p
        idx = rand(1:size(p_float,1))
        idx_val = rand(1:size(p_float,2))
        p_float[idx,idx_val] += randn()
    end
end

function copy_mutation!(p_float::Matrix{Float64},p::Float64)
    if rand() < p
        idx,donor = sample(1:size(p_float,1),2, replace = false)
        if (rand()<=0.5)
            p_float[idx:idx,1:3] = p_float[donor:donor,4:6]
        else
            p_float[donor:donor,1:3] = p_float[idx:idx,4:6]
        end
    end
end

function copy_mutation_env!(p_float::Matrix{Float64},p_int::Vector{Int64},p::Float64)
    if rand() < p
        idx = rand(1:size(p_float,1))
        donor = sample(findall(x->x==4, p_int),1)[1]

        if (rand()<=0.5)
            p_float[idx:idx,1:3] = p_float[donor:donor,4:6]
        else
            p_float[donor:donor,1:3] = p_float[idx:idx,4:6]
        end
    end
end

function copy_receptor_mutation!(p_float::Matrix{Float64},p_int::Vector{Int64},p::Float64)
    if rand() < p
        idx_receptor = rand(1:4)
        idx_env = sample(findall(x->x==4, p_int),1)[1]
        p_float[idx_env:idx_env,4:6] = p_float[idx_receptor:idx_receptor,1:3] # coding 1-3 promoter 4-6

    end
end

function type_mutation!(p_int::Vector{Int64},save::Int,p::Float64)
    if rand() < p
        idx = rand(save+1:size(p_int,1))
        p_int[idx] = rand(1:4)
    end
end

function create_child(p_int1::Vector{Int64},p_float1::Matrix{Float64},p_int2::Vector{Int64},p_float2::Matrix{Float64})

    c_int,c_float = crossover(p_int1,p_float1,p_int2,p_float2)

    return c_int,c_float
end

function create_child(p_int1::Vector{Int64},p_float1::Matrix{Float64},evo_config::Union{POOLConfig,EMAPConfig})

    n_mut = 7

    p = n_mut*rand()

    if p < 1
        c_int,c_float = duplication_mutation(p_float1,p_int1,1.)
    elseif p < 2
        c_int,c_float = deletion_mutation(p_float1,p_int1,evo_config.genome_preserve,1.)
    elseif p < 3
        c_int = copy(p_int1)
        c_float = copy(p_float1)
        creep_mutation!(c_float,1.)
    elseif p < 4
        c_int = copy(p_int1)
        c_float = copy(p_float1)
        copy_receptor_mutation!(c_float,c_int,1.)
    elseif p < 5
        c_int = copy(p_int1)
        c_float = copy(p_float1)
        copy_mutation!(c_float,1.)
    elseif p < 6
        c_int = copy(p_int1)
        c_float = copy(p_float1)
        copy_mutation_env!(c_float,c_int,1.)
    else
        c_int = copy(p_int1)
        c_float = copy(p_float1)
        type_mutation!(c_int,evo_config.genome_preserve,1.)
    end

    return c_int,c_float
end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


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


text:
"""
    SimpleSDMPredictor(CHELSA, BioClim, layer::Integer; left=nothing, right=nothing, bottom=nothing, top=nothing)

Download and prepare bioclim layers from the CHELSA database, and returns
them as an array of `SimpleSDMPredictor`s. Layers are called by their number,
from 1 to 19. The list of available layers is given in a table below.

The keyword argument is `path`, which refers to the path where the function
will look for the geotiff files.

Note that these files are *large* due the fine resolution of the data, and for
this reason this function will return the *integer* version of the layers. Also
note that the bioclim data are only available for the V1 of CHELSA, and are not
from the V2.

It is recommended to *keep* the content of the `path` folder, as it will
eliminate the need to download the tiff files (which are quite large). For
example, calling `bioclim(1:19)` will download and everything, and future
calls will be much faster.
"""
function SimpleSDMPredictor(::Type{CHELSA}, ::Type{BioClim}, layer::Integer=1; kwargs...)
    file = _get_raster(CHELSA, BioClim, layer)
    return geotiff(SimpleSDMPredictor, file; kwargs...)
end

function SimpleSDMPredictor(::Type{CHELSA}, ::Type{BioClim}, layers::AbstractArray; kwargs...)
    @assert eltype(layers) <: Integer
    return [SimpleSDMPredictor(CHELSA, BioClim, l; kwargs...) for l in layers]
end

"""
    SimpleSDMPredictor(CHELSA, BioClim, mod::CMIP5, fut::RepresentativeConcentrationPathway, layer::Integer=1; year="2041-2060", kwargs...)

Returns a *future* layer for a given RCP and model, at a given year (either
2041-2060 or 2061-2080), from the bioclim future climate.
"""
function SimpleSDMPredictor(::Type{CHELSA}, ::Type{BioClim}, mod::CMIP5, fut::RepresentativeConcentrationPathway, layer::Integer=1; year="2041-2060", kwargs...)
    @assert year in ["2041-2060", "2061-2080"]
    file = _get_raster(CHELSA, BioClim, mod, fut, layer, year)
    return geotiff(SimpleSDMPredictor, file; kwargs...)
end

function SimpleSDMPredictor(::Type{CHELSA}, ::Type{BioClim}, mod::CMIP5, fut::RepresentativeConcentrationPathway, layers::AbstractArray; kwargs...)
    @assert eltype(layers) <: Integer
    return [SimpleSDMPredictor(CHELSA, BioClim, mod, fut, l; kwargs...) for l in layers]
end

function SimpleSDMPredictor(::Type{CHELSA}, ::Type{BioClim}, mod::CMIP6, fut::SharedSocioeconomicPathway, layers::AbstractArray; kwargs...)
    @assert eltype(layers) <: Integer
    return [SimpleSDMPredictor(CHELSA, BioClim, mod, fut, l; kwargs...) for l in layers]
end

function SimpleSDMPredictor(::Type{CHELSA}, ::Type{BioClim}, mod::CMIP6, fut::SharedSocioeconomicPathway, layer::Integer=1; year="2041-2060", kwargs...)
    file = _get_raster(CHELSA, BioClim, mod, fut, layer, year)
    return geotiff(SimpleSDMPredictor, file; kwargs...)
end


text:
## MPICommunicator class
immutable MPI_Comm <: FEniCSObject
  pyobject::PyObject
end
mpi_comm_world() = MPI_Comm(dolfin.mpi_comm_world())
mpi_comm_self() = MPI_Comm(dolfin.mpi_comm_world())

export MPI_Comm, mpi_comm_world, mpi_comm_self


text:
# Basically an exact application of
# <https://github.com/JuliaPy/PyCall.jl#using-pycall-from-julia-modules>.
#
# The initialization strategy is similar to the one employed by SymPy.jl, see
# <https://github.com/JuliaPy/SymPy.jl/blob/aea9fb5568022bbd68fe4b69d670ce2a22413a7c/src/SymPy.jl#L124>.

__precompile__()

module Cclib

using PyCall

export bridge
export io
export method
export parser

export ccopen, ccwrite

"""
    bridge

Facilities for moving parsed data to other cheminformatic libraries.
"""
const bridge = PyCall.PyNULL()

"""
    io

Available writers for standard chemical representations.
"""
const io = PyCall.PyNULL()

"""
    method

Analyses and calculations based on data parsed by cclib.
"""
const method = PyCall.PyNULL()

"""
    parser

Parsers for all supported programs.
"""
const parser = PyCall.PyNULL()

"""
Guess the identity of a particular log file and return an instance of it.
"""
const ccopen = PyCall.PyNULL()

"""
Write the parsed data from an outputfile to a standard chemical representation.
"""
const ccwrite = PyCall.PyNULL()

function __init__()
    cclib = PyCall.pyimport_conda("cclib", "cclib")

    # Version of the underlying cclib Python package.
    global __version__ = cclib.__version__

    copy!(bridge, PyCall.pyimport_conda("cclib.bridge", "cclib"))
    copy!(io, PyCall.pyimport_conda("cclib.io", "cclib"))
    copy!(method, PyCall.pyimport_conda("cclib.method", "cclib"))
    copy!(parser, PyCall.pyimport_conda("cclib.parser", "cclib"))

    copy!(ccopen, io.ccopen)
    copy!(ccwrite, io.ccwrite)
end

end


text:
function settrue!(x::BitArray,b)
    fill!(x,false)
    x[b]=true
end

function settrue(N::Integer,b)
    x=BitArray(N)
    x[b]=true
    return x
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# Autogenerated wrapper script for Pangomm_jll for aarch64-apple-darwin
export libpangomm

using Cairomm_jll
using Glibmm_jll
using Pango_jll
JLLWrappers.@generate_wrapper_header("Pangomm")
JLLWrappers.@declare_library_product(libpangomm, "@rpath/libpangomm-2.48.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Cairomm_jll, Glibmm_jll, Pango_jll)
    JLLWrappers.@init_library_product(
        libpangomm,
        "lib/libpangomm-2.48.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SICMModelingToolkit

using ModelingToolkit, Unitful, DifferentialEquations, Optim
using RecursiveArrayTools
using SymbolicUtils
using SymbolicUtils: to_symbolic, substitute
using BenchmarkTools, LinearAlgebra, ForwardDiff, Zygote
using Quadrature, StaticArrays

include("math_base.jl")
include("lagrangians.jl")

export LocalTuple, AbstractLocalTuple, generate_generic_localtuple
export time, position, velocity, acceleration
export Lagrange_equations
export ode_problem_generator
export L_free_particle, L_harmonic_spring, L_kepler

end

max_stars_repo_name: zoemcc/SICMModelingToolkit

max_issues_repo_name: zoemcc/SICMModelingToolkit

max_forks_repo_name: zoemcc/SICMModelingToolkit


text:
using Robotlib
Tworld2base = eye(4)
# """
# `τ = gravity{P}(q::VecOrMat{P}, rm, m, dh::DH, Tworld2base = eye(4))`
# calculates the gravity torques given a joint angles `q`, (distances to center of mass × mass) `rm` and masses `m`. The DH-parameters for the robot must also be provided
# """
function gravity{P}(q::VecOrMat{P}, rm, m, dh::DH, Tworld2base = eye(4))
    Tn  = dh2Tn(dh,q)
    n   = size(Tn,3)-1
    # Find elements that are ≈ 0 and make them zero, smiplifies a lot!
    Tnz = Bool[round(Float64(evalf(subs(Tn[i,j,k],(q[1],1),(q[2],1),(q[3],1),(q[4],1),(q[5],1),(q[6],1),(q[7],1)))),10) == 0 for i=1:4, j=1:4, k=1:n+1]
    Tn[Tnz] = 0
    # Replace the numerical value for π/2 by the SymPi constant PI, this allows for nice simplifications
    Tn = Sym[subs(Tn[i,j,k],("1.5707963267948966",SymPy.PI))  for i=1:4, j=1:4, k=1:n+1]

    Ti = zeros(typeof(q[1]),4,4,n+1)
    Tn = cat(3,Tworld2base,Tn)
    Ti[:,:,1] = Tn[:,:,1]

    for i = 2:n
        Ti[:,:,i] = Ti[:,:,i-1]*Tn[:,:,i]
    end

    gv          = [0, 0, -9.82]
    i           = n
    gi          = Ti[1:3,1:3,i]'*gv
    τ           = Array(Sym,3,n)
    τhat        = Array(Sym,n)
    force       = Array(Sym,3,n)

    τ[:,i]      = skew(rm[:,i]) * gi
    force[:,i]  = gv*m[end]
    τhat[i]     = τ[3,i]

    for i = n-1:-1:1
        Ri      = Tn[1:3,1:3,i+1]
        gi      = Ti[1:3,1:3,i]'*gv # Local frame
        τ[:,i]  = skew(rm[:,i]) * gi + Ri*τ[:,i+1]
        τi      = τ[:,i]
        for k = (i+1):n # This loop formulations seems to be essential, I can not get it to work with accumulated forces
            rk = (trinv(Ti[:,:,i])*Ti[1:4,4,k])[1:3]
            τi += skew(rk)*Ti[1:3,1:3,i]'*force[:,k]
        end

        τhat[i]     = τi[3]
        force[:,i]  = gv*m[i]
    end
    return τhat
end

using Robotlib

using SymPy

# function symbolicGravity(dh::DH)
dh = DHYuMi()
n_joints = size(dh.dhpar,1)

rm  = Sym[symbols("rm$i$j",real=true) for i = 1:3, j=1:n_joints]
m   = Sym[symbols("m$j",real=true) for j=1:n_joints]
q   = Sym[symbols("q$j",real=true) for j=1:n_joints]
tau = Sym[symbols("tau$j",real=true) for j=1:n_joints]

baseAnglesLeft  = [-0.63 , 0.95 , -0.18]
Rbase           = rpy2R(baseAnglesLeft,"xyz")
Tbase           = eye(4)
Tbase[1:3,1:3]  = Rbase

tauhat          = gravity(q, rm, m, dh, Tbase);
w               = [rm[:];m[:]];

AA  = Array(Sym,n_joints, 4n_joints);
AAz = falses(size(AA))

eq  = tauhat;

println("Finding non-zero coefficients")
for i = 1:n_joints
    for j = 1:size(AA,2)
        AA[i,j] = (coeff(expand(eq[i]),w[j])) # This is a nice place to put simplify
        AAz[i,j] = abs(evalf(subs(AA[i,j],(q[1],1),(q[2],1),(q[3],1),(q[4],1),(q[5],1),(q[6],1),(q[7],1)))) < 1e-5
    end
    print(i)
end
println("")

# non_identifiable    = (sum(AA,1) .== 0)[:]
# # non_identifiable = all(AAz,1)[:]
# w                   = w[!non_identifiable];
# AA                  = AA[:,!non_identifiable];
n_joints,n_params   = size(AA)

println("Trying to find expressions common in all terms")
change,res  = cse(AA[:]);
res         = reshape(res,size(AA));

println("Printing results to file")
fid = open("gravityFridaLS.jl","w");
println(fid,"function gravityFridaLS(q)")
for i = 1:n_joints
    println(fid,"q$i = q[$i]")
end
for i = 0:length(change)-1
    sei = change[i+1][2]
    println(fid,"x$i = $sei")
end

println(fid,"A = zeros($n_joints,$n_params)")

for j = 1:n_joints
    for p = 1:n_params
        if res[j,p] != 0 #&& !AAz[j,p]
            println(fid,"A[$j,$p] = $(res[j,p])")
        end
    end
end
println(fid,"return A")
println(fid,"end")
close(fid)
println("Done")
# end


text:
module TiramisuClean

using ..ParallelAccelerator

export remove_unused!

Variable = Union{Slot, SSAValue}

function get_used_def(expr::Expr)
  # TODO support side effects, and nested assignments
  def = nothing
  if expr.head == :(=)
    def = expr.args[1]
  end
  used = union((u for (u, d) in map(get_used_def, expr.args))...)
  delete!(used, def)
  return used, Set{Variable}(def == nothing ? [] : [def])
end

function get_used_def(expr::Variable)
  return Set{Variable}([expr]), Set{Variable}()
end

function get_used_def(expr::Any)
  return Set{Variable}(), Set{Variable}()
end

function remove_unused!(ast::Expr)
  # println("Used, def")
  cleaned = false
  while !cleaned
    used = Set{Variable}()
    def = Set{Variable}()
    for expr in ast.args
      u, d = get_used_def(expr)
      union!(used, u)
      union!(def, d)
    end
    unused = setdiff(def, used)
    # println(unused)
    ast.args = [e for e in ast.args
                if !(isa(e, Expr) && e.head == :(=) && e.args[1] in unused)]
    cleaned = length(unused) == 0
  end
end

end