text:
## --- Hafnium isotopes

    # Calculate the initial Hf ratio and epsilon Hf at time t Ma
    function eHf(Hf176_Hf177, Lu176_Hf177, t; eHfOnly::Bool=true)

        # Lutetium decay constant (Soderlund et al., 2004
        lambda = 1.867E-11

        # Present-day CHUR composition (Bouvier et al., 2008)
        CHUR_Hf176_Hf177 = 0.282785
        CHUR_Lu176_Hf177 = 0.0336

        # Calculate initial Hf ratio at time t
        Hf176_Hf177_t = Hf176_Hf177 .- Lu176_Hf177.*(exp.(t .* 10^6*lambda) .- 1)

        # Calculate CHUR Hf ratio at time t
        CHUR_Hf176_Hf177_t = CHUR_Hf176_Hf177 .- CHUR_Lu176_Hf177.*(exp.(t .* 10^6*lambda) .- 1)

        # Calculate corresponding epsilon Hf
        eHf=(Hf176_Hf177_t ./ CHUR_Hf176_Hf177_t .- 1) .* 10^4

        if eHfOnly
            return eHf
        else
            return (eHf, Hf176_Hf177_t)
        end
    end
    export eHf

    function bin_bsr_eHf(x,Hf176_Hf177,Lu176_Hf177,age,min,max,nbins,x_sigma,Hf176_Hf177_sigma,Lu176_Hf177_sigma,age_sigma,nresamples)
        data = hcat(x,Hf176_Hf177,Lu176_Hf177,age)
        sigma = hcat(x_sigma,Hf176_Hf177_sigma,Lu176_Hf177_sigma,age_sigma)

        means = Array{Float64}(undef,nbins,nresamples)
        c = Array{Float64}(undef,nbins)
        for i=1:nresamples
            dbs = bsresample(data,sigma,length(age))
            eHf_resampled = eHf(dbs[:,2], dbs[:,3], dbs[:,4])
            (c,m,s) = binmeans(dbs[:,1], eHf_resampled, min, max, nbins)
            means[:,i] = m
        end

        m = nanmean(means,dim=2)
        el = m - nanpctile(means,2.5,dim=2)
        eu = nanpctile(means,97.5,dim=2) - m

        return (c, m, el, eu)
    end
    export bin_bsr_eHf

## ---


text:
using LinearAlgebra
using SparseArrays
using Plots

"This code solves Laplace's equation with Dirichlet boundary conditions using a
finite difference method."

m = 100 # number of points
x = LinRange(-1,1,m+2)
xint = x[2:end-1]

f(x) = 5*(Float64((x > -.5) && (x <= 0.0)) - Float64((x < .5) && (x >= 0.0)))
α,β = 1.0,pi

h = x[2]-x[1]
A = (1/h^2) * spdiagm(0=>2*ones(m),-1=>-ones(m-1),1=>-ones(m-1))
b = f.(xint)
b[1] += α/h^2
b[m] += β/h^2

# solve Au = F
uu = A\b
plot(xint,uu,linewidth=2,label="Solution")
plot!(x,uexact.(x),linestyle=:dash)
# plot!(x,.1*f.(x),linecolor=:red,ls=:dash,label="f(x)")
# plot!(leg=:topleft)

# uu - uexact.(xint)

ext: jl


text:
function extract(code::Union{Expr,Symbol,Number}, f::Symbol)

    if typeof(code) == Symbol || typeof(code) <: Number

        return nothing

    elseif code.head == :call && code.args[1] == f

        return code.args[2]

    else

        for i in code.args
            exp = extract(i, f)

            if exp !== nothing
                return exp
            end
        end

        return nothing

    end

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using GLMakie, DynamicalSystems, InteractiveDynamics
using OrdinaryDiffEq

ds = Systems.henonheiles()
diffeq = (alg = Vern9(),)
u0s = [
    [0.0, -0.25, 0.42081, 0.0],
    [0.0, 0.1, 0.5, 0.0],
    [0.0, -0.31596, 0.354461, 0.0591255]
]
trs = [trajectory(ds, 10000, u0; diffeq...)[:, SVector(1,2,3)] for u0 ∈ u0s]
for i in 2:length(u0s)
    append!(trs[1], trs[i])
end

# Inputs:
j = 2 # the dimension of the plane
tr = trs[1]

brainscan_poincaresos(tr, j)


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function words_from_le_bytes_8(
            bytes::AbstractVector{UInt8}
        )::SVector{8, UInt32}
    local data  = view(bytes, 1:(length(bytes) & ~0x3))
    local words = reinterpret(UInt32, data)
    return SVector{8, UInt32}(
        length(words) >= 1 ? ltoh(words[1]) : UInt32(0),
        length(words) >= 2 ? ltoh(words[2]) : UInt32(0),
        length(words) >= 3 ? ltoh(words[3]) : UInt32(0),
        length(words) >= 4 ? ltoh(words[4]) : UInt32(0),
        length(words) >= 5 ? ltoh(words[5]) : UInt32(0),
        length(words) >= 6 ? ltoh(words[6]) : UInt32(0),
        length(words) >= 7 ? ltoh(words[7]) : UInt32(0),
        length(words) >= 8 ? ltoh(words[8]) : UInt32(0)
    )
end

function parent_output(
            left_child_cv::SVector{8, UInt32},
            right_child_cv::SVector{8, UInt32},
            key::SVector{8, UInt32},
            flags::UInt32
        )::Output
    local block_words = vcat(left_child_cv, right_child_cv)

    return Output(key, block_words, 0, UInt32(BLOCK_LEN), PARENT | flags)
end

function parent_cv(
            left_child_cv::SVector{8, UInt32},
            right_child_cv::SVector{8, UInt32},
            key::SVector{8, UInt32},
            flags::UInt32
        )::SVector{8, UInt32}
    return chaining_value(parent_output(left_child_cv, right_child_cv, key, flags))
end

mutable struct Blake3Ctx
    chunk_state::ChunkState
    key::SVector{8, UInt32}
    cv_stack::Vector{SVector{8, UInt32}}
    cv_stack_len::UInt8
    flags::UInt32

    """
        Blake3Ctx(key::AbstractVector{UInt32}=IV, flags::UInt32=UInt32(0))

    Create a new hasher using the specified key.  The key must be 8 values.  Normally this
    method is called with the default values.

    # Examples
    ```
    h = Blake3Ctx()
    update!(h, rand(UInt8, 1024))
    r = digest(h)
    println("Hash => \$(bytes2hex(r))")
    ```
    """
    Blake3Ctx(
                key::AbstractVector{UInt32}=IV,
                flags::UInt32=UInt32(0)
            ) = begin
        length(key) != 8 && throw("Key length must be 8.")
        local skey = SVector{8, UInt32}(key)

        return new(
            ChunkState(skey, UInt64(0), flags),
            skey,
            Vector{SVector{8, UInt32}}(undef, 54),
            0,
            flags
        )
    end

    """
        Blake3Ctx(key::AbstractVector{UInt8})

    Create a Blake3Ctx instance to perform a keyed hash on the data.  Only the first 32 bytes
    are used.

    # Examples
    ```
    h = Blake3Ctx([x for x in UInt8(1):UInt8(32)])
    update!(h, rand(UInt8, 1024))
    r = digest(h)
    println("Hash => \$(bytes2hex(r))")
    ```
    """
    Blake3Ctx(key::AbstractVector{UInt8}) = Blake3Ctx(words_from_le_bytes_8(key), KEYED_HASH)

    """
        Blake3Ctx(content::AbstractString)

    Create a Blake3Ctx instance to perform a derived hash on the data.  The string is hashed
    and the result of the string is then used as the key for the hash.

    # Examples
    ```
    h = Blake3Ctx("My hash context")
    update!(h, rand(UInt8, 1024))
    r = digest(h)
    println("Hash => \$(bytes2hex(r))")
    ```
    """
    Blake3Ctx(content::AbstractString) = begin
        local context_hasher = Blake3Ctx(IV, DERIVE_KEY_CONTEXT)
        update!(context_hasher, Vector{UInt8}(content))

        local context_key = Vector{UInt8}(undef, KEY_LEN)
        digest(context_hasher, context_key)

        local context_key_words = words_from_le_bytes_8(context_key)
        return Blake3Ctx(context_key_words, DERIVE_KEY_MATERIAL)
    end
end

function push_stack(self::Blake3Ctx, cv::SVector{8, UInt32})::Nothing
    self.cv_stack_len += 1
    self.cv_stack[self.cv_stack_len] = cv
    nothing
end

function pop_stack(self::Blake3Ctx)::SVector{8, UInt32}
    self.cv_stack_len -= 1
    return self.cv_stack[self.cv_stack_len+1]
end

# Section 5.1.2 of the BLAKE3 spec explains this algorithm in more detail.
function add_chunk_chaining_value(
            self::Blake3Ctx,
            new_cv::SVector{8, UInt32},
            total_chunks::UInt64
        )::Nothing
    # This chunk might complete some subtrees. For each completed subtree,
    # its left child will be the current top entry in the CV stack, and
    # its right child will be the current value of `new_cv`. Pop each left
    # child off the stack, merge it with `new_cv`, and overwrite `new_cv`
    # with the result. After all these merges, push the final value of
    # `new_cv` onto the stack. The number of completed subtrees is given
    # by the number of trailing 0-bits in the new total number of chunks.
    while total_chunks & 1 == 0
        new_cv = parent_cv(pop_stack(self), new_cv, self.key, self.flags)
        total_chunks = total_chunks >> 1;
    end
    push_stack(self, new_cv)
    nothing
end

"""
    update!(self::Blake3Ctx, input::AbstractVector{UInt8}, datalen=length(input))::Nothing

Hash more data.  This method can be called multiple times until all the data has been
hashed.
"""
function update!(
            self::Blake3Ctx,
            input::AbstractVector{UInt8},
            datalen=length(input)
        )::Nothing
    local output = nothing
    local offset = 1

    while offset <= datalen
        offset, output = update(self.chunk_state, input, offset)
        if output !== nothing
            local chunk_cv     = chaining_value(output)
            local total_chunks = self.chunk_state.chunk_counter + 1
            add_chunk_chaining_value(self, chunk_cv, total_chunks)
            reset(self.chunk_state, self.key, total_chunks, self.flags)
        end
    end
end

"""
    digest(self::Blake3Ctx, out_slice::AbstractArray{UInt8})::Nothing

Generate the hash and write it to the provided array.  The contents of out_slice will be
overwritten.

__Details:__

Normally you would pass in an array of 32 bytes which the method will fill with the hash.
However if you want more (or less) bytes you can pass in a an array of a different size.
"""
function digest(self::Blake3Ctx, out_slice::AbstractArray{UInt8})::Nothing
    # Starting with the Output from the current chunk, compute all the
    # parent chaining values along the right edge of the tree, until we
    # have the root Output.
    local out                    = complete(self.chunk_state)
    local parent_nodes_remaining = UInt(self.cv_stack_len)

    while parent_nodes_remaining > 0
        out = parent_output(
            self.cv_stack[parent_nodes_remaining],
            chaining_value(out),
            self.key,
            self.flags
        )
        parent_nodes_remaining -= 1
    end
    root_output_bytes(out, out_slice)
end

"""
    digest(self::Blake3Ctx, bytes=32)::Vector{UInt8}

Generate the hash and return it.

__Details:__

Normally you would want 32 bytes back (256bits) but you can have it return more or less
bytes based on your need.
"""
function digest(self::Blake3Ctx, bytes=32)::Vector{UInt8}
    local output = Vector{UInt8}(undef, bytes)
    digest(self, output)
    return output
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
Formulation described in:

    Damcı-Kurt, P., Küçükyavuz, S., Rajan, D., & Atamtürk, A. (2016). A polyhedral
    study of production ramping. Mathematical Programming, 158(1), 175-205.
    DOI: https://doi.org/10.1007/s10107-015-0919-9
"""
module DamKucRajAta2016

import ..RampingFormulation

struct Ramping <: RampingFormulation end

end


text:
module FDTD
    using LinearAlgebra
    using Strided
    using DiscreteAxis
    using ProgressMeter
    const EHTuple{T,N} = Tuple{Tuple{Array{T,N},Array{T,N},Array{T,N}},Tuple{Array{T,N},Array{T,N},Array{T,N}}}
    const VecArray{T,N} = Tuple{Array{T,N},Array{T,N},Array{T,N}}

    include("pde_utils.jl")
    include("yee.jl")
    include("fdtd_utils.jl")

    include("update_equations.jl")
    include("cpml.jl")
    include("propagation.jl")

    export FDTD_propagate, radar_propagate, setup_spacetime
    export EMField, VecField3, Medium, Coefficients, PML, pontying, power_density, vec_array
    export PMLAux
    export ε₀, μ₀, η₀, c₀
end


text:
using SeaPearl
using SeaPearlExtras
using ReinforcementLearning
const RL = ReinforcementLearning
using Flux
using GeometricFlux
using BSON: @save, @load
using Random
using Statistics

include("rewards.jl")
include("features.jl")

# -------------------
# Generator
# -------------------
knapsack_generator = SeaPearl.KnapsackGenerator(10, 10, 0.2)

# -------------------
# Internal variables
# -------------------
StateRepresentation = SeaPearl.DefaultStateRepresentation{KnapsackFeaturization, SeaPearl.DefaultTrajectoryState}
numInFeatures = SeaPearl.feature_length(StateRepresentation)

# -------------------
# Experience variables
# -------------------
nbEpisodes = 1000
evalFreq = 100
nbInstances = 3
nbRandomHeuristics = 0

# -------------------
# Agent definition
# -------------------
include("agents.jl")

# -------------------
# Value Heuristic definition
# -------------------
learnedHeuristic = SeaPearl.LearnedHeuristic{StateRepresentation, knapsackReward, SeaPearl.FixedOutput}(agent)
basicHeuristic = SeaPearl.BasicHeuristic((x; cpmodel=nothing) -> SeaPearl.maximum(x.domain)) # Basic value-selection heuristic

# -------------------
# Variable Heuristic definition
# -------------------
struct KnapsackVariableSelection <: SeaPearl.AbstractVariableSelection{false} end

function (::KnapsackVariableSelection)(model::SeaPearl.CPModel)
    i = 1
    while SeaPearl.isbound(model.variables["x[" * string(i) * "]"])
        i += 1
    end
    return model.variables["x[" * string(i) * "]"]
end

valueSelectionArray = [learnedHeuristic, basicHeuristic]

# -------------------
# -------------------
# Core function
# -------------------
# -------------------
function trytrain(nbEpisodes::Int)

    metricsArray, eval_metricsArray = SeaPearl.train!(;
        valueSelectionArray= valueSelectionArray,
        generator=knapsack_generator,
        nbEpisodes=nbEpisodes,
        strategy=SeaPearl.DFSearch(),
        variableHeuristic=KnapsackVariableSelection(),
        out_solver=false,
        verbose=false, #true to print processus
        evaluator=SeaPearl.SameInstancesEvaluator(valueSelectionArray,knapsack_generator; evalFreq=evalFreq, nbInstances=nbInstances),
        restartPerInstances = 1
        )

    return metricsArray, eval_metricsArray
end

metricsArray, eval_metricsArray = trytrain(nbEpisodes)
nothing

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using .ForwardDiff: gradient!
using .DiffResults: GradientResult, value, gradient

struct MALA{D} <: MHSampler
    proposal::D
end

# Create a RandomWalkProposal if we weren't given one already.
MALA(d) = MALA(RandomWalkProposal(d))

# If we were given a RandomWalkProposal, just use that instead.
MALA(d::RandomWalkProposal) = MALA{typeof(d)}(d)

struct GradientTransition{T<:Union{Vector, Real, NamedTuple}, L<:Real, G<:Union{Vector, Real, NamedTuple}} <: AbstractTransition
    params::T
    lp::L
    gradient::G
end

logdensity(model::DensityModel, t::GradientTransition) = t.lp

propose(rng::Random.AbstractRNG, ::MALA, model) = error("please specify initial parameters")
function transition(sampler::MALA, model::DensityModel, params)
    return GradientTransition(params, logdensity_and_gradient(model, params)...)
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::DensityModel,
    sampler::MALA,
    transition_prev::GradientTransition;
    kwargs...
)
    # Extract value and gradient of the log density of the current state.
    state = transition_prev.params
    logdensity_state = transition_prev.lp
    gradient_logdensity_state = transition_prev.gradient

    # Generate a new proposal.
    proposal = sampler.proposal
    candidate = propose(rng, proposal(gradient_logdensity_state), model, state)

    # Compute both the value of the log density and its gradient
    logdensity_candidate, gradient_logdensity_candidate = logdensity_and_gradient(
        model, candidate
    )

    # Compute the log ratio of proposal densities.
    logratio_proposal_density = q(
        proposal(-gradient_logdensity_candidate), state, candidate
    ) - q(proposal(-gradient_logdensity_state), candidate, state)

    # Compute the log acceptance probability.
    logα = logdensity_candidate - logdensity_state + logratio_proposal_density

    # Decide whether to return the previous params or the new one.
    transition = if -Random.randexp(rng) < logα
        GradientTransition(candidate, logdensity_candidate, gradient_logdensity_candidate)
    else
        transition_prev
    end

    return transition, transition
end

"""
    logdensity_and_gradient(model::DensityModel, params)

Return the value and gradient of the log density of the parameters `params` for the `model`.
"""
function logdensity_and_gradient(model::DensityModel, params)
    res = GradientResult(params)
    gradient!(res, model.logdensity, params)
    return value(res), gradient(res)
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto
#see LICENSE.txt
#this work was supported in part by DARPA Contract D15PC00135
#unum-bitwalk.jl

#implements a "bitwalking" functional.  Said functional takes a ulp that isn't
#at maximal fraction length and then breaks it into two ulps that has one extra
#bit of length.  Then it executes the binary function for the ulps and the exacts
#if desired.  Returns an array with the unums that evaluate to true.

function bitwalk{ESS,FSS}(bf, u::Unum{ESS,FSS}, exacts = false, bound_exacts = false)
  is_exact(u) && return [u]
  u.fsize == max_fsize(FSS) && throw(ArgumentError("argument cannot be at max fsize"))
  u_neg = is_negative(u)
  #make sure outer_exact is overridden by the exact parameter.
  bound_exacts = exacts && bound_exacts
  #generate the result array.
  res = Utype[]
  #create a subsidiary function that pushes onto the results array based on the
  #sign of u.

  exact_inner = Unum{ESS,FSS}(UInt16(u.fsize + 1), u.esize, u.flags & UNUM_SIGN_MASK, u.fraction, u.exponent)
  bound_exacts && bf(exact_inner) && (res = vcat(res, exact_inner))

  #generate the inner ulp:  This is just flipping the ubit mask on the inner value.
  ulp_inner = unum_unsafe(exact_inner, exact_inner.flags | UNUM_UBIT_MASK)
  bf(ulp_inner) && (res = vcat(res, ulp_inner))

  #calculate the exact middle.
  new_fraction = u.fraction | __bit_from_top(u.fsize + 2, __frac_cells(FSS))
  exact_middle = Unum{ESS,FSS}(UInt16(u.fsize + 1), u.esize, u.flags & UNUM_SIGN_MASK, new_fraction, u.exponent)
  exacts && bf(exact_middle) && (res = vcat(res, exact_middle))

  #generate the outer ulp:  This is just flipping the ubit mask on the middle value
  ulp_outer = unum_unsafe(exact_middle, exact_middle.flags | UNUM_UBIT_MASK)
  bf(ulp_outer) && (res = vcat(res, ulp_outer))

  if bound_exacts
    exact_outer = __outward_exact(ulp_outer)
    bf(exact_outer) && (res = vcat(res, exact_outer))
  end

  res
end

export bitwalk

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Find Patterns Forming Clumps in a String

const NUCLEOTIDES = ('A', 'C', 'G', 'T')
const NUCLEOTIDE_ORDER = Dict{Char, Int64}('A' => 0, 'C' => 1, 'G' => 2, 'T' => 3)

convert_to_intlist(line) = map(x -> parse(Int, x), split(line))

pattern_to_number(pattern) = foldl((acc, nucleotide) -> 4 * acc + NUCLEOTIDE_ORDER[nucleotide], pattern; init=0)

function number_to_pattern(encoding, k)
    pattern = Char[]
    for _ in 1:k
        encoding, remainder = divrem(encoding, 4)
        push!(pattern, NUCLEOTIDES[remainder+1])
    end
    join(reverse(pattern))
end

function compute_frequencies(text, k)
    frequencies = zeros(Int32, 4 ^ k)
    for ix in 1:(length(text)-k+1)
        pattern = text[ix:ix+k-1]
        code = pattern_to_number(pattern)
        frequencies[code+1] += 1
    end
    frequencies
end

function collect_clump_forming_k_mers(clump_array, k)
    k_mers_in_clump = String[]
    for (ix, isclump) in enumerate(clump_array)
        if isclump
            k_mer = number_to_pattern(ix - 1, k)
            push!(k_mers_in_clump, k_mer)
        end
    end
    k_mers_in_clump
end

function find_clumps_in_text(genome, k, l, t)
    clump_array = falses(4 ^ k)
    frequencies = compute_frequencies(genome[1:l], k)
    for (ix, freq) in enumerate(frequencies)
        if freq >= t
            clump_array[ix] = true
        end
    end

    for ix in 2:(length(genome)-l)
        start_k_mer = genome[ix-1:ix-2+k]
        start_code = pattern_to_number(start_k_mer)
        frequencies[start_code+1] -= 1
        end_k_mer = genome[ix+l-k:ix+l-1]
        end_code = pattern_to_number(end_k_mer)
        frequencies[end_code+1] += 1
        if frequencies[end_code+1] >= t
            clump_array[end_code+1] = true
        end
    end

    collect_clump_forming_k_mers(clump_array, k)
end

function main()
    genome = readline()
    k, l, t = convert_to_intlist(readline())
    result = find_clumps_in_text(genome, k, l, t)
    println(join(result, " "))
end

main()


text:
using DFControl

# CdTe = create_job("CdTe",phd_dir*"CdTe/NSOC",default_inputs[:scf],default_inputs[:bands],server_dir = "CdTe/NSOC")
cdte = load_job(phd_dir*"CdTe/NSOC")
print_data(cdte)
outputs = pull_outputs(cdte)
out = read_qe_output(outputs[6])[:polarization]
using Plots
fermi = read_qe_output(outputs[1])[:fermi]
plot(read_qe_output(outputs[2])[:bands],ylim=[0,5])
size(read_qe_output(outputs[1])[:cell_parameters])
change_atoms!(cdte,read_qe_output(outputs[1])[:atomic_positions])
change_cell_parameters!(cdte,read_qe_output(outputs[1])[:cell_parameters])
add_calculation!(cdte,deepcopy(get_input(cdte,"scf")),1,filename="vcrel.in")
change_flow!(cdte,"vc"=>false)
change_flags!(cdte,"vcrel.in",:calculation => "'vc-relax'")
change_data!(cdte,["nscf1.in","nscf2.in","nscf3.in"],:k_points,[6,6,6,1,1,1])
add_flags!(cdte,["nscf1.in","nscf2.in","nscf3.in"],:control,:nppstr=>6)
change_data_option!(cdte,["nscf1.in","nscf2.in","nscf3.in"],:k_points,:automatic)
print_flags(cdte,"nscf3")
remove_flags!(cdte,Symbol("starting_magnetization(2)"))
cdte.calculations[end].filename = "nscf2.in"
add_flags!(cdte,"nscf1.in",:control,:lberry=>true,:gdir=>1,:nppstr=>12)
add_calculation!(cdte,deepcopy(get_input(cdte,"nscf1")),filename="nscf3.in")
add_flag
print_info(cdte)
change_flags!(cdte,:nppstr => 10)
change_flags!(cdte,"nscf3.in",:gdir => 1)
add_flags!(cdte,"nscf1.in",:control,:nppstr => 10)
change_data!(cdte,["nscf3.in"],:k_points, gen_k_grid(6,6,6,:nscf))
print_flow(cdte)
cdte.server_dir = "CdTe/NSOC/"
print_flow(cdte)
change_flow!(cdte,"nscf1.in"=>true,"nscf2.in"=>true,"scf"=>true)
change_data_option!(cdte,"nscf",:k_points,:crystal)
remove_flags!(cdte,"nscf3.in",:nppstr)
change_flags!(cdte,:occupations=>"'fixed'",:degauss=>0.0f0)
atoms = Dict(:Te => [Point3D{Float32}(0.3333333,0.6666667,0.40000),Point3D{Float32}(0.6666667,.3333333,0.87500)],:Cd =>[Point3D{Float32}(0.3333333,0.6666667,0.0000000),Point3D{Float32}(0.6666667,0.3333333,0.5000000)])
change_atoms!(cdte,atoms,pseudo_set=:pbesol,pseudo_fuzzy="paw")
change_data!(CdTe,"bands",:k_points,[[0.5f0,0.5f0,0.0f0,100],[0.0f0,0.0f0,0.0f0,100f0],[0.5f0,0.0f0,0.0f0,1]])
cell = Float32[0.8660254  -0.5000000   0.0000000;0.0000000   1.0000000   0.0000000;0.0000000   0.0000000   1.6367615]*Float32[4.5700000 0.0 0.0 ;0.0 4.5700000 0.0;0.0 0.0 7.4800000]
change_cell_parameters!(CdTe,cell)
remove_flags!(CdTe,:nbnd)
add_flags!(CdTe,:system,:A => 1.0f0)
change_data_option!(CdTe,:cell_parameters,:alat)
replace_header_word!(cdte,"frontend","defpart")
cdte.server_dir = "CdTe/NSOC/"
change_flow!(cdte,"nscf2.in"=>true)
change_flags!(cdte,:ecutwfc=>35)
submit_job(cdte)
cdte.server_dir

outputs = pull_outputs(CdTe,extras=["*s_j0.5*"])
using  Plots
plot(read_qe_bands_file(outputs[2]),fermi=read_fermi_from_qe_file(outputs[1]),ylims=[-5,5],legend=false)
set_default_input(:scf,get_input(CdTe,"scf"))
set_default_input(:bands,get_input(CdTe,"bands"))

set_default_input(:projwfc,read_qe_input(DFControl.assets_dir*"inputs/qe/projwfc.in",run_command="mpirun -np 24 ~/bin/projwfc.x"))
add_calculation!(CdTe,default_inputs[:projwfc])
change_flow!(CdTe,[("scf",false),("bands",false)])
heatmap(read_qe_kpdos(outputs[3])[1])

CdTe_soc = deepcopy(CdTe)
orig_atoms = Dict(:Te => [Point3D{Float32}(0.3333333,0.6666667,0.375000),Point3D{Float32}(0.6666667,.3333333,0.87500)],:Cd =>[Point3D{Float32}(0.3333333,0.6666667,0.0000000),Point3D{Float32}(0.6666667,0.3333333,0.5000000)])
new_atoms = copy(orig_atoms)
for x=0.1:0.1:0.9
  for (key,val) in orig_atoms
    if key == :Te
      for i=1:length(val)
        new_atoms[key][i] = val[i] .+ Point3D{Float32}(0.0,0.0,x*0.06)
      end
    end
  end
  change_atoms!(cdte,new_atoms,pseudo_set=:pbesol,pseudo_fuzzy="paw")
  cdte.server_dir = "CdTe/NSOC/berryz$x/"
  cdte.local_dir  = phd_dir*"CdTe/NSOC/berryz$x/"
  submit_job(cdte)
end
cdte.local_dir  = phd_dir*"CdTe/NSOC/"
cdte.server_dir = "CdTe/NSOC/"
outputs = pull_outputs(cdte)
polarizations = Array{Tuple{Point3D{Float64},Float64},1}()
for x=0.1:0.1:0.8
  cdte.local_dir  = phd_dir*"CdTe/NSOC/berryz$x/"
  cdte.server_dir = "CdTe/NSOC/berryz$x/"
  files = pull_outputs(cdte)[end-2:end]
  t     = Point3D{Float64}()
  mod   = zero(Float64)
  for f in files
    t_,mod =read_qe_polarization(f,Float64)
    t+=t_
  end
  push!(polarizations,(t,mod))
end
plot(getfield.(getindex.(polarizations,1),:x))
add_flags!(CdTe_soc,:system,:lspinorb=>true,:noncolin=>tru
CdTe_soc.server_dir = "CdTe/SOC/"
CdTe_soc.local_dir  = phd_dir*"CdTe/SOC/"
change_flow!(cdte_soc,"scf"=>false)
submit_job(cdte_soc)
cdte_soc = load_server_job("CdTe/SOC/",phd_dir*"CdTe/SOC/")
print_data(cdte_soc,"bands")
change_atoms!(cdte,orig_atoms)
outputs = pull_outputs(cdte_soc)
plot(read_qe_bands_file(outputs[1]),legend=false,ylims=[10,25])
print_flow(cdte_soc)
submit_job(cdte)

vcrel_cell = Float32[3.550867828  -2.049403922  -0.000000346;0.000598100   4.099843784   0.000000349;-0.000000628   0.000001104  10.094658435]
vcrel_atoms = Dict(:Te => [Point3D{Float32}(0.333293113,0.666706909,0.379372877),Point3D{Float32}(0.666706909,0.333293083,0.879374290)],:Cd =>[Point3D{Float32}(0.333291511,0.666708506,-0.004374027),Point3D{Float32}(0.666708468,0.333291502,0.495626861)])
change_atoms!(cdte,vcrel_atoms)
change_cell_parameters!(cdte,vcrel_cell)

cdte_soc = load_server_job("CdTe/SOC",phd_dir*"CdTe/SOC")
change_flags!(cdte_soc,:verbosity=>"'high'")
change_atoms!(cdte_soc,read_qe_output(outputs[1])[:atomic_positions])
change_cell_parameters!(cdte_soc,read_qe_output(outputs[1])[:cell_parameters])
print_flow(cdte_soc)
submit_job(cdte_soc)


text:
module Pass
using Cassette
using Core: CodeInfo, SlotNumber, SSAValue

Cassette.@context Ctx

function Cassette.execute(ctx::Ctx, callback, f, args...)
    if Cassette.canoverdub(ctx, f, args...)
        _ctx = Cassette.similarcontext(ctx, metadata = callback)
        return Cassette.overdub(_ctx, f, args...) # return result, callback
    else
        return Cassette.fallback(ctx, f, args...), callback
    end
end

function Cassette.execute(ctx::Ctx, callback, ::typeof(println), args...)
    return nothing, () -> (callback(); println(args...))
end

function sliceprintln(::Type{<:Ctx}, ::Type{S}, ir::CodeInfo) where {S}
    @show callbackslotname = gensym("callback")
    @show push!(ir.slotnames, callbackslotname)
    @show push!(ir.slotflags, 0x00)
    @show callbackslot = SlotNumber(length(ir.slotnames))
    @info "getmetadata"
    getmetadata = Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), Expr(:contextslot), QuoteNode(:metadata))

    # insert the initial `callbackslot` assignment into the IR.
    @info "insert"
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                 (stmt, i) -> i == 1 ? 2 : nothing,
                                 (stmt, i) -> [Expr(:(=), callbackslot, getmetadata), stmt])

    # replace all calls of the form `f(args...)` with `callback(f, args...)`, taking care to
    # properly destructure the returned `(result, callback)` into the appropriate statements
    @info "insert more"
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                 (stmt, i) -> begin
                                    i > 1 || return nothing # don't slice the callback assignment
                                    stmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                                    return Base.Meta.isexpr(stmt, :call) ? 3 : nothing
                                 end,
                                 (stmt, i) -> begin
                                     items = Any[]
                                     callstmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                                     push!(items, Expr(:call, callbackslot, callstmt.args...))
                                     push!(items, Expr(:(=), callbackslot, Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), SSAValue(i), 2)))
                                     result = Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), SSAValue(i), 1)
                                     if Base.Meta.isexpr(stmt, :(=))
                                         result = Expr(:(=), stmt.args[1], result)
                                     end
                                     push!(items, result)
                                     return items
                                 end)

    @info "insert a third time"
    # replace return statements of the form `return x` with `return (x, callback)`
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                  (stmt, i) -> Base.Meta.isexpr(stmt, :return) ? 2 : nothing,
                                  (stmt, i) -> begin
                                @info "lambda"
                                @show ex = Expr(:nooverdub, GlobalRef(Core, :tuple))
                                @show a = ex, stmt.args[1]
                                @show typeof(a)
                                @info "SSAV"
                                @show sv = SSAValue(i)
                                @info "return"
                                @show out = [
                                          Expr(:call, a, callbackslot)
                                          Expr(:return, sv)
                                      ]
                                return out
                                end)
    @info "returning $ir"
    return ir
end

const sliceprintlnpass = Cassette.@pass sliceprintln

a = rand(3)
b = rand(3)
function add(a, b)
    # println("I'm about to add $a + $b")
    c = a + b
    # println("c = $c")
    return c
end
add(a, b)
ctx = Ctx(pass=sliceprintlnpass, metadata = () -> nothing);

result, callback = Cassette.overdub(ctx, add, a, b)

@info "calling the callback"
# callback()

end


text:
"""
    function addEventDivision_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions for division events.
"""
function addEventDivision_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if "EventDivision" in keys(abm.declaredUpdates)

        code = abm.declaredUpdates["EventDivision"]

        if !@capture(code,if v_ g__ end)
            error("Erroneos structure of division events. Should be if something update end")
        end

        code = vectorize_(abm,code,p)
        subcode = postwalk(x->@capture(x,if v_ g__ end) ? quote $(g...) end : x, code)
        subcode = unblock(subcode)
        condition = postwalk(x->@capture(x,if v_ g__ end) ? v : x, code)
        condition = unblock(condition)

        #Add atomic_add
        if platform == "cpu"
            pushfirst!(subcode.args, :(ic1New_ = Threads.atomic_add!(NV,Int(1)) + 1))
            pushfirst!(subcode.args, :(idNew_ = Threads.atomic_add!(agentIdMax,Int(2)) + 1))
        elseif platform == "gpu"
            pushfirst!(subcode.args, :(ic1New_ = CUDA.atomic_add!(CUDA.pointer(NV,1),Int32(1)) + 1))
            pushfirst!(subcode.args, :(idNew_ = CUDA.atomic_add!(CUDA.pointer(agentIdMax,1),Int32(2)) + 1))
        end
        for k in ["Local","Identity"]
            vec = Meta.parse(string(lowercase(k),"V"))
            vecCopy = Meta.parse(string(lowercase(k),"VCopy"))
            for (i,j) in enumerate(abm.declaredSymbols[k])
                if j in keys(p.update["EventDivision"])
                    ii = p.update[k][j]
                    #Parse content
                    s = Meta.parse(string(j,"_1"))
                    subcode = postwalk(x -> @capture(x,g_=v__) && g == s ? :($vecCopy[ic1_,$ii] = $(v...)) : x, subcode)
                    s = Meta.parse(string(j,"_2"))
                    subcode = postwalk(x -> @capture(x,g_=v__) && g == s ? :($vecCopy[ic1New_,$ii] = $(v...)) : x, subcode)
                end
            end
        end
        for k in ["Local","Identity"]  #Add other updates of non changed variables
            vec = Meta.parse(string(lowercase(k),"V"))
            vecCopy = Meta.parse(string(lowercase(k),"VCopy"))
            for (i,j) in enumerate(abm.declaredSymbols[k])
                #Add update if doesn't exist
                if !(j in keys(p.update["EventDivision"])) && j != :agentId
                    push!(subcode.args,:($vec[ic1New_,$i] = $vec[ic1_,$i]))
                elseif j != :agentId
                    ii = p.update[k][j]
                    push!(subcode.args,:($vec[ic1_,$i] = $vecCopy[ic1_,$ii]))
                    push!(subcode.args,:($vec[ic1New_,$i] = $vecCopy[ic1New_,$ii]))
                else
                    ii = findfirst(abm.declaredSymbols["Identity"] .== j)
                    push!(subcode.args,:($vec[ic1_,$ii] = idNew_))
                    push!(subcode.args,:($vec[ic1New_,$ii] = idNew_ + 1))
                end
            end
        end
        code = quote
            if $condition
                $subcode
            end
        end

        code = simpleFirstLoopWrapInFunction_(platform,:division_!,code)

        push!(p.declareF.args,code)
        push!(p.args,:(NV))
        if platform == "cpu"
            push!(p.declareVar.args,:(NV = Threads.Atomic{Int}(N)))
            push!(p.declareVar.args,:(agentIdMax = Threads.Atomic{Int}(N)))

            push!(p.execInloop.args,
                    :(begin
                        @platformAdapt division_!(ARGS_)
                        N = NV[]
                    end)
                )
        elseif platform == "gpu"
            push!(p.declareVar.args,:(NV = N .*CUDA.ones(Int,1)))
            push!(p.declareVar.args,:(agentIdMax = N .*CUDA.ones(Int,1)))

            push!(p.execInloop.args,
                    :(begin
                        @platformAdapt division_!(ARGS_)
                        N = Core.Array(NV)[1]
                    end)
                )
        end

    end

    return nothing
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libim_convert_/include/im_convert.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libim_convert_ = "im") : (const libim_convert_ = "libim")  # Name of IM shared lib.

# begin enum imComplex2Real
const imComplex2Real = UInt32
const IM_CPX_REAL = 0
const IM_CPX_IMAG = 1
const IM_CPX_MAG = 2
const IM_CPX_PHASE = 3
# end enum imComplex2Real

# begin enum imGammaFactor
const imGammaFactor = Cint
const IM_GAMMA_LINEAR = 0
const IM_GAMMA_LOGLITE = -10
const IM_GAMMA_LOGHEAVY = -1000
const IM_GAMMA_EXPLITE = 2
const IM_GAMMA_EXPHEAVY = 7
# end enum imGammaFactor

# begin enum imCastMode
const imCastMode = UInt32
const IM_CAST_MINMAX = 0
const IM_CAST_FIXED = 1
const IM_CAST_DIRECT = 2
const IM_CAST_USER = 3
# end enum imCastMode

function imConvertDataType(src_image::Ptr{imImage}, dst_image::Ptr{imImage}, cpx2real::Cint, gamma::Cfloat, abssolute::Cint, cast_mode::Cint)
  ccall( (:imConvertDataType, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}, Cint, Cfloat, Cint, Cint), src_image, dst_image, cpx2real, gamma, abssolute, cast_mode)
end

function imConvertColorSpace(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imConvertColorSpace, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end

function imConvertToBitmap(src_image::Ptr{imImage}, dst_image::Ptr{imImage}, cpx2real::Cint, gamma::Cfloat, abssolute::Cint, cast_mode::Cint)
  ccall( (:imConvertToBitmap, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}, Cint, Cfloat, Cint, Cint), src_image, dst_image, cpx2real, gamma, abssolute, cast_mode)
end

function imImageGetOpenGLData(image::Ptr{imImage}, glformat::Ptr{Cint})
  ccall( (:imImageGetOpenGLData, libim_convert_), Ptr{Cvoid}, (Ptr{imImage}, Ptr{Cint}), image, glformat)
end

function imImageCreateFromOpenGLData(width::Cint, height::Cint, glformat::Cint, gldata::Ptr{Cvoid})
  ccall( (:imImageCreateFromOpenGLData, libim_convert_), Ptr{imImage}, (Cint, Cint, Cint, Ptr{Cvoid}), width, height, glformat, gldata)
end

function imConvertPacking(src_data::Ptr{Cvoid}, dst_data::Ptr{Cvoid}, width::Cint, height::Cint, src_depth::Cint, dst_depth::Cint, data_type::Cint, src_is_packed::Cint)
  ccall( (:imConvertPacking, libim_convert_), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint, Cint, Cint, Cint, Cint), src_data, dst_data, width, height, src_depth, dst_depth, data_type, src_is_packed)
end

function imConvertMapToRGB(data::Ptr{Cuchar}, count::Cint, depth::Cint, packed::Cint, palette::Ptr{Clong}, palette_count::Cint)
  ccall( (:imConvertMapToRGB, libim_convert_), Cvoid, (Ptr{Cuchar}, Cint, Cint, Cint, Ptr{Clong}, Cint), data, count, depth, packed, palette, palette_count)
end

function imConvertRGB2Map(width::Cint, height::Cint, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, map::Ptr{Cuchar}, palette::Ptr{Clong}, palette_count::Ptr{Cint})
  ccall( (:imConvertRGB2Map, libim_convert_), Cint, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Clong}, Ptr{Cint}), width, height, red, green, blue, map, palette, palette_count)
end


text:
"""
    GumbelCopula{d,T}

Fields:
  - θ::Real - parameter

Constructor

    GumbelCopula(d, θ)

The [Gumbel](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Most_important_Archimedean_copulas) copula in dimension ``d`` is parameterized by ``\\theta \\in [1,\\infty)``. It is an Archimedean copula with generator :

```math
\\phi(t) = \\exp{-t^{\\frac{1}{θ}}}
```
"""
struct GumbelCopula{d,T} <: ArchimedeanCopula{d}
    θ::T
end
GumbelCopula(d,θ) = θ >= 1 ? GumbelCopula{d,typeof(θ)}(θ) : @error "Theta must be greater than 1."
ϕ(  C::GumbelCopula,       t) = exp(-t^(1/C.θ))
ϕ⁻¹(C::GumbelCopula,       t) = (-log(t))^C.θ
τ(C::GumbelCopula) = (C.θ-1)/C.θ
τ⁻¹(::Type{GumbelCopula},τ) =1/(1-τ)

function radial_dist(C::GumbelCopula)
    α = 1/C.θ
    β = 1
    γ = cos(π/(2C.θ))^C.θ
    δ = C.θ == 1 ? 1 : 0
    AlphaStable([α,β,γ,δ]...) # for the type promotion...
end

# S(α, β, γ , δ) denotes a stable distribution in
# 1-parametrization [16, p. 8] with characteristic exponent α ∈ (0, 2], skewness β ∈ [−1, 1], scale
# γ ∈ [0,∞), and location δ ∈ R


text:
module Genser

import Base: UUID, ==
using Requires

export
    Tag, GenserDataType, GenserValue, GenserNothingValue, GenserUndefined, GenserNull,
    GenserInt8, GenserUInt8, GenserInt16, GenserUInt16, GenserInt32, GenserUInt32,
    GenserInt64, GenserUInt64, GenserInt128, GenserUInt128,
    GenserFloat16, GenserFloat32, GenserFloat64,
    GenserRational, GenserBigFloat,
    GenserBool, GenserChar, GenserBigInt,
    GenserNumberValue, GenserRealValue,
    GenserIntegerValue, GenserSignedValue, GenserUnsignedValue, GenserFloatValue,
    GenserStringValue, GenserString, GenserURI,
    GenserBinaryValue, GenserBinary, GenserUUID, GenserSymbol,
    GenserSequence, GenserSet, GenserTuple, GenserDict, GenserRecord,
    GenserOptional, GenserVariant, GenserAny,
    gensertypefor, togenser, fromgenser, convert_to_type,
    serialize, deserialize,
    genser_type_for_mime, genser_converter_for_encoding,
    TypeID, Encoding,
    tag, encoding,
    @fieldencoding, @fieldtype

    include("encoding.jl")
    include("datamodel.jl")
    include("derivetype.jl")
    include("convert.jl")
    include("serialize.jl")
    include("deserialize.jl")
    include("registry.jl")
    include("macros.jl")

    function __init__()
        @require JSON="682c06a0-de6a-54ab-a142-c8b1cf79cde6" include("types/json.jl")
    end
end


text:
using PProf
function test()
    for _ in 1:10
        println("Hello World")
    end
end

@pprof test()

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mcke

max_forks_repo_name: permissionx/mcke


text:
include("Backend.jl")

export
	DeviceInfo,
	DType,
	Dim4,
	DimT,
	f32,
	c32,
	f64,
	c64,
	b8,
	s32,
	u32,
	u8,
	s64,
	u64

immutable DeviceInfo
	id
	name
	platform
	toolkit
	compute
end

function toString(info::DeviceInfo)
	"ID: $(info.id)\nName: $(info.name)\nPlatform: $(info.platform)\nToolkit: $(info.toolkit)\nCompute: $(info.compute)"
end

Base.print(io::IOBuffer, info::DeviceInfo) = Base.print(io, toString(info))
Base.print(info::DeviceInfo) = Base.print(toString(info))

typealias DType UInt32

const f32 = DType(0) # 32-bit floating point values
const c32 = DType(1) # 32-bit complex floating point values
const f64 = DType(2) # 64-bit complex floating point values
const c64 = DType(3) # 64-bit complex floating point values
const b8  = DType(4) #  8-bit boolean values
const s32 = DType(5) # 32-bit signed integral values
const u32 = DType(6) # 32-bit unsigned integral values
const u8  = DType(7) #  8-bit unsigned integral values
const s64 = DType(8) # 64-bit signed integral values
const u64 = DType(9) # 64-bit unsigned integral values

const AF_MAT_NONE 		= Int32(0)	# Default.
const AF_MAT_TRANS 		= Int32(1)	# Data needs to be transposed.
const AF_MAT_CTRANS 	= Int32(2)	# Data needs to be conjugate tansposed.
const AF_MAT_CONJ 		= Int32(3)	# Data needs to be conjugate.
const AF_MAT_UPPER 		= Int32(4)	# Matrix is upper triangular.
const AF_MAT_LOWER 		= Int32(5)	# Matrix is lower triangular.
const AF_MAT_DIAG_UNIT 	= Int32(6)	# Matrix diagonal contains unitary values.
const AF_MAT_SYM 		= Int32(7)	# Matrix is symmetric.
const AF_MAT_POSDEF		= Int32(8)	# Matrix is positive definite.
const AF_MAT_ORTHOG 	= Int32(9)	# Matrix is orthogonal.
const AF_MAT_TRI_DIAG 	= Int32(10)	# Matrix is tri diagonal.
const AF_MAT_BLOCK_DIAG = Int32(11)	# Matrix is block diagonal.

asDType(::Type{Float32})          = f32
asDType(::Type{Complex{Float32}}) = c32
asDType(::Type{Float64})          = f64
asDType(::Type{Complex{Float64}}) = c64
asDType(::Type{Bool})             = b8
asDType(::Type{Int32})            = s32
asDType(::Type{UInt32})           = u32
asDType(::Type{UInt8})            = u8
asDType(::Type{Int64})            = s64
asDType(::Type{UInt64})           = u64

asJType(::Type{Val{f32}}) = Float32
asJType(::Type{Val{c32}}) = Complex{Float32}
asJType(::Type{Val{f64}}) = Float64
asJType(::Type{Val{c64}}) = Complex{Float64}
asJType(::Type{Val{b8}}) = Bool
asJType(::Type{Val{s32}}) = Int32
asJType(::Type{Val{u32}}) = UInt32
asJType(::Type{Val{u8}}) = UInt8
asJType(::Type{Val{s64}}) = Int64
asJType(::Type{Val{u64}}) = UInt64

typealias DimT Int64
typealias Dim4 Vector{DimT}

dimsToDim4(dims) =
	if length(dims) == 1
        [dims[1]]
    elseif length(dims) == 2
        [dims[1], dims[2]]
    elseif length(dims) == 3
        [dims[1], dims[2], dims[4]]
    elseif length(dims) == 4
        [dims[1], dims[2], dims[3], dims[4]]
    else
        throw(ArgumentError("Too many dimensions"))
    end

macro sizeRules(len, d)
	quote
		if $len > 3 && $d[4] > 1
			($d[1], $d[2], $d[3], $d[4])
		elseif $len > 2 && $d[3] > 1
			($d[1], $d[2], $d[3])
		elseif $len > 1 && $d[2] > 1
			($d[1], $d[2])
		elseif $len > 0 && $d[1] > 0
			($d[1],)
		else
			()
		end
	end
end

function dimsToSize(d::DimT...)
	len = length(d)
	@sizeRules(len, d)
end

function dimsToSize(d::Dim4)
	len = length(d)
	@sizeRules(len, d)
end

function dimsToSize(d::Tuple)
	len = length(d)
	@sizeRules(len, d)
end

jlLength(dims::Tuple) = length(dimsToSize(dims))

jlLength(dims::DimT...) = length(dimsToSize(dims...))

jlLength(dims::Dim4) = length(dimsToSize(dims))

function afPromote{T,S}(::Type{T},::Type{S})
    if T == S
        return T
    elseif T == Complex{Float64} || S == Complex{Float64}
        return Complex{Float64}
    elseif T == Complex{Float32} || S == Complex{Float32}
        (T == Float64 || S == Float64) && return Complex{Float64}
        return Complex{Float32}
    elseif T == Float64 || S == Float64
        return Float64
    elseif T == Float32 || S == Float32
        return Float32
    elseif T == UInt64 || S == UInt64
        return UInt64
    elseif T == Int64 || S == Int64
        return Int64
    elseif T == UInt32 || S == UInt32
        return UInt32
    elseif T == Int32 || S == Int32
        return Int32
    elseif T == UInt8 || S == UInt8
        return UInt8
    elseif T == Bool || S == Bool
        return Bool
    else
        return Float32
    end
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using wgregseq, DataFrames, CSV, BioSequences, CairoMakie, Statistics, StatsBase

wgregseq.plotting_style.default_makie!()

# Set path
dir = @__DIR__
home_dir = joinpath(split(dir, "/")[1:end-2])

##
filename = "2022-02-15_twist_order.csv"

df = CSV.read(
    "/$home_dir/data/twist_orders/$filename",
    types=Dict(
        "fwd_primer" => String,
        "rev_primer1" => String,
        "rev_primer2" => String,
        "rev_primer3" => String,
    ),
    DataFrame
)
df.fwd_primer = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.fwd_primer)
df.rev_primer1 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer1)
df.rev_primer2 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer2)
df.rev_primer3 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer3)

df.genes = parse.(Vector{String}, df.genes)
df.sequence = [LongDNA{4}(seq) for seq in df.sequence]
##
#wgregseq.quality_control.check_dataframe(df, site_start=27, site_end=186, print_results=false)
wgregseq.quality_control.check_primers_off_target(df)

##
gdf = groupby(df, :rev_primer3)
for _df in gdf
    println([pr[1] for pr in _df.rev_primer3 |> unique])
    println([pr[1] for pr in _df.rev_primer2 |> unique])
    println(nrow(_df))
    println(_df.promoter |> unique |> length)
    println()
end


text:
using Random
using TemporalGPs: RectilinearGrid, SpaceTimeGrid

function FiniteDifferences.to_vec(x::RectilinearGrid)
    v, tup_from_vec = to_vec((x.xl, x.xr))
    function RectilinearGrid_from_vec(v)
        tup = tup_from_vec(v)
        return RectilinearGrid(tup[1], tup[2])
    end
    return v, RectilinearGrid_from_vec
end

@testset "rectilinear_grid" begin
    rng = MersenneTwister(123456)
    Nl = 5
    Nr = 3
    xl = randn(rng, Nl)
    xr = randn(rng, Nr)

    X = RectilinearGrid(xl, xr)
    x = collect(X)
    @test x isa Vector
    @test eltype(x) == eltype(X)
    @test size(X) == (length(x),)
    @test length(x) == length(X)
end


text:
using Base.Test
using MPI

MPI.Init()

function scatter_array(A, root)
    comm = MPI.COMM_WORLD
    T = eltype(A)
    if MPI.Comm_rank(comm) == root
        B = copy(A)
    else
        B = Array(T,1)
    end
    C = MPI.Scatter(B, 1, root, comm)
end

comm = MPI.COMM_WORLD
root = 0

A = collect(1:MPI.Comm_size(comm))
B = scatter_array(A, root)
@test B[1] == MPI.Comm_rank(comm) + 1
for typ in MPI.MPIDatatype.types
    A = convert(Vector{typ},collect(1:MPI.Comm_size(comm)))
    B = scatter_array(A, root)
    @test B[1] == convert(typ,MPI.Comm_rank(comm) + 1)
end

MPI.Finalize()


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using DocOpt

include("image.jl")
using .Image

# Solve optical flow problem at one scale
# Flow problem is formulated as the solution of a linear system where matrix is sparse
# Use an iterative method to
# Pre: size(i1)==size(i2)
function singleScaleOpticalFlow(i1::Matrix{Float32}, i2::Matrix{Float32}, lam::Float32, ni::Int)
# Returns Matrix{Float32} size=(size(i1,1),size(i1,2),2)
  (h, w) = size(i1)
  # Compute the gradient of i1
  # This is basically two stencils, one for x and one for y
  # The stencil is [1 8 0 -8 -1]/12 and leave zeros in the boundary
  Ix = zeros(Float32, h, w)
  Iy = zeros(Float32, h, w)
  for y in 1:h
    for x in 1:w
      if x>2 && x<w-1
        Ix[y,x] = (-i1[y,x+2]-8.0f0*i1[y,x+1]+8.0f0*i1[y,x-1]+i1[y,x-2])/12.0f0
      end
      if y>2 && y<h-1
        Iy[y,x] = (-i1[y+2,x]-8.0f0*i1[y+1,x]+8.0f0*i1[y-1,x]+i1[y-2,x])/12.0f0
      end
    end
  end
  # Compute the time partial derivative
  It = i1-i2

  # Block Jacobi preconditioner
  function blockJacobiPreconditioner(r)
    a = Ix.*Ix.+4.0f0*lam
    b = Ix.*Iy
    c = Iy.*Iy.+4.0f0*lam
    invdet = 1.0f0./(a.*c-b.*b)
    res = zeros(Float32, h, w, 2)
    res[:,:,1] = (c.*r[:,:,1]-b.*r[:,:,2]).*invdet
    res[:,:,2] = (a.*r[:,:,2]-b.*r[:,:,1]).*invdet
    return res
  end

  # Matrix vector product
  function mvm(x)
    # This operation is equivalent to a stencil
    # Want:
    #   -lam for (-1,0,0), (+1,0,0), (0,-1,0), and (0,+1,0)
    #   Ix*Iy for (0,0,flip)
    #   Ix^2+4*lam and Iy^2+4*lam for (0,0,0) respectively (Ix for u, Iy for v)
    # At boundaries just include terms that exist
    y = zeros(Float32, h, w, 2)
    for i in 1:h
      for j in 1:w
        ur = Ix[i,j]*(Ix[i,j]*x[i,j,1]+Iy[i,j]*x[i,j,2])+4.0f0*lam*x[i,j,1]
        vr = Iy[i,j]*(Ix[i,j]*x[i,j,1]+Iy[i,j]*x[i,j,2])+4.0f0*lam*x[i,j,2]
        if i>1
          ur += -lam*x[i-1,j,1]
          vr += -lam*x[i-1,j,2]
        end
        if j>1
          ur += -lam*x[i,j-1,1]
          vr += -lam*x[i,j-1,2]
        end
        if i<h
          ur += -lam*x[i+1,j,1]
          vr += -lam*x[i+1,j,2]
        end
        if j<w
          ur += -lam*x[i,j+1,1]
          vr += -lam*x[i,j+1,2]
        end
        y[i,j,1] = ur
        y[i,j,2] = vr
      end
    end
    return y
  end

  # Initial for iterative solve
  x = zeros(Float32, h, w, 2)
  r = zeros(Float32, h, w, 2)
  r[:,:,1] = -Ix.*It
  r[:,:,2] = -Iy.*It
  z = blockJacobiPreconditioner(r)
  p = z
  rsold = sum(r.*z)

  # Iterate to solution
  # Use preconditioned conjugate gradient algorithm
  for i in 1:ni
    Ap = mvm(p)
    pTAp = sum(p.*Ap)
    alpha = rsold/pTAp
    x += alpha*p
    r -= alpha*Ap
    z = blockJacobiPreconditioner(r)
    rsnew = sum(r.*z)
    beta = rsnew/rsold
    p = z + beta*p
    rsold = rsnew
  end

  return x
end

# Solve optical flow problem at all scales
# Pre: size(i1)==size(i2)
function multiScaleOpticalFlow(i1::Matrix{Float32}, i2::Matrix{Float32}, lam::Float32, ni::Int, ns::Int)
# Returns Matrix{Float32} size=(size(i1,1),size(i1,2),2)
  (h, w) = size(i1)
  if ns==1
    return singleScaleOpticalFlow(i1, i2, lam, ni)
  end
  scale = Float32((Float64(w)/50.0)^(-1.0/ns))
  local u # initialised in the first iteration of the loop
  for i in ns:-1:0
    cw = round(Int, w*scale^Float32(i))
    ch = round(Int, h*scale^Float32(i))
    println("Doing scale: ", cw, "x", ch)
    # Downsample the images
    si1 = Image.downSample(i1, cw, ch)
    si2 = Image.downSample(i2, cw, ch)
    # Interpolate the current flow and warp image 2 by current flow
    if i==ns
      # Coarsest scale, no flow yet, use zeros, warp is image 2
      u = zeros(Float32, ch, cw, 2)
      wi = si2
    else
      u = Image.interpolateFlow(u, cw, ch)
      wi = Image.warpMotion(si2, u, si1)
    end
    # Compute the flow from image 1 to warp
    if i==0
      ni *= 5
    end
    du = singleScaleOpticalFlow(si1, wi, lam, ni)
    # Update the flow
    u += du
  end
  return u
end

function main()
    doc = """opt-flow.jl

Horn-Schunck optical flow estimator.

Usage:
  opt-flow.jl -h | --help
  opt-flow.jl [--image-path-1=<image-path-1>] [--image-path-2=<image-path-2>]

Options:
  -h --help                      Show this screen.
  --image-path-1=<image-path-1>  Specify a path to the first image [default: ../../opt-flow/small_001.dat].
  --image-path-2=<image-path-2>  Specify a path to the second image [default: ../../opt-flow/small_002.dat].
"""
    arguments = docopt(doc)

    dir = dirname(@__FILE__)
    img_file_1 = joinpath(dir, arguments["--image-path-1"])
    img_file_2 = joinpath(dir, arguments["--image-path-2"])

    i1 = Image.readImage(img_file_1)
    i2 = Image.readImage(img_file_2)
    if size(i1)!=size(i2)
        write(STDERR, "images are different sizes")
        exit(-1)
    end

    println("filenames = ", img_file_1, " ", img_file_2)
    println("Image size: ", size(i1, 1), "x", size(i1, 2))

    tic()
    multiScaleOpticalFlow(i1, i2, Float32(0.025), 100, 1)
    println("SELFPRIMED ", toq())

    tic()
    u = multiScaleOpticalFlow(i1, i2, Float32(0.025), 100, 44)
    selftimed = toq()
    Image.writeFlo(u, "out.flo")
    println("SELFTIMED ", selftimed)
end

main()


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

## Initialization

using PyCall, NLopt

@pyimport scipy.io as si

include("displayData.jl")
include("nnCostFunction.jl")
include("randInitializeWeights.jl")
include("checkNNGradients.jl")
include("predict.jl")
# include("fmincg.jl")

## Setup the parameters you will use for this exercise
input_layer_size  = 400  # 20x20 Input Images of Digits
hidden_layer_size = 25   # 25 hidden units
num_labels = 10          # 10 labels, from 1 to 10
                         # (note that we have mapped "0" to label 10)

## =========== Loading and Visualizing Data =============

# Load Training Data
@printf("Loading and Visualizing Data ...\n")

data = si.loadmat("ex4data1.mat")
X = data["X"]
y = data["y"]
m = size(X, 1)

# Randomly select 100 data points to display
sel = randperm(size(X, 1))
sel = sel[1:100]

displayData(X[sel, :])

@printf("Program paused. Press enter to continue.\n")
readline()

## ================  Loading Parameters ================

@printf("\nLoading Saved Neural Network Parameters ...\n")

# Load the weights into variables Theta1 and Theta2
nn = si.loadmat("ex4weights.mat")
Theta1 = nn["Theta1"]
Theta2 = nn["Theta2"]

# Unroll parameters
nn_params = [Theta1[:] ; Theta2[:]]

## ================ Compute Cost (Feedforward) ================
\@printf("\nFeedforward Using Neural Network ...\n")

# Weight regularization parameter (we set this to 0 here).
lambda = 0

J, _ = nnCostFunction(nn_params, input_layer_size, hidden_layer_size,
                      num_labels, X, y, lambda)

println(J)
@printf("Cost at parameters (loaded from ex4weights): %f\n(this value should be about 0.287629)\n", J)

@printf("\nProgram paused. Press enter to continue.\n")
readline()

## =============== Implement Regularization ===============

@printf("\nChecking Cost Function (w/ Regularization) ... \n")

# Weight regularization parameter (we set this to 1 here).
lambda = 1

J, _ = nnCostFunction(nn_params, input_layer_size, hidden_layer_size,
                      num_labels, X, y, lambda)

@printf("Cost at parameters (loaded from ex4weights): %f\n(this value should be about 0.383770)\n", J)

@printf("Program paused. Press enter to continue.\n")
readline()

## ================ Sigmoid Gradient  ================

@printf("\nEvaluating sigmoid gradient...\n")

g = sigmoidGradient([1 -0.5 0 0.5 1])
@printf("Sigmoid gradient evaluated at [1 -0.5 0 0.5 1]:\n  ")
@printf("%s", join(map(x -> @sprintf("%5f", x), g), " "))
@printf("\n\n")

@printf("Program paused. Press enter to continue.\n")
readline()

## ================ Initializing Parameters ================

@printf("\nInitializing Neural Network Parameters ...\n")

initial_Theta1 = randInitializeWeights(input_layer_size, hidden_layer_size)
initial_Theta2 = randInitializeWeights(hidden_layer_size, num_labels)

# Unroll parameters
initial_nn_params = [initial_Theta1[:] ; initial_Theta2[:]]

## =============== Implement Backpropagation ===============

@printf("\nChecking Backpropagation... \n")

#  Check gradients by running checkNNGradients
checkNNGradients()

@printf("\nProgram paused. Press enter to continue.\n")
readline()

## =============== Implement Regularization ===============

@printf("\nChecking Backpropagation (w/ Regularization) ... \n")

#  Check gradients by running checkNNGradients
lambda = 3
checkNNGradients(lambda)

# Also output the costFunction debugging values
debug_J, _  = nnCostFunction(nn_params, input_layer_size,
                             hidden_layer_size, num_labels, X, y, lambda)

@printf("""

  Cost at (fixed) debugging parameters (w/ lambda = 10): %f
  (this value should be about 0.576051)

""", debug_J)

@printf("Program paused. Press enter to continue.\n")

## =================== Training NN ===================
@printf("\nTraining Neural Network... \n")

# options = optimset('MaxIter', 50)
options = Opt(:LD_TNEWTON, length(initial_nn_params))
maxeval!(options, 50)

#  You should also try different values of lambda
lambda = 1

# Create "short hand" for the cost function to be minimized
function costFunction(theta, grad)
  cost, gs = nnCostFunction(theta, input_layer_size, hidden_layer_size,
                            num_labels, X, y, lambda)
  if length(grad) > 0
    grad[:] = gs
  end
  return cost
end

# Now, costFunction is a function that takes in only one argument (the
# neural network parameters)
# [nn_params, cost] = fmincg(costFunction, initial_nn_params, options);

min_objective!(options, costFunction)
(cost, theta, _) = optimize(options, initial_nn_params)

# Obtain Theta1 and Theta2 back from nn_params
Theta1 = reshape(nn_params[1:hidden_layer_size * (input_layer_size + 1)],
                 hidden_layer_size, (input_layer_size + 1))

Theta2 = reshape(nn_params[(1 + (hidden_layer_size * (input_layer_size + 1))):end],
                 num_labels, (hidden_layer_size + 1))

@printf("Program paused. Press enter to continue.\n")
readline()

## =================  Visualize Weights =================

@printf("\nVisualizing Neural Network... \n")

displayData(Theta1[:, 2:end])

@printf("\nProgram paused. Press enter to continue.\n")
readline()

## ================= Implement Predict =================

pred = predict(Theta1, Theta2, X)

@printf("\nTraining Set Accuracy: %f\n", mean(pred .== y) * 100)


text:
# We document types/constructors/functors in one docstring until
# https://github.com/JuliaDocs/Documenter.jl/issues/558 is resolved
"""
    AbstractAggregation

Supertype for any aggregation operator.

See also: [`AggregationStack`](@ref), [`SegmentedSum`](@ref), [`SegmentedMax`](@ref),
    [`SegmentedMean`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
"""
abstract type AbstractAggregation end

@inline _bagnorm(w::Nothing, b) = length(b)
@inline _bagnorm(w::AbstractVector, b) = @views sum(w[b])
@inline _bagnorm(w::AbstractMatrix, b) = @views vec(sum(w[:, b], dims=2))

@inline _weight(w::Nothing, _, _, ::Type{T}) where T = one(T)
@inline _weight(w::AbstractVector, _, j, _) = w[j]
@inline _weight(w::AbstractMatrix, i, j, _) = w[i, j]

@inline _weightsum(ws::Real, _) = ws
@inline _weightsum(ws::AbstractVector, i) = ws[i]

# more stable definitions for r_map and p_map
ChainRulesCore.rrule(::typeof(softplus), x) = softplus.(x), Δ -> (NoTangent(), Δ .* σ.(x))

# our definition of type min for Maybe{...} types
_typemin(t::Type) = typemin(t)
_typemin(::Type{Missing}) = missing
_typemin(::Type{Maybe{T}}) where T = typemin(T)

function _check_agg(a::AbstractAggregation, X::Missing) end
function _check_agg(a::AbstractAggregation, X::AbstractMatrix)
    if size(X, 1) != length(a.ψ)
        DimensionMismatch(
            "Different number of rows in input ($(size(X, 2))) and aggregation ($(length(a.ψ)))"
        ) |> throw
    end
end

include("segmented_sum.jl")
include("segmented_mean.jl")
include("segmented_max.jl")
include("segmented_pnorm.jl")
include("segmented_lse.jl")
# include("transformer.jl")

include("aggregation_stack.jl")

Base.vcat(as::AbstractAggregation...) = reduce(vcat, as |> collect)
function Base.reduce(::typeof(vcat), as::Vector{<:AbstractAggregation})
    AggregationStack(tuple(as...))
end

include("bagcount.jl")

(a::Union{AbstractAggregation, BagCount})(x::ArrayNode, args...) = ArrayNode(a(x.data, args...))

# definitions for mixed aggregations
const names = ["Sum", "Mean", "Max", "PNorm", "LSE"]
for p in filter(p -> length(p) > 1, collect(powerset(collect(1:length(names)))))
    s = Symbol("Segmented", names[p]...)
    @eval begin
        """
            $($(s))([t::Type, ]d::Int)

        Construct [`AggregationStack`](@ref) consisting of $($(
            join("[`Segmented" .* names[p] .* "`](@ref)", ", ", " and ")
       )) operator$($(length(p) > 1 ? "s" : "")).

        $($(
            all(in(["Sum", "Mean", "Max"]), names[p]) ? """
            # Examples
            ```jldoctest
            julia> $(s)(4)
            AggregationStack:
            $(join(" Segmented" .* names[p] .* "(ψ = Float32[0.0, 0.0, 0.0, 0.0])", "\n"))

            julia> $(s)(Float64, 2)
            AggregationStack:
            $(join(" Segmented" .* names[p] .* "(ψ = [0.0, 0.0])", "\n"))
            ```
            """ : ""
        ))

        See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`SegmentedSum`](@ref),
            [`SegmentedMax`](@ref), [`SegmentedMean`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
        """
        function $s(d::Int)
            AggregationStack($((Expr(:call, Symbol("Segmented", n), :d) for n in names[p])...))
        end
    end
    @eval function $s(::Type{T}, d::Int) where T
        AggregationStack($((Expr(:call, Symbol("Segmented", n), :T, :d) for n in names[p])...))
    end
    @eval export $s
end


text:
using GeneralizedCRT
using Test

@testset "GeneralizedCRT" begin

@testset "crt basic" begin
    a, b, p, q = 1, 11, 20, 50
    @test crt((a, p)) == (mod(a, p), p)
    @test crt([(a, p), (b, q)]) == crt(a, b, p, q)
    @test crt(a, b, p, q) == (61, 100)
    @test crt((a, p), (b, q)) == crt(a, b, p, q)
    a = 0
    @test_throws ArgumentError crt(a, b, p, q)
end

@testset "crt_short" begin
    a = Int[]
    m = Int[]
    @test crt(a, m) == (0, 1)

    a = [12]
    m = [10]
    @test crt(a, m) == (2, 10)

    a = [1, 11]
    m = [20, 50]
    @test crt(a, m) == crt(a..., m...)

    a = [1, 1, 1, 1, 1, 0]
    m = [2, 3, 4, 5, 6, 7]
    @test crt(a, m) == (301, 420)

    a = [5, 3, 1]
    m = [6, 8, 10]
    x, = crt(a, m)
    @test a == mod.(x, m)

    a = [2, 5, 7]
    m = [6, 9, 15]
    @test_throws ArgumentError crt(a, m)

    a = vec([2  0  3  3  17  0  0  0  3  1])
    m = vec([5  13  17  12  25  3  21  1  8  2])
    x, lc = crt(a, m)
    @test (x, lc) == (594867, 928200)
    @test a == mod.(x, m)
end

using Random
rng = MersenneTwister(1)

@testset "crt_long" begin
    n = max(GeneralizedCRT.THRESHOLD1, GeneralizedCRT.THRESHOLD2) * 1000 - 100
    m = rand(rng, 1:999, n)
    l = lcm(big.(m))
    x = rand(rng, 1:l)
    a = mod.(x, m)
    @test crt(a, m) == (x, l)
    @test_throws OverflowError crt(zeros(Int, n), m)
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

using FEMBase, AsterReader, Test
using AsterReader: RMEDFile, aster_read_data

datadir = first(splitext(basename(@__FILE__)))

@testset "read nodal field from result rmed file" begin
    rmedfile = joinpath(datadir, "rings.rmed")
    rmed = RMEDFile(rmedfile)
    temp = aster_read_data(rmed, "TEMP")
    @test isapprox(temp[15], 1.0)
    @test isapprox(temp[95], 2.0)
end


text:
function commonpath(ds::Vector{<:AbstractString}, dlm::Char='/')
    0 < length(ds) || return ""
    1 < length(ds) || return String(ds[1])
    p = split(ds[1], dlm)
    mincnt = length(p)
    for d in ds[2:end]
        q = split(d, dlm)
        mincnt = min(mincnt, length(q))
        hits = findfirst(p[1:mincnt] .!= q[1:mincnt])
        if hits != 0 mincnt = hits - 1 end
        if mincnt == 0 return "" end
    end
    1 < mincnt || p[1] != "" || return convert(T, string(dlm))
    return join(p[1:mincnt], dlm)
end

test = ["/home/user1/tmp/coverage/test", "/home/user1/tmp/covert/operator", "/home/user1/tmp/coven/members"]

println("Comparing:\n - ", join(test, "\n - "))
println("for their common directory path yields:\n", commonpath(test))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using HDF5

function test_hdf5_data_layer(backend::Backend, async, T, eps)
  println("-- Testing $(async ? "(Async)" : "") HDF5 Data Layer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  batch_size = 3
  tensor_dim = abs(rand(Int)) % 6 + 1
  data_dim = tuple(rand(1:8, tensor_dim)...)
  println("    > $data_dim")

  data_all = [rand(T, data_dim..., x) for x in [5 1 2]]
  h5fn_all = [string(Mocha.temp_filename(), ".hdf5") for x in 1:length(data_all)]

  for i = 1:length(data_all)
    h5 = h5open(h5fn_all[i], "w")
    h5["data"] = data_all[i]
    close(h5)
  end

  source_fn = string(Mocha.temp_filename(), ".txt")
  open(source_fn, "w") do s
    for fn in h5fn_all
      println(s, fn)
    end
  end

  ############################################################
  # Setup
  ############################################################

  scale = rand(T)
  if async
    layer = AsyncHDF5DataLayer(source = source_fn, tops = [:data], batch_size=batch_size, chunk_size=4,
        transformers=[(:data, DataTransformers.Scale(scale))])
  else
    layer = HDF5DataLayer(source = source_fn, tops = [:data], batch_size=batch_size,
        transformers=[(:data, DataTransformers.Scale(scale))])
  end
  state = setup(backend, layer, Blob[], Blob[])
  @test state.epoch == 0

  data = cat(tensor_dim+1, data_all...)
  data = cat(tensor_dim+1, data, data, data)
  data = data * scale

  data_idx = map(x->1:x, data_dim)
  layer_data = Array{eltype(data)}(tuple(data_dim..., batch_size))
  for i = 1:batch_size:size(data)[end]-batch_size+1
    forward(backend, state, Blob[])
    copy!(layer_data, state.blobs[1])
    @test all(-eps .< layer_data - data[data_idx..., i:i+batch_size-1] .< eps)
  end
  @test state.epoch == 3

  ############################################################
  # Clean up
  ############################################################
  shutdown(backend, state)
  rm(source_fn)
  for fn in h5fn_all
    # workaround for issue #58
    try
      rm(fn)
    catch e
      println(e)
    end
  end
end

function test_hdf5_data_layer(backend::Backend, T, eps)
  test_hdf5_data_layer(backend, false, T, eps)
  test_hdf5_data_layer(backend, true, T, eps)
end

function test_hdf5_data_layer_shuffle(backend::Backend, batch_size, async, n, T)
  println("-- Testing $(async ? "(Async)" : "") HDF5 Data Layer (shuffle,n=$n,b=$batch_size) on $(typeof(backend)){$T}...")

  # To test random shuffling, we generate a dataset containing integer 1:n.
  # Then we run HDF5 layer n times forward, and collect all the output data.
  # For each integer i in 1:n, it should show up exactly b times, where b
  # is the batch size.

  data = reshape(convert(Array{T}, collect(1:n)), 1, 1, 1, n)
  h5fn = string(Mocha.temp_filename(), ".hdf5")
  h5open(h5fn, "w") do file
    file["data"] = data
  end

  source_fn = string(Mocha.temp_filename(), ".txt")
  open(source_fn, "w") do file
    println(file, h5fn)
  end

  if async
    layer = AsyncHDF5DataLayer(source=source_fn, tops=[:data], batch_size=batch_size, shuffle=true, chunk_size=2)
  else
    layer = HDF5DataLayer(source=source_fn, tops=[:data], batch_size=batch_size, shuffle=true)
  end
  state = setup(backend, layer, Blob[], Blob[])

  data_got_all = Int[]
  data_got = zeros(T, 1, 1, 1, batch_size)
  for i = 1:n
    forward(backend, state, Blob[])
    copy!(data_got, state.blobs[1])
    append!(data_got_all, convert(Vector{Int}, data_got[:]))
  end

  for i = 1:n
    @test sum(data_got_all .== i) == batch_size
  end

  if data_got_all[1:n] == collect(1:n)
    println("WARNING: data not shuffled, is today a lucky day or is there a bug?")
  end

  shutdown(backend, state)
  rm(source_fn)
  try
    rm(h5fn)
  catch e
    println(e)
  end
end
function test_hdf5_data_layer_shuffle(backend::Backend, batch_size, n, T)
  # do not run (non-async) HDF5 data layer shuffling on windows, because it is implemented
  # with memmap, which is not working properly on Windows.
  @static is_windows() ? nothing : test_hdf5_data_layer_shuffle(backend, batch_size, false, n, T)

  test_hdf5_data_layer_shuffle(backend, batch_size, true, n, T)
end

function test_hdf5_data_layer_shuffle(backend::Backend, T)
  test_hdf5_data_layer_shuffle(backend, 4, 6, T)
  test_hdf5_data_layer_shuffle(backend, 6, 4, T)
end

function test_hdf5_data_layer(backend::Backend)
  test_hdf5_data_layer(backend, Float32, 1e-3)
  test_hdf5_data_layer(backend, Float64, 1e-7)
  test_hdf5_data_layer_shuffle(backend, Float32)
  test_hdf5_data_layer_shuffle(backend, Float64)
end

if test_cpu
  test_hdf5_data_layer(backend_cpu)
end
if test_gpu
  test_hdf5_data_layer(backend_gpu)
end


text:
##############################   Coulomb  ###################################
mutable struct Coulomb <: EmpiricalPotential
    q_1
    q_2
    ϵ0
    rcutoff
    species :: AbstractVector
end

get_parameters(c::Coulomb) = Parameter{}(())
set_parameters(p::Parameter{}, c::Coulomb) = c

deserialize_parameters(p::Parameter{()}, c::Coulomb) = []
serialize_parameters(p::Vector, c::Coulomb) = Parameter{()}( () )

get_hyperparameters(c::Coulomb) = Parameter{:rcutoff}((c.rcutoff, ))
set_hyperparameters(p::Parameter{(:rcutoff,)}, c::Coulomb) = Coulomb(c.q_1, c.q_2, c.ϵ0, p.rcutoff, c.species)

deserialize_hyperparameters(p::Parameter{(:rcutoff,)}, c::Coulomb) = [p.rcutoff]
serialize_hyperparameters(p::Vector, c::Coulomb) = Parameter{(:rcutoff, )}( (p[1],) )

# ############################# Energies ##########################################

function potential_energy(R::AbstractFloat, c::Coulomb)
    return c.q_1 * c.q_2 / (4.0 * π * c.ϵ0 * R)
end

# ############################### Forces ##########################################

function force(R::AbstractFloat, r::SVector{3,<:AbstractFloat}, c::Coulomb)
    SVector(c.q_1 * c.q_2 / (4.0 * π * c.ϵ0 * R^2) .* r ./ R)
end

# ##############################   Gradients  ###################################

function grad_potential_energy(r::Vector{<:Real}, p::Coulomb)
    error("The Coulomb potential has no trainable parameters")
end
function grad_force(r:: Vector{<:Real}, p::Coulomb)
    error("The Coulomb potential has no trainable parameters")
end

function grad_virial(r::Vector{<:Real}, p::Coulomb)
    error("The Coulomb potential has no trainable parameters")
end


text:
"""
    piecewise(X::Trajectory, [endtime]) -> tt, xx

If X is a jump process with piecewise constant paths and jumps in `X.tt`,
piecewise returns coordinates path for plotting purposes. The second argument
allows to choose the right endtime of the last interval.
"""
function piecewise(X::Trajectory, tend = X.t[end])
    t = [X.t[1]]
    n = length(X)
    append!(t, repeat(X.t[2:n], inner=2))
    push!(t, tend)
    t, repeat(X.x, inner=2)
end


text:
module TestData

const test_cases_data = (
    ([1,0,0,0,99],                    [2,0,0,0,99]),
    ([2,3,0,3,99],                    [2,3,0,6,99]),
    ([2,4,4,5,99,0],                  [2,4,4,5,99,9801]),
    ([1,1,1,4,99,5,6,0,99],           [30,1,1,4,2,5,6,0,99]),
    ([1,9,10,3,2,3,11,0,99,30,40,50], [3500,9,10,70,2,3,11,0,99,30,40,50]),
)

const test_cmp_data = (
    [3,9,8,9,10,9,4,9,99,-1,8],
    [3,9,7,9,10,9,4,9,99,-1,8],
    [3,3,1108,-1,8,3,4,3,99],
    [3,3,1107,-1,8,3,4,3,99],
)

const test_jmp_data = (
    [3,12,6,12,15,1,13,14,13,4,13,99,-1,0,1,9],
    [3,3,1105,-1,9,1101,0,0,12,4,12,99,1],
)

const test_larger_data = (
    [3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,
    1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,
    999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99],
)

end


text:
# test @tr_str
begin
    # for Cladograms
    parsedtree = @tr_str  "(A,B,(C,D)E)F;"
    @test parsedtree == Clado("",[5 1
                                  5 2
                                  5 6
                                  6 3
                                  6 4],2,
                              ["A","B","C","D"],
                              ["F","E"])

    # for Phylograms
    parsedtree = @tr_str "(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;"
    @test parsedtree == Phylo("", [5 1
                                   5 2
                                   5 6
                                   6 3
                                   6 4],2,
                              ["A","B","C","D"],
                              [0.1,0.2,0.5,0.3,0.4],
                              ["F","E"],
                              -1.0)
end

# test writenewick
begin
    # for Phylograms
    @test Phylogenetics.newick(Clado("test",[5 1
                                             5 2
                                             5 6
                                             6 3
                                             6 4],2,
                                     ["A","B","C","D"],
                                     ["F","E"]),false) == "(A,B,(C,D)E)F;"

    # for Cladograms
    @test Phylogenetics.newick(Phylo("test", [5 1
                                              5 2
                                              5 6
                                              6 3
                                              6 4],2,
                                     ["A","B","C","D"],
                                     [0.1,0.2,0.5,0.3,0.4],
                                     ["F","E"],
                                     -1.0),true) == "test(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;"
end


text:
############################################################################################
# Differentiation - Benchmark Model
modelExample = ModelWrapper(ExampleModel(), _val_examplemodel)
objectiveExample = Objective(modelExample, (data1, data2, data3, _idx))

@testset "AutoDiffContainer - Log Objective Results" begin
    tune_fwd = AutomaticDiffTune(:ForwardDiff, objectiveExample)
    fwd = DiffObjective(objectiveExample, tune_fwd)
    theta_unconstrained = randn(_RNG, length(modelExample))
    ## Compute logdensity
    log_density(objectiveExample)
    theta_unconstrained2 = deepcopy(theta_unconstrained)
    #!NOTE: 10th parameter in likelihood for example, so is not compiled away in Reverse Tape
    theta_unconstrained2[10] = Inf
    ## Check if result is finite
    _ld_fin = log_density(objectiveExample, copy(theta_unconstrained))
    _ld_inf = log_density(objectiveExample, copy(theta_unconstrained2))
    _ld_inf_fault = log_density(objectiveExample, copy(theta_unconstrained))
    _ld_inf_fault.θᵤ[10] = Inf

    @test ModelWrappers.checkfinite(_ld_fin)
    @test !ModelWrappers.checkfinite(_ld_inf)
    @test !ModelWrappers.checkfinite(_ld_inf_fault)

    @test ModelWrappers.checkfinite(_ld_inf, _ld_fin) #Infinite to Finite
    @test !ModelWrappers.checkfinite(_ld_fin, _ld_inf) #Finite to Infinite

    @test ModelWrappers.checkfinite(-Inf, 10.0, _ld_fin) #Infinite to Finite
    @test !ModelWrappers.checkfinite(10.0, -Inf, _ld_fin) #Finite to Infinite
    @test !ModelWrappers.checkfinite(-Inf, 10.0, _ld_inf) #Finite to Infinite
end

@testset "AutoDiffContainer - Log Objective AutoDiff compatibility - Vectorized Model" begin
    ## Assign DiffTune
    tune_fwd = AutomaticDiffTune(:ForwardDiff, objectiveExample)
    tune_rd = AutomaticDiffTune(:ReverseDiff, objectiveExample)
    tune_zyg = AutomaticDiffTune(:Zygote, objectiveExample)
    fwd = DiffObjective(objectiveExample, tune_fwd)
    rd = DiffObjective(objectiveExample, tune_rd)
    zyg = DiffObjective(objectiveExample, tune_zyg)
    theta_unconstrained = randn(_RNG, length(modelExample))
    ## Compute logdensity
    log_density(objectiveExample)
    theta_unconstrained2 = deepcopy(theta_unconstrained)
    #!NOTE: 10th parameter in likelihood for example, so is not compiled away in Reverse Tape
    theta_unconstrained2[10] = Inf
    _ld = log_density(objectiveExample, theta_unconstrained2)
    @test isinf(_ld.ℓθᵤ)
    _ld = log_density(fwd, theta_unconstrained2)
    @test isinf(_ld.ℓθᵤ)
    _ld = log_density(rd, theta_unconstrained2)
    @test isinf(_ld.ℓθᵤ)
    _ld = log_density(zyg, theta_unconstrained2)
    @test isinf(_ld.ℓθᵤ)

    ld1 = log_density(fwd, theta_unconstrained)
    ld2 = log_density(rd, theta_unconstrained)
    ld3 = log_density(zyg, theta_unconstrained)
    _ld1 = _log_density(objectiveExample, tune_fwd, theta_unconstrained)
    _ld2 = _log_density(objectiveExample, tune_rd, theta_unconstrained)
    _ld3 = _log_density(objectiveExample, tune_zyg, theta_unconstrained)
    ## Compute Diffresult
    log_density_and_gradient(objectiveExample, tune_fwd, theta_unconstrained2)
    _grad = log_density_and_gradient(fwd, theta_unconstrained2)
    @test isinf(_grad.ℓθᵤ)
    log_density_and_gradient(objectiveExample, tune_rd, theta_unconstrained2)
    _grad = log_density_and_gradient(rd, theta_unconstrained2)
    #!TODO: Need an example where tape does not compile parameter in infi
    @test isinf(_grad.ℓθᵤ)
    log_density_and_gradient(objectiveExample, tune_zyg, theta_unconstrained2)
    _grad = log_density_and_gradient(zyg, theta_unconstrained2)
    @test isinf(_grad.ℓθᵤ)
    _grad1 = _log_density_and_gradient(objectiveExample, tune_fwd, theta_unconstrained)
    _grad2 = _log_density_and_gradient(objectiveExample, tune_rd, theta_unconstrained)
    _grad3 = _log_density_and_gradient(objectiveExample, tune_zyg, theta_unconstrained)
    grad1 = log_density_and_gradient(fwd, theta_unconstrained)
    grad2 = log_density_and_gradient(rd, theta_unconstrained)
    grad3 = log_density_and_gradient(zyg, theta_unconstrained)
    ## Compute manual call ~ Already checked for equality
    ld = objectiveExample(theta_unconstrained)
    grad_mod_fd = ForwardDiff.gradient(objectiveExample, theta_unconstrained)
    grad_mod_rd = ReverseDiff.gradient(objectiveExample, theta_unconstrained)
    grad_mod_zy = Zygote.gradient(objectiveExample, theta_unconstrained)[1]
    ## Compare results
    @test ld - ld1.ℓθᵤ ≈ 0 atol = _TOL
    @test ld - ld2.ℓθᵤ ≈ 0 atol = _TOL
    @test ld - ld3.ℓθᵤ ≈ 0 atol = _TOL
    @test sum(abs.(_grad1[2] - grad1.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(_grad2[2] - grad2.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(_grad3[2] - grad3.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_fd - grad1.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_rd - grad2.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_zy - grad3.∇ℓθᵤ)) ≈ 0 atol = _TOL
    ## Checks
    _output = check_gradients(_RNG, objectiveExample, [:ForwardDiff, :ReverseDiff, :Zygote]; printoutput = false)
    @test sum(abs.(_output.ℓobjective_gradient_diff)) ≈ 0 atol = _TOL
    ## Update DiffTune
    ModelWrappers.update(tune_fwd, objectiveExample)
    ModelWrappers.update(tune_rd, objectiveExample)
    ModelWrappers.update(tune_zyg, objectiveExample)
    ## Config DiffTune
    theta_unconstrained2 = randn(_RNG, length(objectiveExample))
    ModelWrappers._config(ModelWrappers.ADForward(), objectiveExample, theta_unconstrained2)
    ModelWrappers._config(ModelWrappers.ADReverse(), objectiveExample, theta_unconstrained2)
    ModelWrappers._config(ModelWrappers.ADReverseUntaped(), objectiveExample, theta_unconstrained2)
    ModelWrappers._config(ModelWrappers.ADZygote(), objectiveExample, theta_unconstrained2)
end

############################################################################################
# Differentiation - Lower dimensions
modelLowerDim = ModelWrapper(LowerDims(), _val_lowerdims)
objectiveLowerDim = Objective(modelLowerDim, nothing)

@testset "AutoDiffContainer - Log Objective AutoDiff compatibility - Lower dimensions" begin
    ## Assign DiffTune
    autodiff_fd = AutomaticDiffTune(:ForwardDiff, objectiveLowerDim)
    autodiff_rd = AutomaticDiffTune(:ReverseDiff, objectiveLowerDim)
    autodiff_zyg = AutomaticDiffTune(:Zygote, objectiveLowerDim)
    fwd = DiffObjective(objectiveLowerDim, autodiff_fd)
    rd = DiffObjective(objectiveLowerDim, autodiff_rd)
    zyg = DiffObjective(objectiveLowerDim, autodiff_zyg)
    theta_unconstrained = randn(_RNG, length(objectiveLowerDim))
    ## Compute Diffresult
    ld1 = log_density(fwd, theta_unconstrained)
    ld2 = log_density(rd, theta_unconstrained)
    ld3 = log_density(zyg, theta_unconstrained)
    _ld1 = _log_density(objectiveLowerDim, autodiff_fd, theta_unconstrained)
    _ld2 = _log_density(objectiveLowerDim, autodiff_rd, theta_unconstrained)
    _ld3 = _log_density(objectiveLowerDim, autodiff_zyg, theta_unconstrained)
    grad1 = log_density_and_gradient(fwd, theta_unconstrained)
    grad2 = log_density_and_gradient(rd, theta_unconstrained)
    grad3 = log_density_and_gradient(zyg, theta_unconstrained)
    ## Compute manual call ~ Already checked for equality
    ld = objectiveLowerDim(theta_unconstrained)
    grad_mod_fd = ForwardDiff.gradient(objectiveLowerDim, theta_unconstrained)
    grad_mod_rd = ReverseDiff.gradient(objectiveLowerDim, theta_unconstrained)
    grad_mod_zy = Zygote.gradient(objectiveLowerDim, theta_unconstrained)[1]
    ## Compare results
    @test ld - ld1.ℓθᵤ ≈ 0 atol = _TOL
    @test ld - ld2.ℓθᵤ ≈ 0 atol = _TOL
    @test ld - ld3.ℓθᵤ ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_fd - grad1.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_rd - grad2.∇ℓθᵤ)) ≈ 0 atol = _TOL
    @test sum(abs.(grad_mod_zy - grad3.∇ℓθᵤ)) ≈ 0 atol = _TOL
    ## Checks
    _output = check_gradients(_RNG, objectiveLowerDim, [:ForwardDiff, :ReverseDiff, :Zygote]; printoutput = false)
    @test sum(abs.(_output.ℓobjective_gradient_diff)) ≈ 0 atol = _TOL
    ## Results
    ℓDensityResult(objectiveLowerDim, theta_unconstrained)
    ℓDensityResult(objectiveLowerDim)
    ℓDensityResult(fwd)
    ℓDensityResult(rd)
    ℓDensityResult(zyg)
    ℓGradientResult(grad1.θᵤ , grad1.ℓθᵤ , grad1.∇ℓθᵤ)
    ## Update DiffTune
    ModelWrappers.update(autodiff_fd, objectiveLowerDim)
    ModelWrappers.update(autodiff_rd, objectiveLowerDim)
    ModelWrappers.update(autodiff_zyg, objectiveLowerDim)
    ## Config DiffTune
    theta_unconstrained2 = randn(_RNG, length(objectiveLowerDim))
    ModelWrappers._config(ModelWrappers.ADForward(), objectiveLowerDim, theta_unconstrained2)
    ModelWrappers._config(ModelWrappers.ADReverse(), objectiveLowerDim, theta_unconstrained2)
    ModelWrappers._config(ModelWrappers.ADReverseUntaped(), objectiveLowerDim, theta_unconstrained2)
    ModelWrappers._config(ModelWrappers.ADZygote(), objectiveLowerDim, theta_unconstrained2)
end

############################################################################################
# Differentiation - Float32
modelExample2 = ModelWrapper(ExampleModel(), _val_examplemodel, FlattenDefault(; output = Float32))
objectiveExample2 = Objective(modelExample2, (data1, data2, data3, _idx))
objectiveExample2(objectiveExample2.model.val)

fwd = DiffObjective(objectiveExample2, AutomaticDiffTune(:ForwardDiff, objectiveExample2))
rd = DiffObjective(objectiveExample2, AutomaticDiffTune(:ReverseDiff, objectiveExample2))
rd2 = DiffObjective(objectiveExample2, AutomaticDiffTune(:ReverseDiffUntaped, objectiveExample2))
zyg = DiffObjective(objectiveExample2, AutomaticDiffTune(:Zygote, objectiveExample2))

@testset "AutoDiffContainer - Float32 compatibility" begin
    T = Float32
    theta_unconstrained = randn(T, length(objectiveExample2))
    ## Compute Diffresult
    ld1 = log_density(fwd, theta_unconstrained)
    ld2 = log_density(rd, theta_unconstrained)
    ld22 = log_density(rd2, theta_unconstrained)
    ld3 = log_density(zyg, theta_unconstrained)
    grad1 = log_density_and_gradient(fwd, theta_unconstrained)
    grad2 = log_density_and_gradient(rd, theta_unconstrained)
    grad22 = log_density_and_gradient(rd2, theta_unconstrained)
    grad3 = log_density_and_gradient(zyg, theta_unconstrained)
    ## Compare types
    @test ld1.ℓθᵤ isa T && eltype(ld1.θᵤ) == T
    @test ld2.ℓθᵤ isa T && eltype(ld2.θᵤ) == T
    @test ld22.ℓθᵤ isa T && eltype(ld22.θᵤ) == T
    @test ld3.ℓθᵤ isa T && eltype(ld3.θᵤ) == T

    @test grad1.ℓθᵤ isa T && eltype(grad1.θᵤ) == eltype(grad1.∇ℓθᵤ) == T
    @test grad2.ℓθᵤ isa T && eltype(grad2.θᵤ) == eltype(grad2.∇ℓθᵤ) == T
    @test grad22.ℓθᵤ isa T && eltype(grad22.θᵤ) == eltype(grad22.∇ℓθᵤ) == T
    @test grad3.ℓθᵤ isa T && eltype(grad3.θᵤ) == eltype(grad3.∇ℓθᵤ) == T
end

############################################################################################
#Tune Analytic
function fun1(objective::Objective{<:ModelWrapper{M}}, θᵤ::AbstractVector{T}) where {M<:ExampleModel, T<:Real}
    return zeros(size(θᵤ))
end
θᵤ = randn(_RNG, length(objectiveExample))
fun1(objectiveExample, θᵤ)
@testset "AnalyticDiffTune - " begin
    tune_analytic = AnalyticalDiffTune(fun1)
    ModelWrappers.update(tune_analytic, objectiveExample)
    _ld = _log_density(objectiveExample, tune_analytic, θᵤ)
    _ldg =_log_density_and_gradient(objectiveExample, tune_analytic, θᵤ)
    @test _ld == _ldg[1]
    _ldgresult = log_density_and_gradient(objectiveExample, tune_analytic, θᵤ)
    @test _ld == _ldgresult.ℓθᵤ
    @test all(_ldgresult.θᵤ .== θᵤ)
end


text:
module BoundaryTriangulationsTests

using Test
using Gridap.Helpers
using Gridap.Fields
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using FillArrays

domain = (0,4,0,4)
partition = (2,2)
model = simplexify(CartesianDiscreteModel(domain,partition))

btrian = BoundaryTriangulation(model,tags=[7,8])
test_triangulation(btrian)
@test get_background_triangulation(btrian) === get_triangulation(model)

face_s_q = get_cell_ref_map(btrian)

s1 = Point(0.0)
s2 = Point(0.5)
s = [s1,s2]
face_to_s = Fill(s,length(face_s_q))

face_to_q = lazy_map(evaluate,face_s_q,face_to_s)
@test isa(face_to_q,Geometry.FaceCompressedVector)

cell_shapefuns = get_cell_shapefuns(btrian.cell_trian)
cell_grad_shapefuns = lazy_map(Broadcasting(∇),cell_shapefuns)

face_shapefuns = lazy_map(Reindex(cell_shapefuns),get_cell_to_bgcell(btrian))
face_grad_shapefuns = lazy_map(Reindex(cell_grad_shapefuns),get_cell_to_bgcell(btrian))

face_shapefuns_q = lazy_map(evaluate,face_shapefuns,face_to_q)
test_array(face_shapefuns_q,collect(face_shapefuns_q))
@test isa(face_shapefuns_q,Geometry.FaceCompressedVector)

face_grad_shapefuns_q = lazy_map(evaluate,face_grad_shapefuns,face_to_q)
test_array(face_grad_shapefuns_q,collect(face_grad_shapefuns_q))
@test isa(face_grad_shapefuns_q,Geometry.FaceCompressedVector)

face_to_nvec = get_facet_normal(btrian)
face_to_nvec_s = lazy_map(evaluate,face_to_nvec,face_to_s)
test_array(face_to_nvec_s,collect(face_to_nvec_s))

##print_op_tree(face_shapefuns)
#print_op_tree(face_to_q)
#print_op_tree(face_shapefuns_q)
#print_op_tree(face_grad_shapefuns_q)
#print_op_tree(face_to_nvec_s)

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

btrian = BoundaryTriangulation(model)
test_triangulation(btrian)

face_s_q = get_cell_ref_map(btrian)

s1 = Point(0.0)
s2 = Point(0.5)
s = [s1,s2]
face_to_s = Fill(s,length(face_s_q))

face_to_q = lazy_map(evaluate,face_s_q,face_to_s)
@test isa(face_to_q,Geometry.FaceCompressedVector)

cell_shapefuns = get_cell_shapefuns(btrian.cell_trian)
cell_grad_shapefuns = lazy_map(Broadcasting(∇),cell_shapefuns)

face_shapefuns = lazy_map(Reindex(cell_shapefuns),get_cell_to_bgcell(btrian))
face_grad_shapefuns = lazy_map(Reindex(cell_grad_shapefuns),get_cell_to_bgcell(btrian))

face_shapefuns_q = lazy_map(evaluate,face_shapefuns,face_to_q)
test_array(face_shapefuns_q,collect(face_shapefuns_q))
@test isa(face_shapefuns_q,Geometry.FaceCompressedVector)

face_grad_shapefuns_q = lazy_map(evaluate,face_grad_shapefuns,face_to_q)
test_array(face_grad_shapefuns_q,collect(face_grad_shapefuns_q))
@test isa(face_grad_shapefuns_q,Geometry.FaceCompressedVector)

face_to_nvec = get_facet_normal(btrian)
face_to_nvec_s = lazy_map(evaluate,face_to_nvec,face_to_s)
test_array(face_to_nvec_s,collect(face_to_nvec_s))
@test isa(face_to_nvec_s,Geometry.FaceCompressedVector)

#print_op_tree(face_shapefuns_q)
#print_op_tree(face_grad_shapefuns_q)
#print_op_tree(face_to_nvec_s)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))

s2q = get_cell_ref_map(btrian)

q = lazy_map(evaluate,s2q,s)
r = Vector{Point{2,Float64}}[
  [(0.25,0.0),(0.75,0.0)],[(0.0,0.25),(0.0,0.75)],
  [(0.25,0.0),(0.75,0.0)],[(1.0,0.25),(1.0,0.75)],
  [(0.25,1.0),(0.75,1.0)],[(0.0,0.25),(0.0,0.75)],
  [(0.25,1.0),(0.75,1.0)],[(1.0,0.25),(1.0,0.75)]]
test_array(q,r)

q2x = get_cell_map(btrian.cell_trian)
s2x = lazy_map(∘,lazy_map(Reindex(q2x),get_cell_to_bgcell(btrian)),s2q)

x = lazy_map(evaluate,s2x,s)
r = Vector{Point{2,Float64}}[
  [(0.5,0.0),(1.5,0.0)],[(0.0,0.5),(0.0,1.5)],
  [(2.5,0.0),(3.5,0.0)],[(4.0,0.5),(4.0,1.5)],
  [(0.5,4.0),(1.5,4.0)],[(0.0,2.5),(0.0,3.5)],
  [(2.5,4.0),(3.5,4.0)],[(4.0,2.5),(4.0,3.5)]]
test_array(x,r)

nvec = get_facet_normal(btrian)

nvec_s = lazy_map(evaluate,nvec,s)
r = Vector{Point{2,Float64}}[
  [(0.0,-1.0),(0.0,-1.0)],[(-1.0,0.0),(-1.0,0.0)],
  [(0.0,-1.0),(0.0,-1.0)],[(1.0,-0.0),(1.0,-0.0)],
  [(0.0,1.0),(0.0,1.0)],[(-1.0,0.0),(-1.0,0.0)],
  [(0.0,1.0),(0.0,1.0)],[(1.0,-0.0),(1.0,-0.0)]]
test_array(nvec_s,r)

cellids = collect(1:num_cells(model))

face_to_cellid = lazy_map(Reindex(cellids),get_cell_to_bgcell(btrian))
@test face_to_cellid == get_cell_to_bgcell(btrian)

trian = get_background_triangulation(btrian)
q2x = get_cell_map(trian)
s2q = get_cell_ref_map(btrian)
s2x = lazy_map(∘,lazy_map(Reindex(q2x),get_cell_to_bgcell(btrian)),s2q)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))
x = lazy_map(evaluate,s2x,s)
r = Vector{Point{2,Float64}}[
  [(0.5,0.0),(1.5,0.0)],[(0.0,0.5),(0.0,1.5)],
  [(2.5,0.0),(3.5,0.0)],[(4.0,0.5),(4.0,1.5)],
  [(0.5,4.0),(1.5,4.0)],[(0.0,2.5),(0.0,3.5)],
  [(2.5,4.0),(3.5,4.0)],[(4.0,2.5),(4.0,3.5)]]
test_array(x,r)

nvec = get_facet_normal(btrian)
nvec_x = lazy_map(evaluate,nvec,s)

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

btrian = BoundaryTriangulation(model,tags="tag_8")
test_triangulation(btrian)
@test all(1 .== btrian.glue.face_to_lcell)

btrian = BoundaryTriangulation(model,get_face_labeling(model),tags="tag_8")
test_triangulation(btrian)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))
nvec = get_facet_normal(btrian)
nvec_x = lazy_map(evaluate,nvec,s)
s2x = get_cell_map(btrian)
x = lazy_map(evaluate,s2x,s)

oldbtrian = BoundaryTriangulation(model)
bface_to_oldbface = collect(1:4)
btrian = RestrictedTriangulation(oldbtrian,bface_to_oldbface)
test_triangulation(btrian)

nb = get_facet_normal(btrian)
@test length(nb) == num_cells(btrian)

btrian = BoundaryTriangulation(model,tags=["tag_6","tag_7"])
test_triangulation(btrian)

end # module


text:
using SymPy
using Test
using LinearAlgebra
using Plots
using Calculus

# using SymbolicUtils

include("../../lib/KTBC.jl")
using .KTBC: CreateKTBC

using Pkg
Pkg.add("Pluto")
Pluto.run()

# 1 Symbolic path
@vars u₁ u₂ x₁ y₁

u⁺ = [u₁ u₂]
u = [
  u₁
  u₂
]

T₂ = convert(Array{Int}, CreateKTBC(2)[2])

u⁺T₂u = u⁺ * T₂ * u

result₁ = [u₁ * (u₁ - u₂) + u₂ * (-u₁ + 2 * u₂)]

@test u⁺T₂u[1] - result₁[1] == 0

# 3
A = [
  1 -1 0
  0 1 -1
  -1 0 1
]

plot3d(A)

C = (A') * A

plot!(C)

@test_throws PosDefException cholesky(C)

@test A[1:3] + A[4:6] + A[7:9] == [0; 0; 0]

# or
u = ones(Int, 3, 1)
z = zeros(Int, 3, 1)

@test A * u == z


text:
module Simulatr

using StatsBase
using Distributions
using NamedTuples

import StatsBase: sample
import Distributions: Uniform
import NamedTuples: @NT

export simrel

include("simrel.jl")

end # module

max_stars_repo_name: ExtraSpace/JuliaSimulatr

max_issues_repo_name: ExtraSpace/JuliaSimulatr

max_forks_repo_name: ExtraSpace/JuliaSimulatr


text:
using RoadRunner
using Test

# Load the contents of a SBML file into a string variable
sbmlFile = "feedback.xml"
f = open(sbmlFile)
sbmlStr = read(f,String)
close(f)

rr = RoadRunner.createRRInstance()         # Start up roadRunner
RoadRunner.loadSBML(rr, sbmlStr)

@testset "parameters" begin
    @test RoadRunner.getNumberOfGlobalParameters(rr) == 5
    @test RoadRunner.getGlobalParameterIds(rr) == ["J0_VM1", "J0_Keq1", "J0_h", "J4_V4", "J4_KS4"]
end

@testset "species" begin
    @test RoadRunner.getNumberOfFloatingSpecies(rr) == 4
    @test RoadRunner.getFloatingSpeciesIds(rr) == ["S1", "S2", "S3", "S4"]
    @test RoadRunner.getNumberOfBoundarySpecies(rr) == 2
    @test RoadRunner.getBoundarySpeciesIds(rr) == ["X0", "X1"]
end

@testset "steadystates" begin
    @test RoadRunner.steadyState(rr) < 1e-6
end


text:
"""
    tty_has_color()

Check if TTY supports color. This is mainly for lower
Julia version like 1.0.
"""
function tty_has_color()
    if isdefined(Base, :get_have_color)
        return Base.get_have_color()
    else
        return Base.have_color
    end
end

"""
    supports_color256()

Check if the terminal supports color 256.
"""
function supports_color256()
    haskey(ENV, "TERM") || return false
    try
        return parse(Int, readchomp(`tput colors 0`)) == 256
    catch
        return false
    end
end

const ColorType = Union{Int, Symbol}

"""
    ColorPreference

The color preference type.
"""
struct ColorPreference
    fieldname::ColorType
    type::ColorType
    operator::ColorType

    # literal-like
    literal::ColorType
    constant::ColorType
    number::ColorType
    string::ColorType

    # comment-like
    comment::ColorType
    undef::ColorType
    linenumber::ColorType
end

struct PreferenceInvalid <: Exception
    key::String
    type
    got
end

function Base.showerror(io::IO, x::PreferenceInvalid)
    print(io, "preference for ")
    printstyled(io, x.key; color=:light_blue)
    print(io, " is invalid, expect ")
    printstyled(io, x.type; color=:green)
    print(io, " got: ", repr(x.got))
end

Base.show(io::IO, x::ColorPreference) = pprint_struct(io, x)

"""
    ColorPreference(;kw...)

See [`pprint`](@ref) for available keyword configurations.
"""
function ColorPreference(;kw...)
    colors = supports_color256() ? default_colors_256() : default_colors_ansi()
    if color_prefs_toml !== nothing
        merge!(colors, color_prefs_toml)
    end
    colors = merge!(colors, kw)

    args = map(fieldnames(ColorPreference)) do name
        val = colors[string(name)]
        if val isa String
            return Symbol(val)
        elseif val isa Int
            return val
        else
            throw(PreferenceInvalid("GarishPrint.color.$name", Union{String, Int}, val))
        end
    end
    return ColorPreference(args...)
end

"""
    default_colors_ansi()

The default ANSI color theme.
"""
function default_colors_ansi()
    Dict{String, Any}(
        "fieldname" => "light_blue",
        "type" => "green",
        "operator" => "normal",
        "literal" => "yellow",
        "constant" => "yellow",
        "number" => "normal",
        "string" => "yellow",
        "comment" => "light_black",
        "undef" => "normal",
        "linenumber" => "light_black",
    )
end

"""
    default_colors_256()

The default color 256 theme.
"""
function default_colors_256()
    Dict{String, Any}(
        "fieldname" => 039,
        "type" => 037,
        "operator" => 196,
        "literal" => 140,
        "constant" => 099,
        "number" => 140,
        "string" => 180,
        "comment" => 240,
        # undef is actually a constant
        "undef" => 099,
        "linenumber" => 240,
    )
end