text:
############################################################################################
### Types
############################################################################################
"""
Full information maximum likelihood estimation. Can handle observed data with missings.

# Constructor

    SemFIML(;observed, specification, kwargs...)

# Arguments
- `observed::SemObservedMissing`: the observed part of the model
- `specification`: either a `RAMMatrices` or `ParameterTable` object

# Examples
```julia
my_fiml = SemFIML(observed = my_observed, specification = my_parameter_table)
```

# Interfaces
Analytic gradients are available.

# Extended help
## Implementation
Subtype of `SemLossFunction`.
"""
mutable struct SemFIML{INV, C, L, O, M, IM, I, T, U, W} <: SemLossFunction
    inverses::INV #preallocated inverses of imp_cov
    choleskys::C #preallocated choleskys
    logdets::L #logdets of implied covmats

    ∇ind::O

    imp_mean::IM
    meandiff::M
    imp_inv::I

    mult::T

    commutation_indices::U

    interaction::W
end

############################################################################################
### Constructors
############################################################################################

function SemFIML(;observed, specification, kwargs...)

    inverses = broadcast(x -> zeros(x, x), Int64.(pattern_nvar_obs(observed)))
    choleskys = Array{Cholesky{Float64,Array{Float64,2}},1}(undef, length(inverses))

    n_patterns = size(rows(observed), 1)
    logdets = zeros(n_patterns)

    imp_mean = zeros.(Int64.(pattern_nvar_obs(observed)))
    meandiff = zeros.(Int64.(pattern_nvar_obs(observed)))

    nman = Int64(n_man(observed))
    imp_inv = zeros(nman, nman)
    mult = similar.(inverses)

    ∇ind = vec(CartesianIndices(Array{Float64}(undef, nman, nman)))
    ∇ind = [findall(x -> !(x[1] ∈ ind || x[2] ∈ ind), ∇ind) for ind in patterns_not(observed)]

    commutation_indices = get_commutation_lookup(get_n_nodes(specification)^2)

    return SemFIML(
    inverses,
    choleskys,
    logdets,
    ∇ind,
    imp_mean,
    meandiff,
    imp_inv,
    mult,
    commutation_indices,
    nothing
    )
end

############################################################################################
### methods
############################################################################################

function objective!(semfiml::SemFIML, parameters, model)

    if !check_fiml(semfiml, model) return non_posdef_return(parameters) end

    prepare_SemFIML!(semfiml, model)

    objective = F_FIML(rows(observed(model)), semfiml, model, parameters)
    return objective/n_obs(observed(model))
end

function gradient!(semfiml::SemFIML, parameters, model)

    if !check_fiml(semfiml, model) return ones(eltype(parameters), size(parameters)) end

    prepare_SemFIML!(semfiml, model)

    gradient = ∇F_FIML(rows(observed(model)), semfiml, model)/n_obs(observed(model))
    return gradient
end

function objective_gradient!(semfiml::SemFIML, parameters, model)

    if !check_fiml(semfiml, model)
        return non_posdef_return(parameters), ones(eltype(parameters), size(parameters))
    end

    prepare_SemFIML!(semfiml, model)

    objective = F_FIML(rows(observed(model)), semfiml, model, parameters)/n_obs(observed(model))
    gradient = ∇F_FIML(rows(observed(model)), semfiml, model)/n_obs(observed(model))

    return objective, gradient
end

############################################################################################
### Recommended methods
############################################################################################

update_observed(lossfun::SemFIML, observed::SemObserved; kwargs...) =
    SemFIML(;observed = observed, kwargs...)

############################################################################################
### additional functions
############################################################################################

function F_one_pattern(meandiff, inverse, obs_cov, logdet, N)
    F = logdet
    F += meandiff'*inverse*meandiff
    if N > one(N)
        F += dot(obs_cov, inverse)
    end
    F = N*F
    return F
end

function ∇F_one_pattern(μ_diff, Σ⁻¹, S, pattern, ∇ind, N, Jμ, JΣ, model)
    diff⨉inv = μ_diff'*Σ⁻¹

    if N > one(N)
        JΣ[∇ind] .+= N*vec(Σ⁻¹*(I - S*Σ⁻¹ - μ_diff*diff⨉inv))
        @. Jμ[pattern] += (N*2*diff⨉inv)'

    else
        JΣ[∇ind] .+= vec(Σ⁻¹*(I - μ_diff*diff⨉inv))
        @. Jμ[pattern] += (2*diff⨉inv)'
    end

end

function ∇F_fiml_outer(JΣ, Jμ, imply::SemImplySymbolic, model, semfiml)
    G = transpose(JΣ'*∇Σ(imply)-Jμ'*∇μ(imply))
    return G
end

function ∇F_fiml_outer(JΣ, Jμ, imply, model, semfiml)

    Iₙ = sparse(1.0I, size(A(imply))...)
    P = kron(F⨉I_A⁻¹(imply), F⨉I_A⁻¹(imply))
    Q = kron(S(imply)*I_A⁻¹(imply)', Iₙ)
    #commutation_matrix_pre_square_add!(Q, Q)
    Q2 = commutation_matrix_pre_square(Q, semfiml.commutation_indices)

    ∇Σ = P*(∇S(imply) + (Q+Q2)*∇A(imply))

    ∇μ = F⨉I_A⁻¹(imply)*∇M(imply) + kron((I_A⁻¹(imply)*M(imply))', F⨉I_A⁻¹(imply))*∇A(imply)

    G = transpose(JΣ'*∇Σ-Jμ'*∇μ)

    return G
end

function F_FIML(rows, semfiml, model, parameters)
    F = zero(eltype(parameters))
    for i = 1:size(rows, 1)
        F += F_one_pattern(
            semfiml.meandiff[i],
            semfiml.inverses[i],
            obs_cov(observed(model))[i],
            semfiml.logdets[i],
            pattern_n_obs(observed(model))[i])
    end
    return F
end

function ∇F_FIML(rows, semfiml, model)
    Jμ = zeros(Int64(n_man(model)))
    JΣ = zeros(Int64(n_man(model)^2))

    for i = 1:size(rows, 1)
        ∇F_one_pattern(
            semfiml.meandiff[i],
            semfiml.inverses[i],
            obs_cov(observed(model))[i],
            patterns(observed(model))[i],
            semfiml.∇ind[i],
            pattern_n_obs(observed(model))[i],
            Jμ,
            JΣ,
            model)
    end
    return ∇F_fiml_outer(JΣ, Jμ, imply(model), model, semfiml)
end

function prepare_SemFIML!(semfiml, model)
    copy_per_pattern!(semfiml, model)
    batch_cholesky!(semfiml, model)
    #batch_sym_inv_update!(semfiml, model)
    batch_inv!(semfiml, model)
    for i in 1:size(pattern_n_obs(observed(model)), 1)
        semfiml.meandiff[i] .= obs_mean(observed(model))[i] - semfiml.imp_mean[i]
    end
end

function copy_per_pattern!(inverses, source_inverses, means, source_means, patterns)
    @views for i = 1:size(patterns, 1)
        inverses[i] .=
            source_inverses[
                patterns[i],
                patterns[i]]
    end

    @views for i = 1:size(patterns, 1)
        means[i] .=
            source_means[patterns[i]]
    end
end

copy_per_pattern!(
    semfiml,
    model::M where {M <: AbstractSem}) =
    copy_per_pattern!(
        semfiml.inverses,
        Σ(imply(model)),
        semfiml.imp_mean,
        μ(imply(model)),
        patterns(observed(model)))

function batch_cholesky!(semfiml, model)
    for i = 1:size(semfiml.inverses, 1)
        semfiml.choleskys[i] = cholesky!(Symmetric(semfiml.inverses[i]))
        semfiml.logdets[i] = logdet(semfiml.choleskys[i])
    end
    return true
end

function check_fiml(semfiml, model)
    copyto!(semfiml.imp_inv, Σ(imply(model)))
    a = cholesky!(Symmetric(semfiml.imp_inv); check = false)
    return isposdef(a)
end

get_n_nodes(specification::RAMMatrices) = specification.size_F[2]
get_n_nodes(specification::ParameterTable) =
    length(specification.variables[:observed_vars]) +
    length(specification.variables[:latent_vars])


text:
ccall((:FcInit, :libfontconfig), Uint8, ())

# By default fontconfig on OSX does not include user fonts.
if OS_NAME == :Darwin
    ccall((:FcConfigAppFontAddDir, :libfontconfig),
          Uint8, (Ptr{Void}, Ptr{Uint8}),
          C_NULL, b"~/Library/Fonts")
end

const FcMatchPattern = uint32(0)
const FcMatchFont    = uint32(1)
const FcMatchScan    = uint32(2)

function fontconfig_match(desc::String)
    pat = ccall((:FcNameParse, :libfontconfig),
                Ptr{Void}, (Ptr{Uint8},), bytestring(desc))

    ccall((:FcConfigSubstitute, :libfontconfig),
          Uint8, (Ptr{Void}, Ptr{Void}, Int32),
          C_NULL, pat, FcMatchPattern)

    ccall((:FcDefaultSubstitute, :libfontconfig),
          Void, (Ptr{Void},), pat)

    result = Int32[0]
    mat = ccall((:FcFontMatch, :libfontconfig),
                Ptr{Void}, (Ptr{Void}, Ptr{Void}, Ptr{Int32}),
                C_NULL, pat, result)

    matstr = ccall((:FcPatternFormat, :libfontconfig),
                   Ptr{Uint8}, (Ptr{Void}, Ptr{Uint8}),
                   mat, b"%{family} %{style}")
    matstr = bytestring(matstr)

    ccall((:FcPatternDestroy, :libfontconfig),
          Void, (Ptr{Void},), pat)

    ccall((:FcPatternDestroy, :libfontconfig),
          Void, (Ptr{Void},), mat)

    matstr
end

ext: jl


text:
using BrainFlow
using BrainFlowViz

function get_some_board_data(board_shim, nsamples, nchannels=8)
    data = BrainFlow.get_current_board_data(nsamples, board_shim)
    data = transpose(data)
    emg_data = view(data, :, 2:(1+nchannels))
    for chan in 1:nchannels
        emg_channel_data = view(emg_data, :, chan)
        BrainFlow.detrend(emg_channel_data, BrainFlow.CONSTANT)
    end
    return emg_data
end

### Start streaming
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)
params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.GFORCE_DUAL_BOARD, params)
BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)

nchannels = 2
nsamples = 512
data_func = ()->get_some_board_data(board_shim, nsamples, nchannels)

# start a live plotting task
t = @task BrainFlowViz.plot_data(
    data_func,
    nsamples,
    nchannels;
    y_lim = [-200 200],
    theme = :dark,
    color = :lime,
    )

sleep(0.5)
schedule(t)

sleep(4)

### to stop the task:
schedule(t, InterruptException(), error=true)

# this should go into a 'finally' of try/catch/finally
BrainFlow.stop_stream(board_shim)
BrainFlow.release_session(board_shim)


text:
using Test

m = page_model()
include("../src/components/CO2emissions.jl")

add_comp!(m, co2emissions)

set_param!(m, :co2emissions, :e0_baselineCO2emissions, readpagedata(m,"data/e0_baselineCO2emissions.csv"))
set_param!(m, :co2emissions, :er_CO2emissionsgrowth, readpagedata(m, "data/er_CO2emissionsgrowth.csv"))

##running Model
run(m)

emissions= m[:co2emissions,  :e_regionalCO2emissions]

# Recorded data
emissions_compare=readpagedata(m, "test/validationdata/e_regionalCO2emissions.csv")

@test emissions ≈ emissions_compare rtol=1e-3


text:
module SFA

export linsfa
export timediff
export whiten!

include("linsfa.jl")

end


text:
struct BatchedVector{T<:AbstractVector,P<:AbstractVector} <: AbstractVector{T}
    parent::P
    offsets::Vector{Int}
    @inline function BatchedVector{T,P}(parent, offsets) where {T<:AbstractVector,P<:AbstractVector}
        check_offsets(parent, offsets)
        new(parent, offsets)
    end
end

@inline function BatchedVector(parent::P, offsets::AbstractVector{<:Integer}) where {P<:AbstractVector}
    i1 = firstindex(parent)
    T = viewtype(parent, i1:i1)
    return BatchedVector{T,P}(parent, offsets)
end

####
#### Core Array Interface
####

Base.length(A::BatchedVector) = length(A.offsets) - 1

Base.size(A::BatchedVector) = (length(A), )

@propagate_inbounds Base.getindex(A::BatchedVector, i::Int) = view(A.parent, batchrange(A, i))

@propagate_inbounds function Base.setindex!(A::BatchedVector, x, i::Int)
    A.parent[batchrange(A, i)] = x
    return A
end

Base.IndexStyle(::Type{<:BatchedVector}) = IndexLinear()

@propagate_inbounds function batchrange(A::BatchedVector, i::Int)
    j = firstindex(A.parent)
    from = A.offsets[i] + j
    to = A.offsets[i + 1] - 1 + j
    return from:to
end

####
#### Misc
####

Base.copy(A::BatchedVector) = BatchedVector(copy(A.parent), copy(A.offsets))

Base.dataids(A::BatchedVector) = (Base.dataids(A.parent)..., Base.dataids(A.offsets)...)

Base.parent(A::BatchedVector) = A.parent

####
#### Extra
####

"""
    $(SIGNATURES)

View `A` as a vector of batches where `batch(A, batch_lengths)[i]` has
length `batch_lengths[i]`.

# Examples

```jldoctest
julia> A = collect(1:10)
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

julia> B = batch(A, [2, 3, 5])
3-element BatchedVector{SubArray{Int64,1,Array{Int64,1},Tuple{UnitRange{Int64}},true},Array{Int64,1}}:
 [1, 2]
 [3, 4, 5]
 [6, 7, 8, 9, 10]

julia> B[2] == view(A, 3:5)
true
```
"""
function batch(A::AbstractVector, batch_lengths::AbstractVector{<:Integer})
    offsets = Vector{Int}(undef, length(batch_lengths) + 1)
    offsets[1] = cumsum = 0
    for (i, l) in enumerate(batch_lengths)
        cumsum += l
        offsets[i + 1] = cumsum
    end
    return BatchedVector(A, offsets)
end

"""
    $(SIGNATURES)

View `A` as a vector of batches using the same batch lengths as `B`.

See also: [`batch`](@ref).
"""
@inline function similarbatch(A::AbstractVector, B::BatchedVector)
    length(A) == length(B.parent) || argerror("length(A) != length(parent(B))")
    return BatchedVector(A, copy(B.offsets))
end

"""
    $(SIGNATURES)

Returns the unbatched parent array wrapped by `B`.
"""
SpecialArrays.flatten(B::BatchedVector) = B.parent

####
#### 3rd Party
####

@inline function UnsafeArrays.unsafe_uview(A::BatchedVector)
    BatchedVector(uview(A.parent), A.offsets)
end

####
#### Util
####

check_offsets(A::BatchedVector) = check_offsets(A.parent, A.offsets)

function check_offsets(parent::AbstractVector, offsets::AbstractVector{<:Integer})
    length(offsets) >= 1 || argerror("offsets cannot be empty")
    first(offsets) == 0 || argerror("First offset is non-zero")
    for i in LinearIndices(offsets)[2:end]
        o1 = offsets[i-1]
        o2 = offsets[i]
        # TODO 0 length slice
        o2 > o1 || argerror("Overlapping indices found in offsets")
    end
    # TODO support batching < length(parent)
    offsets[end] == length(parent) || argerror("offsets[end] does not equal to length(parent)")
    return nothing
end


text:
"""
Defines functions and constants managing the Kroki service the rest of the
package uses to render diagrams. These services can be either local or remote.

This module also enables management of a local service instance, provided
[Docker](https://docker.com) and [Docker
Compose](https://docs.docker.com/compose/) are available on the system.
"""
module Service

using ..Documentation
@setupDocstringMarkup()

"""
The default [`ENDPOINT`](@ref) to use, i.e. the [publicly hosted
version](https://kroki.io).
"""
const DEFAULT_ENDPOINT = "https://kroki.io"

"""
A specialized `Exception` to include reporting instructions for specific types
of errors that may occur while trying to execute `docker-compose`.
"""
struct DockerComposeExecutionError <: Exception
  message::String
end
Base.showerror(io::IO, error::DockerComposeExecutionError) = print(
  io,
  """
An error occurred while executing `docker-compose`.

This may be caused by a change in its interface. If you believe this error to
be caused by Kroki.jl itself instead of a configuration error on the system,
please file an issue through https://github.com/bauglir/Kroki.jl/issues/new
including the error message below and a description of when the error occurred.

The reported error was:

$(error.message)
""",
)

"The currently active Kroki service endpoint being used."
const ENDPOINT = Ref{String}("https://kroki.io")

"Path to the Docker Compose definitions for running a local Kroki service."
const SERVICE_DEFINITION_FILE = realpath("$(@__DIR__)/../../support/docker-services.yml")

"""
Helper function for executing Docker Compose commands.

Returns captured stdout.

Throws an `ErrorException` if Docker and/or Docker Compose aren't available.
Throws a `DockerComposeExecutionError` if any other exception occurs during
execution.
"""
function executeDockerCompose(cmd::Vector{String})
  captured_stderr = IOBuffer()
  captured_stdout = IOBuffer()

  try
    run(pipeline(
      `docker-compose --file $(SERVICE_DEFINITION_FILE) --project-name krokijl $cmd`;
      stderr = captured_stderr,
      stdout = captured_stdout,
    ))
  catch exception
    exception isa Base.IOError &&
    throw(ErrorException("Missing dependencies! Docker and/or Docker Compose do not appear to be available"))

    throw(DockerComposeExecutionError(String(take!(captured_stderr))))
  end

  String(take!(captured_stdout))
end
executeDockerCompose(cmd::String) = executeDockerCompose([cmd])
# The function should be called through indirection in cases where it needs to
# be mocked out in tests
const EXECUTE_DOCKER_COMPOSE = Ref{Any}(executeDockerCompose)

"""
Sets the [`ENDPOINT`](@ref) using a fallback mechanism if no `endpoint` is
provided.

The fallback mechanism checks for a `KROKI_ENDPOINT` environment variable
specifying an endpoint (e.g. to be used across Julia instances). If this
environment variable is also not present the [`DEFAULT_ENDPOINT`](@ref) is
used.

This can, for instance, be used in cases where a [privately hosted
instance](https://docs.kroki.io/kroki/setup/install/) is available or when a
local service has been [`start!`](@ref)ed.

Returns the value that [`ENDPOINT`](@ref) got set to.

# Examples

- `setEndpoint!()`
- `setEndpoint!("http://localhost:8000")`
"""
function setEndpoint!(
  endpoint::AbstractString = get(ENV, "KROKI_ENDPOINT", DEFAULT_ENDPOINT),
)
  ENDPOINT[] != endpoint && @info "Setting Kroki service endpoint to $(endpoint)."
  ENDPOINT[] = String(endpoint)
end

"""
Starts the Kroki service components on the local system, optionally, ensuring
[`ENDPOINT`](@ref) points to them.

Pass `false` to the function to prevent the [`ENDPOINT`](@ref) from being
updated. The default behavior is to update.
"""
function start!(update_endpoint::Bool = true)
  @info "Starting Kroki service components."
  EXECUTE_DOCKER_COMPOSE[](["up", "--detach"])
  update_endpoint && setEndpoint!("http://localhost:8000")
  return
end

"""
Returns a `NamedTuple` where the keys are the names of the service components
and the values their corresponding 'running' state.

# Examples
```
julia> status()
(core = true, mermaid = false)
```
"""
function status()
  states_to_bool = Pair{String, Bool}["running" => true, "stopped" => false]

  services_by_state = map(states_to_bool) do (state, running)::Pair{String, Bool}
    services_with_state =
      EXECUTE_DOCKER_COMPOSE[](["ps", "--filter", "status=$state", "--services"])

    service_symbols = Symbol.(split(services_with_state, '\n'; keepempty = false))

    (; zip(service_symbols, Iterators.repeated(running))...)
  end

  merge(services_by_state...)
end

"""
Stops any running Kroki service components ensuring [`ENDPOINT`](@ref) no
longer points to the stopped service.

Cleans up left-over containers by default. This behavior can be turned off by
passing `false` to the function.
"""
function stop!(perform_cleanup::Bool = true)
  @info "Stopping Kroki service components."
  EXECUTE_DOCKER_COMPOSE[]("stop")
  perform_cleanup && EXECUTE_DOCKER_COMPOSE[](["rm", "--force"])
  setEndpoint!()
  return
end

"Updates the Docker images for the individual Kroki service components."
function update!()
  EXECUTE_DOCKER_COMPOSE[](["pull", "--quiet"])
  return
end

# Ensure `ENDPOINT` has a valid value
__init__() = setEndpoint!()

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

pathToFMU = joinpath(dirname(@__FILE__), "../model/SpringPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

@test fmiInstantiate!(myFMU; loggingOn=false) != 0

dt = 0.01
t_start = 0.0
t_stop = 8.0

data = fmiSimulateCS(myFMU, dt, t_start, t_stop, ["mass.s", "mass.v"])
@test length(data.dataPoints) > 0
@test length(data.dataPoints[1]) > 0
# ToDo: Time series comparision

fmiUnload(myFMU)


text:
"A fixed-length vector whose data is stored in some manually managed region of memory."
struct ManualVector{T} <: AbstractArray{T, 1}
    ptr::Manual{T}
    length::Int64

    function ManualVector{T}(ptr::Manual{T}, length::Int64) where {T}
        @assert isbits(T)
        new(ptr, length)
    end
end

@inline function get_address(pv::ManualVector{T}, i::Int) where {T}
    (0 < i <= pv.length) || throw(BoundsError(pv, i))
    Manual{T}(pv.ptr.ptr + (i-1)*sizeof(T))
end

function get_address(pv::Manual{ManualVector{T}}, i::Int) where {T}
    get_address(unsafe_load(pv), i)
end

# array interface

function Base.size(pv::ManualVector)
    (pv.length,)
end

@inline function Base.getindex(pv::ManualVector{T}, i::Int) where {T}
    unsafe_load(get_address(pv, i))
end

@inline function Base.setindex!(pv::ManualVector{T}, v, i::Int) where {T}
    unsafe_store!(get_address(pv, i), v)
end

function Base.IndexStyle(_::Type{ManualVector{T}}) where {T}
    Base.IndexLinear()
end

# copying, with correct handling of overlapping regions
# TODO use memcopy
function Base.copy!(dest::ManualVector{T}, doff::Int, src::ManualVector{T},
    soff::Int, n::Int) where {T}
    if doff < soff
        for i in 0:n-1 dest[doff+i] = src[soff+i] end
    else
        for i in n-1:-1:0 dest[doff+i] = src[soff+i] end
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 36511846-ff42-11ea-1c58-a1739934bbea
using Pkg, DrWatson

# ╔═╡ 36515568-ff42-11ea-3e8d-1955d756177a
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, GLM
	using StatisticalRethinking
end

# ╔═╡ 3656a862-ff42-11ea-34ea-39c9fb1daca2
include(projectdir("models", "06", "m6.11s.jl"));

# ╔═╡ b1fb477a-ff40-11ea-1283-e7c61d1ea1d5
md"## Clip-06-28s.jl"

# ╔═╡ 36604d54-ff42-11ea-00cf-8f348cf3a0a9
begin
	N = 200
	b_GP = 1                               # Direct effect of G on P
	b_GC = 0                               # Direct effect of G on C
	b_PC = 1                               # Direct effect of P on C
	b_U = 2                                # Direct effect of U on P and C
	df = DataFrame(
	  :u => 2 * rand(Bernoulli(0.5), N) .- 1,
	  :g => rand(Normal(), N)
	)
	df[!, :p] = [rand(Normal(b_GP * df[i, :g] + b_U * df[i, :u]), 1)[1] for i in 1:N]
	df[!, :c] = [rand(Normal(b_PC * df[i, :p] + b_GC * df[i, :g] + b_U * df[i, :u]), 1)[1] for i in 1:N]
	Text(precis(df; io=String))
end

# ╔═╡ 366b3adc-ff42-11ea-0b42-fd1d057bd787
stan6_12 = "
data {
  int <lower=0> N;
  vector[N] C;
  vector[N] P;
  vector[N] G;
  vector[N] U;
}
parameters {
  real <lower=0> sigma;
  real a;
  real b_PC;
  real b_GC;
  real b_U;
}
model {
  vector[N] mu;
  sigma ~ exponential(1);
  a ~ normal(0, 1);
  b_PC ~ normal(0, 1);
  b_GC ~ normal(0, 1);
  b_U ~ normal(0, 1);
  mu = a + b_PC * P + b_GC * G + b_U * U;
  C ~ normal(mu, sigma);
}
";

# ╔═╡ 366c251e-ff42-11ea-39e5-c3e32f574c75
begin
	m6_12s = SampleModel("m6.12s", stan6_12)
	m6_12_data = Dict(
	  :N => nrow(df),
	  :C => df[:, :c],
	  :P => df[:, :p],
	  :G => df[:, :g],
	  :U => df[:, :u]
	)
	rc6_12s = stan_sample(m6_12s, data=m6_12_data)
	post6_12s_df = read_samples(m6_12s, :dataframe)
	Text(precis(post6_12s_df; io=String))
end

# ╔═╡ 3676ee36-ff42-11ea-200c-594878417287
if success(rc6_12s)
	(s6_12s, p6_12s) = plot_model_coef([m6_11s, m6_12s], [:a, :b_PC, :b_GC, :b_U])
	p6_12s
end

# ╔═╡ 42995640-ff42-11ea-280f-0dbbffc082e5
s6_12s

# ╔═╡ 367794a8-ff42-11ea-17d8-01eebed2403d
if success(rc6_12s)
  scale!(df, [:g, :c, :p])
  q = quantile(df[:, :p_s], [0.45, 0.60])

  function is_in(v::Vector, q::Vector)
    findall(x -> q[1] < x < q[2], v)
  end

  v = is_in(df[:, :p], q)

  df1 = DataFrame(
    :c_s => df[v, :c_s],
    :g_s => df[v, :g_s]
  )

  ols = lm(@formula(c_s ~ g_s), df1)

  plot(xlab="grandparent education (g)", ylab="grandchild education (c)", leg=false)
  for i in 1:nrow(df)
    if df[i, :u] == -1
      if df[i, :p] >= q[1] && df[i, :p] <= q[2]
        scatter!([df[i, :g_s]], [df[i, :c_s]], color=:blue)
      else
        scatter!([df[i, :g_s]], [df[i, :c_s]], color=:lightgrey)
      end
    else
      if df[i, :p] >= q[1] && df[i, :p] <= q[2]
        scatter!([df[i, :g_s]], [df[i, :c_s]], color=:blue)
      else
        scatter!([df[i, :g_s]], [df[i, :c_s]], color=:lightblue)
      end
    end
  end

  x = -3:0.01:3
  plot!(x, coef(ols)[1] .+ coef(ols)[2] * x, color=:black)

end

# ╔═╡ 69592a8e-ff43-11ea-142f-4bdbe33da88f
md"## End of clip-06-28s.jl"

# ╔═╡ Cell order:
# ╟─b1fb477a-ff40-11ea-1283-e7c61d1ea1d5
# ╠═36511846-ff42-11ea-1c58-a1739934bbea
# ╠═36515568-ff42-11ea-3e8d-1955d756177a
# ╠═3656a862-ff42-11ea-34ea-39c9fb1daca2
# ╠═36604d54-ff42-11ea-00cf-8f348cf3a0a9
# ╠═366b3adc-ff42-11ea-0b42-fd1d057bd787
# ╠═366c251e-ff42-11ea-39e5-c3e32f574c75
# ╠═3676ee36-ff42-11ea-200c-594878417287
# ╠═42995640-ff42-11ea-280f-0dbbffc082e5
# ╠═367794a8-ff42-11ea-17d8-01eebed2403d
# ╟─69592a8e-ff43-11ea-142f-4bdbe33da88f


text:
# probably should be rewritten in a more fuctional way

const FILE_TYPE = "fio.file"

function AE_open(line::Ref{UInt}, ns::Namespace, filename::AetherlangObject{String}, mode::AetherlangObject{String})::AetherlangObject
    try
        AetherlangObject{IO}(open(dirname(ARGS[1])*"/"*filename.dataref, mode.dataref), Ref(FILE_TYPE))
    catch
        throw(AetherlangError("Error when opening the file $(filename.dataref) for $(mode.dataref)"))
    end
end

function AE_close(line::Ref{UInt}, ns::Namespace, file::AetherlangObject)::AetherlangObject
    close(file.dataref)
    AetherlangObject()
end

function AE_fprint(line::Ref{UInt}, ns::Namespace, file::AetherlangObject, args...)::AetherlangObject
    try
        for i in 1:length(args)-1
            write(file.dataref, string(args[i])*' ')
        end
        write(file.dataref, string(args[end]))
        AetherlangObject()
    catch
        throw(AetherlangError("Error when writing to file"))
    end
end
function AE_fprintln(line::Ref{UInt}, ns::Namespace, file::AetherlangObject, args...)::AetherlangObject
    try
        for i in 1:length(args)-1
            write(file.dataref, string(args[i])*' ')
        end
        write(file.dataref, string(args[end])*'\n')
        AetherlangObject()
    catch
        throw(AetherlangError("Error when writing to file"))
    end
end
function AE_finput(line::Ref{UInt}, ns::Namespace, file::AetherlangObject)::AetherlangObject
    AetherlangObject(readline(file.dataref))
end

fio_namespace_modify = Namespace(
    "open" => AetherlangObject(AE_open),
    "close" => AetherlangObject(AE_close),
    "fprintln" => AetherlangObject(AE_fprintln),
    "fprint" => AetherlangObject(AE_fprint),
    "finput" => AetherlangObject(AE_finput)
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
# adjacency list

#################################################
#
#   GenericAdjacencyList{V, VList, AList}
#
#   V:          vertex type
#   VList:      the type of vertex list
#   AdjList:    adjacency list
#
#################################################

type GenericAdjacencyList{V, VList, AdjList} <: AbstractGraph{V, Edge{V}}
    is_directed::Bool
    vertices::VList
    nedges::Int
    adjlist::AdjList
end

typealias SimpleAdjacencyList GenericAdjacencyList{Int, Range1{Int}, Vector{Vector{Int}}}
typealias AdjacencyList{V} GenericAdjacencyList{V, Vector{V}, Vector{Vector{V}}}

@graph_implements GenericAdjacencyList vertex_list vertex_map adjacency_list

## construction

simple_adjlist(nv::Int; is_directed::Bool=true) = SimpleAdjacencyList(is_directed, 1:nv, 0, multivecs(Int, nv))

adjlist{V}(vs::Vector{V}; is_directed::Bool=true) = AdjacencyList{V}(is_directed, vs, 0, Vector{V}[])
adjlist{V}(::Type{V}; is_directed::Bool=true) = adjlist(V[]; is_directed=is_directed)

## required interfaces

is_directed(g::GenericAdjacencyList) = g.is_directed

num_vertices(g::GenericAdjacencyList) = length(g.vertices)
vertices(g::GenericAdjacencyList) = g.vertices
vertex_index{V<:ProvidedVertexType}(v::V, g::GenericAdjacencyList{V}) = vertex_index(v)

num_edges(g::GenericAdjacencyList) = g.nedges

out_degree{V}(v::V, g::GenericAdjacencyList{V}) = length(g.adjlist[vertex_index(v,g)])
out_neighbors{V}(v::V, g::GenericAdjacencyList{V}) = g.adjlist[vertex_index(v,g)]

## mutation

function add_vertex!{V}(g::GenericAdjacencyList{V}, v::V)
    push!(g.vertices, v)
    push!(g.adjlist, Array(V,0))
    v
end
add_vertex!(g::GenericAdjacencyList, x) = add_vertex!(g, make_vertex(g, x))

function add_edge!{V}(g::GenericAdjacencyList{V}, u::V, v::V)
    nv::Int = num_vertices(g)
    iu = vertex_index(u, g)::Int
    push!(g.adjlist[iu], v)
    g.nedges += 1
    if !g.is_directed
        iv = vertex_index(v, g)::Int
        push!(g.adjlist[iv], u)
    end
end

## constructing from matrices

function simple_adjlist{T<:Number}(A::AbstractMatrix{T}; is_directed::Bool=true)
    n = size(A, 1)
    size(A, 2) == n || error("A must be square")
    alist = multivecs(Int, n)
    m = 0

    @inbounds if is_directed
        for i = 1:n
            nbs = alist[i]::Vector{Int}
            for j = 1:n
                if isnz(A[i,j])
                    push!(nbs, j)
                    m += 1
                end
            end
        end
    else
        for i = 1:n
            nbs = alist[i]::Vector{Int}
            if isnz(A[i,i])
                push!(nbs, i)
                m += 1
            end
            for j = i+1:n
                if isnz(A[i,j])
                    push!(nbs, j)
                    push!(alist[j], i)
                    m += 1
                end
            end
        end
    end
    return SimpleAdjacencyList(is_directed, 1:n, m, alist)
end


text:
Maybe{T} = Union{T,Nothing}

issubmodule(m::Module, s) = s isa Module && parentmodule(s) == m && m != s

function submodules(m::Module)
    nms = filter!(names(m, all=true)) do y
        Base.isdefined(m, y) && !Base.isdeprecated(m, y)
    end
    symbols = Core.eval.(Ref(m), nms)
    filter!(x -> issubmodule(m, x), symbols)
end

# list of recursive submodules of m, including m itself
function recsubmodules(m::Module)
    subs = submodules(m)
    if isempty(subs)
        recs = subs
    else
        recs = mapreduce(recsubmodules, vcat, subs)
    end
    pushfirst!(recs, m)
    recs
end

function allequal(xs)
    local val
    for x in xs
        if !@isdefined(val)
            val = x
        else
            isequal(val, x) || return false
        end
    end
    true
end

# test if a module or one of its parents was replaced
function is_replaced(mod::Module)
    par = parentmodule(mod)
    while par != mod
        getfield(par, nameof(mod)) != mod && return true
        mod = par
        par = parentmodule(par)
    end
    false
end

# list of recursive parent modules of mod, starting with mod
function parentmodules(mod::Module)
    pars = [mod]
    while true
        par = parentmodule(mod)
        par == mod && return pars
        push!(pars, par)
        mod = par
    end
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
# pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
# const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
# includet(joinpath(LOOPVECBENCHDIR, "driver.jl"))

using Distributed, LoopVectorization, JLD2

const LOOPVECBENCHDIR = joinpath(pkgdir(LoopVectorization), "benchmark")
include(joinpath(LOOPVECBENCHDIR, "benchmarkflops.jl"))
include(joinpath(LOOPVECBENCHDIR, "plotbenchmarks.jl"))

nprocs_to_add() = (Sys.CPU_THREADS >> 1) - 1
start_worker(wid) = remotecall(include, wid, joinpath(LOOPVECBENCHDIR, "setup_worker.jl"))
function start_workers()
    addprocs(nprocs_to_add())
    foreach(wait, map(start_worker, workers()))
end
stop_workers() = rmprocs(workers())

blastests() = [
    "LoopVectorization",
    "Julia", "Clang",
    "GFortran", "icc", "ifort",
    "g++ & Eigen-3", "clang++ & Eigen-3",
    "GFortran-builtin", "ifort-builtin",
    "OpenBLAS", "MKL"
]
function benchmark_AmulB(sizes)
    tests = blastests()
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> A_mul_B_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_AmulBt(sizes)
    tests = blastests()
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> A_mul_Bt_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_AtmulB(sizes)
    tests = blastests()
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> At_mul_B_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_AtmulBt(sizes)
    tests = blastests()
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> At_mul_Bt_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_dot(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3"]#, "OpenBLAS"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> dot_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_selfdot(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3"]#, "OpenBLAS"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> selfdot_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_Amulvb(sizes)
    tests = blastests()
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> A_mul_vb_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_Atmulvb(sizes)
    tests = blastests()
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> At_mul_vb_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_dot3(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3", "LinearAlgebra" ]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> dot3_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_sse(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3", "MKL"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> sse_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_exp(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> exp_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_aplusBc(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> aplusBc_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_AplusAt(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3", "GFortran-builtin", "ifort-builtin"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> AplusAt_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_random_access(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> randomaccess_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_logdettriangle(sizes)
    # tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "g++ & Eigen-3", "clang++ & Eigen-3", "LinearAlgebra"]
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort", "LinearAlgebra"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> logdettriangle_bench!(sm, is[2], is[1]), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_filter2d(sizes, K)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    pmap(is -> filter2d_bench_run!(sm, is[2], is[1], K), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end
function benchmark_filter2ddynamic(sizes)
    K = OffsetArray(rand(Float64, 3, 3), -1:1, -1:1)
    benchmark_filter2d(sizes, K)
end
function benchmark_filter2d3x3(sizes)
    K = SizedOffsetMatrix{Float64,-1,1,-1,1}(rand(3,3))
    benchmark_filter2d(sizes, K)
end
function benchmark_filter2dunrolled(sizes)
    tests = ["LoopVectorization", "Julia", "Clang", "GFortran", "icc", "ifort"]
    start_workers()
    sm = SharedMatrix(Matrix{Float64}(undef, length(tests), length(sizes)))
    K = SizedOffsetMatrix{Float64,-1,1,-1,1}(rand(3,3))
    pmap(is -> filter2dunrolled_bench_run!(sm, is[2], is[1], K), enumerate(sizes))
    br = BenchmarkResult(Matrix(sm), tests, sizes)
    stop_workers()
    br
end

# sizes = 23:23
sizes = 256:-1:2
longsizes = 1024:-1:2

logdettriangle_bench = benchmark_logdettriangle(sizes); println("logdet(LowerTriangular(A)) benchmark results:"); println(logdettriangle_bench)
dot3_bench = benchmark_dot3(sizes); println("x' * A * y benchmark results:"); println(dot3_bench)

AmulB_bench = benchmark_AmulB(sizes); println("A * B benchmark results:"); println(AmulB_bench)
AmulBt_bench = benchmark_AmulBt(sizes); println("A * B' benchmark results:"); println(AmulBt_bench)
AtmulBt_bench = benchmark_AtmulBt(sizes); println("A' * B' benchmark results:"); println(AtmulBt_bench)
AtmulB_bench = benchmark_AtmulB(sizes); println("A' * B benchmark results:"); println(AtmulB_bench)

Amulvb_bench = benchmark_Amulvb(sizes); println("A * b benchmark results:"); println(Amulvb_bench)
Atmulvb_bench = benchmark_Atmulvb(sizes); println("A' * b benchmark results:"); println(Atmulvb_bench)

dot_bench = benchmark_dot(longsizes); println("a' * b benchmark results:"); println(dot_bench)
selfdot_bench = benchmark_selfdot(longsizes); println("a' * a benchmark results:"); println(selfdot_bench)

sse_bench = benchmark_sse(sizes); println("Benchmark resutls of summing squared error:"); println(sse_bench)
aplusBc_bench = benchmark_aplusBc(sizes); println("Benchmark results of a .+ B .* c':"); println(aplusBc_bench)
AplusAt_bench = benchmark_AplusAt(sizes); println("Benchmark results of A * A':"); println(AplusAt_bench)

filter2d_dynamic_bench = benchmark_filter2ddynamic(sizes); println("Benchmark results for dynamically sized 3x3 convolution:"); println(filter2d_dynamic_bench)
filter2d_3x3_bench = benchmark_filter2d3x3(sizes); println("Benchmark results for statically sized 3x3 convolution:"); println(filter2d_3x3_bench)
filter2d_unrolled_bench = benchmark_filter2dunrolled(sizes); println("Benchmark results for unrolled 3x3 convolution:"); println(filter2d_unrolled_bench)

vexp_bench = benchmark_exp(sizes); println("Benchmark results of exponentiating a vector:"); println(vexp_bench)
randomaccess_bench = benchmark_random_access(sizes); println("Benchmark results from using a vector of indices:"); println(randomaccess_bench)

const v = 1
using Cairo, Fontconfig
const PICTURES = joinpath(pkgdir(LoopVectorization), "docs", "src", "assets")
saveplot(f, br) = draw(PNG(joinpath(PICTURES, f * "$v.png"), 12inch, 8inch), plot(br))

# If only rerunning a few, remove them from load.
# @load "benchmarkresults.jld2" logdettriangle_bench filter2d_dynamic_bench filter2d_3x3_bench filter2d_unrolled_bench dot_bench selfdot_bench dot3_bench sse_bench aplusBc_bench AplusAt_bench vexp_bench randomaccess_bench AmulB_bench AmulBt_bench AtmulB_bench AtmulBt_bench Amulvb_bench Atmulvb_bench

saveplot("bench_dot3_v", dot3_bench);
saveplot("bench_dot_v", dot_bench);
saveplot("bench_selfdot_v", selfdot_bench);
saveplot("bench_sse_v", sse_bench);
saveplot("bench_aplusBc_v", aplusBc_bench);
saveplot("bench_AplusAt_v", AplusAt_bench);
saveplot("bench_AmulB_v", AmulB_bench);
saveplot("bench_AmulBt_v", AmulBt_bench);
saveplot("bench_AtmulB_v", AtmulB_bench);
saveplot("bench_AtmulBt_v", AtmulBt_bench);
saveplot("bench_Amulvb_v", Amulvb_bench);
saveplot("bench_Atmulvb_v", Atmulvb_bench);

saveplot("bench_logdettriangle_v", logdettriangle_bench);
saveplot("bench_filter2d_dynamic_v", filter2d_dynamic_bench);
saveplot("bench_filter2d_3x3_v", filter2d_3x3_bench);
saveplot("bench_filter2d_unrolled_v", filter2d_unrolled_bench);
saveplot("bench_exp_v", vexp_bench);
saveplot("bench_random_access_v", randomaccess_bench);

@save "benchmarkresults.jld2" logdettriangle_bench filter2d_dynamic_bench filter2d_3x3_bench filter2d_unrolled_bench dot_bench selfdot_bench dot3_bench sse_bench aplusBc_bench AplusAt_bench vexp_bench randomaccess_bench AmulB_bench AmulBt_bench AtmulB_bench AtmulBt_bench Amulvb_bench Atmulvb_bench


text:
function metropolis_sampler(;target_pdf, known_pdf, known_sampler, starting_state, n_iter, n_burn = 100, n_thinning = 100)

    function transition_probability(target_pdf, known_pdf, new_state, previous_state)

        """metropolis_independant_sampler transition probability, see eq 1 in https://arxiv.org/abs/1705.00686"""
        min(1, target_pdf(new_state) * known_pdf(previous_state)/(target_pdf(previous_state) * known_pdf(new_state)))
    end

    function do_with_probability(p)

        """returns true with probability p, false with (1-p)"""

        rand() < p ? true : false

    end

    function update_state_parameters(target_pdf, known_pdf, new_state, previous_state)

        if do_with_probability(transition_probability(target_pdf, known_pdf, new_state, previous_state))

            return new_state
        else
            return previous_state
        end

    end

    function metropolis_iteration!(target_pdf, known_pdf, known_sampler, new_state, previous_state)
        new_state[:] = known_sampler()[:]
        new_state[:] = update_state_parameters(target_pdf, known_pdf, new_state, previous_state)[:]
        previous_state[:] = new_state[:]
    end

    samples = Array{eltype(starting_state)}(undef, n_iter, size(starting_state,1))

    previous_state = similar(starting_state) #x
    new_state = similar(previous_state)
    previous_state = starting_state

    if n_burn == n_thinning
        @showprogress for iter in 1:n_iter
            for i in 1:n_thinning
                metropolis_iteration!(target_pdf, known_pdf, known_sampler, new_state, previous_state)
            end
            samples[iter,:] = new_state[:]
        end
    else
        throw(Exception("n_burn != n_thinning not implemented"))
    end

    samples

end


text:
module test_Plot6

import ModiaMath

# Desired:
#   using Test
#
# In order that Test needs not to be defined in the user environment, it is included via ModiaMath:
using ModiaMath.Test
t = range(0.0, stop=10.0, length=100)

# Define signals
result = Dict{AbstractString,Any}("time" => t      , "phi1" => sin.(t)     , "phi2" => 0.5 * sin.(t),
                                  "w1"   => cos.(t), "w2"   => 0.6 * cos.(t))

# Plots
ModiaMath.plot(result, :phi1)                                        # 1 signal in one diagram
ModiaMath.plot(result, (:phi1, :phi2, :w1), figure=2)                # 3 signals in one diagram
ModiaMath.plot(result, [:phi1, :phi2, :w1], figure=3)                # 3 diagrams in form of a vector (every diagram has one signal)
ModiaMath.plot(result, ["phi1" "phi2";
                        "w1"   "w2"   ], figure=4)                   # 4 diagrams in form of a matrix (every diagram has one signal)
ModiaMath.plot(result, [ (:phi1, :phi2), (:w1) ], figure=5)           # 2 diagrams in form of a vector
ModiaMath.plot(result, [ (:phi1,)           (:phi2, :w1);
                         (:phi1, :phi2, :w1)  (:w2,)     ], figure=6)   # 4 diagrams in form of a matrix
ModiaMath.plot(result, :w1, xAxis=:phi1, figure=7)                   # Plot w1=f(phi1) in one diagram

# Print result variables
println("\n... result variables = ", ModiaMath.resultTable(result))

end


text:
using FFTW,LinearAlgebra,Plots,SpecialFunctions; function pic()
N=64;P=16N;ϵ=1/10;dt=ϵ^2/N/2;T=4048;
W=400;w=W/P*N;Z=zeros;ξ=im*Z(N);x=rand(P);a=2π*rand(P);v=randn(P,3);
O()=(ξ[1]*=false;ξ);r=Z(N);
Ex=Z(N);Ey=Z(N);Ez=Z(N);Bx=Z(N);By=Z(N);Bz=Z(N).+1;X=copy(x);V=copy(v);
Fx=Z(N);Fy=Z(N);Fz=Z(N);Cy=Z(N);Cz=Z(N);
ik=2π*im*vcat(0,1:N/2,-N/2+1:-1);D=Z(T,4);irn(c)=Int(mod1(round(c*N),N));
f(g,c)=erf((g-c)*N)/2;ff(i,c)=f((i+0.5)/N,c)-f((i-0.5)/N,c);η=1e-8;
s(c,u=1)=((mod1(i,N),u*ff(i,c)) for i∈(-6:6).+irn(c));
d(y,u=0*y.+1)=(r.*=0;for j∈s.(y,u);for k∈j;r[k[1]]+=k[2]*w;end;end;r);
fields = zeros(N,T,12);rif(x)=real.(ifft(x));
for t ∈1:T;
  X.=x;V.=v; Ex.=real.(ifft((ξ.=fft(d(X))./ik;O())));
  x.=X.+v[:,1]*dt; x.=mod.(x,1);p=sortperm(x);x=x[p];v=v[p,:];
  Cy.=By; Cz.=Bz; By.-=dt.*rif(-ik.*fft(Ez)); Bz.-=dt.*rif(ik.*fft(Ey));
  by=(By.+Cy)./2; bz=(Bz.+Cz)./2;
  for j∈1:P;
    θ=s(x[j]); a=sum(k->[Ex[k[1]],Ey[k[1]],Ez[k[1]]]*k[2],θ);
    v[j,:].+=a*dt/2
    v[j,:].+=cross(v[j,:],sum(k->[Bx[k[1]],by[k[1]],bz[k[1]]]*k[2],θ))*dt;
    v[j,:].+=a*dt/2
  end
  Ey.+=dt.*rif(-ik.*fft(Bz)/ϵ^2 .-fft(d((X.+x)/2,v[:,2])));
  Ez.+=dt.*rif(+ik.*fft(By)/ϵ^2 .-fft(d((X.+x)/2,v[:,3])));
  fields[:,t,1].=Ex;fields[:,t,2].=Ey;fields[:,t,3].=Ez;fields[:,t,4].=By;
  fields[:,t,5].=Bz;fields[:,t,6].=d(x);fields[:,t,7].=d(x,v[:,1]);
  fields[:,t,8].=d(x, v[:,2]);fields[:,t,9].=d(x, v[:,3]);
  fields[:,t,10].=d(x, v[:,1].^2);fields[:,t,11].=d(x, v[:,2].^2);
  fields[:,t,12].=d(x, v[:,3].^2);
end;return fields;end


text:
function thomas(a::Vector{Float64}, b::Vector{Float64}, c::Vector{Float64},
                d::Vector{Float64}, n::Int64)

    x = copy(d)
    c_prime = copy(c)

    # Setting initial elements
    c_prime[1] /= b[1]
    x[1] /= b[1]

    for i = 2:n
        # Scale factor is for c_prime and x
        scale = 1.0 / (b[i] - c_prime[i-1]*a[i])
        c_prime[i] *= scale
        x[i] = (x[i] - a[i] * x[i-1]) * scale
    end

    # Back-substitution
    for i = n-1:-1:1
        x[i] -= (c_prime[i] * x[i+1])
    end

    return x

end

function main()
    a = [0.0, 2.0, 3.0]
    b = [1.0, 3.0, 6.0]
    c = [4.0, 5.0, 0.0]
    d = [7.0, 5.0, 3.0]

    println(
        """The system
        $(join((b[1], c[1], "",   "|", d[1]), "\t"))
        $(join((a[2], b[2], c[2], "|", d[2]), "\t"))
        $(join(("",   a[3], b[3], "|", d[3]), "\t"))
        Has the solution:"""
    )

    soln = thomas(a, b, c, d, 3)

    println(soln)
end

main()


text:
function solve()
    k = parse(Int, readline())
    s = readline()

    if length(s) ≤ k
        s
    else
        s[1:k] * "..."
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
include("Include.jl")

function update_model_dictionary(parameter_guess_array,model_data_dictionary)

    # what is the host_type?
    host_type = :cell_free

     # Phase 1: parameter update =========================================================================== #
    # update the paramaters in the model data dictionary -
    # for now - lets only search over dG's -
    R = model_data_dictionary["R"]
    T_K = model_data_dictionary["T_K"]

    # compute W -
    tmp_W_array = Float64[]
    for index = 1:8
        parameter_guess = parameter_guess_array[index]
        value = exp(-1*parameter_guess/(R*T_K))
        push!(tmp_W_array,value)
    end

    # update the control W's -
    control_parameter_dictionary = model_data_dictionary["control_parameter_dictionary"]
    control_parameter_dictionary["W_clssrA_RNAP"] = tmp_W_array[1]          # 1
	control_parameter_dictionary["W_clssrA_sigma_28"] = tmp_W_array[2]      # 2
	control_parameter_dictionary["W_deGFPssrA_RNAP"] = tmp_W_array[3]       # 3
	control_parameter_dictionary["W_deGFPssrA_sigma_70"] = tmp_W_array[4]   # 4
	control_parameter_dictionary["W_deGFPssrA_clssrA"] = tmp_W_array[5]     # 5
	control_parameter_dictionary["W_sigma_28_RNAP"] = tmp_W_array[6]        # 6
	control_parameter_dictionary["W_sigma_28_sigma_70"] = tmp_W_array[7]    # 7
	control_parameter_dictionary["W_sigma_28_clssrA"] = tmp_W_array[8]      # 8
    model_data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary

    binding_parameter_dictionary = model_data_dictionary["binding_parameter_dictionary"]
    binding_parameter_dictionary["n_clssrA_sigma_28"] = parameter_guess_array[9]
	binding_parameter_dictionary["K_clssrA_sigma_28"] = parameter_guess_array[10]
	binding_parameter_dictionary["n_deGFPssrA_sigma_70"] = parameter_guess_array[11]
	binding_parameter_dictionary["K_deGFPssrA_sigma_70"] = parameter_guess_array[12]
	binding_parameter_dictionary["n_deGFPssrA_clssrA"] = parameter_guess_array[13]
	binding_parameter_dictionary["K_deGFPssrA_clssrA"] = parameter_guess_array[14]
	binding_parameter_dictionary["n_sigma_28_sigma_70"] = parameter_guess_array[15]
	binding_parameter_dictionary["K_sigma_28_sigma_70"] = parameter_guess_array[16]
	binding_parameter_dictionary["n_sigma_28_clssrA"] = parameter_guess_array[17]
    binding_parameter_dictionary["K_sigma_28_clssrA"] = parameter_guess_array[18]
    model_data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary

    # time constant modifier -
    time_constant_modifier_array = [
		0.0	                          ;	  # 1	clssrA
		0.0	                          ;	  # 2	deGFPssrA
		0.0	                          ;	  # 3	sigma_28
		0.0                           ;	  # 4	sigma_70
		parameter_guess_array[19]     ;	  # 5	mRNA_clssrA
		parameter_guess_array[20]     ;	  # 6	mRNA_deGFPssrA
		parameter_guess_array[21]     ;	  # 7	mRNA_sigma_28
		1.0 	                      ;	  # 8	mRNA_sigma_70
		parameter_guess_array[22]	  ;	  # 9	protein_clssrA
		parameter_guess_array[23]	  ;	  # 10	protein_deGFPssrA
		parameter_guess_array[24]	  ;	  # 11	protein_sigma_28
		1.0	                          ;	  # 12	protein_sigma_70
	]
    model_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

    # setup degradation_modifier_array -
    degradation_modifier_array = [
		0.0	                            ;	# 1	clssrA
		0.0	                            ;	# 2	deGFPssrA
		0.0	                            ;	# 3	sigma_28
		0.0	                            ;	# 4	sigma_70
		parameter_guess_array[25]	    ;	# 5	mRNA_clssrA 0.0085
		parameter_guess_array[26]		;	# 6	mRNA_deGFPssrA 0.05
		parameter_guess_array[27]		;	# 7	mRNA_sigma_28
		0.0                             ;   # 8 mRNA_sigma_70
		parameter_guess_array[28]		;	# 9	protein_clssrA
		parameter_guess_array[29]		;	# 10	protein_deGFPssrA
		parameter_guess_array[30]		;	# 11	protein_sigma_28
	    parameter_guess_array[31]		;	# 12	protein_sigma_70
    ]

    # update the translation time -
    model_data_dictionary["half_life_translation_capacity"] = parameter_guess_array[32]

    # lastly, update KL -
    biophysical_constants_dictionary = model_data_dictionary["biophysical_constants_dictionary"]
    biophysical_constants_dictionary["translation_saturation_constant"] = parameter_guess_array[33]
    model_data_dictionary["biophysical_constants_dictionary"] = biophysical_constants_dictionary

    # grab defaults -
    species_symbol_type_array = model_data_dictionary["species_symbol_type_array"]
    protein_coding_length_array = model_data_dictionary["protein_coding_length_array"]
    gene_coding_length_array = model_data_dictionary["gene_coding_length_array"]
    time_constant_modifier_array = model_data_dictionary["time_constant_modifier_array"]
    initial_condition_array = model_data_dictionary["initial_condition_array"]

    # # get gene IC -
    idx_gene = findall(x->x==:gene,species_symbol_type_array)
    gene_abundance_array = initial_condition_array[idx_gene]

    # Precompute the translation parameters -
	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
    model_data_dictionary["translation_parameter_array"] = translation_parameter_array

	# Precompute the kinetic limit of transcription -
	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
    model_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

    # Dilution degrdation matrix -
    dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)
    model_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix
    # ===================================================================================================== #

    # return -
    return model_data_dictionary
end

# Function to set some values for the weights, gene lengths etc -
function customize_data_dictionary(default_data_dictionary::Dict{String,Any}, host_type::Symbol)::Dict{String,Any}

    # make a deepcopy -
    customized_data_dictionary = deepcopy(default_data_dictionary)

    # set the gene lengths -
    gene_length_array = customized_data_dictionary["gene_coding_length_array"]
    gene_length_array[1] = 744.0    # clssrA
    gene_length_array[2] = 711.0    # deGFPssrA
    gene_length_array[3] = 720.0    # S28
    gene_length_array[4] = 720.0    # S70 => doesn't matter, gene expression = 0

    # set the protein lengths -
    protein_coding_length_array = customized_data_dictionary["protein_coding_length_array"]
    protein_coding_length_array[1] = 248.0  # clssrA
    protein_coding_length_array[2] = 237.0  # deGFPssrA
    protein_coding_length_array[3] = 240.0  # S28
    protein_coding_length_array[4] = 240.0  # S70 => doesn't matter, translation = 0

    # set ICs -
    initial_condition_array = customized_data_dictionary["initial_condition_array"]
    initial_condition_array[1] = 0.001  # clssrA gene
    initial_condition_array[2] = 0.008  # deGFP gene
    initial_condition_array[3] = 0.0015 # S28 gene
    initial_condition_array[4] = 0.0    # no S70 gene, just S70 protein in the extract
    #initial_condition_array[11] = 0.020 # sigma28 protein in muM
    initial_condition_array[12] = 0.035 # S70 protein in muM

    # setup the W's -
    control_parameter_dictionary = customized_data_dictionary["control_parameter_dictionary"]
    control_parameter_dictionary["W_clssrA_RNAP"] = 0.000014
	control_parameter_dictionary["W_clssrA_sigma_28"] = 1.0
	control_parameter_dictionary["W_deGFPssrA_RNAP"] = 0.000014
	control_parameter_dictionary["W_deGFPssrA_sigma_70"] = 10.0
	control_parameter_dictionary["W_deGFPssrA_clssrA"] = 100.0
	control_parameter_dictionary["W_sigma_28_RNAP"] = 0.00014
	control_parameter_dictionary["W_sigma_28_sigma_70"] = 100.0
	control_parameter_dictionary["W_sigma_28_clssrA"] = 0.1
	control_parameter_dictionary["W_sigma_70_RNAP"] = 0.0

    # setup the binding parameters -
    bP = readdlm("bP-Opt.dat")
    binding_parameter_dictionary = customized_data_dictionary["binding_parameter_dictionary"]
    bp_symbol_name_array = String[]
    push!(bp_symbol_name_array, "n_clssrA_sigma_28")
    push!(bp_symbol_name_array, "K_clssrA_sigma_28")
    push!(bp_symbol_name_array, "n_deGFPssrA_sigma_70")
    push!(bp_symbol_name_array, "K_deGFPssrA_sigma_70")
    push!(bp_symbol_name_array, "n_deGFPssrA_clssrA")
    push!(bp_symbol_name_array, "K_deGFPssrA_clssrA")
    push!(bp_symbol_name_array, "n_sigma_28_sigma_70")
    push!(bp_symbol_name_array, "K_sigma_28_sigma_70")
    push!(bp_symbol_name_array, "n_sigma_28_clssrA")
    push!(bp_symbol_name_array, "K_sigma_28_clssrA")
    for (index,key_text) in enumerate(bp_symbol_name_array)

        # get value -
        value = bP[index]

        # cache -
        binding_parameter_dictionary[key_text] = value
    end

    # setup degradation_modifier_array -
    degradation_modifier_array = [
		0.0	    ;	# 1	clssrA
		0.0	    ;	# 2	deGFPssrA
		0.0	    ;	# 3	sigma_28
		0.0	    ;	# 4	sigma_70
		0.05    ;	# 5	mRNA_clssrA 0.0085
		0.05	;	# 6	mRNA_deGFPssrA 0.05
		0.05	;	# 7	mRNA_sigma_28
		1.0	    ;	# 8	mRNA_sigma_70
		30.0	;	# 9	protein_clssrA
		28.0	;	# 10	protein_deGFPssrA
		1.5	    ;	# 11	protein_sigma_28
		1.0	    ;	# 12	protein_sigma_70
	]
    customized_data_dictionary["degradation_modifier_array"] = degradation_modifier_array
    biophysical_constants_dictionary = customized_data_dictionary["biophysical_constants_dictionary"]
    species_symbol_type_array = customized_data_dictionary["species_symbol_type_array"]
    protein_coding_length_array = customized_data_dictionary["protein_coding_length_array"]
    gene_coding_length_array = customized_data_dictionary["gene_coding_length_array"]

    # get gene IC -
    idx_gene = findall(x->x==:gene,species_symbol_type_array)
    gene_abundance_array = initial_condition_array[idx_gene]

    # Dilution degrdation matrix -
	dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)
    customized_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix

    # time constant modifier -
    time_constant_modifier_array = [
		0.0	                    ;	  # 1	clssrA
		0.0	                    ;	  # 2	deGFPssrA
		0.0	                    ;	  # 3	sigma_28
		0.0                     ;	  # 4	sigma_70
        0.5549555127597618      ;	  # 5	mRNA_clssrA
        0.9961201331095983      ;	  # 6	mRNA_deGFPssrA
		0.9573079006791378      ;	  # 7	mRNA_sigma_28
        2.47595196225463	    ;	  # 8	mRNA_sigma_70
		1.2340314042937286	    ;	  # 9	protein_clssrA
		2.0833904410438735 	    ;	  # 10	protein_deGFPssrA
		4.99999996459193	    ;	  # 11	protein_sigma_28
        2.47595196225463	    ;	  # 12	protein_sigma_70
	]
    customized_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

    # Precompute the translation parameters -
	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
    customized_data_dictionary["translation_parameter_array"] = translation_parameter_array

	# Precompute the kinetic limit of transcription -
	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
    customized_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

    # return -
    return customized_data_dictionary
end

function main(path_to_ensemble_file::String, path_to_sim_dir::String; sample_fraction::Float64=1.0)

    # load the default data_dictionary -
    time_start = 0.0
    time_stop = 16.0
    time_step_size = 0.01

    # what is the host_type?
    host_type = :cell_free

    # path to parameters -
    path_to_biophysical_constants_file = "./CellFree.json"

    # Load the data dictionary (uses the default biophysical_constants file)
    default_data_dictionary = build_data_dictionary((time_start,time_stop,time_step_size), path_to_biophysical_constants_file, host_type)

    # Set some values for the weights, gene lengths etc -
    custom_data_dictionary = customize_data_dictionary(default_data_dictionary,host_type)

    # load the ensemble file -
    parameter_ensemble_array = readdlm(path_to_ensemble_file)

    # sort the array -
    parameter_ensemble_array = sort_parameter_ensemble_array(parameter_ensemble_array)

    # what is the size of the ensemble -
    (number_of_parameters, number_of_samples) = size(parameter_ensemble_array)

    # take the top sample fraction -
    number_of_samples = round(Int,sample_fraction*number_of_samples)
    for sample_index = 1:number_of_samples

        # grab the parameter array -
        local_parameter_array = parameter_ensemble_array[:,sample_index]

        # update the default data_dictionary to reflect the new parameters -
        model_data_dictionary = update_model_dictionary(local_parameter_array, custom_data_dictionary)

        # solve the model equations -
        (T,X) = SolveBalances(time_start,time_stop,time_step_size,model_data_dictionary)

        # dump -
        data_array = [T X]
        filename = "$(path_to_sim_dir)/simulation-P$(sample_index).dat"
        writedlm(filename, data_array)

        # give user some notification -
        @show sample_index
    end
end

# setup paths -
path_to_sim_file = "$(pwd())/simulations"
path_to_ensemble_file = "$(pwd())/Ensemble-c1-restriction-T20.dat"

# call -
main(path_to_ensemble_file, path_to_sim_file; sample_fraction = 1.0)

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
function _make_dense_block(
    num_layers,
    bn_size,
    growth_rate,
    dropout,
    stage_index,
)
    result = Sequential()
    for _ = 1:num_layers
        add(result, _make_dense_layer(growth_rate, bn_size, dropout))
    end
    return result
end

function _make_dense_layer(growth_rate, bn_size, dropout)
    new_features = Sequential()
    add(new_features, BatchNorm(relu))
    add(new_features, Conv2D(bn_size * growth_rate, (1, 1), use_bias = false))
    add(new_features, BatchNorm(relu))
    add(
        new_features,
        Conv2D(growth_rate, (3, 3), padding = 1, use_bias = false),
    )
    if dropout > 0
        add(new_features, Dropout(dropout))
    end
    out = Concurrent()
    add(out, identity)
    add(out, new_features)
    return out
end

function _make_transition(num_output_features)
    out = Sequential()
    add(out, BatchNorm(relu))
    add(out, Conv2D(num_output_features, (1, 1), use_bias = false))
    add(out, AvgPool2D())
    return out
end

mutable struct DenseNet <: Photon.Layer
    features
    output
end

function DenseNet(
    num_init_features,
    growth_rate,
    block_config;
    bn_size = 4,
    dropout = 0,
    classes = 1000,
)

    features = Sequential()
    add(
        features,
        Conv2D(
            num_init_features,
            (7, 7),
            strides = 2,
            padding = 3,
            use_bias = false,
        ),
    )
    add(features, BatchNorm(relu))
    add(features, MaxPool2D(3, strides = 2, padding = 1))
    # Add dense blocks
    num_features = num_init_features
    for (i, num_layers) in enumerate(block_config)
        add(
            features,
            _make_dense_block(num_layers, bn_size, growth_rate, dropout, i + 1),
        )
        num_features = num_features + num_layers * growth_rate
        if i != length(block_config) - 1
            num_features = floor(Int, num_features / 2)
            add(features, _make_transition(num_features))
        end
    end
    add(features, BatchNorm(relu))
    add(features, AvgPool2D(7))
    add(features, Flatten())

    output = Dense(classes)
    DenseNet(features, output)
end

function (dn::DenseNet)(x)
    x = dn.features(x)
    dn.output(x)
end

densenet_spec = Dict(
    121 => (64, 32, [6, 12, 24, 16]),
    161 => (96, 48, [6, 12, 36, 24]),
    169 => (64, 32, [6, 12, 32, 32]),
    201 => (64, 32, [6, 12, 48, 32]),
)

DenseNet121(; bn_size = 4, dropout = 0, classes = 1000) =
    DenseNet(densenet_spec[121]..., bn_size = bn_size, dropout=dropout, classes=classes)

DenseNet161(; bn_size = 4, dropout = 0, classes = 1000) =
    DenseNet(densenet_spec[161]..., bn_size = bn_size, dropout=dropout, classes=classes)

DenseNet169(; bn_size = 4, dropout = 0, classes = 1000) =
    DenseNet(densenet_spec[169]..., bn_size = bn_size, dropout=dropout, classes=classes)

DenseNet201(; bn_size = 4, dropout = 0, classes = 1000) =
    DenseNet(densenet_spec[201]..., bn_size = bn_size, dropout=dropout, classes=classes)


text:
# Include some methods for processing the IMDB dataset. The dataset should be available at
#   https://www.imdb.com/interfaces/
# These methods assume that the data has been downloaded and extracted. You should pass the folder
# containing these files to the methods for processing them.
using CSV
using DataFrames

TITLE_BASICS_FILENAME = "title.basics.tsv"
NAME_BASICS_FILENAME = "name.basics.tsv"
CREDITS_FILENAME = "title.principals.tsv"

struct Imdb
    # film ID -> film Title
    films::Dict{String, String}

    # actor ID -> actor name
    actors::Dict{String, String}

    # actor ID -> cluster ID
    clusters::Dict{String, Int64}

    # film ID -> array of credited people's IDs
    credits::Dict{String, Array{String}}
end

"""
    loadimdb(directory)

Load the imdb dataset into memory, returning an Imdb object.
"""
function loadimdb(directory)
    # First, we will construct the films dictionary.
    #
    # Select the 'movies' from the titles file, and store the ID with the Movie title.
    filename = joinpath(directory, TITLE_BASICS_FILENAME)
    df = DataFrame(CSV.File(filename, select=[:tconst, :titleType, :primaryTitle]))
    moviesdf = df[df.titleType .== "movie", :]
    films = Dict(Pair.(moviesdf.tconst, moviesdf.primaryTitle))

    # Clear the dataframe to free up memory
    moviesdf = nothing

    # Now, by a similar method, we construct the actors dictionary.
    filename = joinpath(directory, NAME_BASICS_FILENAME)
    df = DataFrame(CSV.File(filename, select=[:nconst, :primaryName]))
    actors = Dict(Pair.(df.nconst, df.primaryName))

    # Now, construct the credits dictionary
    filename = joinpath(directory, CREDITS_FILENAME)
    df = DataFrame(CSV.File(filename, select=[:tconst, :nconst, :category]))

    # Select only the rows which correspond to films in our film dictionary.
    df = df[in(keys(films)).(df.tconst), :]

    # Select only the rows which correspond to actors, actresses and directors.
    df = df[in.(df.category, Ref(["actor", "actress", "director"])), :]

    # Grouping by movie, construct the credits dictionary.
    groups = groupby(df, :tconst)
    credits = Dict()
    clusters = Dict()
    for group in groups
        this_tconst = group[1, :tconst]
        this_actors = group[group.category .== "actor", :].nconst
        this_actresses = group[group.category .== "actress", :].nconst
        this_directors = group[group.category .== "director", :].nconst
        this_principals = []

        # Update the clusters dictionary
        for actor in this_actors
            clusters[actor] = 1
        end
        for actress in this_actresses
            clusters[actress] = 2
        end
        for director in this_directors
            clusters[director] = 3
        end

        # Keep at most the first k of each category
        k = 1
        append!(this_principals, this_actresses[1 : min(k, length(this_actresses))])
        append!(this_principals, this_actors[1 : min(k, length(this_actors))])
        append!(this_principals, this_directors[1 : min(k, length(this_directors))])
        credits[this_tconst] = this_principals
    end

    return Imdb(films, actors, clusters, credits)
end

"""
    writeCreditEdgelist(imdbdata, filename)
    writeCreditEdgelist(imdbDirectory)

Store the IMDB credit hyperaph in the edgelist format.

Create five files:
  - filename.edgelist
  - filename.edges
  - filename.vertices
  - filename.gt
  - filename.clusters

If only the imdb directory is specified, then the data is loaded, and then saved, with
filename="credit".

In the filename.edgelist file, there will be a single line for each edge in the hypergraph,
containing a space seperated list of nodes.

The filename.edges file will contain the name of the film corresponding to each edge.

The filename.vertices file will contain the name of the person corresponding to each vertex.
"""
function writeCreditEdgelist(imdbData::Imdb, filename::String)
    # Start by opening the ouput files
    edgelistFile = open("$filename.edgelist", "w")
    edgesFile = open("$filename.edges", "w")
    verticesFile = open("$filename.vertices", "w")
    gtFile = open("$filename.gt", "w")
    clusterFile = open("$filename.clusters", "w")

    # We will construct a dictionary of person IDs -> node IDs as we go along.
    person_node = Dict()

    vᵢ = 1
    for (movie_id, credit_list) in imdbData.credits
        # Add this edge to the edges file
        println(edgesFile, imdbData.films[movie_id])

        for person_id in credit_list
            # Add this person to the person_node dictionary if they are not there already, and
            # save to the vertices file.
            if person_id ∉ keys(person_node)
                person_node[person_id] = vᵢ
                vᵢ += 1

                # Write the name of the vertex
                if person_id ∈ keys(imdbData.actors)
                    println(verticesFile, imdbData.actors[person_id])
                else
                    println(verticesFile, "missing")
                end

                # Write the cluster of the vertex
                if person_id ∈ keys(imdbData.clusters)
                    println(gtFile, imdbData.clusters[person_id])
                else
                    println(verticesFile, "-1")
                end
            end

            print(edgelistFile, "$(person_node[person_id]) ")
        end

        # Add a newline at the end of each line
        print(edgelistFile, "\n")
    end

    # Write the cluster ID file
    print(clusterFile, "actor\nactress\ndirector")

    close(edgelistFile)
    close(edgesFile)
    close(verticesFile)
    close(gtFile)
    close(clusterFile)
end

function writeCreditEdgelist(imdbDirectory::String)
    imdbData = loadimdb(imdbDirectory)
    writeCreditEdgelist(imdbData, joinpath(imdbDirectory, "credit"))
end

max_stars_repo_name: pmacg/graph-data-processing

max_issues_repo_name: pmacg/graph-data-processing

max_forks_repo_name: pmacg/graph-data-processing


text:
function ueberpruefeMatrixVektora(A,m)

f = [2; 2; 2; 5; 1; 1; 6; 5; 0];

if norm(A*f-transpose(m)) <= 0.01
    # homogene Lösungen berechnen
    homogen = nullspace(A);

    # inhomogene Lösung berechnen
    inhomogen = pinv(A)*m';

    if size(homogen,2) > 0
        println("Euer Gleichungssystem stimmt zwar, allerdings betrachtet ihr noch nicht genug Strahlen.")
    else
        println("Ihr habt genug Strahlen verwendet. Die Lösung eures Gleichungssystems lautet:");
        print("f1 = ", round(inhomogen[1];digits = 2),"\n")
        print("f2 = ", round(inhomogen[2];digits = 2),"\n")
        print("f3 = ", round(inhomogen[3];digits = 2), "\n")
        print("f4 = ", round(inhomogen[4];digits = 2), "\n")
        print("f5 = ", round(inhomogen[5];digits = 2), "\n")
        print("f6 = ", round(inhomogen[6];digits = 2), "\n")
        print("f7 = ", round(inhomogen[7];digits = 2),"\n")
        print("f8 = ", round(inhomogen[8];digits = 2),"\n")
        print("f9 = ", round(inhomogen[9];digits = 2));
    end
else
    println("Schaut euch noch einmal den Verlauf der Strahlen an und überdenkt eure Matrix.")
end

end


text:
"""
    `rotstress(stress,theta)`

Rotates `stress` `theta` degrees. `stress` = [sigma1,sigma2,tau12]
"""
function rotstress!(stress::AbstractArray{<:Real,1},theta::Real)
  c = cosd(theta)
  s = sind(theta)
  c2 = c^2
  s2 = s^2
  cs = c*s

  rot11 = c2
  rot12 = s2
  rot13 = 2*cs
  rot21 = s2
  rot22 = c2
  rot23 = -2*cs
  rot31 = -cs
  rot32 = cs
  rot33 = c2-s2

  tmp1 = stress[1]
  tmp2 = stress[2]
  tmp3 = stress[3]

  stress[1] = rot11*tmp1 + rot12*tmp2 + rot13*tmp3
  stress[2] = rot21*tmp1 + rot22*tmp2 + rot23*tmp3
  stress[3] = rot31*tmp1 + rot32*tmp2 + rot33*tmp3

  return stress
end

function rotstress(stress::AbstractArray{<:Real,1}, theta::Real)
  tmp = copy(stress)
  return rotstress!(tmp, theta)
end

"""
    `rotstrain(strain,theta)`

Rotates `strain` `theta` degrees. `strain` = [eps1,eps2,gamma12]
"""
function rotstrain!(strain::AbstractArray{<:Real,1},theta::Real)
  c = cosd(theta)
  s = sind(theta)
  c2 = c^2
  s2 = s^2
  cs = c*s

  rot11 = c2
  rot12 = s2
  rot13 = cs
  rot21 = s2
  rot22 = c2
  rot23 = -cs
  rot31 = -2*cs
  rot32 = 2*cs
  rot33 = c2-s2

  tmp1 = strain[1]
  tmp2 = strain[2]
  tmp3 = strain[3]

  strain[1] = rot11*tmp1 + rot12*tmp2 + rot13*tmp3
  strain[2] = rot21*tmp1 + rot22*tmp2 + rot23*tmp3
  strain[3] = rot31*tmp1 + rot32*tmp2 + rot33*tmp3

  return strain
end

function rotstrain(strain::AbstractArray{<:Real,1},theta::Real)
  tmp = copy(strain)
  return rotstrain!(tmp, theta)
end

"""
    `getplystrain(nply::AbstractArray{<:Integer,1}, tply::AbstractArray{<:Real,1},
        theta::AbstractArray{<:Real,1}, resultantstrain::AbstractArray{<:Real,1})`

    `getplystrain(lam::Laminate, resultantstrain::AbstractArray{<:Real,1})`

Calculates strains in each ply aligned with principal material direction.
"""
function getplystrain(nply::AbstractArray{<:Integer,1}, tply::AbstractArray{<:Real,1},
    theta::AbstractArray{<:Real,1}, resultantstrain::AbstractArray{<:Real,1})

    # number of laminas
    nlam = length(nply)

    # coordinates of top and bottom of plies
    z = getz(tply, nply)

    # local strains at top and bottom of each lamina
    localstrain = [resultantstrain[1:3]+z[k]*resultantstrain[4:6] for k=1:nlam+1]

    # rotate to align with material principal axes
    lowerplystrain = [rotstrain(localstrain[k], theta[k]) for k = 1:nlam]
    upperplystrain = [rotstrain(localstrain[k+1],theta[k]) for k = 1:nlam]

  return lowerplystrain, upperplystrain
end
getplystrain(lam::Laminate, resultantstrain::AbstractArray{<:Real,1}) =
    getplystrain(lam.nply, lam.tply, lam.theta, resultantstrain)

"""
    `getplystress(plystrain::AbstractArray{<:AbstractArray{<:Real,1},1},
        q::AbstractArray{<:AbstractArray{<:Real,2},1},
        matid::AbstractArray{<:Integer,1})`

    `getplystress(plystrain::AbstractArray{<:AbstractArray{<:Real,1},1},
        q::AbstractArray{<:AbstractArray{<:Real,2},1}, lam::Laminate)`

Calculates ply stresses from ply strains and material stiffness matrices (`q`)
"""
function getplystress(plystrain::AbstractArray{<:AbstractArray{<:Real,1},1},
    q::AbstractArray{<:AbstractArray{<:Real,2},1},
    matid::AbstractArray{<:Integer,1})

    return [q[matid[i]]*plystrain[i] for i in 1:length(matid)]
end
getplystress(plystrain::AbstractArray{<:AbstractArray{<:Real,1},1},
    q::AbstractArray{<:AbstractArray{<:Real,2},1}, lam::Laminate) =
    getplystress(plystrain, q, lam.matid)


text:
"""
    gs(leaf::LeafModels{I,E,A,<:AbstractGsModel,S},gs_mod)
    gs(leaf::LeafModels{I,E,A,<:AbstractGsModel,S},meteo<:Atmosphere)
    gs!(leaf::LeafModels{I,E,A,<:AbstractGsModel,S},gs_mod)
    gs!(leaf::LeafModels{I,E,A,<:AbstractGsModel,S},meteo<:Atmosphere)

Default method to compute the stomatal conductance for CO₂ (mol m-2 s-1), it takes the form:

`leaf.stomatal_conductance.g0 + gs_closure(leaf,meteo) * leaf.status.A`

where gs_closure(leaf,meteo) computes the stomatal closure, and must be
implemented for the type of `leaf.stomatal_conductance`. The stomatal conductance is not
allowed to go below `leaf.stomatal_conductance.gs_min`.

# Arguments

- `leaf::LeafModels{I,E,A,<:AbstractGsModel,S}`: A leaf struct holding the parameters for the model. See
[`LeafModels`](@ref), and [`Medlyn`](@ref) or [`ConstantGs`](@ref) for the conductance models.
- `gs_mod`: the output from a [`gs_closure`](@ref) implementation (the conductance models
generally only implement this function)
- `meteo<:Atmosphere`: meteo data, see [`Atmosphere`](@ref)

# Examples

```julia
meteo = Atmosphere(T = 22.0, Wind = 0.8333, P = 101.325, Rh = 0.4490995)

# Using a constant value for Gs:

leaf = LeafModels(stomatal_conductance = Medlyn(0.03,12.0), # Instance of a Medlyn type
            A = 20.0, Cₛ = 380.0, Dₗ = meteo.VPD)

# Computing the stomatal conductance using the Medlyn et al. (2011) model:
gs(leaf,meteo)
```
"""
function gs(leaf::LeafModels{I,E,A,Gs,S},gs_mod) where {I,E,A,Gs<:AbstractGsModel,S}
    max(leaf.stomatal_conductance.gs_min,
        leaf.stomatal_conductance.g0 + gs_mod * leaf.status.A)
end

function gs(leaf::LeafModels{I,E,A,Gs,S},meteo::M) where {I,E,A,Gs<:AbstractGsModel,S,M<:Atmosphere}
    max(leaf.stomatal_conductance.gs_min,
        leaf.stomatal_conductance.g0 + gs_closure(leaf,meteo) * leaf.status.A)
end

function gs!(leaf::LeafModels{I,E,A,Gs,S},gs_mod) where {I,E,A,Gs<:AbstractGsModel,S}
    leaf.status.Gₛ = max(leaf.stomatal_conductance.gs_min,
        leaf.stomatal_conductance.g0 + gs_mod * leaf.status.A)
end

function gs!(leaf::LeafModels{I,E,A,Gs,S},meteo::M) where {I,E,A,Gs<:AbstractGsModel,S,M<:Atmosphere}
    leaf.status.Gₛ = max(leaf.stomatal_conductance.gs_min,
        leaf.stomatal_conductance.g0 + gs_closure(leaf,meteo) * leaf.status.A)
end


text:
"""
    dataset(name::String)::DataFrame

Load the dataset associated with `name`.

Some datasets used by Coles (2001) are available using the following names:
 - `portpirie`: annual maximum sea-levels in Port Pirie,
 - `glass`: breaking strengths of glass fibers
 - `fremantle`: annual maximum sea-levels in Fremantle
 - `rain`: daily rainfall accumulations in south-west England
 - `wooster`: daily minimum temperatures recorded in Wooster
 - `dowjones`: daily closing prices of the Dow Jones Index
 These datasets have been retrieved using the R package [*ismev*](https://cran.r-project.org/web/packages/ismev/index.html).

# Examples
```julia-repl
julia> Extremes.dataset("portpirie")
```
"""
function dataset(name::String)::DataFrame

    filename = joinpath(dirname(@__FILE__), "..", "data", string(name, ".csv"))
    if isfile(filename)
        return CSV.read(filename, DataFrame)
    end
    error("There is no dataset with the name '$name'")

end


text:
function get_input_reg(nbits::Int, samples::AbstractArray)
    get_input_reg(nbits, Vector(packbits(samples)))
end

function get_input_reg(nbits::Int, samples::Vector{<:Integer})
    nbatch = length(samples)
    config = zeros(ComplexF64, nbatch, 1<<nbits)
    for i = 1:nbatch
        config[i,Int(samples[i])+1] = 1
    end
    # transpose the register for better performance
    return arrayreg(transpose(config); nbatch=nbatch)
end

function VAN.free_energy(β::Real, H::AbstractBlock, sampler::AbstractSampler, circuit::AbstractBlock, samples)
    mean(free_energy_local(β, H, sampler, circuit, samples))
end

function energy(H::AbstractBlock, sampler, circuit::AbstractBlock, samples)
    reg = get_input_reg(nqubits(H), samples)
    mean(real.(expect(H, reg|>circuit)))
end

function energy2(H::AbstractBlock, sampler::AbstractSampler, circuit::AbstractBlock, samples)
    reg = get_input_reg(nqubits(H), samples)
    mean(real.(expect(H*H, reg|>circuit)))
end

function entropy(H::AbstractBlock, sampler::AbstractSampler, circuit::AbstractBlock, samples)
    -mean(get_logp(sampler, samples))
end

function purity(H::AbstractBlock, sampler::AbstractSampler, circuit::AbstractBlock, samples)
    mean(exp.(get_logp(sampler, samples)))
end

function spectra(H::AbstractBlock, sampler::AbstractSampler, circuit::AbstractBlock, samples)
    reg = get_input_reg(nqubits(H), samples)
    real.(expect(H, reg|>circuit))
end

function free_energy_local(β::Real, H::AbstractBlock, sampler, circuit::AbstractBlock, samples)
   reg = get_input_reg(nqubits(H), samples)
   return transpose(real.(expect(H, reg|>circuit))) + get_logp(sampler, samples)/β
end

# E(x) = <x|U'HU|x>
# d([E(x)]{x∼p} - TS)/dθ = [d(E(x))/dθ]{x∼p}
function grad_θ(H, circuit::AbstractBlock, samples)
    reg = get_input_reg(nqubits(H), samples)
    _, paramsδ = expect'(H, reg=>circuit)
    return paramsδ /size(samples)[end] # divide by batch size since there is a mean in loss
end

function loss_reinforce(β, H, sampler, circuit::AbstractBlock, samples)
    f = free_energy_local(β, H, sampler, circuit, samples)
    logp = get_logp(sampler, samples)
    b = mean(f)
    return mean(logp .* (f .- b))
end

function grad_sampler(β, H, sampler, circuit::AbstractBlock, samples)
    return Zygote.gradient(loss_reinforce, β, H, sampler, circuit, samples)[3]
end


text:
# module SearchLightEx

using SearchLight

struct SQLFunction <: SQLType
    name::Symbol
    fields::Tuple
end

struct SQLFunctionName <: SQLType
    name::Symbol
end
(f::SQLFunctionName)(args...) = SQLFunction(f.name, args)

function sql_functions(s)
    :(($(s...),) = $(map(x->SQLFunctionName(x), s)))
end

macro sql_functions(args...)
    esc(sql_functions(args))
end

# function AS(a::SQLFunction, b::Symbol)

# end

###### SearchLight ###############
function SearchLight.to_select_part(m::Type{T}, cols::Vector{Union{SQLColumn, SQLFunction}}, joins = SQLJoin[])::String where {T<:AbstractModel}
    SearchLight.Database.to_select_part(m, cols, joins)
end

function SearchLight.to_select_part(m::Type{T}, c::SQLFunction)::String where {T<:AbstractModel}
    SearchLight.to_select_part(m, [c])
end

########## SearchLight.Database ##############
function SearchLight.Database.prepare_column_name(fn_column::SQLFunction, _m::T) where {T<:AbstractModel}
    function process_sql_function(fn_column::SQLFunction)
        "$(fn_column.name)($(process_fields(fn_column.fields)))"
    end

    function process_fields(fields::Tuple)
        result = []

        for field in fields
            field_processed = if isa(field, SQLFunction)
                process_sql_function(field)
            elseif isa(field, Symbol)
                column_data::Dict{Symbol,Any} = SearchLight.from_literal_column_name(field |> String)

                if ! haskey(column_data, :table_name)
                    column_data[:table_name] = table_name(_m)
                end
                if ! haskey(column_data, :alias)
                    column_data[:alias] = ""
                end

                "$(SearchLight.to_fully_qualified(column_data[:column_name], column_data[:table_name]))"
            else
                "$field"
            end
            push!(result, field_processed)
        end

        """$(join(result, ","))"""
    end

    process_sql_function(fn_column)
end

function SearchLight.Database.to_select_part(m::Type{T}, cols::Vector{Union{SQLColumn, SQLFunction}}, joins = SQLJoin[])::String where {T<:AbstractModel}
    SearchLight.Database.DatabaseAdapter.to_select_part(m, cols, joins)
end

function SearchLight.Database._to_select_part(m::Type{T}, cols::Vector{Union{SQLColumn, SQLFunction}}, joins = SQLJoin[])::String where {T<:AbstractModel}
    _m::T = m()

    joined_tables = []

    if has_relation(_m, RELATION_HAS_ONE)
      rels = _m.has_one
      joined_tables = vcat(joined_tables, map(x -> is_lazy(x) ? nothing : (x.model_name)(), rels))
    end

    if has_relation(_m, RELATION_HAS_MANY)
      rels = _m.has_many
      joined_tables = vcat(joined_tables, map(x -> is_lazy(x) ? nothing : (x.model_name)(), rels))
    end

    if has_relation(_m, RELATION_BELONGS_TO)
      rels = _m.belongs_to
      joined_tables = vcat(joined_tables, map(x -> is_lazy(x) ? nothing : (x.model_name)(), rels))
    end

    filter!(x -> x != nothing, joined_tables)

    if ! isempty(cols)
      table_columns = []
      cols = vcat(cols, columns_from_joins(joins))

      for column in cols
        push!(table_columns, prepare_column_name(column, _m))
      end

      return join(table_columns, ", ")
    else
      table_columns = join(to_fully_qualified_sql_column_names(_m, persistable_fields(_m), escape_columns = true), ", ")
      table_columns = isempty(table_columns) ? String[] : vcat(table_columns, map(x -> prepare_column_name(x, _m), columns_from_joins(joins)))

      related_table_columns = String[]
      for rels in map(x -> to_fully_qualified_sql_column_names(x, persistable_fields(x), escape_columns = true), joined_tables)
        for col in rels
          push!(related_table_columns, col)
        end
      end

      return join([table_columns ; related_table_columns], ", ")
    end
end

####### SearchLight.Database.DatabaseAdapter ###########

# end

max_stars_repo_name: huaxk/rlservice

max_issues_repo_name: huaxk/rlservice

max_forks_repo_name: huaxk/rlservice


text:
using Pkg; Pkg.activate("..")
using ArgParse, Benchmarks, Serialization

function main(args)
    s = ArgParseSettings()
    @add_arg_table s begin
        "--nthreads"
            arg_type = Int
            required = true
            help = "Select the number of threads"
        "--datafile"
            arg_type = String
            required = true
            help = "File where the serialized data is held"
        "--refresh"
            action = :store_true
            help = "Start a new data collection"
    end

    parsed_args = parse_args(args, s)

    # Run the benchmark
    data = Benchmarks.benchmark(Kernel, threads = parsed_args["nthreads"])

    # If this is not a refresh, deserialize the datafile
    # Apparently "datafile" is stored as an array ...
    datafile = parsed_args["datafile"]
    if !parsed_args["refresh"]
        past_data = deserialize(datafile)
        push!(past_data, data)
        data_to_save = past_data
    else
        data_to_save = [data]
    end

    serialize(datafile, data_to_save)
    return 0
end

main(ARGS)

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
module TestCLIWrapper

using Test
using TSML

using TSML.CLIWrappers:tsmlrun,rawstat,aggregatedstat,aggregatedoutput,imputedstat,imputedoutput

const inputname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
const outputname = joinpath(dirname(pathof(TSML)),"../data/testdata_output.csv")

function test_cli()
  op1=tsmlrun(inputname,"")
  op2=imputedoutput(inputname,"")
  @test mean(op1.Value .== op2.Value) == 1.0
  op5=tsmlrun(inputname,"","dd/mm/yyyy HH:MM","stat")
  op6=imputedstat(inputname,"")
  @test (op5[1:1,:] .== op6[1:1,:]) |> Matrix |> sum == 20
end
@testset "CLIWrapper: check cliwrapper" begin
  test_cli()
end

end


text:
function Base.show(io::IO, x::EXPR, offset=0, d=0, er=false)
    T = typof(x)
    c =  T === ErrorToken || er ? :red : :normal
    # Print span as 1-based range of the source string. This presentation is
    # simple to understand when strings are presented to CSTParser.parse().
    print(io, lpad(offset + 1, 3), ":", rpad(offset + x.fullspan, 3), " ")
    if isidentifier(x)
        printstyled(io, " "^d, typof(x) == NONSTDIDENTIFIER ? valof(x.args[2]) : valof(x), color=:yellow)
        x.meta !== nothing && show(io, x.meta)
        println(io)
    elseif isoperator(x)
        printstyled(io, " "^d, "OP: ", kindof(x), "\n", color=c)
    elseif iskw(x)
        printstyled(io, " "^d, kindof(x), "\n", color=:magenta)
    elseif ispunctuation(x)
        if kindof(x) === Tokens.LPAREN
            printstyled(io, " "^d, "(\n", color=c)
        elseif kindof(x) === Tokens.RPAREN
            printstyled(io, " "^d, ")\n", color=c)
        elseif kindof(x) === Tokens.LSQUARE
            printstyled(io, " "^d, "[\n", color=c)
        elseif kindof(x) === Tokens.RSQUARE
            printstyled(io, " "^d, "]\n", color=c)
        elseif kindof(x) === Tokens.COMMA
            printstyled(io, " "^d, ",\n", color=c)
        else
            printstyled(io, " "^d, "PUNC: ", kindof(x), "\n", color=c)
        end
    elseif isliteral(x)
        printstyled(io, " "^d, "$(kindof(x)): ", valof(x), "\n", color=c)
    else
        printstyled(io, " "^d, T, color=c)
        if x.meta !== nothing
            print(io, "( ")
            show(io, x.meta)
            print(io, ")")
        end
        println(io)
        x.args === nothing && return
        for a in x.args
            show(io, a, offset, d + 1, er)
            offset += a.fullspan
        end
    end
end

struct CSTInfiniteLoop <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::CSTInfiniteLoop)
    print(io, ex.msg)
end


text:
module GraphModels
using LightGraphs

using SemanticModels.ModelTools
using SemanticModels.ModelTools.CategoryTheory
import SemanticModels.ModelTools.CategoryTheory: ⊔, FinSetMorph, dom, codom, func
import SemanticModels.ModelTools: model
export GraphModel, model, ⊔, GraphMorph, dom, codom, func, verify

struct GraphModel <: AbstractModel
  model::AbstractGraph
end

function model(::Type{GraphModel}, g::AbstractGraph)
    return GraphModel(g)
end

⊔(g::GraphModel, h::GraphModel) = GraphModel(blockdiag(g.model,h.model))

function (f::FinSetMorph)(gModel::G) where G <: GraphModel
    g = gModel.model
    dom(f) == vertices(g) || throw(DomainError(vertices(g), "dom(f) = $(dom(f)) but nv(g) = $(nv(g))"))
    ϕ = func(f)
    GraphModel(map(edges(g)) do e
        s,t = e.src, e.dst
        LightGraphs.Edge(ϕ(s), ϕ(t))
      end |> Graph)
end

"""    GraphMorph{T,F} <: Morph

morphisms in the category of Finite Graphs. The objects must be a subtype of AbstractGraph.

You can take a `FinSetMorph` and lift it to a graph homomorphism. This is the functor that
takes the finite set `1:n`, to the empty graph with `n` vertices.

"""
struct GraphMorph{T, F} <: AbstractMorph
    dom::T
    codom::T
    fun::F
end

"""    GraphMorph(g::AbstractGraph, f::FinSetMorph)

is defined to be the graph homomorphism you get by functorially lifting `f`.
That is, `f` acts on the vertex set of `g` as an `Int->Int` function, and then
must act on the edges consistently.
"""
GraphMorph(g::GraphModel, f::FinSetMorph) = GraphMorph(g, f(g), f)

dom(m::GraphMorph) = m.dom
codom(m::GraphMorph) = m.codom
func(m::GraphMorph) = begin
    f = func(m.fun)
    return i->f(i)
end

"""    verify(m::GraphMorph)

validate a graph homomorphism by checking that all the edges in `dom(m)` and map to edges in `codom(m)`.
"""
verify(m::GraphMorph) = begin
    dom(m.fun) == vertices(dom(m)) || return false
    codom(m.fun) == vertices(codom(m)) || return false
    E = Set(edges(codom(m)))
    f = func(m)
    map(edges(dom(m))) do e
        u,v = f(e.src), f(e.dst)
        if u > v
            u,v = v,u
        end
        LightGraphs.Edge(u,v) in E
    end |> all
end

end


text:
using Revise
using AlgebraicPetri
using OrdinaryDiffEq
using DiffEqFlux, Flux
using Catalyst
using ModelExploration

SIR =  LabelledPetriNet([:S, :I, :R],
                        :inf => ((:S, :I)=>(:I, :I)),
                        :rec => (:I=>:R));

initial_conc=[1e7,1e5,1e-10]
true_rates = [1e-8,1e-4]
data = generate_data(SIR, true_rates,initial_conc)

initial_states = initial_conc;
p = [-10.,-10.,-10.];
g = SIR;
g=SIR

prob = ODEProblem(ReactionSystem(g),
                    zeros(ns(g)), # initial state
                    (0.0,100.0),
                    zeros(ns(g)+nt(g)));

cur_p = exp.(p)
prob′ = remake(prob, p=cur_p, u0=initial_states, tspan=(0.0,100.0))
sol = solve(prob′, Tsit5())
sol_steps = [sol(i) for i in 1:100]
prediction = [sum_infected_states(g, s) for s in sol_steps]
sum(abs2, data .- prediction), sol

#x = train(SIR, data)

function test_train(g::LabelledPetriNet, data::Vector{Float64})
  param_guess = [1e-7,1e-5,1e7,1e5,1e-10] # correct [1e-8,1e-4,1e7,1e5,1e-10]
  prob = ODEProblem(ReactionSystem(g),
                    zeros(ns(g)), # initial state
                    (0.0,100.0),
                    zeros(ns(g)+nt(g))) # current variables in model

  l_func = make_loss(g, prob, data, states=Dict([1=>1e7,2=>1e5,3=>1e-10]),
                      rates=Dict([2=>1e-5]))

  p = DiffEqFlux.sciml_train(l_func,log.(param_guess),ADAM(0.1),maxiters = 1000)
  return p, l_func
end

x, l = test_train(SIR, data)