text:
function kFoldClassPropnPrep(X::Array, y::Vector, numFolds::Int)

  classes = sort(unique(y))                                     # Vector containing names of class outcomes
  numClasses = length(classes)                                  # Number of output classes
  numPerClass = Int[sum(y .== k) for k in classes]              # Vector containing number of instances of each class

  # Check that there are at least 'numFolds' instances of each class
  for k = 1:numClasses
      numPerClass[k] < numFolds && display("ERROR! Label $k has fewer than $numFolds (number of folds) samples.")
  end

  # FoldSizes will be a matrix holding the number of vectors to place in each fold for each category.
  # The number of folds may not divide evenly into the number of vectors so we need to distribute the remainder.
  FoldSizes = zeros(Int, numClasses, numFolds)                  # Initiate empty FoldSizes matrix
  for k = 1:numClasses                                          # For each class ...
    num_k = numPerClass[k]                                      #  record the number of instances in that class
    for f = 1:numFolds                                          #  Then for each fold ...
      foldSize = int(ceil(num_k / (numFolds - f + 1)))          #   calculate the allocation from this class to this fold
      FoldSizes[k, f] = foldSize                                #   record that number in the FoldSizes matrix
      num_k -= foldSize                                         #   reduce the number to be allocated to the remaining folds
    end
  end

  # Check the fold sizes sum up correctly by comparing to the total number for each class
  any(sum(FoldSizes, 2) .!= numPerClass) && display("ERROR! The sum of fold sizes does not equal the number of vectors.")

  # Randomly sort the matrix X and vector y and then group by class for proportional sampling
  m, n = size(X)                                                # m = # of examples, n = # of features
  randVecs = [X y][randperm(m), :]                              # Randomly shuffle the input vectors and output classes
  groupVecs = sortrows(randVecs, by=x->x[n + 1])                # Group the vectors by class for proportional sampling
  X_sorted = groupVecs[:, 1:n]                                  # Retrieve shuffled and sorted X
  y_sorted = groupVecs[:, n + 1]                                # Retrieve shuffled and sorted y

  # Calculate arrays containing starting and ending indices of each class for each fold
  cvStart = cumsum(numPerClass) .- numPerClass .+ 1 .+ cumsum(FoldSizes, 2) .- FoldSizes
  cvEnd = cvStart .+ FoldSizes - 1

  # Return everything needed for the k-fold CV loop
  return X_sorted, y_sorted, classes, FoldSizes, cvStart, cvEnd
end

function kFoldClassPropnArrays(X_sorted::Array, y_sorted::Vector, classes::Vector,
                               cvStart::Array{Int,2}, cvEnd::Array{Int,2}, roundNumber::Int)

  numClasses = length(classes)                                  # Number of output classes
  numVectors = length(y_sorted)                                 # Number of vectors

  cv_inds = Int[]                                               # Ititialise empty vector to contain cv set indices
  for k = 1 : numClasses                                        # For each class ...
    append!(cv_inds, cvStart[k, roundNumber] : cvEnd[k, roundNumber])    #  fetch the cv set indices for this round
  end
  train_inds = setdiff(1:numVectors, cv_inds)                   # Vector of training set indices != cv set indices

  X_train = X_sorted[train_inds, :]                             # Matrix containing training set vectors
  y_train = y_sorted[train_inds]                                # Vector containing training set classes
  X_cv = X_sorted[cv_inds, :]                                   # Matrix containing cross validation set vectors
  y_cv = y_sorted[cv_inds]                                      # Vector containing cross validation set classes

  return X_train, y_train, X_cv, y_cv                           # Return training and cross validation arrays
end

function shufflePartition2(X::Array, y::Vector, p_cv::Float64=0.2)

  # Function
  # Inputs
  # Outputs

  num_cv = int(length(y) * p_cv)                                # Number in cross validation set
  cv_inds = shuffle([1:length(y)] .<= num_cv)                   # Vector of shuffled indices

  X_train = X[!cv_inds, :]                                      # Matrix containing training set vectors
  y_train = y[!cv_inds]                                         # Vector containing training set classes
  X_cv = X[cv_inds, :]                                          # Matrix containing cross validation set vectors
  y_cv = y[cv_inds]                                             # Vector containing cross validation set classes

  return X_train, y_train, X_cv, y_cv                           # Return training and cross validation arrays
end

function shufflePartition3(X::Array, y::Vector, p_cv::Float64=0.2, p_test::Float64=0.2)

  num_cv = int(length(y) * p_cv)                                # Number in cross validation set
  num_test = int(length(y) * p_test)                            # Number in test set
  num_train = int(length(y)) - num_cv - num_test                # Number in training set

  # Vector of shuffled indices
  inds = shuffle([ones(Int, num_train); ones(Int, num_cv)*2; ones(Int, num_test)*3])

  X_train = X[inds .== 1, :]                                    # Matrix containing training set vectors
  y_train = y[inds .== 1]                                       # Vector containing training set classes
  X_cv = X[inds .== 2, :]                                       # Matrix containing cross validation set vectors
  y_cv = y[inds .== 2]                                          # Vector containing cross validation set classes
  X_test = X[inds .== 3, :]                                     # Matrix containing test set vectors
  y_test = y[inds .== 3]                                        # Vector containing test set classes

  return X_train, y_train, X_cv, y_cv, X_test, y_test           # Return training, cv and test arrays
end

function propnPartition2(X::Array, y::Vector, p_cv::Float64=0.2)

  classes = sort(unique(y))                                     # Vector containing names of class outcomes
  numClasses = length(classes)                                  # Number of output classes
  numPerClass = Int[sum(y .== k) for k in classes]              # Vector containing number of instances of each class
  num_cv = int(length(y) * p_cv)                                # Number in cross validation set

  # Randomly sort the matrix X and vector y and then group by class for proportional sampling
  m, n = size(X)                                                # m = # of examples, n = # of features
  randVecs = [X y][randperm(m), :]                              # Randomly shuffle the input Matrix and output Vector
  groupVecs = sortrows(randVecs, by=x->x[n + 1])                # Group the vectors by class for proportional sampling
  X_sorted = groupVecs[:, 1:n]                                  # Retrieve shuffled and sorted X Matrix
  y_sorted = groupVecs[:, n + 1]                                # Retrieve shuffled and sorted y Vector

  cv_size = int(numPerClass * p_cv)                             # Vector of size of each class in the cv set
  cv_start = cumsum(numPerClass) .- numPerClass + 1             # Vector of start indices of each class
  cv_end = cv_start .+ cv_size - 1                              # Vector of end indices in cv set of each class

  cv_inds = Int[]                                               # Ititialise empty vector to contain cv set indices
  for k = 1 : numClasses                                        # For each class ...
    append!(cv_inds, cv_start[k] : cv_end[k])                   #  fetch the cv set indices for this round
  end
  train_inds = setdiff(1:m, cv_inds)                            # Vector of training set indices != cv set indices

  X_train = X_sorted[train_inds, :]                             # Matrix containing training set vectors
  y_train = y_sorted[train_inds]                                # Vector containing training set classes
  X_cv = X_sorted[cv_inds, :]                                   # Matrix containing cross validation set vectors
  y_cv = y_sorted[cv_inds]                                      # Vector containing cross validation set classes

  return X_train, y_train, X_cv, y_cv                           # Return training and cross validation arrays
end


text:
import Base: size, getindex

"""
`DelayEmbedding` wraps any vector and exposes a view with vectors of
`dim` consecutive elements, with differences between indices being `lag`.
For `dim` equal to `1` it reduces to a regular `dim`-sized slidig window.
"""
immutable DelayEmbedding{T, Vec <: AbstractVector} <: AbstractVector{AbstractVector{T}}
  vec :: Vec
  dim :: Int
  lag :: Int

  function DelayEmbedding(vec :: Vec, dim :: Int, lag :: Int)
    if dim < 1
      throw(ArgumentError("Dimension has to be positive"))
    elseif lag < 1
      throw(ArgumentError("Lag has to be positive"))
    else
      new(vec, dim, lag)
    end
  end
end

function DelayEmbedding{Vec <: AbstractVector}(vec :: Vec, dim :: Int, lag :: Int = 1)
  DelayEmbedding{eltype(vec), Vec}(vec, dim, lag)
end

size(de :: DelayEmbedding) = max(0, size(de.vec, 1) - (de.dim - 1) * de.lag)

function getindex(de :: DelayEmbedding, i :: Int)
  slice(de.vec, i : de.lag : (i + (de.dim - 1) * de.lag))
end

ext: jl


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "ParamBitCircuit test" begin

    function test_integrity(pbc)
        bc = pbc.bitcircuit
        params = pbc.params
        @test num_elements(pbc) == size(bc.elements, 2)
        for el in 1:size(bc.elements,2)
            d = bc.elements[1,el]
            @test bc.nodes[1,d] <= el
            @test bc.nodes[2,d] >= el
            p = bc.elements[2,el]
            @test el ∈ bc.parents[bc.nodes[3,p]:bc.nodes[4,p]]
            s = bc.elements[3,el]
            @test el ∈ bc.parents[bc.nodes[3,s]:bc.nodes[4,s]]
        end
        for node in 1:size(bc.nodes,2)
            first_el = bc.nodes[1,node]
            last_el = bc.nodes[2,node]
            if first_el != 0
                for i = first_el:last_el
                    @test bc.elements[1,i] == node
                end
            end
            first_par = bc.nodes[3,node]
            last_par = bc.nodes[3,node]
            if first_par != 0
                for i = first_par:last_par
                    par = bc.parents[i]
                    @test bc.elements[2,par] == node || bc.elements[3,par] == node
                end
            else
                @test node == num_nodes(pbc) || node <= num_leafs(pbc)
            end
        end
        @test sum(length, bc.layers) == size(bc.nodes,2)
    end

    r = fully_factorized_circuit(PlainProbCircuit, 10)
    test_integrity(ParamBitCircuit(r, 10))

    r = fully_factorized_circuit(LogisticCircuit, 10; classes=2)
    test_integrity(ParamBitCircuit(r, 2, 10))

    vtree = PlainVtree(10, :balanced)
    r = fully_factorized_circuit(StructProbCircuit, vtree)
    test_integrity(ParamBitCircuit(r, 10))

end


text:
# Example model ZIPoisson provided by Chris Fisher

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking

import Distributions: logpdf, rand

struct ZIPoisson{T1,T2} <: DiscreteUnivariateDistribution
    logλ::T1
    w::T2
end

function logpdf(d::ZIPoisson, y::Int)
    LL = 0.0
    if y == 0
        LLs = zeros(typeof(d.logλ), 2)
        LLs[1] = log(d.w)
        LLs[2] = log(1 - d.w) - exp(d.logλ)
        LL = logsumexp(LLs)
    else
        LL = log(1 - d.w) + logpdf(LogPoisson(d.logλ), y)
    end
    return LL
end

function rand(d::ZIPoisson)
    return rand() <= d.w ? 0 : rand(Poisson(exp(d.logλ)))
end

rand(d::ZIPoisson, N::Int) = map(_->rand(d), 1:N)

inv_logit(x) = 1/(1 + exp(-x))

@model ppl12_3bt(data) = begin
    a1 ~ Normal(1, .5)
    logλ ~ Normal(-1.5, 1)
    w = inv_logit(a1)
    data .~ ZIPoisson(logλ, w)
end

#Random.seed!(74591)
data = rand(ZIPoisson(-1.5, .70), 1000)

n_samples = 3000
n_adapt = 1500
config = NUTS(n_adapt, .65)
n_chains = 4
chns12_3bt = sample(ppl12_3bt(data), config, MCMCThreads(), n_samples, n_chains, progress=true)
chns12_3bt |> display


text: AbstractDTM = Union{DTM}


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

"""
    CompleteSnapshot()

Seals and completes the snapshot after all of the required blocks of data have been written to it. Completing the snapshot changes the status to completed. You cannot write new blocks to a snapshot after it has been completed.

# Required Parameters
- `snapshotId`: The ID of the snapshot.
- `x-amz-ChangedBlocksCount`: The number of blocks that were written to the snapshot.

# Optional Parameters
- `x-amz-Checksum`: An aggregated Base-64 SHA256 checksum based on the checksums of each written block. To generate the aggregated checksum using the linear aggregation method, arrange the checksums for each written block in ascending order of their block index, concatenate them to form a single string, and then generate the checksum on the entire string using the SHA256 algorithm.
- `x-amz-Checksum-Aggregation-Method`: The aggregation method used to generate the checksum. Currently, the only supported aggregation method is LINEAR.
- `x-amz-Checksum-Algorithm`: The algorithm used to generate the checksum. Currently, the only supported algorithm is SHA256.
"""
complete_snapshot(snapshotId, x_amz_ChangedBlocksCount; aws_config::AWSConfig=global_aws_config()) = ebs("POST", "/snapshots/completion/$(snapshotId)", Dict{String, Any}("headers"=>Dict{String, Any}("x-amz-ChangedBlocksCount"=>x_amz_ChangedBlocksCount)); aws_config=aws_config)
complete_snapshot(snapshotId, x_amz_ChangedBlocksCount, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ebs("POST", "/snapshots/completion/$(snapshotId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("headers"=>Dict{String, Any}("x-amz-ChangedBlocksCount"=>x_amz_ChangedBlocksCount)), args)); aws_config=aws_config)

"""
    GetSnapshotBlock()

Returns the data in a block in an Amazon Elastic Block Store snapshot.

# Required Parameters
- `blockIndex`: The block index of the block from which to get data. Obtain the BlockIndex by running the ListChangedBlocks or ListSnapshotBlocks operations.
- `blockToken`: The block token of the block from which to get data. Obtain the BlockToken by running the ListChangedBlocks or ListSnapshotBlocks operations.
- `snapshotId`: The ID of the snapshot containing the block from which to get data.

"""
get_snapshot_block(blockIndex, blockToken, snapshotId; aws_config::AWSConfig=global_aws_config()) = ebs("GET", "/snapshots/$(snapshotId)/blocks/$(blockIndex)", Dict{String, Any}("blockToken"=>blockToken); aws_config=aws_config)
get_snapshot_block(blockIndex, blockToken, snapshotId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ebs("GET", "/snapshots/$(snapshotId)/blocks/$(blockIndex)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("blockToken"=>blockToken), args)); aws_config=aws_config)

"""
    ListChangedBlocks()

Returns information about the blocks that are different between two Amazon Elastic Block Store snapshots of the same volume/snapshot lineage.

# Required Parameters
- `secondSnapshotId`: The ID of the second snapshot to use for the comparison.  The SecondSnapshotId parameter must be specified with a FirstSnapshotID parameter; otherwise, an error occurs.

# Optional Parameters
- `firstSnapshotId`: The ID of the first snapshot to use for the comparison.  The FirstSnapshotID parameter must be specified with a SecondSnapshotId parameter; otherwise, an error occurs.
- `maxResults`: The number of results to return.
- `pageToken`: The token to request the next page of results.
- `startingBlockIndex`: The block index from which the comparison should start. The list in the response will start from this block index or the next valid block index in the snapshots.
"""
list_changed_blocks(secondSnapshotId; aws_config::AWSConfig=global_aws_config()) = ebs("GET", "/snapshots/$(secondSnapshotId)/changedblocks"; aws_config=aws_config)
list_changed_blocks(secondSnapshotId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ebs("GET", "/snapshots/$(secondSnapshotId)/changedblocks", args; aws_config=aws_config)

"""
    ListSnapshotBlocks()

Returns information about the blocks in an Amazon Elastic Block Store snapshot.

# Required Parameters
- `snapshotId`: The ID of the snapshot from which to get block indexes and block tokens.

# Optional Parameters
- `maxResults`: The number of results to return.
- `pageToken`: The token to request the next page of results.
- `startingBlockIndex`: The block index from which the list should start. The list in the response will start from this block index or the next valid block index in the snapshot.
"""
list_snapshot_blocks(snapshotId; aws_config::AWSConfig=global_aws_config()) = ebs("GET", "/snapshots/$(snapshotId)/blocks"; aws_config=aws_config)
list_snapshot_blocks(snapshotId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ebs("GET", "/snapshots/$(snapshotId)/blocks", args; aws_config=aws_config)

"""
    PutSnapshotBlock()

Writes a block of data to a snapshot. If the specified block contains data, the existing data is overwritten. The target snapshot must be in the pending state. Data written to a snapshot must be aligned with 512-byte sectors.

# Required Parameters
- `BlockData`: The data to write to the block. The block data is not signed as part of the Signature Version 4 signing process. As a result, you must generate and provide a Base64-encoded SHA256 checksum for the block data using the x-amz-Checksum header. Also, you must specify the checksum algorithm using the x-amz-Checksum-Algorithm header. The checksum that you provide is part of the Signature Version 4 signing process. It is validated against a checksum generated by Amazon EBS to ensure the validity and authenticity of the data. If the checksums do not correspond, the request fails. For more information, see  Using checksums with the EBS direct APIs in the Amazon Elastic Compute Cloud User Guide.
- `blockIndex`: The block index of the block in which to write the data. A block index is a logical index in units of 512 KiB blocks. To identify the block index, divide the logical offset of the data in the logical volume by the block size (logical offset of data/524288). The logical offset of the data must be 512 KiB aligned.
- `snapshotId`: The ID of the snapshot.
- `x-amz-Checksum`: A Base64-encoded SHA256 checksum of the data. Only SHA256 checksums are supported.
- `x-amz-Checksum-Algorithm`: The algorithm used to generate the checksum. Currently, the only supported algorithm is SHA256.
- `x-amz-Data-Length`: The size of the data to write to the block, in bytes. Currently, the only supported size is 524288. Valid values: 524288

# Optional Parameters
- `x-amz-Progress`: The progress of the write process, as a percentage.
"""
put_snapshot_block(BlockData, blockIndex, snapshotId, x_amz_Checksum, x_amz_Checksum_Algorithm, x_amz_Data_Length; aws_config::AWSConfig=global_aws_config()) = ebs("PUT", "/snapshots/$(snapshotId)/blocks/$(blockIndex)", Dict{String, Any}("BlockData"=>BlockData, "headers"=>Dict{String, Any}("x-amz-Checksum"=>x_amz_Checksum, "x-amz-Checksum-Algorithm"=>x_amz_Checksum_Algorithm, "x-amz-Data-Length"=>x_amz_Data_Length)); aws_config=aws_config)
put_snapshot_block(BlockData, blockIndex, snapshotId, x_amz_Checksum, x_amz_Checksum_Algorithm, x_amz_Data_Length, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ebs("PUT", "/snapshots/$(snapshotId)/blocks/$(blockIndex)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("BlockData"=>BlockData, "headers"=>Dict{String, Any}("x-amz-Checksum"=>x_amz_Checksum, "x-amz-Checksum-Algorithm"=>x_amz_Checksum_Algorithm, "x-amz-Data-Length"=>x_amz_Data_Length)), args)); aws_config=aws_config)

"""
    StartSnapshot()

Creates a new Amazon EBS snapshot. The new snapshot enters the pending state after the request completes.  After creating the snapshot, use  PutSnapshotBlock to write blocks of data to the snapshot.

# Required Parameters
- `VolumeSize`: The size of the volume, in GiB. The maximum size is 16384 GiB (16 TiB).

# Optional Parameters
- `ClientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request. Idempotency ensures that an API request completes only once. With an idempotent request, if the original request completes successfully. The subsequent retries with the same client token return the result from the original successful request and they have no additional effect. If you do not specify a client token, one is automatically generated by the AWS SDK. For more information, see  Idempotency for StartSnapshot API in the Amazon Elastic Compute Cloud User Guide.
- `Description`: A description for the snapshot.
- `Encrypted`: Indicates whether to encrypt the snapshot. To create an encrypted snapshot, specify true. To create an unencrypted snapshot, omit this parameter. If you specify a value for ParentSnapshotId, omit this parameter. If you specify true, the snapshot is encrypted using the CMK specified using the KmsKeyArn parameter. If no value is specified for KmsKeyArn, the default CMK for your account is used. If no default CMK has been specified for your account, the AWS managed CMK is used. To set a default CMK for your account, use  ModifyEbsDefaultKmsKeyId. If your account is enabled for encryption by default, you cannot set this parameter to false. In this case, you can omit this parameter. For more information, see  Using encryption in the Amazon Elastic Compute Cloud User Guide.
- `KmsKeyArn`: The Amazon Resource Name (ARN) of the AWS Key Management Service (AWS KMS) customer master key (CMK) to be used to encrypt the snapshot. If you do not specify a CMK, the default AWS managed CMK is used. If you specify a ParentSnapshotId, omit this parameter; the snapshot will be encrypted using the same CMK that was used to encrypt the parent snapshot. If Encrypted is set to true, you must specify a CMK ARN.
- `ParentSnapshotId`: The ID of the parent snapshot. If there is no parent snapshot, or if you are creating the first snapshot for an on-premises volume, omit this parameter. If your account is enabled for encryption by default, you cannot use an unencrypted snapshot as a parent snapshot. You must first create an encrypted copy of the parent snapshot using CopySnapshot.
- `Tags`: The tags to apply to the snapshot.
- `Timeout`: The amount of time (in minutes) after which the snapshot is automatically cancelled if:   No blocks are written to the snapshot.   The snapshot is not completed after writing the last block of data.   If no value is specified, the timeout defaults to 60 minutes.
"""
start_snapshot(VolumeSize; aws_config::AWSConfig=global_aws_config()) = ebs("POST", "/snapshots", Dict{String, Any}("VolumeSize"=>VolumeSize, "ClientToken"=>string(uuid4())); aws_config=aws_config)
start_snapshot(VolumeSize, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ebs("POST", "/snapshots", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("VolumeSize"=>VolumeSize, "ClientToken"=>string(uuid4())), args)); aws_config=aws_config)


text:
###
### Conversion & Promotion
###
###
### Conversion methods for LongSequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

###
### Promotion
###
for alph in (DNAAlphabet, RNAAlphabet)
    @eval function Base.promote_rule(::Type{LongSequence{A}}, ::Type{LongSequence{B}}) where {A<:$alph,B<:$alph}
        return LongSequence{promote_rule(A, B)}
    end
end

###
### Conversion
###

function Base.convert(::Type{T}, seq::LongSequence{<:NucleicAcidAlphabet}) where
         {T<:LongSequence{<:NucleicAcidAlphabet}}
    return T(seq)
end

# Convert from a LongSequence to to a DNA or RNA vector
function Base.convert(::Type{LongSequence{A}}, seq::Vector) where A<:Alphabet
    return LongSequence{A}(seq)
end

Base.parse(::Type{LongSequence{A}}, seq::AbstractString) where A = LongSequence{A}(seq)


text:
function markminima(J, tol=1e-4; debug=false)
    N = length(J)
    M = zeros(Bool,N)
    P = zeros(N)
    # find all minima
    M[1] = J[1] < J[2]
    prev = J[2] - J[1]
    curr = 0.0
    P[1] = prev
    for i=2:(N-1)
        curr = J[i+1] - J[i]
        curr = (abs(curr) < tol) ? 0 : curr
        debug && println("$i: $prev<=0 && $curr>=0")
        M[i] = prev<=0 && curr>=0
        prev=curr
        P[i] = prev
    end
    M[end] = J[end] < J[end-1]
    debug && println(M)
    debug && println(P[M])
    return (M, P)
end

function findglobalmin(J, tol=1e-4; debug=false)
    N = length(J)
    N == 1 && return 1
    # find all minima
    M, P = markminima(J; debug=false)
    # zero minima below tolerance (skip first & last)
    minima, MI = if sum(M[2:end-1]) > 0
        Mnew = copy(M)
        Mnew[1] = Mnew[end] = false
        mval, mi = findmin(J[Mnew])
        map(v->(abs(v) < tol) ? 0 : v, J[Mnew].-mval), Mnew
    else
        J[M], M
    end
    debug && println(minima)
    # get minima index
    findall(MI)[(last∘findmin)(minima)]
end

function findglobalmin2(J; debug=false)
    N = length(J)
    N == 1 && return 1
    f(x) = J[round(Int,x)]
    r = optimize(f, 1, N)
    debug && println(r)
    idxs = sortperm(J)
    si = findfirst(x-> x>=minimum(r), J[idxs])
    si === nothing ? N : si
    idxs[si]
end

function findlocalmin(J)
    findmin(J)[2]
end

function findlocalminrev(J)
    c = length(J)
    Jmax = Inf
    midx = c
    for i in c:-1:1
        if J[c] < Jmax
            Jmax = J[c]
        else
            break
        end
        midx = i
    end
    return idx
end

function findfirstminimum(J; debug=false)
    N = length(J)
    N == 1 && return 1
    # find all minima
    M, P = markminima(J; debug=debug)
    m, midx = Inf, 0
    for i in findall(M)
        J[i] |> println
        if J[i] < m
            m = J[i]
            midx = i
        else
            break
        end
    end
    return midx
end

function mixlogpdf(logps::AbstractMatrix{T}, ps::Vector{Vector{T}},
                   idxs::Vector{Vector{Int}}) where {T <: AbstractFloat}
    # using the formula below for numerical stability
    #
    # logpdf(d, x) = log(sum_i pri[i] * pdf(cs[i], x))
    #              = log(sum_i pri[i] * exp(logpdf(cs[i], x)))
    #              = log(sum_i exp(logpri[i] + logpdf(cs[i], x)))
    #              = m + log(sum_i exp(logpri[i] + logpdf(cs[i], x) - m))
    #
    #  m is chosen to be the maximum of logpri[i] + logpdf(cs[i], x)
    #  such that the argument of exp is in a reasonable range
    #
    M = length(idxs)
    N, L = size(logps)

    @assert all(e->length(e[1]) == length(e[2]), zip(idxs, ps)) "Sizes of indexes and priors must be equal"

    lps = Matrix{T}(undef, N, M)

    @inbounds for k in eachindex(idxs)
        p = ps[k]
        ids = idxs[k]
        K = length(p)
        lp = Vector{T}(undef, K)

        j = 1
        for logp in eachrow(logps)
            m = -Inf   # m <- the maximum of log(p(cs[i], x)) + log(pri[i])
            for i in eachindex(p)
                pi = p[i]
                if pi > 0.0
                    # lp[i] <- log(p(cs[i], x)) + log(pri[i])
                    lp_i = logp[ids[i]] + log(pi)
                    lp[i] = lp_i
                    if lp_i > m
                        m = lp_i
                    end
                end
            end
            v = 0.0
            for i = 1:K
                if p[i] > 0.0
                    v += exp(lp[i] - m)
                end
                lp[i] = 0
            end
            lps[j, k] = m + log(v)
            j+=1
        end
    end
    return lps
end

function score(fn, cl::C, data::AbstractMatrix) where {C <: ClusteringResult}
    assign = assignments(cl)
    clusters = [[i] for i in 1:nclusters(cl)]
    models = [model(data, findall(i .∈  assign)) for i in 1:nclusters(cl)]
    mcr = ModelClusteringResult(models, assign)
    fn(clusters, mcr, data)
end
score(fn, cl::PLMClusteringResult, data::AbstractMatrix) = fn(cl.clusters, cl.models, data)

max_stars_repo_name: wildart/PLMC

max_issues_repo_name: wildart/PLMC

max_forks_repo_name: wildart/PLMC


text:
"""
    FeatureDescriptor

A feature, in this context, refers to a distinct attribute manifested
by an atomic system.
A FeatureDescriptor describes a feature/class of features - i.e., its name,
possible values, etc.

A FeatureDescriptor does NOT store any actual instances of the value(s) of the
feature it describes.
Simply put, it can be understood to be "features of a feature".

All FeatureDescriptors MUST also describe an encoding and decoding scheme.
This can (and should) be easily done using a Codec.
"""
module FeatureDescriptor

using Base: Int16
using ..ChemistryFeaturization.AbstractType:
    AbstractAtoms, AbstractFeatureDescriptor, AbstractCodec
using ..ChemistryFeaturization.Codec: OneHotOneCold

import ..ChemistryFeaturization.encodable_elements
encodable_elements(fd::AbstractFeatureDescriptor) =
    throw(MethodError(encodable_elements, fd))
export encodable_elements

"""
    get_value(fd::AbstractFeatureDescriptor, atoms::AbstractAtoms)
Get the value(s) of feature corresponding to feature descriptor `fd` for structure `atoms`.
This function computes and returns the value that would actually get encoded by [`encode`](@ref).
"""
get_value(fd::AbstractFeatureDescriptor, atoms::AbstractAtoms) =
    throw(MethodError(fd, atoms))
(fd::AbstractFeatureDescriptor)(atoms::AbstractAtoms) = get_value(fd, atoms)

import ..ChemistryFeaturization.encode
"""
    encode(fd::AbstractFeatureDescriptor, atoms::AbstractAtoms)
Encode features for `atoms` using the feature descriptor `fd`.
"""
encode(fd::AbstractFeatureDescriptor, atoms::AbstractAtoms) =
    encode(fd.encoder_decoder, get_value(fd, atoms))
export encode
export get_value

import ..ChemistryFeaturization.decode
"""
    decode(fd::AbstractFeatureDescriptor, encoded_feature)
Decode `encoded_feature` using the feature descriptor `fd`.
"""
decode(fd::AbstractFeatureDescriptor, encoded_feature) =
    decode(fd.encoder_decoder, encoded_feature)
export decode

include("abstractfeatures.jl")
encode(efd::AbstractAtomFeatureDescriptor, atoms::AbstractAtoms) =
    hcat(encode(efd.encoder_decoder, get_value(efd, atoms))...)

import ..ChemistryFeaturization.output_shape
output_shape(afd::AbstractFeatureDescriptor) = output_shape(afd, afd.encoder_decoder)
output_shape(::AbstractAtomFeatureDescriptor, ed::OneHotOneCold) = output_shape(ed)
export output_shape

include("bondfeatures.jl")

include("elementfeature.jl")
export ElementFeatureDescriptor, encode, decode, output_shape

include("orbitalfeature.jl")
export OrbitalFeatureDescriptor

include("pairfeature.jl")
export PairFeatureDescriptor

include("speciesfeature.jl")
export SpeciesFeatureDescriptor

end


text:
function transpose_strings(input::AbstractArray)

end


text:
"""
    Ciphertext

A ciphertext element, consisting of two or more polynomials. It can be created from a `Plaintext`
element by encrypting it with an appropriate `Encryptor` instance. `Ciphertext` instances should
usually not be modified directly by the user but only through the corresponding functions of
`Evaluator`. Decryption is performed via a `Decryptor` instance, which converts a `Ciphertext` back
to a `Plaintext` instance.

See also: [`Plaintext`](@ref), [`Encryptor`](@ref), [`Decryptor`](@ref), [`Evaluator`](@ref)
"""
mutable struct Ciphertext <: SEALObject
  handle::Ptr{Cvoid}

  function Ciphertext()
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Ciphertext_Create1, libsealc), Clong,
                   (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   C_NULL, handleref)
    @check_return_value retval
    return Ciphertext(handleref[])
  end

  function Ciphertext(context)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Ciphertext_Create3, libsealc), Clong,
                   (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   context, C_NULL, handleref)
    @check_return_value retval
    return Ciphertext(handleref[])
  end

  function Ciphertext(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::Ciphertext)
  if isallocated(object)
    @check_return_value ccall((:Ciphertext_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function scale(encrypted::Ciphertext)
  value = Ref{Cdouble}(0)
  retval = ccall((:Ciphertext_Scale, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Cdouble}),
                 encrypted, value)
  @check_return_value retval
  return Float64(value[])
end

function scale!(encrypted::Ciphertext, value::Float64)
  retval = ccall((:Ciphertext_SetScale, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Cdouble}),
                 encrypted, value)
  @check_return_value retval
  return encrypted
end

function parms_id(encrypted::Ciphertext)
  parms_id_ = zeros(UInt64, 4)
  retval = ccall((:Ciphertext_ParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 encrypted, parms_id_)
  @check_return_value retval
  return parms_id_
end

function Base.size(encrypted::Ciphertext)
  sizeref = Ref{UInt64}(0)
  retval = ccall((:Ciphertext_Size, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 encrypted, sizeref)
  @check_return_value retval
  return (Int(sizeref[]),)
end
Base.length(encrypted::Ciphertext) = size(encrypted)[1]

function save_size(compr_mode, encrypted::Ciphertext)
  result = Ref{Int64}(0)
  retval = ccall((:Ciphertext_SaveSize, libsealc), Clong,
                 (Ptr{Cvoid}, UInt8, Ref{Int64}),
                 encrypted, compr_mode, result)
  @check_return_value retval
  return Int(result[])
end
save_size(encrypted::Ciphertext) = save_size(ComprModeType.default, encrypted)

function save!(buffer::DenseVector{UInt8}, length::Integer,
               compr_mode::ComprModeType.ComprModeTypeEnum, encrypted::Ciphertext)
  out_bytes = Ref{Int64}(0)
  retval = ccall((:Ciphertext_Save, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}, UInt64, UInt8, Ref{Int64}),
                 encrypted, buffer, length, compr_mode, out_bytes)
  @check_return_value retval
  return Int(out_bytes[])
end
function save!(buffer::DenseVector{UInt8}, length::Integer, encrypted::Ciphertext)
  return save!(buffer, length, ComprModeType.default, encrypted)
end
function save!(buffer::DenseVector{UInt8}, encrypted::Ciphertext)
  return save!(buffer, length(buffer), encrypted)
end

function load!(encrypted::Ciphertext, context::SEALContext, buffer::DenseVector{UInt8}, length)
  in_bytes = Ref{Int64}(0)
  retval = ccall((:Ciphertext_Load, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt8}, UInt64, Ref{Int64}),
                 encrypted, context, buffer, length, in_bytes)
  @check_return_value retval
  return Int(in_bytes[])
end
function load!(encrypted::Ciphertext, context::SEALContext, buffer::DenseVector{UInt8})
  return load!(encrypted, context, buffer, length(buffer))
end

function reserve!(encrypted::Ciphertext, capacity)
  retval = ccall((:Ciphertext_Reserve3, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64),
                 encrypted, capacity)
  @check_return_value retval
  return encrypted
end


text:
module MLLabelUtils

using StatsBase
using LearnBase
using MappedArrays
import Base.Broadcast: broadcastable

export

    ind2label,
    label2ind,

    labeltype,

    label,
    nlabel,

    poslabel,
    neglabel,
    isposlabel,
    isneglabel,

    ObsDim,
    classify,
    classify!,

    convertlabel,
#    convertlabel!,

    labelmap,
    labelmap!,
    labelfreq,
    labelfreq!,
    labelmap2vec,

    LabelEnc,
    labelenc,
    islabelenc,

    convertlabelview

include("learnbase.jl")
include("labelencoding.jl")
include("classify.jl")
include("convertlabel.jl")
include("labelmap.jl")

end # module


text:
using Printf

# isolate eigvec tests for use with special cases
function vectest(A::Matrix{T}, S::Schur{T2}, vtol; normal=false) where {T<:Complex, T2<:Complex}
    n = size(A,1)
    ulp = eps(real(T))
    VR = eigvecs(S)
    if verbosity[] > 1
        vec_err = norm(A * VR - VR * diagm(0 => S.values)) / (n * norm(A) * ulp)
        println("r.eigenvector error: $vec_err, vtol = $vtol")
        @test vec_err < 1000
    else
        @test norm(A * VR - VR * diagm(0 => S.values)) / (n * norm(A) * ulp) < vtol
    end
    VL = eigvecs(S, left=true)
    if verbosity[] > 1
        vec_err = norm(A' * VL - VL * diagm(0 => conj.(S.values))) / (n * norm(A) * ulp)
        println("l.eigenvector error: $vec_err, vtol = $vtol")
        @test vec_err < 1000
    else
        @test norm(A' * VL - VL * diagm(0 => conj.(S.values))) / (n * norm(A) * ulp) < vtol
    end
    if normal
        # check orthonormality of vectors where appropriate
        if (T == ComplexF16) && (n > 10)
            # is this simply too many ops for this type?
            @test_broken norm(VR*VR'-I) / (n * ulp) < vtol
            @test_broken norm(VL*VL'-I) / (n * ulp) < vtol
        else
            @test norm(VR*VR'-I) / (n * ulp) < vtol
            @test norm(VL*VL'-I) / (n * ulp) < vtol
        end
    end
end

function schurtest(A::Matrix{T}, tol; normal=false) where {T<:Complex}
    n = size(A,1)
    ulp = eps(real(T))
    # schur() uses eigtype(), which promotes ComplexF16 to ComplexF32
    if (real(T) <: BlasFloat) || (LinearAlgebra.eigtype(T) != T)
        S = GenericSchur.gschur(A)
    else
        S = schur(A)
    end
    # test 1: S.T is upper triangular
    @test all(tril(S.T,-1) .== 0)
    # test 2: norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    if verbosity[] > 1
        decomp_err = norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp)
        println("decomp. error: ", decomp_err)
        @test decomp_err < tol
    else
        @test norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    end
    # test 3: S.Z is unitary: norm(I - S.Z * S.Z') / (n * ulp) < tol
    @test norm(I - S.Z * S.Z') / (n * ulp) < tol
    # test 4: S.values are e.v. of T
    @test all(csort(S.values) .== csort(diag(S.T)))

    # It is tempting to check eigenvalues against LAPACK eigvals(A)
    # for suitable types, but the comparison is misleading without
    # detailed analysis because many of the test matrices are severly
    # non-normal.

    if normal
        # verify that Schur diagonalizes normal matrices
        @test norm(triu(S.T,1)) / (n * norm(A) * ulp) < tol
    end
    vtol = tol
    vectest(A, S, tol, normal=normal)
end

function hesstest(A::Matrix{T}, tol) where {T<:Complex}
    n = size(A,1)
    ulp = eps(real(T))
    H = invoke(GenericSchur._hessenberg!,
               Tuple{StridedMatrix{Ty},} where Ty,
               copy(A))
    Q = GenericSchur._materializeQ(H)
    if VERSION < v"1.3"
        HH = triu(H.data, -1)
    else
        HH = H.H
        # test 1: H.H is upper triangular
        @test all(tril(HH,-2) .== 0)
    end
    # test 2: norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    decomp_err = norm(A - Q * HH * Q') / (n * norm(A) * ulp)
    @test decomp_err < tol
    # test 3: S.Z is orthogonal: norm(I - S.Z * S.Z') / (n * ulp) < tol
    orth_err = norm(I - Q * Q') / (n * ulp)
    @test orth_err < tol
    # test 4: subdiagonal is real
    @test all(isreal.(diag(HH,-1)))
end

function hesstest(A::Hermitian{T}, tol) where {T<:Complex}
    n = size(A,1)
    ulp = eps(real(T))
    H = invoke(GenericSchur._hessenberg!,
               Tuple{Hermitian{Ty},} where {Ty <: Complex},
               copy(A))
    Q = GenericSchur._materializeQ(H)
    HH = H.H
    # test 1: H.H is formally symtridiagonal
    @test isa(HH, SymTridiagonal)
    @test eltype(HH) == real(eltype(A))
    # test 2: norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    decomp_err = norm(A - Q * HH * Q') / (n * norm(A) * ulp)
    @test decomp_err < tol
    # test 3: S.Z is orthogonal: norm(I - S.Z * S.Z') / (n * ulp) < tol
    orth_err = norm(I - Q * Q') / (n * ulp)
    @test orth_err < tol
end

"""
generate a random unitary matrix

Unitary matrices are normal, so Schur decomposition is diagonal for them.
(This is not only another test, but sometimes actually useful.)
"""
function randu(::Type{T},n) where {T<:Complex}
    if real(T) ∈ [Float16,Float32,Float64]
        A = randn(T,n,n)
        F = qr(A)
        return Matrix(F.Q) * Diagonal(sign.(diag(F.R)))
    else
        # don't have normal deviates for other types, but need appropriate QR
        A = randn(ComplexF64,n,n)
        F = qr(convert.(T,A))
        return Matrix(F.Q) * Diagonal(sign.(diag(F.R)))
    end
end

Random.seed!(1234)

@testset "Hessenberg $T" for T in [ComplexF64, Complex{Float16}]
    n = 32
    tol = 10
    A = rand(T,n,n)
    hesstest(A, tol)
    s = 100 * floatmin(real(T))
    Asmall = s * A
    hesstest(Asmall, tol)
end

if VERSION >= v"1.3"
    @testset "Hermitian Hessenberg $T" for T in [ComplexF64, Complex{Float16}]
        n = 32
        tol = 10
        A = Hermitian(rand(T,n,n), :U)
        hesstest(A, tol)
        s = 100 * floatmin(real(T))
        Asmall = s * A
        hesstest(Asmall, tol)
        A = Hermitian(rand(T,n,n), :L)
        hesstest(A, tol)
        s = 100 * floatmin(real(T))
        Asmall = s * A
        hesstest(Asmall, tol)
    end
end

@testset "group $T" for T in [Complex{BigFloat},Complex{Float16}]
unfl = floatmin(real(T))
ovfl = one(real(T)) / unfl
ulp = eps(real(T))
ulpinv = one(real(T)) / ulp
rtulp = sqrt(ulp)
rtulpi = one(real(T)) / rtulp

ens = [4,32]
@testset "random general" begin
    for n in ens
        A = rand(T,n,n)
        schurtest(A,10)
    end
end
@testset "random unitary" begin
    for n in ens
        A = randu(T,n)
        schurtest(A,20,normal=true)
    end
end

end # group testset

@testset "Godunov (complex)" begin
    # Set precision fine enough for eigval condition to apply,
    # but not fine enough to make it a slam-dunk.
    setprecision(BigFloat, 80) do
        A,v,econd = godunov(Complex{BigFloat})
        S = schur(A)
        δ = norm(S.Z*S.T*S.Z'-A)
        @test δ < 100 * eps(big(1.0)) * norm(A)
        t = 3*δ*econd
        @test isapprox(csort(S.values),v,atol=t)
        vnew = eigvals(A)
        @test isapprox(csort(vnew),v,atol=t)
    end
end

# v0.4.0 breaks on this matrix
function _example1(T=Complex{Float32})
    iseed = [4066,2905,502,2389]
    ulpinv = 1 / eps(real(T))
    A,d,ds = latme!(zeros(T,5,5),1.0,4,ulpinv,1.0,iseed=iseed)
    return A
end
@testset "short QR sweep" begin
    A = _example1(ComplexF32)
    schurtest(A,20)
end

@testset "group $T" for T in [ComplexF64, ComplexF32]

unfl = floatmin(real(T))
ovfl = one(real(T)) / unfl
ulp = eps(real(T))
ulpinv = one(real(T)) / ulp
rtulp = sqrt(ulp)
rtulpi = one(real(T)) / rtulp

# sizes of matrices to test
ens = [4,32]

itypemax = 10
# On author's dev system, tests pass with tol=5 for most and 6 for all,
# but generalization is risky.
# (LAPACK uses 20, FWIW.)
tols = fill(20,itypemax)

simrconds = [1.0, rtulpi, 0.0]
magns = [1.0, ovfl*ulp, unfl*ulpinv]

# test matrix from zdrves:
# conds=0 means index 3 for simrconds (i.e. irrelevant)
# type magn mode conds
# 1 1 0 0
# 2 1 0 0
# 3 1 0 0
# 4 1 4 0
# 4 1 3 0
# 4 1 1 0
# 4 2 4 0
# 4 3 4 0
# 6 1 4 1
# 6 1 3 1
# 6 1 1 1
# 6 1 5 1
# 6 1 4 2
# 6 1 3 2
# 6 1 1 2
# 6 1 5 2
# 6 2 5 2
# 6 3 5 2
# 9 1 4 0
# 9 2 3 0
# 9 3 1 0

# Other stuff in zdrves: pointless or not even used
# 1- zero
# 2- identity
# 4- diagonal, ev specified
# 5- symmetric, ev specified, not used in zdrves

@testset "Jordan" begin
    itype=3
    kmagn = [1]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
            anorm = magns[kmagn[j]]
            A = diagm(1 => ones(T,n-1)) + anorm * I
            schurtest(A,tols[itype])
        end
    end
end

@testset "general, ev specified" begin
    itype=6
    kmagn =  [1,1,1,1,1,1,1,1,2,3]
    kmode =  [4,3,1,5,4,3,1,5,5,5]
    kconds = [1,1,1,1,2,2,2,2,2,2]
    nj = length(kmode)
    for n in ens
        for j in 1:nj
#               CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE,
#     $                      'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
#     $                      A, LDA, WORK( 2*N+1 ), IINFO )
            A = zeros(T,n,n)
            simrcond = simrconds[kconds[j]]
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = ulpinv
            verbosity[] > 1 &&
                println("latme: n=$n $anorm, $imode, $rcond, $simrcond")
            latme!(A,anorm,imode,rcond,simrcond)
            schurtest(A,tols[itype])
        end
    end
end

@testset "diagonal, random ev" begin
    itype=7
    kmode = [6]
    kmagn = [1]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            A = zeros(T,n,n)
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = 1.0
            latmr!(A,anorm,imode,rcond,kl=0,ku=0)
            schurtest(A,tols[itype])
        end
    end
end

@testset "symmetric, random ev" begin
    itype=8
    kmagn = [1]
    kmode = [6]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     # note sym='H'
     #           CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            A = zeros(T,n,n)
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = 1.0
            latmr!(A,anorm,imode,rcond,sym='H',kl=n,ku=n)
            schurtest(A,tols[itype])
        end
    end
end

@testset "general, random ev" begin
    itype=9
    kmagn = [1,2,3]
    # kmode = [4,3,1] # data statement is a lie
    kmode = [6,6,6]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            # zdrves zeros first 2 rows/cols & last row w/o explanation
            # presumably to test balancing, which we don't do, so skip it.
            A = zeros(T,n,n)
            imode = kmode[j]
            anorm = magns[kmagn[j]]
            rcond = 1.0
            (verbosity[] > 1) &&
                println("latmr: n=$n $anorm, $imode, $rcond")
            latmr!(A,anorm,imode,rcond)
            schurtest(A,tols[itype])
        end
    end
end

@testset "triangular, specified ev" begin
    # This is done peculiarly to exercise branches in the eigenvector solver.
    # In particular, we stuff the UT directly into a Schur object so we know exactly
    # what's on the diagonal.
    # (Probably happens for normal schur(), but I didn't promise that anywhere.)
    kconds =   [2,2,2,2,2,2,2,2,2,2,2,2]
    kmagn =    [1,1,1,1,2,2,2,2,3,3,3,3]
    kmode =    [5,0,0,0,5,0,0,0,5,0,0,0]
    mydmodes = [0,0,1,2,0,0,1,2,0,0,1,2]
    nj = length(kmagn)
    itype = 6
    for n in ens
        for j in 1:nj
            A = zeros(T,n,n)
            imode = kmode[j]
            rcond = ulpinv
            anorm = magns[kmagn[j]]
            d_mode = mydmodes[j]
            (verbosity[] > 1) &&
                println("latmr: n=$n $anorm, $imode, $rcond")
            if imode == 0
                d = rand(T,n)
                if d_mode == 0
                    d[2] = d[n-1]
                elseif d_mode == 1
                    d[2] = d[n-1] * nextfloat(one(real(T)))
                else
                    d[2] = d[n-1] * prevfloat(one(real(T)))
                end
                latmr!(A, anorm, imode, rcond, kl=0, d=d)
            else
                latmr!(A, anorm, imode, rcond, kl=0)
            end
            S = LinearAlgebra.Schur(A, Diagonal(ones(T,n)), diag(A))
            vectest(A, S, tols[itype])
        end
    end
end

@testset "triangular, random ev" begin
    # trivial, not even run in zdrves
    itype=10
    kconds = [1]
    kmagn = [1]
    kmode = [1]
    nj = length(kmagn)
    for n in ens
        for j in 1:nj
     # note: kl=0
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
#
            A = zeros(T,n,n)
            imode = 6
            anorm = magns[kmagn[j]]
            rcond = 1.0
            latmr!(A,anorm,imode,rcond,kl=0)
            schurtest(A,tols[itype])
        end
    end
end

end # group testset


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
using LaTeXStrings
theme(:wong)
gr(aspect_ratio=1,legend=nothing,axis=nothing,border=nothing,ticks=nothing,
   markerstrokewidth=0,markersize=1)

Q_GLegendre_rq_N3 = [-0.861136311594054,-0.339981043584857,0.339981043584856,0.861136311594053,-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.861136311594054,-0.339981043584857,0.339981043584856,0.861136311594053,0.168314227951316,-0.542461986333870,0.168314227951316,-0.625852241617446,-0.542461986333870,-0.625852241617446];
Q_GLegendre_sq_N3 = [-1.000000000000000,-1.000000000000000,-1.000000000000000,-1.000000000000000,-0.861136311594054,-0.339981043584857,0.339981043584856,0.861136311594053,0.861136311594054,0.339981043584857,-0.339981043584856,-0.861136311594053,-0.542461986333870,0.168314227951316,-0.625852241617446,0.168314227951316,-0.625852241617446,-0.542461986333870];

r1D, w1D = gauss_lobatto_quad(0,0,4)

plot([-1;-1],[-1;1],seriescolor=:black,lw=3,label = nothing)
scatter!([-0.542461986333870],[0.168314227951316],markersize=100,size=(500,500),markerstrokestyles=:dot,markerstrokecolor=:grey,markerstrokewidth=2,markercolor=:white)
plot!([-1;1],[1;-1],seriescolor=:black,lw=3,label = nothing)
plot!([1;-1],[-1;-1],seriescolor=:black,lw=3,label = nothing)
#plot!([r1D[2];r1D[4]],[0;0],seriescolor=:grey,linestyle=:dash,lw=3,label = nothing)
#plot!([0;0],[r1D[2];r1D[4]],seriescolor=:grey,linestyle=:dash,lw=3,label = nothing)
scatter!(Q_GLegendre_rq_N3, Q_GLegendre_sq_N3,markersize=5,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!([-0.542461986333870;-0.625852241617446;-0.339981043584857], [0.168314227951316;0.168314227951316;0.339981043584857],markersize=5,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
plot!([-0.54246198633387,-0.625852241617446],[0.168314227951316;0.168314227951316],seriescolor=:grey,linestyle=:dash,lw=3,label = nothing)
plot!([-0.54246198633387,-0.339981043584857],[0.168314227951316;0.339981043584857],seriescolor=:grey,linestyle=:dash,lw=3,label = nothing)
# scatter!(r1D,ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(r1D,r1D[2]*ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(r1D,zeros(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(r1D,r1D[4]*ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(r1D,-1*ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!([0;r1D[2];0;r1D[4];0],[r1D[2];0;0;0;r1D[4]],markersize=7,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
#scatter!(r1D[2:4],ones(3),markersize=7,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
savefig("SBP_quad.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using RecipesBase
using TypedTables: Table

@recipe function f(table::Table, quantity::Symbol)
    xguide --> "t"
    yguide --> String(quantity)

    for i in eachindex(eval(:($table.$quantity[1])))
        @series begin
            legend --> :false
            table.t, eval(:([value[$i] for value in $table.$quantity]))
        end
    end
end


text:
using SnpArrays, Test, Statistics, Random, TraitSimulation
using DataFrames, LinearAlgebra, OrdinalMultinomialModels, GLM, DelimitedFiles
Random.seed!(1234)

# demo genetic data can be found in Snparrays package directory
filepath = SnpArrays.datadir()
filename = "EUR_subset"
EUR = SnpArray(SnpArrays.datadir(filename * ".bed"));

EUR_data = SnpData(SnpArrays.datadir(filename));
bimfile = EUR_data.snp_info # store the snp_info with the snp names
snpid  = bimfile[:, :snpid] # store the snp names in the snpid vector
causal_snp_index = findfirst(x -> x ==  "rs62057050" , snpid) # find the index of the snp of interest by snpid
locus = convert(Vector{Float64}, @view EUR[:, causal_snp_index])
famfile = EUR_data.person_info
sex = map(x -> strip(x) == "F" ? 0.0 : 1.0, famfile[!, :sex])
X = [sex locus]
β_full =
 [0.52
0.05]

# first ordinal power for the single snp
θ = [1.0, 1.2, 1.4]
trait_ordered = OrderedMultinomialTrait(X, β_full, θ, LogitLink())

n_sim = 10

writepheno(filename, filepath, trait_ordered, n_sim)

pheno_ordered = readdlm("pheno.txt", Any)

# test to see if the outputted phenotype file has the correct number of columns, the first two columns are family ID and Individual ID respectively. The remaining columns are the # of simulated univariate traits
@test size(pheno_ordered, 2) == n_sim + 2
# first column FID
@test famfile[:, 1] == string.(pheno_ordered[2:end, 1])
# second column IID
@test famfile[:, 2] == string.(pheno_ordered[2:end, 2])

## VCM trait bivariate test
β = [β_full β_full]
GRM = grm(EUR, minmaf = 0.05);
I_n = Matrix(I, size(GRM))
ΣE = [1.0 0.0; 0.0 1.0]
ΣA = [0.3 0.0; 0.0 0.3]
vc = @vc ΣA ⊗ 2GRM + ΣE ⊗ I_n;
trait_vcm = VCMTrait(X, β, vc)

n_sim = 5

writepheno(filename, filepath, trait_vcm, n_sim)

pheno_ordered = readdlm("pheno.txt", Any)

# test to see if the outputted phenotype file has the correct number of columns, the first two columns are family ID and Individual ID respectively. The remaining columns are the # of simulated univariate traits
n_traits = ntraits(trait_vcm)

@test size(pheno_ordered, 2) == n_traits * n_sim + 2
# first column FID
@test famfile[:, 1] == string.(pheno_ordered[2:end, 1])
# second column IID
@test famfile[:, 2] == string.(pheno_ordered[2:end, 2])


text:
@defcomp n2oemissions begin
    region = Index()

    e_globalN2Oemissions = Variable(index=[time], unit="Mtonne/year")
    e0_baselineN2Oemissions = Parameter(index=[region], unit="Mtonne/year")
    e_regionalN2Oemissions = Variable(index=[time,region], unit="Mtonne/year")
    er_N2Oemissionsgrowth = Parameter(index=[time,region], unit="%")

    function run_timestep(p, v, d, t)
        # note that Hope (2009) states that Equations 1-12 for methane also apply to N2O

        # eq.4 in Hope (2006) - regional N2O emissions as % change from baseline
        for r in d.region
            v.e_regionalN2Oemissions[t,r] = p.er_N2Oemissionsgrowth[t,r] * p.e0_baselineN2Oemissions[r] / 100
        end

        # eq. 5 in Hope (2006) - global N2O emissions are sum of regional emissions
        v.e_globalN2Oemissions[t] = sum(v.e_regionalN2Oemissions[t,:])
    end
end


text:
# CREATE
"""
    create_entity!(new_object::IEntity,
                   dbconn::LibPQ.Connection;
                   creator::Union{IEntity,Missing} = missing)

Persists an instance to the database.
Properties of type vector of IEntity are not persisted
(see function `update_vector_property!` for that usage).
The instance gets enriched by the values generated by the database in the process
(eg. autogenerated primary keys).
If the optional argument 'creator' is given and that the ORM module of the class
is configured to track the creator then a entry will be inserted in the
modification table
"""
function create_entity!(new_object::IEntity,
                       dbconn::LibPQ.Connection;
                       creator::Union{IEntity,Missing} = missing)

   Controller.create_entity!(new_object,
                             dbconn;
                             creator = creator)
end

"""
    create_in_bulk_using_copy(entities::Vector{T},
                                       dbconn::LibPQ.Connection) where T <: IEntity

Persists many instances to the database using PostgreSQL 'COPY'
"""
function create_in_bulk_using_copy(entities::Vector{T},
                                   dbconn::LibPQ.Connection) where T <: IEntity

   Controller.create_in_bulk_using_copy(entities,
                                        dbconn)
end

# RETRIEVE
"""
    retrieve_entity(filter_object::IEntity,
                    retrieve_complex_props::Bool,
                    dbconn::LibPQ.Connection)

Retrieves a vector or T instances that match the filter.
If there is no match, it returns an empty vector.
If `retrieve_complex_props` is true then a complex prop (i.e. a property of
type IEntity) is fully loaded using an additional query to the database ;
if not, an instance of IEntity is simply created using the values of the foreign
key, the other properties are left to missing.
"""
function retrieve_entity(filter_object::IEntity,
                         retrieve_complex_props::Bool,
                         dbconn::LibPQ.Connection)

     Controller.retrieve_entity(filter_object,
                               retrieve_complex_props,
                               dbconn)
end

"""
    retrieve_one_entity(filter_object::T,
                        retrieve_complex_props::Bool,
                        dbconn::LibPQ.Connection)

Retrieves an instance ot T that matches the filter.
If there is no match then the function returns missing.
If there are more than one match, an error is thrown.
If `retrieve_complex_props` is true then a complex prop (i.e. a property of
type IEntity) is fully loaded using an additional query to the database ;
if not, an instance of IEntity is simply created using the values of the foreign
key, the other properties are left to missing.
"""
function retrieve_one_entity(filter_object::T,
                             retrieve_complex_props::Bool,
                             dbconn::LibPQ.Connection) where T <: IEntity
     Controller.retrieve_one_entity(filter_object,
                                     retrieve_complex_props,
                                     dbconn)
end

"""
    update_entity!(updated_object::IEntity,
                   dbconn::LibPQ.Connection;
                   editor::Union{IEntity,Missing} = missing)

Updates an instance to the database.
Properties of type vector of IEntity are not updated
(see function `update_vector_property!` for that usage).
The instance gets enriched by the values generated by the database in the process
(eg. values generated by triggers).
If the optional argument 'editor' is given and that the ORM module of the class
is configured to track the editor then a entry will be inserted in the
modification table
"""
function update_entity!(updated_object::IEntity,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )
    Controller.update_entity!(updated_object,
                              dbconn;
                              editor = editor
                              )
end

"""
    update_vector_property!(updated_object::IEntity,
                            updated_property::Symbol,
                            dbconn::Union{Missing, LibPQ.Connection};
                            editor::Union{Missing, IEntity} = missing)

Updates a property of type vector of IEntity

"""
function update_vector_property!(updated_object::IEntity,
                                 updated_property::Symbol,
                                 dbconn::Union{Missing, LibPQ.Connection};
                                 editor::Union{Missing, IEntity} = missing)

      Controller.update_vector_property!(updated_object,
                                        updated_property,
                                        dbconn;
                                        editor = editor)

end

"""
    delete_entity(deleted_object::IEntity,
                  dbconn::LibPQ.Connection;
                  editor::Union{IEntity,Missing} = missing)

Deletes an entity on the basis of its ID
"""
function delete_entity(deleted_object::IEntity,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )
    Controller.delete_entity(deleted_object,
                             dbconn;
                             editor = editor
                             )
end

"""
    delete_entity_alike(filter_object::IEntity,
                        dbconn::LibPQ.Connection;
                        editor::Union{IEntity,Missing} = missing)

Deletes all entities matching the the filter
"""
function delete_entity_alike(filter_object::IEntity,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )
    Controller.delete_entity_alike(filter_object,
                                   dbconn;
                                   editor = editor
                                   )
end

"""
   execute_plain_query(query_string::String,
                       query_args::Union{Vector,Missing},
                       dbconn::LibPQ.Connection)

Execute a query and with the given arguments returns a dataframe
"""
function execute_plain_query(query_string::String,
                             query_args::Union{Vector,Missing},
                             dbconn::LibPQ.Connection)
   Controller.execute_plain_query(query_string,
                                  query_args,
                                  dbconn)
end

"""
   execute_query_and_handle_result(query_string::String,
                                   data_type::DataType,
                                   query_args::Union{Vector,Missing},
                                   retrieve_complex_props::Bool,
                                   dbconn::LibPQ.Connection)

Execute a query and with the given arguments and convert the rows to the given
data_type in the same way that `retrieve_entity` would do it.
"""
function execute_query_and_handle_result(query_string::String,
                                         data_type::DataType,
                                         query_args::Union{Vector,Missing},
                                         retrieve_complex_props::Bool,
                                         dbconn::LibPQ.Connection)

    Controller.execute_query_and_handle_result(query_string,
                                               data_type,
                                               query_args,
                                               retrieve_complex_props,
                                               dbconn)
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestWiringRDF
using Test

using Serd, Serd.RDF
using Catlab, Catlab.Theories, Catlab.WiringDiagrams

using SemanticFlowGraphs

const R = Resource

A, B, C, D = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C, :D)
f = Hom(:f, A, B)
g = Hom(:g, B, C)
h = Hom(:h, D, D)
diagram = to_wiring_diagram(otimes(compose(f,g), h))
_, stmts = wiring_diagram_to_rdf(diagram)
#write_rdf(stdout, stmts)

vs = box_ids(diagram)
vin, vout = input_id(diagram), output_id(diagram)
find_vertex = value -> vs[findfirst(v -> box(diagram, v).value == value, vs)]
fv, gv, hv = find_vertex(:f), find_vertex(:g), find_vertex(:h)

# Check RDF type of boxes and ports.
@test Triple(Blank("n$fv"), R("rdf","type"), R("monocl","Box")) in stmts
@test Triple(Blank("n$fv:in1"), R("rdf","type"), R("monocl","Port")) in stmts
@test Triple(Blank("n$fv:out1"), R("rdf","type"), R("monocl","Port")) in stmts

# Check values of boxes and ports.
@test Triple(Blank("n$fv"), R("monocl","value"), Literal("f")) in stmts
@test Triple(Blank("n$fv:in1"), R("monocl","value"), Literal("A")) in stmts
@test Triple(Blank("n$fv:out1"), R("monocl","value"), Literal("B")) in stmts

# Check wires, in both reified and non-reified form.
i = first([ i for (i, wire) in enumerate(wires(diagram))
            if wire.source.box == fv && wire.target.box == gv ])
@test Triple(Blank("e$i"), R("rdf","type"), R("monocl","Wire")) in stmts
@test Triple(Blank("e$i"), R("monocl","source"), Blank("n$fv:out1")) in stmts
@test Triple(Blank("e$i"), R("monocl","target"), Blank("n$gv:in1")) in stmts
@test Triple(Blank("n$fv:out1"), R("monocl","wire"), Blank("n$gv:in1")) in stmts

end

max_stars_repo_name: IBM/semanticflowgraph

max_issues_repo_name: epatters/semanticflowgraph

max_forks_repo_name: IBM/semanticflowgraph


text:
#InspectDR: Event handling
#-------------------------------------------------------------------------------

#=TODO:
Use Reactive?  Don't want to risk slowing down InspectDR by adding another
module for the simple interaction needed here.
=#

#==Types
===============================================================================#
struct HandlerInfo
	listener::Any
	fn::Function
end
#Raises event:
#fn(listener, source, args...)
#    where source is the signal generator.

#==Functions
===============================================================================#
raiseevent(::Nothing, source, args...) = nothing

function raiseevent(info::HandlerInfo, source, args...)
	eventhandler = info.fn
	eventhandler(info.listener, source, args...)
end

function raiseevent(infolist::Vector{HandlerInfo}, source, args...)
	for info in infolist
		raiseevent(info, source, args...)
	end
end

#Last line


text:
######### CmdStan program example  ###########

using CmdStan, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  bernoullidata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
  inittheta = Dict("theta" => 0.60)

  global stanmodel, sdf
  stanmodel = Stanmodel(name="bernoulli", model=bernoullimodel,
    random=CmdStan.Random(seed=-1),
    num_warmup=1000, printsummary=false);

  global rc, sim
  rc, a3d, cnames = stan(stanmodel, bernoullidata, ProjDir,
    init=inittheta, CmdStanDir=CMDSTAN_HOME)

  if rc == 0
    sdf  = read_summary(stanmodel)
    @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.34 rtol=0.1
  end

end # cd


text:
function SL_wait(ta::Array{Union{Task,Nothing}, 1}, t_interval::Int64)
  t = 0
  while t < 300
    if any([!istaskdone(t) for t in ta])
      sleep(t_interval)
      t += t_interval
    elseif t ≥ 60
      println("      one or more queries incomplete after 60 s; skipping test.")
      for i = 1:4
        ta[i] = nothing
        GC.gc()
      end
      break
    else
      tf1 = fetch(ta[1])
      tf2 = fetch(ta[2])
      tf3 = fetch(ta[3])
      tf4 = fetch(ta[4])
      @test tf1 == tf2 == tf3 == tf4
      break
    end
  end
  return nothing
end

function wait_on_data!(S::GphysData, tmax::Real)
  τ = 0.0
  t = 5.0
  printstyled(string("      (sleep up to ", tmax, " s)\n"), color=:green)
  redirect_stdout(out) do
    show(S)

    # Here we actually wait for data to arrive
    while isempty(S)
      if any(isopen.(S.c)) == false
        break
      end
      if τ > tmax-t
        show(S)
        break
      end
      sleep(t)
      τ += t
    end

    # Close the connection cleanly (write & close are redundant, but
    # write should close it instantly)
    for q = 1:length(S.c)
      if isopen(S.c[q])
        if q == 3
          show(S)
        end
        close(S.c[q])
      end
    end
    sleep(t)
    τ += t
  end

  # Synchronize (the reason we used d0,d1 in our test sessions)
  prune!(S)
  if !isempty(S)
    sync!(S, s="first")
    str = string("time elapsed = ", τ, " s")
    @info(str)
    printstyled("      "*str*"\n", color=:green)
  else
    str = (if τ < tmax
      string("connection closed after ", τ, " s")
    else
      string("no data after ", τ, " s...is server up?")
    end)
    @warn(str)
    printstyled("      "*str*"\n", color=:green)
  end
  return nothing
end


text: global const BIOMART_URL = "http://www.ensembl.org/biomart/martservice"


text:
"""
Load the columns of the JDF as vectors
"""

load_columns(jdf::JDFFile; args...) = load_columns(path(jdf); args...)

function load_columns(indir; cols = Symbol[], verbose = false)
    # starting from DataFrames.jl 0.21 the colnames are strings
    cols = collect(string.(cols))

    if verbose
        println("loading $indir in parallel")
    end

    metadatas = open(joinpath(indir, "metadata.jls")) do io
        deserialize(io)
    end

    if length(cols) == 0
        cols =collect(string.(metadatas.names))
    else
        scmn = setdiff(cols, collect(string.(metadatas.names)))
        if length(scmn) > 0
            throw("columns $(reduce((x,y) -> string(x) * ", " * string(y), scmn)) are not available, please ensure you have spelt them correctly")
        end
    end

    results = Vector{Any}(undef, length(cols))
    names = Vector{String}(undef, length(cols))

    # rate limit channel
    c1 = Channel{Bool}(Threads.nthreads())
    atexit(() -> close(c1))

    i = 1
    for (name, metadata) in zip(metadatas.names, metadatas.metadatas)
        name_str = string(name)
        if name_str in cols
            put!(c1, true)
            results[i] = @spawn begin
                io = BufferedInputStream(open(joinpath(indir, string(name)), "r"))
                new_result = column_loader(metadata.type, io, metadata)
                close(io)
                (name = name_str, task = new_result)
            end
            take!(c1)
            i += 1
        end
    end

    # run the collection of results this serially
    result_vectors = Vector{Any}(undef, length(cols))

    for (i, result) in enumerate(results)
        if verbose
            println("Extracting $(fetch(result).name)")
        end

        new_result = fetch(result).task
        colname = fetch(result).name
        names[i] = colname
        if new_result === nothing
            result_vectors[i] = Vector{Missing}(missing, metadatas.rows)
        else
            result_vectors[i] = new_result
        end
    end

    return names, result_vectors
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
### For debugging @async crahshes
import Base.sync_varname
import Base.@async

macro async(expr)

    tryexpr = quote
        try
            $expr
        catch err
            @warn "error within async" exception=err # line $(__source__.line):
            @show stacktrace(catch_backtrace())
        end
    end

    thunk = esc(:(()->($tryexpr)))

    var = esc(sync_varname)
    quote
        local task = Task($thunk)
        if $(Expr(:isdefined, var))
            push!($var, task)
        end
        schedule(task)
        task
    end
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

hpc = Plasm.mkpol(V,EV)
Plasm.view(hpc)

Plasm.mkpol(V,CV)

Plasm.view(V,CV)

Plasm.view(Lar.cuboid([1,1,1],true))

Plasm.view( (V,FV) )


text:
#=
methods for handling backtest results of strategy objects
=#

mutable struct Results
    trades::Dict{String,TS}
    backtest::Dict{String,TS{Float64}}
    optimization::Matrix{Float64}
    function Results(trades::Dict{String,TS}=Dict{String,TS}(),
                     backtest::Dict{String,TS{Float64}}=Dict{String,TS{Float64}}(),
                     optimization::Matrix{Float64}=Matrix{Float64}(0,0))
        return new(trades, backtest, optimization)
    end
end


text:
export PARSDMM_precompute_distribute_Minkowski

function PARSDMM_precompute_distribute_Minkowski(
                                      TD_OP_c1       ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
                                      TD_OP_c2        ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
                                      TD_OP_sum        ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
                                      set_Prop_c1,
                                      set_Prop_c2,
                                      set_Prop_sum,
                                      comp_grid,
                                      options
                                      ) where {TF<:Real,TI<:Integer}

N = prod(comp_grid.n)
p = length(TD_OP_c1)
q = length(TD_OP_c2)
r = length(TD_OP_sum)

if options.feasibility_only==true
  s=p+q+r
else
  s=p+q+r+1
end

AtA=Vector{SparseMatrixCSC{TF,TI}}(undef,s)
ZERO_MAT = spzeros(TF,N,N)
Id_MAT   = SparseMatrixCSC{TF}(LinearAlgebra.I,N,N)

for i=1:p #first component
  if set_Prop_c1.dense[i]==true
    if set_Prop_c1.AtA_diag[i]==true
      AtA[i]=[ Id_MAT ZERO_MAT ; ZERO_MAT ZERO_MAT ]
      #set_Prop_c1.AtA_offsets[i]=[0] #do not need to add additional offset info
    else
      error("provided a dense non orthogoal transform-domain operator")
    end
  else
    AtA[i] =  [ TD_OP_c1[i]'*TD_OP_c1[i] ZERO_MAT ; ZERO_MAT ZERO_MAT ]
    #set_Prop_c1.AtA_offsets[i] #do not need to add additional offset info
  end
end
for i=1:q #second component
  if set_Prop_c2.dense[i]==true
    if set_Prop_c2.AtA_diag[i]==true
      AtA[p+i]= [ ZERO_MAT ZERO_MAT ; ZERO_MAT Id_MAT ]
      #set_Prop_c1.AtA_offsets[p+i]=[0] #do not need to add additional offset info
    else
      error("provided a dense non orthogoal transform-domain operator")
    end
  else
    AtA[p+i] = [ ZERO_MAT ZERO_MAT ; ZERO_MAT TD_OP_c2[i]'*TD_OP_c2[i] ]
    #do not need to add additional offset info
  end
end
for i=1:r #total image
  if set_Prop_sum.dense[i]==true
    if set_Prop_sum.AtA_diag[i]==true
      AtA[i+p+q]=[ Id_MAT Id_MAT ; Id_MAT Id_MAT ]
      #set_Prop_sum.AtA_offsets[i]=[-N 0 N]
    else
      error("provided a dense non orthogoal transform-domain operator")
    end
  else
    AtA[i+p+q] = [ TD_OP_sum[i]'*TD_OP_sum[i] TD_OP_sum[i]'*TD_OP_sum[i] ; TD_OP_sum[i]'*TD_OP_sum[i] TD_OP_sum[i]'*TD_OP_sum[i] ]
    #set_Prop_sum.AtA_offsets[i]=[set_Prop_sum.AtA_offsets[i]-N set_Prop_sum.AtA_offsets[i] set_Prop_sum.AtA_offsets[i]+N]
  end
end

#modify the transform-domain operators to include two blocks per operator, i.e.,
# TD_OP = [A] -> TD_OP = [A 0] or [0 A] or [A A]. This is one row of \tilde{A}
for i=1:length(TD_OP_c1)
  TD_OP_c1[i] = [ TD_OP_c1[i] spzeros(TF,size(TD_OP_c1[i],1),N) ];
end
for i=1:length(TD_OP_c2)
  TD_OP_c2[i] = [ spzeros(TF,size(TD_OP_c2[i],1),N) TD_OP_c2[i] ];
end
for i=1:length(TD_OP_sum)
  TD_OP_sum[i] = [ TD_OP_sum[i] TD_OP_sum[i] ];
end

#add the identity matrix as the operator for the squared distance from the point we want to project
if options.feasibility_only==false
  push!(TD_OP_sum,[ Id_MAT Id_MAT ]);
  push!(set_Prop_sum.TD_n,comp_grid.n)
  push!(set_Prop_sum.AtA_offsets,[0]) #push a dummy value (doesn't matter because the offsets will bet detected automatically)
  push!(set_Prop_sum.banded,true)
  push!(set_Prop_sum.AtA_diag,false)
  push!(set_Prop_sum.dense,false)
  push!(set_Prop_sum.ncvx,false)
  push!(set_Prop_sum.tag,("distance squared","identity","matrix",""))
  AtA[s] = [ Id_MAT Id_MAT ; Id_MAT Id_MAT ]
end

#create one set_Prop for all operators
set_Prop = deepcopy(set_Prop_c1)

append!(set_Prop.AtA_diag,set_Prop_c2.AtA_diag)
append!(set_Prop.AtA_offsets,set_Prop_c2.AtA_offsets)
append!(set_Prop.TD_n,set_Prop_c2.TD_n)
append!(set_Prop.banded,set_Prop_c2.banded)
append!(set_Prop.dense,set_Prop_c2.dense)
append!(set_Prop.ncvx,set_Prop_c2.ncvx)
append!(set_Prop.tag,set_Prop_c2.tag)

append!(set_Prop.AtA_diag,set_Prop_sum.AtA_diag)
append!(set_Prop.AtA_offsets,set_Prop_sum.AtA_offsets)
append!(set_Prop.TD_n,set_Prop_sum.TD_n)
append!(set_Prop.banded,set_Prop_sum.banded)
append!(set_Prop.dense,set_Prop_sum.dense)
append!(set_Prop.ncvx,set_Prop_sum.ncvx)
append!(set_Prop.tag,set_Prop_sum.tag)

#if all AtA are banded -> convert to compressed diagonal storage (CDS/DIA) format
if sum(set_Prop.banded[1:s].=true)==s
  for i=1:s
    (AtA[i],set_Prop.AtA_offsets[i]) = mat2CDS(AtA[i])
    set_Prop.AtA_offsets[i]=convert(Vector{TI},set_Prop.AtA_offsets[i])
  end
   AtA=convert(Vector{Array{TF,2}},AtA);
   set_Prop.AtA_offsets=set_Prop.AtA_offsets[1:s]
end

#allocate arrays of vectors
y       = Vector{Vector{TF}}(undef,s);
l       = Vector{Vector{TF}}(undef,s);

#create one TD_OP that contains all transform-domain operators
TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef,s)
TD_OP[1:p]     = TD_OP_c1
TD_OP[1+p:p+q] = TD_OP_c2
TD_OP[1+p+q:s] = TD_OP_sum

for ii=1:s #initialize all rho's, gamma's, y's and l's
    y[ii]       = zeros(TF,size(TD_OP[ii],1))
    l[ii]       = zeros(TF,length(y[ii]));
end

if options.parallel==true
  TD_OP = distribute(TD_OP)
  y     = distribute(y)
  l     = distribute(l)
end

return TD_OP,set_Prop,AtA,l,y
end


text:
# Packages. Note: make sure that all the packages that you use here are
# (1) available in either your global environment or your current page environment
# (2) are installed in the `.github/workflows/deploy.yml` file. See both the
# Project.toml and the `deploy.yml` file here as examples.
#
using DelimitedFiles

# ========================================================================

###########
### 001 ###
###########

# approach 1; uses DelimitedFiles
function hfun_members_table(params::Vector{String})::String
    path_to_csv = params[1]
    members = readdlm(path_to_csv, ',', skipstart=1)
    # write a simple table
    io = IOBuffer()
    write(io, "
")
    write(io, "
Name	GitHub alias
")
    for (name, alias) in eachrow(members)
        write(io, "
")
        write(io, "$name	")
        write(io, """$alias	""")
        write(io, "
")
    end
    write(io, "
")
    return String(take!(io))
end

###########
### 007 ###
###########

# case 1
hfun_case_1() =
    """
var read from foo is $(pagevar("foo", "var"))
"""

# case 2, note the `@delay`
@delay function hfun_case_2()
    all_tags = globvar("fd_page_tags")
    isnothing(all_tags) && return ""
    all_tags = union(values(all_tags)...)
    tagstr = strip(prod("$t " for t in all_tags))
    return """
tags: { $tagstr }
"""
end

###########
### 008 ###
###########

function lx_capa(com, _)
    # this first line extracts the content of the brace
    content = Franklin.content(com.braces[1])
    output = replace(content, "a" => "A")
    return "**$output**"
end

function env_cap(com, _)
    option = Franklin.content(com.braces[1])
    content = Franklin.content(com)
    output = replace(content, option => uppercase(option))
    return "~~~~~~$output~~~~~~"
end


text:
import MQTT: topic_wildcard_len_check, filter_wildcard_len_check, MQTTException

@testset "topic_wildcard_len_check" begin
    @test_throws MQTTException topic_wildcard_len_check("+")
    @test topic_wildcard_len_check("foo") == nothing
    @test_throws MQTTException topic_wildcard_len_check("#")
    @test_throws MQTTException topic_wildcard_len_check("")
end;

@testset "filter_wildcard_len_check" begin
    @test_throws MQTTException filter_wildcard_len_check("")
    @test_throws MQTTException filter_wildcard_len_check("#/")
    @test_throws MQTTException filter_wildcard_len_check("f+oo/bar/more")
    @test_throws MQTTException filter_wildcard_len_check("f#oo/bar/more")
    @test filter_wildcard_len_check("foo/bar/more") == nothing
    @test filter_wildcard_len_check("foo/bar/more/#") == nothing
    @test filter_wildcard_len_check("foo/+/bar/more") == nothing
end;

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
module GoogleAnalytics

# for interaction with API endpoints
using JSON
using JSON3
using HTTP
import MbedTLS
using Base64
#using Sockets
# core
import Base: show
using Dates
using Lazy
# for data processing
using JuliaDB
# for data storage
using CSVFiles

export authenticate,
       JWTAuth

include("utils/request_data_types.jl")
include("utils/response_data_types.jl")
include("utils/authentication.jl")
include("utils/api_comms.jl")
include("utils/data_formatting.jl")

end # module


text:
using ADCME
using JLD2

@load "data.jld2" z

n = 64
noise = parse(Float64, ARGS[1])

τ = Variable(1.0)
x = placeholder(Float64, shape = [n,])
y = placeholder(Float64, shape = [n,])
w = placeholder(Float64, shape = [n,])
σ, κ = 0.08, 0.5
α = 0.5
Δt = 0.01
yhat = (
    x + κ * (τ - α*x)*Δt + σ * sqrt(x) * sqrt(Δt) * w  + 1/4*σ^2*Δt*(w^2-1)
)/(1+(1-α)*κ*Δt)

θ = Variable(fc_init([2,20,20,20,1]))
hDξ = sigmoid(squeeze(fc([x yhat], [20,20,20,1], θ)))
Dξ = sigmoid(squeeze(fc([x y], [20,20,20,1], θ)))
# LF = mean(log((1-hDξ)/hDξ))
LF = -mean(log(hDξ))
LD = -mean(log(Dξ) + log(1-hDξ))

function generate_batch_data()
    idx = rand(1:9999, n)
    abs.(z[idx] + noise * randn(n)), abs.(z[idx.+1] + noise * randn(n)), randn(n)
end

optF = RMSPropOptimizer(1e-3).minimize(LF, var_list=[τ])
optD = RMSPropOptimizer(1e-3).minimize(LD, var_list=[θ])

sess = Session(); init(sess)
x0, y0, w0 = generate_batch_data()
    fd = Dict(
        x => x0, y => y0, w => w0
    )
@info run(sess, [LF, LD], feed_dict = fd)

fiter = 1
diter = 5
l1 = []
l2 = []
τs = []

for i = 1:50000
    x0, y0, w0 = generate_batch_data()

    fd = Dict(
        x => x0, y => y0, w => w0
    )

    ld = 0
    for k = 1:diter
        _, ld = run(sess, [optD, LD], feed_dict = fd)
    end

    lf = 0
    for k = 1:fiter
        _, lf = run(sess, [optF, LF], feed_dict = fd)
    end

    τ0 = run(sess, τ)

    # if length(loss)==0
    #     loss = [l1[end]]
    # end
    push!(l1, ld)
    push!(l2, lf)
    push!(τs, τ0)
    if mod(i, 100)==0
        @info i, ld, lf, τ0
    end
end

kappa = join(string.(τs), ',')
db = Database("noise.db")
execute(db, """
CREATE TABLE IF NOT EXISTS tau (
    noise real,
    kappa text
)
""")
execute(db, """
INSERT INTO tau VALUES ($noise, \"$kappa\")
""")
close(db)
# close("all")
# plot(τs)
# savefig("test$noise.png")


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains definitions for VertexNodes in the query dag.

################################################# IMPORT/EXPORT ############################################################

###
# A Part of the DAG. The execution of a node returns
# a value which may be used higher up in the dag
###
abstract Node

###
# Denotes a Graph
###
abstract GraphNode <: Node

###
# Denotes a vector of values
###
abstract VectorNode <: Node

###
# Denotes a vertex or edge property
###
abstract Property <: Node

################################################# GRAPHNODE ################################################################

###
# Leaf node containing a symbol denoting a graph object.
###
immutable SimpleGraphNode <: GraphNode
   gs::Symbol
end

Base.show(io::IO, x::SimpleGraphNode) = write(io, "Graph($(x.gs))")

###
# Node denoting a filter on vertices, determined by an array of booleans.
###
immutable FilterNode <: GraphNode
   graph::GraphNode
   bools::VectorNode
end

Base.show(io::IO, x::FilterNode) = write(io, "filter($(x.graph), $(x.bools))")

###
# Node denoting a subset of either vertex properties, or edge properties, or subsets of both.
###
immutable SelectNode <: GraphNode
   graph::GraphNode
   props::Vector{Property}
end

Base.show(io::IO, x::SelectNode) = write(io, "select($(x.props),$(x.graph))")

################################################# VECTORNODE ###############################################################

###
# Denotes a single literal value to be broadcast.
###
immutable LiteralNode{T} <: VectorNode
   val::T
end

Base.show{T}(io::IO, x::LiteralNode{T}) = write(io, "$(x.val)")

###
# Denotes a column from the vertex or edge dataframes
###
immutable TableNode <: VectorNode
   graph::GraphNode
   prop::Property
end

Base.show(io::IO, x::TableNode) = write(io, "Property($(x.graph), $(x.prop))")

###
# Denotes an operation to be performed on two vectors
###
immutable VectorOperation <: VectorNode
   op::Function
   args::Vector{VectorNode}
end

Base.show(io::IO, x::VectorOperation) = write(io, "{$(x.op)($(x.args))}")

################################################# PROPERTY #################################################################

###
# A vertex property
###
immutable VertexProperty <: Property
   prop::Symbol
end

Base.show(io::IO, x::VertexProperty) = write(io, "v.$(x.prop)")

###
# An edge property
###
immutable EdgeProperty <: Property
   prop::Symbol
end

Base.show(io::IO, x::EdgeProperty) = write(io, "e.$(x.prop)")

###
# An edge source vertex property
###
immutable EdgeSourceProperty <: Property
   prop::Symbol
end

Base.show(io::IO, x::EdgeSourceProperty) = write(io, "s.$(x.prop)")

###
# An edge target vertex property
###
immutable EdgeTargetProperty <: Property
   prop::Symbol
end

Base.show(io::IO, x::EdgeTargetProperty) = write(io, "t.$(x.prop)")


text:
@testset "Arc-Swift" begin

    # See paper:
    # Qi & Manning 2017: "Arc-swift: A Novel Transition System for Dependency Parsing"
    # https://nlp.stanford.edu/pubs/qi2017arcswift.pdf

    # note that the transition sequences tested here are compared with
    # the arc-swift implementation at https://github.com/qipeng/arc-swift

    # @test projective_only(ArcSwift())

    chopsticks_tree = test_sentence("chopsticks.conll")
    ketchup_tree = test_sentence("ketchup.conll")

    @testset "Arc Eager Reduce" begin
        oracle = Oracle(ArcEager(), static_oracle, typed)
        @test last.(oracle(chopsticks_tree)) ==
            [Shift(), LeftArc("nsubj"), RightArc("root"), RightArc("dobj"), Reduce(),
             Shift(), LeftArc("case"), RightArc("inst"), Reduce(), RightArc(".")]

        @test last.(oracle(ketchup_tree)) ==
            [Shift(), LeftArc("nsubj"), RightArc("root"), RightArc("dobj"),
             RightArc("nmod"), RightArc("case"), Reduce(), Reduce(), Reduce(),
             RightArc(".")]
    end

    @testset "Arc Eager Shift" begin
        oracle = Oracle(ArcEager(), static_oracle_prefer_shift, typed)

        @test last.(oracle(chopsticks_tree)) ==
            [Shift(), LeftArc("nsubj"), RightArc("root"), RightArc("dobj"),
             Shift(), LeftArc("case"), Reduce(), RightArc("inst"), Reduce(),
             RightArc(".")]

        # @test last.(oracle(ketchup_tree)) ==
        @test last.(oracle(ketchup_tree)) ==
            [Shift(), LeftArc("nsubj"), RightArc("root"),
             RightArc("dobj"), RightArc("nmod"), RightArc("case"),
             Reduce(), Reduce(), Reduce(), RightArc(".")]
    end

    @testset "Arc Swift" begin

        oracle = Oracle(ArcSwift(), static_oracle, typed)

        @test last.(oracle(chopsticks_tree)) ==
            [Shift(), LeftArc(1, "nsubj"), RightArc(1, "root"), RightArc(1, "dobj"),
             Shift(), LeftArc(1, "case"), RightArc(2, "inst"), RightArc(2, ".")]

        @test last.(oracle(ketchup_tree)) ==
            [Shift(), LeftArc(1, "nsubj"), RightArc(1, "root"), RightArc(1, "dobj"),
             RightArc(1, "nmod"), RightArc(1, "case"), RightArc(4, ".")]

        for tree in (chopsticks_tree, ketchup_tree)
            c1 = initconfig(ArcSwift(), [t.form for t in tree])
            c2 = initconfig(ArcSwift(), tree)
            @test stack(c1) == stack(c2)
            @test buffer(c1) == buffer(c2)
            @test [t.form for t in c1.A] == [t.form for t in c2.A]
            @test startswith(showstr(c1), "ArcSwiftConfig")
        end
    end

end


text:
using Brownies
using Test

@testset "Brownies.jl" begin
    include("typestest.jl")
    include("structuretests.jl")
    include("thermodynamictests.jl")
    include("dynamicstests.jl")
    include("utilstests.jl")
end


text:
using SymmetricPolynomials
using Test
import SymmetricPolynomials: semi_elementary_polynomial, dim, is_elementary, to_string
import SymmetricPolynomials: push!, to_polynomial

@testset "Monomials" begin
    include("monomialTests.jl")
end

@testset "Polynomials" begin
    include("polynomialTests.jl")
end

@testset "Decompose" begin
    include("decomposeTests.jl")
end

@testset "Evaluate" begin
    include("evaluateTests.jl")
end

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
function normal_wrench(b::TypedElasticBodyBodyCache{T}) where {T}
    wrench_lin = zeros(SVector{3,T})
    wrench_ang = zeros(SVector{3,T})
    @inbounds begin
    for k_trac = 1:length(b.TractionCache)
        trac = b.TractionCache[k_trac]
        p_dA = calc_p_dA(trac)
        λ_s = p_dA * trac.n̂
        wrench_lin += λ_s
        wrench_ang += cross(trac.r_cart, λ_s)
    end
    end
    return Wrench(b.mesh_2.FrameID, wrench_ang, wrench_lin)
end

function normal_wrench_cop(b::TypedElasticBodyBodyCache{T}) where {T}
    wrench_lin = zeros(SVector{3,T})
    wrench_ang = zeros(SVector{3,T})
    int_p_dA_cop = zeros(SVector{3,T})
    int_p_dA = zero(T)
    @inbounds begin
    for k_trac = 1:length(b.TractionCache)
        trac = b.TractionCache[k_trac]
        p_dA = calc_p_dA(trac)
        λ_s = p_dA * trac.n̂
        wrench_lin += λ_s
        wrench_ang += cross(trac.r_cart, λ_s)
        int_p_dA += p_dA
        int_p_dA_cop += p_dA * trac.r_cart
    end
    end
    return Point3D(b.mesh_2.FrameID, int_p_dA_cop / int_p_dA), Wrench(b.mesh_2.FrameID, wrench_ang, wrench_lin)
end


text:
############### Display Function

import Base.Multimedia.display;
import Base.Multimedia.print;

function print(p::Union{AbstractMonteCarloProcess, AbstractPayoff})
    fldnames = collect(fieldnames(typeof(p)))
    print(typeof(p), "(")
    if (length(fldnames) > 0)
        print(fldnames[1], "=", getfield(p, fldnames[1]))
        popfirst!(fldnames)
        for name in fldnames
            print(",", name, "=", getfield(p, name))
        end
    end
    print(")")
end

function display(p::Union{AbstractMonteCarloProcess, AbstractPayoff})
    fldnames = collect(fieldnames(typeof(p)))
    print(typeof(p), "(")
    if (length(fldnames) > 0)
        print(fldnames[1], "=", getfield(p, fldnames[1]))
        popfirst!(fldnames)
        for name in fldnames
            print(",", name, "=", getfield(p, name))
        end
    end
    println(")")
end

# Agnostic getters and setters, useful for testing and calibration

function get_parameters(model::XProcess) where {XProcess <: AbstractPayoff}
    fields_ = fieldnames(XProcess)
    N1 = length(fields_)
    param_ = []
    for i = 1:N1
        tmp_ = getproperty(model, fields_[i])
        #typeof(tmp_) <: Bool ? continue : nothing;
        append!(param_, tmp_)
    end
    if (N1 > 0)
        typecheck_ = typeof(sum(param_) + prod(param_))
        param_ = convert(Array{typecheck_}, param_)
    end

    return param_
end

function get_parameters(model::XProcess) where {XProcess <: BaseProcess}
    fields_ = fieldnames(XProcess)
    N1 = length(fields_)
    param_ = []
    for i = 1:N1
        tmp_ = getproperty(model, fields_[i])
        typeof(tmp_) <: AbstractUnderlying ? continue : nothing
        append!(param_, tmp_)
    end
    typecheck_ = typeof(sum(param_) + prod(param_))
    param_ = convert(Array{typecheck_}, param_)

    return param_
end

function get_parameters(model::XProcess) where {XProcess <: VectorialMonteCarloProcess}
    fields_ = fieldnames(XProcess)
    N1 = length(fields_)
    param_ = []
    for i = 1:N1
        tmp_ = getproperty(model, fields_[i])
        typeof(tmp_) <: Tuple ? continue : nothing
        append!(param_, tmp_)
    end
    typecheck_ = typeof(sum(param_) + prod(param_))
    param_ = convert(Array{typecheck_}, param_)

    return param_
end

function set_parameters!(model::XProcess, param::Array{num}) where {XProcess <: BaseProcess, num <: Number}
    fields_ = fieldnames(XProcess)
    N1 = length(fields_) - 1
    if N1 != length(param)
        error("Check number of parameters of the model")
    end

    for (symb, nval) in zip(fields_, param)
        if (symb != :underlying)
            setproperty!(model, symb, nval)
        end
    end
end

function set_parameters!(model::LogNormalMixture, param::Array{num}) where {num <: Number}
    N1 = div(length(param) + 1, 2)
    if N1 * 2 != length(param) + 1 || N1 < 2
        error("Check number of parameters of the model")
    end
    eta = param[1:N1]
    lam = param[(N1+1):end]
    fields_ = fieldnames(LogNormalMixture)
    setproperty!(model, fields_[1], eta)
    setproperty!(model, fields_[2], lam)
    return
end

function set_parameters!(model::ShiftedLogNormalMixture, param::Array{num}) where {num <: Number}
    N1 = div(length(param), 2)
    if N1 * 2 != length(param) || N1 < 2
        error("Check number of parameters of the model")
    end
    eta = param[1:N1]
    lam = param[(N1+1):(2*N1-1)]
    alfa = param[end]
    fields_ = fieldnames(ShiftedLogNormalMixture)
    setproperty!(model, fields_[1], eta)
    setproperty!(model, fields_[2], lam)
    setproperty!(model, fields_[3], alfa)
    return
end