text:
# GOAL: refactor the code in such a way so that verifier could be written outsied the main module.

# verify(proposition, proof, verifier) -> Bool
# verify(simulator) -> Bool

# The proposition of course can be verified knowing the secret (secret key or randomization factors) alone:
# verify(proposition, secret) -> Bool

# The challnege is last step obtained executing finite state machine. Does not need to contain responses.
# verify(proposition, proof, challenge) -> Bool

# prove(proposition, secret, verifier) -> Simulator

# shuffle(𝔀, g, pk) -> Tuple{Proposition, Secret}
# shuffle(𝔀, g, pk, verifier) -> Simulator

# decrypt(A, g, sk) -> Proposition # secret key is the one which is possible to be used
# decrypt(A, g, sk, verifier) -> Simulator

# Derives a public key for a given secret key
# Proving knowldge of sk such that: x = g^sk
# power(g, sk) -> Proposition
# power(g, sk, verifier) -> Simulator

abstract type Proposition end
abstract type Proof end
abstract type Verifier end

function step end
function challenge end

struct Simulator
    proposition::Proposition # Proposition type is the one which selects a type of prover being used
    proof::Proof
    verifier::Verifier
end

struct Shuffle <: Proposition
    g
    pk
    𝐞
    𝐞′
end

struct ShuffleSecret
    𝛙
    𝐫′
end

struct PoSProof <: Proof
    𝐜
    𝐜̂
    t
    s
end

struct PoSChallenge
    𝐡 # Independent set of generators
    𝐮 # PoS commitment challenge
    c # Last bit of a challenge
end

function verify(proposition::Shuffle, secret::ShuffleSecret)
    error("TODO")
end

function gen_shuffle(enc::Enc, e::ElGamal, r::Vector{T}) where T <: Integer

    e_enc = enc(e, r)
    ψ = sortperm(e_enc)

    sort!(e_enc)

    (; g, pk) = enc

    proposition = Shuffle(g, pk, e, e_enc)
    secret = ShuffleSecret(ψ, r)

    return proposition, secret
end

function gen_commitment(g::G, 𝐡::Vector{G}, b::Vector, r::Integer) where G <: Generator

    #(; g, h) = crs
    com = g^r * prod(𝐡 .^ b)

    return com
end

function gen_perm_commitment(g::G, 𝐡::Vector{G}, 𝛙::Vector, 𝐫::Vector) where G <: Generator

    #(; g, 𝐡) = crs

    commitments = [g^𝐫[j] * 𝐡[i] for (i, j) in enumerate(𝛙)]
    sorted_commitments = commitments[𝛙]

    return sorted_commitments
end

function gen_commitment_chain(g::Generator, c0::T, 𝐮::Vector, 𝐫::Vector) where T

    N = length(𝐮)

    𝐜 = Vector{T}(undef, N)

    𝐜[1] = g^𝐫[1] * c0^𝐮[1]

    for i in 2:N
        𝐜[i] = g^𝐫[i] * 𝐜[i-1]^𝐮[i]
    end

    return 𝐜
end

∑(𝐱, q) = mod(sum(𝐱), q) ### Need to improve
∏(𝐱) = prod(𝐱)
∏(f, 𝐱) = prod(f, 𝐱)

_a(x) = b(x)
_b(x) = a(x)

using Random: default_rng, rand

function prove(proposition::Shuffle, secret::ShuffleSecret, verifier::Verifier; rng = default_rng())

    (; 𝛙, 𝐫′) = secret
    (; g, pk, 𝐞, 𝐞′) = proposition

    v1 = step(verifier, proposition) # So I could keep a proposition in the coresponding state machine in the end
    𝐡, h = challenge(v1)

    @assert length(𝛙) == length(𝐞)

    N = length(𝛙)
    q = order(g)

    𝐫 = rand(rng, 2:q-1, N)
    𝐫̂ = rand(rng, 2:q-1, N)
    𝛚 = rand(rng, 2:q-1, 4)
    𝛚̂ = rand(rng, 2:q-1, N)
    𝛚̂′ = rand(rng, 2:q-1, N)

    𝐚′ = _a(𝐞′)
    𝐛′ = _b(𝐞′)

    𝐜 = gen_perm_commitment(g, 𝐡, 𝛙, 𝐫)

    v2 = step(v1, 𝐜)
    𝐮 = challenge(v2)

    𝐮′ = 𝐮[𝛙]

    𝐜̂ = gen_commitment_chain(g, h, 𝐮′, 𝐫̂)

    𝐯 = Vector(undef, N)
    𝐯[N] = 1
    for i in N-1:-1:1
        𝐯[i] = 𝐮′[i+1] * 𝐯[i+1]
    end

    r̄ = ∑(𝐫, q)
    r̂ = ∑(𝐫̂ .* 𝐯, q)
    r̃ = ∑(𝐫 .* 𝐮, q)
    r′ = ∑(𝐫′ .* 𝐮, q)

    t₁ = g^𝛚[1]
    t₂ = g^𝛚[2]

    t₃ = g^𝛚[3] * ∏(𝐡 .^ 𝛚̂′)

    t₄₁ = pk^(-𝛚[4]) * ∏(𝐚′ .^ 𝛚̂′)
    t₄₂ = g^(-𝛚[4]) * ∏(𝐛′ .^ 𝛚̂′)

    𝐭̂ = Vector(undef, N)
    𝐭̂[1] = g^𝛚̂[1] * h^𝛚̂′[1]
    for i in 2:N
        𝐭̂[i] = g^𝛚̂[i] * 𝐜̂[i-1]^𝛚̂′[i]
    end

    y = (𝐞, 𝐞′, 𝐜, 𝐜̂, pk)
    t = (t₁, t₂, t₃, (t₄₁, t₄₂), 𝐭̂)

    v3 = step(v2, 𝐜̂, t)
    c = challenge(v3)

    s₁ = mod(𝛚[1] + c * r̄, q)
    s₂ = mod(𝛚[2] + c * r̂, q)
    s₃ = mod(𝛚[3] + c * r̃, q)
    s₄ = mod(𝛚[4] + c * r′, q)

    𝐬̂ = mod.(𝛚̂ .+ c .* 𝐫̂, q) ### What can I do if I have a 0 as one of the elements?
    𝐬′ = mod.(𝛚̂′ .+ c .* 𝐮′, q)

    s = (s₁, s₂, s₃, s₄, 𝐬̂, 𝐬′) # Do I need to ensure that `s` are without 0 elements

    proof = PoSProof(𝐜, 𝐜̂, t, s)

    simulator = Simulator(proposition, proof, verifier)

    return simulator
end

function verify(proposition::Shuffle, proof::PoSProof, verifier::Verifier)

    v1 = step(verifier, proposition)

    (; 𝐜) = proof
    v2 = step(v1, 𝐜)

    (; 𝐜̂, t) = proof
    v3 = step(v2, 𝐜̂, t)

    chg = PoSChallenge(v3)
    return verify(proposition, proof, chg)
end

function verify(proposition::Shuffle, proof::PoSProof, challenge::PoSChallenge)

    (; g, pk, 𝐞, 𝐞′) = proposition
    (; 𝐜, 𝐜̂, t, s) = proof
    (; 𝐡, 𝐮, c) = challenge

    #(; g, pk, 𝐞, 𝐞′, 𝐡, 𝐜, 𝐮, 𝐜̂, t, c, s) = v

    (s₁, s₂, s₃, s₄, 𝐬̂, 𝐬′) = s # 𝐬̂ = 𝐤_B, 𝐬′ = 𝐤_E, s₃ = k_A
    (t₁, t₂, t₃, (t₄₁, t₄₂), 𝐭̂) = t # 𝐭̂ = 𝐁′, t₃ = A′, t₄ = F

    𝐚 = _a(𝐞)
    𝐛 = _b(𝐞)
    𝐚′ = _a(𝐞′)
    𝐛′ = _b(𝐞′)
    h = 𝐡[1]

    q = order(g)
    N = length(𝐞)

    c̄ = ∏(𝐜) / ∏(𝐡)
    u = mod(∏(𝐮), q)

    ĉ = 𝐜̂[N] / h^u
    c̃ = ∏(𝐜 .^ 𝐮)

    a′ = ∏(𝐚 .^ 𝐮)
    b′ = ∏(𝐛 .^ 𝐮)

    t₁′ = c̄^(-c) * g^s₁
    t₂′ = ĉ^(-c) * g^s₂
    t₃′ = c̃^(-c) * g^s₃ * ∏(𝐡 .^ 𝐬′)

    t₄₁′ = a′^(-c) * pk^(-s₄) * ∏(𝐚′ .^ 𝐬′)
    t₄₂′ = b′^(-c) * g^(-s₄) * ∏(𝐛′ .^ 𝐬′)

    𝐭̂′ = Vector(undef, N)

    𝐭̂′[1] = 𝐜̂[1]^(-c) * g^𝐬̂[1] * h^𝐬′[1]    #ĉ0 = h

    for i in 2:N
        𝐭̂′[i] = 𝐜̂[i]^(-c) * g^𝐬̂[i] * 𝐜̂[i-1]^𝐬′[i]
    end

    @show t₁ == t₁′
    @show t₂ == t₂′
    @show t₃ == t₃′
    @show t₄₁ == t₄₁′
    @show t₄₂ == t₄₂′

    for i in 1:N
        @show 𝐭̂[i] == 𝐭̂′[i]
    end

    return true # ToDo
end

verify(simulator::Simulator) = verify(simulator.proposition, simulator.proof, simulator.verifier)


text:
abstract type AbstractBikeSharing <: Name end
struct BikeSharingHour <: AbstractBikeSharing end
struct BikeSharingDay <: AbstractBikeSharing end

# mandatory methods
format(::Type{<:AbstractBikeSharing}) = TabularData
task(::Type{<:AbstractBikeSharing}) = Regression
source(::Type{<:AbstractBikeSharing}) = "https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset"

function downloadlink(N::Type{<:AbstractBikeSharing})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/00275/Bike-Sharing-Dataset.zip"
end

function preprocess(N::Type{BikeSharingHour})
    istarget = 17
    iscategorical = vcat(3:10, 15:16)
    toremove = 1

    return function (path)
        unzip(path; skip = f -> f.name in ["day.csv", "Readme.txt"])
        file = joinpath(datadir(N), "hour.csv")
        preprocess_csv(N, file, :train; istarget, iscategorical, toremove, header = true)
        return
    end
end

function preprocess(N::Type{BikeSharingDay})
    istarget = 16
    iscategorical = vcat(3:9, 14:15)
    toremove = 1

    return function (path)
        unzip(path; skip = f -> f.name in ["hour.csv", "Readme.txt"])
        file = joinpath(datadir(N), "day.csv")
        preprocess_csv(N, file, :train; istarget, iscategorical, toremove, header = true)
        return
    end
end

nattributes(::Type{BikeSharingHour}) = (16, )
nattributes(::Type{BikeSharingDay}) = (15, )
ninstances(::Type{BikeSharingHour}) = (17389, 0, 0)
ninstances(::Type{BikeSharingDay}) = (731, 0, 0)

# optional methods
function checksum(::Type{<:AbstractBikeSharing})
    return "b70182d0d0508e9abbb79306ce5c0cec34869000f8220175ac83d11dbe845401"
end

name(::Type{<:AbstractBikeSharing}) = "BikeSharing Compressive Strength Data Set"
author(::Type{<:AbstractBikeSharing}) = ["Hadi Fanaee-T", "Joao Gama"]

function citation(::Type{<:AbstractBikeSharing})
    return """
    @article{
        year={2013},
        issn={2192-6352},
        journal={Progress in Artificial Intelligence},
        doi={10.1007/s13748-013-0040-3},
        title={Event labeling combining ensemble detectors and background knowledge},
        url={[Web Link]},
        publisher={Springer Berlin Heidelberg},
        keywords={Event labeling; Event detection; Ensemble learning; Background knowledge},
        author={Fanaee-T, Hadi and Gama, Joao},
        pages={1-15}
    }
    """
end

ext: jl

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: identitystore
using AWS.Compat
using AWS.UUIDs

"""
    DescribeGroup()

Retrieves the group metadata and attributes from GroupId in an identity store.

# Required Parameters
- `GroupId`: The identifier for a group in the identity store.
- `IdentityStoreId`: The globally unique identifier for the identity store, such as d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated string which contains number and lower case letters. This value is generated at the time that a new identity store is created.

"""
describe_group(GroupId, IdentityStoreId; aws_config::AWSConfig=global_aws_config()) = identitystore("DescribeGroup", Dict{String, Any}("GroupId"=>GroupId, "IdentityStoreId"=>IdentityStoreId); aws_config=aws_config)
describe_group(GroupId, IdentityStoreId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = identitystore("DescribeGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GroupId"=>GroupId, "IdentityStoreId"=>IdentityStoreId), args)); aws_config=aws_config)

"""
    DescribeUser()

Retrieves the user metadata and attributes from UserId in an identity store.

# Required Parameters
- `IdentityStoreId`: The globally unique identifier for the identity store, such as d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated string which contains number and lower case letters. This value is generated at the time that a new identity store is created.
- `UserId`: The identifier for a user in the identity store.

"""
describe_user(IdentityStoreId, UserId; aws_config::AWSConfig=global_aws_config()) = identitystore("DescribeUser", Dict{String, Any}("IdentityStoreId"=>IdentityStoreId, "UserId"=>UserId); aws_config=aws_config)
describe_user(IdentityStoreId, UserId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = identitystore("DescribeUser", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IdentityStoreId"=>IdentityStoreId, "UserId"=>UserId), args)); aws_config=aws_config)

"""
    ListGroups()

Lists the attribute name and value of the group that you specified in the search. We only support DisplayName as a valid filter attribute path currently, and filter is required. This API returns minimum attributes, including GroupId and group DisplayName in the response.

# Required Parameters
- `IdentityStoreId`: The globally unique identifier for the identity store, such as d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated string which contains number and lower case letters. This value is generated at the time that a new identity store is created.

# Optional Parameters
- `Filters`: A list of Filter objects, which is used in the ListUsers and ListGroups request.
- `MaxResults`: The maximum number of results to be returned per request, which is used in the ListUsers and ListGroups request to specify how many results to return in one page. The length limit is 50 characters.
- `NextToken`: The pagination token used for the ListUsers and ListGroups APIs. This value is generated by the identity store service and is returned in the API response if the total results are more than the size of one page, and when this token is used in the API request to search for the next page.
"""
list_groups(IdentityStoreId; aws_config::AWSConfig=global_aws_config()) = identitystore("ListGroups", Dict{String, Any}("IdentityStoreId"=>IdentityStoreId); aws_config=aws_config)
list_groups(IdentityStoreId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = identitystore("ListGroups", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IdentityStoreId"=>IdentityStoreId), args)); aws_config=aws_config)

"""
    ListUsers()

Lists the attribute name and value of the user that you specified in the search. We only support UserName as a valid filter attribute path currently, and filter is required. This API returns minimum attributes, including UserId and UserName in the response.

# Required Parameters
- `IdentityStoreId`: The globally unique identifier for the identity store, such as d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated string which contains number and lower case letters. This value is generated at the time that a new identity store is created.

# Optional Parameters
- `Filters`: A list of Filter objects, which is used in the ListUsers and ListGroups request.
- `MaxResults`: The maximum number of results to be returned per request, which is used in the ListUsers and ListGroups request to specify how many results to return in one page. The length limit is 50 characters.
- `NextToken`: The pagination token used for the ListUsers and ListGroups APIs. This value is generated by the identity store service and is returned in the API response if the total results are more than the size of one page, and when this token is used in the API request to search for the next page.
"""
list_users(IdentityStoreId; aws_config::AWSConfig=global_aws_config()) = identitystore("ListUsers", Dict{String, Any}("IdentityStoreId"=>IdentityStoreId); aws_config=aws_config)
list_users(IdentityStoreId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = identitystore("ListUsers", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IdentityStoreId"=>IdentityStoreId), args)); aws_config=aws_config)


text:
type CoupledSchroedinger1D <: TSSM.TimeSplittingSpectralMethodComplex1D
    m1::Schroedinger1D
    m2::Schroedinger1D
    delta
    e
    function CoupledSchroedinger1D(nx::Integer, xmin::Real, xmax::Real,
                        delta::Real, e::Real)
       m1 = Schroedinger1D(nx, xmin, xmax)
       m2 = Schroedinger1D(nx, xmin, xmax)
       ev1 = get_eigenvalues(m1, true)
       ev2 = get_eigenvalues(m2, true)
       f = 2*2*pi/(xmax-xmin)
       ev1[1:end] = ev1 - f*delta*[0:nx/2-1;-nx/2:-1]
       ev2[1:end] = ev2 + f*delta*[0:nx/2-1;-nx/2:-1]
       new(m1, m2, delta, e)
    end
end

type WfCoupledSchroedinger1D <: TSSM.WaveFunctionComplex1D
   psi1::WfSchroedinger1D
   psi2::WfSchroedinger1D
   m::CoupledSchroedinger1D
   function WfCoupledSchroedinger1D( m::CoupledSchroedinger1D )
      new( WfSchroedinger1D(m.m1), WfSchroedinger1D(m.m2), m )
   end
end

function TSSM.wave_function(m::CoupledSchroedinger1D )
    WfCoupledSchroedinger1D(m)
end

function TSSM.propagate_A!(psi::WfCoupledSchroedinger1D, dt::Number)
    propagate_A!(psi.psi1, dt)
    propagate_A!(psi.psi2, dt)
end

function TSSM.copy!(target::WfCoupledSchroedinger1D, source::WfCoupledSchroedinger1D)
    copy!(target.psi1, source.psi1)
    copy!(target.psi2, source.psi2)
end

function TSSM.scale!(psi::WfCoupledSchroedinger1D, factor::Number)
    scale!(psi.psi1, factor)
    scale!(psi.psi2, factor)
end

function TSSM.axpy!(psi1::WfCoupledSchroedinger1D, psi2::WfCoupledSchroedinger1D, factor::Number)
    axpy!(psi1.psi1, psi2.psi1, factor)
    axpy!(psi1.psi2, psi2.psi2, factor)
end

type Function2
    f1::Function
    f2::Function
end

function TSSM.set!(psi::WfCoupledSchroedinger1D, f::Function2)
    set!(psi.psi1, f.f1)
    set!(psi.psi2, f.f2)
end

function TSSM.set!(psi::WfCoupledSchroedinger1D, f::Function2, t::Real)
    set!(psi.psi1, f.f1, t)
    set!(psi.psi2, f.f2, t)
end

function TSSM.norm(psi::WfCoupledSchroedinger1D )
   sqrt(norm(psi.psi1)^2+norm(psi.psi2)^2)
end

function TSSM.distance(psi1::WfCoupledSchroedinger1D, psi2::WfCoupledSchroedinger1D)
   sqrt(distance(psi1.psi1, psi2.psi1)^2+distance(psi1.psi2,psi2.psi2)^2)
end

function TSSM.propagate_B!(psi::WfCoupledSchroedinger1D, dt::Number)
   u1 = get_data(psi.psi1, true)
   u2 = get_data(psi.psi2, true)
   to_real_space!(psi.psi1)
   to_real_space!(psi.psi2)
   u1[1:end] = exp(1im*dt*(abs(u1).^2 + psi.m.e*abs(u2).^2)).*u1
   u2[1:end] = exp(1im*dt*(psi.m.e*abs(u1).^2 + abs(u2).^2)).*u2
end


text:
# using Revise

using Test

using IncrementalInference, Distributions
using DistributedFactorGraphs
using Statistics
# going to introduce two new constraint types
import Base: convert
import IncrementalInference: getSample

mutable struct DevelopPrior <: FunctorSingleton
  x::Distribution
end
getSample(dpl::DevelopPrior, N::Int=1) = (reshape(rand(dpl.x, N),1,N), )

mutable struct DevelopLikelihood <: FunctorPairwise
  x::Distribution
end
getSample(dpl::DevelopLikelihood, N::Int=1) = (reshape(rand(dpl.x, N),1,N), )
function (vv::DevelopLikelihood)(res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple,
            wXi::Array{Float64,2},
            wXj::Array{Float64,2}  )::Nothing
  #
  res[1] = meas[1][idx] - (wXj[1,idx] - wXi[1,idx])
  nothing
end

global N  = 100
global fg = initfg()

@testset "test populate factor graph with a multi-hypothesis factor..." begin

global v1 = addVariable!(fg, :x1, ContinuousScalar, N=N)

global pr = DevelopPrior(Normal(10.0,1.0))
global f1 = addFactor!(fg,[:x1],pr)

ensureAllInitialized!(fg)

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl",121)

global pts = evalFactor2(fg, f1, v1.label, N=N)

@test sum(abs.(pts .- 1.0) .< 5) < 30
@test sum(abs.(pts .- 10.0) .< 5) > 30

global v2 = addVariable!(fg, :x2, ContinuousScalar, N=N)
global pp = DevelopLikelihood(Normal(100.0,1.0))
global f2 = addFactor!(fg, [:x1;:x2], pp)

ensureAllInitialized!(fg)

@test abs(Statistics.mean(getVal(fg, :x2))-110.0) < 10.0

global v3 = addVariable!(fg, :x3, ContinuousScalar, N=N)
global v4 = addVariable!(fg, :x4, ContinuousScalar, N=N)

global ppMH = DevelopLikelihood(Normal(90.0,1.0))
global f3 = addFactor!(fg, [:x2;:x3;:x4], ppMH, multihypo=[1.0;0.5;0.5])

# @test getData(f3).fnc.hypoverts == [:x3, :x4]
@test sum(abs.(solverData(f3).fnc.hypotheses.p[1] .- 0.0)) < 0.1  # 1.0 becomes 0.0 for computational convenience
@test sum(abs.(solverData(f3).fnc.hypotheses.p[2:3] .- 0.5)) < 0.1

manualinit!(fg, :x2, 1*ones(1,100))
manualinit!(fg, :x3, 2*ones(1,100))
manualinit!(fg, :x4, 3*ones(1,100))

end

@testset "Test multi-hypothesis factor convolution exploration" begin

global pts = approxConv(fg, :x2x3x4f1, :x2, N=N)

@test 99 < sum(pts .<= -70.0)

global pts = approxConv(fg, :x2x3x4f1, :x3, N=N)

@test 25 < sum(pts .== 3.0) < 75

global pts = approxConv(fg, :x2x3x4f1, :x4, N=N)

@test 25 < sum(pts .== 2.0) < 75

end

println("Packing converters")

mutable struct PackedDevelopPrior <: PackedInferenceType
  x::String
  PackedDevelopPrior() = new()
  PackedDevelopPrior(x) = new(x)
end
function convert(::Type{PackedDevelopPrior}, d::DevelopPrior)
  PackedDevelopPrior(string(d.x))
end
function convert(::Type{DevelopPrior}, d::PackedDevelopPrior)
  DevelopPrior(IncrementalInference.normalfromstring(d.x))
end

mutable struct PackedDevelopLikelihood <: PackedInferenceType
  x::String
  PackedDevelopLikelihood() = new()
  PackedDevelopLikelihood(x) = new(x)
end
function convert(::Type{PackedDevelopLikelihood}, d::DevelopLikelihood)
  PackedDevelopLikelihood(string(d.x))
end
function convert(::Type{DevelopLikelihood}, d::PackedDevelopLikelihood)
  DevelopLikelihood(IncrementalInference.extractdistribution(d.x))
end

@testset "test packing and unpacking the data structure" begin

  global topack = solverData(f1)
  global dd = convert(PackedFunctionNodeData{PackedDevelopPrior},topack)
  global unpacked = convert(FunctionNodeData{CommonConvWrapper{DevelopPrior}},dd)

  @test abs(unpacked.fnc.usrfnc!.x.μ - 10.0) < 1e-10
  @test abs(unpacked.fnc.usrfnc!.x.σ - 1.0) < 1e-10

  global topack = solverData(f3)
  global dd = convert(PackedFunctionNodeData{PackedDevelopLikelihood},topack)
  global unpacked = convert(FunctionNodeData{CommonConvWrapper{DevelopLikelihood}},dd)

  # @test unpacked.fnc.hypoverts == Symbol[:x3; :x4]
  @test sum(abs.(unpacked.fnc.hypotheses.p[1] .- 0.0)) < 0.1
  @test sum(abs.(unpacked.fnc.hypotheses.p[2:3] .- 0.5)) < 0.1
  # str = "Symbol[:x3, :x4];[0.5, 0.5]"
  # IncrementalInference.parsemultihypostr(str)

end

# start a new factor graph
global N = 200
global fg = initfg()

@testset "test tri-modal factor..." begin

global v1 = addVariable!(fg, :x1, ContinuousScalar, N=N)

global pr = DevelopPrior(Normal(10.0,1.0))
global f1 = addFactor!(fg,[:x1],pr)

ensureAllInitialized!(fg)

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl",121)

global pts = approxConv(fg, Symbol(f1.label), :x1, N=N)

@test sum(abs.(pts .- 1.0) .< 5) < 30
@test sum(abs.(pts .- 10.0) .< 5) > 30

global v2 = addVariable!(fg, :x2, ContinuousScalar, N=N)
global pp = DevelopLikelihood(Normal(100.0,1.0))
global f2 = addFactor!(fg, [:x1;:x2], pp)

ensureAllInitialized!(fg)

@test abs(Statistics.mean(getVal(fg, :x2))-110.0) < 10.0

global v3 = addVariable!(fg, :x3, ContinuousScalar, N=N)
global v4 = addVariable!(fg, :x4, ContinuousScalar, N=N)
global v5 = addVariable!(fg, :x5, ContinuousScalar, N=N)

global ppMH = DevelopLikelihood(Normal(90.0,1.0))
global f3 = addFactor!(fg, [:x2;:x3;:x4;:x5], ppMH, multihypo=(1.0,0.333,0.333,0.334))

# @test getData(f3).fnc.hypoverts == [:x3, :x4]
@test sum(abs.(solverData(f3).fnc.hypotheses.p[1] .- 0.0)) < 0.1  # 1.0 becomes 0.0 for computational convenience
@test sum(abs.(solverData(f3).fnc.hypotheses.p[2] .- 0.333)) < 0.001
@test sum(abs.(solverData(f3).fnc.hypotheses.p[3] .- 0.333)) < 0.001
@test sum(abs.(solverData(f3).fnc.hypotheses.p[4] .- 0.334)) < 0.001

manualinit!(fg, :x2 ,1*ones(1,100))
manualinit!(fg, :x3 ,2*ones(1,100))
manualinit!(fg, :x4 ,3*ones(1,100))
manualinit!(fg, :x5 ,4*ones(1,100))

# solve for certain idx
global pts = approxConv(fg, :x2x3x4x5f1, :x2, N=N)

@test 0.95*N < sum(pts .<= -70.0)

# solve for one of uncertain variables
global pts = approxConv(fg, :x2x3x4x5f1, :x3, N=N)

@test 0.1*N < sum(80 .< pts .< 100.0) < 0.5*N
@test 0.1*N < sum(pts .== 3.0) < 0.5*N
@test 0.1*N < sum(pts .== 4.0) < 0.5*N

@test sum(70 .< pts .< 110.0) + sum(pts .== 3.0) + sum(pts .== 4.0) == N

# solve for one of uncertain variables
global pts = approxConv(fg, :x2x3x4x5f1, :x4, N=N)

@test 0.1*N < sum(80 .< pts .< 100.0) < 0.5*N
@test 0.1*N < sum(pts .== 2.0) < 0.5*N
@test 0.1*N < sum(pts .== 4.0) < 0.5*N

@test sum(80 .< pts .< 100.0) + sum(pts .== 2.0) + sum(pts .== 4.0) == N

# solve for one of uncertain variables
global pts = approxConv(fg, :x2x3x4x5f1, :x5, N=N)

@test 0.1*N < sum(80 .< pts .< 100.0) < 0.5*N
@test 0.1*N < sum(pts .== 2.0) < 0.5*N
@test 0.1*N < sum(pts .== 3.0) < 0.5*N

@test sum(80 .< pts .< 100.0) + sum(pts .== 2.0) + sum(pts .== 3.0) == N

end

#


text:
using JuMP, Gurobi, Distributions, StatsFuns

const gurobi_env = Gurobi.Env()

# DIMENSIONS
const numVariables = 100
const numJCC = 100
const conRHS = numVariables^2

# CASES
# number of samples, replicates for scenario approximations
const NumSamples = ceil.(Int64,logspace(1,5,50))
const NumReplicates = 20

# FILENAMES
const objFile = "objectiveValue.txt"
const solnFile = "solution.txt"
const riskFile = "riskLevel.txt"
const solnTimeFile = "solutionTime.txt"
const checkingTimeFile = "checkingTime.txt"
const sortingTimeFile = "sortingTime.txt"
const callbackFile = "numCallbacks.txt"
const numConstraintsFile = "numConstraintsEnforced.txt"

# compute negative probability of all constraints being satisfied
# for given scenarios of demands d and a given decision vector x
function computeRiskLevel(x::Array{Float64})

	x_sqr = Array{Float64}(0)
	for j = 1:numVariables
		if(x[j] > 0)
			push!(x_sqr,(x[j])^2)
		end
	end

	n::Int64 = size(x_sqr,1)
	x_sqr = sort(x_sqr)

	if(x_sqr[n] > x_sqr[1])
		baseFactor1::Float64 = 1.0
		for i = 1:n-1
			baseFactor1 *= sqrt(x_sqr[n]/x_sqr[i])
		end

		k::Int64 = 1
		error_tol::Float64 = 1E-09/numJCC
		error::Float64 = 1.0
		while(error > error_tol)
			k *= 2
			baseFactor2 = exp(lgamma(n/2.0 + k) - lgamma(n/2.0))
			baseFactor3 = exp(k*log(x_sqr[n]/x_sqr[1] - 1.0) - k*log(2) - lfact(k))
			baseFactor4 = chisqcdf(n+2*k,conRHS/x_sqr[n])
			error = baseFactor1*baseFactor2*baseFactor3*baseFactor4
		end

		p::Float64 = 2.0/(1.0/x_sqr[1] + 1.0/x_sqr[n])
		g_base = 1.0 - p./x_sqr
		g = zeros(k-1)
		for i = 1:k-1
			g[i] = sum(g_base.^i)
		end

		c = zeros(Float64,k)
		c[1] = 1.0
		for j = 1:n
			c[1] *= sqrt(p/x_sqr[j])
		end
		for j = 1:k-1
			for r = 1:j
				c[j+1] += g[j+1-r]*c[r]
			end
			c[j+1] /= (2.0*j)
		end

		prob = 0.0
		for j = 0:k-1
			prob += c[j+1]*chisqcdf(n+2*j,conRHS/p)
		end

	else
		prob = chisqcdf(n,conRHS/x_sqr[1])
	end

	riskLevel::Float64 = 1.0 - (prob)^(numJCC)

	return riskLevel
end

# generate random samples of the demands
function generateRandomSamples(numSamples::Int64)

	xi = rand(Normal(0.0,1.0), numJCC, numVariables, numSamples)
	return xi
end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function getScenarioConstraintViolations(x::Array{Float64},xi::SubArray{Float64,2,Array{Float64,3},Tuple{Int64,Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true},pickFirstSetOfConstraints::Bool,maxNumConstraintsPerIteration::Int64)

	const numSamples = size(xi,2)
	numViolated::Int64 = 0
	constraintViolations = Float64[]
	constraintViolationIndices = Int64[]
	constraintViolationTolerance::Float64 = 1E-06

	for iter = 1:numSamples
		@views con::Float64 = norm(xi[:,iter].*x) - sqrt(conRHS)

		if(con >= constraintViolationTolerance)
			numViolated += 1
			push!(constraintViolations,con)
			push!(constraintViolationIndices,iter)
			if(pickFirstSetOfConstraints && numViolated == maxNumConstraintsPerIteration)
				break
			end
		end
	end

	return numViolated, constraintViolations, constraintViolationIndices
end

# write solution to file
function initializeFiles(details_file::String)

	writeMode::String = "w"

	open(details_file, writeMode) do f
		write(f,"Model: $modelName \n")
		write(f,"numVariables: $numVariables \n")
		write(f,"numJCC: $numJCC \n")
		write(f,"conRHS: $conRHS \n")
		write(f,"NumSamples: $NumSamples \n")
		write(f,"NumReplicates: $NumReplicates \n")
		write(f,"Method: Scenario approximations \n")
	end

end

# write solution to file
function writeSolutionToFile(dirName::String, status::Symbol, objvar_opt::Float64, x_opt::Array{Float64}, riskLevel::Float64, solutionTime::Float64, checkingTime::Float64, sortingTime::Float64, numCallbacks::Int64, numConstraintsAdded::Int64)

	if(status == Symbol("Optimal"))
		writeMode::String = "a"

		obj_file::String = dirName * objFile
		open(obj_file, writeMode) do f
			write(f,"$objvar_opt \n")
		end

		solution_file::String = dirName * solnFile
		open(solution_file, writeMode) do f
			for i = 1:numVariables
				write(f,"$(x_opt[i]) ")
			end
			write(f,"\n")
		end

		risk_file::String = dirName * riskFile
		open(risk_file, writeMode) do f
			write(f,"$riskLevel \n")
		end

		solntime_file::String = dirName * solnTimeFile
		open(solntime_file, writeMode) do f
			write(f,"$solutionTime \n")
		end

		checktime_file::String = dirName * checkingTimeFile
		open(checktime_file, writeMode) do f
			write(f,"$checkingTime \n")
		end

		sorttime_file::String = dirName * sortingTimeFile
		open(sorttime_file, writeMode) do f
			write(f,"$sortingTime \n")
		end

		callback_file::String = dirName * callbackFile
		open(callback_file, writeMode) do f
			write(f,"$numCallbacks \n")
		end

		constraints_file::String = dirName * numConstraintsFile
		open(constraints_file, writeMode) do f
			write(f,"$numConstraintsAdded \n")
		end
	end

end

# solve the scenario approximation model
function solveScenarioApproximationModel(numScenarios::Int64,dirName::String)

	tic()

	# sample the random variables
	xi = generateRandomSamples(numScenarios)

	maxNumConstraintsPerIteration::Int64 = 10
	numConstraintsAdded::Int64 = 0

	pickFirstSetOfConstraints::Bool = false

	sortingTime::Float64 = 0.0
	checkingTime::Float64 = 0.0
	numCallbacks::Int64 = 0

	objvar_opt::Float64 = Inf
	y_opt = zeros(Float64,numVariables)
	status::Symbol = Symbol("Error")

	numViolations::Int64 = numScenarios

	constraintIndices = [Int64[] for i=1:numJCC]

	while(numViolations > 0)

		numCallbacks += 1
		numViolations = 0

		numConstraintsConsidered = zeros(Int64,numJCC)
		for i = 1:numJCC
			numConstraintsConsidered[i] = size(constraintIndices[i])[1]
		end

		mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

		@variable(mod, objvar)
		@variable(mod, 0 <= y[1:numVariables] <= 10^6)

		@objective(mod, Min, objvar)
		@constraint(mod, -sum(y[j] for j = 1:numVariables) <= objvar)
		@constraint(mod, [i=1:numJCC, scen=1:numConstraintsConsidered[i]], norm(xi[i,:,constraintIndices[i][scen]].*y) <= sqrt(conRHS))

		status = solve(mod)
		objvar_opt = getvalue(objvar)
		y_opt = getvalue(y)

		for i = 1:numJCC

			numViolatedConstraints::Int64 = 0
			constraintViolations = Float64[]
			constraintViolationIndices = Int64[]

			tic()

			@views numViolatedConstraints, constraintViolations, constraintViolationIndices = getScenarioConstraintViolations(y_opt,xi[i,:,:],pickFirstSetOfConstraints,maxNumConstraintsPerIteration)

			checkingTime += toq()

			numConstraintsEnforced::Int64 = min(size(constraintViolations,1),maxNumConstraintsPerIteration)

			tic()
			relevantIndices = 1:1:numConstraintsEnforced
			if(!pickFirstSetOfConstraints)
				relevantIndices = sortperm(constraintViolations, rev=true)
			end
			sortingTime += toq()

			for iter = 1:numConstraintsEnforced
				push!(constraintIndices[i],constraintViolationIndices[relevantIndices[iter]])
			end

			numConstraintsAdded += numConstraintsEnforced

			numViolations += numViolatedConstraints
		end

	end

	riskLevel::Float64 = computeRiskLevel(y_opt)

	solutionTime::Float64 = toq()

	@printf "  Objective: %2.5f,  Risklevel: %.6f,  Time: %4.2f,  NumCallbacks: %d,  NumConstraints: %d \n" objvar_opt riskLevel solutionTime numCallbacks numConstraintsAdded

	writeSolutionToFile(dirName, status, objvar_opt, y_opt, riskLevel, solutionTime, checkingTime, sortingTime, numCallbacks, numConstraintsAdded)

	return status, objvar_opt, y_opt, riskLevel, solutionTime, checkingTime, sortingTime, numCallbacks, numConstraintsAdded
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
using Pkg

# add dependencies
# ----------------

package_dir = normpath(joinpath(@__DIR__, ".."))
atomjl_file = joinpath(package_dir, "src", "Atom.jl")
project_file = joinpath(package_dir, "Project.toml")
test_file = joinpath(package_dir, "test", "runtests.jl")
precompile_file = joinpath(package_dir, "src", "precompile.jl")

lines = readlines(atomjl_file; keep = true)

toml = Pkg.TOML.parsefile(project_file)
test_deps = get(toml, "extras", nothing)

@info "Adding temporary dependencies ..."
test_deps !== nothing && Pkg.add([PackageSpec(; name = name, uuid = uuid) for (name, uuid) in test_deps])
Pkg.add("SnoopCompile")

# generate and assert
# -------------------

using SnoopCompile

try
    @info "Generating `precompile` statements ..."
    try
        open(atomjl_file, "w") do io
            for line in lines
                if occursin("_precompile_()", line)
                    write(io, "# _precompile_()\n") # comment out
                else
                    write(io, line)
                end
            end
        end
        @debug "Commented out `_precompile_` call in $atomjl_file for `precompile` statement generation"

        inf_timing = @snoopi include(test_file)
        pc = SnoopCompile.parcel(inf_timing; blacklist=["Main"]) # NOTE: don't include functions used in test

        if (stmts = get(pc, :Atom, nothing)) !== nothing
            open(precompile_file, "w") do io
                println(io, "# This file is mostly generated by `scripts/generate_precompile.jl`\n")
                if any(str->occursin("__lookup", str), stmts)
                    println(io, SnoopCompile.lookup_kwbody_str)
                end
                println(io, "function _precompile_()")
                println(io, "    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing")
                for stmt in sort(stmts)
                    if startswith(stmt, "isdefined")
                        println(io, "    try; $(stmt); catch err; @debug err; end") # don't assert on this
                    else
                        println(io, "    try; @assert($(stmt)); catch err; @debug err; end")
                    end
                end
                println(io, "end")
            end
        end
    catch e
        printstyled(
            "`precompile` statement generation failed with the following error:\n";
            bold = true, color = :lightred
        )
        @error e
    end

    @info "Asserting generated statements ..."
    try
        open(atomjl_file, "w") do io
            for line in lines
                if occursin("# _precompile_()", line)
                    write(io, "_precompile_()\n") # comment in
                else
                    write(io, line)
                end
            end
        end
        @debug "Commented in `_precompile_` call in $atomjl_file for `precompile` statement assertion"

        run(pipeline(`julia --project=. --color=yes -e '
            ENV["JULIA_DEBUG"] = "Atom"
            using Atom # should invoke precompile
        '`))
    catch e
        printstyled(
            "`precompile` statement assertion failed with the following error:\n";
            bold = true, color = :lightred
        )
        @error e
    finally
        # let lines = readlines(precompile_file; keep = true)
        #     open(precompile_file, "w") do io
        #         for line in lines
        #             write(io, replace(line, "@assert " => ""))
        #         end
        #     end
        # end
        # @info "Removed `@assert`s in $precompile_file"
    end
catch e
    printstyled(
        "Unexpected error happened:\n";
        bold = true, color = :lightred
    )
    @error e
finally
    open(atomjl_file, "w") do io
        for line in lines
            write(io, line)
        end
    end
    @info "Restored the original state in $atomjl_file"
end

# remove dependencies
# -------------------

@info "Removing temporary dependencies ..."
test_deps === nothing || Pkg.rm([PackageSpec(; name = name, uuid = uuid) for (name, uuid) in test_deps])
Pkg.rm("SnoopCompile")


text:
module Runtests

using Test, FaultDetectionTools

@testset "Test FaultDetectionTools" begin
# test FDI tools
include("test_fdiutils.jl")
# test FDI analysis tools
include("test_fdigenspec.jl")
# test FDI synthesis functions
include("test_efdsyn.jl")
include("test_efdisyn.jl")
include("test_afdsyn.jl")
include("test_afdisyn.jl")
include("test_emmsyn.jl")
include("test_ammsyn.jl")
# test model detection functions
include("test_emdsyn.jl")
end

end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
#!/usr/bin/env julia

using DataFrames
import JLD

import Celeste.AccuracyBenchmark
import Celeste.Log
import Celeste.ParallelRun

data_rows = DataFrame[]
for jld_path in ARGS
    Log.info("Reading $jld_path...")
    sources = JLD.load(jld_path)["results"]
    Log.info("Found $(length(sources)) sources")
    for source in sources
        push!(data_rows,
              AccuracyBenchmark.variational_parameters_to_data_frame_row(
                  source.vs)
              )
    end
end
data = vcat(data_rows...)

if length(ARGS) == 1
    csv_path = string(splitext(splitdir(ARGS[1])[2])[1], ".csv")
else
    csv_path = "celeste_catalog.csv"
end
Log.info("Writing $csv_path...")
writetable(csv_path, data)


text:
  push!(LOAD_PATH, "F:\\julia_files\\lfecon\\src")
  using lfecon
  using Distributions
  using PyPlot

  include("mw_lfst_pv.jl")
  include("mw_lrstat_vec.jl")
  include("mw_marg_lf.jl")
  include("mw_rw_cov.jl")
  include("mw_s_cov_ha.jl")
  include("mw_h_cov_ha_i0_i1_i2.jl") # ok
  include("mw_h_test_i0_i1_pv.jl") # ok
  #include("mw_stest_pv.jl")
  param_lfst_m = 10
  param_lfst_t = 20
  param_lfur_m = 14
  param_h = 6
  param_s = 5
  nsim = 100
  n_wsim = 100
  n_dsim = 100
  t1 = 500; t2 = 1000
  rnd_delta = randn(t1, t2)
  rmatrix_w_string = rnd_delta
  t1 = 1100; t2 = 1000
  rnd_w = randn(t1, t2)
  rmatrix_d_string =  rnd_w
  itrend = 1
  dat = Simul_ARMA(T = 123)
  x_t = psi_compute(length(dat), 12)'dat
  # Htest using I(0) Null
  ind_int = 0
  pv_h_i0, h_i0_stat=  mw_h_test_i0_i1_pv(x_t,param_h,itrend,ind_int,n_wsim,rmatrix_w_string,nsim)
  # Htest using UR Null @
  ind_int = 1
  pv_h_i1, h_i1_stat = mw_h_test_i0_i1_pv(x_t,param_h,itrend,ind_int,n_wsim,rmatrix_w_string,nsim)
# Stest
  q = size(x_t, 1)
  psi = psi_compute(100, q+1)#, itrend)
  psi = psi[:, 1:q]
  cov_h0_i0 = psi'psi
  cov_u = mw_rw_cov(100)
  cov_h0_i1 = psi'cov_u*psi

  pv_s_i0=mw_stest_pv(x_t,cov_h0_i0,param_s,n_dsim,rmatrix_d_string,nsim)


text:
Fs = 16000
Nf = 40 # input feat dim
# feature extraction (saved to data)
get_feats(x) = identity(x)
# on the fly feature processing
function feats_post(x)
  Fs_original=20000
  Fs=16000
  x = load(x)
  if typeof(x) <: Tuple
    x = x[1]
  else
    x = x.data
  end
  x = x[:]
  x = resample(x,Fs//Fs_original)
  x = mfcc(x, float(Fs);
    dither=true, numcep=40, nbands=40, minfreq=20.0, maxfreq=7600.0)[1]
  return Float32.(x)
end
subsample = 3 # out subsampling


text:
abstract type AbstractFeasibility end

initialize_feasibility!(lshaped::AbstractLShapedSolver) = initialize_feasibility!(lshaped, lshaped.feasibility)
update_objective!(lshaped::AbstractLShapedSolver, cut::AbstractHyperPlane) = update_objective!(lshaped, cut, lshaped.feasibility)
active_model_objectives(lshaped::AbstractLShapedSolver) = active_model_objectives(lshaped, lshaped.feasibility)

# Fallback for objective update
update_objective!(lshaped::AbstractLShapedSolver, cut::AbstractHyperPlane, ::AbstractFeasibility) = nothing

struct IgnoreFeasibility <: AbstractFeasibility end

initialize_feasibility!(::AbstractLShapedSolver, ::IgnoreFeasibility) = nothing
active_model_objectives(lshaped::AbstractLShapedSolver, ::IgnoreFeasibility) = fill(true, nthetas(lshaped))
handle_feasibility(::IgnoreFeasibility) = false
initial_theta_coefficient(::IgnoreFeasibility) = 1.0
ncuts(::IgnoreFeasibility) = 0

struct HandleFeasibility{T <: AbstractFloat} <: AbstractFeasibility
    active_model_objectives::Vector{Bool}
    cuts::Vector{SparseFeasibilityCut{T}}

    function HandleFeasibility(::Type{T}) where T <: AbstractFloat
        return new{T}(Vector{Bool}(), Vector{SparseFeasibilityCut{T}}())
    end
end

function initialize_feasibility!(lshaped::AbstractLShapedSolver, feasibility::HandleFeasibility)
    append!(feasibility.active_model_objectives, fill(false, nthetas(lshaped)))
end

function active_model_objectives(::AbstractLShapedSolver, feasibility::HandleFeasibility)
    return feasibility.active_model_objectives
end

function update_objective!(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}, feasibility::HandleFeasibility)
    # Ensure that θi is included in minimization if feasibility cuts are used
    c = MPB.getobj(lshaped.mastersolver.lqmodel)
    if c[length(lshaped.x) + cut.id] == 0.0
        c[length(lshaped.x) + cut.id] = 1.0
        feasibility.active_model_objectives[cut.id] = true
        MPB.setobj!(lshaped.mastersolver.lqmodel,c)
    end
    return nothing
end
function update_objective!(lshaped::AbstractLShapedSolver, cut::AggregatedOptimalityCut, feasibility::HandleFeasibility)
    # Ensure that θis are included in minimization if feasibility cuts are used
    c = MPB.getobj(lshaped.mastersolver.lqmodel)
    ids = filter(i->c[length(lshaped.x)+i] == 0.0, cut.ids)
    if !isempty(ids)
        c[length(lshaped.x) .+ ids] .= 1.0
        feasibility.active_model_objectives[ids] .= true
        MPB.setobj!(lshaped.mastersolver.lqmodel,c)
    end
    return nothing
end

handle_feasibility(::HandleFeasibility) = true
initial_theta_coefficient(::HandleFeasibility) = 0.0
ncuts(feasibility::HandleFeasibility) = length(feasibility.cuts)


text:
import Base.StatusActive, Base.StatusOpen, Base.StatusPaused

abstract type RedisConnectionBase end

mutable struct RedisConnection <: RedisConnectionBase
    host::AbstractString
    port::Integer
    password::AbstractString
    db::Integer
    socket::TCPSocket
end

function RedisConnection(; host="127.0.0.1", port=6379, password= "", db=0)
    try
        socket = connect(host, port)
        connection = RedisConnection(host, port, password, db, socket)
        on_connect(connection)
    catch
        throw("Failed to connect to Redis server")
    end
end

function RedisConnection(conn::RedisConnection)
    try
        socket = connect(conn.host, conn.port)
        connection = RedisConnection(conn.host, conn.port, conn.password, conn.db, socket)
        on_connect(connection)
    catch
        throw("Failed to connect to Redis server")
    end
end

@inline function reconnect(conn::RedisConnectionBase)
      conn.socket = connect(conn.host, conn.port)
      on_connect(conn)
end

@inline function on_connect(conn::RedisConnectionBase)
    conn.password != "" && auth(conn, conn.password)
    conn.db != 0        && select(conn, conn.db)
    conn
end

@inline disconnect(conn::RedisConnectionBase) = close(conn.socket)
@inline is_connected(conn::RedisConnectionBase) =conn.socket.status in [StatusActive ,StatusOpen,StatusPaused]
@inline send_command(conn::RedisConnectionBase, command::AbstractString) =write(conn.socket, command)

# redis 回复消息解析
struct redisreply{T} end

function reply(conn::TCPSocket)
    tmp = readline(conn)::String
    syms, value = tmp[1] , tmp[2:end]
    reply(redisreply{Symbol(syms)}, value, conn)
end

function reply(::Type{redisreply{:*}}, value::AbstractString, conn::TCPSocket)
    num = parse(Int, value)
    num == -1 && return nothing
    num == 0 && return nothing
    res = Vector{Union{Int, String, Nothing}}(undef, num)
    for i in 1:num
        @inbounds res[i] = reply(conn)
    end
    return res
end

function reply(::Type{redisreply{:$}}, value::AbstractString, conn::TCPSocket)
    num = parse(Int, value)
    num == -1 ? nothing : read(conn, num +2) |> q -> view(q, 1:num) |> String
end

reply(::Type{redisreply{:(:)}}, value::AbstractString, conn::TCPSocket) = parse(Int, value)
reply(::Type{redisreply{:+}}, value::AbstractString, conn::TCPSocket)   = value
reply(::Type{redisreply{:-}}, value::AbstractString, conn::TCPSocket)   = throw(value)
reply(conn::TCPSocket , num::Int) = num >=1 ? reply(redisreply{:*}, string(num), conn) : nothing

# 输入redis 命令拼接函数
execute_send(conn::RedisConnectionBase, command::AbstractVector) = send_command(conn, pack_command(command))
execute_send(conn::RedisConnectionBase, command::AbstractString) = send_command(conn, command)

@inline function execute_reply(conn::RedisConnectionBase, command::Vector{String})
    execute_send(conn, command)
    reply(conn.socket)
end

function pack_command(command::Vector{String})::String
    packed_command = "*$(length(command))\r\n"
    for token in command
        packed_command = string(packed_command, "\$$(length(token))\r\n", token, "\r\n")
    end
    packed_command
end

Merge_parameters(command...) = vcat(map(Merge_parameter, command)...)::Vector{String}

Merge_parameter(token::Symbol) = string(token)
Merge_parameter(token::Number) = string(token)
Merge_parameter(token::AbstractString) = string(token)
Merge_parameter(token::Vector)  = map(Merge_parameter, token)::Vector{String}
Merge_parameter(token::Set) = json(token)
Merge_parameter(token::Dict) = json(token)
Merge_parameter(token::Tuple)  = json(token)

# 生成函数 宏

extra(d::Expr) = d.head == :(::) ? d.args[1] : d
extra(d::Symbol) = d
extra(d::AbstractString) = d

function genfunction(  kw::Vector )

    func = Expr(:call , Symbol(kw[1]) ,Expr(:(::) , :conn , :RedisConnection))
    func1 = Expr(:call , Symbol(kw[1]))

    length(kw) > 1 && begin
                        append!(func.args, kw[2:end] )
                        append!(func1.args, kw[2:end] )
                        end
    tmp = [ extra(i) for i in kw ]

    block = Expr(:block , Expr(:call, :execute_reply, :conn, Expr(:call, :Merge_parameters ,tmp... )) )
    block1 = Expr(:block ,  Expr(:call, :pack_command , Expr(:call, :Merge_parameters ,tmp... ) ))

   esc(Expr(:block , Expr(:function , func, block), Expr(:function , func1, block1)))

end

macro genfunction( kw... )
     genfunction( collect(kw) )
end

function pipe_trans(conn::RedisConnectionBase, comms::Vector{ <: AbstractString}, num::Int)
     execute_send(conn, join(comms) )
     reply(conn.socket, num )
end

function pipe_trans(conn::RedisConnectionBase, comms::AbstractString, num::Int)
     execute_send(conn, comms  )
     reply(conn.socket, num )
end

pipelines(conn::RedisConnectionBase, fun... ) = pipe_trans(conn, collect(fun), length(collect(fun)))
pipelines(conn::RedisConnectionBase, fun::AbstractString, num::Int ) = pipe_trans(conn, fun,  num )

function transactions(conn::RedisConnectionBase, fun... )
     comms = [multi(), collect(fun)..., exec()]
     pipe_trans(conn, comms, length(comms))
end

function transactions(conn::RedisConnectionBase, fun::AbstractString, num::Int )
     comms = join[multi(), fun , exec()]
     pipe_trans(conn, comms, num )
end

max_stars_repo_name: frankill/jredis

max_issues_repo_name: frankill/jredis

max_forks_repo_name: frankill/jredis


text:
#! /usr/bin/env julia

# Copyright 2017 John Hanley.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# The software is provided "AS IS", without warranty of any kind, express or
# implied, including but not limited to the warranties of merchantability,
# fitness for a particular purpose and noninfringement. In no event shall
# the authors or copyright holders be liable for any claim, damages or
# other liability, whether in an action of contract, tort or otherwise,
# arising from, out of or in connection with the software or the use or
# other dealings in the software.

# usage:
#     src/file_source /some/path [ more paths ...]
#
# Output shall be one of "DEB pkg_name", "URL url", or "UNKNOWN".

using SQLite

function file_source(db, qpath)
    sel = "select src, pkg  from file_source  where pathname = ?"
    return get(SQLite.query(db, sel; values=[qpath])[1, 2])
end

qpath, = ARGS  # the filepath we are querying
db = SQLite.DB("/tmp/file_source/pkg_contents.sqlite")
src = file_source(db, qpath)
println(src)


text:
using TrainView

using DataFrames, CSV

# Some track and camera properties
    focal_length = 40e-3
    track_gauge = 1.435

    # horizontal curvature
    β = 0.0

    # vertical curvature
    α = 0.0

    left_track  = [[x, -track_gauge/2.0 + β*x^2, α*x^2] for x = LinRange(4.0,40.0,40)]
    right_track = [[x, track_gauge/2.0 + β*x^2, α*x^2] for x = LinRange(4.0,40.0,40)]

# The reference track image
    cameraposition_reference = [0.0,0.0,2.4];
    camera_reference = VideoCamera(cameraposition_reference;
        focal_length = focal_length
        # , ψθφ = [0.0,-5.0,0.0] .* (pi/180.0)
    )

    left_ref_uvs = camera_image(camera_reference, left_track);
    right_ref_uvs = camera_image(camera_reference, right_track);

# The disturbed track image
    cameraposition = cameraposition_reference + [0.0,-0.06,0.05];
    camera = VideoCamera(cameraposition;
        focal_length = focal_length
        , ψθφ = [1.05,-1.2,0.001] .* (pi/180.0)
    )

    uvs_left = camera_image(camera, left_track);
    uvs_right = camera_image(camera, right_track);

# Save the results
    data_ref = DataFrame(
        left_ref_u = [uv[1] for uv in left_ref_uvs],
        left_ref_v = [uv[2] for uv in left_ref_uvs],
        right_ref_u = [uv[1] for uv in right_ref_uvs],
        right_ref_v = [uv[2] for uv in right_ref_uvs]
    )
    data = DataFrame(
        left_u = [uv[1] for uv in uvs_left],
        left_v = [uv[2] for uv in uvs_left],
        right_u = [uv[1] for uv in uvs_right],
        right_v = [uv[2] for uv in uvs_right]
    )

# CSV.write("reference_image.csv", data_ref)
# CSV.write("disturbed_image.csv", data)

# Plot the results
using Plots

w = camera.opticalproperties.sensor_width;
h = camera.opticalproperties.sensor_height;

plot([-w/2.0,w/2.0,w/2.0,-w/2.0,-w/2.0],[-h/2.0,-h/2.0,h/2.0,h/2.0,-h/2.0], color=:yellow, lab = "")

w = camera_reference.opticalproperties.sensor_width;
h = camera_reference.opticalproperties.sensor_height;

plot!([-w/2.0,w/2.0,w/2.0,-w/2.0,-w/2.0],[-h/2.0,-h/2.0,h/2.0,h/2.0,-h/2.0], color=:red, lab = "")

scatter!([uv[1] for uv in left_ref_uvs],[uv[2] for uv in left_ref_uvs])
scatter!([uv[1] for uv in right_ref_uvs],[uv[2] for uv in right_ref_uvs])
scatter!([uv[1] for uv in left_uvs],[uv[2] for uv in left_uvs])
scatter!([uv[1] for uv in right_uvs],[uv[2] for uv in right_uvs])


text:
#!/usr/bin/env julia

function packagelt(a::String, b::String)
    a == "julia" && b != "julia" && return true
    a != "julia" && b == "julia" && return false
    return lowercase(a) < lowercase(b)
end

≲(v::VersionNumber, t::NTuple{0,Integer}) = true
≲(v::VersionNumber, t::NTuple{1,Integer}) = v.major ≤ t[1]
≲(v::VersionNumber, t::NTuple{2,Integer}) = v.major < t[1] ||
                                            v.major ≤ t[1] && v.minor ≤ t[2]
≲(v::VersionNumber, t::NTuple{3,Integer}) = v.major < t[1] ||
                                            v.major ≤ t[1] && v.minor < t[2] ||
                                            v.major ≤ t[1] && v.minor ≤ t[2] && v.patch ≤ t[3]

≲(t::NTuple{0,Integer}, v::VersionNumber) = true
≲(t::NTuple{1,Integer}, v::VersionNumber) = t[1] ≤ v.major
≲(t::NTuple{2,Integer}, v::VersionNumber) = t[1] < v.major ||
                                            t[1] ≤ v.major && t[2] ≤ v.minor
≲(t::NTuple{3,Integer}, v::VersionNumber) = t[1] < v.major ||
                                            t[1] ≤ v.major && t[2] < v.minor ||
                                            t[1] ≤ v.major && t[2] ≤ v.minor && t[3] ≤ v.patch

function compress_versions(inc::Vector{VersionNumber}, from::Vector{VersionNumber})
    issorted(inc) || (inc = sort(inc))
    exc = sort!(setdiff(from, inc))
    pairs = []
    if isempty(exc)
        lo, hi = first(inc), last(inc)
        push!(pairs, (lo.major, lo.minor) => (hi.major, hi.minor))
    else
        for v in inc
            t = (v.major, v.minor)
            if any(t ≲ w ≲ t for w in exc)
                t = (v.major, v.minor, v.patch)
            end
            if isempty(pairs) || any(pairs[end][1] ≲ w ≲ t for w in exc)
                push!(pairs, t => t) # need a new interval
            else
                pairs[end] = pairs[end][1] => t # can be merged with last
            end
        end
    end
    @assert all(any(p[1] ≲ v ≲ p[2] for p ∈ pairs) for v ∈ inc)
    @assert all(!any(p[1] ≲ v ≲ p[2] for p ∈ pairs) for v ∈ exc)
    return pairs
end
compress_versions(f::Function, from::Vector{VersionNumber}) =
    compress_versions(filter(f, from), from)
compress_versions(vi::VersionInterval, from::Vector{VersionNumber}) =
    compress_versions(v->v in vi, from)
compress_versions(inc, from) = compress_versions(inc, collect(from))

versions_string(p::Pair) = versions_string(p...)
versions_string(a::Tuple{}, b::Tuple{}) = "*"
versions_string(a::NTuple{m,Integer}, b::NTuple{n,Integer}) where {m,n} =
    a == b ? join(a, '.') : "$(join(a, '.'))-$(join(b, '.'))"

versions_repr(x) = repr(versions_string(x))
versions_repr(v::Vector) = length(v) == 1 ? repr(versions_string(v[1])) :
    "[" * join(map(repr∘versions_string, v), ", ") * "]"

## compress per-version data ##

function compress_versions_data(
    d::Dict{VersionNumber,Dict{String,String}},
    versions::Vector{VersionNumber},
)
    kvs = Dict{Pair{String,String},Vector{VersionNumber}}()
    for (ver, x) in d, (k, v) in x
        push!(get!(kvs, k => v, VersionNumber[]), ver)
    end
    tx = Tuple{Any,String,String}[]
    for (kv, vs) in kvs
        for v in compress_versions(sort!(vs), versions)
            push!(tx, (v, kv...))
        end
    end
    sort!(tx, by=t->(t[1][1], t[1][2], t[2], t[3]))
    return Tuple{String,String,String}[
        (versions_string(t[1]), t[2], t[3]) for t in tx
    ]
end

## dict utility functions ##

function invert_map(fwd::Dict{K,V}) where {K,V}
    rev = Dict{V,Vector{K}}()
    for (k, v) in fwd
        push!(get!(rev, v, K[]), k)
    end
    return rev
end

function invert_map(fwd::Dict{Vector{K},V}) where {K,V}
    rev = Dict{V,Vector{K}}()
    for (k, v) in fwd
        append!(get!(rev, v, K[]), k)
    end
    return rev
end

flatten_keys(d::Dict{Vector{K},V}) where {K,V} =
    isempty(d) ? Dict{K,V}() : Dict{K,V}(k => v for (ks, v) in d for k in ks)

## TOML generation functions ##

function write_toml(f::Function, names::String...)
    path = joinpath(names...) * ".toml"
    mkpath(dirname(path))
    open(path, "w") do io
        f(io)
    end
end

toml_key(str::String) = contains(str, r"[^\w-]") ? repr(str) : str
toml_key(strs::String...) = join(map(toml_key, [strs...]), '.')

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
"""
    echotime(f)

Get echo time `TE` of the time profile `f`, which is the end of the last characteristic
interval.
"""
echotime(f::TimeProfile) = intervals(f)[end]

"""
    echotime(gradient)

Get echo time `TE` of `gradient`.
"""
echotime(gradient::GeneralGradient) = gradient.TE
echotime(gradient::ScalarGradient) = echotime(gradient.profile)


text:
include("/usr/share/julia/base/latex_symbols.jl")

# find(latex_symbols, "⟺")  # returns  "\\Longleftrightarrow"
open(joinpath(dirname(@__FILE__), "..", "latex.ts"), "w") do f
  println(f, "export const latexSymbols = {")
  # sort by name length and name content
  symbols = sort(collect(latex_symbols))
  symbols = sort(symbols, by=x->length(x[1]))
  for (name, sym) in symbols
    println(f, "    '" * name[2:length(name)] * "': '$sym',")
    # println(f, "    '\\$name': '$sym',")
  end
  println(f, "};")
end


text:
"""
    bus_dist(b1::Bus, b2::Bus)

OBSOLETE. Compute geometrical distance between two buses. Obsolete now that
buses are located by latitude and longitude.
"""
function bus_dist(b1::Bus, b2::Bus)
    return sqrt((b1.coords[1] - b2.coords[1])^2 + (b1.coords[2] - b2.coords[2])^2)
end

"""
    haversine{T<:Real}(c1::Tuple{T,T}, c2::Tuple{T,T})

Compute distance between two pairs of latitude and longitude using the
haversine formula.
"""
function haversine{T<:Real}(c1::Tuple{T,T}, c2::Tuple{T,T})
    const AVG_RADIUS = 6369.783 # Earth radius at 38.8 latitude (in km).
    d = 2 * AVG_RADIUS * asin(sqrt(sind((c2[1] - c1[1]) / 2)^2
        + cosd(c2[1]) * cosd(c1[1]) * sind((c2[2] - c1[2]) / 2)^2))
    return d
end

function haversine(c1::LatLon, c2::LatLon)
    return haversine((c1.lat, c1.lon), (c2.lat, c2.lon))
end

"""
    haversine(b1::Bus, b2::Bus)

Compute distance between two buses from their latitude and longitude using the
haversine formula.
"""
function haversine(b1::Bus, b2::Bus)
    return haversine(b1.coords, b2.coords)
end

function mean{T <: Real}(v::Vector{LatLon{T}})
    xmean = mean([c.lat for c in v])
    ymean = mean([c.lon for c in v])
    return LatLon(xmean, ymean)
end

"""
    subs_dist(inds1, inds2, buses, b_dist)

Compute distance between two substations using equation 1 from the REFERENCE.

REFERENCE: Birchfield, Adam B., et al. "Grid Structural Characteristics as
Validation Criteria for Synthetic Networks."
IEEE Transactions on Power Systems (2016).
"""
function subs_dist{T <: Real}(
    inds1::Vector{Int},
    inds2::Vector{Int},
    buses::Vector{LoadBus},
    b_dist::Matrix{T},
)
    dist = 0
    totw = 0
    for i1 in inds1, i2 in inds2
        totw += (buses[i1].population + buses[i2].population)
    end
    for i1 in inds1, i2 in inds2
        dist += ((buses[i1].population + buses[i2].population)
                * b_dist[i2,i1] / totw)
    end
    return dist
end

function load_sub_dist(subs, buses, b_dist)
    b_inds = sizehint!([], length(subs)) # indexes of buses in each substations
    for sub in subs
        temp = sizehint!(Int[],length(sub.grouping))
        for load in sub.grouping
          push!(temp,load.id)
        end
        push!(b_inds, temp)
    end
    n = length(subs)
    dist_mat = fill(Inf, (n, n))
    for s1 in 1:(n - 1), s2 in (s1 + 1):n
        dist_mat[s1, s2] = subs_dist(b_inds[s1], b_inds[s2], buses, b_dist)
        dist_mat[s2, s1] = dist_mat[s1, s2]
    end
    return dist_mat
end

function distance{T<:Bus}(buses::Vector{T})
    n = length(buses)
    dist_mat = fill(Inf, (n, n))
    for b1 in 1:(n - 1), b2 in (b1 + 1):n
        dist_mat[b1, b2] = haversine(buses[b1], buses[b2])
        dist_mat[b2, b1] = dist_mat[b1, b2]
    end
    return dist_mat
end


text:
function drop_clock(cur_player::AbstractPlayer)

  ( cur_player.game.in_focus ) || return

  isnull(cur_player.clock.push) &&
    push_clock(cur_player)

  cur_uuid = UUIDs.uuid4()

  cur_player.clock.drop = cur_uuid

  cur_player.clock.lock =
    Nullable{UUIDs.UUID}()

  cur_level = cur_player.level

  cur_sleep_time = 0.85

  cur_sleep_time ^= cur_level

  cur_sleep_time *= 725

  cur_sleep_time += cur_level

  cur_sleep_time /= 1000

  @async while check_clock(cur_player, :drop, cur_uuid)
    did_step = step(cur_player)

    if did_step
      cur_player.clock.lock =
        Nullable{UUIDs.UUID}()
    else
      isnull(cur_player.clock.lock) &&
        lock_clock(cur_player)

      isnull(cur_player.clock.drag) &&
        drag_clock(cur_player)
    end

    sleep(cur_sleep_time)
  end

end


text:
np = nr * (ns + 3) + 1
p = randn(Float64, np) .* 0.05;
p[1:nr] .+= 0.8;  # w_b
p[nr * (ns + 1) + 1:nr * (ns + 2)] .+= 0.8;  # w_out
p[nr * (ns + 2) + 1:end - 1] .+= 0.1;  # w_b | w_Ea
p[end] = 0.1;  # slope

function p2vec(p)
    slope = p[end] .* 1.e1
    w_b = p[1:nr] .* (slope * 10.0)
    w_b = clamp.(w_b, 0, 50)

    w_out = reshape(p[nr + 1:nr * (ns + 1)], ns, nr)
    @. w_out[1, :] = clamp(w_out[1, :], -1.1, 0.0)
    @. w_out[end, :] = clamp(abs(w_out[end, :]), 0.0, 1.1)

    if p_cutoff > 0.0
        w_out[findall(abs.(w_out) .< p_cutoff)] .= 0.0
    end

    w_out[ns - 1:ns - 1, :] .=
        -sum(w_out[1:ns - 2, :], dims=1) .- sum(w_out[ns:ns, :], dims=1)

    w_in_Ea = abs.(p[nr * (ns + 1) + 1:nr * (ns + 2)] .* (slope * 100.0))
    w_in_Ea = clamp.(w_in_Ea, 0.0, 300.0)

    w_in_b = abs.(p[nr * (ns + 2) + 1:nr * (ns + 3)])

    # w_in_ocen = abs.(p[nr*(ns+3)+1:nr*(ns+4)])
    # w_in_ocen = clamp.(w_in_ocen, 0.0, 1.5)

    # if p_cutoff > 0.0
    #     w_in_ocen[findall(abs.(w_in_ocen) .< p_cutoff)] .= 0.0
    # end

    w_in = vcat(clamp.(-w_out, 0.0, 1.1), w_in_Ea', w_in_b')
    return w_in, w_b, w_out
end

function display_p(p)
    w_in, w_b, w_out = p2vec(p)
    println("\n species (column) reaction (row)")
    println("w_in | Ea | b | n_ocen | lnA | w_out")
    show(stdout, "text/plain", round.(hcat(w_in', w_b, w_out'), digits=2))
    # println("\n w_out")
    # show(stdout, "text/plain", round.(w_out', digits=3))
    println("\n")
end
display_p(p);

function getsampletemp(t, T0, beta)
    T = T0 + beta / 60 * t  # K/min to K/s
    return T
end

const R = -1.0 / 8.314e-3  # universal gas constant, kJ/mol*K
@inbounds function crnn!(du, u, p, t)
    logX = @. log(clamp(u, lb, 2.0))
    T = getsampletemp(t, T0, beta)
    w_in_x = w_in' * vcat(logX, R / T, log(T))
    du .= w_out * (@. exp(w_in_x + w_b))
end

tspan = [0.0, 1.0];
u0 = zeros(ns);
u0[1] = 1.0;
prob = ODEProblem(crnn!, u0, tspan, p, abstol=lb)

condition(u, t, integrator) = u[1] < lb * 5.0
affect!(integrator) = terminate!(integrator)
_cb = DiscreteCallback(condition, affect!)

alg = TRBDF2();
sense = ForwardSensitivity(autojacvec=true)
# sense = BacksolveAdjoint()
# sense = ForwardDiffSensitivity()
function pred_n_ode(p, i_exp, exp_data)
    global T0, beta, ocen
    global w_in, w_b, w_out
    T0, beta, ocen = l_exp_info[i_exp, :]
    w_in, w_b, w_out = p2vec(p)

    ts = @view(exp_data[:, 1])
    tspan = [ts[1], ts[end]]
    sol = solve(
        prob,
        alg,
        tspan=tspan,
        p=p,
        saveat=ts,
        sensalg=sense,
        maxiters=maxiters,
        # callback=_cb,
    )

    if sol.retcode == :Success
        nothing
    else
        @sprintf("solver failed beta: %.0f", beta)
    end
    if length(sol.t) > length(ts)
        return  sol[:, 1:length(ts)]
    else
        return sol
    end
end

function loss_neuralode(p, i_exp)
    exp_data = l_exp_data[i_exp]
    pred = Array(pred_n_ode(p, i_exp, exp_data))
    masslist = sum(clamp.(@view(pred[1:end - 1, :]), 0.0, Inf), dims=1)'
    gaslist = clamp.(@views(pred[end, :]), 0.0, Inf)

    s = sample(1:length(masslist), batchsize)

    loss = mae(masslist[s], @view(exp_data[s, 3])) + mae(gaslist[s], 1 .- @view(exp_data[s, 3]))
    # - sum(clamp.(pred, -Inf, -lb)) .* 1.e-3
    return loss
end
@time loss = loss_neuralode(p, 1)
# using BenchmarkTools
# @benchmark loss = loss_neuralode(p, 1)
# @benchmark grad = ForwardDiff.gradient(x -> loss_neuralode(x, 1), p)

max_stars_repo_name: DENG-MIT/CRNN-Pyrolysis

max_issues_repo_name: DENG-MIT/CRNN-Pyrolysis

max_forks_repo_name: DENG-MIT/CRNN-Pyrolysis


text:
function RCCSDpT()
    prec = select_precision(Fermi.CurrentOptions["precision"])
    RCCSDpT{prec}()
end

function RCCSDpT{T}() where T <: AbstractFloat
    refwfn = Fermi.HartreeFock.RHF()

    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    @output "Transforming Integrals..."

    ccsd_wfn = RCCSD{T}()
    RCCSDpT{T}(ccsd_wfn, refwfn.ints)
end

function RCCSDpT{T}(ccsd::RCCSD, ints::IntegralHelper) where T <: AbstractFloat

    @output "\n   • Perturbative Triples Started\n\n"

    T1 = ccsd.T1.data
    T2 = ccsd.T2.data

    Vvvvo = permutedims(ints["OVVV"], (4,2,3,1))
    Vvooo = permutedims(ints["OOOV"], (4,2,1,3))
    Vvovo = permutedims(ints["OOVV"], (3,1,4,2))

    o,v = size(T1)
    Et::T = 0.0

    fo = diag(ints["FOO"])
    fv = diag(ints["FVV"])

    # Pre-allocate Intermediate arrays
    W  = Array{T}(undef, v,v,v)
    V  = Array{T}(undef, v,v,v)

    # Pre-allocate Views
    Vvvvo_4k =    Array{T}(undef,v,v,v)
    Vvooo_2k_3j = Array{T}(undef,v,o)
    Vvooo_2k_3i = Array{T}(undef,v,o)
    Vvooo_2i_3k = Array{T}(undef,v,o)
    Vvooo_2j_3k = Array{T}(undef,v,o)
    T2_1k_2j    = Array{T}(undef,v,v)
    T2_1k_2i    = Array{T}(undef,v,v)
    T2_1i_2k    = Array{T}(undef,v,v)
    T2_1j_2k    = Array{T}(undef,v,v)
    T2_1k       = Array{T}(undef,o,v,v)
    T1_1k       = Array{T}(undef,v)
    Vvovo_2j_4k = Array{T}(undef,v,v)
    Vvovo_2i_4k = Array{T}(undef,v,v)
    Vvvvo_4i    = Array{T}(undef,v,v,v)
    T2_1i       = Array{T}(undef,o,v,v)
    T1_1i       = Array{T}(undef,v)
    Vvvvo_4j    = Array{T}(undef,v,v,v)
    Vvooo_2i_3j = Array{T}(undef,v,o)
    Vvooo_2j_3i = Array{T}(undef,v,o)
    T2_1i_2j    = Array{T}(undef,v,v)
    T2_1j_2i    = Array{T}(undef,v,v)
    T2_1j       = Array{T}(undef,o,v,v)
    Vvovo_2i_4j = Array{T}(undef,v,v)
    T1_1j       = Array{T}(undef,v)
    @output "Computing energy contribution from occupied orbitals:\n"
    for i in 1:o
        @output "→ Orbital {} of {}\n" i o
        Vvvvo_4i .= view(Vvvvo, :,:,:,i)
        T2_1i    .= view(T2, i,:,:,:)
        T1_1i    .= view(T1, i, :)
        for j in 1:i
            Vvvvo_4j    .= view(Vvvvo, :,:,:,j)
            Vvooo_2i_3j .= view(Vvooo,:,i,j,:)
            Vvooo_2j_3i .= view(Vvooo,:,j,i,:)
            T2_1i_2j    .= view(T2, i,j,:,:)
            T2_1j_2i    .= view(T2, j,i,:,:)
            T2_1j       .= view(T2, j,:,:,:)
            Vvovo_2i_4j .= view(Vvovo, :,i,:,j)
            T1_1j       .= view(T1, j, :)
            δij = Int(i == j)
            for k in 1:j
                Vvvvo_4k    .= view(Vvvvo, :,:,:,k)
                Vvooo_2k_3j .= view(Vvooo,:,k,j,:)
                Vvooo_2k_3i .= view(Vvooo,:,k,i,:)
                Vvooo_2i_3k .= view(Vvooo,:,i,k,:)
                Vvooo_2j_3k .= view(Vvooo,:,j,k,:)

                T2_1k_2j .= view(T2, k,j,:,:)
                T2_1k_2i .= view(T2, k,i,:,:)
                T2_1i_2k .= view(T2, i,k,:,:)
                T2_1j_2k .= view(T2, j,k,:,:)

                T2_1k .= view(T2, k,:,:,:)
                T1_1k .= view(T1, k, :)

                Vvovo_2j_4k .= view(Vvovo, :,j,:,k)
                Vvovo_2i_4k .= view(Vvovo, :,i,:,k)

                # Build intermediates W and V for the set i j k
                @tensoropt begin
                    W[a,b,c]  = (Vvvvo_4i[b,d,a]*T2_1k_2j[c,d] - Vvooo_2k_3j[c,l]*T2_1i[l,a,b]  # ijk abc
                              +  Vvvvo_4i[c,d,a]*T2_1j_2k[b,d] - Vvooo_2j_3k[b,l]*T2_1i[l,a,c]  # ikj acb
                              +  Vvvvo_4k[a,d,c]*T2_1j_2i[b,d] - Vvooo_2j_3i[b,l]*T2_1k[l,c,a]  # kij cab
                              +  Vvvvo_4k[b,d,c]*T2_1i_2j[a,d] - Vvooo_2i_3j[a,l]*T2_1k[l,c,b]  # kji cba
                              +  Vvvvo_4j[c,d,b]*T2_1i_2k[a,d] - Vvooo_2i_3k[a,l]*T2_1j[l,b,c]  # jki bca
                              +  Vvvvo_4j[a,d,b]*T2_1k_2i[c,d] - Vvooo_2k_3i[c,l]*T2_1j[l,b,a]) # jik bac

                    V[a,b,c]  = W[a,b,c] + Vvovo_2j_4k[b,c]*T1_1i[a] + Vvovo_2i_4k[a,c]*T1_1j[b] + Vvovo_2i_4j[a,b]*T1_1k[c]
                end

                # Compute Energy contribution
                δjk = Int(j == k)
                for a in 1:v
                    for b in 1:a
                        δab = Int(a == b)
                        for c in 1:b
                            Dd = fo[i] + fo[j] + fo[k] - fv[a] - fv[b] - fv[c]
                            δbc = Int(b == c)
                            @inbounds X = (W[a,b,c]*V[a,b,c] + W[a,c,b]*V[a,c,b] + W[b,a,c]*V[b,a,c] + W[b,c,a]*V[b,c,a] + W[c,a,b]*V[c,a,b] + W[c,b,a]*V[c,b,a])
                            @inbounds Y = (V[a,b,c] + V[b,c,a] + V[c,a,b])
                            @inbounds Z = (V[a,c,b] + V[b,a,c] + V[c,b,a])
                            E = (Y - 2*Z)*(W[a,b,c] + W[b,c,a] + W[c,a,b]) + (Z - 2*Y)*(W[a,c,b]+W[b,a,c]+W[c,b,a]) + 3*X
                            Et += E*(2-δij-δjk)/(Dd*(1+δab+δbc))
                        end
                    end
                end
            end
        end
    end
    @output "Final (T) contribution: {:15.10f}\n" Et
    @output "CCSD(T) energy:         {:15.10f}\n" Et+ccsd.CorrelationEnergy
    return RCCSDpT{T}(ccsd, Et)
end


text:
# How to run:
#    julia --project=. scripts/parallel.jl
#
# To use a different batch size, set environment variable:
#    export BATCH_SIZE=1000
#    julia --project=. scripts/parallel.jl

using Distributed
using TypeDBClient

@info "Delete database biograkn, create it again and define the schema"
rootdir = joinpath(@__DIR__, "..")
include("$rootdir/src/define_schema.jl")

TypeDBClient.dbconnect("127.0.0.1") do client
    delete_db_create_db(client, "biograkn")
end

@everywhere ENV["BATCH_SIZE"] = 1000

@info "Preparing worker processes"
procs = addprocs(21)
@everywhere rootdir = joinpath(@__DIR__, "..")
@everywhere import Pkg
@everywhere Pkg.activate(rootdir)
@everywhere include("$rootdir/src/load_dataset.jl")

@info "Loading function in master process"
include("$rootdir/src/load_dataset.jl")

@info "Loading data to master process"
queries = make_queries(load_data())

batches = collect(Iterators.partition(queries, 1000))
@info "Data is divided into: $(length(batches))"

failure_arr = String[]

const remote_channel = RemoteChannel(()->Channel{Tuple}(10))
const failure_channel = RemoteChannel(()->Channel(10))
const result_channel = RemoteChannel(()->Channel(10))

@info "Start failure collecting"
@async begin
    while isopen(failure_channel)
        yield()
        if isready(failure_channel)
            push!(failure_arr, take!(failure_channel))
        end
    end
end

@info "Warming up"
@sync for n in procs
    @info "Warming up worker $n"
    @async remotecall_wait(warmup, n)
end

@info "Start worker processes"
for n in procs
    @spawnat n parallel_channel(remote_channel, failure_channel, result_channel)
end

@info "Put work on the workers"
for i in 1:length(batches)
    @async put!(remote_channel, (i,batches[i]))
end

for _ in 1:length(batches)
    take!(result_channel)
end

length(failure_arr)
# put!(failure_channel, "hier")

# filter!(x->x != "hier", failure_arr)

# client = CoreClient("127.0.0.1", 1729)
# sess = CoreSession(client, "biograkn", Proto.Session_Type.DATA, request_timout = Inf, error_time = 20)
# trans = transaction(sess, Proto.Transaction_Type.WRITE)

# for query in failure_arr
#    @info insert(trans, query)
# end

# commit(trans)
# @everywhere batch_size = get_env_batch_size()
# @info "Batch size = $(get_env_batch_size())"

# @info "Warming up"
# @time try
#     parallel_run()
# catch ex
#     @info "Got exception... be brave, just continue!"
# end

# @info "Test runs"
# run_results = parallel_run()
# run_reuslts = parallel_run()

# # @time parallel_run()

# rmprocs(procs)
# length(filter(x->isempty(x),run_results))


text:
# # Variational problems
#
# In this example, we will numerically simulate an entropy-regularised Wasserstein gradient flow
# approximating the Fokker-Planck and porous medium equations.
#
# The connection between Wasserstein gradient flows and (non)-linear PDEs is due to Jordan, Kinderlehrer and Otto [^JKO98], and
# an easy-to-read overview of the topic is provided in Section 9.3 [^PC19]
#
# [^JKO98]: Jordan, Richard, David Kinderlehrer, and Felix Otto. "The variational formulation of the Fokker--Planck equation." SIAM journal on mathematical analysis 29.1 (1998): 1-17.
# [^PC19]: Peyré, Gabriel, and Marco Cuturi. "Computational optimal transport: With applications to data science." Foundations and Trends® in Machine Learning 11.5-6 (2019): 355-607.
#
# ## Fokker-Planck equation as a $W_2$ gradient flow
# For a potential function $\Psi$ and noise level $\sigma^2$, the Fokker-Planck equation (FPE) is
# ```math
# \partial_t \rho_t = \nabla \cdot (\rho_t \nabla \Psi) + \frac{\sigma^2}{2} \Delta \rho_t,
# ```
# and we take no-flux (Neumann) boundary conditions.
#
# This describes the evolution of a massless particle undergoing both diffusion (with diffusivity $\sigma^2$) and drift (along potential $\Psi$) according to the stochastic differential equation
# ```math
# dX_t = -\nabla \Psi(X_t) dt + \sigma dB_t.
# ```
# The result of Jordan, Kinderlehrer and Otto (commonly referred to as the JKO theorem) states that
# $\rho_t$ evolves following the 2-Wasserstein gradient flow of the Gibbs free energy functional
# ```math
#   F(\rho) = \int \Psi d\rho + \int \log(\rho) d\rho.
# ```
#
# ## Implicit schemes for gradient flows
# In an Euclidean space, the gradient flow of a functional $F$ is simply the solution of an ordinary differential equation
# ```math
#  \dfrac{dx(t)}{dt} = -\nabla F(x(t)).
# ```
# Of course, there is a requirement that $F$ is smooth. A more general formulation of a gradient flow that allows
# $F$ to be non-smooth is the implicit scheme
# ```math
#   x_{t+\tau} = \operatorname{argmin}_x \frac{1}{2} \| x - x_t \|_2^2 + \tau F(x).
# ```
# As the timestep $\tau$ shrinks, $x_t$ becomes a better and better approximation to the gradient flow of $F$.
#
# ## Wasserstein gradient flow
# In the context of the JKO theorem, we seek $\rho_t$ that is the gradient flow of $F$ with
# respect to the 2-Wasserstein distance. This can be achieved by choosing the $W_2$ metric in the implicit step:
# ```math
#   \rho_{t + \tau} = \operatorname{argmin}_{\rho} d_{W_2}^2(\rho_{t}, \rho) + \tau F(\rho).
# ```
# Finally, a numerical scheme for computing this gradient flow can be developed by using the entropic regularisation
# of optimal transport on a discretised domain
# ```math
#   \rho_{t + \tau} = \operatorname{argmin}_{\rho} \operatorname{OT}_\varepsilon(\rho_{t}, \rho) + \tau F(\rho),
# ```
# where
# ```math
#   \operatorname{OT}_\varepsilon(\alpha, \beta) = \min_{\gamma \in \Pi(\alpha, \beta)} \sum_{i,j} \frac{1}{2} \| x_i - x_j \|_2^2 \gamma_{ij} + \varepsilon \sum_{i, j} \gamma_{ij} \log(\gamma_{ij}).
# ```
# Each step of this problem is a minimisation problem with respect to $\rho$.
# Since we use entropic optimal transport which is differentiable, this can be solved using gradient-based methods.

# ## Problem setup
#
using OptimalTransport
using Distances
using LogExpFunctions
using Optim
using Plots
using StatsBase
using ReverseDiff

using LinearAlgebra
using Logging

# Here, we set up the computational domain that we work on - we discretize the interval $[-1, 1]$.
# The natural boundary conditions to use will be Neumann (zero flux), see e.g. [^Santam2017]
#
# [^Santam2017]: Santambrogio, Filippo. "{Euclidean, metric, and Wasserstein} gradient flows: an overview." Bulletin of Mathematical Sciences 7.1 (2017): 87-154.

support = range(-1, 1; length=64)
C = pairwise(SqEuclidean(), support');

# Now we set up various functionals that we will use.
#
# We define the generalised entropy (Equation (4.4) of [^Peyre2015]) as follows. For $m = 1$ this is just the "regular" entropy, and $m = 2$ this is squared $L_2$.
#
# [^Peyre2015]: Peyré, Gabriel. "Entropic approximation of Wasserstein gradient flows." SIAM Journal on Imaging Sciences 8.4 (2015): 2323-2351.
function E(ρ; m=1)
    if m == 1
        return sum(xlogx.(ρ)) - sum(ρ)
    elseif m > 1
        return dot(ρ, @. (ρ^(m - 1) - m) / (m - 1))
    end
end;

# Now define $\psi(x)$ to be a potential energy function that has two potential wells at $x = ± 0.5$.
ψ(x) = 10 * (x - 0.5)^2 * (x + 0.5)^2;
plot(support, ψ.(support); color="black", label="Scalar potential")

# Having defined $\psi$, this induces a potential energy functional $\Psi$ on probability distributions $\rho$:
# ```math
#    \Psi(\rho) = \int \psi(x) \rho(x) dx = \langle \psi, \rho \rangle .
# ```
Ψ = ψ.(support);

# Define the time step $\tau$ and entropic regularisation level $\varepsilon$, and form the associated Gibbs kernel $K = e^{-C/\varepsilon}$.
τ = 0.05
ε = 0.01
K = @. exp(-C / ε);

# We define the (non-smooth) initial condition $\rho_0$ in terms of step functions.
H(x) = x > 0
ρ0 = @. H(support + 0.25) - H(support - 0.25)
ρ0 = ρ0 / sum(ρ0)
plot(support, ρ0; label="Initial condition ρ0", color="blue")

# `G_fpe` is the objective function for the implicit step scheme
# ```math
# G_\mathrm{fpe}(\rho) = \operatorname{OT}_\varepsilon(\rho_{t}, \rho) + \tau F(\rho),
# ```
# and we seek to minimise in $\rho$.
function G_fpe(ρ, ρ0, τ, ε, C)
    return sinkhorn2(ρ, ρ0, C, ε; regularization=true, maxiter=250) + τ * (dot(Ψ, ρ) + E(ρ))
end;

# `step` solves the implicit step problem to produce $\rho_{t + \tau}$ from $\rho_t$.
function step(ρ0, τ, ε, C, G)
    ## only print error messages
    obj = u -> G(softmax(u), ρ0, τ, ε, C)
    opt = with_logger(SimpleLogger(stderr, Logging.Error)) do
        optimize(
            obj,
            ones(size(ρ0)),
            LBFGS(),
            Optim.Options(; iterations=50, g_tol=1e-6);
            autodiff=:forward,
        )
    end
    return softmax(Optim.minimizer(opt))
end
# Now we simulate `N = 10` iterates of the gradient flow and plot the result.

N = 10
ρ = similar(ρ0, size(ρ0, 1), N)
ρ[:, 1] = ρ0
for i in 2:N
    @info i
    ρ[:, i] = step(ρ[:, i - 1], τ, ε, C, G_fpe)
end
colors = range(colorant"red"; stop=colorant"blue", length=N)
plot(
    support,
    ρ;
    title=raw"$F(\rho) = \langle \psi, \rho \rangle + \langle \rho, \log(\rho) \rangle$",
    palette=colors,
    legend=nothing,
)

# ## Porous medium equation
#
# The porous medium equation (PME) is the nonlinear PDE
# ```math
# \partial_t \rho = \nabla \cdot (\rho \nabla \Psi) + \Delta \rho^m,
# ```
# again with Neumann boundary conditions. The value of $m$ in the PME corresponds to picking $m$ in the generalised entropy functional.
# Now, we will solve the PME with $m = 2$ as a Wasserstein gradient flow.
#
function G_pme(ρ, ρ0, τ, ε, C)
    return sinkhorn2(ρ, ρ0, C, ε; regularization=true, maxiter=250) +
           τ * (dot(Ψ, ρ) + E(ρ; m=2))
end;

# set up as previously
N = 10
ρ = similar(ρ0, size(ρ0, 1), N)
ρ[:, 1] = ρ0
for i in 2:N
    ρ[:, i] = step(ρ[:, i - 1], τ, ε, C, G_pme)
end
plot(
    support,
    ρ;
    title=raw"$F(\rho) = \langle \psi, \rho \rangle + \langle \rho, \rho - 1\rangle$",
    palette=colors,
    legend=nothing,
)

# ## Exploiting duality
#
# The previous examples solved the minimisation problem for the implicit gradient flow step directly, involving automatic differentiation through the Sinkhorn iterations used to compute $\operatorname{OT}_\varepsilon(\rho_t, \rho)$ each time a gradient needs to be evaluated.
# While this is straightforward to implement, it is computationally costly.
# An alternative approach for convex variational problems is to proceed via the [dual problem](https://en.wikipedia.org/wiki/Duality_(optimization)).
# The benefit of proceeding via the dual problem is that the part of the dual minimisation problem corresponding to the (entropy-regularised) optimal transport loss is typically available in closed form. This is in contrast to the primal problem, where evaluation of the objective and its gradients requires potentially many Sinkhorn iterations.
#
# Consider a general convex and unconstrained problem. Under (usually satisfied) conditions for strong duality to hold, we have
# ```math
# \begin{aligned}
# &\min_{\rho} \operatorname{OT}_{\varepsilon}(\rho_0, \rho) + \mathcal{F}(\rho)  \\
# &= \min_{\rho} \sup_{u}\left[\langle \rho, u \rangle - \operatorname{OT}^*_{\varepsilon}(\rho_0, u)\right] + \mathcal{F}(\rho)  \\
# &= \sup_{u} \min_{\rho} \langle \rho, u \rangle - \operatorname{OT}^*_{\varepsilon}(\rho_0, u) + \mathcal{F}(\rho) \\
# &= \sup_{u} - \operatorname{OT}^*_{\varepsilon}(\rho_0, u) + \min_{\rho} \langle \rho, u \rangle + \mathcal{F}(\rho) \\
# &= \sup_{u} - \operatorname{OT}^*_{\varepsilon}(\rho_0, u) - \sup_{\rho} \langle \rho, -u \rangle - \mathcal{F}(\rho) \\
# &= \sup_{u} - \operatorname{OT}^*_{\varepsilon}(\rho_0, u) - \mathcal{F}^*(-u).
# \end{aligned}
# ```
# Thus, the dual problem is
# ```math
# \min_{u} \operatorname{OT}^*_{\varepsilon}(\rho_0, u) + \mathcal{F}^*(-u).
# ```
#
# The upshot here is that $u \mapsto \operatorname{OT}^*_{\varepsilon}(\rho_0, u)$ and its gradient is available in closed form. This is a known result in the literature [^CP18].
#
# [^CP18]: Cuturi, Marco, and Gabriel Peyré. “Semi-Dual Regularized Optimal Transport.” ArXiv: Learning, 2018.
#
# The formulas we state below are lifted from statements in [^Z21].
#
# [^Z21]: Zhang, Stephen Y. “A Unified Framework for Non-Negative Matrix and Tensor Factorisations with a Smoothed Wasserstein Loss.” ArXiv: Machine Learning, 2021.
#
# ```math
# \begin{aligned}
# \operatorname{OT}^*_{\varepsilon}(\rho_0, u) &= -\varepsilon \left\langle \rho_0, \log\left( \dfrac{\rho_0}{K e^{u/\varepsilon}} \right) - 1\right\rangle, \\
# \nabla_u \operatorname{OT}^*_{\varepsilon}(\rho_0, u) &= K^\top \left( \dfrac{\rho_0}{K e^{u/\varepsilon}} \right) \odot e^{u/\varepsilon}.
# \end{aligned}
# ```
# At optimality, we can recover the primal optimal point $\rho^\star$ from the dual optimal point $u^\star$ following the formula
# ```math
# \rho^\star = e^{u^\star/\varepsilon} \odot K^\top \dfrac{\rho_0}{K e^{u^\star/\varepsilon}}.
# ```
#
# When $\mathcal{F}^*(\cdot)$ is also available in closed form (this is not always the case), the dual problem has a closed form objective and can generally be solved much more efficiently than the primal problem.
#
# In the setting of the Fokker-Planck and porous medium equations, the function $\mathcal{F}$ can be identified with
#
# ```math
# \mathcal{F}(\rho) = \tau \left[ \langle \psi, \rho \rangle + E_m(\rho) \right].
# ```
#
# A straightforward computation shows that
# ```math
# \mathcal{F}^*(u) = \tau E_m^*\left( \frac{u}{\tau}-\psi \right),
# ```
# where
# ```math
#     E_m^*(u) = \begin{cases}
#     \langle e^u, \mathbf{1} \rangle, & m = 1 \\
#     \sum_i \left[ \left( u_i + \frac{m}{m-1} \right) \left( \frac{m-1}{m} u_i + 1 \right)^{\frac{1}{m-1}} - \frac{1}{m-1} \left( \frac{m-1}{m} u_i + 1 \right)^{\frac{m}{m-1}} \right], & m > 1.
#     \end{cases}
# ```
# In particular, for $m = 2$ we have a simpler formula
# ```math
# E_2^*(u) = \left\| \frac{u}{2} + 1 \right\|_2^2
# ```
#
# We now implement $E_m^*$ for $m = 1, 2$.
E_dual(u, m::Val{1}) = sum(exp.(u))
function E_dual(u, m::Val{2})
    return dot(u / 2 .+ 1, u / 2 .+ 1)
end;
#
# So, the dual problem we are dealing with reads
# ```math
# \min_{u} \operatorname{OT}^*_{\varepsilon}(\rho_0, u) + \tau E_m^*\left( \frac{-u}{\tau}-\psi \right),
# ```
# and we can thus set up `G_dual_fpe`, the dual objective.
#
function G_dual_fpe(u, ρ0, τ, ε, K)
    return OptimalTransport.Dual.ot_entropic_semidual(ρ0, u, ε, K) +
           τ * E_dual(-u / τ - Ψ, Val(1))
end;
#
# Now we set up `step` as previously, except this time we need to convert from the optimal dual variable $u^\star$ to the primal variable $\rho^\star$. In the code, this is handled by `getprimal_ot_entropic_semidual`. We use `ReverseDiff` in this problem.
#
function step(ρ0, τ, ε, K, G)
    obj = u -> G(u, ρ0, τ, ε, K)
    opt = optimize(
        obj,
        (∇, u) -> ReverseDiff.gradient!(∇, obj, u),
        zeros(size(ρ0)),
        LBFGS(),
        Optim.Options(; iterations=250, g_tol=1e-6),
    )
    return OptimalTransport.Dual.getprimal_ot_entropic_semidual(
        ρ0, Optim.minimizer(opt), ε, K
    )
end;
#
# Now we can solve the dual problem as previously, and we note that the dual formulation is solved an order of magnitude faster than the primal formulation.
#
ρ = similar(ρ0, size(ρ0, 1), N)
ρ[:, 1] = ρ0
for i in 2:N
    ρ[:, i] = step(ρ[:, i - 1], τ, ε, K, G_dual_fpe)
end
colors = range(colorant"red"; stop=colorant"blue", length=N)
plot(
    support,
    ρ;
    title=raw"$F(\rho) = \langle \psi, \rho \rangle + \langle \rho, \log(\rho) \rangle$",
    palette=colors,
    legend=nothing,
)

# Setting `m = 2`, we can simulate instead the porous medium equation.
#
function G_dual_pme(u, ρ0, τ, ε, K)
    return OptimalTransport.Dual.ot_entropic_semidual(ρ0, u, ε, K) +
           τ * E_dual(-u / τ - Ψ, Val(2))
end
ρ = similar(ρ0, size(ρ0, 1), N)
ρ[:, 1] = ρ0
for i in 2:N
    @info i
    ρ[:, i] = step(ρ[:, i - 1], τ, ε, K, G_dual_pme)
end
colors = range(colorant"red"; stop=colorant"blue", length=N)
plot(
    support,
    ρ;
    title=raw"$F(\rho) = \langle \psi, \rho \rangle + \langle \rho, \rho - 1\rangle$",
    palette=colors,
    legend=nothing,
)


text:
#===============================================================================
# Test suite for PIDLib
#
# Author: Carl Cortright
# Date: 1/2/2017
#
===============================================================================#
workspace()
include("PID.jl")

using Base.Test
using PIDLib
using PyPlot

controller = PID(0.1, 0.9, 0.0, 70)
process(t) = 2.0t + 60

plantVal = process(0)
plantVals = Array(Float64, 25)
for i in 1:25
  plantVal = process(i) - Update(controller, plantVal)
  plantVals[i] = plantVal
  println(plantVal)
end

# Plot the results
x = linspace(0,1,25)
plot(x, plantVals)

max_issues_repo_name: carlcortright/Playground


text:
export align

struct State
    i1::Int
    j1::Int
    i2::Int
    j2::Int
    isfirst::Bool
end

isvalid(s::State) = s.i1 > 0 && s.i1 <= s.j1 && s.i2 > 0 && s.i2 <= s.j2

function match(t1::Vector{Int}, t2::Vector{Int})
    states = [State(1,length(t1),1,length(t2),true)]
    aligns = Tuple[]
    while !isempty(states)
        st = pop!(states)
        newst = getlcs(t1[st.i1:st.j1], t2[st.i2:st.j2], st.isfirst)
        isvalid(newst) || continue

        newst = (newst[1]+st.i1-1, newst.j1+st.i1-1, newst.i2+st.i2-1, newst.j2+st.i2-1)
        isvalid(newst) || continue
        push!(aligns, newst)

        left = State(st.i1, newst.i1-1, st.i2, newst.i2-1, false)
        right = State(newst.j1+1, st.j1, newst.j2+1, st.j2, true)
        isvalid(left) && push!(states, left)
        isvalid(right) && push!(states, right)
    end

    array = zeros(Int, length(t1))
    for a in aligns
        for i = 0:a.j1-a.i1
            array[a.i1+i] = a.i2 + i
        end
    end
    array
end

function align(t1::Vector{Int}, t2::Vector{Int})
    isvalid(st::NTuple{4,Int}) = st[1] > 0 && st[1] <= st[2] && st[3] > 0 && st[3] <= st[4]
    states = [((1,length(t1),1,length(t2)),true)]
    aligns = NTuple{4,Int}[]
    while !isempty(states)
        st,isfirst = pop!(states)
        newst = getlcs(t1[st[1]:st[2]], t2[st[3]:st[4]], isfirst)
        isvalid(newst) || continue

        newst = (newst[1]+st[1]-1, newst[2]+st[1]-1, newst[3]+st[3]-1, newst[4]+st[3]-1)
        isvalid(newst) || continue
        push!(aligns, newst)

        left = st[1], newst[1]-1, st[3], newst[3]-1
        right = newst[2]+1, st[2], newst[4]+1, st[4]
        isvalid(left) && push!(states, (left,false))
        isvalid(right) && push!(states, (right,true))
    end

    array = zeros(Int, length(t1))
    #dict = Dict{Int,Int}()
    for a in aligns
        for i = 0:a[2]-a[1]
            array[a[1]+i] = a[3] + i
            #dict[a[1]+i] = a[3]+i
        end
    end
    array
    #dict
end

function align(s1::String, s2::String)
    t1 = Int[Int(s1[i]) for i=1:length(s1)]
    t2 = Int[Int(s2[i]) for i=1:length(s2)]
    align(t1, t2)
end

function getlcs(t1::Vector{Int}, t2::Vector{Int}, isfirst::Bool, minlen=1)
    m = length(t1)
    text = copy(t1)
    append!(text, t2)
    sa = sais(text)
    lcps = lcparray(sa, text)

    maxinds, maxlcp = Int[], 0
    for i = 2:length(lcps)
        k1, k2 = sa[i-1], sa[i]
        (k1 <= m) == (k2 <= m) && continue
        lcp = min(lcps[i], m-min(k1,k2)+1)
        if lcp >= maxlcp
            lcp > maxlcp && (maxinds = Int[])
            push!(maxinds, i)
            maxlcp = lcp
        end
    end

    maxk1, maxk2 = 0, 0
    for i in maxinds
        k1, k2 = sa[i-1], sa[i]
        lcp = min(lcps[i], m-min(k1,k2)+1)
        for j = i-1:-1:2
            (sa[j-1] <= m) == (sa[i-1] <= m) || break
            lcps[j] < maxlcp && break
            if isfirst
                sa[j-1] < k1 && (k1 = sa[j-1])
            else
                sa[j-1] > k1 && (k1 = sa[j-1])
            end
        end
        for j = i+1:length(lcps)
            (sa[j] <= m) == (sa[i] <= m) || break
            lcps[j] < maxlcp && break
            if isfirst
                sa[j] < k2 && (k2 = sa[j])
            else
                sa[j] > k2 && (k2 = sa[j])
            end
        end
        k1, k2 = k1 > k2 ? (k2,k1) : (k1,k2)
        k2 -= m
        if isfirst
            if maxk1 == 0 || k1 < maxk1
                maxk1,maxk2 = k1,k2
            end
        else
            if maxk1 == 0 || k1 > maxk1
                maxk1,maxk2 = k1,k2
            end
        end
    end
    maxlcp >= minlen ? (maxk1,maxk1+maxlcp-1,maxk2,maxk2+maxlcp-1) : (0,0,0,0)
end

function getlcs2(t1::Vector{Int}, t2::Vector{Int}, isfirst::Bool, minlen=1)
    m = length(t1)
    text = copy(t1)
    append!(text, t2)
    sa = sais(text)
    lcps = lcparray(sa, text)

    maxk1, maxk2, maxlcp = 0, 0, 0
    for i = 2:length(lcps)
        k1, k2 = sa[i-1], sa[i]
        (k1 <= m) == (k2 <= m) && continue
        k1, k2 = k1 > k2 ? (k2,k1) : (k1,k2)
        lcp = min(lcps[i], m-k1+1)
        cond = isfirst ? lcp > maxlcp : lcp >= maxlcp
        cond || continue

        if cond
            maxk1, maxk2, maxlcp = k1, k2-m, lcp
        end
    end
    maxlcp >= minlen ? (maxk1,maxk1+maxlcp-1,maxk2,maxk2+maxlcp-1) : (0,0,0,0)
end


text:
t = Template(;
    user="Student2Pro",
    dir="C:/Users/Cloud Johnson/.julia/dev",
    authors="Alex",
    julia=v"1.5",
    plugins=[
        License(),
        Git(;name="CloudJohnson", email="cbtk_5ifth@qq.com", manifest=true),
        GitHubActions(; x86=true),
        Codecov(),
        Documenter{GitHubActions}(),
        TravisCI(),
    ],
)

#t("PkgName")


text:
module AGRN

using Logging
using Printf

include("config.jl")
include("chromosome.jl")
include("mutation.jl")
include("crossover.jl")
include("evolution.jl")

end


text:
#examples of using sac2jld2 functions to convert sac files to jld2 format.
using SeisIO, SeisConvert
# Saves individual sac file (in) to JLD2 file.
# The output file name will be the same as the input, with new extension *.jld2.

infile = "BK.HUMO.BHZ.sac"
sac2jld2(infile)
#

# """
# Saves individual sac file (in) to JLD2 file (out), which is explicitly specified.
# """

sac2jld2(infile,"BK.HUMO.BHZ_younameit.jld2")

# """
# Saves sac files based on a filelist generated using ls() in Julia to
# JLD2 file (out), which is explicitly specified.
# """

filelist = ls("sac/Raw/Event_2010_001/*.sac")

sac2jld2(filelist,"2010.001T00.00.00","Event_2010_001.jld2")

# """
# Searches for all sac files scanning throgh the sacdirlist (Array{String,1}). The sac files
#     in each sac directory are assigned to one time stamp based on timestamplist (Array{String,1}).
#     The two Arrays must be equal in size. The list can be generated from a file list (using ls in shell),
#     using a utility function in this package: readfilelist(infilelist::String).
#
# All sac files will be saved into one single JLD2 file (out), which is explicitly specified as outfile.
# """
sacdirlist = ls("sac/Raw/Event_201*")
timestamplist=["2010.001T00.00.00","2010.002T00.00.00","2011.041T00.00.00"]
sac2jld2(sacdirlist,timestamplist,"testoutput_3groups.jld2")


text:
 using Test
 using LogicCircuits

include("../helper/little_circuits.jl")

@testset "Nnf c2d IO test" begin

    c2d_example = """
        nnf 15 17 4
        L -3
        L -2
        L 1
        A 3 2 1 0
        L 3
        O 3 2 4 3
        L -4
        A 2 6 5
        L 4
        A 2 2 8
        A 2 1 4
        L 2
        O 2 2 11 10
        A 2 12 9
        O 4 2 13 7
        """

    function test_circuit(c)
        @test c isa LogicCircuit
        @test is⋁gate(c)
        @test num_nodes(c) == 15
        @test num_edges(c) == 17
        @test num_variables(c) == 4
    end

    circuit = parse(LogicCircuit, c2d_example, NnfFormat())

    test_circuit(circuit)

    mktempdir() do tmp
        temp_path = "$tmp/test.nnf"
        write(temp_path, circuit)

        circuit2 = read(temp_path, LogicCircuit)

        test_circuit(circuit2)

        @test prob_equiv(circuit, circuit2, 10)

        write("$temp_path.gz", circuit)

        circuit2 = read("$temp_path.gz", LogicCircuit)

        test_circuit(circuit2)

        @test prob_equiv(circuit, circuit2, 10)

    end
end

@testset "Nnf SDD IO test" begin

    circuit = readme_sdd()

    mktempdir() do tmp
        temp_path = "$tmp/test.nnf"
        write(temp_path, circuit)

        circuit2 = read(temp_path, LogicCircuit)

        @test num_nodes(circuit) == num_nodes(circuit2)
        @test num_edges(circuit) == num_edges(circuit2)
        @test prob_equiv(circuit, circuit2, 10)

        write("$temp_path.gz", circuit)

        circuit2 = read("$temp_path.gz", LogicCircuit)

        @test num_nodes(circuit) == num_nodes(circuit2)
        @test num_edges(circuit) == num_edges(circuit2)
        @test prob_equiv(circuit, circuit2, 10)
    end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using Plots
# gr()
data = rand(21, 100)
heatmap(1:size(data, 1),
    1:size(data,2), data,
    c = cgrad([:blue, :white, :red, :yellow]),
    xlabel = "x values", ylabel = "y values",
    title = "My title")

sbst = se[:, se.coldata.condition .== "ocean"]
data = assay(sbst, "sim")

heatmap(1:20, 1:size(sbst, 2), assay(sbst, "sim"),
        c = cgrad([:blue, :white, :red, :yellow]),
        xlabel = "x values", ylabel = "y values",
        title = "My title")


text:
export create_missing!

using Missings: disallowmissing

"""
	create_missing!(df, col::Symbol)

Create a new column for where `col` is missing
"""
create_missing!(df, col::Symbol; prefix="", suffix = "_missing") = begin
	df[!, prefix*string(col)*suffix] = ismissing.(df[!, col])
	if eltype(df[!, col]) <: Union{String, Missing}
		df[!, col] .= disallowmissing.(coalesce.(df[!, col], "JULIA.MISSING"))
	else
		df[!, col] .= disallowmissing.(coalesce.(df[!, col], zero(eltype(df[!, col]))))
	end
	df
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
"""
    struct ConstantChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}

Charge density model that assumes a constant charge density everywhere.

## Fields
* `ρ::T`: the constant value of the charge density.

## Definition in Configuration File

A `ConstantChargeDensity` is defined in the configuration file through the field `charge_density`
(of a `passive` or `surrounding`) with `name: constant` and a `value` for `ρ`.

An example definition of a constant charge density looks like this:
```yaml
charge_density:
  name: constant
  value: 1.0e-10 # C/m³
```
"""
struct ConstantChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}
    ρ::T
end

function get_charge_density(cdm::ConstantChargeDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    return cdm.ρ
end

function ChargeDensity(T::DataType, t::Val{:constant}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, input_units::NamedTuple)
    ρ::T = haskey(dict, "value") ? _parse_value(T, dict["value"], input_units.length^(-3)) : T(0)
    ConstantChargeDensity{T}( ρ )
end


text:
module RedirectRequest

import ..Layer, ..request
using ..URIs
using ..Messages
using ..Pairs: setkv
using ..Header
import ..@debug, ..DEBUG_LEVEL

"""
    request(RedirectLayer, method, ::URI, headers, body) -> HTTP.Response

Redirects the request in the case of 3xx response status.
"""

abstract type RedirectLayer{Next <: Layer} <: Layer end
export RedirectLayer

function request(::Type{RedirectLayer{Next}},
                 method::String, url::URI, headers, body;
                 redirect_limit=3, forwardheaders=false, kw...) where Next
    count = 0
    while true

        res = request(Next, method, url, headers, body; kw...)

        if (count == redirect_limit
        ||  !isredirect(res)
        ||  (location = header(res, "Location")) == "")
            return res
        end

        @static if VERSION > v"0.7.0-DEV.2338"
        kw = merge(merge(NamedTuple(), kw), (parent = res,))
        else
        setkv(kw, :parent, res)
        end
        url = absuri(location, url)
        if forwardheaders
            headers = filter(h->!(h[1] in ("Host", "Cookie")), headers)
        else
            headers = Header[]
        end

        @debug 1 "➡  Redirect: $url"

        count += 1
    end

    @assert false "Unreachable!"
end

end # module RedirectRequest


text:
using Test
using Cleaner: CleanTable, compact_columns, compact_rows, compact_rows!, compact_table,
    cols, size, names

@testset "Compact functions working fine" begin
    testRM1 = CleanTable(
        [:A, :B, :C], [[missing, missing, missing], [1, missing, 3], ["x", "", "z"]]
    )
    testRM2 = CleanTable(
        [:A, :B, :C], [["", missing, ""], [1, missing, 3], ["x", missing, "z"]]
    )

    @test size(compact_columns(testRM1)) == (3, 2)
    @test size(compact_rows(testRM2)) == (2, 3)
    @test size(compact_table(testRM1)) == (3, 2)

    @test size(compact_columns(testRM2; empty_values=[""])) == (3, 2)
    @test size(compact_rows(testRM1; empty_values=[""])) == (2, 3)
    @test size(compact_table(testRM2; empty_values=[""])) == (2, 2)

    @test names(compact_columns(testRM1)) == [:B, :C]

    if Threads.nthreads() > 1
        testRM1 = CleanTable([:A, :B], [collect(1:1_000_000), repeat([missing], 1_000_000)])
        testRM2 = CleanTable([:A, :B],
        [append!(Union{Missing, Int64}[missing], collect(1:1_000_000)),
        append!(Union{Missing, Int64}[missing], collect(1:1_000_000))])

        @test size(compact_columns(testRM1)) == (1_000_000, 1)
        @test size(compact_rows!(testRM2; empty_values=[1, 3, 10])) == (999_997, 2)
        @test !in(10, cols(testRM2)[1]) && !in(10, cols(testRM2)[2])
    end
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner