text:
#apply_weights(A::Array, W::Array) = sum(A .* W)
apply_weights(A::Array, W::Missing) = sum(A .* (1/length(A)))

@inline function orthogonal_sine(t::Real, j::Real, block_num::Integer, block_width::Real)
    # Note for simplicity and speed I am going to assume that we will never call this after (k+1)h_n. This is due to the way it will be calculated in the loop.
    start_of_block = block_num*block_width
    return  t > start_of_block ? ((sqrt(2 * block_width))/(j*pi)) *sin((j * pi * (t-start_of_block))/block_width) : 0.0
end

function spectral_lmm_array(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts); regularisation::Union{Missing,Symbol} = :covariance_default, regularisation_params::Dict = Dict(),
                          only_regulise_if_not_PSD::Bool = false, numJ::Integer = 100, num_blocks::Integer = 10, block_width::Real = (maximum(ts.df[:,ts.time])) / num_blocks,
                          microstructure_noise_var::Dict{Symbol,<:Real} = two_scales_volatility(ts, assets)[2])
    numAssets = length(assets)
    SS = zeros(num_blocks, numJ, numAssets)
    indices = Dict(assets .=> collect(1:numAssets))
    lasttimes = zeros(numAssets)
    lastvalues = zeros(numAssets)
    H_kn = zeros(num_blocks, numAssets)
    for i in 1:nrow(ts.df)
        asset = ts.df[i, ts.grouping]
        if (asset in assets) == false continue end
        assetslot = indices[asset]
        newtime = ts.df[i, ts.time]
        newval  = ts.df[i, ts.value]
        if lasttimes[assetslot] < eps()
            lasttimes[assetslot]  = newtime
            lastvalues[assetslot] = newval
            continue
        end
        timegap = (newtime - lasttimes[assetslot])
        centertime = lasttimes[assetslot] + 0.5*timegap
        ret = simple_differencing([newval], [lastvalues[assetslot]])[1]
        lasttimes[assetslot]  = newtime
        lastvalues[assetslot] = newval

        block_num = min(num_blocks-1,Int(floor(newtime/block_width)))

        # Looking over j values for the S Array.
        for j in 1:numJ
             SS[block_num+1, j, assetslot] += (pi * j/block_width) * ret * orthogonal_sine(centertime, j, block_num, block_width)
        end
        # Now doing the H_k^n matrix which is independent of j.
        H_kn[block_num+1, assetslot] += (timegap^2) * (microstructure_noise_var[asset]/block_width)
    end
    # Now we assemble the covariance matrices for every j,k.
    R = eltype(SS)
    uncorrected_covar_matrices = Array{Array{Hermitian{R},1},1}(undef,num_blocks)
    for block_num in 1:num_blocks
        hermitian_array = Array{Hermitian{R},1}(undef, numJ)
        for j in 1:numJ
            hermitian_array[j] = Hermitian( SS[block_num,j,:] * transpose(SS[block_num,j,:])  )
        end
        uncorrected_covar_matrices[block_num] = hermitian_array
    end

    # The efficient weights (see page 1324 of Bibinger, Hautsch, Malec & Reiss 2014) are not implemented here. It is highly error proone as
    # many inversions are required which is 1) computationally intensive and 2) can cause singularity errors as there is no guarantee of psd.
    # So going to fall back to equal weighting.
    weights = repeat([missing],num_blocks)

    corrected_matrices = Array{CovarianceMatrix{R},1}(undef,num_blocks)
    for block_num in 1:num_blocks
        weights_for_block = weights[block_num]
        block_H = Diagonal(H_kn[block_num,:])
        adjustments = ((pi .*  collect(1:numJ)) ./block_width).^2
        block_adjustments = Array{Diagonal,1}(undef,numJ)
        for nn in 1:numJ block_adjustments[nn] = adjustments[nn] * block_H end
        corrected_mats = uncorrected_covar_matrices[block_num] .- block_adjustments
        mat = Hermitian(apply_weights(corrected_mats, weights_for_block))

        # Regularisation
        dont_regulise = ismissing(regularisation) || (only_regulise_if_not_PSD && is_psd_matrix(mat))
        mat = dont_regulise ? mat : regularise(mat, ts, assets, regularisation; regularisation_params... )

        # In some cases we get negative terms on the diagonal with this algorithm.
        negative_diagonals = findall(diag(mat) .< eps())
        covar = make_nan_covariance_matrix(assets, ts.time_period_per_unit)
        if length(negative_diagonals) < 1
            corr, vols = cov2cor_and_vol(mat, 1) # No adjustment for duration as we already have the spot matrix.
            covar.correlation = corr
            covar.volatility = vols
         end
        corrected_matrices[block_num] = covar
    end
    return corrected_matrices
end

"""
    spectral_covariance(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                        regularisation::Union{Missing,Symbol} = :covariance_default,
                        regularisation_params::Dict = Dict(), only_regulise_if_not_PSD::Bool = false,
                        numJ::Integer = 100, num_blocks::Integer = 10,
                        block_width::Real = (maximum(ts.df[:,ts.time]) - minimum(ts.df[:,ts.time])) / num_blocks,
                        microstructure_noise_var::Dict{Symbol,<:Real} = two_scales_volatility(ts, assets)[2])

Estimation of a CovarianceMatrix using the spectral covariance method.
### Inputs
* `ts` - The tick data.
* `assets` - The assets you want to estimate volatilities for.
* `regularisation` - A symbol representing what regularisation technique should be used. If missing no regularisation is performed.
* `regularisation_params` - keyword arguments to be consumed in the regularisation algorithm.
* `only_regulise_if_not_PSD` - Should regularisation only be attempted if the matrix is not psd already.
* `numJ` - The number of J values. See the paper for details.
* `num_blocks` - The number of blocks to split the time frame into. See the paper for details.
* `block_width` - The width of each block to split the time frame into.
* `microstructure_noise_var` - Estimates of microstructure noise variance for each asset.
### Returns
* A `CovarianceMatrix`.

### References
Bibinger M, Hautsch N, Malec P, Reiss M (2014). “Estimating the quadratic covariation matrix from noisy observations: Local method of moments and efficiency.” The Annals of Statistics, 42(4), 1312–1346. doi:10.1214/14-AOS1224.
"""
function spectral_covariance(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts);
                             regularisation::Union{Missing,Symbol} = :covariance_default,
                             regularisation_params::Dict = Dict(), only_regulise_if_not_PSD::Bool = false,
                             numJ::Integer = 100, num_blocks::Integer = 10,
                             block_width::Real = (maximum(ts.df[:,ts.time]) - minimum(ts.df[:,ts.time])) / num_blocks,
                             microstructure_noise_var::Dict{Symbol,<:Real} = two_scales_volatility(ts, assets)[2])
    corrected_matrices = spectral_lmm_array(ts, assets; regularisation = regularisation, regularisation_params = regularisation_params, only_regulise_if_not_PSD = only_regulise_if_not_PSD, numJ = numJ, num_blocks = num_blocks,
                                            block_width = block_width, microstructure_noise_var = microstructure_noise_var)
    cor_weights = repeat([1.0], length(corrected_matrices))
    spectral_estimate = combine_covariance_matrices(corrected_matrices, cor_weights)
    return spectral_estimate
end


text:
include("./tower.jl")

function sample_solver(n)
    if n == 5
        const top = [2 3 2 1 2]
        const right = [2 4 3 4 1]
        const bottom = [2 2 3 4 1]
        const left = [3 1 3 2 2]
        solve_tower(n, top, right, bottom, left)
    elseif n == 3
        const top = [2 2 1]
        const right = [1 2 2]
        const bottom = [2 1 3]
        const left = [3 1 2]
        solve_tower(n, top, right, bottom, left)
    end
end

ext: jl

max_stars_repo_name: roscroft/solver

max_issues_repo_name: roscroft/solver

max_forks_repo_name: roscroft/solver


text:
using IncrementalInference
using Test

@testset "Test delete tree clique" begin

fg = generateCanonicalFG_lineStep(3;
                                   poseEvery=1,
                                   landmarkEvery=3,
                                   posePriorsAt=[],
                                   landmarkPriorsAt=[0],
                                   sightDistance=2,
                                   solverParams=SolverParams(algorithms=[:default, :parametric]))

# getSolverParams(fg).graphinit = false
# getSolverParams(fg).treeinit = true
# getSolverParams(fg).dbg = true
# getSolverParams(fg).useMsgLikelihoods = true
# getSolverParams(fg).dbg = true
# getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = true

smtasks = Task[]
oldtree, smt, hists = solveTree!(fg; smtasks=smtasks, verbose=false, recordcliqs=ls(fg));

@test IIF.isRoot(oldtree, IIF.CliqueId(1))
@test IIF.isRoot(oldtree, IIF.getClique(oldtree,1))
@test !IIF.isRoot(oldtree, IIF.CliqueId(2))
@test !IIF.isRoot(oldtree, IIF.CliqueId(3))

IIF.deleteClique!(oldtree, IIF.CliqueId(1))

@test IIF.isRoot(oldtree, IIF.CliqueId(2))
@test IIF.isRoot(oldtree, IIF.CliqueId(3))

@test IIF.getClique(oldtree, :x0) == IIF.getClique(oldtree, IIF.CliqueId(3)) == IIF.getClique(oldtree, 1)
@test IIF.getClique(oldtree, :x3) == IIF.getClique(oldtree, IIF.CliqueId(2)) == IIF.getClique(oldtree, 2)

# drawTree(oldtree, show=true)

tree, smt, hists = solveTree!(fg, oldtree; smtasks=smtasks, verbose=false, recordcliqs=ls(fg));

# csmAnimate(tree, hists, frames=1)

end


text:
export combineLatest
export PushEach, PushEachBut, PushNew, PushNewBut, PushStrategy

import Base: show

"""
    combineLatest(sources...; strategy = PushEach())
    combineLatest(sources::S, strategy::G = PushEach()) where { S <: Tuple, U }

Combines multiple Observables to create an Observable whose values are calculated from the latest values of each of its input Observables.
Accept optinal update strategy object.

# Arguments
- `sources`: input sources
- `strategy`: optional update strategy for batching new values together

# Examples
```jldoctest
using Rocket

latest = combineLatest(of(1), from(2:5))

subscribe!(latest, logger())
;

# output

[LogActor] Data: (1, 2)
[LogActor] Data: (1, 3)
[LogActor] Data: (1, 4)
[LogActor] Data: (1, 5)
[LogActor] Completed
```

```
using Rocket

latest = combineLatest(of(1) |> async(0), from(2:5) |> async(0), strategy = PushNew())

subscribe!(latest, logger())
;

# output

[LogActor] Data: (1, 2)
[LogActor] Completed
```

See also: [`Subscribable`](@ref), [`subscribe!`](@ref), [`PushEach`](@ref), [`PushEachBut`](@ref), [`PushNew`](@ref), [`PushNewBut`](@ref), [`PushStrategy`](@ref)
"""
function combineLatest end

"""
    PushEach

`PushEach` update strategy specifies combineLatest operator to emit new value each time an inner observable emit a new value

See also: [`combineLatest`](@ref), [`PushEachBut`](@ref), [`PushNew`](@ref), [`PushNewBut`](@ref), [`PushStrategy`](@ref)
"""
struct PushEach end

"""
    PushEachBut{I}

`PushEachBut` update strategy specifies combineLatest operator to emit new value if and only if an inner observable with index `I` have a new value

See also: [`combineLatest`](@ref), [`PushEach`](@ref), [`PushNew`](@ref), [`PushNewBut`](@ref), [`PushStrategy`](@ref)
"""
struct PushEachBut{I} end

"""
    PushNew

`PushNew` update strategy specifies combineLatest operator to emit new value if and only if all inner observables have a new value

See also: [`combineLatest`](@ref), [`PushEach`](@ref), [`PushEachBut`](@ref), [`PushNewBut`](@ref), [`PushStrategy`](@ref)
"""
struct PushNew end

"""
    PushNewBut{I}

`PushNewBut{I}` update strategy specifies combineLatest operator to emit new value if and only if all inner observables except with index `I` have a new value

See also: [`combineLatest`](@ref), [`PushEach`](@ref), [`PushEachBut`](@ref), [`PushNew`](@ref), [`PushStrategy`](@ref)
"""
struct PushNewBut{I} end

"""
    PushStrategy(strategy::BitArray{1})

`PushStrategy` update strategy specifies combineLatest operator to emit new value if and only if all inner observables with index such that `strategy[index] = false` have a new value

See also: [`combineLatest`](@ref), [`PushEach`](@ref), [`PushEachBut`](@ref), [`PushNew`](@ref), [`PushNewBut`](@ref), [`collectLatest`](@ref)
"""
struct PushStrategy
    strategy :: BitArray{1}
end

function push_update!(::Int, ::BitArray{1}, ::BitArray{1}, ::PushEach)
    return nothing
end

function push_update!(::Int, vstatus::BitArray{1}, ::BitArray{1}, ::PushEachBut{I}) where I
    @inbounds vstatus[I] = false
    return nothing
end

function push_update!(nsize::Int, vstatus::BitArray{1}, cstatus::BitArray, ::PushNew)
    unsafe_copyto!(vstatus, 1, cstatus, 1, nsize)
    return nothing
end

function push_update!(nsize::Int, vstatus::BitArray{1}, cstatus::BitArray, ::PushNewBut{I}) where I
    push_update!(nsize, vstatus, cstatus, PushNew())
    @inbounds vstatus[I] = true
    return nothing
end

function push_update!(nsize::Int, vstatus::BitArray{1}, cstatus::BitArray, strategy::PushStrategy)
    push_update!(nsize, vstatus, cstatus, PushNew())
    map!(|, vstatus, vstatus, strategy.strategy)
    return nothing
end

combineLatest(; strategy = PushEach())                                       = error("combineLatest operator expects at least one inner observable on input")
combineLatest(args...; strategy = PushEach())                                = combineLatest(args, strategy)
combineLatest(sources::S, strategy::G = PushEach()) where { S <: Tuple, G }  = CombineLatestObservable{combined_type(sources), S, G}(sources, strategy)
combineLatest(sources::V, strategy::G = PushEach()) where { V <: Vector, G } = CombineLatestObservable{combined_type(sources), V, G}(sources, strategy)

##

struct CombineLatestInnerActor{L, W, I} <: Actor{L}
    wrapper :: W
end

Base.show(io::IO, inner::CombineLatestInnerActor{L, W, I}) where { L, W, I } = print(io, "CombineLatestInnerActor($L, $I)")

on_next!(actor::CombineLatestInnerActor{L, W, I}, data::L) where { L, W, I } = next_received!(actor.wrapper, data, Val{I}())
on_error!(actor::CombineLatestInnerActor{L, W, I}, err)    where { L, W, I } = error_received!(actor.wrapper, err, Val{I}())
on_complete!(actor::CombineLatestInnerActor{L, W, I})      where { L, W, I } = complete_received!(actor.wrapper, Val{I}())

##

struct CombineLatestActorWrapper{S, A, G}
    storage :: S
    actor   :: A

    nsize    :: Int
    strategy :: G           # Push update strategy
    cstatus  :: BitArray{1} # Completion status
    vstatus  :: BitArray{1} # Values status
    ustatus  :: BitArray{1} # Updates status

    subscriptions :: Vector{Teardown}

    CombineLatestActorWrapper{S, A, G}(storage::S, actor::A, strategy::G) where { S, A, G } = begin
        nsize   = length(storage)
        cstatus = falses(nsize)
        vstatus = falses(nsize)
        ustatus = falses(nsize)
        subscriptions = fill!(Vector{Teardown}(undef, nsize), voidTeardown)
        return new(storage, actor, nsize, strategy, cstatus, vstatus, ustatus, subscriptions)
    end
end

push_update!(wrapper::CombineLatestActorWrapper) = push_update!(wrapper.nsize, wrapper.vstatus, wrapper.cstatus, wrapper.strategy)

cstatus(wrapper::CombineLatestActorWrapper, index) = @inbounds wrapper.cstatus[index]
vstatus(wrapper::CombineLatestActorWrapper, index) = @inbounds wrapper.vstatus[index]
ustatus(wrapper::CombineLatestActorWrapper, index) = @inbounds wrapper.ustatus[index]

dispose(wrapper::CombineLatestActorWrapper) = begin fill!(wrapper.cstatus, true); foreach(s -> unsubscribe!(s), wrapper.subscriptions) end

function next_received!(wrapper::CombineLatestActorWrapper, data, index::Val{I}) where I
    setstorage!(wrapper.storage, data, index)
    @inbounds wrapper.vstatus[I] = true
    @inbounds wrapper.ustatus[I] = true
    if all(wrapper.vstatus) && !all(wrapper.cstatus)
        push_update!(wrapper)
        next!(wrapper.actor, snapshot(wrapper.storage))
    end
end

function error_received!(wrapper::CombineLatestActorWrapper, err, index::Val{I}) where I
    if !(@inbounds wrapper.cstatus[I])
        dispose(wrapper)
        error!(wrapper.actor, err)
    end
end

function complete_received!(wrapper::CombineLatestActorWrapper, ::Val{I}) where I
    if !all(wrapper.cstatus)
        @inbounds wrapper.cstatus[I] = true
        if ustatus(wrapper, I)
            @inbounds wrapper.vstatus[I] = true
        end
        if all(wrapper.cstatus) || (@inbounds wrapper.vstatus[I] === false)
            dispose(wrapper)
            complete!(wrapper.actor)
        end
    end
end

##

struct CombineLatestObservable{T, S, G} <: Subscribable{T}
    sources  :: S
    strategy :: G
end

function on_subscribe!(observable::CombineLatestObservable{T, S, G}, actor::A) where { T, S, G, A }
    storage = getmstorage(T)
    wrapper = CombineLatestActorWrapper{typeof(storage), A, G}(storage, actor, observable.strategy)

    for (index, source) in enumerate(observable.sources)
        @inbounds wrapper.subscriptions[index] = subscribe!(source, CombineLatestInnerActor{eltype(source), typeof(wrapper), index}(wrapper))
        if cstatus(wrapper, index) === true && vstatus(wrapper, index) === false
            dispose(wrapper)
            break
        end
    end

    if all(wrapper.cstatus)
        dispose(wrapper)
    end

    return CombineLatestSubscription(wrapper)
end

##

struct CombineLatestSubscription{W} <: Teardown
    wrapper :: W
end

as_teardown(::Type{ <: CombineLatestSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::CombineLatestSubscription)
    dispose(subscription.wrapper)
    return nothing
end

Base.show(io::IO, ::CombineLatestObservable{D}) where D  = print(io, "CombineLatestObservable($D)")
Base.show(io::IO, ::CombineLatestSubscription)           = print(io, "CombineLatestSubscription()")

##


text:
@info "Starting PointCloud tests..."
@testset "PointCloud" begin
    for (T) in [
            (Float32),
            (Float64),
            (Int64),
        ]
        points = rand(T,1024,3)
        rep = PointCloud(points)

        @test rep.npoints == 1024
        @test size(rep.points) == (1024,3)
        @test rep.points isa Array{Float32,2}
    end
end


text:
using Libdl

depsfile = joinpath(dirname(@__FILE__),"deps.jl")

if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    f = open(depsfile,"w")
    if Sys.iswindows()
        path = replace(path, "\\" => "\\\\")
    end
    println(f,"const libxprs = \"$(path)\"")
    close(f)
end

libname = string(Sys.iswindows() ? "" : "lib", "xprs", ".", Libdl.dlext)
paths_to_try = String[]

push!(paths_to_try, libname)

if haskey(ENV, "XPRESSDIR")
  push!(paths_to_try, joinpath(ENV["XPRESSDIR"], Sys.iswindows() ? "bin" : "lib", libname))
end

global found_xprs = false
for l in paths_to_try
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        global found_xprs = true
        @info("Found $l")
        write_depsfile(l)
        break
    end
end

if !found_xprs
    error("Unable to locate Xpress installation, please check your enviroment variable XPRESSDIR . Note that Xpress must be obtained separately from fico.com")
end


text:
# we order the time dimension last because Julia is column-major (see:
# https://github.com/Alexander-Barth/NCDatasets.jl/issues/87#issuecomment-636098859)

using NCDatasets
using OrderedCollections

mutable struct NetCDFWriter <: AbstractWriter
    filename::Union{Nothing, String}

    NetCDFWriter() = new(nothing)
end

function full_name(writer::NetCDFWriter, filename = nothing)
    if !isnothing(filename)
        return filename * ".nc"
    else
        return writer.filename * ".nc"
    end
end

function init_data(nc::NetCDFWriter, filename, dims, vars)
    Dataset(full_name(nc, filename), "c") do ds
        # define spatial and time dimensions
        for (dn, (dv, da)) in dims
            defDim(ds, dn, length(dv))
        end
        defDim(ds, "time", Inf) # Inf sets UNLIMITED dimension

        # include dimensions as variables
        for (dn, (dv, da)) in dims
            defVar(ds, dn, dv, (dn,), attrib = da)
        end
        defVar(
            ds,
            "time",
            Float64,
            ("time",),
            attrib = OrderedDict(
                "units" => "seconds since 1900-01-01 00:00:00",
                "long_name" => "time",
            ),
        )

        # define variables
        for (vn, (vd, vt, va)) in vars
            defVar(ds, vn, vt, (vd..., "time"), attrib = va)
        end
    end
    nc.filename = filename
    return nothing
end

function append_data(nc::NetCDFWriter, varvals, simtime)
    Dataset(full_name(nc), "a") do ds
        timevar = ds["time"]
        t = length(timevar) + 1
        timevar[t] = simtime
        for (vn, vv) in varvals
            dsvar = ds[vn]
            dsvar[ntuple(_ -> Colon(), ndims(vv))..., t] = vv
        end
    end
    return nothing
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
abstract type AbstractSigma end

get(sigma::AbstractSigma) = sigma.σ

@with_kw struct ConstSigma <: AbstractSigma
    σ::Real
end

ConstSigma() = ConstSigma(1)

get(sigma::ConstSigma) = [sigma.σ]

@with_kw struct VarSigma <: AbstractSigma
    σ::Array{Float64,1}
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

"""
     parse_scene(
          inputstream::InputStream,
          variables::Dict{String, Float64} = Dict()
          ) :: Scene

Read a scene description from the given `inputstream` and
return a `Scene` object.
Throws `GrammarError` if an error occurs.

Call internally the following parsing functions:
- [`read_token`](@ref)
- [`StopToken`](@ref)
- [`expect_identifier`](@ref)
- [`expect_symbol`](@ref)
- [`expect_number`](@ref)
- [`parse_sphere`](@ref)
- [`parse_plane`](@ref)
- [`parse_camera`](@ref)
- [`parse_material`](@ref)

Call internally the following functions and structs of the program:
- [`add_shape!`](@ref)
- [`add_light!`](@ref)

See also: [`InputStream`](@ref), [`Scene`](@ref)
"""
function parse_scene(inputstream::InputStream, variables::Dict{String, Float64} = Dict{String, Float64}())
#function parse_scene(inputstream::InputStream, variables::Dict{T1, T2} = Dict{T1, T2}()) where {T1, T2}
     scene = Scene()
     scene.float_variables = copy(variables)
     scene.overridden_variables = keys(variables)

     while true
          what = read_token(inputstream)
          isa(what.value, StopToken) && (break)

          if !isa(what.value, KeywordToken)
               throw(GrammarError(what.location, "expected a keyword instead of '$(what)'"))
          end

          if what.value.keyword == FLOAT
               variable_name = expect_identifier(inputstream)

               # Save this for the error message
               variable_loc = inputstream.location

               expect_symbol(inputstream, "(")
               variable_value = expect_number(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    # Only define the variable if it was not defined by the user *outside* the scene file
                    # (e.g., from the command line)
                    scene.float_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == STRING
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = expect_string(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.string_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == BOOL
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = expect_bool(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.bool_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == VECTOR
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = parse_vector(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.vector_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == COLOR
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = parse_color(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.color_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == BRDFS
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = parse_brdf(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.brdf_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == PIGMENT
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = parse_pigment(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.pigment_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == TRANSFORMATION
               variable_name = expect_identifier(inputstream)
               variable_loc = inputstream.location
               expect_symbol(inputstream, "(")
               variable_value = parse_transformation(inputstream, scene)
               expect_symbol(inputstream, ")")

               if (variable_name ∈ scene.variable_names) && !(variable_name ∈ scene.overridden_variables)
                    throw(GrammarError(variable_loc, "variable «$(variable_name)» cannot be redefined"))
               end

               if variable_name ∉ scene.overridden_variables
                    scene.transformation_variables[variable_name] = variable_value
                    push!(scene.variable_names, variable_name)
               end

          elseif what.value.keyword == POINTLIGHT
               add_light!(scene.world, parse_pointlight(inputstream, scene))

          elseif what.value.keyword == CAMERA
               if !isnothing(scene.camera)
                    throw(GrammarError(what.location, "You cannot define more than one camera"))
               end
               scene.camera = parse_camera(inputstream, scene)

          elseif what.value.keyword ==  MATERIAL
               name, material = parse_material(inputstream, scene)
               scene.materials[name] = material

          elseif what.value.keyword == PRINT
               println(inputstream, scene)

          elseif what.value.keyword == ASSERT
               assert(inputstream, scene)

          elseif what.value.keyword == SPHERE
               add_shape!(scene.world, parse_sphere(inputstream, scene))

          elseif what.value.keyword == PLANE
               add_shape!(scene.world, parse_plane(inputstream, scene))

          elseif what.value.keyword == CUBE
               add_shape!(scene.world, parse_cube(inputstream, scene))

          elseif what.value.keyword == TRIANGLE
               add_shape!(scene.world, parse_triangle(inputstream, scene))
          end
     end

     return scene
end


text:
let

        # A �~F~R C �~F~P B
        bn = BayesNet()
        push!(bn, StaticCPD(:a, Categorical([1.0,0.0])))
        push!(bn, StaticCPD(:b, Categorical([0.0,1.0])))
        push!(bn, CategoricalCPD{Bernoulli}(:c, [:a, :b], [2,2], [Bernoulli(0.1), Bernoulli(0.2), Bernoulli(1.0), Bernoulli(0.4)]))

        t5 = gibbs_sample(bn, 5, 100; thinning=5, consistent_with=Assignment(),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test size(t5) == (5, 3)
        @test t5[!,:a] == [1,1,1,1,1]
        @test t5[!,:b] == [2,2,2,2,2]
        @test t5[!,:c] == [1,1,1,1,1]

        config = GibbsSampler(burn_in=100, thinning=5)
        t5 = rand(bn, config, 5)
        @test size(t5) == (5, 3)
        @test t5[!,:a] == [1,1,1,1,1]
        @test t5[!,:b] == [2,2,2,2,2]
        @test t5[!,:c] == [1,1,1,1,1]

        bn2 = BayesNet()
        push!(bn2, StaticCPD(:a, Categorical([0.5,0.5])))
        push!(bn2, StaticCPD(:b, Categorical([0.5,0.5])))
        push!(bn2, CategoricalCPD{Categorical}(:c, [:a, :b], [2,2], [
               Categorical([1.0, 0, 0, 0]), Categorical([0, 1.0, 0, 0]),
               Categorical([0, 0, 1.0, 0]), Categorical([0, 0, 0, 1.0])]))

        t6 = gibbs_sample(bn2, 5, 100; thinning=5, consistent_with=Assignment(:c=>1),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test t6[!,:a] == [1,1,1,1,1]
        @test t6[!,:b] == [1,1,1,1,1]
        @test t6[!,:c] == [1,1,1,1,1]

	config.evidence = Assignment(:c=>1)
	t6 = rand(bn2, config, 5)
	@test t6[!,:a] == [1,1,1,1,1]
        @test t6[!,:b] == [1,1,1,1,1]
        @test t6[!,:c] == [1,1,1,1,1]

        t7 = gibbs_sample(bn2, 5, 100; thinning=5, consistent_with=Assignment(:c=>2),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test t7[!,:a] == [2,2,2,2,2]
        @test t7[!,:b] == [1,1,1,1,1]
        @test t7[!,:c] == [2,2,2,2,2]

        rand()
	config.evidence = Assignment(:c=>2)
        t7 = rand(bn2, config, 5)
        @test t7[!,:a] == [2,2,2,2,2]
        @test t7[!,:b] == [1,1,1,1,1]
        @test t7[!,:c] == [2,2,2,2,2]

        bn3 = BayesNet()
        push!(bn3, StaticCPD(:a, Normal(2.5, 1.5)))
        push!(bn3, StaticCPD(:b, Categorical([0.3,0.3,0.4])))
        push!(bn3, LinearGaussianCPD(:c, [:a, :b], [0.5, 0.5], 1.0, 0.5))

        t8 = gibbs_sample(bn3, 5, 100; thinning=5, consistent_with=Assignment(),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test size(t8) == (5, 3)

        # unlikely c
        t9 = gibbs_sample(bn3, 5, 100; thinning=5, consistent_with=Assignment(:c=>1.0),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test size(t9) == (5, 3)

        # use optional parameters and border cases for other parameters
        t10 = gibbs_sample(bn3, 5, 0; thinning=0, consistent_with=Assignment(:c=>2.0, :b=>1),
             variable_order=(Vector{Symbol}([:c, :a, :b])), time_limit=1000000,
             error_if_time_out=false, initial_sample=Assignment(:a=>1, :b=>1, :c=>2.0))
        @test size(t10)[2] == 3 && size(t10)[1] <= 5

	config.burn_in = 0
	config.evidence = Assignment(:c=>2.0, :b=>1)
	config.variable_order = Vector{Symbol}([:c, :a, :b])
	config.time_limit = 1000000
	config.error_if_time_out = false
	config.initial_sample = Assignment(:a=>1, :b=>1, :c=>2.0)
	t10 = rand(bn3, config, 5)
        @test size(t10)[2] == 3 && size(t10)[1] <= 5

        # test early return from time limit
        t11 = gibbs_sample(bn3, 25000, 100; thinning=4, consistent_with=Assignment(),
             variable_order=nothing, time_limit= 1000,
             error_if_time_out=false, initial_sample=nothing)
        @test size(t10)[2] == 3 && size(t10)[1] <= 25000

	config = GibbsSampler(burn_in=100, thinning=4, time_limit = 1000, error_if_time_out=false)
	t11 = rand(bn3, config, 25000)
        @test size(t10)[2] == 3 && size(t10)[1] <= 25000

        d_bn = DiscreteBayesNet()
        push!(d_bn, DiscreteCPD(:a, [1.0,0.0]))
        push!(d_bn, DiscreteCPD(:b, [0.0,1.0]))
        push!(d_bn, CategoricalCPD{Categorical{Float64,Vector{Float64}}}(:c, [:a, :b], [2,2],
                        [Categorical([0.1, 0.9]), Categorical([0.2, 0.8]), Categorical([1.0, 0.0]), Categorical([0.4, 0.6])]))

	t12 = gibbs_sample(d_bn, 5, 100; thinning=5, consistent_with=Assignment(),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test size(t12) == (5, 3)
        @test t12[!,:a] == [1,1,1,1,1]
        @test t12[!,:b] == [2,2,2,2,2]
        @test t12[!,:c] == [1,1,1,1,1]

        t13 = gibbs_sample(d_bn, 5, 100; thinning=5, consistent_with=Assignment(:c=>1),
             variable_order=nothing, time_limit=nothing,
             error_if_time_out=true, initial_sample=nothing)
        @test size(t13) == (5, 3)
        @test t13[!,:a] == [1,1,1,1,1]
        @test t13[!,:b] == [2,2,2,2,2]
        @test t13[!,:c] == [1,1,1,1,1]

        # test bad parameters given to gibbs_sample
        pass = false
        try
            error_test = gibbs_sample(bn, 0, 100; thinning=5, consistent_with=Assignment(),
                 variable_order=nothing, time_limit=nothing,
                 error_if_time_out=true, initial_sample=nothing)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            error_test = gibbs_sample(bn, 5, -1; thinning=5, consistent_with=Assignment(),
                 variable_order=nothing, time_limit=nothing,
                 error_if_time_out=true, initial_sample=nothing)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            error_test = gibbs_sample(bn, 5, 100; thinning=-1, consistent_with=Assignment(),
                 variable_order=nothing, time_limit=nothing,
                 error_if_time_out=true, initial_sample=nothing)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            v_order = Vector{Symbol}([:a, :c])
            error_test = gibbs_sample(bn, 5, 100; thinning=5, consistent_with=Assignment(),
                 variable_order=v_order, time_limit=nothing,
                 error_if_time_out=true, initial_sample=nothing)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            v_order = Vector{Symbol}([:a, :c])
            config = GibbsSampler(burn_in=100, variable_order=v_order, thinning=5)
            error_test = rand(bn, config, 5)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            v_order = Vector{Symbol}([:a, :c, :b, :d])
            error_test = gibbs_sample(bn, 5, 100; thinning=5, consistent_with=Assignment(),
                 variable_order=v_order, time_limit=nothing,
                 error_if_time_out=true, initial_sample=nothing)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            error_test = gibbs_sample(bn, 5, 100; thinning=5, consistent_with=Assignment(),
                 variable_order=nothing, time_limit=0,
                 error_if_time_out=true, initial_sample=nothing)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            error_test = gibbs_sample(bn, 5, 100; thinning=5, consistent_with=Assignment(),
                 variable_order=nothing, time_limit=nothing,
                 error_if_time_out=true, initial_sample=Assignment(:a => 1, :b => 1))
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            config = GibbsSampler(burn_in=100, thinning=5, initial_sample=Assignment(:a => 1, :b => 1))
            error_test = rand(bn, config, 5)
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            error_test = gibbs_sample(bn, 5, 100; thinning=5, consistent_with=Assignment(:c => 2),
                 variable_order=nothing, time_limit=nothing,
                 error_if_time_out=true, initial_sample=Assignment(:a => 1, :b => 1, :c => 1))
        catch e
            pass = true
        end
        @test pass

        pass = false
        try
            error_test = gibbs_sample(bn2, 5, 100; initial_sample=Assignment(:a => 1, :b => 1, :c => 3))
        catch e
            pass = true
        end
        @test pass
end


text:
function _get_time_series(
    container::OptimizationContainer,
    component::PSY.Component,
    attributes::TimeSeriesAttributes{T},
) where {T <: PSY.TimeSeriesData}
    initial_time = get_initial_time(container)
    time_steps = get_time_steps(container)
    forecast = PSY.get_time_series(
        T,
        component,
        get_time_series_name(attributes);
        start_time=initial_time,
        count=1,
    )
    ts_vector = IS.get_time_series_values(
        component,
        forecast,
        initial_time;
        len=length(time_steps),
        ignore_scaling_factors=true,
    )
    return ts_vector
end

function _get_time_series(
    container::OptimizationContainer,
    component::PSY.HybridSystem,
    subcomponent::S,
    attributes::TimeSeriesAttributes{T},
) where {S <: PSY.Component, T <: PSY.TimeSeriesData}
    initial_time = get_initial_time(container)
    time_steps = get_time_steps(container)
    forecast = PSY.get_time_series(
        T,
        component,
        make_subsystem_time_series_name(subcomponent, get_time_series_name(attributes));
        start_time=initial_time,
        count=1,
    )
    ts_vector = IS.get_time_series_values(
        component,
        forecast,
        initial_time;
        len=length(time_steps),
        ignore_scaling_factors=true,
    )
    return ts_vector
end

function get_time_series(
    container::OptimizationContainer,
    component::T,
    parameter::TimeSeriesParameter,
    meta=CONTAINER_KEY_EMPTY_META,
) where {T <: PSY.Component}
    parameter_container = get_parameter(container, parameter, T, meta)
    return _get_time_series(container, component, parameter_container.attributes)
end

function get_time_series(
    container::OptimizationContainer,
    component::S,
    subcomponent_type::Type{T},
    parameter::TimeSeriesParameter,
    meta=CONTAINER_KEY_EMPTY_META,
) where {S <: PSY.HybridSystem, T <: PSY.Component}
    parameter_container = get_parameter(container, parameter, S, meta)
    subcomponent = get_subcomponent(component, subcomponent_type)
    return _get_time_series(
        container,
        component,
        subcomponent,
        parameter_container.attributes,
    )
end

# This is just for temporary compatibility with current code. Needs to be eliminated once the time series
# refactor is done.
function get_time_series(
    container::OptimizationContainer,
    component::PSY.Component,
    forecast_name::String,
)
    return _get_time_series(
        container,
        component,
        TimeSeriesAttributes(PSY.Deterministic, forecast_name),
    )
end


text:
#= -------------------------------------------------------------------------
# @file problema2_PID.jl
#
# @date 12/08/17 01:24:52
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
using Sundials
using PyPlot
#=------------------------------------------------------------------------------
                              constants
------------------------------------------------------------------------------=#
const  A      =  2.0
const  R      =  2.390057361376673e-5
const  cₑ     =  4.184e6
const  V      =  0.5
const  θ_ref  =  60.0 + 273.0
const θ₀_inicial     =  rand(273.0:290.0)
const  G_inv  = 0.4 * 836800.0
const  θₐ    =  293.0
q_nec = ((θ_ref - θₐ) / (R)) - A * G_inv
const  τ      = R * cₑ * V
α = cₑ * V
f(u) = u[2]
# funcion para generar una perturbacion aleatoria escalon
θ_a_f(t) = (100.0 < t < 300.0) ? θₐ - rand(1.0:10.0): θₐ

"""
Funcion que modela el problema 2
"""
function heat_panel(t, θ, dθ)
   dθ[1] = (f(u) / (α)) + ( ((A * G_inv)) / α - ((θ[1] - θₐ) / τ ) + q_nec / α)
end

function heat_panel_pert(t, θ, dθ)
   dθ[1] = (f(u) / (α)) + ( ((A * G_inv)) / α - ((θ[1] - θ_a_f(t)) / τ ) + q_nec / α)
end
#=------------------------------------------------------------------------------
                        constants of the algorithm
------------------------------------------------------------------------------=#
K_p = 1 / R
K_i = 837.0
tauI = 110
tauD = 0.0
t_sim = 500.0
δt  = .01
time_sim = linspace(0, 500, 10000)
set_point = zeros(time_sim) * θ_ref
error = zeros(time_sim) # zeros like t
θ = zeros(time_sim) # zeros like t
#= θ₀ = zeros(time_sim) # zeros like t =#
dθ = zeros(time_sim) # zeros like t
error_int = zeros(time_sim) # zeros like t
P = zeros(time_sim)
I = zeros(time_sim)
D = zeros(time_sim)
pv = zeros(time_sim)
controller_out = zeros(time_sim)
# Upper and Lower limits on OP
controller_max = 3e6
controller_min = 0.0

θ₀ = zeros(1)
θ₀[1] = θ₀_inicial
controller_min = 0
pv[1] = θ_ref
u = ones(time_sim) * θ_ref
for sample in 0:length(time_sim) - 2
   Δt = time_sim[sample+1] - time_sim[sample+2]
   error[sample+1] = set_point[sample+1] - pv[sample+1]
   if sample >= 1  # calculate starting on second cycle
     error_int[sample+1] = error_int[sample] + error[sample+1] * Δt
   end
   P[sample+1] = K_p * error[sample+1]
   I[sample+1] = K_i* error_int[sample+1]
   controller_out[sample+1] = controller_out[1] + P[sample+1] + I[sample+1]
   if controller_out[sample+1] > controller_max  # check upper limit
      controller_out[sample+1] = controller_max
      error_int[sample+1] = error_int[sample+1] - error[sample+1] * Δt # anti-reset windup
   end
   if controller_out[sample+1] < controller_min  # check lower limit
      controller_out[sample+1] = controller_min
      error_int[sample+1] = error_int[sample+1] - error[sample+1] * Δt # anti-reset windup
   end
   ts = [time_sim[sample+1],time_sim[sample+2]]
   #= @show ts =#
   u[sample+2]= controller_out[sample+1]
   sol = Sundials.cvode(heat_panel_pert, θ₀,ts)
   #= sol = integrate.odeint(heat_panel, θ₀[1], ts) =#
   θ[sample+2] = sol[end]
   θ₀[1] = θ[sample+2]
   pv[sample+2] = θ[sample+2]
end

plot(time_sim[2:end], θ[2:end] - 273.0)
xlabel(L"t\,[seg]")
ylabel(L"T\,[^{\circ}C]")
grid("on")
matplotlib["rcParams"][:update](["font.size" => 14, "font.family" => "serif"])
text(0, θ₀_inicial - 273.0, latexstring("\\theta_{0}="*"$(θ₀_inicial-273.0)"), fontsize=10)
legend()


text:
export epgRotation, epgRelaxation, epgDephasing, rfRotation

"""
    epgRotation(alpha::Float64, F::Vector{T}, Z::Vector{T}; statesConsidered=nothing, phi::Float64=0.0)

applies Bloch-rotation (<=> RF pulse) to a set of EPG states.

# Arguments
* `alpha::Float64`           - flip angle of the RF pulse
* `F::Vector{T}`             - transverse dephasing stats
* `Z::Vector{T}`             - longitudinal dephasing stats
* `statesConsidered=nothing` - number of dephasing states to consider (nothing means all states are taken into account)
* `phi::Float64=0.0`         - phase of the RF pulse
"""
function epgRotation(alpha::Float64, F::Vector{T}, Z::Vector{T}; statesConsidered=nothing, phi::Float64=0.0) where T
  # apply rotation to all states per default
  numStates = length(Z)

  n = numStates
  statesConsidered!=nothing && (n=statesConsidered)

  R = rfRotation(alpha, phi)

  rotatedF = fill!(copy(F), 0)
  rotatedZ = fill!(copy(Z), 0)

  for i = 0:n-1
    posIdx = numStates + i
    negIdx = numStates - i

    rotStates = R*[F[posIdx]; conj(F[negIdx]); Z[i+1]]
    rotatedF[posIdx] = rotStates[1]
    rotatedF[negIdx] = conj(rotStates[2])
    rotatedZ[i+1] = rotStates[3]
  end

  return rotatedF, rotatedZ
end

"""
    epgRelaxation( R1::Float64, R2::Float64, t::Float64, F::Vector{T}, Z::Vector{T}) where T

applies relaxation matrices to a set of EPG states.

# Arguments
* `R1::Float64`   - R1
* `R2::Float64`   - R2
* `t::Float64`    - length of time interval in s
* `F::Vector{T}`  - transverse dephasing stats
* `Z::Vector{T}`  - longitudinal dephasing stats
"""
function epgRelaxation( R1::Float64, R2::Float64, t::Float64, F::Vector{T}, Z::Vector{T}) where T

  relaxedF = exp(-R2*t).*F
  relaxedZ = fill!(copy(Z), 0)
  relaxedZ[2:end] = exp(-R1*t)*Z[2:end]
  relaxedZ[1] = exp(-R1*t)*(Z[1]-1.0) + 1.0

  return relaxedF, relaxedZ
end

"""
    epgDephasing(F::Vector{T}, n=1) where T = circshift(F[:],n)

shifts the transverse dephasing states `F` corresponding to n dephasing-cycles.
"""
epgDephasing(F::Vector{T}, n=1) where T = circshift(F[:],n)

"""
    rfRotation(alpha, phi=0.)

returns the rotation matrix for a pulse with flip angle `alpha` and phase `phi`.
"""
function rfRotation(alpha, phi=0.)
  R = [ cos(alpha/2.)^2   exp(2*im*phi)*sin(alpha/2.)^2   -im*exp(im*phi)*sin(alpha);
        exp(-2*im*phi)*sin(alpha/2.)^2   cos(alpha/2.)^2   im*exp(-im*phi)*sin(alpha);
        -im/2 .*exp(-im*phi)*sin(alpha)   im/2 .*exp(im*phi)*sin(alpha)   cos(alpha) ]
end


text:
export BasicDQNLearner

using Random
using Flux

struct BasicDQNLearner{Q,F,R} <: AbstractLearner
    approximator::Q
    loss_func::F
    γ::Float32
    batch_size::Int
    min_replay_history::Int
    rng::R
end

(learner::BasicDQNLearner)(obs) =
    obs |> get_state |>
    x ->
        send_to_device(device(learner.approximator), x) |> learner.approximator |>
        send_to_host

function BasicDQNLearner(;
    approximator::Q,
    loss_func::F,
    γ = 0.99f0,
    batch_size = 32,
    min_replay_history = 32,
    seed = nothing,
) where {Q,F}
    rng = MersenneTwister(seed)
    BasicDQNLearner{Q,F,typeof(rng)}(
        approximator,
        loss_func,
        γ,
        batch_size,
        min_replay_history,
        rng,
    )
end

function RLBase.update!(learner::BasicDQNLearner, batch)
    Q, γ, loss_func, batch_size =
        learner.approximator, learner.γ, learner.loss_func, learner.batch_size
    s, r, t, s′ = map(
        x -> send_to_device(device(Q), x),
        (batch.state, batch.reward, batch.terminal, batch.next_state),
    )
    a = CartesianIndex.(batch.action, 1:batch_size)

    gs = gradient(params(Q)) do
        q = batch_estimate(Q, s)[a]
        q′ = vec(maximum(batch_estimate(Q, s′); dims = 1))
        G = r .+ γ .* (1 .- t) .* q′
        loss_func(G, q)
    end

    update!(Q, gs)
end

function RLBase.extract_experience(
    t::CircularCompactSARTSATrajectory,
    learner::BasicDQNLearner,
)
    if length(t) > learner.min_replay_history
        inds = rand(learner.rng, 1:length(t), learner.batch_size)
        map(get_trace(t, :state, :action, :reward, :terminal, :next_state)) do x
            consecutive_view(x, inds)
        end
    else
        nothing
    end
end


text:
export tri_indices

"""
Iterator over Cartesian indices of the stored region of a `TriMatrix` in an
order corresponding to the elements of its data array. Created with the
[`tri_indices`](@ref) function.
"""
struct TriIndexIterator{L<:TriLayout}
	n::Int

	TriIndexIterator(::L, n::Int) where {L<:TriLayout} = new{L}(n)
end

Base.IteratorSize(::Type{<:TriIndexIterator}) = Base.HasLength()
Base.length(ti::TriIndexIterator{L}) where L = nelems(L, ti.n)
Base.IteratorEltype(::Type{<:TriIndexIterator}) = Base.HasEltype()
Base.eltype(::Type{<:TriIndexIterator}) = Tuple{Int, Int}
TriLayout(::TriIndexIterator{L}) where L = L()

Base.iterate(e::TriIndexIterator) = iterate(e, (1, 1))

function Base.iterate(e::TriIndexIterator{L}, state::Tuple{Int, Int}) where L
	i, j = state

	# Reached end of column
	if j > i
		i += 1
		j = 1
	end

	# Finished last row
	i > (hasdiag(L) ? e.n : e.n - 1) && return nothing

	# Advance column
	next = (i, j + 1)

	value = cartesian_from_tril(L(), i, j)
	return value, next
end

"""
	$(FUNCTIONNAME)(layout::TriLayout, n::Integer)

Iterate over the cartesian indices of a `TriMatrix` corresponding to the linear
indices of its data array, in order.

This should be equivalent to the following generator expression:

    (lin2car(layout, i) for i in 1:nelems(layout, n))
"""
tri_indices(layout::TriLayout, n::Integer) = TriIndexIterator(layout, n)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct BgpSettings <: SwaggerModel
    asn::Any # spec type: Union{ Nothing, Int64 } # spec name: asn
    bgpPeeringAddress::Any # spec type: Union{ Nothing, String } # spec name: bgpPeeringAddress
    peerWeight::Any # spec type: Union{ Nothing, Int32 } # spec name: peerWeight
    bgpPeeringAddresses::Any # spec type: Union{ Nothing, Vector{IPConfigurationBgpPeeringAddress} } # spec name: bgpPeeringAddresses

    function BgpSettings(;asn=nothing, bgpPeeringAddress=nothing, peerWeight=nothing, bgpPeeringAddresses=nothing)
        o = new()
        validate_property(BgpSettings, Symbol("asn"), asn)
        setfield!(o, Symbol("asn"), asn)
        validate_property(BgpSettings, Symbol("bgpPeeringAddress"), bgpPeeringAddress)
        setfield!(o, Symbol("bgpPeeringAddress"), bgpPeeringAddress)
        validate_property(BgpSettings, Symbol("peerWeight"), peerWeight)
        setfield!(o, Symbol("peerWeight"), peerWeight)
        validate_property(BgpSettings, Symbol("bgpPeeringAddresses"), bgpPeeringAddresses)
        setfield!(o, Symbol("bgpPeeringAddresses"), bgpPeeringAddresses)
        o
    end
end # type BgpSettings

const _property_map_BgpSettings = Dict{Symbol,Symbol}(Symbol("asn")=>Symbol("asn"), Symbol("bgpPeeringAddress")=>Symbol("bgpPeeringAddress"), Symbol("peerWeight")=>Symbol("peerWeight"), Symbol("bgpPeeringAddresses")=>Symbol("bgpPeeringAddresses"))
const _property_types_BgpSettings = Dict{Symbol,String}(Symbol("asn")=>"Int64", Symbol("bgpPeeringAddress")=>"String", Symbol("peerWeight")=>"Int32", Symbol("bgpPeeringAddresses")=>"Vector{IPConfigurationBgpPeeringAddress}")
Base.propertynames(::Type{ BgpSettings }) = collect(keys(_property_map_BgpSettings))
Swagger.property_type(::Type{ BgpSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BgpSettings[name]))}
Swagger.field_name(::Type{ BgpSettings }, property_name::Symbol) =  _property_map_BgpSettings[property_name]

function check_required(o::BgpSettings)
    true
end

function validate_property(::Type{ BgpSettings }, name::Symbol, val)
    if name === Symbol("asn")
        Swagger.validate_param(name, "BgpSettings", :maximum, val, 4294967295, false)
        Swagger.validate_param(name, "BgpSettings", :minimum, val, 0, false)
    end
end


text:
"""
    shannon(se::SummarizedExperiment, assay_name::String)

Evaluates the shannon diversity index for each sample and returns it as a vector.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
"""
function shannon(se::SummarizedExperiment, assay_name::String)

    # create empty vector with length = n° samples
    shannon_vector = zeros(size(assay(se, assay_name))[2])

    # evaluate shannon diversity index for each sample / column
    for sample in 1:size(assay(se, assay_name))[2]

        shannon_vector[sample] = shannon(assay(se, assay_name)[:, sample])

    end

    # output vector with results for each sample
    return shannon_vector

end

"""
    shannon!(se::SummarizedExperiment, assay_name::String)

Evaluates the shannon diversity index for each sample and stores it into `coldata(se)`.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
"""
function shannon!(se::SummarizedExperiment, assay_name::String)

    shannon_vector = zeros(size(assay(se, assay_name))[2])

    for sample in 1:size(assay(se, assay_name))[2]

        shannon_vector[sample] = shannon(assay(se, assay_name)[:, sample])

    end

    se.coldata[!, :shannon] = shannon_vector
    return se

end

"""
    ginisimpson(se::SummarizedExperiment, assay_name::String)

Evaluates the ginisimpson diversity index for each sample and returns it as a vector.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
"""
function ginisimpson(se::SummarizedExperiment, assay_name::String)

    ginisimpson_vector = zeros(size(assay(se, assay_name))[2])

    for sample in 1:size(assay(se, assay_name))[2]

        ginisimpson_vector[sample] = ginisimpson(assay(se, assay_name)[:, sample])

    end

    return ginisimpson_vector

end

"""
    ginisimpson!(se::SummarizedExperiment, assay_name::String)

Evaluates the ginisimpson diversity index for each sample and stores it into `coldata(se)`.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
"""
function ginisimpson!(se::SummarizedExperiment, assay_name::String)

    ginisimpson_vector = zeros(size(assay(se, assay_name))[2])

    for sample in 1:size(assay(se, assay_name))[2]

        ginisimpson_vector[sample] = ginisimpson(assay(se, assay_name)[:, sample])

    end

    se.coldata[!, :ginisimpson] = ginisimpson_vector
    return se

end


text:
using Millboard
using Base.Test

board = ["한글"]
@test """
+===+======+
|   |    1 |
+===+======+
| 1 | 한글 |
+---+------+""" == board |> table |> string

board = ["한\n글"; "abc"]
@test """
+===+=====+
|   |   1 |
+===+=====+
|   |  한 |
| 1 |  글 |
+---+-----+
| 2 | abc |
+---+-----+""" == board |> table |> string


text:
@testset "Test LAStopping" begin
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
#
# The following examples illustrate a specialized Stopping for the
# LinearOperators type.
# https://github.com/JuliaSmoothOptimizers/LinearOperators.jl
#
using LinearAlgebra, SparseArrays, LinearOperators #v.1.0.1

##############################################################
#TEST

n = 5
A = rand(n,n)
xref = 100 * rand(n)
x0 = zeros(n)
b    = A * xref
sA  = sparse(A)
opA = LinearOperator(A)

mLO = LinearSystem(A, b)
sLO = LLSModel(sA, b)
opLO = LinearSystem(opA, b)
meta = StoppingMeta(max_cntrs = init_max_counters_linear_operators())
mLOstp_meta = LAStopping(mLO, meta, GenericState(x0)) #this is different because of optimality_check
mLOstp = LAStopping(mLO, GenericState(x0), max_cntrs = init_max_counters_linear_operators())
sLOstp = LAStopping(sLO, GenericState(x0))
short_stop = LAStopping(A,b, sparse = true) #note that sparse is true by default
maxcn = init_max_counters_linear_operators(nprod = 1)
opLOstp = LAStopping(opLO, GenericState(x0), max_cntrs = maxcn)

opLOstp2 = LAStopping(opLO, GenericState(x0), optimality_check = linear_system_check)
@test opLOstp2.meta.max_cntrs == init_max_counters_linear_operators()
@test opLOstp2.meta.optimality_check == opLOstp.meta.optimality_check

mLOstp_src = LAStopping(mLO, meta, StopRemoteControl(), GenericState(x0))

@test typeof(mLOstp.pb.counters) <: LACounters
@test typeof(sLOstp.pb.counters) <: NLSCounters
@test typeof(opLOstp.pb.counters) <: LACounters
@test typeof(short_stop.pb.counters) <: NLSCounters
@test typeof(short_stop.pb) <: LLSModel

@test start!(mLOstp) == false
@test start!(sLOstp) == false
@test start!(opLOstp) == false

update_and_stop!(mLOstp, x = xref, res = Float64[])
@test status(mLOstp) == :Optimal
@test mLOstp.meta.resources == false
update_and_stop!(sLOstp, x = xref, res = Float64[])
@test status(sLOstp) == :Optimal
@test sLOstp.meta.resources == false
update_and_stop!(opLOstp, x = xref, res = Float64[])
@test status(opLOstp) == :Optimal
@test opLOstp.meta.resources == true

include("_additional_linear_operator_algebra.jl")

"""
Randomized block Kaczmarz
"""
function RandomizedBlockKaczmarz(stp :: AbstractStopping; kwargs...)

    b   = stp.pb.b
    x0  = stp.current_state.x
    m, n = length(b), length(x0)
    A   = try
           stp.pb.A #LSSModel.A does not exist
       catch
           sparse(stp.pb.Arows, stp.pb.Acols, stp.pb.Avals, m, n)
       end

    xk  = x0

    OK = start!(stp)

    while !OK

     i  = Int(floor(rand() * m)+1) #rand a number between 1 and m
     Ai = A[i,:]
     xk = Ai == 0 ? x0 : x0 - (dot(Ai,x0)-b[i])/dot(Ai,Ai) * Ai
     #xk  = Ai == 0 ? x0 : x0 - (Ai' * x0-b[i])/(Ai' * Ai) * Ai

     if xtype(stp.current_state) <: SparseVector && !(typeof(xk) <: SparseVector)
        OK = update_and_stop!(stp, x = sparse(xk), res = sparse(Float64[]))
     elseif xtype(stp.current_state) <: Vector && typeof(xk) <: SparseVector
        OK = update_and_stop!(stp, x = convert(Vector, xk), res = Float64[])
     else
        OK = update_and_stop!(stp, x = xk, res = Stopping._init_field(typeof(stp.current_state.res)))
     end
     x0  = xk

    end

 return stp
end

m, n = 200, 100 #size of A: m x n
A    = 100 * rand(m, n) #It's a dense matrix :)
xref = 100 * rand(n)
b    = A * xref

#Our initial guess
x0 = zeros(n)

la_stop = LAStopping(A, b, GenericState(x0),
                     max_iter = 150000,
                     rtol = 1e-6,
                     max_cntrs = init_max_counters_NLS(residual = 150000))
#Be careful using GenericState(x0) would not work here without forcing convert = true
#in the update function. As the iterate will be a SparseVector to the contrary of initial guess.
#Tangi: maybe start! should send a Warning for such problem !?
sa_stop = LAStopping(sparse(A), b, GenericState(sparse(x0)),
                     max_iter = 150000,
                     rtol = 1e-6,
                     max_cntrs = init_max_counters_NLS(residual = 150000))
op_stop = LAStopping(LinearSystem(LinearOperator(A), b),
                     GenericState(x0),
                     max_iter = 150000,
                     rtol = 1e-6,
                     max_cntrs = init_max_counters_linear_operators(nprod = 150000))
opbis_stop = LAStopping(LinearOperator(A), b)

try
 @timed RandomizedBlockKaczmarz(la_stop)
 @test status(la_stop) == :Optimal
 @timed RandomizedBlockKaczmarz(sa_stop)
 @test status(sa_stop) == :Optimal
catch
    @warn "If LSSModel.A does not exist consider [la_stop.pb.Avals[i,j] for (i) in la_stop.pb.Arows, j in la_stop.pb.Acols]"
    #https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/lls_model.jl
end

@timed RandomizedBlockKaczmarz(op_stop)
@test status(op_stop) == :Optimal

update!(la_stop.current_state, x = xref)
@test normal_equation_check(la_stop.pb, la_stop.current_state) <= la_stop.meta.atol
update!(op_stop.current_state, x = xref)
@test normal_equation_check(op_stop.pb, op_stop.current_state) <= la_stop.meta.atol

end


text:
module ThermoState

using Unitful
import Unitful: @u_str

include("utilities.jl")
include("docs.jl")
include("types.jl")
include("spec.jl")
include("state_type.jl")
include("units.jl")
include("show.jl")
include("specs_model.jl")
include("fromspecs_props.jl")

#Spec types
export Spec,ThermodynamicState,VariableSpec

#spec functions
export spec,state

#spec utilities
export value,specification,values,get_spec
#FromState Model
export FromState

#unitful utilities
export normalize_units,default_units,convert_unit

#property functions
export pressure,temperature,mass,moles,molar_mass
export mass_volume, mol_volume, total_volume
export mass_density, mol_density
export mass_enthalpy, mol_enthalpy,total_enthalpy
export mass_gibbs, mol_gibbs, total_gibbs
export mass_helmholtz, mol_helmholtz, total_helmholtz
export mass_internal_energy, mol_internal_energy, total_internal_energy
export mass_entropy, mol_entropy, total_entropy
export mass_fraction, mol_fraction
export mass_number, mol_number
export options, phase, quality, mass_vapor_fraction,mol_vapor_fraction
export molecular_weight

#macros
export @to_units, @spec_str

#state type for dispach
export state_type

#modules
export Types,QuickStates,StatePoints

end # module


text:
@inline _eta(z) = sqrt(1+z*z)+log(z/(1+sqrt(1+z*z)))
@inline _p(z)   = inv(sqrt(1+z*z))

function Uk_polynomials(max_order)
  P0   = Polynomial([1.0])
  out  = [P0]
  mul_int  = Polynomial([1.0, 0.0, -5.0])
  mul_frnt = Polynomial([0.0, 0.0, 1.0, 0.0, -1.0])/2
  for j in 1:max_order
    Pjm1 = out[end]
    Pjm1_int = integrate(mul_int*Pjm1)/8
    Pjm1_drv = derivative(Pjm1)
    newP     = mul_frnt*Pjm1_drv + Pjm1_int - Pjm1_int(0.0)/8
    push!(out, newP)
  end
  out
end

# Should I compute these and hard-code them into this file?
if !(@isdefined UK_COEF_TUPLES)
  const UK_POLYS = Uk_polynomials(20)
end

function _besselk_vz_asv(v, z, maxorder, modify=false)
  ez  = _eta(z)
  pz  = _p(z)
  if modify
    # same trick: (x^v)*exp(-v*...) in the log.
    mulval = exp(v*log(z*v)-v*ez)
  else
    mulval = exp(-v*ez)
  end
  out = sqrt(pi/(2*v))*mulval/sqrt(sqrt(1+z*z))
  # series part:
  ser = zero(z)
  sgn = one(z)
  _v  = one(v)
  @inbounds for j in 1:maxorder
    ser += sgn*UK_POLYS[j](pz)/_v
    sgn *= -one(z)
    _v  *= v
  end
  out*ser
end

_besselk_asv(v, z, maxorder=10, modify=false) = _besselk_vz_asv(v, z/v, maxorder, modify)


text:
using BinaryBuilder

# Collection of sources required to build Pixman
name = "Pixman"
version = v"0.40.0"
sources = [
    ArchiveSource("https://www.cairographics.org/releases/pixman-$(version).tar.gz",
                  "6d200dec3740d9ec4ec8d1180e25779c00bc749f94278c8b9021f5534db223fc"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pixman-*/

./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libpixman-1", :libpixman)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# truncation.jl
#
# Implements truncation schemes for truncating a tensor with svd, leftorth or rightorth

abstract type TruncationScheme end

struct NoTruncation <: TruncationScheme
end
notrunc() = NoTruncation()

struct TruncationError <: TruncationScheme
    epsilon::Real
end
truncerr(epsilon::Real) = TruncationError(epsilon)
Base.eps(t::TruncationError) = t.epsilon

struct TruncationDimension <: TruncationScheme
    D::Integer
    epsilon::Real
end
truncdim(D::Int) = TruncationDimension(D, 0)
dim(t::TruncationDimension) = t.D
Base.eps(t::TruncationDimension) = t.epsilon

struct TruncationSpace{S<:ElementarySpace} <: TruncationScheme
    space::S
    epsilon::Real
end
truncspace(space::ElementarySpace) = TruncationSpace(space, 0)
space(s::TruncationSpace) = s.space
dim(s::TruncationSpace) = dim(s.space)
Base.eps(t::TruncationSpace) = t.epsilon

#
#
# # in order not having to specify all possibilities:
# Base.:|(trunc1::TruncationScheme, trunc2::TruncationScheme) = |(trunc2, trunc1)
# Base.:&(trunc1::TruncationScheme, trunc2::TruncationScheme) = &(trunc2, trunc1)
# Base.:|(trunc1::NoTruncation, trunc2::NoTruncation) = trunc1
# Base.:|(trunc1::NoTruncation, trunc2::TruncationScheme) = trunc2
#
# Base.:&(trunc1::NoTruncation, trunc2::NoTruncation) = trunc1
# Base.:&(trunc1::NoTruncation, trunc2::TruncationScheme) = trunc1
#
# Base.:|(trunc1::TruncationError, trunc2::TruncationError) = TruncationError(max(eps(trunc1), eps(trunc2)))
# Base.:&(trunc1::TruncationError, trunc2::TruncationError) = TruncationError(min(eps(trunc1), eps(trunc2)))
#
#
# Base.:|(trunc1::TruncationDimension, trunc2::TruncationDimension) = TruncationDimension(min(dim(trunc1), dim(trunc2)), max(eps(trunc1), eps(trunc2)))
# Base.:&(trunc1::TruncationDimension, trunc2::TruncationDimension) = TruncationDimension(max(dim(trunc1), dim(trunc2)), min(eps(trunc1), eps(trunc2)))
#
# Base.:|(trunc1::TruncationDimension, trunc2::TruncationError) = TruncationDimension(dim(trunc1), max(eps(trunc1), eps(trunc2)))
# Base.:&(trunc1::TruncationDimension, trunc2::TruncationError) = TruncationError(min(eps(trunc1), eps(trunc2)))
#
# Base.&(trunc1::TruncationSpace, trunc2::TruncationError) = TruncationSpace(space(trunc1), max(eps(trunc1), eps(trunc2)))
# Base.&(trunc2::TruncationError, trunc1::TruncationSpace) = TruncationSpace(space(trunc1), max(eps(trunc1), eps(trunc2)))


text:
w=readdlm("W.txt")

x=readdlm("x.txt")

y=readdlm("y.txt")

yb=1974
tmax,jmax=size(x)

cst=cohorts_t1(yb,tmax,jmax,x,y,w);

eq_csmtx(csv2cst(cst2csv2(cst)),cst)

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
export PetscMat, PetscMatSetType, PetscSetUp, PetscMatSetValues, PetscMatAssemblyBegin, PetscMatAssemblyEnd, PetscMatSetSizes, PetscMatGetSize, PetscMatGetValues, PetscMatGetOwnershipRange, PetscMatXAIJSetPreallocation, PetscMatMPIAIJSetPreallocation, PetscMatSetFromOptions, PetscMatGetInfo, PetscMatMatMult, PetscMatNorm, PetscMatZeroEntries, PetscMatSetValuesBlocked, MatSetOption, MatCreateShell, MatShellSetOperation, MatShellSetContext, MatShellGetContext, MatGetType, MatCreateTranspose, MatTranspose

type PetscMat  <: AbstractArray{PetscScalar, 2}
  pobj::Ptr{Void}
  function PetscMat(comm::MPI_Comm)
#    comm = PETSC_COMM_SELF();
    vec = Array(Ptr{Void},1)
    err = ccall( (:MatCreate,  libpetsclocation), PetscErrorCode, (comm_type, Ptr{Ptr{Void}}),comm,vec);
    vec = new(vec[1])
#    finalizer(vec,PetscDestroy)
    # does not seem to be called immediately when vec is no longer visible, is it called later during garbage collection?
    return vec
  end

  function PetscMat(pobj::Ptr{Void})  # default constructor
    return new(pobj)
  end
end

function MatCreateShell(arg1::MPI_Comm,arg2::Integer,arg3::Integer,arg4::Integer,arg5::Integer, arg6::Ptr{Void})
  # arg6 is the user provided context
    arg7 = Ref{Ptr{Void}}()
    ccall((:MatCreateShell,petsc),PetscErrorCode,(MPI_Comm,PetscInt,PetscInt,PetscInt,PetscInt,Ptr{Void}, Ref{Ptr{Void}}),arg1,arg2,arg3,arg4,arg5,arg6,arg7)
    return PetscMat(arg7[])
end

function MatShellSetOperation(arg1::PetscMat,arg2::MatOperation,arg3::Ptr{Void})
# arg3 is a function pointer, and must have the signature:
# void fname(Mat, vec, vec) for MATOP_MULT

    ccall((:MatShellSetOperation,petsc),PetscErrorCode,(Ptr{Void},MatOperation,Ptr{Void}),arg1.pobj,arg2,arg3)
end

function MatShellGetContext(arg1::PetscMat)
# get the user provided context for the matrix shell
    arg2 = Ref{Ptr{Void}}()
    ccall((:MatShellGetContext,petsc),PetscErrorCode,(Ptr{Void},Ref{Ptr{Void}}),arg1.pobj,arg2)
    return arg2[]  # turn it into a julia object here?
end

function MatShellSetContext(arg1::PetscMat, arg2::Ptr{Void})
# get the user provided context for the matrix shell
    ccall((:MatShellSetContext,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),arg1.pobj,arg2)
    return nothing  # turn it into a julia object here?
end

function MatGetType(arg1::PetscMat)
    arg2 = Ref{Ptr{UInt8}}()
    ccall((:MatGetType,petsc),PetscErrorCode,(Ptr{Void}, Ref{Ptr{UInt8}}),arg1.pobj,arg2)
    return bytestring(arg2[])
end

function MatCreateTranspose(arg1::PetscMat)
  arg2 = Ref{Ptr{Void}}()

  ccall((:MatCreateTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Ptr{Void}}),arg1.pobj, arg2)

  return PetscMat(arg2[])
end

"""
  Constructs the transpose of a given matrix.  Can be in place or out of place.

  Inputs:
    A: matrix to take the transpose of

  Keywords:
    inplace: whether or not to do the transpose in place

  Outputs:
    A PetscMat object, either the original object A if the transpose was done
    in place, or a new matrix object if the transpose was done out of place
"""
function MatTranspose(A::PetscMat; inplace::Bool=false, mat_initialized=false)
  if inplace
    B = Ref{Ptr{Void}}(A.pobj)
    reuse = MAT_REUSE_MATRIX
  else
    B = Ref{Ptr{Void}}()
    reuse = MAT_INITIAL_MATRIX
  end

  println("before, A.pobj = ", A.pobj)
    ccall((:MatTranspose,petsc),PetscErrorCode,(Ptr{Void},MatReuse,Ptr{Ptr{Void}}),A.pobj, reuse, B)

  println("after, A.pobj = ", A.pobj)
  println("typeof(B) = ", typeof(B))
  println("after, B[] = ", B[])

  if inplace
    return A
  else
    return PetscMat(B[])
  end
end

  function PetscDestroy(vec::PetscMat)
    if (vec.pobj != 0)
      err = ccall( (:MatDestroy,  libpetsclocation), PetscErrorCode, (Ptr{Ptr{Void}},), &vec.pobj);
    end
#    vec.pobj = 0

#    println("Petsc Mat Destroy called")
#    sleep(5)
  end

immutable PetscMatInfo
    block_size::PetscLogDouble
    nz_allocated::PetscLogDouble
    nz_used::PetscLogDouble
    nz_unneeded::PetscLogDouble
    memory::PetscLogDouble
    assemblies::PetscLogDouble
    mallocs::PetscLogDouble
    fill_ratio_given::PetscLogDouble
    fill_ratio_needed::PetscLogDouble
    factor_mallocs::PetscLogDouble

    function PetscMatInfo()  # incomplete initialization
      return new()
    end
end

function show(io::IO, obj::PetscMatInfo)
# print the fields of PetscMatInfo
#  println("PetscMatInfo:")
  println(io, "  block_size : ", obj.block_size)
  println(io, "  nz_allocated : ", obj.nz_allocated)
  println(io, "  nz_used : ", obj.nz_used)
  println(io, "  nz_unneeded : ", obj.nz_unneeded)
  println(io, "  memory : ", obj.memory)
  println(io, "  assemblies : ", obj.assemblies)
  println(io, "  mallocs : ", obj.mallocs)
  println(io, "  fill_ratio_given : ", obj.fill_ratio_given)
  println(io, "  fill_ratio_needed : ", obj.fill_ratio_needed)
  println(io, "  factor_mallocs : ", obj.factor_mallocs)
end

  function PetscMatSetFromOptions(mat::PetscMat)
    ccall((:MatSetFromOptions,petsc),PetscErrorCode,(Ptr{Void},), mat.pobj)

  end

  function PetscMatSetType(vec::PetscMat,name)
    err = ccall( (:MatSetType,  libpetsclocation), PetscErrorCode,(Ptr{Void}, Cstring), vec.pobj,name);

  end

  function PetscSetUp(vec::PetscMat)
    err = ccall( ( :MatSetUp,  libpetsclocation), PetscErrorCode, (Ptr{Void},), vec.pobj);

    MatSetOption(vec, MAT_ROW_ORIENTED, PETSC_TRUE)  # julia data is column-major
  end

#=
  PETSC_MAT_FLUSH_ASSEMBLY = 1;
  PETSC_MAT_FINAL_ASSEMBLY = 0
=#

  function PetscMatSetValues(vec::PetscMat,idi::Array{PetscInt},idj::Array{PetscInt},array::Array{PetscScalar},flag::Integer)
    # remember, only matrices can be inserted into a Petsc matrix
    # if array is a 3 by 3, then idi and idj are vectors of length 3
#    idi = idi
#    idj = idj

    @assert length(idi)*length(idj) == length(array)

    # do check here to ensure array is the right shape (remember tranpose)
    err = ccall( ( :MatSetValues,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar},Int32), vec.pobj,length(idi), idi, length(idj), idj,array,flag);
#    idi = idi
#    idj = idj
    return err
  end

  function PetscMatSetValuesBlocked(mat::PetscMat, idi::Array{PetscInt}, idj::Array{PetscInt}, v::Array{PetscScalar}, flag::Integer)

      err = ccall((:MatSetValuesBlocked, libpetsclocation), PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},PetscInt,Ptr{PetscInt},Ptr{PetscScalar}, Int32),mat.pobj, length(idi), idi, length(idj), idj, v,flag)

      return err
  end

  function MatSetOption(mat::PetscMat,arg2::MatOption,arg3::Union{Integer, Bool})
      println("mat.pobj = ", mat.pobj)
      ierr = ccall((:MatSetOption, libpetsclocation),PetscErrorCode, (Ptr{Void},MatOption,PetscBool), mat.pobj, arg2, arg3)

      if ierr != 0
	println(STDERR, "Error: MatSetOption returned non zero exit status")
      end
  end

  # 1-based indexing for both Petsc and regular matrices
  #----------------------------------------------------------------------------
  """
    1-based indexing for both regular and Pets matrices.
    The flag must be either PETSC_INSERT_VALUES or PETSC_ADD_VALUES.

    idxm and idxn cannot alias
  """
  function set_values1!(mat::PetscMat, idxm::Array{PetscInt}, idxn::Array{PetscInt},
                        vals::Array{PetscScalar}, flag::Integer=PETSC_INSERT_VALUES)

    for i=1:length(idxm)
      idxm[i] -= 1
    end

    for i=1:length(idxn)
      idxn[i] -= 1
    end

    err = PetscMatSetValues(mat, idxm, idxn, vals, flag)

    for i=1:length(idxm)
      idxm[i] += 1
    end

    for i=1:length(idxn)
      idxn[i] += 1
    end

    return err
  end

  function set_values1!{T}(mat::AbstractMatrix, idxm::Array{PetscInt}, idxn::Array{PetscInt},
                           vals::Array{T}, flag::Integer=PETSC_INSERT_VALUES)

    if flag == PETSC_INSERT_VALUES
      for i=1:length(idxn)
        for j=1:length(idxm)
          mat[idxm[j], idxn[i]] = vals[j, i]
        end
      end
    elseif flag == PETSC_ADD_VALUES
      for i=1:length(idxn)
        for j=1:length(idxm)
          mat[idxm[j], idxn[i]] += vals[j, i]
        end
      end
    end

    return PetscErrorCode(0)
  end

  function set_values1!{T}(mat::SparseMatrixCSC, idxm::Array{PetscInt}, idxn::Array{PetscInt},
                           vals::Array{T}, flag::Integer=PETSC_INSERT_VALUES)

    if flag == PETSC_INSERT_VALUES
      println("STRANGE: inside PETSC_INSERT_VALUES")
      for i=1:length(idxn)
        for j=1:length(idxm)
          mat[idxm[j], idxn[i]] = vals[j, i]
        end
      end
    elseif flag == PETSC_ADD_VALUES  # optimized += implementation
      # hoist
      colptr = mat.colptr
      rowval = mat.rowval
      nzval = mat.nzval
      for i=1:length(idxn)
        row_start = colptr[idxn[i]]
        row_end = colptr[idxn[i]+1] - 1
        rowvals_extract = unsafe_view(rowval, row_start:row_end)
        for j=1:length(idxm)
          idx = searchsortedfirst(rowvals_extract, idxm[j])
          idx = row_start + idx - 1
          nzval[idx] += vals[j, i]
        end
      end
    end

    return PetscErrorCode(0)
  end

  function get_values1!(mat::PetscMat, idxm::Array{PetscInt}, idxn::Array{PetscInt},
                        vals::Array{PetscScalar})

    for i=1:length(idxm)
      idxm[i] -= 1
    end

    for i=1:length(idxn)
      idxn[i] -= 1
    end

    err = PetscMatGetValues(mat, idxm, idxn, vals)

    for i=1:length(idxm)
      idxm[i] += 1
    end

    for i=1:length(idxn)
      idxn[i] += 1
    end

    return err
  end

  function get_values1!(mat::AbstractMatrix, idxm::Array{PetscInt}, idxn::Array{PetscInt},
                        vals::Array)

    for i=1:length(idxn)
      for j=1:length(idxm)
        vals[j, i] = mat[idxm[j], idxn[i]]
      end
    end

    return PetscErrorCode(0)
  end

  function PetscMatAssemblyBegin(obj::AbstractArray, flg::Integer=0)
  end

  function PetscMatAssemblyBegin(obj::PetscMat,flg::Integer)
    err = ccall( ( :MatAssemblyBegin,  libpetsclocation), PetscErrorCode,(Ptr{Void},Int32), obj.pobj,flg);
  end

  function PetscMatAssemblyBegin(obj::PetscMat)
    return PetscMatAssemblyBegin(obj,PETSC_MAT_FINAL_ASSEMBLY);
  end

  function PetscMatAssemblyEnd(obj::AbstractArray, flag=0)
  end

  function PetscMatAssemblyEnd(obj::PetscMat,flg::Integer)
    err = ccall( ( :MatAssemblyEnd,  libpetsclocation), PetscErrorCode,(Ptr{Void},Int32), obj.pobj,flg);
  end

  function PetscMatAssemblyEnd(obj::PetscMat)
    return PetscMatAssemblyEnd(obj,PETSC_MAT_FINAL_ASSEMBLY);
  end

  function PetscMatSetSizes(vec::PetscMat,m::PetscInt, n::PetscInt, M::PetscInt, N::PetscInt)
    err = ccall( ( :MatSetSizes,  libpetsclocation), PetscErrorCode, (Ptr{Void},PetscInt, PetscInt, PetscInt, PetscInt), vec.pobj,m,n,M,N);

    return err
  end

  function PetscView(obj::PetscMat,viewer)
    err = ccall( (:MatView,  libpetsclocation), PetscErrorCode, (Ptr{Void}, Int64),obj.pobj,0);
  end

  function PetscMatGetSize(obj::PetscMat)
    m = Array(PetscInt, 1)
    n = Array(PetscInt, 1)
    err = ccall(Libdl.dlsym(libpetsc, :MatGetSize), PetscErrorCode,(Ptr{Void}, Ptr{PetscInt},Ptr{PetscInt}), obj.pobj,m,n);
    return m[1],n[1]
  end

export PetscMatGetLocalSize
function PetscMatGetLocalSize(mat::PetscMat)
    m = Array(PetscInt, 1)
    n = Array(PetscInt, 1)
    ccall((:MatGetLocalSize,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt}),mat.pobj, m, n)
    return m[1], n[1]
end

# Petsc populates v row major, so Julia will see the returned array as being transposed
function PetscMatGetValues(obj::PetscMat, idxm::Array{PetscInt, 1}, idxn::Array{PetscInt, 1}, v::Array{PetscScalar, 2})
    # do check here to ensure v is the right shape

    ccall((:MatGetValues,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},PetscInt,Ptr{PetscInt},Ptr{PetscScalar}), obj.pobj, length(idxm), idxm, length(idxn), idxn, v)

end

function PetscMatGetOwnershipRange(mat::PetscMat)
    low = Array(PetscInt,1)
    high = Array(PetscInt,1)
    ccall((:MatGetOwnershipRange,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt}),mat.pobj, low, high)

    return low[1], high[1]
end

function PetscMatZeroEntries(mat::PetscMat)
    ccall((:MatZeroEntries,petsc),PetscErrorCode,(Ptr{Void},),mat.pobj)
end

function PetscMatZeroEntries(mat::AbstractMatrix)
  fill!(mat, 0.0)
end

function PetscMatZeroEntries(mat::SparseMatrixCSC)
  fill!(mat.nzval, 0.0)
end

### new function ###

function PetscMatNorm(mat::PetscMat, ntype::NormType)
    nrm = Array(PetscReal, 1)
    ccall((:MatNorm,petsc),PetscErrorCode,(Ptr{Void}, NormType,Ptr{PetscReal}), mat.pobj, ntype, nrm)
    return nrm[1]
end

export PetscMatAXPY, PetscMatAYPX, PetscMatScale, PetscMatShift
function PetscMatAXPY(Y::PetscMat, a::PetscScalar, X::PetscMat, str::PetscMatStructure)
    ccall((:MatAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void},PetscMatStructure), Y.pobj, a, X.pobj, str)
end

function PetscMatAYPX( Y::PetscMat, a::PetscScalar, X::PetscMat, str::PetscMatStructure)
    ccall((:MatAYPX,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void}, PetscMatStructure), Y.pobj, a, X.pobj, str)
end

function PetscMatScale(mat::PetscMat, a::Number)
    ccall((:MatScale,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), mat.pobj, PetscScalar(a))
end

function PetscMatScale(mat::AbstractArray, a::Number)
  scale!(mat, PetscScalar(a))
end

function PetscMatScale(mat::SparseMatrixCSC, a::Number)
  scale!(mat.nzval, PetscScalar(a))
end

function PetscMatShift(mat::PetscMat, a::PetscScalar)
    ccall((:MatShift,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), mat.pobj, a)
end

function PetscMatShift(mat::AbstractMatrix, a::PetscScalar)

  for i=1:size(mat, 1)
    mat[i, i] += a
  end

  return nothing
end

export PetscMatMult, PetscMatMultAdd, PetscMatMultTranspose, PetscMatMultHermitianTranspose

function PetscMatMult(mat::PetscMat, x::PetscVec, y::PetscVec)
    ccall((:MatMult,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), mat.pobj, x.pobj, y.pobj)
end

function PetscMatMultAdd(A::PetscMat, x::PetscVec, y::PetscVec, z::PetscVec)
    ccall((:MatMultAdd,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}), A.pobj, x.pobj, y.pobj, z.pobj)
end

function PetscMatMultTranspose(A::PetscMat, x::PetscVec, y::PetscVec)
    ccall((:MatMultTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), A.pobj, x.pobj, y.pobj)
end

function PetscMatMultHermitianTranspose(mat::PetscMat, x::PetscVec, y::PetscVec)
    ccall((:MatMultHermitianTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}),mat.pobj, x.pobj, y.pobj)
end

function PetscMatMatMult(A::PetscMat, B::PetscMat, scall::MatReuse, fill::PetscReal, C::PetscMat)
    arr = [C.pobj]
    ccall((:MatMatMult,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}, MatReuse, PetscReal, Ptr{Ptr{Void}}), A.pobj, B.pobj, scall, fill, arr)
    if C.pobj != arr[1]
      PetscDestroy(C)
      C.pobj = arr[1]
    end

end

function PetscMatXAIJSetPreallocation(mat::PetscMat, bs::PetscInt, dnnz::AbstractArray{PetscInt, 1}, onnz::AbstractArray{PetscInt,1}, dnnzu::AbstractArray{PetscInt, 1}, onnzu::AbstractArray{PetscInt, 1})
# this is a unified interface for matrix preallocation for the Petsc built in
# matrix types: Aij, Bij, and their respective symmetric forms SAij, SBij
# for a non symmetric format matrix, dnnzu and onnzu are not required
# and vice versa for a non symmetric format matrix

    ccall((:MatXAIJSetPreallocation,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},Ptr{PetscInt},Ptr{PetscInt},Ptr{PetscInt}), mat.pobj, bs, dnnz, onnz, dnnzu, onnzu)
end

# Matrix preallocation
function PetscMatMPIAIJSetPreallocation(mat::PetscMat, d_nz::PetscInt, d_nnz::PetscInt_arr_or_null, o_nz::PetscInt, o_nnz::PetscInt_arr_or_null)

    ccall((:MatMPIAIJSetPreallocation,petsc),PetscErrorCode,(Ptr{Void}, PetscInt,Ptr{PetscInt},PetscInt,Ptr{PetscInt}), mat.pobj, d_nz, d_nnz, o_nz, o_nnz)
end

function PetscMatGetInfo(mat::PetscMat, info_type::Int32)
    matinfo = PetscMatInfo()  # create uninitialized struct
    ref_matinfo = Ref{PetscMatInfo}(matinfo)
    ccall((:MatGetInfo,petsc),PetscErrorCode,(Ptr{Void}, Int32,Ref{PetscMatInfo}), mat.pobj ,info_type, ref_matinfo)
    return ref_matinfo[]
end

# some functions to make a PetscMat look like a Julia array
# these only work after a matrix has been assembled
# some of them only work for uniprocess matrices
import Base.size
import Base.getindex

function size(A::PetscMat)

  return PetscMatGetLocalSize(A)
end

function size(A::PetscMat, dim::Integer)
  dims = PetscMatGetLocalSize(A)
  return dims[dim]
end

function getindex(A::PetscMat, i::Integer, j::Integer)
  # not efficient

  i_ = [PetscInt(i - 1)]
  j_ = [PetscInt(j - 1)]
  val = Array(PetscScalar, 1, 1)

  PetscMatGetValues(A, i_, j_, val)

  return val[1]
end

function getindex(A::PetscMat, idx::Integer)
  # linear indexing
  m,n = PetscMatGetLocalSize(A)
  i = idx % m
  j = div(idx, m)

  # zero based indexing
  i -= 1
  j -= 1

  i_ = [PetscInt(i)]
  j_ = [PetscInt(j)]
  val = Array(PetscScalar, 1, 1)

  PetscMatGetValues(A, i_, j_, val)
  return val[1]
end


text:
"""
    term(t::AbstractTheory, root::Σ, args::Vector{AbstractTerm}) -> AbstractTerm

Produce a term in theory `t` with root `root` and arguments `args`. If an unexpected input
is provided, raise an exception.
"""
function term end

"""
    theory(T::Type{<:AbstractTerm}) -> AbstractTheory

Produce the theory which contains type `T`.
"""
function theory end
theory(t::T) where {T<:AbstractTerm} = theory(T)

"""
    vars(t::AbstractTerm) -> Set{Variable}

Produce the set of variables which appear as subterms of `t`.
"""
function vars end

"""
    priority(::Type{<:AbstractTerm}) -> Int

Priority of a term type, used to produce a total ordering over terms.

!!! note

    The `priority` function must be injective; in other words, every type must map to a
    distinct priority.
"""
function priority end
priority(t::T) where {T<:AbstractTerm} = priority(T)

"""
    s::AbstractTerm >ₜ t::AbstractTerm -> Bool

Total ordering on terms.
"""
function >ₜ end
>ₜ(s::AbstractTerm, t::AbstractTerm) = priority(s) > priority(t)

"""
    matcher(t::AbstractTerm, V::Set{Variable}) -> Tuple{AbstractMatcher,Set{Variable}}

Generate a matcher for `t`, given that variables `V` will already be matched. Produce a set
of variables which are guaranteed to be fixed during matching.
"""
function matcher end

"""
    match!(σ, pattern::AbstractMatcher, term::AbstractTerm) -> Union{AbstractSubproblem,Nothing}

Match `term` against `pattern` given the partial substitution `σ`, mutating `σ` and
producing a subproblem to solve or producing `nothing` if a match is impossible.

!!! note

    Unless `pattern` and `term` are derived from the same theory, `match!` should
    necessarily produce `nothing`.
"""
match!(::Any, ::AbstractMatcher, ::AbstractTerm) = nothing

"""
    compile(pattern::AbstractMatcher, V) -> Tuple{Symbol,Expr}

Compile `pattern` to a native `Expr` containing a function for matching against `pattern`,
given that the variables in `V` are fixed.
"""
compile(pattern::AbstractMatcher) = compile(pattern, Set{Variable}())

"""
    map(f, t::AbstractTerm) -> AbstractTerm

Call `f` on each subterm of `t`.
"""
Base.map(f, t::AbstractTerm)

"""
    rewriter(t::AbstractTheory) -> AbstractRewriter

Produce a fresh rewriter for theory `t`.
"""
function rewriter end

"""
    rewrite(rw::AbstractRewriter, t::AbstractTerm) -> AbstractTerm

Rewrite `t` using `rw`, producing `nothing` if the process fails.
"""
function rewrite end

"""
    compile(rw::AbstractRewriter) -> Tuple{Symbol,Expr}

Compile `rw` to a native `Expr` containing a function for rewriting with respect to `rw`.
If `name` is provided, the function will be defined as `name`.
"""
compile(rw::AbstractRewriter)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Defines a set of pods (namely those matching the labelSelector relative to the given namespace(s)) that this pod should be co-located (affinity) or not co-located (anti-affinity) with, where co-located is defined as running on a node whose value of the label with key <topologyKey> matches that of any node on which a pod of the set of pods is running

    IoK8sApiCoreV1PodAffinityTerm(;
        labelSelector=nothing,
        namespaces=nothing,
        topologyKey=nothing,
    )

    - labelSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : A label query over a set of resources, in this case pods.
    - namespaces::Vector{String} : namespaces specifies which namespaces the labelSelector applies to (matches against); null or empty list means \"this pod's namespace\"
    - topologyKey::String : This pod should be co-located (affinity) or not co-located (anti-affinity) with the pods matching the labelSelector in the specified namespaces, where co-located is defined as running on a node whose value of the label with key topologyKey matches that of any node on which any of the selected pods is running. Empty topologyKey is not allowed.
"""
mutable struct IoK8sApiCoreV1PodAffinityTerm <: SwaggerModel
    labelSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: labelSelector
    namespaces::Any # spec type: Union{ Nothing, Vector{String} } # spec name: namespaces
    topologyKey::Any # spec type: Union{ Nothing, String } # spec name: topologyKey

    function IoK8sApiCoreV1PodAffinityTerm(;labelSelector=nothing, namespaces=nothing, topologyKey=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodAffinityTerm, Symbol("labelSelector"), labelSelector)
        setfield!(o, Symbol("labelSelector"), labelSelector)
        validate_property(IoK8sApiCoreV1PodAffinityTerm, Symbol("namespaces"), namespaces)
        setfield!(o, Symbol("namespaces"), namespaces)
        validate_property(IoK8sApiCoreV1PodAffinityTerm, Symbol("topologyKey"), topologyKey)
        setfield!(o, Symbol("topologyKey"), topologyKey)
        o
    end
end # type IoK8sApiCoreV1PodAffinityTerm

const _property_map_IoK8sApiCoreV1PodAffinityTerm = Dict{Symbol,Symbol}(Symbol("labelSelector")=>Symbol("labelSelector"), Symbol("namespaces")=>Symbol("namespaces"), Symbol("topologyKey")=>Symbol("topologyKey"))
const _property_types_IoK8sApiCoreV1PodAffinityTerm = Dict{Symbol,String}(Symbol("labelSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("namespaces")=>"Vector{String}", Symbol("topologyKey")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PodAffinityTerm }) = collect(keys(_property_map_IoK8sApiCoreV1PodAffinityTerm))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodAffinityTerm }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodAffinityTerm[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodAffinityTerm }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodAffinityTerm[property_name]

function check_required(o::IoK8sApiCoreV1PodAffinityTerm)
    (getproperty(o, Symbol("topologyKey")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodAffinityTerm }, name::Symbol, val)
end


text:
# This file includes testset for sources.

@testset "ClockTestSet" begin
    @info "Running ClockTestSet ..."

    # Clock construction
    Clock(0 : 1)
    Clock(0. : 2. : 10.)
    Clock(sort(rand(10)))
    Clock(
        Channel(0) do ch
            for t in 1 : 10
                put!(ch, t)
            end
        end
    )
    Clock(0, 2, 10)

    # Check Clock defaults.
    clk = Clock(0. : 5.)
    vals = Float64[]
    for val in clk
        push!(vals, val)
    end
    @test vals == collect(0. : 5.)

    # Pausing Clock
    clk = Clock(0 : 1 : 10)
    item, state = iterate(clk)
    @test (item, state) == (0, 0)
    item, state = iterate(clk, state)
    @test (item, state) == (1, 1)
    pause!(clk)
    iter = iterate(clk)
    @test iter === nothing

    @info "Done ClockTestSet."
end  # testset

@testset "GeneratorsTestSet" begin
    @info "Running GeneratorsTestSet ..."
    # FunctionGenerator construction
    gen = SinewaveGenerator()
    @test typeof(gen.trigger) == Inpin{Float64}
    @test typeof(gen.handshake) == Outpin{Bool}
    @test !hasfield(typeof(gen), :input)
    @test typeof(gen.output) == Outport{Outpin{Float64}}

    @def_source struct  Mygen{OP, RO} <: AbstractSource
        output::OP = Outport(2)
        readout::RO = t -> [sin(t), cos(t)]
    end
    gen = Mygen()
    @test length(gen.output) == 2

    # Driving FunctionGenerator
    gen = Mygen(readout = t -> t, output=Outport(1))
    trg = Outpin()
    hnd = Inpin{Bool}()
    ip = Inport()
    connect!(gen.output, ip)
    connect!(trg, gen.trigger)
    connect!(gen.handshake, hnd)
    task = launch(gen)
    task2 = @async while true
        all(take!(ip) .=== NaN) && break
    end
    for t in 1. : 10.
        put!(trg, t)
        take!(hnd)
    end
    @test content(ip[1].link.buffer) == collect(1:10)
    @test !istaskdone(task)
    put!(trg, NaN)
    @test istaskdone(task)
    put!(gen.output, [NaN])
    @test istaskdone(task2)

    # Construction of other generators
    sinegen = SinewaveGenerator()
    dampedsinegen = DampedSinewaveGenerator()
    sqauregen = SquarewaveGenerator()
    trigen = TriangularwaveGenerator()
    congen = ConstantGenerator()
    rampgen = RampGenerator()
    stepgen = StepGenerator()
    expgen = ExponentialGenerator()
    dampedexpgen = DampedExponentialGenerator()

    # Mutaton of generators
    @test_throws Exception sinegen.amplitude = 5.
    @test_throws Exception sqauregen.high = 10.

    # Test redefining new source types.
    @test_throws Exception @eval @def_source struct Mygen{RO,OP}
        readout::RO = t -> sin(t)
        output::OP = Outport()
    end

    @test_throws Exception @eval @def_source struct Mygen{RO,OP} <: SomeDummyType
        readout::RO = t -> sin(t)
        output::OP = Outport()
    end

    @info "Done GeneratorsTestSet ..."
end  # testset


text:
#=

    This is a port of the R2 model ClinicalTrial.cs

    with the datafiles
    - clinical_trial_control.csv  (size 1000)
    - clinical_trial_treatment.csv (size 1000)

    length :1000 mean(control): 0.513 mean(treatment): 0.51

    Output from the R2 model ClinicalTrial.cs (for isEffective)
    ```
    Mean: 0.113
    Variance: 0.100331
    Number of accepted samples = 318
    0:00:14:048.362.900
    ```

    Output from this Turing.jl model:

Summary Statistics
     parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

    isEffective    0.0659    0.2481     0.0025    0.0140   246.7994    1.0080       25.5512
    probControl    0.5152    0.2910     0.0029    0.0184   182.4067    0.9999       18.8846
  probTreatment    0.5278    0.2789     0.0028    0.0180   173.0546    1.0028       17.9164
        probAll    0.4985    0.0829     0.0008    0.0050   223.0815    1.0112       23.0957

Quantiles
     parameters      2.5%     25.0%     50.0%     75.0%     97.5%
         Symbol   Float64   Float64   Float64   Float64   Float64

    isEffective    0.0000    0.0000    0.0000    0.0000    1.0000
    probControl    0.0225    0.2196    0.5287    0.7653    0.9847
  probTreatment    0.0196    0.3207    0.5222    0.8098    0.9720
        probAll    0.1196    0.5019    0.5118    0.5190    0.5348

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function clinical_trial(control,treatment)
    n = length(control)

    isEffective   ~ flip(0.5)
    probControl   ~ Beta(1,1)
    probTreatment ~ Beta(1,1)
    probAll       ~ Beta(1,1)

    if isEffective
        control   ~ filldist(Bernoulli(probControl),n)
        treatment ~ filldist(Bernoulli(probTreatment),n)
    else
        control   ~ filldist(Bernoulli(probAll),n)
        treatment ~ filldist(Bernoulli(probAll),n)

    end
end

control = parse.(Int,split(readline("clinical_trial_control.csv"),","))
treatment = parse.(Int,split(readline("clinical_trial_treatment.csv"),","))
println("len(control):$(length(control)) len(treatment):$(length(treatment)) mean(control):$(mean(control)) mean(treatment): $(mean(treatment))")
# println("treatment:$treatment")

model = clinical_trial(control, treatment)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns[[:isEffective,:probControl,:probTreatment,:probAll]])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AccessUriOutput <: SwaggerModel
    output::Any # spec type: Union{ Nothing, AccessUriRaw } # spec name: output

    function AccessUriOutput(;output=nothing)
        o = new()
        validate_property(AccessUriOutput, Symbol("output"), output)
        setfield!(o, Symbol("output"), output)
        o
    end
end # type AccessUriOutput

const _property_map_AccessUriOutput = Dict{Symbol,Symbol}(Symbol("output")=>Symbol("output"))
const _property_types_AccessUriOutput = Dict{Symbol,String}(Symbol("output")=>"AccessUriRaw")
Base.propertynames(::Type{ AccessUriOutput }) = collect(keys(_property_map_AccessUriOutput))
Swagger.property_type(::Type{ AccessUriOutput }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_AccessUriOutput[name]))}
Swagger.field_name(::Type{ AccessUriOutput }, property_name::Symbol) =  _property_map_AccessUriOutput[property_name]

function check_required(o::AccessUriOutput)
    true
end

function validate_property(::Type{ AccessUriOutput }, name::Symbol, val)
end


text:
@testset "eachleaf!" begin
    # test with/without threads
    for node in (Node{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2, 3}(),
                 Node{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2, 1}(),
                 HashNode{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2, 3}(),
                 HashNode{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2, 1}(),
                 DynamicNode{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2}(8,8),
                 DynamicNode{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2}(2,2),
                 DynamicHashNode{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2}(8,8),
                 DynamicHashNode{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2}(2,2),)
        A = TreeView(node)
        A .= 1

        # eachleaf! for all
        eachleaf!(x -> 3x, A)
        @test all(==(3), A)

        # eachleaf! for part matrix
        eachleaf!(x -> 2x, A, :, 1:3)
        @test all(==(6), A[:, 1:3])
        @test all(==(3), A[:, 4:end])
    end
end