text:
@testset "EmbedOnehot" begin
    multi(x, n) = collect(Iterators.repeated(x, n))

    after_enc = [2,2,4,5,6,1,1]
    multi_after_enc = hcat(multi(after_enc, 5)...)

    e = Embed(10, 6)

    @test e(after_enc) == hcat(map(i->e.embedding[:,i], after_enc)...)
    @test eltype(e(after_enc, 0.5)) == Float32

    @test size(e) == (10, 6)
    @test size(e, 1) == 10
    @test size(e, 2) == 6
end

ext: jl


text:
# Implements the APIs required by an application master to communicate with the resource manager and the node managers.
# The underlying methods are implementations of the following protobuf services:
# - applicationmaster_protocol.proto
# - containermanagement_protocol.proto

# In managed mode, the AMRMToken available in file CONTAINER_TOKEN_FILE_ENV_NAME is in Java serialized format (https://apache.googlesource.com/hadoop-common/+/HADOOP-6685/src/java/org/apache/hadoop/security/Credentials.java).
# Since that is unusable, we are forced to operate only in unmanaged mode, where we get it from the application report.

const YARN_CONTAINER_MEM_DEFAULT = 128
const YARN_CONTAINER_CPU_DEFAULT = 1
const YARN_CONTAINER_LOCATION_DEFAULT = "*"
const YARN_CONTAINER_PRIORITY_DEFAULT = 1
const YARN_NM_CONN_KEEPALIVE_SECS = 5*60

@doc doc"""
YarnAppMaster is a skeleton application master. It provides the generic scafolding methods which can be used to create specific
application masters for different purposes.
""" ->
type YarnAppMaster
    amrm_conn::YarnAMRMProtocol

    host::AbstractString
    port::Integer
    tracking_url::AbstractString

    max_mem::Int32
    max_cores::Int32
    available_mem::Int32
    available_cores::Int32
    nodes::YarnNodes
    containers::YarnContainers

    queue::AbstractString
    managed::Bool

    response_id::Int32

    registration::Nullable{RegisterApplicationMasterResponseProto}
    lck::Lock

    function YarnAppMaster(rmhost::AbstractString, rmport::Integer, ugi::UserGroupInformation,
                amhost::AbstractString="", amport::Integer=0, amurl::AbstractString="")
        amrm_conn = YarnAMRMProtocol(rmhost, rmport, ugi)
        lck = makelock()
        put!(lck, 1)

        new(amrm_conn, amhost, amport, amurl,
            Int32(0), Int32(0), Int32(0), Int32(0),
            YarnNodes(ugi), YarnContainers(),
            "", true, 1,
            Nullable{RegisterApplicationMasterResponseProto}(), lck)
    end
end

macro locked(yam, ex)
    quote
        try
            take!(yam.lck)
            $(esc(ex))
        finally
            put!(yam.lck, 1)
        end
    end
end

function show(io::IO, yam::YarnAppMaster)
    show(io, yam.amrm_conn)
    println(io, "    Memory: available:$(yam.available_mem), max:$(yam.max_mem), can schecule:$(can_schedule_mem(yam))")
    println(io, "    Cores: available:$(yam.available_cores), max:$(yam.max_cores), can schedule:$(can_schedule_cores(yam))")
    println(io, "    Queue: $(yam.queue)")
    show(io, yam.nodes)
    show(io, yam.containers)
end

callback(yam::YarnAppMaster, on_container_alloc::Nullable, on_container_finish::Nullable) =
    callback(yam.containers, on_container_alloc, on_container_finish)

function submit(client::YarnClient, unmanagedappmaster::YarnAppMaster)
    @logmsg("submitting unmanaged application")
    clc = launchcontext()
    app = submit(client, clc, YARN_CONTAINER_MEM_DEFAULT, YARN_CONTAINER_CPU_DEFAULT; unmanaged=true)

    # keep the am_rm token
    tok = am_rm_token(app)
    channel = unmanagedappmaster.amrm_conn.channel
    add_token(channel.ugi, token_alias(channel), tok)

    # register the unmanaged appmaster
    unmanagedappmaster.managed = false
    register(unmanagedappmaster)
    wait_for_attempt_state(app, Int32(1), YarnApplicationAttemptStateProto.APP_ATTEMPT_RUNNING) || throw(YarnException("Application attempt could not be launched"))

    # initialize complete node list for appmaster
    nodes(client; nodelist=unmanagedappmaster.nodes)

    # start the am_rm processing task once the app attempt starts running if it is an unmanaged application master
    @async(process_am_rm(unmanagedappmaster))
    app
end

function register(yam::YarnAppMaster)
    inp = RegisterApplicationMasterRequestProto()
    if !isempty(yam.host) && (yam.port > 0)
        set_field!(inp, :host, yam.host)
        set_field!(inp, :port, convert(Int32,yam.port))
    end
    !isempty(yam.tracking_url) && set_field!(inp, :tracking_url, yam.tracking_url)

    resp = @locked(yam, registerApplicationMaster(yam.amrm_conn, inp))
    yam.registration = Nullable(resp)

    if isfilled(resp, :maximumCapability)
        yam.max_mem = resp.maximumCapability.memory
        yam.max_cores = resp.maximumCapability.virtual_cores
    end
    @logmsg("max capability: mem:$(yam.max_mem), cores:$(yam.max_cores)")
    if isfilled(resp, :queue)
        yam.queue = resp.queue
    end

    # start the am_rm processing task on registration if it is a managed application master
    yam.managed && @async(process_am_rm(yam))
    nothing
end

can_schedule_mem(yam::YarnAppMaster) = can_schedule(yam, SchedulerResourceTypes.MEMORY)
can_schedule_cores(yam::YarnAppMaster) = can_schedule(yam, SchedulerResourceTypes.CPU)
function can_schedule(yam::YarnAppMaster, restype::Int32)
    reg = yam.registration
    !isfilled(reg, :scheduler_resource_types) && (return false)
    restype in reg.scheduler_resource_types
end

function _unregister(yam::YarnAppMaster, finalstatus::Int32, diagnostics::AbstractString)
    inp = protobuild(FinishApplicationMasterRequestProto, Dict(:final_application_status => finalstatus))
    !isempty(yam.tracking_url) && set_field!(inp, :tracking_url, yam.tracking_url)
    !isempty(diagnostics) && set_field!(inp, :diagnostics, diagnostics)

    resp = @locked(yam, finishApplicationMaster(yam.amrm_conn, inp))
    resp.isUnregistered && (yam.registration = Nullable{RegisterApplicationMasterResponseProto}())
    resp.isUnregistered
end

unregister(yam::YarnAppMaster, success::Bool, diagnostics::AbstractString="") = _unregister(yam, success ? FinalApplicationStatusProto.APP_SUCCEEDED : FinalApplicationStatusProto.APP_FAILED, diagnostics)
kill(yam::YarnAppMaster, diagnostics::AbstractString="") = _unregister(yam, FinalApplicationStatusProto.APP_KILLED, diagnostics)

container_allocate(yam::YarnAppMaster, numcontainers::Int; opts...) = request_alloc(yam.containers, numcontainers; opts...)
container_release(yam::YarnAppMaster, cids::ContainerIdProto...) = request_release(yam.containers, cids...)

container_start(yam::YarnAppMaster, cid::ContainerIdProto, container_spec::ContainerLaunchContextProto) = container_start(yam, yam.containers.containers[cid], container_spec)
function container_start(yam::YarnAppMaster, container::ContainerProto, container_spec::ContainerLaunchContextProto)
    @logmsg("starting container $(container)")
    req = protobuild(StartContainerRequestProto, Dict(:container_launch_context => container_spec, :container_token => container.container_token))
    inp = protobuild(StartContainersRequestProto, Dict(:start_container_request => [req]))

    nodeid = container.nodeId
    nm_conn = get_connection(yam.nodes, nodeid)
    success = false
    resp = StartContainersResponseProto()
    try
        resp = startContainers(nm_conn, inp)
        success = true
    finally
        release_connection(yam.nodes, nodeid, nm_conn, success)
    end
    success || throw(YarnException("Error starting container"))
    isfilled(resp, :succeeded_requests) || throw(YarnException(isfilled(resp,:failed_requests) ? resp.failed_requests[1] : "Error starting container"))
    cid = resp.succeeded_requests[1]
    (cid == container.id) || throw(YarnException("Unexpected container id mismatch"))
    set_busy(yam.containers, cid)
    cid
end

container_stop(yam::YarnAppMaster, cid::ContainerIdProto) = container_stop(yam, yam.containers.containers[cid])
function container_stop(yam::YarnAppMaster, container::ContainerProto)
    @logmsg("stopping container $container")

    inp = protobuild(StopContainersRequestProto, Dict(:container_id => [container.id]))
    nodeid = container.nodeId
    nm_conn = get_connection(yam.nodes, nodeid)
    success = false
    resp = StopContainersResponseProto()
    try
        resp = stopContainers(nm_conn, inp)
        success = true
    finally
        release_connection(yam.nodes, nodeid, nm_conn, success)
    end
    success || throw(YarnException("Error stopping container"))
    isfilled(resp, :succeeded_requests) || throw(YarnException(isfilled(resp,:failed_requests) ? resp.failed_requests[1] : "Error stopping container"))
    cid = resp.succeeded_requests[1]
    (cid == container.id) || throw(YarnException("Unexpected container id mismatch"))
    set_free(yam.containers, cid)
    cid
end

function _update_rm(yam::YarnAppMaster)
    @logmsg("started processing am-rm messages")
    inp = AllocateRequestProto()

    # allocation and release requests
    (alloc_pending,release_pending) = torequest(yam.containers)
    @logmsg("alloc pending: $alloc_pending")
    @logmsg("release pending: $release_pending")
    !isempty(alloc_pending) && set_field!(inp, :ask, alloc_pending)
    !isempty(release_pending) && set_field!(inp, :release, release_pending)

    # send one-up response id
    if yam.response_id == typemax(Int32)
        yam.response_id = 1
    else
        yam.response_id += 1
    end
    set_field!(inp, :response_id, yam.response_id)

    @logmsg(inp)
    resp = @locked(yam, allocate(yam.amrm_conn, inp))
    @logmsg(resp)

    # store/update tokens
    channel = yam.amrm_conn.channel
    ugi = channel.ugi
    isfilled(resp, :am_rm_token) && add_token(ugi, token_alias(channel), resp.am_rm_token)
    if isfilled(resp, :nm_tokens)
        for nmtok in resp.nm_tokens
            add_token(ugi, token_alias(nmtok.nodeId), nmtok.token)
        end
    end

    # update available headroom
    if isfilled(resp, :limit)
        yam.available_mem = resp.limit.memory
        yam.available_cores = resp.limit.virtual_cores
    end

    # update node and container status
    update(yam.nodes, resp)
    update(yam.containers, resp)
    @logmsg("finished processing am-rm messages")
    @logmsg(yam)
    nothing
end

haverequests(yam::YarnAppMaster) = haverequests(yam.containers)

function process_am_rm(yam::YarnAppMaster)
    @logmsg("started am-rm processor task")
    stopped = ()->isnull(yam.registration)
    stopwaiting = ()->(haverequests(yam) || isnull(yam.registration))
    waittime = 10.
    while !stopped()
        t1 = time()
        _update_rm(yam)
        t2 = time()
        interval = t1 + 60 - t2
        if interval > 0
            timedwait(stopwaiting, interval; pollint=5.0)
            t3 = time()
            waittime = (waittime*5 + (t3 - t2))/6
            (waittime < 10) && sleep(10-waittime)
        end
    end
    @logmsg("stopped am-rm processor task")
    nothing
end


text:
# Decision trees as found in DecisionTree Julia package.
module DecisionTreeLearners

import DecisionTree
DT = DecisionTree

# standard included modules
using DataFrames
using Random
using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils

import AutoMLPipeline.AbsTypes: fit!, transform!
export fit!, transform!

export PrunedTree, RandomForest, Adaboost

# Pruned CART decision tree.

"""
    PrunedTree(
      Dict(
        :purity_threshold => 1.0,
        :max_depth => -1,
        :min_samples_leaf => 1,
        :min_samples_split => 2,
        :min_purity_increase => 0.0
      )
    )

Decision tree classifier.
See [DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparmeters:
- `:purity_threshold` => 1.0 (merge leaves having >=thresh combined purity)
- `:max_depth` => -1 (maximum depth of the decision tree)
- `:min_samples_leaf` => 1 (the minimum number of samples each leaf needs to have)
- `:min_samples_split` => 2 (the minimum number of samples in needed for a split)
- `:min_purity_increase` => 0.0 (minimum purity needed for a split)

Implements `fit!`, `transform!`
"""
mutable struct PrunedTree <: Learner
  name::String
  model::Dict
  args::Dict

  function PrunedTree(args=Dict())
    default_args = Dict(
      :name => "prunetree",
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_args => Dict(
        # Merge leaves having >= purity_threshold CombineMLd purity.
        :purity_threshold => 1.0,
        # Maximum depth of the decision tree (default: no maximum).
        :max_depth => -1,
        # Minimum number of samples each leaf needs to have.
        :min_samples_leaf => 1,
        # Minimum number of samples in needed for a split.
        :min_samples_split => 2,
        # Minimum purity needed for a split.
        :min_purity_increase => 0.0
      )
    )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

"""
    fit!(tree::PrunedTree, features::DataFrame, labels::Vector)

Optimize the hyperparameters of `PrunedTree` instance.
"""
function fit!(ptree::PrunedTree, features::DataFrame, labels::Vector)
  instances=convert(Matrix,features)
  args = ptree.args[:impl_args]
  btreemodel = DT.build_tree(
    labels,
    instances,
    0, # num_subfeatures (keep all)
    args[:max_depth],
    args[:min_samples_leaf],
    args[:min_samples_split],
    args[:min_purity_increase])
  btreemodel = DT.prune_tree(btreemodel, args[:purity_threshold])
  ptree.model = Dict(
                    :dtmodel => btreemodel,
                    :impl_args => args
                   )
end

"""
    transform!(ptree::PrunedTree, features::DataFrame)

Predict using the optimized hyperparameters of the trained `PrunedTree` instance.
"""
function transform!(ptree::PrunedTree, features::DataFrame)
  instances=convert(Matrix,features)
  model = ptree.model[:dtmodel]
  return DT.apply_tree(model, instances)
end

# Random forest (CART).

"""
    RandomForest(
      Dict(
        :output => :class,
        :num_subfeatures => 0,
        :num_trees => 10,
        :partial_sampling => 0.7,
        :max_depth => -1
      )
    )

Random forest classification.
See [DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparmeters:
- `:num_subfeatures` => 0  (number of features to consider at random per split)
- `:num_trees` => 10 (number of trees to train)
- `:partial_sampling` => 0.7 (fraction of samples to train each tree on)
- `:max_depth` => -1 (maximum depth of the decision trees)
- `:min_samples_leaf` => 1 (the minimum number of samples each leaf needs to have)
- `:min_samples_split` => 2 (the minimum number of samples in needed for a split)
- `:min_purity_increase` => 0.0 (minimum purity needed for a split)

Implements `fit!`, `transform!`
"""
mutable struct RandomForest <: Learner
  name::String
  model::Dict
  args::Dict
  function RandomForest(args=Dict())
    default_args = Dict(
      :name => "rf",
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_args => Dict(
        # Number of features to train on with trees (default: 0, keep all).
        :num_subfeatures => 0,
        # Number of trees in forest.
        :num_trees => 10,
        # Proportion of trainingset to be used for trees.
        :partial_sampling => 0.7,
        # Maximum depth of each decision tree (default: no maximum).
        :max_depth => -1
      )
    )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

"""
    fit!(forest::RandomForest, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}

Optimize the parameters of the `RandomForest` instance.
"""
function fit!(forest::RandomForest, features::DataFrame, labels::Vector)
  instances=convert(Matrix,features)
  # Set training-dependent options
  impl_args = forest.args[:impl_args]
  # Build model
  model = DT.build_forest(
    labels,
    instances,
    impl_args[:num_subfeatures],
    impl_args[:num_trees],
    impl_args[:partial_sampling],
    impl_args[:max_depth]
  )
  forest.model = Dict(
                      :dtmodel => model,
                      :impl_args => impl_args
                     )
end

"""
    transform!(forest::RandomForest, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Predict using the optimized hyperparameters of the trained `RandomForest` instance.
"""
function transform!(forest::RandomForest, features::DataFrame)
  instances = features
  instances=convert(Matrix,features)
  model = forest.model[:dtmodel]
  return DT.apply_forest(model, instances)
end

# Adaboosted decision stumps.

"""
    Adaboost(
      Dict(
        :output => :class,
        :num_iterations => 7
      )
    )

Adaboosted decision tree stumps. See
[DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:
- `:num_iterations` => 7 (number of iterations of AdaBoost)

Implements `fit!`, `transform!`
"""
mutable struct Adaboost <: Learner
  name::String
  model::Dict
  args::Dict
  function Adaboost(args=Dict())
    default_args = Dict(
      :name => "adaboost",
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_args => Dict(
        # Number of boosting iterations.
        :num_iterations => 7
      )
    )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

"""
    fit!(adaboost::Adaboost, features::DataFrame, labels::Vector)

Optimize the hyperparameters of `Adaboost` instance.
"""
function fit!(adaboost::Adaboost, features::DataFrame, labels::Vector)
  instances = convert(Matrix,features)
  # NOTE(svs14): Variable 'model' renamed to 'ensemble'.
  #              This differs to DecisionTree
  #              official documentation to avoid confusion in variable
  #              naming within CombineML.
  ensemble, coefficients = DT.build_adaboost_stumps(
    labels, instances, adaboost.args[:impl_args][:num_iterations]
  )
  adaboost.model = Dict(
    :ensemble => ensemble,
    :coefficients => coefficients
  )
end

"""
    transform!(adaboost::Adaboost, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Predict using the optimized hyperparameters of the trained `Adaboost` instance.
"""
function transform!(adaboost::Adaboost, features::DataFrame)::Vector{<:Any}
  instances = convert(Matrix,features)
  return DT.apply_adaboost_stumps(
    adaboost.model[:ensemble], adaboost.model[:coefficients], instances
  )
end

end # module


text:
# # Hydrostatic Equilibrium test for Richards Equation

# This tutorial shows how to use `ClimateMachine` code to solve
# Richards equation in a column of soil. We choose boundary
# conditions of zero flux at the top and bottom of the column,
# and then run the simulation long enough to see that the system
# is approaching hydrostatic equilibrium, where the gradient of the
# pressure head is equal and opposite the gradient of the
# gravitational head. Note that the [`SoilWaterModel`](@ref
# ClimateMachine.Land.SoilWaterModel) includes
# a prognostic equation for the volumetric ice fraction,
# as ice is a form of water that must be accounted for to ensure
# water mass conservation. If freezing and thawing are not turned on
# (the default), the amount of ice in the model is zero for all space and time
# (again by default).

# The equations are:

# ``
# \frac{ ∂ ϑ_l}{∂ t} = ∇ ⋅ K (T, ϑ_l, θ_i; ν, ...) ∇h( ϑ_l, z; ν, ...).
# ``

# ``
# \frac{ ∂ θ_i}{∂ t} = 0
# ``

# Here

# ``t`` is the time (s),

# ``z`` is the location in the vertical (m),

# ``T`` is the temperature of the soil (K),

# ``K`` is the hydraulic conductivity (m/s),

# ``h`` is the hydraulic head (m),

# ``ϑ_l`` is the augmented volumetric liquid water fraction,

# ``θ_i`` is the volumetric ice fraction, and

# ``ν, ...`` denotes parameters relating to soil type, such as porosity.

# We will solve this equation in an effectively 1-d domain with ``z ∈ [-10,0]``,
# and with the following boundary and initial conditions:

# ``- K ∇h(t, z = 0) = 0 ẑ ``

# `` -K ∇h(t, z = -10) = 0 ẑ``

# `` ϑ(t = 0, z) = ν-0.001 ``

# `` θ_i(t = 0, z) = 0.0. ``

# where ``\nu`` is the porosity.

# A word about the hydraulic conductivity: please see the
# [`hydraulic functions`](./hydraulic_functions.md) tutorial
# for options regarding this function. The user can choose to make it depend
# on the temperature and the amount of ice in the soil; the default, which we use
# here, is that `K` only depends on the liquid moisture content.

# Lastly, our formulation of this equation allows for a continuous solution in both
# saturated and unsaturated areas, following [1].

# # Preliminary setup

# - Load external packages

using MPI
using OrderedCollections
using StaticArrays
using Statistics

# - Load CLIMAParameters and ClimateMachine modules

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state

# - Define the float type desired (`Float64` or `Float32`)
const FT = Float64;

# - Initialize ClimateMachine for CPU
ClimateMachine.init(; disable_gpu = true);

# Load plot helpers:
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(clima_dir, "docs", "plothelpers.jl"));

# # Set up the soil model

# We want to solve Richards equation alone, without simultaneously
# solving the heat equation. Because of that, we choose a
# [`PrescribedTemperatureModel`](@ref
# ClimateMachine.Land.PrescribedTemperatureModel).
# The user can supply a function for temperature,
# depending on time and space; if this option is desired, one could also
# choose to model the temperature dependence of viscosity, or to drive
# a freeze/thaw cycle, for example. If the user simply wants to model
# Richards equation for liquid water, the defaults will allow for that.
# Here we ignore the effects of temperature and freezing and thawing,
# using the defaults.

soil_heat_model = PrescribedTemperatureModel();

# Define the porosity, Ksat, and specific storage values for the soil. Note
# that all values must be given in mks units. The soil parameters chosen
# roughly correspond to Yolo light clay.
soil_param_functions = SoilParamFunctions{FT}(
    porosity = 0.495,
    Ksat = 0.0443 / (3600 * 100),
    S_s = 1e-3,
);

# Define the boundary conditions. The user can specify two conditions,
# either at the top or at the bottom, and they can either be Dirichlet
# (on `ϑ_l`) or Neumann (on `-K∇h`). Note that fluxes are supplied as
# scalars, inside the code they are multiplied by ẑ. The two conditions
# not supplied must be set to `nothing`.

surface_flux = (aux, t) -> eltype(aux)(0.0)
bottom_flux = (aux, t) -> eltype(aux)(0.0)
surface_state = nothing
bottom_state = nothing
bc = GeneralBoundaryConditions(
    Dirichlet(surface_state = surface_state, bottom_state = bottom_state),
    Neumann(surface_flux = surface_flux, bottom_flux = bottom_flux),
)

# Define the initial state function. The default for `θ_i` is zero.
ϑ_l0 = (aux) -> eltype(aux)(0.494);

# Create the SoilWaterModel. The defaults are a temperature independent
# viscosity, and no impedance factor due to ice. We choose to make the
# hydraulic conductivity a function of the moisture content `ϑ_l`,
# and employ the vanGenuchten hydraulic model with `n` = 2.0. The van
# Genuchten parameter `m` is calculated from `n`, and we use the default
# value for `α`.
soil_water_model = SoilWaterModel(
    FT;
    moisture_factor = MoistureDependent{FT}(),
    hydraulics = vanGenuchten{FT}(n = 2.0),
    initialϑ_l = ϑ_l0,
    boundaries = bc,
);

# Create the soil model - the coupled soil water and soil heat models.
m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model);

# We are ignoring sources and sinks here, like runoff or freezing and thawing.
sources = ();

# Define the function that initializes the prognostic variables. This
# in turn calls the functions supplied to `soil_water_model`.
function init_soil_water!(land, state, aux, localgeo, time)
    state.soil.water.ϑ_l = eltype(state)(land.soil.water.initialϑ_l(aux))
    state.soil.water.θ_i = eltype(state)(land.soil.water.initialθ_i(aux))
end

# Create the land model - in this tutorial, it only includes the soil.
m = LandModel(
    param_set,
    m_soil;
    source = sources,
    init_state_prognostic = init_soil_water!,
);

# # Specify the numerical configuration and output data.

# Specify the polynomial order and vertical resolution.
N_poly = 5;
nelem_vert = 20;

# Specify the domain boundaries.
zmax = FT(0);
zmin = FT(-10);

# Create the driver configuration.
driver_config = ClimateMachine.SingleStackConfiguration(
    "LandModel",
    N_poly,
    nelem_vert,
    zmax,
    param_set,
    m;
    zmin = zmin,
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
);

# Choose the initial and final times, as well as a timestep.
t0 = FT(0)
timeend = FT(60 * 60 * 24 * 4)
dt = FT(5);

# Create the solver configuration.
solver_config =
    ClimateMachine.SolverConfiguration(t0, timeend, driver_config, ode_dt = dt);

# Determine how often you want output.
const n_outputs = 5;

const every_x_simulation_time = ceil(Int, timeend / n_outputs);

# Create a place to store this output.
state_types = (Prognostic(), Auxiliary(), GradientFlux())
dons_arr = Dict[dict_of_nodal_states(solver_config, state_types; interp = true)]
time_data = FT[0] # store time data

callback = GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
    dons = dict_of_nodal_states(solver_config, state_types; interp = true)
    push!(dons_arr, dons)
    push!(time_data, gettime(solver_config.solver))
    nothing
end;

# # Run the integration
ClimateMachine.invoke!(solver_config; user_callbacks = (callback,));

# Get the final state and create plots:
dons = dict_of_nodal_states(solver_config, state_types; interp = true)
push!(dons_arr, dons)
push!(time_data, gettime(solver_config.solver));

# Get z-coordinate
z = get_z(solver_config.dg.grid; rm_dupes = true);

# # Create some plots
# We'll plot the moisture content vs depth in the soil, as well as
# the expected slope of `ϑ_l` in hydrostatic equilibrium when the soil
# is saturated. For `ϑ_l` values above porosity, the soil is
# saturated, and the pressure head changes from being equal to the matric
# potential to the pressure generated by compression of water and the soil
# matrix. The pressure head is continuous at porosity, but the derivative
# is not.
slope = -1e-3

output_dir = @__DIR__;

t = time_data ./ (60 * 60 * 24);

plot(
    dons_arr[1]["soil.water.ϑ_l"],
    dons_arr[1]["z"],
    label = string("t = ", string(t[1]), "days"),
    xlim = [0.47, 0.501],
    ylabel = "z",
    xlabel = "ϑ_l",
    legend = :bottomleft,
    title = "Equilibrium test",
);
plot!(
    dons_arr[4]["soil.water.ϑ_l"],
    dons_arr[4]["z"],
    label = string("t = ", string(t[4]), "days"),
);
plot!(
    dons_arr[6]["soil.water.ϑ_l"],
    dons_arr[6]["z"],
    label = string("t = ", string(t[6]), "days"),
);
plot!(
    (dons_arr[1]["z"] .+ 10.0) .* slope .+ dons_arr[6]["soil.water.ϑ_l"][1],
    dons_arr[1]["z"],
    label = "expected",
);

plot!(
    1e-3 .+ dons_arr[1]["soil.water.ϑ_l"],
    dons_arr[1]["z"],
    label = "porosity",
);
# save the output.
savefig(joinpath(output_dir, "equilibrium_test_ϑ_l_vG.png"))
# ![](equilibrium_test_ϑ_l_vG.png)
# # References
# [1] Woodward, C. S., and C. N. Dawson (2000), Analysis of expanded mixed
# finite element methods for a nonlinear parabolic equation modeling
# flow into variably saturated porous media, SIAM J. Numer. Anal., 37, 701–724

max_forks_repo_name: jm-c/CLIMA


text:
export
    tsp_hobo

# functions below assumes that hobos accepts ony one or zero

function are_same(hobos::Vector{HOBO{T}}, n::Int) where T<:Real
    bit_length = length(hobos)
    bits = reverse(digits(n, base=2, pad=bit_length))
    @assert length(bits) == bit_length
    reduce(and, [b == 0 ? one(T) - h : h  for (h,b)=zip(hobos, bits)])
end

function are_different(v1::Vector{HOBO{T}}, v2::Vector{HOBO{T}}) where T<:Real
    neg(reduce(or, xor.(v1, v2)))
end

# assumes numbering from 0 and that the vectors of HOBO
# are of the form var[...]
function smaller_equal_than(hobos::Vector{HOBO{T}}, n::Int) where T<:Real
    k = length(hobos)
    # no point for n = 2^k-1 or n = 0
    @assert n+1 == 2^k || 1 <= n <= 2^k-2
    bits = reverse(digits(n, base=2))
    zero_pos = findall(x->x==0, bits)
    @assert bits[1] == 1
    @assert length(bits) == k

    result = zero(T)*hobos[1]
    for ind=zero_pos
        result_tmp = prod(b==0 ? one(T) - h : h for (h,b)=zip(hobos, bits[1:ind-1]))
        result += result_tmp*hobos[ind]
    end
    result
end

function _hamiltonian_cycle_complete_hobo(vars::Variables)
    n, logn = indexcapacities(vars)
    hobo = HOBO(vars)

    # prohibit numbers
    for t = 1:n
        hobo += smaller_equal_than([vars[t,i] for i=1:logn], n-1)
    end

    # different for different time
    for t1 = 1:n, t2 = (t1+1):n
        hobo += are_different([vars[t1,i] for i=1:logn], [vars[t2,i] for i=1:logn])
    end
    hobo
end

function _hobo_tsp_part(w::Matrix{T}, vars::Variables) where T<:Real
    n = size(w, 1)
    logn = ceil(Int, log(2, n))

    @assert size(w, 2) == n
    hobo = HOBO(vars)

    # prohibit numbers
    for t=1:n, u=1:n, v=1:n
        if u != v
            hobo_tmp = are_same([vars[t,i] for i=1:logn], u-1)
            hobo_tmp *= are_same([vars[t%n+1,i] for i=1:logn], v-1)
            hobo += w[u,v]*hobo_tmp
        end
    end
    hobo
end

function tsp_hobo(w::Matrix{T}, vars::Variables; B::Real=1., A::Real=2*B*maximum(w)) where T<:Real
    n = size(w, 1)
    logn = ceil(Int, log(2, n))

    @assert size(w, 2) == n
    @assert indexcapacities(vars) == [n, logn]
    @assert A >= zero(T) && B >= zero(T)

    A*_hamiltonian_cycle_complete_hobo(vars) + B*_hobo_tsp_part(w, vars)
end

function tsp_hobo(w::Matrix; name::String="TSP", B::Real=1., A::Real=2*B*maximum(w))
    n = size(w, 1)
    tsp_hobo(w, Variables(name, [n, ceil(Int, log(2, n))], ["timestep", "city_binary"]), B=B, A=A)
end


text:
# functions definition

function continued_fraction(x, terms)
    val = x
    for i = 1:terms
        if i==terms
            val = 1/(1-val)
        else
            val = x/(1-val)
        end
    end
    return val
   end

function continued_exp(x, terms)
    val = 1
    for i = 1:terms
        val = exp(val*x)
    end
    return val
end

# continued fraction call section

N = 200
conv = 100
df = 0.01
s = zeros(N, N)
pt = zeros(Complex{Float64},conv)
for x = 1:N , y = 1:N
  dx = df * x - N * df/2
  dy = df * y - N * df/2
  for i = 1:conv
    pt[i] = continued_fraction(dx+dy*im,i)
  end
  s[x,y] = length(Set(pt))
end

using ImageView
using Images
ImageView.imshow(s)

# continued exponential call section

N = 2000
conv = 100
df = 0.01
s = zeros(N, N)
pt = zeros(Complex{Float64},conv)
for x = 1:N , y = 1:N
  dx = df * x - N * df/2
  dy = df * y - N * df/2
  for i = 1:conv
    pt[i] = continued_exp(dx+dy*im,i)
  end
  s[x,y] = length(Set(pt))
end

using ImageView
using Images
ImageView.imshow(s)


text:
struct Margins
    # critical frequency
    ω_c::Float64
    # gain cross-over frequency
    ω_g::Float64
    # gain margin
    gain_margin::Float64
    # phase margin (radians)
    phase_margin::Float64
end

"""
    margins = gain_phase_margins(g_ol, ω_c_guess=0.001, ω_g_guess=0.001)

compute critical frequency (radians / time), gain crossover frequency (radians / time),
gain margin, and phase margin (radians) of a closed loop, given its
closed loop transfer function `g_ol::TransferFunction`.

if ω_c or ω_g is not found (i.e. if either are `NaN`), but the `bode_plot` clearly shows
a critical/gain crossover frequency, adjust `ω_c_guess` or `ω_g_guess` to find the root.

# Example
```
g_ol = 2 * exp(-s) / (5 * s + 1)
margins = gain_phase_margins(g_ol)
margins.ω_c # critical freq. (radians / time)
margins.ω_g # gain crossover freq. (radians / time)
margins.gain_margin # gain margin
margins.phase_margin # phase margin (radians)
```
"""
function gain_phase_margins(g_ol::TransferFunction; ω_c_guess::Float64=0.001, ω_g_guess::Float64=0.001)
    # critical frequency ∠ G(i ω_c) = -π
    ω_c = NaN
    try
        ω_c = fzero(ω -> angle(evaluate(g_ol, im * ω)) + π, ω_c_guess, atol=2*sqrt(eps()))
    catch da_error
        if isa(da_error, Roots.ConvergenceFailed)
            ω_c = NaN
        else
            println(da_error)
            error("something went wrong when computing ω_c")
        end
    end
    # check for spurious solution (often ω_c d.n.e.) TODO check issue on Roots.jl
    if ! isapprox(angle(evaluate(g_ol, im * ω_c)), -π, rtol=0.01)
        ω_c = NaN
    end

    # gain cross over frequency | G(i ω_g) | = 1
    ω_g = NaN
    try
        ω_g = fzero(ω -> abs(evaluate(g_ol, im * ω)) - 1.0, ω_g_guess, atol=2*sqrt(eps()))
    catch da_error
        if isa(da_error, Roots.ConvergenceFailed)
            ω_g = NaN
        else
            println(da_error)
            error("something went wrong when computing ω_g")
        end
    end
    # check for spurious solution (often ω_c d.n.e.) TODO check issue on Roots.jl
    if ! isapprox(abs(evaluate(g_ol, im * ω_g)), 1.0, rtol=0.01)
        ω_g = NaN
    end
    # gain margin = 1 / | G(i ω_c) |
    gm = 1 / abs(evaluate(g_ol, im * ω_c))
    # phase margin = ∠ G(i ω_g) + π
    pm = π + angle(evaluate(g_ol, im * ω_g))
    return Margins(ω_c, ω_g, gm, pm)
end

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
using SwJ1644
using CSV

#df = ReadData("../lc.csv")
df = CSV.read("../lc.csv")
#@show data

#plotlc(data)
using DataFrames
#using Test

df2 = sort!(df)
#@show plot(df2[1],df2[4])
#A1 = convert(Array,df2)
t = convert(Array{Float64},df2[1])
y = convert(Array{Float64},df2[4])
#Peaks = peaks(t,y) # An easy way to find the peaks
Peaks,Dips = peaks2(t,y) # A more proper way.

N = Int64(length(Peaks)/ndims(Peaks))
tPeaks1 = Peaks[:,1]
tPeaks2 = ones(N-1) #initialize an array
DeltatPeak = ones(N-1) #initialize an array
for i=1:N-1
    tPeaks2[i] = (tPeaks1[i+1]+tPeaks1[i])/2
    DeltatPeak[i] = tPeaks1[i+1]-tPeaks1[i]
end
println("The observed number of peaks are: ", N-1)

### To get the Period changing theoretically

function eptFun(ep0::Float64,t0::Float64,t::Float64,nTidal::Float64)
    ep1 = ep0*(t/t0)^nTidal
end

#function ept(t::Float64)
#    ep2 = eptFun(ep0,t0,t,nTidal)
#end #Trying to use this to reduce the number of dumbs. New version of Julia forbids it?

function calcul()
    # parameters
    DL = 3.7e27 # Luminosity distance in unit of cm
    c = 2.9979e10
    eta = 10 # The efficiency converting mass to X-rays
    mStar = 5*2e33 # Mass of the star
    mm = mStar
    ep0 = 1e-3 # ep0: a small value, short for 'epsilon', a part being accreted
    t0 = 1.0 # Right now, the ti is just the number i
    nTidal = -0.2
    MBH = 7e6 # Mass of the central BH, in unit of M_sun
    Rp = 1*1e14 # Rp: the distance of pericenter to the BH, in unit of cm
    Rp = Rp/1.49597871E13 # convert to AU
    N = 4*Int64(floor(1/ep0)) # the integer part of a float
    dm = ones(N) # Initialize the mass
    e = zeros(N)
    e[1] = 0.6
    G = 6.67242E-8
    P = zeros(N)
    P[1] = 2*pi/sqrt(G*MBH) * (Rp/(1-e[1]))^1.5
    N2 = 0
    for i in 2:N
        ep = eptFun(ep0,t0,Float64(i),nTidal) #ept(Float64(i))
        tmp1 = sqrt(1+e[i-1])-ep
        e[i] = (tmp1/(1-ep))^2-1.0
        dm[i] = ep * mm
        mm = mm - dm[i] # why global? weird. Ans: When the whole block is inside a function, you don't need.
        if e[i] >= 1.0
            println("This is the last orbit! eccentricity:", e[i], " mass:", mm)
            break
        end
        P[i] = 2*pi/sqrt(G*MBH) * (Rp/(1-e[i]))^1.5
        #println(i, " tmp1: ", tmp1, " e[i]: ", e[i], " P[i]:", P[i])
        N2 = N2 + 1
    end
    e2 = e[1:N2-1] #truncate the tails, which was set to be larger
    P2 = P[1:N2-1]
    tp3 = P2[1]
    tp4 = ones(N2-1)
    tp4[1] = tp3
    for i in 2:N2-1
        tp3 = tp3 + P2[i]
        tp4[i] = tp3
    end
    dm2 = dm[1:N2-1]
    dm2[1] = dm2[2] # tmp use
    #Fx = 1e6 * eta .* dm2 * c^2 / (4*pi*DL^2) ./ P2 #(time)
    z = 0.35 #redshift
    alpha = 0.6 #spectral index
    Fx = (1+z) * ((1-alpha)*(10/0.3)^(-alpha)/((10/0.3)^(1-alpha)-1)) * 1e23 * eta .* dm2 * c^2 / (4*pi*DL^2) /(1e4*2.418e14) ./ P2 #flux at 10 keV, in Jy
    println("The predicted number of orbits are: ", N2-1)
    return tp4, P2, Fx
end

tp4, P2, Fx = calcul()

#= using Plots
#using Gadfly
#plot(tPeaks2,DeltatPeak)
plot(log10.(tPeaks2),log10.(DeltatPeak))
savefig("fig.png")
=#

# To plot the period-t diagram
using RCall
@rlibrary ggplot2
dat = DataFrame(tP=tPeaks2, DP=DeltatPeak)
dat2 = DataFrame(tp2=tp4, DP2=P2)
ggplot(data=dat, aes(x=log10.(dat[:tP]), y=log10.(dat[:DP]))) +
    geom_point(size=2) +
    geom_line(data=dat2, aes(x=log10.(dat2.tp2),y=log10.(dat2.DP2))) +
    geom_point(data=dat2, aes(x=log10.(dat2.tp2),y=log10.(dat2.DP2)), size=1, color="red") +
    xlab("log10(time) (s)") + ylab("log10(Period) (s)") +
    theme_bw()
#print(p)
ggsave(file="Period-t.pdf")

# To plot the light curves, together with the peaks, dips and predictions.
tt = DataFrame(t=t, y=y)
dfPeaks = DataFrame(t=Peaks[:,1], p=Peaks[:,2])
dfDips = DataFrame(t=Dips[:,1], d=Dips[:,2])
dfPredicts = DataFrame(tp=tp4)
dfFx = DataFrame(t=tp4, Fx=Fx)
p = ggplot(data=tt, aes(x=tt[:t],y=tt[:y]))+
    geom_point(size=0.1, color="grey")+
    #geom_vline(xintercept=dfPredicts[:tp], size=0.1, color="grey")+
    geom_point(data=dfPeaks, aes(x=dfPeaks[:t],y=dfPeaks[:p]), size=0.5, color="red")+
    #geom_point(data=dfDips, aes(x=dfDips[:t],y=dfDips[:d]), size=0.5, color="blue")+
    geom_point(data=dfFx, aes(x=dfFx[:t],y=dfFx[:Fx]), size=0.2, color="orange")+
    scale_x_log10() + scale_y_log10()+
    xlab("log10(time) (s)") + ylab("log10(Flux) (erg/cm^2/s)") +
    theme_bw()
#p = p + theme(legend.title=element_blank())
#print(p)
ggsave(file="light-curves.pdf")

# Just plot part of the light curves to have a better look
##p = p + xlim(1.5e6,1.6e6)
#p = p + xlim(1.e3,5e3)
#ggsave(file="light-curves-zoom-in.pdf")


text:
# Functions for data collection

# mymean(a) = isempty(a) ? missing : mean(a)
mymean(a) = isempty(a) ? 0.0 : mean(a)
mystd(a) = isempty(a) ? 0.0 : std(a)
mymedian(a) = isempty(a) ? 0.0 : median(a)
mysum(a) = isempty(a) ? 0.0 : sum(a)

"""
knowledge output of the researcher given the problems it has worked on and their publication success.
"""
function mean_knowledge_output(researcher::Researcher)
	produced_knowledge = zeros(Float64, length(researcher.problem_history))
	published_indices = findall(a-> !ismissing(a) && a == true, researcher.publication_success)
	for index in published_indices
		problem = researcher.problem_history[index]
		produced_knowledge[index] = problem.information
	end
	return mymean(produced_knowledge)
end

submissions_per_author(r::Researcher) = length(r.publication_success)

publications_per_author(r::Researcher) = length(findall(a-> !ismissing(a) && a==true, r.publication_success))

isgranted(r::Researcher) = r.grant > 0 ? true : false

function totalcitations(r::Researcher)
	j = sum(r.publication_citations)
	if !isnan(j)
		return 0
	else
		return j
	end
end

npapers(r::Researcher) = length(findall(a-> !ismissing(a) && a==true, r.publication_success))

nprojects(r::Researcher) = length(r.publication_success)

pop(r::Researcher) = 1

adata = [
	(mean_knowledge_output, mysum),
	(mean_knowledge_output, mymean),
	# (mean_knowledge_output, mystd),
	(mean_knowledge_output, mymedian),
	(submissions_per_author, mymean),
	# (submissions_per_author, mystd),
	(submissions_per_author, mymedian),
	(publications_per_author, mymean),
	# (publications_per_author, mystd),
	(publications_per_author, mymedian),
	(:experience, mymean),
	(:experience, mymedian),
	# (:experience, mystd),
	(isgranted, count),
	(:grant, mymean),
	(totalcitations, mymean),
	# (totalcitations, mystd),
	(totalcitations, mymedian),
	(:risk_taking, mymean),
	# (:risk_taking, mystd),
	(:risk_taking, mymedian),
	(:generality, mymean),
	# (:generality, mystd),
	(:generality, mymedian),
	(npapers, mymean),
	# (npapers, mystd),
	(npapers, mymedian),
	(nprojects, mymean),
	(nprojects, mymedian),
	# (nprojects, mystd),
	(pop, mysum)
]


text:
using InteractBulma, InteractBase
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "ijulia" begin
    @test !InteractBase.isijulia()
end

@testset "deps" begin
    libs = InteractBase.libraries(Bulma())
    @test all(isfile, libs)
end


text:
struct Load <: AbstractPipeline
    table::Observable{SimpleTable}
    files::Observable{Union{Nothing, Dict}}
    value::Observable{SimpleTable}
end

function get_vertices(l::Load)
    names = collect(Tables.columnnames(l.table[]))
    return [Vertex(:Load, Symbol[], names)]
end

get_vertex_names(::Load) = [:Load]

function Load(table::Observable{SimpleTable})
    files = Observable{Union{Nothing, Dict}}(nothing)
    value = Observable(table[])
    return Load(table, files, value)
end

function jsrender(session::Session, l::Load)

    ui = scrollable_component(FilePicker(l.files))

    tryon(session, l.files) do files
        !isempty(files) && (l.table[] = readfiles(files))
    end

    tryon(session, l.table) do table
        l.value[] = table
    end

    return jsrender(session, ui)
end


text:
@with_kw struct Param

    #Resume previously run sumulation?
    resume_simulation   ::Bool      = :false
    input_file_name     ::String    = ""

    #Temporal variables for simulation
    t_sim_init          :: Float64   = 0.0
    t_sim_final         :: Float64   = 0.0
    out_every_t         :: Float64   = 1.0
    deltat              :: Float64   = 0.05

    #Max possible sim
    max_runtime         ::String    = "70:00:00"

    #Spacial variables
    xmin            :: Float64      = -30.0
    xmax            :: Float64      =  30.0
    xnodes          :: Integer      =  1024

    ymin            :: Float64      = -30.0
    ymax            :: Float64      =  30.0
    ynodes          :: Integer      =  1024

    #Initial Conditions
    A0              ::Float64       = 1.0
    σ               ::Float64       = 2.0
    r0              ::Float64       = 5.0
    ω               ::Float64       = 5.0
    m               ::Integer       = 0.0

    #Klein Gordon Mass
    μ               ::Float64       = 1.0

    #Dissipation terms
    dissipation     ::Bool          =  :false

    R_orbit         ::Float64       = 10
    R1              ::Float64       = 1
    R2              ::Float64       = 1

    Ω               ::Float64       = 0.2
    alpha1          ::Float64       = 1.0
    alpha2          ::Float64       = 1.0

    #Periodic BC?
    Boundaries      ::Symbol        = :radial # :square :periodic

    #Data variables
    folder          :: String    = "data/"
    fname           :: String    = "default_file_name"

end


text:
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("FVheat1d");

useLog("FVheat1dlog")

# Configuration setup
domain(1)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 40 # number of elements
mesh(LINEMESH, elsperdim=n)

# Variables and BCs
u = variable("u", SCALAR, CELL)
boundary(u, 1, FLUX, "0")

# Time interval and initial condition
T = 0.1;
timeInterval(T)
initial(u, "x<0.5 ? sin(pi*x)^4 : 0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
coefficient("D", 0.1) # Diffusion rate
fluxAndSource(u, "-D*dot(grad(u),normal())", "0")

#exportCode("fvheat1dcode") # uncomment to export generated code to a file
#importCode("fvheat1dcode") # uncomment to import code from a file

solve(u)

@finalize()

##### Uncomment below to plot
# # The initial condition
# u0 = zeros(n);
# x = Femshop.fv_info.cellCenters[:]
# for i=1:n
#     u0[i] = x[i]<0.5 ? sin(pi*x[i])^4 : 0
# end

# using Plots
# pyplot();
# display(plot([x x], [u0 u.values[:]], markershape=:circle, label=["initial" "t="*string(T)]))


text:
# ==============================================================================
# Modele JuMP pour calculer la relaxation linéaire du 2SPA sur un objectif donne, activant eventuellement une ϵ-contrainte

function computeLinearRelax2SPA(  nbvar::Int,
                                  nbctr::Int,
                                  A::Array{Int,2},
                                  c1::Array{Int,1},
                                  c2::Array{Int,1},
                                  epsilon,
                                  obj::Int
)
  model = Model(GLPK.Optimizer)
  @variable(model, 0.0 <= x[1:nbvar] <= 1.0 )
  @constraint(model, [i=1:nbctr],(sum((x[j]*A[i,j]) for j in 1:nbvar)) == 1)
  if obj == 1
    @objective(model, Min, sum((c1[i])*x[i] for i in 1:nbvar))
    @constraint(model, sum((c2[i])*x[i] for i in 1:nbvar) <= epsilon)
  else
    @objective(model, Min, sum((c2[i])*x[i] for i in 1:nbvar))
    @constraint(model, sum((c1[i])*x[i] for i in 1:nbvar) <= epsilon)
  end
  optimize!(model)
  return objective_value(model), value.(x)
end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
function formpanel(;
  panelname::String = PARAM[:panelname],
  savepath::String = PARAM[:workingpath],
  testdate::Date = PARAM[:testdate],
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM,
  binextension::String = BIN_EXTENSION,
  includecomp::Bool = PARAM[:compincludecomp]
  )::DataFrame

  #housekeeping
  local panel::DataFrame
  local crsp::DataFrame
  local panelpath::String = "$savepath\\$panelname.$binextension"

  if PARAM[:crspfrequencysuffix] == "w"
    @info "loading weekly data"
  elseif PARAM[:crspfrequencysuffix] == "m"
    @warn "Using monthly data"
  else
    @assert false
  end

  #The below block directs the panel creation
  if PARAM[:refreshpanel]

    #quick way to filter the dates
    if PARAM[:crspfrequencysuffix] == "w"
      crsp = makeweekly(selectsummarydata)
    elseif PARAM[:crspfrequencysuffix] == "m"
      crsp = makemonthly(selectsummarydata)
    end

    if includecomp #allows us to short circuit the merge if not using comp data
      comp = prepcomp()
      panel = mergecrspcomp(;crsp, comp)
    else
      panel = crsp
    end
    @assert includecomp == (:gvkey ∈ propertynames(panel))
    @info "$(includecomp ? "comp rows included! unmerged rows dropped" :
      "comp excluded, so rows that would have been unmerged were kept.")"

    panel = processpanel(panel)

    if PARAM[:testpanel]

      #time series test output
      testpermno::Int = PARAM[:testpermnomult]
      testview::SubDataFrame = view(panel, (panel.permno .% testpermno) .== 0, :)
      if size(testview,1)>20_000 #prevents scenarios where the file is larger than needed
        testview=view(testview, 1:20_000,:)
      end
      testview |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_panel.csv")

      #cross-sectional test output
      testview = view(panel, panel.date .== testdate,:)
      testview |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_panel_$testdate.csv")

    end

    @info "Panel formed of size $(size(panel))."

    outbinstream(panelpath, panel)
  else

    panel = inbinstream(panelpath)
  end

  return panel
end

function boundpanelbydate!(panel::DataFrame;
  lowerdatebound::NDate=PARAM[:lowerdatebound],
  upperdatebound::NDate=PARAM[:upperdatebound])
  if !(lowerdatebound === nothing)
    @info "filtering panel to date≥$(lowerdatebound) (current lowest: $(minimum(panel.date)))"
    filter!(:date=>d->d≥lowerdatebound, panel)
  end
  if !(upperdatebound === nothing)
    @info "filtering panel to date≤$(upperdatebound) (current highest: $(maximum(panel.date)))"
    filter!(:date=>d->d≤upperdatebound, panel)
  end

end

function placebostrategies!(panel::DataFrame)
  N=size(panel,1)
  panel.placebo = rand(N) |> Vector{MFloat64}
  panel.placebofixed = Vector{MFloat64}(undef, N)
  panel.placebofixeddiscrete = Vector{MFloat64}(undef, N)

  discretevector::Vector{Float64} = collect(0.5:0.5:5.0)
  for spanel ∈ groupby(panel, :permno)
    spanel.placebofixed .= rand()
    spanel.placebofixeddiscrete .= rand(discretevector)
  end

  return nothing
end

#iterate through the momentums trategies and put them into the dataframe
function fundamentalstrategies!(panel::DataFrame;
    strategies=PARAM[:fundamentalstrategies],
    bookequityfloor::Float64 = PARAM[:compbookequityfloor])

    #book to market
    panel.bm = finiteormissing.(panel.bkequity .* 1_000_000 ./ panel.mc)

    ######testing code for comparison to Lewellen
    assetfloor(x) = max(x,bookequityfloor)#max(x,bookequityfloor)
    assetfloor(::Missing) = missing
    panel.fyear = year.(panel.fdate)
    transform!(groupby(panel, [:permno, :fyear]), [:Lmc, :bkequity] |> AsTable =>
      ((t)->log(assetfloor(t.bkequity[1]) * 1_000_000 / t.Lmc[1])) =>
      :lbLmquarterly)
    @assert length(panel.lbLmquarterly) == size(panel,1)
    panel.lbLm = log.(panel.bkequity .* 1_000_000 ./ panel.Lmc) .|> finiteormissing
    #######

    winsorizewithin!(panel, Fsource=:ret, prop=PARAM[:qwinsorize], Fgroup=:date,
      twosided=true, Fnew=:WXret)

    for Fstrategy ∈ strategies
      FWstrategy::Symbol = Symbol(:W, Fstrategy)
      FWXstrategy::Symbol = Symbol(:WX, Fstrategy)
      panel[!,FWstrategy] = winsorizequantile(panel[!, Fstrategy], PARAM[:qwinsorize], twosided=true)

      winsorizewithin!(panel, Fsource=Fstrategy, prop=PARAM[:qwinsorize], Fgroup=:date,
        twosided=true, Fnew=FWXstrategy)
      lagwithin2!(panel, [Fstrategy, FWstrategy, FWXstrategy], :permno, date=:date) #this could be parallelized

      #NOTE- fairly sure we don't need the below
      #FWiivolWstrategy = Symbol(:Wiivol, FWstrategy)
      #panel[!,Symbol(:Wiivol, FWstrategy)] =  panel[!, :Wiivol] .* panel[!, FWstrategy]
      #panel[!,Symbol(:LWiivol, FWstrategy)] =  panel[!, :LWiivol] .* panel[!, Symbol(:L, FWstrategy)]
    end

end

function processpanel(panel; momentumstrategies=PARAM[:momentumstrategies],
    includecomp::Bool = PARAM[:compincludecomp])
  #winsorize and compute the nomalized volume
  panel.Wdvol = (winsorizequantile(
    panel.dvol, PARAM[:qwinsorize], twosided=true) ./ PARAM[:volscale])

  panel.dvolscaled = panel.dvol ./ PARAM[:volscale]
  differencewithin2!(panel, [:Wdvol], :permno, date=:date)
  #computenvol!(panel, :Wdvol)

  #this forms a market cap control and the basis of the size characteristic
  panel.lmc = log.(panel.mc) .|> finiteormissing
  panel.mcscaled = panel.mc ./ PARAM[:volscale]
  lagwithin2!(panel, [:mc, :lmc, :mcscaled], :permno, date=:date)
  lagwithin2!(panel, [:Lmc, :Llmc, :Lmcscaled], :permno, date=:date)
  #transform!(groupby(panel, :permno),
  #  )
  transform!(groupby(panel, :permno),
    :lmc => ((lmc)->mean(skipmissing(lmc))) => :fixedlmc,
    :mc => ((mc)->mean(skipmissing(mc))) => :fixedmc)

  panel.vecof1s = ones(size(panel,1)) #used for fixed (constant)

  #winsorize idiosyncratic volatiltiy
  if (panel.ivol .!== missing ) |> all
    panel.Wivol = winsorizequantile(panel.ivol, PARAM[:qwinsorize], twosided=true)
    panel.Wiivol = 1.0 ./ panel.Wivol
    lagwithin2!(panel, [:Wivol, :Wiivol], :permno, date=:date)
  else
    select!(panel, Not(:ivol))
  end

  #form the momentum characteristic
  #panel=panel[panel.ndays .≥ 10, :] #WARNING- comment this
  if PARAM[:momentumstrategyfrequency] === :aggregate
    select!(panel, Not([:plret; keys(PARAM[:momentumstrategies])|>collect]))
    panel.lret = (log).(panel.ret .+ 1)
    panel.plret = adjustindustry(panel, Fmc=:mc)
    panel.lplret = (log).(panel.plret .+ 1)
    momentumstrategies!(panel)
  elseif PARAM[:momentumstrategyfrequency] === :day
    @assert setdiff([:plret; keys(PARAM[:momentumstrategies])|>collect], propertynames(panel)) |> isempty
  elseif PARAM[:momentumstrategyfrequency] === :test
    rename!(panel, Dict(s=>Symbol(:d_,s) for s ∈ [:plret; keys(PARAM[:momentumstrategies])|>collect]))
    panel.lret = (log).(panel.ret .+ 1)
    panel.plret = adjustindustry(panel, Fmc=:mc)
    panel.lplret = (log).(panel.plret .+ 1)
    momentumstrategies!(panel)
    throw("aggregate momentum strategies computed")
  else
    throw("unrecognized momentum strategy frequency (:day, :aggregate) are acceptable")
  end

  #use this for within-winsorization
  panel.yearmonth = Symbol.(:Y, year.(panel.date), :M, month.(panel.date))
  #momentumstrategies!(panel)
  conditionmomentumstrategies!(panel)

  if includecomp
    fundamentalstrategies!(panel)
  end

  #form characteristics normalized by market value

  #form market value characteristic
  panel.mcthorn = finiteormissing.(panel.mc) ./ 10.0^9
  lagwithin2!(panel, [:mcthorn], :permno, date=:date)
  lagwithin2!(panel, [:Lmcthorn], :permno, date=:date)

  boundpanelbydate!(panel) #restrict the dates
  placebostrategies!(panel) #create placebo characteristics

  return panel
end

#=function ensurecompleteness(crsp::DataFrame; Ftimegroup::Symbol,
  datefrequency::Symbol)

  #drop all time groups (month or week) with insufficient data
  #use a loose absolute criteria and a stricter criteria based on a fraction of the median
  mindaysmedianfrac::Float64 = PARAM[:mindaysmedianfrac][datefrequency]
  absolutemindays::Int = PARAM[:absolutemindays][datefrequency]
  crsp.tokeep = trues(size(crsp,1))
  scrsps::GroupedDataFrame = groupby(crsp, Ftimegroup)
  Threads.@threads for i ∈ 1:length(scrsps)
    scrsp::SubDataFrame = scrsps[i]
    #now employ some heuristics to make sure we don't have too much missing data
    medianndays::Int = median(scrsp.ndays)
    #minndays::Int = medianndays ≥ 2 ? medianndays - 1 : medianndays
    minndays::Int = max(ceil(medianndays*mindaysmedianfrac) |> Int, absolutemindays)
    scrsp.tokeep .= scrsp.ndays .≥ minndays
  end

  crsp = crsp[crsp.tokeep, :]
  select!(crsp, Not(:tokeep))

  return crsp
end=#

function selectsummarydata(crsp::DataFrame;
  #NOTE: Probably ok for now, but consider doing these in quantiles by cross-section instead
  minprice::Float64 = PARAM[:minprice],
  minmarketcap::Float64 = PARAM[:minmarketcap],
  minnysequantile::Float64 = PARAM[:minnysequantile],
  Ftimegroup::Symbol,
  datefrequency::Symbol)

  #@assert all((!ismissing).(crsp.mc))
  #@assert all((!ismissing).(crsp.exchcd))
  crsp.mc = Vector{Float64}(crsp.mc)
  crsp.exchcd = Vector{Int}(crsp.exchcd)

  function assignnysequantile(t)
    nysemc = t.mc[t.exchcd .== 1]
    nysecdf = ecdf(nysemc)
    return nysecdf.(t.mc)
  end

  @assert issorted(crsp, [:permno, :date])
  @assert issorted(crsp, [:permno, Ftimegroup])
  filtercols = [:permno, Ftimegroup, :mc, :exchcd, :price]
  @assert setdiff(filtercols, propertynames(crsp)) |> isempty
  crspeop::DataFrame = combine(groupby(crsp[!,filtercols], [:permno, Ftimegroup]),
    filtercols .=> last,
    nrow => :ndays, renamecols=false)
  #debugmsg("dbg2.01", crspeop)
  @assert issorted(crspeop, [:permno, Ftimegroup])
  transform!(groupby(crspeop, Ftimegroup), [:mc, :exchcd] |> AsTable =>
    assignnysequantile => :nysequantile)
  #debugmsg("dbg2.02", crspeop)
  nysebreakpoint(t) = minimum(t.mc[t.nysequantile .≥ minnysequantile])
  transform!(groupby(crspeop, Ftimegroup), [:mc, :nysequantile] |> AsTable =>
    nysebreakpoint => :breakpoints)
  #debugmsg("dbg2.03", crspeop)
  #compute quantiles based on nyse thresholds
  averagebreakpoint::Float64 = mean(combine(groupby(crspeop, Ftimegroup), :breakpoints =>
    last, renamecols=false).breakpoints)

  @info "Average NYSE breakpoint: $averagebreakpoint"

  #question- should I filter out firms
  filter!(AsTable([:price, :mc, :nysequantile]) =>
    (r ->
    (r.price ≥ minprice) &&
    (r.mc ≥ minmarketcap) &&
    (r.nysequantile ≥ minnysequantile)), crspeop)

  @assert all(crspeop.mc .≥ crspeop.breakpoints)
  crspeop = crspeop[!, [:permno, Ftimegroup, :ndays, :nysequantile]]
  #debugmsg("dbg2.04", crspeop)
  #drop all time groups (month or week) with insufficient data
  #use a loose absolute criteria and a stricter criteria based on a fraction of the median
  #mindaysmedianfrac::Float64 = PARAM[:mindaysmedianfrac][datefrequency]
  absolutemindays::Int = PARAM[:absolutemindays][datefrequency]

  medianndays(ndays) = max(floor(median(ndays)) - 1 |> Int,1)
  transform!(groupby(crspeop, Ftimegroup), :ndays => medianndays => :medianndays)
  crspeop = crspeop[(crspeop.ndays .≥ crspeop.medianndays) .& (crspeop.medianndays .≥ absolutemindays), :]
  select!(crspeop, Not([:medianndays]))

  #debugmsg("dbg2.05", crspeop)

  return crspeop
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using DEBmicroTrait
using CSV, DataFrames, Statistics
using GLM

################################################
# I/O
df = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/dethlefsen_schmidt_2008.csv", DataFrame, missingstring="N/A")
gdf = groupby(df, :Species_Name)
df = combine(gdf, :Protein_Quantity_Measurement => mean, :Size_Measurement => mean, :DNA_Quantity_Measurement => mean, :RNA_Quantity_Measurement => mean, :Specific_Growth_Rate_Measurement => mean, )
df_out = DataFrame()
################################################

################################################
# Dethlefsen & Schmidt (2008): translational_power = μ*P/R
# Protein data
d_p = 1.37e6
P = df.Protein_Quantity_Measurement_mean
V_P = P/d_p
df_out.V_P = V_P
# RNA data
d_r = 1.79e6
R = df.RNA_Quantity_Measurement_mean
V_R = R/d_r
df_out.V_R = V_R
# Growth rate data
gmax = df.Specific_Growth_Rate_Measurement_mean
df_out.gmax = gmax
# Cell size data
V_cell = df.Size_Measurement_mean
df_out.V_cell = V_cell
# k_E calculation
k_E = gmax.*V_P./V_R
df_out.k_E = k_E
################################################

################################################
#Kempes et al. (2016): allometric predictions
# Cell size prediction
d_DNA = 2e6
DNA = df.DNA_Quantity_Measurement_mean
V_DNA = DNA/d_DNA
v_N = 1.47e-27
L_DNA = V_DNA/v_N
V_cell_model = DEBmicroTrait.genome_size_to_cell_volume(L_DNA)
df_out.V_cell_model = V_cell_model
# Protein prediction
V_P_model = DEBmicroTrait.cell_volume_to_protein_volume(V_cell)
df_out.V_P = V_P
df_out.V_P_model = V_P_model
# RNA prediction
V_R_model = DEBmicroTrait.cell_volume_to_ribosome_volume(V_cell, gmax)
df_out.V_R = V_R
df_out.V_R_model = V_R_model
# tRNA predictions
V_tRNA_model = DEBmicroTrait.cell_volume_to_tRNA_volume(V_cell, gmax)
df_out.tRNA_model = V_tRNA_model
# mRNA predictions
V_mRNA_model = DEBmicroTrait.cell_volume_to_mRNA_volume(V_cell, gmax)
df_out.mRNA_model = V_mRNA_model
# total RNA
V_mtR_model = V_R_model .+ V_tRNA_model
df_out.V_mtR_model = V_mtR_model
# k_E prediction
k_E_model= DEBmicroTrait.translation_power(V_P_model, V_R_model, log(2)./gmax)
df_out.k_E_model = k_E_model
################################################

################################################
# Statistics
linearRegressor = lm(@formula(log(k_E) ~ log(k_E_model)), df_out)
r2(linearRegressor)
# test for significant difference to 1-1 line
df_out.diff = k_E_model-k_E
linearRegressor = lm(@formula(diff ~ k_E_model), df_out)
# RMSD
RMSD = sqrt((1/(12-1)*sum((k_E_model - k_E).^2)))
################################################

################################################
# I/O
CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/manuscript/files/benchmark_protein_synthesis.csv", df_out)
################################################

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using Test
using FWFTables
using FixedSizeStrings

specs = [
         Varspec("var1", FixedSizeString, 1, 2, 0),
         Varspec("var2", FixedSizeString, 3, 2, 0),
        ]
tb = FWFTables.File(IOBuffer("0102\n0304\n0506\n"), specs)

@test length(tb) == 3
@test tb.var1[1] == FixedSizeString{2}("01")
@test tb.var1[2] == FixedSizeString{2}("03")
@test tb.var1[3] == FixedSizeString{2}("05")
@test tb.var2[1] == FixedSizeString{2}("02")
@test tb.var2[2] == FixedSizeString{2}("04")
@test tb.var2[3] == FixedSizeString{2}("06")

# no newline at end
tb = FWFTables.File(IOBuffer("0102\n0304\n0506"), specs)
@test length(tb) == 3

# empty input
tb = FWFTables.File(IOBuffer(""), specs)
@test length(tb) == 0

# empty input
tb = FWFTables.File(IOBuffer("\n"), specs)
@test length(tb) == 0

max_stars_repo_name: HenricoWitvliet/FWFTables


text:
module TestsStackBuilding

using Test
using TMLEEpistasis
using MLJBase
using TOML
using TMLE
using MLJGLMInterface
using MLJLinearModels
using MLJModels

include("helper_fns.jl")

@testset "Categorical target TMLE built from configuration file" begin
    tmle_config = joinpath("config", "tmle_config.toml")
    build_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    tmle_bin = TMLEEpistasis.estimator_from_toml(TOML.parsefile(tmle_config), queries, Bool, adaptive_cv=false)

    @test tmle_bin.threshold == 0.001
    # Test binary target TMLE's Qstack
    @test tmle_bin.Q̅.measures == [log_loss]
    @test tmle_bin.F isa LinearBinaryClassifier
    ## Checking Qstack.metalearner
    @test tmle_bin.Q̅.metalearner isa LogisticClassifier
    @test tmle_bin.Q̅.metalearner.fit_intercept == false
    ## Checking Qstack.resampling
    @test tmle_bin.Q̅.resampling isa StratifiedCV
    @test tmle_bin.Q̅.resampling.nfolds == 2
    ## Checking Qstack EvoTree models
    @test tmle_bin.Q̅.EvoTreeClassifier_1.nrounds == 10
    ## Checking Qstack  Interaction Logistic models
    @test tmle_bin.Q̅.InteractionLMClassifier_1 isa TMLEEpistasis.InteractionLMClassifier
    @test tmle_bin.Q̅.InteractionLMClassifier_1.interaction_transformer.column_pattern == r"^RS_"
    ## Checking Qstack HAL model
    @test tmle_bin.Q̅.HALClassifier_1.lambda == 10
    @test tmle_bin.Q̅.HALClassifier_1.smoothness_orders == 1
    @test tmle_bin.Q̅.HALClassifier_1.cv_select == false
    @test tmle_bin.Q̅.HALClassifier_1.num_knots == [10, 5]

    # Test binary target TMLE Qstack
    tmle_cont = TMLEEpistasis.estimator_from_toml(TOML.parsefile(tmle_config), queries, Real, adaptive_cv=false)
    @test tmle_cont.Q̅.measures == [rmse]
    @test tmle_cont.F isa MLJGLMInterface.LinearRegressor
    ## Checking Qstack.metalearner
    @test tmle_cont.Q̅.metalearner isa MLJLinearModels.LinearRegressor
    @test tmle_cont.Q̅.metalearner.fit_intercept == false

    ## Checking Qstack.resampling
    @test tmle_cont.Q̅.resampling isa CV
    @test tmle_cont.Q̅.resampling.nfolds == 2
    ## Checking Qstack EvoTree models
    @test tmle_cont.Q̅.EvoTreeRegressor_1.nrounds == 10
    @test tmle_cont.Q̅.EvoTreeRegressor_2.nrounds == 20
    ## Checking Qstack Interaction Linear model
    @test tmle_cont.Q̅.InteractionLMRegressor_1.interaction_transformer.column_pattern == r"^RS_"
    ## Checking Qstack HAL model
    @test tmle_cont.Q̅.HALRegressor_1.lambda == 10
    @test tmle_cont.Q̅.HALRegressor_1.smoothness_orders == 1
    @test tmle_cont.Q̅.HALRegressor_1.cv_select == false
    @test tmle_cont.Q̅.HALRegressor_1.num_knots == [10, 5]

    # Both TMLE have the same G Stack
    expected_queries = [
        Query(case=(RSID_10="AG", RSID_100="AG"), control=(RSID_10="GG", RSID_100="GG"), name="QUERY_1"),
        Query(case=(RSID_10="AG", RSID_100="AA"), control=(RSID_10="GG", RSID_100="GG"), name="QUERY_2"),
        Query(case=(RSID_10="AG", RSID_100="AA"), control=(RSID_10="TT", RSID_100="GG"), name="QUERY_3")
    ]
    for tmle in [tmle_bin, tmle_cont]
        @test tmle.G.model.measures == [log_loss]
        test_queries(tmle.queries, expected_queries)
        # Checking Gstack
        @test tmle.G isa FullCategoricalJoint
        ## Checking Gstack.metalearner
        @test tmle.G.model.metalearner isa LogisticClassifier
        @test tmle.G.model.metalearner.fit_intercept == false
        ## Checking Gstack.resampling
        @test tmle.G.model.resampling isa StratifiedCV
        @test tmle.G.model.resampling.nfolds == 2
        ## Checking Gstack models
        @test tmle.G.model.LogisticClassifier_1.lambda == 1.0
        @test tmle.G.model.EvoTreeClassifier_1.nrounds == 10
    end
    rm(queryfile)
end

@testset "Test standard ATE TMLE build" begin
    tmle_config = joinpath("config", "tmle_config.toml")
    build_ate_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    tmle_bin = TMLEEpistasis.estimator_from_toml(TOML.parsefile(tmle_config), queries, Bool, adaptive_cv=true)
    tmle_cont = TMLEEpistasis.estimator_from_toml(TOML.parsefile(tmle_config), queries, Real, adaptive_cv=true)
    expected_queries = [
        Query(case=(RSID_10="AG",), control=(RSID_10="GG",), name="QUERY_1"),
        Query(case=(RSID_10="AA",), control=(RSID_10="GG",), name="QUERY_2"),
    ]
    for (type, tmle) in [("binary", tmle_bin), ("continuous", tmle_cont)]
        test_queries(tmle.queries, expected_queries)
        # Checking Gstack
        @test tmle.G isa Stack
        @test tmle.G.resampling isa TMLEEpistasis.AdaptiveCV
        @test tmle.G.resampling.cv isa StratifiedCV

        @test tmle.Q̅.resampling isa TMLEEpistasis.AdaptiveCV
        if type =="binary"
            @test tmle.Q̅.resampling.cv isa StratifiedCV
        else
            @test tmle.Q̅.resampling.cv isa CV
        end
    end
    rm(queryfile)
end

end;

true


text:
# To implement joints between two rigid bodies.
# Location of joint specified initially.

# include("RigidBody.jl")
# include("OrientationConversion.jl")

using LinearAlgebra

mutable struct Joint
    RB1::RigidBody
    pos1::Vector{T} where T<:Real # Coordinates of joint in body-fixed frame of body 1

    RB2::RigidBody
    pos2::Vector{T} where T<:Real# Coordinates of joint in body-fixed frame of body 2

    type::String # "Revolute", "Revolute2", "Spherical", "Weld", "Free", "Spring"

    # Revolute Joint
    axis::Vector{T} where T<:Real # Axis for revolute joint

    # Spring-Damper
    k::T where T<:Real
    restLen::T where T<:Real
    c::T where T<:Real

    # Joint Actuation: Defined in the body frame of the parent link
    jointF::Vector{T} where T<:Real # Length: 6 [F;τ]

    function Joint(body1::RigidBody, body2::RigidBody, rj1::Vector{T},
        rj2::Vector{T}; type::String="Free",
        axis::Vector{T}=[0.0;0.0;1.0], k::T=0.0, rL::T=0.0, c::T=0.0,
        jointForce::Vector{T} = zeros(T,3), jointTorque::Vector{T} = zeros(T,3)) where T<:Real
        # Default values for type, axis, k, and restLen provided.
        allowedTypes = ["Revolute", "Revolute2", "Spherical", "Weld", "Free", "Spring"]
        if !in(type,allowedTypes)
            error("Unknown joint specified.")
        end

        this = new()
        this.RB1 = body1;
        this.RB2 = body2
        this.pos1 = rj1 # Position of joint in body frame of body1 wrt body1 cm
        this.pos2 = rj2 # Position of joint in body frame of body2 wrt body2 cm
        this.type = type
        this.axis = axis # Revolute joint axis
        this.k = k # Spring constant
        this.c = c # Damper
        this.restLen = rL # Rest length of spring
        this.jointF = [jointForce;jointTorque]

        return this
    end
end

# Number of constraint equations for all kinds of joints
JointDict = Dict(["Revolute"=>7, "Revolute2"=>6, "Spherical" => 5, "Weld"=>9, "Spring" => 2, "Free" => 2]);
JointDictIn = Dict(["Revolute"=>6, "Revolute2"=>5, "Spherical" => 4, "Weld"=>8 , "Spring" => 1, "Free" => 1]); # When the first body is the inertial frame

function body_in_jointVec(jVec::Vector{Joint}, b_id::Int64)
# returns the joint id which contains b_id in the vector of joints
    # j = 0;
    # j_i = 0;
    RB = InertialFrameAsRB();
    for i=1:length(jVec)
        if (b_id == jVec[i].RB1.bodyID)
            # j_i = i;
            RB = jVec[i].RB1;
            break;
        elseif (b_id ==  jVec[i].RB2.bodyID)
            # j_i = i;
            # j = 2;
            RB = jVec[i].RB2;
            break;
        else continue;
        end
    end
    return RB;
    # j_i is the position of the first joint in the vector that contains corresponding body
    # j is either 1 or 2, i.e., which body of the corresponding joint
end

function filterInertialJoint!(jVec::Vector{Joint})
    # Remove the inertial body/frame occurring twice in a joint tree
    for i=2:length(jVec)
        if (jVec[i].RB1.m == 0 || jVec[i].RB2.m == 0)
            deleteat!(jVec,i);
            break;
        end
        # RB1 = jVec[i].RB1; RB2 = jVec[i].RB2;
        # for j=1:length(jVec)
        #
        # end
    end
end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
# Source code from this PR
# https://github.com/JuliaStats/HypothesisTests.jl/pull/124
using HypothesisTests
using LinearAlgebra
using StatsFuns: norminvcdf, normccdf
export ShapiroWilkTest, SWCoeffs

#=
From:
PATRICK ROYSTON
Approximating the Shapiro-Wilk W-test for non-normality
*Statistics and Computing* (1992) **2**, 117-119
DOI: [10.1007/BF01891203](https://doi.org/10.1007/BF01891203)
=#

# TODO: Rerun simulation and polynomial fitting

const ROYSTON_COEFFS = Dict{String, Vector{Float64}}(
"C1" => [0.0E0, 0.221157E0, -0.147981E0, -0.207119E1, 0.4434685E1, -0.2706056E1],
"C2" => [0.0E0, 0.42981E-1, -0.293762E0, -0.1752461E1, 0.5682633E1, -0.3582633E1],
"C3" => [0.5440E0, -0.39978E0, 0.25054E-1, -0.6714E-3],
"C4" => [0.13822E1, -0.77857E0, 0.62767E-1, -0.20322E-2],
"C5" => [-0.15861E1, -0.31082E0, -0.83751E-1, 0.38915E-2],
"C6" => [-0.4803E0, -0.82676E-1, 0.30302E-2],
"C7" => [0.164E0, 0.533E0],
"C8" => [0.1736E0, 0.315E0],
"C9" => [0.256E0, -0.635E-2],
"G"  => [-0.2273E1, 0.459E0]
)

for (s,c) in ROYSTON_COEFFS
    @eval $(Symbol("_"*s))(x) = Base.Math.@horner(x, $(c...))
end

#=
The following hardcoded constants has been replaced by more precise values:

SQRTH = sqrt(2.0)/2.0 # 0.70711E0
TH = 3/8 # 0.375E0
SMALL = eps(1.0) # 1E-19
PI6 = π/6 # 0.1909859E1
STQR = asin(sqrt(0.75)) # 0.1047198E1
=#

struct SWCoeffs <: AbstractVector{Float64}
    N::Int
    A::Vector{Float64}
end

Base.size(SWc::SWCoeffs) = (SWc.N,)
Base.IndexStyle(::Type{SWCoeffs}) = IndexLinear()

function Base.getindex(SWc::SWCoeffs, i::Int)
    if i <= lastindex(SWc.A)
        return SWc.A[i]
    elseif i <= length(SWc)
        if isodd(SWc.N) && i == div(SWc.N, 2) + 1
            return 0.0
        else
            return -SWc.A[SWc.N+1-i]
        end
    else
        throw(BoundsError(SWc, i))
    end
end

function SWCoeffs(N::Int)
    if N < 3
        throw(ArgumentError("N must be greater than or equal to 3: $N"))
    elseif N == 3 # exact
        return SWCoeffs(N, [sqrt(2.0)/2.0])
    else
        # Weisberg&Bingham 1975 statistic; store only positive half of m:
        # it is (anti-)symmetric; hence '2' factor below
        m = [-norminvcdf((i - 3/8)/(N + 1/4)) for i in 1:div(N,2)]
        mᵀm = 2sum(abs2, m)

        x = 1/sqrt(N)

        a₁ = m[1]/sqrt(mᵀm) + _C1(x) # aₙ = cₙ + (...)

        if N ≤ 5
            ϕ = (mᵀm - 2m[1]^2)/(1 - 2a₁^2)
            m = m/sqrt(ϕ) # A, but reusing m to save allocs
            m[1] = a₁
        else
            a₂ = m[2]/sqrt(mᵀm) + _C2(x) # aₙ₋₁ = cₙ₋₁ + (...)
            ϕ = (mᵀm - 2m[1]^2 - 2m[2]^2)/(1 - 2a₁^2 - 2a₂^2)
            m = m/sqrt(ϕ) # A, but reusing m to save allocs
            m[1], m[2] = a₁, a₂
        end

        return SWCoeffs(N, m)
    end
end

function swstat(X::AbstractArray{T}, A::SWCoeffs) where T<:Real

    if X[end] - X[1] < lastindex(X)*eps(eltype(X))
        throw("Data seems to be constant!")
    end

    AX = dot(A,X)
    S² = sum(abs2, X.-mean(X))

    return AX^2/S²
end

function pvalue(W::T, A::SWCoeffs, N1=A.N) where T<:Real

    if A.N == 3 # exact by Shapiro&Wilk 1965
        return π/6 * (asin(sqrt(W)) - asin(sqrt(0.75)))
    elseif A.N ≤ 11
        γ = _G(A.N)
        if γ ≤ log(1 - W)
            return zero(W)
        end
        w = -log(γ - log(1 - W))
        μ = _C3(A.N)
        σ = exp(_C4(A.N))
    else
        w = log(1 - W)
        μ = _C5(log(A.N))
        σ = exp(_C6(log(A.N)))
    end

    if (A.N - N1) > 0
        throw("Not implemented yet!")
    end

    return normccdf((w - μ)/σ)
end

struct ShapiroWilkTest <: HypothesisTests.HypothesisTest
    SWc::SWCoeffs            # Expectation of order statistics for Shapiro-Wilk test
    W::Float64            # test statistic
    N1::Int               #
end

HypothesisTests.testname(::ShapiroWilkTest) = "Shapiro-Wilk normality test"
HypothesisTests.population_param_of_interest(t::ShapiroWilkTest) =
("Squared correlation of data and SWCoeffes (W)", 1.0, t.W)

HypothesisTests.default_tail(::ShapiroWilkTest) = :left

function HypothesisTests.show_params(io::IO, t::ShapiroWilkTest, ident)
    println(io, ident, "number of observations:         ", t.SWc.N)
    println(io, ident, "censored ratio:                 ", (t.SWc.N - t.N1)/t.SWc.N)
    println(io, ident, "W-statistic:                    ", t.W)
end

HypothesisTests.pvalue(t::ShapiroWilkTest) = pvalue(t.W, t.SWc, t.N1)

function ShapiroWilkTest(
    X::AbstractArray{T}, SWc::SWCoeffs=SWCoeffs(length(X)), N1=length(X)) where {T<:Real}

    N = length(X)
    #fatal errors
    if N < 3
        throw("Need at least 3 samples.")
    elseif N1 > N
        throw("N1 must be less than or equal to length(X)")
    elseif length(SWc) ≠ length(X)
        throw("Length of the sample differs from Shapiro-Wilk coefficients!")
    end

    #non-fatal errors
    if N > 5000
        @warn("p-value may be unreliable for samples larger than 5000 points")
    elseif (N1 < N) && (N < 20)
        @warn("Number of samples is < 20. Censoring may produce unreliable p-value.")
    elseif (N - N1)/N > 0.8
        @warn("(N - N1)/N > 0.8. Censoring too much data may produce unreliable p-value.")
    end

    if !issorted(view(X, 1:N1))
        @warn("Shapiro-Wilk requires sorted data")
        W = swstat(sort(X[1:N1]), SWc)
    else
        W = swstat(view(X, 1:N1), SWc)
    end

    return ShapiroWilkTest(SWc, W, N1)
end

#=
# To compare with the standard ALGORITHM AS R94 fortran subroutine
#  * grab scipys (swilk.f)[https://github.com/scipy/scipy/blob/master/scipy/stats/statlib/swilk.f];
#  * compile
#  ```
#  gfortran -shared -fPIC -o swilk.so swilk.f
#  gfortran -fdefault-integer-8 -fdefault-real-8 -shared -fPIC swilk.f -o swilk64.so
#  ```

for (lib, I, F) in (("./swilk64.so", Int64, Float64),
                    ("./swilk.so"  , Int32, Float32))
    @eval begin
        function swilkfort!(X::AbstractVector{$F}, A::AbstractVector{$F}, computeA=true)

            w, pval = Ref{$F}(0.0), Ref{$F}(0.0)
            ifault = Ref{$I}(0)

            ccall((:swilk_, $lib),
                Void,
                (
                Ref{Bool},  # INIT if false compute SWCoeffs in A, else use A
                Ref{$F},    # X    sample
                Ref{$I},    # N    samples length
                Ref{$I},    # N1   (upper) uncensored data length
                Ref{$I},    # N2   length of A
                Ref{$F},    # A    A
                Ref{$F},    # W    W-statistic
                Ref{$F},    # PW   p-value
                Ref{$I},    # IFAULT error code (see swilk.f for meaning)
                ),
                !computeA, X, length(X), length(X), length(A), A, w, pval, ifault)
            return (w[], pval[], ifault[], A)
        end
        swilkfort(X::Vector{$F}) = swilkfort!(X, zeros($F, div(length(X),2)))
    end
end
=#


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

# UniversalFallback is needed for starting values
const CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())

import SCS

for T in [SCS.Direct, SCS.Indirect]
    optimizer = SCS.Optimizer(linear_solver=T, eps=1e-6)
    MOI.set(optimizer, MOI.Silent(), true)

    @testset "SolverName" begin
        @test MOI.get(optimizer, MOI.SolverName()) == "SCS"
    end

    @testset "supports_allocate_load" begin
        @test MOIU.supports_allocate_load(optimizer, false)
        @test !MOIU.supports_allocate_load(optimizer, true)
    end

    MOI.empty!(CACHE)
    cached = MOIU.CachingOptimizer(CACHE, optimizer)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)
    config = MOIT.TestConfig(atol=1e-5)

    @testset "Unit" begin
        MOIT.unittest(bridged, config, [
            # FIXME `NumberOfThreads` not supported.
            "number_threads",
            # `TimeLimitSec` not supported.
            "time_limit_sec",
            # ArgumentError: The number of constraints in SCSModel must be greater than 0
            "solve_unbounded_model",
            # Integer and ZeroOne sets are not supported
            "solve_integer_edge_cases", "solve_objbound_edge_cases",
            "solve_zero_one_with_bounds_1",
            "solve_zero_one_with_bounds_2",
            "solve_zero_one_with_bounds_3"])
    end

    @testset "Continuous linear problems with $T" begin
        MOIT.contlineartest(bridged, config)
    end

    @testset "ADMMIterations attribute with $T" begin
        MOIT.linear1test(bridged, config)
        @test MOI.get(bridged, SCS.ADMMIterations()) > 0
    end

    @testset "Continuous quadratic problems with $T" begin
        MOIT.qcptest(bridged, config)
    end

    @testset "Continuous conic problems with $T" begin
        MOIT.contconictest(bridged, config, ["rootdets", "logdets"])
    end
end


text:
# Index
# =====
#
# Index types for genomic intervals.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

import BGZFStreams: BGZFStream, VirtualOffset

include("chunk.jl")
include("bgzfindex.jl")
include("tabix.jl")


text:
using Flux, Test
using Tracker: gradcheck

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)

@testset "Tracker" begin

@test gradtest(Flux.mse, rand(5,5), rand(5, 5))
@test gradtest(Flux.crossentropy, rand(5,5), rand(5, 5))

@test gradtest(x -> Flux.normalise(x), rand(4,3))
@test gradtest(x -> Flux.normalise(x, dims = 2), rand(3,4))

end


text:
using CoupledDipole, Revise, Plots, DelimitedFiles

### --- INITIAL CONDITIONS ARE CHANGED TO AGREE WITH QUTIP SAVED DATA ---
function CoupledDipole.get_initial_conditions(problem::SimulationScalar)
    return zeros(ComplexF64, problem.atoms.N)
end

function CoupledDipole.get_initial_conditions(problem::SimulationMeanField)
    β₀ = zeros(ComplexF64, problem.atoms.N)
    z₀ = -ones(problem.atoms.N)
    u₀ = vcat(β₀, z₀)
    return u₀
end

### ------------ FIXED INPUTS ---------------------
r = [ -1.01015   -0.455099     0.543082;
      -2.96935   -0.0781967   -0.00764304;
      -0.714778  -0.248556     0.683987;
      -1.20512    0.376068    -1.85059;
      -1.89765    0.84797      0.248158;
      0.664911   0.00503915  -0.177228]

atoms = Cube(r, 10);

Δ = -2.0
s = 0.0022
laser = PlaneWave_3D(:z, s, Δ);

### ------------ CREATE PROBLEM AND EVOLVE ---------------------
simulationScalar = ScalarProblem(atoms, laser);
simulationMeanField = MeanFieldProblem(atoms, laser);

aS = time_evolution(simulationScalar, time_max=100);
aMF = time_evolution(simulationMeanField, time_max=100);

### ------------ GET INTENSITIES ---------------------
iS = zeros( length(aS.time_array) )
for t ∈ eachindex(iS)
    iS[t] = get_scattered_intensity(simulationScalar, aS.βₜ[t], deg2rad(35))
end

iMF = zeros( length(aMF.time_array) )
for t ∈ eachindex(iMF)
    iMF[t] = get_scattered_intensity(simulationMeanField, vcat(aMF.βₜ[t], aMF.zₜ[t]), deg2rad(35))
end

matrixData = readdlm(pwd()*"/test/benchmarks/Field_qutip.txt", '\t', Float64, '\n')
time_Qutip = matrixData[:,1]
iQutip = matrixData[:,2]

### ------------ STEADY STATES -----------------------
ss_S = get_steady_state(simulationScalar)
ss_MF = get_steady_state(simulationMeanField; time_max=1000)

i_ss_S = get_scattered_intensity(simulationScalar, ss_S, deg2rad(35))
i_ss_MF = get_scattered_intensity(simulationMeanField, ss_MF , deg2rad(35))

### ------------ PLOTS ---------------------
plot(aS.time_array, iS, label="Scalar", lw=5, ylabel="intensity", xlabel="time")
plot!(aMF.time_array, iMF, label="Mean Field", lw=4, linestyle=:dash, size=(800,600))
plot!( time_Qutip, iQutip, label="Qutip", lw=4, linestyle=:dot  )
plot!(guidefont=17, tickfont=15, legendfontsize=15, legend=:bottomright)

hline!([i_ss_S], label="i_ss_S", lw=5, linestyle=:dash)
hline!([i_ss_MF], label="i_ss_MF", lw=5, linestyle=:dot)


text:
include("POD_experiment/nous1.jl")
include("basic/gamsworld.jl")

@testset "POD instances" begin

@testset "full strong branching" begin
    println("==================================")
    println("full strong branching")
    println("==================================")

    m = batch_problem()

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            strong_branching_nsteps= 100,
            strong_branching_perc = 100,
            strong_restart = false,
            debug = true)...
    ))

    status = solve(m)
    @test status == MOI.LOCALLY_SOLVED

    juniper_val = JuMP.objective_value(m)

    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 285506.5082, atol=opt_atol, rtol=opt_rtol)

    bm = JuMP.backend(m)
    innermodel = bm.optimizer.model.inner
    debugDict = innermodel.debugDict
    counter_test(debugDict,Juniper.getnbranches(innermodel))
end

@testset "break strong branching time limit" begin
    println("==================================")
    println("break strong branching time limit")
    println("==================================")

    m = batch_problem()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            strong_branching_time_limit=0.01,
            time_limit = 4,
            strong_restart = false)...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT

    juniper_val = JuMP.objective_value(m)
    best_bound_val = JuMP.objective_bound(m)
    gap_val = getobjgap(m)
    println("Best bound: ", best_bound_val)
    println("Objective: ", juniper_val)

    # minimization problem
    @test best_bound_val <= juniper_val || isnan(juniper_val)
end

@testset "nous1 restart" begin
    println("==================================")
    println("nous1 restart")
    println("==================================")

    m = get_nous1()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            strong_restart = true,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED
end

@testset "nous1 no restart" begin
    println("==================================")
    println("nous1 no restart")
    println("==================================")

    m = get_nous1()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            strong_restart=false,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED
end

@testset "reliability" begin
    println("==================================")
    println("reliability")
    println("==================================")

    m = batch_problem()

    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:Reliability,
            reliability_branching_perc = 50,
            reliability_branching_threshold = 5,
            strong_restart = true,
            debug = true)...
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED

    juniper_val = JuMP.objective_value(m)
    best_bound_val = JuMP.objective_bound(m)
    gap_val = getobjgap(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)
    println("best_bound_val: ", best_bound_val)
    println("gap_val: ", gap_val)

    objval = 285506.5082
    @test isapprox(juniper_val, objval, atol=1e0)
    @test isapprox(best_bound_val, objval, atol=1e0)
    @test isapprox(gap_val, 0, atol=1e-2)

    bm = JuMP.backend(m)
    innermodel = bm.optimizer.model.inner
    debugDict = innermodel.debugDict
    counter_test(debugDict,Juniper.getnbranches(innermodel))
    end

end


text:
using LinearAlgebra, Arpack, FastRounding, ValidatedNumerics
using ..DynamicDefinition, ..BasisDefinition

export invariant_vector

"""
Return the (hopefully unique) invariant vector of the dynamic with discretized operator P.

The vector is normalized so that integral_covector(B)*w ≈ 1
"""
function invariant_vector(B::Basis, Q::DiscretizedOperator; tol = 0.0)
	mQ = mid(Q)
	n = size(Q)[1]
	# setting a larger nev seems to slow things down
	F = eigs(mQ; tol=tol, nev=1, ritzvec=true, v0=ones((n,)))
	w = F[2][:, 1]
	@assert imag(w) ≈ zeros(n)
	w = real(w) # this seems a pretty safe assumption.
	# In the Ulam case, in principle we could enforce w >= 0, but in practice
	# it will hardly ever be relevant.
	w = w ./ (mid.(integral_covector(B))*w) #normalization
	return w
end

"""
Return an upper bound to Q_h*w - w in the given norm
"""
function residualbound(N::Type{<:NormKind}, Q::DiscretizedOperator, w::AbstractVector)
	return normbound(N, Q*w - w)
end

"""
Bounds rigorously the distance of w from the fixed point of Q (normalized with integral = 1),
using a vector of bounds norms[k] ≥ ||Q_h^k|_{U_h^0}||.
"""
function distance_from_invariant(B::Basis, D::Dynamic, Q::DiscretizedOperator, w::AbstractVector, norms::Vector; ε₁::Float64 = residualbound(weak_norm(B), Q, w), ε₂::Float64 = mag(integral_covector(B) * w - 1), normQ::Float64 = opnormbound(weak_norm(B), Q))
	if ε₂ > 1e-8
		@error "w does not seem normalized correctly"
	end
	us = BasisDefinition.invariant_measure_strong_norm_bound(B, D)
	Cs = infinite_sum_norms(norms)
	Kh =  BasisDefinition.weak_projection_error(B)
	normw = normbound(weak_norm(B), w)

	return Cs ⊗₊ (2. ⊗₊ Kh ⊗₊ (1. ⊕₊ normQ) ⊗₊ us ⊕₊ ε₁ ⊘₊ (1. ⊖₋ ε₂)) ⊕₊ ε₂ ⊘₊ (1. ⊖₋ ε₂) ⊗₊ normw
end

# """
# This function returns a sequence of Cᵢ, \\tilde{C}ᵢ for a matrix P
# on a subspace V such that ||P^i|_V||_1\\leq C_i and
# ||P^i|_V||_{\\infty}\\leq \\tilde{C}_i, with respect to the
# """
# function contractmatrix(B::Basis, P::AbstractMatrix{Interval{T}}, m) where {T}
# 	# vector of the Cᵢ for P
# 	C = zeros(m)
# 	S = zeros((length(B), m))
# 	tilde_C = zeros(m)
#
# 	PP = mid.(P)
#
# 	for v in BasisDefinition.AverageZero(B)
# 		λ₁, λ₂ = norm(v, 1), norm(v, Inf)
# 		for i in 1:m
# 			v = PP*v
# 			C[i] = max(C[i], norm(v, 1)/λ₁)
# 			S[:, i]+=abs.(v)/λ₂
# 		end
# 	end
#
#
# 	for k in 1:m
# 		tilde_C[k] = maximum(S[:,k])
# 	end
#
# 	# we keep track of the error due to the basis we have chosen
# 	η₁ = BasisDefinition.spaceconstant(B, Val(:L1))
# 	η₂ = BasisDefinition.spaceconstant(B, Val(:L∞))
#
#  	return η₁*C, η₂*tilde_C
# end
#
# @deprecate contractmatrix norms_of_powers

"""
This function returns the bound on the weak norm of the discretized operator
"""
function boundnorm(B::Basis, P::AbstractMatrix{Interval{T}}, m) where {T}
	W₁, W₂ = BasisDefinition.bound_weak_norm_from_linalg_norm(B)
	α₁ = BasisDefinition.bound_linalg_norm_L1_from_weak(B)
	α₂ = BasisDefinition.bound_linalg_norm_L∞_from_weak(B)
	C, tilde_C = contractmatrix(B, P, m)
	return (W₁/α₁)*C+(W₂/α₂)*tilde_C
end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
######### Stan batch program example  ###########

using StanSample, MCMCChains

ProjDir = @__DIR__

dyes ="
data {
  int BATCHES;
  int SAMPLES;
  real y[BATCHES, SAMPLES];
  //vector[SAMPLES] y[BATCHES];
}

parameters {
  real<lower=0> tau_between;
  real<lower=0> tau_within;
  real theta;
  real mu[BATCHES];
}

transformed parameters {
  real sigma_between;
  real sigma_within;
  sigma_between <- 1/sqrt(tau_between);
  sigma_within <- 1/sqrt(tau_within);
}

model {
  theta ~ normal(0.0, 1E5);
  tau_between ~ gamma(.001, .001);
  tau_within ~ gamma(.001, .001);

  mu ~ normal(theta, sigma_between);
  for (n in 1:BATCHES)
    y[n] ~ normal(mu[n], sigma_within);
}

generated quantities {
  real sigmasq_between;
  real sigmasq_within;

  sigmasq_between <- 1 / tau_between;
  sigmasq_within <- 1 / tau_within;
}
"

dyesdata = Dict(:BATCHES => 6,
    :SAMPLES => 5,
    :y => [
        1545 1540 1595 1445 1595;
        1520 1440 1555 1550 1440;
        1630 1455 1440 1490 1605;
        1595 1515 1450 1520 1560;
        1510 1465 1635 1480 1580;
        1495 1560 1545 1625 1445]
  )

data = (BATCHES=6, SAMPLES=5, y=dyesdata[:y])

tmpdir = joinpath(@__DIR__, "tmp")
sm = SampleModel("dyes", dyes, tmpdir);

rc = stan_sample(sm; use_json=true, data)

if success(rc)
  chns = read_samples(sm, :mcmcchains; include_internals=true)

  chn = set_section(chns,
    Dict(
      :parameters => ["theta",
        "tau_between", "tau_within",
        "sigma_between", "sigma_within",
        "sigmasq_between", "sigmasq_within"],
      :mu => ["mu.$i" for i in 1:6],
      :internals => names(chns, [:internals])
    )
  )
  chn |> display
  println()
  describe(chn, sections=[:mu]) |> display
  println()
  describe(chn, sections=[:internals]) |> display
end


text:
# This file is a part of ArraysOfArrays.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "Package ArraysOfArrays" begin
    include("functions.jl")
    include("array_of_similar_arrays.jl")
    include("vector_of_arrays.jl")
    include("test_statsbase_support.jl")
end


text:
using PlotlyJS, Rsvg

exe_time = [2.0, 6.0, 4.2, 3.0, 23.5, 51.2]
code_line = [124, 125, 80, 75, 84, 55]
text = ["C++", "C", "Scala","Julia","Python","R"]
col = [:goldenrod, "#7F25C4", :mediumseagreen, :crimson, :lightgrey, :darkcyan]

s = scatter(;x=code_line, y=exe_time,
            text=text,
            marker=attr(color=col, size=20),
            line=attr(color=:grey, width=1));
layout1 = Layout(width=670, height=470,
                 #title="Speed vs. Code Length Tradeoff",
                 margin=attr(r=20,l=50,b=50,t=20),
                 xaxis=attr(title="Lines of Code"),
                 yaxis=attr(title="Execution Time (seconds)",zeroline=false))
p = plot(s, layout1)
savefig(p, "scatter.html")
run(`sed -i.bak '8,13!d' scatter.html`)
run(`sed -i.bak 's/newPlot('\''[a-z|0-9|-]*'\''/newPlot('\''scatplot'\''/' scatter.html`)

b1 = bar(x=text, y=code_line, marker=attr(color=col))
a = [Dict("text"=>code_line[i], "x"=>text[i], "y"=>code_line[i], "showarrow"=>false,
          "xanchor"=>"center", "yanchor"=>"bottom") for i in 1:length(exe_time)]
layout2 = Layout(width=670, height=300,
                 text=code_line,
                 annotations=a,
                 margin=attr(t=0),
                 yaxis=attr(zeroline=false,showticklabels=false,showgrid=false))
p2 = plot(b1,layout2)
savefig(p2, "bar_line.html")
run(`sed -i.bak '8,13!d' bar_line.html`)
run(`sed -i.bak 's/newPlot('\''[a-z|0-9|-]*'\''/newPlot('\''conciseness'\''/' bar_line.html`)

b2 = bar(x=text, y=exe_time, marker=attr(color=col))
a2 = [Dict("text"=>exe_time[i], "x"=>text[i], "y"=>exe_time[i], "showarrow"=>false,
          "xanchor"=>"center", "yanchor"=>"bottom") for i in 1:length(exe_time)]
layout3 = Layout(width=670, height=300,
                 text=code_line,
                 annotations=a2,
                 margin=attr(t=0),
                 yaxis=attr(zeroline=false,showticklabels=false,showgrid=false))
p3 = plot(b2,layout3)
savefig(p3, "bar_time.html")
run(`sed -i.bak '8,13!d' bar_time.html`)
run(`sed -i.bak 's/newPlot('\''[a-z|0-9|-]*'\''/newPlot('\''speed'\''/' bar_time.html`)

#plot([Plot(b1), Plot(b2)])
#=
include("graph.jl")
=#

max_stars_repo_name: luiarthur/padawan

max_issues_repo_name: luiarthur/padawan

max_forks_repo_name: luiarthur/padawan


text:
#! /Users/malcolm/bin/julia -q
#
nargs = size(ARGS)[1]
if nargs > 2
  println("usage: ftop [dir [nfiles]]")
  exit(-1)
else
  dir = (nargs > 0) ? ARGS[1] : pwd()
  nf  = (nargs > 1) ? ARGS[2] : 20
end
efind =  `find $dir -type f -iname "*" -exec du -sh "{}" + `
println("Directory: $dir")
try
  run(pipeline(efind,`sort -rh`,`head -$nf`))
  println("Done.")
catch
  println("No files found.")
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type BAT.AbstractSamplingAlgorithm

Abstract type for BAT sampling algorithms. See [`bat_sample`](@ref).
"""
abstract type AbstractSamplingAlgorithm end

"""
    bat_sample(
        [rng::AbstractRNG],
        target::BAT.AnySampleable,
        [algorithm::BAT.AbstractSamplingAlgorithm]
    )::DensitySampleVector

Draw samples from `target` using `algorithm`.

Depending on sampling algorithm, the samples may be independent or correlated
(e.g. when using MCMC).

Returns a NamedTuple of the shape

```julia
(result = X::DensitySampleVector, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_sample`, add methods to
    `bat_sample_impl` instead.
"""
function bat_sample end
export bat_sample

function bat_sample_impl end

@inline function bat_sample(rng::AbstractRNG, target::AnySampleable, algorithm::AbstractSamplingAlgorithm; kwargs...)
    r = bat_sample_impl(rng, target, algorithm; kwargs...)
    result_with_args(r, (rng = rng, algorithm = algorithm), kwargs)
end

@inline function bat_sample(target::AnySampleable; kwargs...)
    rng = bat_default_withinfo(bat_sample, Val(:rng), target)
    algorithm = bat_default_withinfo(bat_sample, Val(:algorithm), target)
    bat_sample(rng, target, algorithm; kwargs...)
end

@inline function bat_sample(target::AnySampleable, algorithm::AbstractSamplingAlgorithm; kwargs...)
    rng = bat_default_withinfo(bat_sample, Val(:rng), target)
    bat_sample(rng, target, algorithm; kwargs...)
end

@inline function bat_sample(rng::AbstractRNG, target::AnySampleable; kwargs...)
    algorithm = bat_default_withinfo(bat_sample, Val(:algorithm), target)
    bat_sample(rng, target, algorithm; kwargs...)
end

function argchoice_msg(::typeof(bat_sample), ::Val{:rng}, x::AbstractRNG)
    "Initializing new RNG of type $(typeof(x))"
end

function argchoice_msg(::typeof(bat_sample), ::Val{:algorithm}, x::AbstractSamplingAlgorithm)
    "Using sampling algorithm $x"
end

"""
    abstract type AbstractSampleGenerator

*BAT-internal, not part of stable public API.*

Abstract super type for sample generators.
"""
abstract type AbstractSampleGenerator end
export AbstractSampleGenerator