text:
mutable struct ConstantMean{T<:Real} <: MeanPrior{T}
    C::T
    opt::Optimizer
end

"""
**ConstantMean**
```julia
ConstantMean(c::T=1.0;opt::Optimizer=Adam(α=0.01))
```

Construct a prior mean with constant `c`
Optionally set an optimizer `opt` (`Adam(α=0.01)` by default)
"""
function ConstantMean(c::T=1.0;opt::Optimizer=Adam(α=0.01)) where {T<:Real}
    ConstantMean{T}(c,opt)
end

function update!(μ::ConstantMean{T},grad::AbstractVector{T}) where {T<:Real}
    μ.C += update(μ.opt,sum(grad))
end

Base.:+(x::Real,y::ConstantMean{<:Real}) = x+y.C
Base.:+(x::AbstractVector{<:Real},y::ConstantMean{<:Real}) = x.+y.C
Base.:+(x::ConstantMean{<:Real},y::Real) = y+x.C
Base.:+(x::ConstantMean{<:Real},y::AbstractVector{<:Real}) = y.+x.C
Base.:+(x::ConstantMean{<:Real},y::ConstantMean{<:Real}) = ConstantMean(x.C+y.C)
Base.:-(x::Real,y::ConstantMean) = x - y.C
Base.:-(x::AbstractVector{<:Real},y::ConstantMean) = x .- y.C
Base.:-(x::ConstantMean{<:Real},y::Real) = x.C - y
Base.:-(x::ConstantMean{<:Real},y::AbstractVector{<:Real}) = x.C .- y
Base.:-(x::ConstantMean{<:Real},y::ConstantMean{<:Real}) = ConstantMean(x.C-y.C)
Base.:*(A::AbstractMatrix{<:Real},y::ConstantMean{T}) where {T<:Real} = y.C*A*ones(T,size(A,2))
Base.:*(y::ConstantMean{T},A::AbstractMatrix{<:Real}) where {T<:Real} = y.C*ones(T,1,size(A,1))*A
Base.:convert(::T1,x::ConstantMean{T2}) where {T1<:Real,T2<:Real} = T1(x.C)


text:
import ThorinDistributions
using DynamicPolynomials, Random
Random.seed!(123)

function get_hankels_truth(s,a,b)
    m = length(s)-1
    if iseven(m)
        # m = 2n
        n = Int(m//2)
        H_low = [s[i+j-1] for i in 1:n+1,j in 1:n+1]
        H_up = [(a+b)*s[i+j] - s[i+j+1] - a*b*s[i+j-1] for i in 1:n,j in 1:n]
    else
        # m = 2n+1
        n = Int((m-1)//2)
        H_low = [s[i+j] - a*s[i+j-1] for i in 1:n+1,j in 1:n+1]
        H_up = [b*s[i+j-1] - s[i+j] for i in 1:n+1,j in 1:n+1]
    end
    return H_low,H_up
end

@testset "HankelMatrices.jl tests" begin

    m=20
    @polyvar x[1:m+1]
    H_low, H_up = vec.(ThorinDistributions.get_hankels(x,2,3))
    H_low_op, H_up_op = ThorinDistributions.vec_H.(ThorinDistributions.hankel_operators(m,2,3))

    @testset "coherence of hankel_operator" begin
        @test H_low_op*x == H_low
        @test H_up_op*x == H_up
    end

    a = rand()
    b = rand()
    a,b = min(a,b),max(a,b)

    @testset "coherence of get_hankels for random bounds" begin
        @test get_hankels_truth(x,3,4) == ThorinDistributions.get_hankels(x,3,4)
        @test get_hankels_truth(x,a,b) == ThorinDistributions.get_hankels(x,a,b)
    end

    x = randn(100)

    @testset "coherence for random values" begin
        @test get_hankels_truth(x,0,1) == ThorinDistributions.get_hankels(x,0,1)
        @test get_hankels_truth(x,a,b) == ThorinDistributions.get_hankels(x,a,b)
    end
end

ext: jl


text:
using TextAnalysis, Languages, Clustering

filename = joinpath(dirname(@__FILE__), "data", "highlights.csv")

function get_headers(source)
    headers = [ String(s) for s in source[1,:] ]
    source = source[2:end,:] # clear the headers
    headers, source
end

source = readcsv(filename)
headers, source = get_headers(source)

sample = source[:,5]
length(sample)

#fd = FileDocument(sample)
#sd = StringDocument(fd)

text = join(sample, ",")
sd = StringDocument(text)
crps = Corpus([sd])

remove_punctuation!(sd)
update_lexicon!(crps)
lexicon(crps)

update_inverse_index!(crps)
inverse_index(crps)
hash_function!(crps, TextHashFunction())

crps["Handelskrieg"]

m = DocumentTermMatrix(crps)
D = dtm(m, :dense)
T = tf_idf(D)
cl = kmeans(T, 5)

m = DocumentTermMatrix(crps)
k = 2            # number of topics
iteration = 1000 # number of gibbs sampling iterations
alpha = 0.1      # hyper parameter
beta = 0.1       # hyber parameter
l = lda(m, k, iteration, alpha, beta) # l is k x word matrix.
                                      # value is probablity of occurrence of a word in a topic.

max_stars_repo_name: davidemorotti/newsradar

max_issues_repo_name: davidemorotti/newsradar

max_forks_repo_name: davidemorotti/newsradar


text:
abstract type ExecutionContext end

# These are "soft" interfaces, not all of these methods apply to every subtype of ExecutionContext.
increment!(ctx::T, w::Float64) where T <: ExecutionContext = ctx.weight += w
get_subselection(ctx::T, addr) where T <: ExecutionContext = get_sub(ctx.select, addr)
get_subparameters(ctx::T, addr) where T <: ExecutionContext = get_sub(ctx.params, addr)
visit!(ctx::T, addr) where T <: ExecutionContext = visit!(ctx.visited, addr)
get_prev(ctx::T, addr) where T <: ExecutionContext = get_sub(ctx.prev, addr)
function add_choice!(ctx::T, addr, cs::ChoiceSite) where T <: ExecutionContext
    ctx.score += get_score(cs)
    add_choice!(ctx.tr, addr, cs)
end
function add_call!(ctx::T, addr, cs::CallSite) where T <: ExecutionContext
    ctx.score += get_score(cs)
    add_call!(ctx.tr, addr, cs)
end
function add_call!(ctx::T, addr, cs::CallSite, sc::Float64) where T <: ExecutionContext
    ctx.score += get_score(cs) + sc
    add_call!(ctx.tr, addr, cs)
end
function add_call!(ctx::T, cs::CallSite) where T <: ExecutionContext
    ctx.score += get_score(cs)
    # TODO: should only work for VectorizedTraces - make error explicit here.
    add_call!(ctx.tr, cs)
end

@dynamo function (mx::ExecutionContext)(a...)
    ir = IR(a...)
    ir == nothing && return
    recur!(ir)
    return ir
end

function (mx::ExecutionContext)(::typeof(Core._apply_iterate), f, c::typeof(rand), args...)
    return mx(c, flatten(args)...)
end

# ------------ includes ------------ #

include("contexts/generate.jl")
include("contexts/simulate.jl")
include("contexts/propose.jl")
include("contexts/score.jl")
include("contexts/update.jl")
include("contexts/regenerate.jl")
include("contexts/backpropagate.jl")


text:
import SchumakerSpline.evaluate

function evaluate(f::Piecewise_Function, point::Real)
    which_function = searchsortedlast(f.starts_, point)
    return evaluate(f.functions_[which_function], point)
end
function (s::Piecewise_Function)(x::Union{Real,Date,DateTime,DatePeriod})
    return evaluate(s, x)
end

function derivative(f::Piecewise_Function)
    derivatives = derivative.( f.functions_)
    return Piecewise_Function(f.starts_, derivatives)
end

function indefinite_integral(f::Piecewise_Function)
    indef_integrals = indefinite_integral.( f.functions_)
    return Piecewise_Function(f.starts_, indef_integrals)
end

function +(f::Piecewise_Function,number::Real)
    functions_with_addition = f.functions_ .+ number
    return Piecewise_Function(f.starts_, functions_with_addition)
end
function -(f::Piecewise_Function,number::Real)
    functions_with_subtraction = f.functions_ .- number
    return Piecewise_Function(f.starts_, functions_with_subtraction)
end
function *(f::Piecewise_Function,number::Real)
    functions_with_multiplication = f.functions_ .* number
    return Piecewise_Function(f.starts_, functions_with_multiplication)
end
function /(f::Piecewise_Function,number::Real)
    functions_with_division = f.functions_ ./ number
    return Piecewise_Function(f.starts_, functions_with_division)
end

function create_common_pieces(f1::Piecewise_Function,f2::Piecewise_Function)
    starts_ = unique(sort!(vcat(f1.starts_,f2.starts_)))
    functions1_ = Array{UnivariateFunction}(undef, length(starts_))
    functions2_ = Array{UnivariateFunction}(undef, length(starts_))
    for i in 1:length(starts_)
        point = starts_[i]
        which_1 = searchsortedlast(f1.starts_, point)
        functions1_[i] = f1.functions_[which_1]
        which_2 = searchsortedlast(f2.starts_, point)
        functions2_[i] = f2.functions_[which_2]
    end
    return Piecewise_Function(starts_, functions1_) , Piecewise_Function(starts_, functions2_)
end

function +(f1::Piecewise_Function,f2::Piecewise_Function)
    c_f1, c_f2 = create_common_pieces(f1,f2)
    starts_ = c_f1.starts_
    number_of_pieces = length(starts_)
    functions_ = Array{UnivariateFunction}(undef, number_of_pieces)
    for i in 1:number_of_pieces
        functions_[i] = c_f1.functions_[i] + c_f2.functions_[i]
    end
    return Piecewise_Function(starts_, functions_)
end
function *(f1::Piecewise_Function,f2::Piecewise_Function)
    c_f1, c_f2 = create_common_pieces(f1,f2)
    starts_ = c_f1.starts_
    number_of_pieces = length(starts_)
    functions_ = Array{UnivariateFunction}(undef, number_of_pieces)
    for i in 1:number_of_pieces
        functions_[i] = c_f1.functions_[i] * c_f2.functions_[i]
    end
    return Piecewise_Function(starts_, functions_)
end
function -(f1::Piecewise_Function,f2::Piecewise_Function)
    c_f1, c_f2 = create_common_pieces(f1,f2)
    starts_ = c_f1.starts_
    number_of_pieces = length(starts_)
    functions_ = Array{UnivariateFunction}(undef, number_of_pieces)
    for i in 1:number_of_pieces
        functions_[i] = c_f1.functions_[i] - c_f2.functions_[i]
    end
    return Piecewise_Function(starts_, functions_)
end


text:
# Copied from CUDA.jl/src/pool.jl

"""
    NonReentrantLock

Simple non-reentrant lock that errors when trying to reenter on the same task.
"""
struct NonReentrantLock # <: Base.AbstractLock
    rl::ReentrantLock
    NonReentrantLock() = new(ReentrantLock())
end

function Base.lock(nrl::NonReentrantLock)
    @assert !islocked(nrl.rl) || nrl.rl.locked_by !== current_task()
    lock(nrl.rl)
end

function Base.trylock(nrl::NonReentrantLock)
    @assert !islocked(nrl.rl) || nrl.rl.locked_by !== current_task()
    trylock(nrl.rl)
end

Base.unlock(nrl::NonReentrantLock) = unlock(nrl.rl)

# NonReentrantLock may be taken around code that might call the GC, which might
# reenter through finalizers.  Avoid that by temporarily disabling finalizers
# running concurrently on this thread.
enable_finalizers(on::Bool) = ccall(:jl_gc_enable_finalizers, Cvoid,
                                    (Ptr{Cvoid}, Int32,), Core.getptls(), on)
macro safe_lock(l, ex)
    quote
        temp = $(esc(l))
        lock(temp)
        enable_finalizers(false)
        try
            $(esc(ex))
        finally
            unlock(temp)
            enable_finalizers(true)
        end
    end
end

# If we actually want to acquire a lock from a finalizer, we can't cause a task
# switch. As a NonReentrantLock can only be taken by another thread that should
# be running, and not a concurrent task we'd need to switch to, we can safely
# spin.
macro safe_lock_spin(l, ex)
    quote
        temp = $(esc(l))
        while !trylock(temp)
            # we can't yield here
        end
        enable_finalizers(false) # retains compatibility with non-finalizer callers
        try
            $(esc(ex))
        finally
            unlock(temp)
            enable_finalizers(true)
        end
    end
end


text:
using DashBase
@testset "lazy + dynamic" begin
    test_resource = Resource(
        relative_package_path = "test.js",
        dynamic = true
    )
    @test test_resource.async == :lazy
    @test DashBase.isdynamic(test_resource, true)
    @test DashBase.isdynamic(test_resource, false)

    test_resource = Resource(
        relative_package_path = "test.js",
        dynamic = false
    )
    @test test_resource.async == :none
    @test !DashBase.isdynamic(test_resource, true)
    @test !DashBase.isdynamic(test_resource, false)

    @test_throws ArgumentError test_resource = Resource(
        relative_package_path = "test.js",
        dynamic = false,
        async = :false
    )

    test_resource = Resource(
        relative_package_path = "test.js",
        async = :lazy
    )

    @test test_resource.async == :lazy
    @test DashBase.isdynamic(test_resource, true)
    @test DashBase.isdynamic(test_resource, false)

    test_resource = Resource(
        relative_package_path = "test.js",
        async = :eager
    )

    @test test_resource.async == :eager
    @test !DashBase.isdynamic(test_resource, true)
    @test DashBase.isdynamic(test_resource, false)
end

@testset "register" begin
    DashBase.register_package(
        ResourcePkg(
            "dash_html_components",
            "path",
            version = "1.2.3",
            [
                Resource(
                    relative_package_path = "dash_html_components.min.js",
                    external_url = "https://unpkg.com/dash-html-components@1.0.1/dash_html_components/dash_html_components.min.js",
                    dynamic = nothing,
                    async = nothing,
                    type = :js
                ),
                Resource(
                    relative_package_path = "dash_html_components.min.js.map",
                    external_url = "https://unpkg.com/dash-html-components@1.0.1/dash_html_components/dash_html_components.min.js.map",
                    dynamic = true,
                    async = nothing,
                    type = :js
                )
            ]
        )
    )
    @test length(DashBase.main_registry().components) == 1
end

max_stars_repo_name: waralex/DashBase

max_issues_repo_name: waralex/DashBase

max_forks_repo_name: waralex/DashBase


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

@doc raw"""RawExtension is used to hold extensions in external versions.  To use this, make a field which has RawExtension as its type in your external, versioned struct, and Object in your internal struct. You also need to register your various plugin types.  // Internal package: type MyAPIObject struct {  runtime.TypeMeta `json:\",inline\"`  MyPlugin runtime.Object `json:\"myPlugin\"` } type PluginA struct {  AOption string `json:\"aOption\"` }  // External package: type MyAPIObject struct {  runtime.TypeMeta `json:\",inline\"`  MyPlugin runtime.RawExtension `json:\"myPlugin\"` } type PluginA struct {  AOption string `json:\"aOption\"` }  // On the wire, the JSON will look something like this: {  \"kind\":\"MyAPIObject\",  \"apiVersion\":\"v1\",  \"myPlugin\": {   \"kind\":\"PluginA\",   \"aOption\":\"foo\",  }, }  So what happens? Decode first uses json or yaml to unmarshal the serialized data into your external MyAPIObject. That causes the raw JSON to be stored, but not unpacked. The next step is to copy (using pkg/conversion) into the internal struct. The runtime package's DefaultScheme has conversion functions installed which will unpack the JSON stored in RawExtension, turning it into the correct object type, and storing it in the Object. (TODO: In the case where the object is of an unknown type, a runtime.Unknown object will be created and stored.)

    IoK8sApimachineryPkgRuntimeRawExtension(;
    )

"""
mutable struct IoK8sApimachineryPkgRuntimeRawExtension <: SwaggerModel

    function IoK8sApimachineryPkgRuntimeRawExtension(;)
        o = new()
        o
    end
end # type IoK8sApimachineryPkgRuntimeRawExtension

const _property_map_IoK8sApimachineryPkgRuntimeRawExtension = Dict{Symbol,Symbol}()
const _property_types_IoK8sApimachineryPkgRuntimeRawExtension = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }) = collect(keys(_property_map_IoK8sApimachineryPkgRuntimeRawExtension))
Swagger.property_type(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgRuntimeRawExtension[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgRuntimeRawExtension[property_name]

function check_required(o::IoK8sApimachineryPkgRuntimeRawExtension)
    true
end

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


text:
# common args for all hdmr tests
testargs = Dict(:verbose=>false,:showwarnings=>true)

@testset "hdmr" begin

hdmrmodels = [
    @model(h ~ x + z + cat + y, c ~ x + z + cat + y)
    @model(h ~ x + z + cat + y, c ~ z + cat + y)
    @model(h ~ x + z + cat, c ~ x + z + cat + y)
    @model(h ~ x + z + cat, c ~ z + cat + y)
    ]

@testset "$M" for M in (InclusionRepetition,Hurdle)

includel = (M == InclusionRepetition)

@testset "$(replace(string(f),"\n" => " "))" for f in hdmrmodels

# # parse and merge rhs terms
trmszero = HurdleDMR.getrhsterms(f, :h)
trmspos = HurdleDMR.getrhsterms(f, :c)
trms, inzero, inpos = HurdleDMR.mergerhsterms(trmszero,trmspos)

# create model matrix
covarsb, countsb, as = modelcols(trms, covarsdf, counts; model=M)

# inzero and inpos may be different in mm with factor variables
inzero, inpos = HurdleDMR.mapins(inzero, inpos, as)
pzero = length(inzero)
ppos = length(inpos)

# resolve projdirs
projdir = HurdleDMR.ixprojdir(as, :y)
dirpos = something(findfirst(isequal(projdir), inpos), 0)
dirzero = something(findfirst(isequal(projdir), inzero), 0)

# hurdle dmr parallel local cluster
hdmrcoefs = hdmr(covars, counts, M; inzero=inzero, inpos=inpos, parallel=true, testargs...)
coefsHppos, coefsHpzero = coef(hdmrcoefs)
@test size(coefsHppos) == (ppos+1, d)
@test size(coefsHpzero) == (pzero+1, d)
@test_throws MethodError coef(hdmrcoefs, AllSeg())

# test Int matrix for counts
hdmrcoefsint = hdmr(covars, countsint, M; inzero=inzero, inpos=inpos, parallel=true, testargs...)
coefsHpposint, coefsHpzeroint = coef(hdmrcoefsint)
@test coefsHpposint == coefsHppos
@test coefsHpzeroint == coefsHpzero

# test standardize=false
hdmrcoefsunstd = hdmr(covars, counts, M; standardize=false, inzero=inzero, inpos=inpos, parallel=true, testargs...)
coefsHpposunstd, coefsHpzerounstd = coef(hdmrcoefsunstd)
@test coefsHpposunstd != coefsHppos
@test coefsHpzerounstd != coefsHpzero
hdmrcoefsunstd2 = hdmr(covars, counts, M; standardize=false, inzero=inzero, inpos=inpos, local_cluster=false, parallel=true, testargs...)
coefsHpposunstd2, coefsHpzerounstd2 = coef(hdmrcoefsunstd2)
@test coefsHpposunstd ≈ coefsHpposunstd2
@test coefsHpzerounstd ≈ coefsHpzerounstd2

hdmrcoefsb = fit(HDMRCoefs{M}, covars, counts; inzero=inzero, inpos=inpos, parallel=true, testargs...)
@test coef(hdmrcoefsb)[1] == coefsHppos
@test coef(hdmrcoefsb)[2] == coefsHpzero
@test n == nobs(hdmrcoefs)
@test d == ncategories(hdmrcoefs)
@test ppos == ncovarspos(hdmrcoefs)
@test pzero == ncovarszero(hdmrcoefs)

# select=MinBIC()
hdmrcoefsb = fit(HDMRCoefs{M}, covars, counts; inzero=inzero, inpos=inpos, select=MinBIC(), testargs...)
@test coef(hdmrcoefsb)[1] != coefsHppos
# @test coef(hdmrcoefsb)[2] != coefsHpzero
hdmrb = fit(HDMR{M}, covars, counts; inzero=inzero, inpos=inpos, select=MinBIC(), testargs...)
@test coef(hdmrb)[1] == coef(hdmrcoefsb)[1]
@test coef(hdmrb)[2] == coef(hdmrcoefsb)[2]

# hurdle dmr parallel remote cluster
hdmrcoefs2 = hdmr(covars, counts, M; inzero=inzero, inpos=inpos, parallel=true, local_cluster=false, testargs...)
coefsHppos2, coefsHpzero2 = coef(hdmrcoefs2)
@test coefsHppos ≈ coefsHppos2
@test coefsHpzero ≈ coefsHpzero2
hdmrcoefs2 = fit(HDMRCoefs{M}, covars, counts; inzero=inzero, inpos=inpos, parallel=true, local_cluster=false, testargs...)
@test coef(hdmrcoefs2)[1] ≈ coefsHppos2
@test coef(hdmrcoefs2)[2] ≈ coefsHpzero2
@test_throws ErrorException predict(hdmrcoefs2,newcovars)

hdmrpaths3 = fit(HDMRPaths{M}, covars, counts; inzero=inzero, inpos=inpos, parallel=true, testargs...)
coefsHppos3, coefsHpzero3 = coef(hdmrpaths3)
@test coefsHppos3 ≈ coefsHppos
@test coefsHpzero3 ≈ coefsHpzero
@test n == nobs(hdmrpaths3)
@test d == ncategories(hdmrpaths3)
@test ppos == ncovarspos(hdmrpaths3)
@test pzero == ncovarszero(hdmrpaths3)
@test size(hdmrpaths3.nlpaths,1) == d
η = predict(hdmrpaths3,newcovars)
@test sum(η, dims=2) ≈ ones(size(η, 1))
coefsallpos, coefsallzero = coef(hdmrpaths3, AllSeg())
@test size(coefsallpos,1) > 1
@test size(coefsallpos,2) == ppos+1
@test size(coefsallpos,3) == d
@test size(coefsallzero,1) > 1
@test size(coefsallzero,2) == pzero+1
@test size(coefsallzero,3) == d

# hurdle dmr serial
hdmrcoefs3 = hdmr(covars, counts, M; inzero=inzero, inpos=inpos, parallel=false, testargs...)
coefsHspos, coefsHszero = coef(hdmrcoefs3)
@test coefsHppos ≈ coefsHspos
@test coefsHpzero ≈ coefsHszero
hdmrcoefs3 = fit(HDMRCoefs{M}, covars, counts; inzero=inzero, inpos=inpos, parallel=false, testargs...)
@test coef(hdmrcoefs3)[1] ≈ coefsHspos
@test coef(hdmrcoefs3)[2] ≈ coefsHszero

# serial paths
hdmrcoefs4 = hdmrpaths(covars, counts, M; inzero=inzero, inpos=inpos, parallel=false, testargs...)
coefsHspos4, coefsHszero4 = coef(hdmrcoefs4)
@test coefsHspos == coefsHspos4
@test coefsHszero == coefsHszero4

# using a dataframe and formula
hdmrcoefsdf = fit(HDMRCoefs{M}, f, covarsdf, counts; inzero=inzero, inpos=inpos, parallel=true, testargs...)
@test coef(hdmrcoefsdf)[1] == coefsHppos
@test coef(hdmrcoefsdf)[2] ≈ coefsHpzero
@test_throws ErrorException predict(hdmrcoefsdf,newcovars)

hdmrpathsdf = fit(HDMRPaths{M}, f, covarsdf, counts; inzero=inzero, inpos=inpos, parallel=true, testargs...)
@test coef(hdmrpathsdf)[1] == coefsHppos3
@test coef(hdmrpathsdf)[2] ≈ coefsHpzero3
@test predict(hdmrpathsdf,newcovars) ≈ η

if M == InclusionRepetition
    srprojcounts = counts - posindic(counts)
else
    srprojcounts = counts
end

zHpos = srproj(coefsHppos, srprojcounts)
@test size(zHpos) == (n,ppos+1)

zHzero = srproj(coefsHpzero, posindic(counts))
@test size(zHzero) == (n,pzero+1)

Z1 = srproj(coefsHppos, coefsHpzero, counts, dirpos, dirzero;
    includel=includel, intercept=true)

if dirpos > 0
    z1pos = srproj(coefsHppos, srprojcounts, dirpos)
    @test z1pos ≈ zHpos[:,[dirpos,ppos+1]]
    @test hascol(Z1,z1pos[:,1])
    @test !includel || hascol(Z1,z1pos[:,2])
end

if dirzero > 0
    z1zero = srproj(coefsHpzero, posindic(counts), dirzero)
    @test z1zero ≈ zHzero[:,[dirzero,pzero+1]]
    @test hascol(Z1,z1zero[:,1])
    @test !includel || hascol(Z1,z1zero[:,2])
end

Z1b = srproj(hdmrcoefs, counts, dirpos, dirzero; intercept=true)
@test Z1 == Z1b

@test_throws ErrorException srproj(coefsHppos, coefsHpzero, counts, 0, 0; intercept=true)

X1, X1_nocounts, inz = srprojX(coefsHppos,coefsHpzero,counts,covars,projdir;
    inzero=inzero, inpos=inpos, includel=includel, includem=true)
ix = filter!(x->x!=projdir,collect(1:p))
@test X1_nocounts ≈ [ones(n) covars[:,ix]] rtol=1e-8
@test X1 ≈ [X1_nocounts Z1] rtol=1e-8
X1b, X1_nocountsb, inzb = srprojX(hdmrcoefs,counts,covars,projdir;
    inzero=inzero, inpos=inpos, includem=true)
@test X1 ≈ X1b rtol=1e-8
@test X1_nocountsb ≈ X1_nocountsb rtol=1e-8
@test inz == inzb

X2, X2_nocounts, inz = srprojX(coefsHppos,coefsHpzero,counts,covars,projdir;
    inzero=inzero, inpos=inpos, includel=includel, includem=false)
@test X2_nocounts ≈ X1_nocounts rtol=1e-8
if includel && dirzero > 0
    @test X2 ≈ [X1[:,1:end-2] X1[:,end]] rtol=1e-8
else
    @test X2 ≈ X1[:,1:end-1] rtol=1e-8
end
X2b, X2_nocountsb, inzb = srprojX(hdmrcoefs,counts,covars,projdir;
    inzero=inzero, inpos=inpos, includem=false)
@test X2 ≈ X2b rtol=1e-8
@test X2_nocountsb ≈ X2_nocountsb rtol=1e-8
@test inz == inzb

X3, X3_nocounts, inz3 = @test_logs (:info, r"rank\(X\)") srprojX(coefsHppos,coefsHpzero,zero(counts),covars,projdir;
    inzero=inzero, inpos=inpos, includel=includel, includem=true)
@test X3_nocounts == [ones(n) covars[:,setdiff(1:p,[projdir])]]
@test inz3 == [2]

X3, X3_nocounts, inz3b = @test_logs (:info, "includezpos == false. dropping zpos.") srprojX(coefsHppos,coefsHpzero,zero(counts),covars,projdir;
    inzero=inzero, inpos=inpos, includel=includel, includem=true, inz=inz3)
@test X3_nocounts == [ones(n) covars[:,setdiff(1:p,[projdir])]]
@test inz3 == inz3b

# HIR
hir = fit(CIR{HDMR{M}, LinearModel},covars,counts,projdir;
    inzero=inzero, inpos=inpos, nocounts=true, testargs...)
@test coefbwd(hir)[1] ≈ coef(hdmrcoefs)[1]
@test coefbwd(hir)[2] ≈ coef(hdmrcoefs)[2]

hirglm = fit(CIR{HDMR{M}, GeneralizedLinearModel},covars,counts,projdir,Gamma();
    inzero=inzero, inpos=inpos, nocounts=true, testargs...)
@test coefbwd(hirglm)[1] ≈ coef(hdmrcoefs)[1]
@test coefbwd(hirglm)[2] ≈ coef(hdmrcoefs)[2]
@test !(coeffwd(hirglm)[1] ≈ coeffwd(hir)[1])
@test !(coeffwd(hirglm)[2] ≈ coeffwd(hir)[2])

hirdf = fit(CIR{HDMR{M}, LinearModel},f,covarsdf,counts,:y;
    inzero=inzero, inpos=inpos, nocounts=true, testargs...)
@test coefbwd(hirdf)[1] ≈ coef(hdmrcoefs)[1]
@test coefbwd(hirdf)[2] ≈ coef(hdmrcoefs)[2]
@test coeffwd(hirdf) ≈ coeffwd(hir)
hirglmdf = fit(CIR{HDMR{M}, GeneralizedLinearModel},f,covarsdf,counts,:y,Gamma();
    inzero=inzero, inpos=inpos, nocounts=true, testargs...)
@test coefbwd(hirglmdf)[1] ≈ coef(hdmrcoefs)[1]
@test coefbwd(hirglmdf)[2] ≈ coef(hdmrcoefs)[2]
@test !(coeffwd(hirglmdf)[2] ≈ coeffwd(hirdf)[2])
@test !(coeffwd(hirglmdf)[1] ≈ coeffwd(hirdf)[1])

# select=MinBIC()
hirdfb = fit(CIR{HDMR{M}, LinearModel},f,covarsdf,counts,:y;
    inzero=inzero, inpos=inpos, nocounts=true, select=MinBIC(), testargs...)
@test coefbwd(hirdf)[1] != coefbwd(hirdfb)[1]
# @test coefbwd(hirdf)[2] != coefbwd(hirdfb)[2]
@test coeffwd(hirdf)[1] != coeffwd(hirdfb)[1]
@test coeffwd(hirdf)[2] != coeffwd(hirdfb)[2]

zlm = lm(hcat(ones(n,1),Z1,covars[:,1:4]),covars[:,projdir])
@test r2(zlm) ≈ r2(hir)
@test adjr2(zlm) ≈ adjr2(hir)
@test predict(zlm,hcat(ones(10,1),Z1[1:10,:],covars[1:10,1:4])) ≈ predict(hir,covars[1:10,:],counts[1:10,:])

zlmnocounts = lm(hcat(ones(n,1),covars[:,1:4]),covars[:,projdir])
@test r2(zlmnocounts) ≈ r2(hir; nocounts=true)
@test adjr2(zlmnocounts) ≈ adjr2(hir; nocounts=true)
@test predict(zlmnocounts,hcat(ones(10,1),covars[1:10,1:4])) ≈ predict(hir,covars[1:10,:],counts[1:10,:]; nocounts=true)

end

@testset "HDMR{$M} degenerate cases" begin

@info("Testing hdmr degenerate cases. The 2 following warnings by workers are expected ...")

# column j is never zero, so hj=1 for all observations
zcounts = deepcopy(counts)
zcounts[:,2] = zeros(n)
zcounts[:,3] = ones(n)

# make sure we are not adding all zero obseravtions
m = sum(zcounts, dims=2)
@test sum(m .== 0) == 0

# hurdle dmr parallel local cluster
hdmrcoefs = fit(HDMR{M}, covars, zcounts; parallel=true, testargs...)
coefsHppos, coefsHpzero = coef(hdmrcoefs)
@test size(coefsHppos) == (p+1, d)
@test size(coefsHpzero) == (p+1, d)
@test coefsHppos[:,2] == zeros(p+1)
@test coefsHpzero[:,2] == zeros(p+1)
@test coefsHppos[:,3] == zeros(p+1)
@test coefsHpzero[:,3] == zeros(p+1)
@test HurdleDMR.includelinX(hdmrcoefs) == includel

if M == Hurdle
    η3 = 0.36
    η4 = 0.6
else
    η3 = 0.32
    η4 = 0.64
end

# hurdle dmr parallel remote cluster
hdmrcoefs2 = fit(HDMRPaths{M},covars, zcounts; parallel=true, local_cluster=false, testargs...)
coefsHppos2, coefsHpzero2 = coef(hdmrcoefs)
@test coefsHppos ≈ coefsHppos2
@test coefsHpzero ≈ coefsHpzero2
η = predict(hdmrcoefs2,newcovars)
@test sum(η, dims=2) ≈ ones(size(newcovars, 1))
@test η[:,2] == zeros(size(newcovars,1))
@test η[:,3] ≈ ones(size(newcovars,1))*η3 rtol=0.05
@test η[:,4] ≈ ones(size(newcovars,1))*η4 rtol=0.06
@test HurdleDMR.includelinX(hdmrcoefs2) == includel

# hurdle dmr serial paths
rx = Regex("fit\\($M...\\) failed for countsj")
hdmrcoefs3 = @test_logs (:warn, rx) fit(HDMRPaths{M},covars, zcounts; parallel=false, testargs...)
coefsHppos3, coefsHpzero3 = coef(hdmrcoefs3)
@test coefsHppos ≈ coefsHppos3
@test coefsHpzero ≈ coefsHpzero3
η = predict(hdmrcoefs3,newcovars)
@test sum(η, dims=2) ≈ ones(size(newcovars, 1))
@test η[:,2] == zeros(size(newcovars,1))
@test η[:,3] ≈ ones(size(newcovars,1))*η3 rtol=0.05
@test η[:,4] ≈ ones(size(newcovars,1))*η4 rtol=0.06
@test HurdleDMR.includelinX(hdmrcoefs3) == includel

# hurdle dmr serial coefs
hdmrcoefs4 = @test_logs (:warn, r"failed on count dimension 2") fit(HDMR{M},covars, zcounts; parallel=false, testargs...)
coefsHppos4, coefsHpzero4 = coef(hdmrcoefs4)
@test size(coefsHppos4) == (p+1, d)
@test size(coefsHpzero4) == (p+1, d)
@test coefsHppos4[:,2] == zeros(p+1)
@test coefsHpzero4[:,2] == zeros(p+1)
@test coefsHppos4[:,3] == zeros(p+1)
@test coefsHpzero4[:,3] == zeros(p+1)
@test HurdleDMR.includelinX(hdmrcoefs4) == includel

if M == InclusionRepetition
    @info("Testing hdmr degenerate cases. The 12 following warnings by workers are expected ...")
end

zcounts = Matrix(deepcopy(counts))
Random.seed!(13)
for I = eachindex(zcounts)
    if iszero(zcounts[I])
        zcounts[I] = rand(1:10)
    end
end

# make sure we are not adding all zero obseravtions
m = sum(zcounts, dims=2)
@test sum(m .== 0) == 0
if M == InclusionRepetition
    l = sum(posindic(zcounts), dims=2)
    @test sum((m .- l) .== 0) == 0
    @test sum((d .- l) .== 0) == n
end

# hurdle dmr parallel local cluster
hdmrcoefs = fit(HDMR{M}, covars, zcounts; parallel=true, testargs...)
coefsHppos, coefsHpzero = coef(hdmrcoefs)
@test size(coefsHppos) == (p+1, d)
@test size(coefsHpzero) == (p+1, d)
if M == Hurdle
	@test coefsHppos[:,2] != zeros(p+1)
	@test coefsHppos[:,3] != zeros(p+1)
else
	@test coefsHpzero[:,2] == zeros(p+1)
	@test coefsHpzero[:,3] == zeros(p+1)
end
@test HurdleDMR.includelinX(hdmrcoefs) == includel

# hurdle dmr parallel remote cluster
hdmrcoefs2 = fit(HDMRPaths{M}, covars, zcounts; parallel=true, local_cluster=false, testargs...)
coefsHppos2, coefsHpzero2 = coef(hdmrcoefs2)
@test coefsHppos ≈ coefsHppos2
@test coefsHpzero ≈ coefsHpzero2
@test HurdleDMR.includelinX(hdmrcoefs2) == includel

# just checking the fit(HDMRPaths...) ignore local_cluster
hdmrcoefs3 = fit(HDMRPaths{M},covars, zcounts; parallel=true, local_cluster=true, testargs...)
coefsHppos3, coefsHpzero3 = coef(hdmrcoefs3)
@test coefsHppos2 == coefsHppos3
@test coefsHpzero2 == coefsHpzero3
@test HurdleDMR.includelinX(hdmrcoefs3) == includel

# hurdle dmr parallel remote cluster
if M == Hurdle
    hdmrcoefs4 = fit(HDMRPaths{M},covars, zcounts; parallel=false, testargs...)
else
    rx = Regex("fit\\($M...\\) failed for countsj")
    warnings = [(:warn, rx) for i=1:4]
    hdmrcoefs4 = @test_logs(warnings..., fit(HDMRPaths{M},covars, zcounts; parallel=false, testargs...))
end
coefsHppos4, coefsHpzero4 = coef(hdmrcoefs4)
@test coefsHppos ≈ coefsHppos4
@test coefsHpzero ≈ coefsHpzero4
@test HurdleDMR.includelinX(hdmrcoefs4) == includel

end

end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2df34aa4-a3c7-11ec-3508-61f14f475394
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("PlutoUI")
		Pkg.add("Statistics")
		Pkg.add("StatsBase")
		Pkg.add("ImageMorphology")
		Pkg.add("ImageFiltering")
		Pkg.add("CSV")
		Pkg.add("DataFrames")
		Pkg.add("GLM")
		Pkg.add(url="https://github.com/JuliaHealth/DICOM.jl")
		Pkg.add(url="https://github.com/Dale-Black/DICOMUtils.jl")
		Pkg.add(url="https://github.com/Dale-Black/Phantoms.jl")
		Pkg.add(url="https://github.com/Dale-Black/CalciumScoring.jl")
		Pkg.add("CairoMakie")
	end

	using PlutoUI
	using Statistics
	using StatsBase: quantile!
	using ImageMorphology
	using ImageFiltering
	using CSV
	using DataFrames
	using GLM
	using DICOM
	using DICOMUtils
	using Phantoms
	using CalciumScoring
	using CairoMakie
end

# ╔═╡ 61d478bf-7bde-4184-9231-fd4a8bb1e6a1
begin
	VENDER = "GE_Revolution"
	BASE_PATH = "/Users/daleblack/Google Drive/Datasets/"
end

# ╔═╡ 64eb1515-c69e-41f0-8fdb-f93389a7003f
md"""
# Save results
"""

# ╔═╡ f5509a75-4a0d-484d-9b9d-7fbe72876cfe
if ~isdir(string(cd(pwd, "..") , "/data/output/", VENDER, "_agatston_mass_score_script"))
	mkdir(string(cd(pwd, "..") , "/data/output/", VENDER, "_agatston_mass_score_script",))
end

# ╔═╡ 09b33976-a24e-4bd7-b8a6-7fc47efb6f94
output_path = string(cd(pwd, "..") , "/data/output/", VENDER, "_agatston_mass_score_script", "/output.csv")

# ╔═╡ a1c582c1-85c5-498f-9ecb-e70cd7819715
function create_mask(array, mask)
	@assert size(array) == size(mask)
	idxs = findall(x -> x == true, mask)
	overlayed_mask = zeros(size(array))
	for idx in idxs
		overlayed_mask[idx] = array[idx]
	end
	return overlayed_mask
end

# ╔═╡ b8afde9f-a02a-4165-908e-41afe98109e8
begin
	dfs = []
	for s in 1:10
		SCAN_NUMBER = s
		root_path = string(BASE_PATH, VENDER)
		dcm_path_list = dcm_list_builder(root_path)
		pth = dcm_path_list[SCAN_NUMBER]
		scan = basename(pth)
		header, dcm_array, slice_thick_ori1 = dcm_reader(pth)

		# Segment Heart
		masked_array, center_insert, mask = mask_heart(header, dcm_array, size(dcm_array, 3)÷2)

		# Segment Calcium Rod
		calcium_image, slice_CCI, quality_slice, cal_rod_slice = mask_rod(masked_array, header)

		# Segment Calcium Inserts
		mask_L_HD, mask_M_HD, mask_S_HD, mask_L_MD, mask_M_MD, mask_S_MD, mask_L_LD, mask_M_LD, mask_S_LD = mask_inserts(
	            dcm_array, masked_array, header, slice_CCI, center_insert
		)

		# Mask Calibration Factor
		output = calc_output(masked_array, header, slice_CCI, 130, trues(3, 3))
	    insert_centers = calc_centers(dcm_array, output, header, center_insert, slice_CCI)
		rows, cols = Int(header[tag"Rows"]), Int(header[tag"Columns"])
		pixel_size = DICOMUtils.get_pixel_size(header)
		mass_cal_factor, angle_0_200HA, water_rod_metrics = mass_calibration(masked_array, insert_centers[:Large_LD], center_insert, cal_rod_slice, rows, cols, pixel_size)

		# Score Large Inserts
		arr = masked_array[:, :, slice_CCI-2:slice_CCI+2]
		single_arr = masked_array[:, :, slice_CCI]

		## High Density
		mask_L_HD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_L_HD_3D[:, :, z] = mask_L_HD
		end
		dilated_mask_L_HD = dilate(dilate(mask_L_HD_3D))
		alg = Agatston()
		overlayed_mask_l_hd = create_mask(arr, dilated_mask_L_HD)
		agat_l_hd, mass_l_hd = score(overlayed_mask_l_hd, pixel_size, mass_cal_factor, alg)

		## Medium Density
		mask_L_MD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_L_MD_3D[:, :, z] = mask_L_MD
		end
		dilated_mask_L_MD = dilate(dilate(mask_L_MD_3D))
		overlayed_mask_l_md = create_mask(arr, dilated_mask_L_MD)
		agat_l_md, mass_l_md = score(overlayed_mask_l_md, pixel_size, mass_cal_factor, alg)

		## Low Density
		mask_L_LD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_L_LD_3D[:, :, z] = mask_L_LD
		end
		dilated_mask_L_LD = dilate(dilate(mask_L_LD_3D))
		overlayed_mask_l_ld = create_mask(arr, dilated_mask_L_LD)
		agat_l_ld, mass_l_ld = score(overlayed_mask_l_ld, pixel_size, mass_cal_factor, alg)

		# Score Medium Inserts
		## High Density
		mask_M_HD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_M_HD_3D[:, :, z] = mask_M_HD
		end
		dilated_mask_M_HD = dilate(dilate(dilate(dilate(mask_M_HD_3D))))
		overlayed_mask_m_hd = create_mask(arr, dilated_mask_M_HD)
		agat_m_hd, mass_m_hd = score(overlayed_mask_m_hd, pixel_size, mass_cal_factor, alg)

		## Medium Density
		mask_M_MD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_M_MD_3D[:, :, z] = mask_M_MD
		end
		dilated_mask_M_MD = dilate(dilate(dilate(dilate(mask_M_MD_3D))))
		overlayed_mask_m_md = create_mask(arr, dilated_mask_M_MD)
		agat_m_md, mass_m_md = score(overlayed_mask_m_md, pixel_size, mass_cal_factor, alg)

		## Low Density
		mask_M_LD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_M_LD_3D[:, :, z] = mask_M_LD
		end
		dilated_mask_M_LD = dilate(dilate(dilate(dilate(mask_M_LD_3D))))
		overlayed_mask_m_ld = create_mask(arr, dilated_mask_M_LD)
		agat_m_ld, mass_m_ld = score(overlayed_mask_m_ld, pixel_size, mass_cal_factor, alg)

		# Score Small Inserts
		## High Density
		mask_S_HD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_S_HD_3D[:, :, z] = mask_S_HD
		end
		dilated_mask_S_HD = dilate(dilate(dilate(dilate(dilate((mask_S_HD_3D))))))
		overlayed_mask_s_hd = create_mask(arr, dilated_mask_S_HD)
		agat_s_hd, mass_s_hd = score(overlayed_mask_s_hd, pixel_size, mass_cal_factor, alg)

		## Medium Density
		mask_S_MD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_S_MD_3D[:, :, z] = mask_S_MD
		end
		dilated_mask_S_MD = dilate(dilate(dilate(dilate(dilate(mask_S_MD_3D)))))
		overlayed_mask_s_md = create_mask(arr, dilated_mask_S_MD)
		agat_s_md, mass_s_md = score(overlayed_mask_s_md, pixel_size, mass_cal_factor, alg)

		## Low Density
		mask_S_LD_3D = Array{Bool}(undef, size(arr))
		for z in 1:size(arr, 3)
			mask_S_LD_3D[:, :, z] = mask_S_LD
		end
		dilated_mask_S_LD = dilate(dilate(dilate(dilate(dilate(mask_S_LD_3D)))))
		overlayed_mask_s_ld = create_mask(arr, dilated_mask_S_LD)
		agat_s_ld, mass_s_ld = score(overlayed_mask_s_ld, pixel_size, mass_cal_factor, alg)

		# Results
		density_array = [0, 200, 400, 800]
		inserts = [
			"Low Density",
			"Medium Density",
			"High Density"
		]

		## Agatston
		calculated_agat_large = [
			agat_l_ld,
			agat_l_md,
			agat_l_hd
		]
		calculated_agat_medium = [
			agat_m_ld,
			agat_m_md,
			agat_m_hd
		]
		calculated_agat_small = [
			agat_s_ld,
			agat_s_md,
			agat_s_hd
		]

		## Mass
		ground_truth_mass_large = [
			19.6,
			39.3,
			78.5
		] # mg
		calculated_mass_large = [
			mass_l_ld,
			mass_l_md,
			mass_l_hd
		]
		ground_truth_mass_medium = [
			4.2,
			8.5,
			17.0
		]
		calculated_mass_medium = [
			mass_m_ld,
			mass_m_md,
			mass_m_hd
		]
		ground_truth_mass_small = [
			0.2,
			0.3,
			0.6
		]
		calculated_mass_small = [
			mass_s_ld,
			mass_s_md,
			mass_s_hd
		]

		df = DataFrame(
			scan = scan,
			inserts = inserts,
			calculated_agat_large = calculated_agat_large,
			calculated_agat_medium = calculated_agat_medium,
			calculated_agat_small = calculated_agat_small,
			ground_truth_mass_large = ground_truth_mass_large,
			calculated_mass_large = calculated_mass_large,
			ground_truth_mass_medium = ground_truth_mass_medium,
			calculated_mass_medium = calculated_mass_medium,
			ground_truth_mass_small = ground_truth_mass_small,
			calculated_mass_small = calculated_mass_small,
			mass_cal_factor = mass_cal_factor
		)
		push!(dfs, df)
	end
end

# ╔═╡ 94bd4f44-cd03-4602-8dd5-b3f0ba99e2ed
new_df = vcat(dfs[1:10]...)

# ╔═╡ ef536a67-d206-4c42-ae09-bfc8c9e1abfc
CSV.write(output_path, new_df)

# ╔═╡ Cell order:
# ╠═2df34aa4-a3c7-11ec-3508-61f14f475394
# ╠═61d478bf-7bde-4184-9231-fd4a8bb1e6a1
# ╠═b8afde9f-a02a-4165-908e-41afe98109e8
# ╟─64eb1515-c69e-41f0-8fdb-f93389a7003f
# ╠═94bd4f44-cd03-4602-8dd5-b3f0ba99e2ed
# ╠═f5509a75-4a0d-484d-9b9d-7fbe72876cfe
# ╠═09b33976-a24e-4bd7-b8a6-7fc47efb6f94
# ╠═ef536a67-d206-4c42-ae09-bfc8c9e1abfc
# ╟─a1c582c1-85c5-498f-9ecb-e70cd7819715

max_stars_repo_name: MolloiLab/CAC-stanford-data

max_issues_repo_name: MolloiLab/CAC-stanford-data

max_forks_repo_name: MolloiLab/CAC-stanford-data


text:
struct MemoryLocation{D}
    address         ::Address{D}
    other_addresses ::Vector{Address{D}}
    offsets         ::Vector{Address{D}}
end

function MemoryLocation(a::T, b::Vector{T}, c::Vector{T}) where T <: NTuple{D,<:Integer} where D
    return MemoryLocation{D}(Address(a), Address.(b), Address.(c))
end

function MemoryLocation(a::T, b::T, c::Vector{T}) where T <: NTuple{D,<:Integer} where D
    return MemoryLocation{D}(Address(a), [Address(b)], Address.(c))
end

function MemoryLocation(a::T, c::Vector{T}) where T <: NTuple{D,<:Integer} where D
    return MemoryLocation{D}(Address(a), Address{D}[], Address.(c))
end

function mem_addresses(v::Vector{MemoryLocation{D}}) where D
    a = Set{Address{D}}()
    for i in v
        push!(a, i.address)
        for j in i.other_addresses
            push!(a, j)
        end
    end
    return a
end

mem_locations(v::Vector{MemoryLocation{D}}) where D = [i.address for i in v]

function mem_neighbors(v::Vector{MemoryLocation{D}}) where D
    a = Set{Address{D}}()
    for i in v
        base = i.address
        for o in i.offsets
            push!(a, base + o)
        end
    end
    return a
end

function mem_neighbor(v::Vector{MemoryLocation{D}}) where D
    a = Array{Address{D},1}()
    for i in v
        base = i.address
        for o in i.offsets
            push!(a, base + o)
        end
    end
    return a
end

function asap_cluster(num_links, A)
    arch = TopLevel{A,2}("asap_cluster")

    # Extra parameters
    num_fifos = 2
    memories = [
                # Memory Cluster 1
                MemoryLocation((8, 8),   (8, 9)  , [(-1,0),(-1,1)]),
                MemoryLocation((8, 10),  (8, 11) , [(-1,0),(-1,1)]),
                MemoryLocation((9, 8),   (9, 9)  , [(1,0),(1,1)]),
                MemoryLocation((9, 10),  (9, 11) , [(1,0),(1,1)]),
                # Memory Cluster 2
                MemoryLocation((8, 24),  (8, 25) , [(-1,0),(-1,1)]),
                MemoryLocation((8, 26),  (8, 27) , [(-1,0),(-1,1)]),
                MemoryLocation((9, 24),  (9, 25) , [(1,0),(1,1)]),
                MemoryLocation((9, 26),  (9, 27) , [(1,0),(1,1)]),
                # Memory Cluster 3
                MemoryLocation((24, 8),  (24, 9) , [(-1,0),(-1,1)]),
                MemoryLocation((24, 10), (24, 11), [(-1,0),(-1,1)]),
                MemoryLocation((25, 8),  (25, 9) , [(1,0),(1,1)]),
                MemoryLocation((25, 10), (25, 11), [(1,0),(1,1)]),
                # Memory Cluster 4
                MemoryLocation((24, 24), (24, 25), [(-1,0),(-1,1)]),
                MemoryLocation((24, 26), (24, 27), [(-1,0),(-1,1)]),
                MemoryLocation((25, 24), (25, 25), [(1,0),(1,1)]),
                MemoryLocation((25, 26), (25, 27), [(1,0),(1,1)]),
        ]

    mem_addrs = mem_addresses(memories)
    mems = mem_locations(memories)
    mem_neighbor_addrs = mem_neighbors(memories)

    ####################
    # Normal Processor #
    ####################
    processor = build_processor_tile(num_links)
    for r in 2:33, c in 2:33
        addr = Address(r,c)
        # Skip if this is an address occupied by a memory or memory neighbor.
        if in(addr, mem_addrs) || in(addr, mem_neighbor_addrs)
            continue
        end
        add_child(arch, processor, addr)
    end

    ####################
    # Memory Processor #
    ####################
    # Instantiate a memory processor at every address neighboring a memory.
    memory_processor = build_processor_tile(num_links, include_memory = true)
    for mem_neighbor_addr in mem_neighbor_addrs
        add_child(arch, memory_processor, mem_neighbor_addr)
    end

    #################
    # 2 Port Memory #
    #################
    memory = build_memory(2)
    for mem_addr in mems
        add_child(arch, memory, mem_addr)
    end

    #################
    # Input Handler #
    #################
    input_handler = build_input_handler(1)
    for (r,c) ∈ zip((2,13,15,19), (1, 34, 1, 34))
        add_child(arch, input_handler, Address(r,c))
    end

    ##################
    # Output Handler #
    ##################
    output_handler = build_output_handler(1)
    for (r,c) ∈ zip((2,13,15,19), (34, 1, 34, 1))
        add_child(arch, output_handler, Address(r,c))
    end

    #######################
    # Global Interconnect #
    #######################
    connect_processors(arch, num_links)
    connect_io(arch, num_links)
    connect_memories_cluster(arch, memories)
    return arch
end

function connect_memories_cluster(arch, memories)
    # Create metadata dictionary for the memory links.
    request_metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "memory_request_link",
   )

    response_metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "memory_response_link",
   )

    ###########################
    # Connect 1 port memories #
    ###########################
    proc_rule = x -> search_metadata!(x, "attributes", "memory_processor", in)
    mem_rule = x -> search_metadata!(x, "attributes", "memory_2port", in)

    for mem in memories
        address = mem.address
        for (i, offset) in enumerate(mem.offsets)
            # Build Request Link.
            offset_rule = [(-offset, "memory_out", "in[$(i-1)]")]
            connection_rule(arch, offset_rule, proc_rule, mem_rule,
                            metadata = request_metadata,
                            valid_addresses = (address + offset,))

            # Build response link.
            offset_rule = [(offset, "out[$(i-1)]", "memory_in")]
            connection_rule(arch, offset_rule, mem_rule, proc_rule,
                            metadata = response_metadata,
                            valid_addresses = (address,))
        end
    end

    return nothing
end


text:
using Turing: Tracker

function turing_inference(
    prob::DiffEqBase.DEProblem,
    alg,
    t,
    data,
    priors;
    likelihood_dist_priors = [InverseGamma(2, 3)],
    likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
    num_samples=1000, num_chains=1,
    sampler = Turing.NUTS(0.65),
    syms = [Turing.@varname(theta[i]) for i in 1:length(priors)],
    sample_u0 = false,
    save_idxs = nothing,
    progress = false,
    kwargs...,
)
    N = length(priors)
    Turing.@model function mf(x, ::Type{T} = Float64) where {T <: Real}
        theta = Vector{T}(undef, length(priors))
        for i in 1:length(priors)
            theta[i] ~ NamedDist(priors[i], syms[i])
        end
        σ = Vector{T}(undef, length(likelihood_dist_priors))
        for i in 1:length(likelihood_dist_priors)
            σ[i] ~ likelihood_dist_priors[i]
        end
        nu = save_idxs === nothing ? length(prob.u0) : length(save_idxs)
        u0 = convert.(T, sample_u0 ? theta[1:nu] : prob.u0)
        p = convert.(T, sample_u0 ? theta[(nu + 1):end] : theta)
        if length(u0) < length(prob.u0)
            # assumes u is ordered such that the observed variables are in the begining, consistent with ordered theta
            for i in length(u0):length(prob.u0)
                push!(u0, convert(T,prob.u0[i]))
            end
        end
        _saveat = t === nothing ? Float64[] : t
        sol = solve(prob, alg; u0=u0, p=p, saveat = _saveat, progress = progress, save_idxs = save_idxs, kwargs...)
        failure = size(sol, 2) < length(_saveat)

        if failure
            Turing.acclogp!(_varinfo, -Inf)
            return
        end
        if ndims(sol) == 1
            x ~ likelihood(Array(sol), theta, Inf, σ)
        else
            for i = 1:length(t)
                x[:, i] ~ likelihood(sol[:, i], theta, sol.t[i], σ)
            end
        end
        return
    end false

    # Instantiate a Model object.
    model = mf(data)
    chn = sample(model, sampler, MCMCThreads(), num_samples, num_chains; progress = progress)
    return chn
end


text:
module Day11

export day11

function parseinput(filename)
    grid = Vector{Char}()
    line = readline(filename)
    xsize = length(line)
    for line in readlines(filename)
        push!(grid, collect(line)...)
    end
    reshape(grid,(xsize, div(length(grid), xsize)))
end

function adjacent_neighbours(grid, x, y)
    count = 0
    for dx in -1:1
        for dy in -1:1
            a = x+dx
            b = y+dy
            if 1 <= a <= size(grid,1) && 1 <= b <= size(grid,2) && grid[a,b] == '#'
                count += 1
            end
        end
    end
    grid[x,y] == '#' ? count-1 : count
end

function occupied_in_direction(grid, x, y, dx, dy)
    x += dx
    y += dy
    while 1 <= x <= size(grid,1) && 1 <= y <= size(grid,2)
        if grid[x,y] == 'L'
            return false
        elseif grid[x,y] == '#'
            return true
        end
        x += dx
        y += dy
    end
    false
end

function seen_neighbours(grid, x, y)
    count = 0
    for dx in -1:1
        for dy in -1:1
            if dx == 0 && dy == 0
                continue
            end
            if occupied_in_direction(grid, x, y, dx, dy)
                count += 1
            end
        end
    end
    count
end

function next(grid::Array{Char, 2}, neighbours, threshold)
    newgrid = copy(grid)
    changes = false
    for x in axes(grid,1)
        for y in axes(grid,2)
            n = neighbours(grid, x, y)
            state = grid[x,y]
            if state == 'L' && n == 0
                newgrid[x,y] = '#'
                changes = true
            elseif state == '#' && n >= threshold
                newgrid[x,y] = 'L'
                changes = true
            end
        end
    end
    (newgrid, changes)
end

function day11(::Val{:a}, filename)
    grid = parseinput(filename)
    changes = true
    while changes
        (grid, changes) = next(grid, adjacent_neighbours, 4)
    end
    count(grid) do ch
        ch == '#'
    end
end

function day11(::Val{:b}, filename)
    grid = parseinput(filename)
    changes = true
    while changes
        (grid, changes) = next(grid, seen_neighbours, 5)
    end
    count(grid) do ch
        ch == '#'
    end
end

end


text:
using ZygoteRules, Test
using ZygoteRules: typeless, named

foo(x) = 2x
@adjoint foo(x) = 3x, ȳ -> (3,)

# using Zygote

# @test gradient(foo, x) == (3,)

"""
Extract the first parameter declaration from a function definition

This macro is necessary because the same syntax is parsed differently if it is part of a
parameter list, see the first test.
"""
macro param(def)
  @assert def.head == :function
  head = def.args[1]
  call = head.head == :where ? head.args[1] : head
  QuoteNode(call.args[2])
end

@testset "Test utilities" begin
  @test @param(function f(a::Int = 5) end) != :(a::Int = 5)

  @test @param(function f(x) end) == :(x)
  @test @param(function f(a::Int = 5) end) == Expr(:kw, :(a::Int), 5)
end

@testset "Macro utilities" begin
  @testset "extracting parameter names" begin
    @test typeless(@param function f(x) end) == :x
    @test typeless(@param function f(b::Float64 = π) end) == :b
  end

  @testset "naming anonymous parameters" begin
    @test typeless(named(@param function f(w) end)) == :w
    @test typeless(named(@param function f(σ::Float32 = 1.0) end)) == :σ

    @test typeless(named(@param function f(::Int) end)) isa Symbol
    @test typeless(named(@param function f(::Type{Val{x}}) where {x} end)) isa Symbol
  end
end


text:
"""
   IO functions for the HARPS spectrograph
    https://www.eso.org/sci/facilities/lasilla/instruments/harps.html
Author: Alex Wise and collaborators
Created: April 2021
"""

"""
Create Dataframe containing filenames and key data for all folders YYYY-MM-DDThh:mm:ss.mss (year-month-day T hour:minute:second:milisecond in directory
We assume each YYYY-MM-DDThh:mm:ss.mss folder contains an e2ds.fits spectrum, which is the format for HARPS spectra from the ESO archive.
"""
function make_manifest(data_path::String ; max_spectra_to_use::Int = 1000 )
    df_foldernames = EchelleInstruments.make_manifest(data_path,r"^20\d{2}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}.\d{3}$")

    #assume each HARPS spectrum is named e2ds.fits, as they come formatted from the ESO archive.
    df_files = DataFrame(read_metadata(joinpath(df_foldernames.Filename[1],"e2ds.fits")))
    keys = propertynames(df_files)
    allowmissing!(df_files, keys[map(k->k∉[:Filename, :bjd, :target, :airmass],keys)] )

    if length(df_foldernames.Filename) >= 2
        map(fn->add_metadata_from_fits!(df_files,joinpath(fn,"e2ds.fits")),df_foldernames.Filename[2:end])
    end

    #load binaryRV corrections
    if all(df_files[:,:target] .== df_files[1,:target]) && isfile(joinpath(data_path,string(df_files[1,:target],"_binaryRV.csv")))
        @info "HARPS.make_manifest found a binary RV correction file. Loading binary RVs."
        binaryRV = CSV.read(joinpath(data_path,string(df_files[1,:target],"_binaryRV.csv")), DataFrame, header=["binaryRV"])
        df_files = hcat(df_files,binaryRV) #TODO: edit binaryRV.csv to contain timestamps, check that the times match
    end

    df_files

end

"""Create Dict containing filename and default metadata from file."""
function read_metadata(fn::Union{String,FITS})
    fields_to_save = metadata_symbols_default(HARPS2D())
    fields_str_to_save = metadata_strings_default(HARPS2D())
    dict = read_metadata_from_fits(fn,fields=fields_to_save,fields_str=fields_str_to_save)
    check_metadata_fields_expected_present(dict,fields_to_save)
    dict[:Filename] = fn
    if haskey(dict,:airmass) && typeof(dict[:airmass]) == String
       dict[:airmass] = parse(Float64,dict[:airmass])
    end
    dict[:wfile] = joinpath(fn[1:end-33],"waves",dict[:wfile])
    return dict
end

function add_metadata_from_fits!(df::DataFrame, fn::String)
    println("# Reading metadata from ", fn)
    metadata_keep = read_metadata(fn)
    #reject file from df if it is a calibration (only keep spectra)
    #if metadata_keep[:target] != "Solar" return df end
    #calibration_target_substrings = ["Etalon","LFC","ThArS","LDLS","FlatBB"]
    #if any(map(ss->contains( metadata_keep[:target],ss),calibration_target_substrings))
    #   return df
    #end

    for k in propertynames(df)
        if typeof(metadata_keep[k]) == Missing continue end
        if typeof(metadata_keep[k]) == Nothing
            metadata_keep[k] = missing
            continue
        end
        if (eltype(df[!,k]) <: Union{Real,Union{<:Real,Missing}}) && !(typeof(metadata_keep[k]) <: Union{Real,Union{<:Real,Missing}})
            metadata_keep[k] = missing
        end
    end
    #println("metadata_keep = ", metadata_keep)
    if haskey(metadata_keep,:airmass) && typeof(metadata_keep[:airmass]) == String
       metadata_keep[:airmass] = parse(Float64,metadata_keep[:airmass])
    end

    metadata_keep[:wfile] = joinpath(fn[1:end-33],"waves",metadata_keep[:wfile])

    push!(df, metadata_keep)
    return df
end

""" Read HARPS wavelength file.
    Note:HARPS wavelengths are stored in Float-32 format, which reduces wavelength precision.
    This function restores that precision using a smooth polynomial fit """
function read_λ(fn::String)
    λ_fits = FITS(fn)
    λ_original = convert(Array{Float64,2},FITSIO.read(λ_fits[1]))
    nx, norder = size(λ_original)
    λ_new = zeros(nx,norder)
    for i in 1:norder
        waves_fit = Polynomials.fit(1:nx,λ_original[:,i],2)
        λ_new[:,i] = waves_fit.(1:nx)
    end
    return λ_new
end

function read_λ(fn::String, orders_to_read::AR ) where  { AR<:AbstractRange }
    λ_fits = FITS(fn)
    λ_original = convert(Array{Float64,2},FITSIO.read(λ_fits[1],:,orders_to_read))
    nx, norder = size(λ_original)
    λ_new = zeros(nx,norder)
    for i in 1:norder
        waves_fit = Polynomials.fit(1:nx,λ_original[:,i],2)
        λ_new[:,i] = waves_fit.(1:nx)
    end
    return λ_new
end

""" Read HARPS data from FITS file, and return in a Spectra2DBasic object."""
function read_data   end

function read_data(f::FITS, metadata::Dict{Symbol,Any} )
    @assert length(f) >= 1
    #@assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    λ = read_λ(metadata[:wfile])
    flux = FITSIO.read(f[1])
    # std = NaNMath.sqrt.(flux)
    var = deepcopy(flux)
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, HARPS2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(f::FITS, metadata::Dict{Symbol,Any}, orders_to_read::AR ) where  { AR<:AbstractRange }
    @assert length(f) >= 1
    #@assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    λ = read_λ(metadata[:wfile], orders_to_read)
    flux = FITSIO.read(f[1],:,orders_to_read)
    # std = NaNMath.sqrt.(flux)
    var = deepcopy(flux)
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, HARPS2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    read_data(f, metadata)
end

function read_data(fn::String)
    f = FITS(fn)
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    read_data(f, metadata)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index})
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata)
end

function read_data(fn::String, orders_to_read::AR ) where  { AR<:AbstractRange }
    f = FITS(fn)
    metadata = read_metadata(f)
    read_data(f,metadata, orders_to_read)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index}, orders_to_read::AR ) where  { AR<:AbstractRange }
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata, orders_to_read)
end

function apply_barycentric_correction!(λ::AbstractArray{T,1}, z::Real) where { T<:Real }

end

""" Read normalization function for approximate continuum+blaze correction
TODO: Figure out mapping between these orders and what's stored in the FITS file before using. """
function read_normalization(fn::String = joinpath(pkgdir(EchelleInstruments),"data","neid","HARPS_normalization_function.csv") )
    CSV.read(fn, DataFrame)
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using Revise
using ZonalFlow

tspan   = (0.0,1000.0);
tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            saveat=50
           );

domain  = Domain(extent=(2π,2π),res=(8,8));
# domain  = Domain(extent=(2π,2π),res=(12,12));
coeffs  = Coefficients(Ω=2π,θ=0.0,μ=0.05,ν=0.0,ν₄=0.0);
forcing = PointJet(Ξ=1.0,Δθ=0.1,τ=20.0);
prob    = BetaPlane(domain,coeffs,forcing);

dn      = "validate/pointjet/8x8/"
# dn      = "validate/pointjet/12x12/"

eqs     = [NL(),GQL(Λ=prob.d.nx-1),GCE2(Λ=prob.d.nx-1),GQL(0),CE2()];
labels  = ["nl","gql","gce2","ql","ce2"];
sols    = integrate(prob,eqs,tspan;tsargs...);
write(prob,eqs,sols,dn=dn,labels=labels)

gce2_0    = integrate(prob,GCE2(0),tspan;tsargs...);
write(prob,GCE2(0),gce2_0,dn=dn,fn="gce2_0")

gql    = integrate(prob,GQL(1),tspan;tsargs...);
write(prob,GQL(1),gql,dn=dn,fn="gql_1")

gce2    = integrate(prob,GCE2(Λ=1,poscheck=true,poscheckat=20),tspan;tsargs...);
write(prob,GCE2(1),gce2,dn=dn,fn="gce2_1")

gql    = integrate(prob,GQL(3),tspan;tsargs...);
write(prob,GQL(3),gql,dn=dn,fn="gql_3")

gce2    = integrate(prob,GCE2(Λ=3,poscheck=true,poscheckat=20),tspan;tsargs...);
write(prob,GCE2(3),gce2,dn=dn,fn="gce2_3")

gql    = integrate(prob,GQL(5),tspan;tsargs...);
write(prob,GQL(5),gql,dn=dn,fn="gql_5")

gce2    = integrate(prob,GCE2(5),tspan;tsargs...);
write(prob,GCE2(5),gce2,dn=dn,fn="gce2_5")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule BiCePS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("BiCePS")
JLLWrappers.@generate_main_file("BiCePS", UUID("7196b6e2-84d2-5468-aff3-a0acff4d3b58"))
end  # module BiCePS_jll


text:
using .Zygote: Zygote, @adjoint, @nograd, pullback

@adjoint istraining() = true, _ -> nothing
@nograd Bijectors._debug

@adjoint function mapvcat(f, args...)
    g(f, args...) = map(f, args...)
    return pullback(g, f, args...)
end
@adjoint function eachcolmaphcat(f, x1, x2)
    function g(f, x1, x2)
        init = reshape(f(view(x1, :, 1), x2[1]), :, 1)
        return reduce(hcat, [f(view(x1, :, i), x2[i]) for i in 2:size(x1, 2)]; init = init)
    end
    return pullback(g, f, x1, x2)
end
@adjoint function eachcolmaphcat(f, x)
    function g(f, x)
        init = reshape(f(view(x, :, 1)), :, 1)
        return reduce(hcat, [f(view(x, :, i)) for i in 2:size(x, 2)]; init = init)
    end
    return pullback(g, f, x)
end
@adjoint function _sum(f, args...)
    g(f, args...) = sum(map(f, args...))
    return pullback(g, f, args...)
end
@adjoint function _sumeachcol(f, x1, x2)
    g(f, x1, x2) = sum([f(view(x1, :, i), x2[i]) for i in 1:size(x1, 2)])
    return pullback(g, f, x1, x2)
end

@adjoint function logabsdetjac(b::Log{1}, x::AbstractVector)
    return -sum(log, x), Δ -> (nothing, -Δ ./ x)
end
@adjoint function logabsdetjac(b::Log{1}, x::AbstractMatrix)
    return -vec(sum(log, x; dims = 1)), Δ -> (nothing, .- Δ' ./ x)
end

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:ZygoteAD}, Inverse{<:ADBijector{<:ZygoteAD}}},
    x::Real
)
    return Zygote.gradient(b, x)[1]
end
function jacobian(
    b::Union{<:ADBijector{<:ZygoteAD}, Inverse{<:ADBijector{<:ZygoteAD}}},
    x::AbstractVector{<:Real}
)
    return Zygote.jacobian(b, x)
end
@adjoint function _logabsdetjac_scale(a::Real, x::Real, ::Val{0})
    return _logabsdetjac_scale(a, x, Val(0)), Δ -> (inv(a) .* Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0})
    J = fill(inv.(a), length(x))
    return _logabsdetjac_scale(a, x, Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{0})
    J = fill(size(x, 1) / a, size(x, 2))
    return _logabsdetjac_scale(a, x, Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::AbstractVector, x::AbstractVector, ::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    J = inv.(a)
    return _logabsdetjac_scale(a, x, Val(1)), Δ -> (J .* Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Val{1})
    Jᵀ = repeat(inv.(a), 1, size(x, 2))
    return _logabsdetjac_scale(a, x, Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end

## Positive definite matrices

@adjoint function replace_diag(::typeof(log), X)
    f(i, j) = i == j ? log(X[i, j]) : X[i, j]
    out = f.(1:size(X, 1), (1:size(X, 2))')
    out, ∇ -> begin
        g(i, j) = i == j ? ∇[i, j] / X[i, j] : ∇[i, j]
        (nothing, g.(1:size(X, 1), (1:size(X, 2))'))
    end
end
@adjoint function replace_diag(::typeof(exp), X)
    f(i, j) = ifelse(i == j, exp(X[i, j]), X[i, j])
    out = f.(1:size(X, 1), (1:size(X, 2))')
    out, ∇ -> begin
        g(i, j) = ifelse(i == j, ∇[i, j] * exp(X[i, j]), ∇[i, j])
        (nothing, g.(1:size(X, 1), (1:size(X, 2))'))
    end
end

@adjoint function _logpdf_with_trans_pd(
    d,
    X::AbstractMatrix{<:Real},
    transform::Bool,
)
    return pullback(_logpdf_with_trans_pd_zygote, d, X, transform)
end
function _logpdf_with_trans_pd_zygote(
    d,
    X::AbstractMatrix{<:Real},
    transform::Bool,
)
    T = eltype(X)
    Xcf = cholesky(X, check = false)
    if !issuccess(Xcf)
        Xcf = cholesky(X + max(eps(T), eps(T) * norm(X)) * I, check = true)
    end
    lp = getlogp(d, Xcf, X)
    if transform && isfinite(lp)
        U = Xcf.U
        @inbounds for i in 1:dim(d)
            lp += (dim(d) - i + 2) * log(U[i, i])
        end
        lp += dim(d) * log(T(2))
    end
    return lp
end

# Simplex adjoints

@adjoint function _simplex_bijector(X::AbstractVector, b::SimplexBijector)
    return _simplex_bijector(X, b), Δ -> (simplex_link_jacobian(X)' * Δ, nothing)
end
@adjoint function _simplex_inv_bijector(Y::AbstractVector, b::SimplexBijector)
    return _simplex_inv_bijector(Y, b), Δ -> (simplex_invlink_jacobian(Y)' * Δ, nothing)
end

@adjoint function _simplex_bijector(X::AbstractMatrix, b::SimplexBijector)
    return _simplex_bijector(X, b), Δ -> begin
        maphcat(eachcol(X), eachcol(Δ)) do c1, c2
            simplex_link_jacobian(c1)' * c2
        end, nothing
    end
end
@adjoint function _simplex_inv_bijector(Y::AbstractMatrix, b::SimplexBijector)
    return _simplex_inv_bijector(Y, b), Δ -> begin
        maphcat(eachcol(Y), eachcol(Δ)) do c1, c2
            simplex_invlink_jacobian(c1)' * c2
        end, nothing
    end
end

@adjoint function logabsdetjac(b::SimplexBijector, x::AbstractVector)
    return logabsdetjac(b, x), Δ -> begin
        (nothing, simplex_logabsdetjac_gradient(x) * Δ)
    end
end
@adjoint function logabsdetjac(b::SimplexBijector, x::AbstractMatrix)
    return logabsdetjac(b, x), Δ -> begin
        (nothing, maphcat(eachcol(x), Δ) do c, g
            simplex_logabsdetjac_gradient(c) * g
        end)
    end
end

# LocationScale fix

@adjoint function _clamp(x::T, dist::LocationScale) where {T <: Real}
    function f(x, dist)
        ϵ = _eps(T)
        bounds = (minimum(dist) + ϵ, maximum(dist) - ϵ)
        if x < bounds[1]
            clamped_x = bounds[1]
        elseif x > bounds[2]
            clamped_x = bounds[2]
        else
            clamped_x = x
        end
        DEBUG && _debug("x = $x, bounds = $bounds, clamped_x = $clamped_x")
        return clamped_x
    end
    return pullback(f, x, dist)
end


text:
using MLJBase:
    CV,
    evaluate!,
    fit,
    fit!,
    machine,
    predict,
    rms
using MLJTuringInterface
using Test
using Turing:
    Diagonal,
    I,
    MvNormal,
    Normal,
    NUTS,
    @model

"""
    _allclose(A, B; atol=0.1)

Helper function to return whether the elements in `A` and `B` are approximately equal.
Also throws an informative error if the elements are not approximately equal.
"""
function _allclose(A, B; atol=0.1)
    @assert length(A) == length(B)
    for (i, t) in enumerate(zip(A, B))
        a, b = t
        if !isapprox(a, b; atol)
            error("!($a ≈ $b) at position $i")
        end
    end
    return true
end

@testset "turinginterface" begin
    # Using centered data to speed up sampling.
    X = (; X = collect(-3:3))
    y = collect(-0.9:0.3:0.9)

    @model function linear_regression(X, y)
        σ = 0.2
        intercept ~ Normal(0, σ)
        coef ~ Normal(0, σ)

        mu = intercept .+ X.X * coef
        y ~ MvNormal(mu, σ^2 * I)
    end

    model = linear_regression
    n_samples = 100
    sampler = NUTS()
    tm = TuringModel(model, n_samples, sampler)

    @testset "interface" begin
        verbosity = 1
        fitresult, _, _ = fit(tm, verbosity, X, y)
        chns = fitresult

        @test MLJTuringInterface._parameter_mean(chns, :intercept) ≈ 0 atol=0.05
        @test MLJTuringInterface._parameter_mean(chns, :coef) ≈ 0.3 atol=0.05

        predictions = predict(tm, chns, X)
        @test _allclose(y, predictions)
    end

    @testset "cross-validation" begin
        mach = machine(tm, X, y)
        measure = rms
        nfolds = 2
        resampling = CV(; nfolds)
        result = evaluate!(mach; measure, resampling)
        # Smoke test.
        @test length(result.fitted_params_per_fold) == nfolds
    end
end

@testset "named_data" begin
    A = collect(-3:3)
    n = length(A)
    B = fill(1, n)
    X = (; A, B)
    y = collect(-0.9:0.3:0.9)

    @model function multivariate_regression(X, y)
        σ = 0.2
        intercept ~ Normal(0, σ)
        d = size(X, 2)
        coef ~ MvNormal(Diagonal(fill(σ ^ 2, d)))

        mu = intercept .+ X * coef
        y ~ MvNormal(mu, σ^2 * I)
    end

    model = multivariate_regression
    n_samples = 100
    sampler = NUTS()
    renamer = Dict(["coef[$i]" => "coef[$key]" for (i, key) in enumerate(keys(X))])
    tm = TuringModel(model, n_samples, sampler; renamer)
    mach = machine(tm, X, y)
    fit!(mach)

    chns = mach.fitresult
    @test MLJTuringInterface._parameter_mean(chns, "coef[A]") ≈ 0.3 atol=0.05
end


text:
function add_kepler_contribution(obj::TimeDerivativesEquations, arg0::jdouble)
    return jcall(obj, "addKeplerContribution", void, (jdouble,), arg0)
end

function add_mass_derivative(obj::TimeDerivativesEquations, arg0::jdouble)
    return jcall(obj, "addMassDerivative", void, (jdouble,), arg0)
end

function add_non_keplerian_acceleration(obj::TimeDerivativesEquations, arg0::Vector3D)
    return jcall(obj, "addNonKeplerianAcceleration", void, (Vector3D,), arg0)
end


text:
module Gmm

export gmmObj!, gmmObjOnly, gmmObjGradOnly!, gmmMoment
using TypeDef, Helper, Estim, UserFileIO
using LinearAlgebra

function gmmObj!(g::Vector{Float64}, θ₂::Vector{Float64},
  dat::datType, dim::dimType, fxp::fxpType, mat::matType,
  indx::indxType, pred::predType, diags::diagsType)

  m = gmmMoment(θ₂,dat,dim,fxp,mat,indx,pred,diags)
  gmmGrad!(g,m,dim,fxp,mat,indx,diags)
  return dot(m, mat.invPhi*m)
end

function gmmObjOnly(
  θ₂::Vector{Float64},
  dat::datType,
  dim::dimType,
  fxp::fxpType,
  mat::matType,
  indx::indxType,
  pred::predType,
  diags::diagsType,
  update::Bool=false)

  m = gmmMoment(θ₂,dat,dim,fxp,mat,indx,pred,diags,update)
  return dot(m, mat.invPhi*m)
end
function gmmObjOnly(
  θ₂::Float64,
  dat::datType,
  dim::dimType,
  fxp::fxpType,
  mat::matType,
  indx::indxType,
  pred::predType,
  diags::diagsType,
  update::Bool=false)

  m = gmmMoment([θ₂],dat,dim,fxp,mat,indx,pred,diags,update)
  return dot(m, mat.invPhi*m)
end

function gmmObjGradOnly!(g::Vector{Float64}, θ₂::Vector{Float64},
  dat::datType, dim::dimType, fxp::fxpType, mat::matType,
  indx::indxType, pred::predType, diags::diagsType)

  m = gmmMoment(θ₂,dat,dim,fxp,mat,indx,pred,diags)
  gmmGrad!(g,m,dim,fxp,mat,indx,diags)
end

function gmmGrad!(g::Vector{Float64}, m::Vector{Float64},
  dim::dimType, fxp::fxpType, mat::matType,
  indx::indxType, diags::diagsType)

  ∂δ∂θ = unpackStackedArray(fxp.∂δ∂θ, dim.obs,dim.θ₂,dim.market,indx.mkt)
  g[:] = 2*∂δ∂θ'*mat.Z*mat.invPhi*m
end

function gmmMoment(
  θ₂::Vector{Float64},
  dat::datType,
  dim::dimType,
  fxp::fxpType,
  mat::matType,
  indx::indxType,
  pred::predType,
  diags::diagsType,
  update::Bool=false)

  #setupMu!(fxp.eμ,pred.eμ₀,pred.adj, θ₂,dim,indx,mat.L,mat.Lν,mat.Xᵢ,mat.νᵢ)
  setupMu!(fxp.eμ,pred.eμ₀,pred.adj, θ₂,mat.Xᵢ,mat.Dᵢ,mat.νᵢ,dim,indx)
  fixPoint!(fxp.eδ,fxp.∂δ∂θ, fxp.eμ,dat,dim,mat,indx,pred,diags,θ₂)
  #fixPoint!(fxp.eδ,fxp.∂δ∂θ, fxp.eμ,dat,dim,mat,indx,pred,diags)
  for t=1:dim.market
    if any(isnan.(fxp.eδ[t]))
      # println(" Error: fxp.eδ NaN")
      return NaN*ones(size(mat.Z,2),)
    end
  end
  if update==true
    optimIV!(mat,θ₂,dat,dim,fxp,indx,pred,diags)
  end
  linearIV!(mat.ξ, fxp.eδ,fxp.eδ1d, dim,mat,indx;output=3)
  m  = mat.Z'*mat.ξ
  # m  = view(mat.Z,mat.mask,:)'*view(mat.ξ,mat.mask,1)
  return m
end
# function gmmMoment(θ₂::Float64,
#   dat::datType, dim::dimType, fxp::fxpType, mat::matType,
#   indx::indxType, pred::predType, diags::diagsType, update::Bool=false)
#
#   setupMu!(fxp.eμ,pred.eμ₀,pred.adj, θ₂,dim,mat,indx)
#   fixPoint!(fxp.eδ,fxp.∂δ∂θ, fxp.eμ,dat,dim,mat,indx,pred,diags,θ₂)
#   if update==true
#     optimIV!(mat,θ₂,dat,dim,fxp,indx,pred,diags)
#   end
#   linearIV!(mat.ξ, fxp.eδ,fxp.eδ1d, dim,mat,indx;output=3)
#   m  = mat.Z'*mat.ξ
#   # m  = view(mat.Z,mat.mask,:)'*view(mat.ξ,mat.mask,1)
#   return m
# end

end # end of module

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((im+cot(x))^-1*sin(x)^4, x) == :(-3//32*(im+cot(x))^-2+1//32*(im+-1*cot(x))^-2+-1*im*(-8*cot(x)+8im)^-1+3*im*(16im+16*cot(x))^-1+-5//16*im*x+-1//24*im*(im+cot(x))^-3)
@test integrate((im+cot(x))^-1*sin(x)^3, x) == :(-4//15*im*cos(x)^3+4//5*im*cos(x)+im*(5im+5*cot(x))^-1*sin(x)^3)
@test integrate((im+cot(x))^-1*sin(x)^2, x) == :(-1//8*(im+cot(x))^-2+im*(4im+4*cot(x))^-1+-1*im*(-8*cot(x)+8im)^-1+-3//8*im*x)
@test integrate((im+cot(x))^-1*sin(x), x) == :(2//3*im*cos(x)+im*(3im+3*cot(x))^-1*sin(x))
@test integrate((im+cot(x))^-1*csc(x), x) == :(im*(im+cot(x))^-1*csc(x))
@test integrate((im+cot(x))^-1*csc(x)^2, x) == :(-1*im*x+log(sin(x)))
@test integrate((im+cot(x))^-1*csc(x)^3, x) == :(-1*csc(x)+im*arctanh(cos(x)))
@test integrate((im+cot(x))^-1*csc(x)^4, x) == :(-1//2*cot(x)^2+im*cot(x))
@test integrate((im+cot(x))^-1*csc(x)^5, x) == :(-1//3*csc(x)^3+(1/2)*im*arctanh(cos(x))+(1/2)*im*cot(x)*csc(x))
@test integrate((im+cot(x))^-1*csc(x)^6, x) == :(-1//2*cot(x)^2+-1//4*cot(x)^4+im*cot(x)+1//3*im*cot(x)^3)
@test integrate((im+cot(x))^-1*csc(x)^7, x) == :(-1//5*csc(x)^5+3//8*im*arctanh(cos(x))+1//4*im*csc(x)^3*cot(x)+3//8*im*cot(x)*csc(x))
@test integrate((a+b*cot(x))^-1*csc(x)^6, x) == :(-1//4*b^-1*cot(x)^4+-1*b^-5*(a^2+b^2)^2*log(a+b*cot(x))+-1//2*b^-3*cot(x)^2*(a^2+2*b^2)+1//3*a*b^-2*cot(x)^3+a*b^-4*(a^2+2*b^2)*cot(x))
@test integrate((a+b*cot(x))^-1*csc(x)^4, x) == :(-1//2*b^-1*cot(x)^2+a*b^-2*cot(x)+-1*b^-3*(a^2+b^2)*log(a+b*cot(x)))
@test integrate((a+b*cot(x))^-1*csc(x)^2, x) == :(-1*b^-1*log(a+b*cot(x)))
@test integrate((a+b*cot(x))^-1*sin(x)^2, x) == :(-1*b^3*(a^2+b^2)^-2*log(a*sin(x)+b*cos(x))+-1*(2*a^2+2*b^2)^-1*sin(x)^2*(b+a*cot(x))+(1/2)*a*x*(a^2+b^2)^-2*(a^2+3*b^2))
@test integrate((a+b*cot(x))^-1*sin(x)^4, x) == :(-1*b^5*(a^2+b^2)^-3*log(a*sin(x)+b*cos(x))+-1*(4*a^2+4*b^2)^-1*sin(x)^4*(b+a*cot(x))+-1//8*(a^2+b^2)^-2*sin(x)^2*(4*b^3+a*(3*a^2+7*b^2)*cot(x))+1//8*a*x*(a^2+b^2)^-3*(3*a^4+15*b^4+10*a^2*b^2))
@test integrate((a+b*cot(x))^-1*csc(x)^5, x) == :(-1//3*b^-1*csc(x)^3+b^-4*(a^2+b^2)^3//2*arctanh((a^2+b^2)^-1//2*(b+-1*a*cot(x))*sin(x))+(1/2)*a*b^-2*arctanh(cos(x))+-1*b^-3*(a^2+b^2)*csc(x)+a*b^-4*(a^2+b^2)*arctanh(cos(x))+(1/2)*a*b^-2*cot(x)*csc(x))
@test integrate((a+b*cot(x))^-1*csc(x)^3, x) == :(-1*b^-1*csc(x)+a*b^-2*arctanh(cos(x))+b^-2*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(b+-1*a*cot(x))*sin(x)))
@test integrate((a+b*cot(x))^-1*csc(x), x) == :(-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(-1b+a*cot(x))*sin(x)))
@test integrate((a+b*cot(x))^-1*sin(x), x) == :(b^2*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(b+-1*a*cot(x))*sin(x))+-1*a*(a^2+b^2)^-1*cos(x)+-1*b*(a^2+b^2)^-1*sin(x))
@test integrate((a+b*cot(x))^-1*sin(x)^3, x) == :(a*(3*a^2+3*b^2)^-1*cos(x)^3+b^4*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(b+-1*a*cot(x))*sin(x))+-1*a*(a^2+b^2)^-1*cos(x)+-1*b*(3*a^2+3*b^2)^-1*sin(x)^3+-1*b^3*(a^2+b^2)^-2*sin(x)+-1*a*b^2*(a^2+b^2)^-2*cos(x))
@test integrate((a+b*cot(x))^-2*csc(x)^2, x) == :(b^-1*(a+b*cot(x))^-1)
@test integrate((a+b*cot(x))^n*csc(x)^2, x) == :(-1*b^-1*(1+n)^-1*(a+b*cot(x))^(1+n))

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.1.2 (d csc)^m (a+b cot)^n.jl


text:
include(joinpath("..", "src_PolynomialOptim", "PolynomialOptim.jl"))
using LightXML, MAT

import Base.copy

"""
    Link(orig, dest)

Represents a link from a bus origin `orig` to a bus destination `dest`
Util for stock or treat data about transmission lines or transformers
### Examples
```
julia > Link("BUS_1","BUS_4") #represents a link between BUS_1 and BUS_4
```
"""
struct Link
  orig::String
  dest::String
end

"""
    Abstract type AbstractNodeLabel

Abstract type for any element about a node
### Examples
```julia
struct MatpowerLoad <: AbstractNodeLabel
struct GOCVolt <: AbstractNodeLabel
struct IIDMGenerator <: AbstractNodeLabel
struct GOCShunt <: AbstractNodeLabel

```
"""
abstract type AbstractNodeLabel end

"""
    abstract type AbstractLinkLabel

Abstract type for any element about a link
### Examples
```
struct MatpowerLine_π <: AbstractLinkLabel
struct GOCLineπ_withtransformer <: AbstractLinkLabel
struct GOCNullImpedance_notransformer <: AbstractLinkLabel
```
"""
abstract type AbstractLinkLabel end

"""
    abstract type AbstractInput
Abstract type for any entry format
### Examples
```
type MatpowerInput <: AbstractInput end
type GOCInput <:AbstractInput end
type IIDMInput <: AbstractInput end
```
"""
abstract type AbstractInput end

"""
    DataSource(bus::Dict{String, Dict{String, Any}}, link::Dict{Link, Dict{String, Any}})

Store network data either in an attribute `bus` or in an attribute `link` depending on data type.

### Attributes
- `bus` : Dict{String, Dict{String, Any}}
- `link`: Dict{Link, Dict{String, Any}}

### Example
```
Instance Matpower WB2 with two nodes
julia > databus1 = Dict("Volt" => MatpowerVolt("BUS_1",1,0.95, 1.05), "Gen_1"=> MatpowerGenerator(-400 im, 600+400im, 3, [0 2 0], true))
julia > databus2 = Dict("Volt" => MatpowerVolt("BUS_2",2,0.95, 1.05), "Load" => MatpowerLoad(350-350im))
julia > busdata = Dict("BUS_1"=> databus1, "BUS_2" => databus2)
julia > link1to2 = Dict("LinkMP_1" => MatpowerLine_π(100,0.04,0.2,0.0,0.0,0.0))
julia > linkdata = Dict( Link("BUS_1","BUS_2") => link1to2)
julia > ds = DataSource(busdata,linkdata)
```
"""
struct DataSource
  # baseMVA::Number
  bus::Dict{String, Dict{String, Any}}
  link::Dict{Link, Dict{String, Any}}
end

"""
    struct GridStructure(scenario::String, node_linksin::Dict{String, Set{Link}}, node_linksout::Dict{String, Set{Link}})

Store information about network to construct power balances.

### Attributes
- `scenario::String` : scenario to study
- `generator_types::Set{Type}` : specify which types are generators
- `node_linksin::Dict{String, Set{Link}}` : for each node, set of links having this node for destination
- `node_linksout::Dict{String, Set{Link}}` : for each node, set of links having this node for origin

### Example
```
Instance Matpower WB2 with two nodes
julia > scenario = "BaseCase"
julia > node_linksin = Dict( "BUS_2" => Set{Link}([Link("BUS_1","BUS_2")]) )
julia > node_linksout = Dict( "BUS_1" => Set{Link}([Link("BUS_1, "BUS_2")]) )
julia > gs = GridStructure(scenario,generator_types, node_linksin, node_linksout)
```
"""
mutable struct GridStructure
  scenario::String
  node_linksin::Dict{String, Set{Link}}
  node_linksout::Dict{String, Set{Link}}
end

"""
    struct MathematicalProgramming(node_formulations::Dict{String, Dict{String, Symbol}}, link_formulations::Dict{String, Dict{String, Symbol}}, node_vars::Dict{String, Dict{String, Variable}}, link_vars::Dict{Link, Dict{String, Variable}})

Store information about mathematical formulation.
A node element or a link element can be formulated with a minimal number of variables or with variables representing all the quantities linked to the element.
This formulation is represented by symbols in the dictonaries `node_formulations` or `link_formulations` depending on the type of element.
Store created variables in `node_vars` or `link_vars`.

### Example
```
Instance Matpower WB2 with two nodes
julia > node_formulations =
julia > link_formulations = Dict( Link("BUS_1","BUS_2") => Dict("LinkMP_1" => :NbMinVar))
julia > node_vars =
julia > link_vars =

```
"""
mutable struct MathematicalProgramming
  node_formulations::Dict{String, Dict{String, Symbol}} # :NbMinVar, :NewVar, :None
  link_formulations::Dict{Link, Dict{String, Symbol}}
  node_vars::Dict{String, Dict{String, Variable}}
  link_vars::Dict{Link, Dict{String, Variable}}
end

"""
    struct Scenario(ds::DataSource, gs::GridStructure, mp::MathematicalProgramming)

### Example
```
Instance Matpower WB2 with two nodes
julia > basecase_scenario = Scenario(ds, gs, mp)
```
"""
mutable struct Scenario
  ds::DataSource
  gs::GridStructure
  mp::MathematicalProgramming
end

"""
    OPFproblems

### Examples
```

```
"""
const OPFProblems = Dict{String, Scenario}

# files to include for each entry type
include(joinpath("..","src_GOC","GOC_files.jl"))
include(joinpath("..","src_Matpower","Matpower_files.jl"))
include(joinpath("..","src_IIDM","IIDM_files.jl"))

# common files
include("add_generator_variables.jl")
include("build_globalpb.jl")
include("build_Problem.jl")
include("check_feasibility.jl")
include("default_behaviour.jl")
include("load_OPFproblems.jl")
include("naming_conventions.jl")
include("PowSysMod_accessors.jl")
include("PowSysMod_operators.jl")
include("utils_plots.jl")

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
include("experience_based_sampling_model.jl")
include("time_based_sample_model.jl")
include("prioritized_sweeping_sampling_model.jl")


text:
using DataFrames, DataFramesMeta, CSV
using PGFPlots
using Statistics, HypothesisTests

""" Functions for results of Fig.1B and Fig.1C """
# Get mean RT per subject and averaged across subjects
# for surprise and non-surprise trials
# Check statistical significance of difference
function processreactiontime(;
        filepath =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/reactiontime_persubj/"),
        )

    savepath = makehomesavepath("processreactiontime")
    mkdir(savepath)

    subjectfiles = readdir(filepath)

    dfresults = DataFrame(subjectname = String[],
                    mean_rt_surprise = Float64[],
                    std_rt_surprise = Float64[],
                    se_rt_surprise = Float64[],
                    mean_rt_nonsurprise = Float64[],
                    std_rt_nonsurprise = Float64[],
                    se_rt_nonsurprise = Float64[],
                    mean_rt_diff = Float64[],
                    )

    dfresults_after4 = DataFrame(subjectname = String[],
                    mean_rt_surprise = Float64[],
                    std_rt_surprise = Float64[],
                    se_rt_surprise = Float64[],
                    mean_rt_nonsurprise = Float64[],
                    std_rt_nonsurprise = Float64[],
                    se_rt_nonsurprise = Float64[],
                    mean_rt_diff = Float64[],
                    )
    dfresult_scatter = DataFrame()
    dfresult_scatter_after4 = DataFrame()
    dfresult_scatter_diff = DataFrame()
    dfresult_scatter_diff_after4 = DataFrame()
    for i in subjectfiles
        @show i
        # subjectname = split(split(i, "_")[end], ".")[end-1]
        dfalgo = CSV.read(joinpath(filepath, i), header=false)
        rt_surp = @where(dfalgo, :Column2 .== 1)[!, :Column1]
        rt_nonsurp = @where(dfalgo, :Column2 .== 0)[!, :Column1]

        # Remove the first one cause some subjects had their eyes closed
        # and did not realize the experiment has started:
        rt_nonsurp = rt_nonsurp[2:end]

        rt_surp = rt_removeoutliers(rt_surp, showstr = "surp")
        rt_nonsurp = rt_removeoutliers(rt_nonsurp, showstr = "non_surp")
        # ---------
        #
        mean_rt_surprise, std_rt_surprise, se_rt_surprise = calc_stats_nan(reshape(rt_surp, :, 1)') # reshape to make it 2-d (datapoints x 1) and transpose
        mean_rt_nonsurprise, std_rt_nonsurprise, se_rt_nonsurprise = calc_stats_nan(reshape(rt_nonsurp, :, 1)') # reshape to make it 2-d (datapoints x 1) and transpose

        push!(dfresults, [i, mean_rt_surprise[1], std_rt_surprise[1], se_rt_surprise[1],
                                        mean_rt_nonsurprise[1], std_rt_nonsurprise[1], se_rt_nonsurprise[1],
                                        mean_rt_surprise[1] - mean_rt_nonsurprise[1]
                                        ])
        # In columns: first nonsurprise, then surprise
        dfresult_scatter[!, Symbol(i*"_mean_rt")] = [mean_rt_nonsurprise[1], mean_rt_surprise[1]]
        dfresult_scatter[!, Symbol(i*"_std_rt")] = [std_rt_nonsurprise[1], std_rt_surprise[1]]
        dfresult_scatter[!, Symbol(i*"_se_rt")] = [se_rt_nonsurprise[1], se_rt_surprise[1]]

        # ----------------------------------------------------------------------
        # Remove the data points before the 5th episode (introduction of surprise trials)
        # ----------------------------------------------------------------------
        rt_surp_after4 = @where(dfalgo, :Column2 .== 1, :Column3 .> 4)[!, :Column1]
        rt_nonsurp_after4 = @where(dfalgo, :Column2 .== 0, :Column3 .> 4)[!, :Column1]

        # ----- Remove outliers ----
        rt_surp_after4 = rt_removeoutliers(rt_surp_after4, showstr = "surp - after4")
        rt_nonsurp_after4 = rt_removeoutliers(rt_nonsurp_after4, showstr = "non_surp - after4")
        # ---------

        mean_rt_surprise_after4, std_rt_surprise_after4, se_rt_surprise_after4 = calc_stats_nan(reshape(rt_surp_after4, :, 1)') # reshape to make it 2-d (datapoints x 1) and transpose
        mean_rt_nonsurprise_after4, std_rt_nonsurprise_after4, se_rt_nonsurprise_after4 = calc_stats_nan(reshape(rt_nonsurp_after4, :, 1)') # reshape to make it 2-d (datapoints x 1) and transpose

        push!(dfresults_after4, [i, mean_rt_surprise_after4[1], std_rt_surprise_after4[1], se_rt_surprise_after4[1],
                                        mean_rt_nonsurprise_after4[1], std_rt_nonsurprise_after4[1], se_rt_nonsurprise_after4[1],
                                        mean_rt_surprise_after4[1] - mean_rt_nonsurprise_after4[1]
                                        ])
        # In columns: first nonsurprise, then surprise
        dfresult_scatter_after4[!, Symbol(i*"_mean_rt")] = [mean_rt_nonsurprise_after4[1], mean_rt_surprise_after4[1]]
        dfresult_scatter_after4[!, Symbol(i*"_std_rt")] = [std_rt_nonsurprise_after4[1], std_rt_surprise_after4[1]]
        dfresult_scatter_after4[!, Symbol(i*"_se_rt")] = [se_rt_nonsurprise_after4[1], se_rt_surprise_after4[1]]

    end
    CSV.write(joinpath(savepath, "reactiontime_stats.csv"), dfresults, delim = " ");
    CSV.write(joinpath(savepath, "reactiontime_stats_scatter.csv"), dfresult_scatter, delim = " ");

    CSV.write(joinpath(savepath, "reactiontime_stats_after4.csv"), dfresults_after4, delim = " ");
    CSV.write(joinpath(savepath, "reactiontime_stats_scatter_after4.csv"), dfresult_scatter_after4, delim = " ");

    println("All non-surprise trials:")
    @show OneSampleTTest(dfresults[!,:mean_rt_surprise], dfresults[!,:mean_rt_nonsurprise])
    pval = pvalue(OneSampleTTest(dfresults[!,:mean_rt_surprise], dfresults[!,:mean_rt_nonsurprise]))
    println(" ------ ")
    println("After removing trials before the 5th episode:")
    @show OneSampleTTest(dfresults_after4[!,:mean_rt_surprise], dfresults_after4[!,:mean_rt_nonsurprise])
    pval_after4 = pvalue(OneSampleTTest(dfresults_after4[!,:mean_rt_surprise], dfresults_after4[!,:mean_rt_nonsurprise]))

    @show length(dfresults[:,:mean_rt_surprise])
    @show sqrt(length(dfresults[:,:mean_rt_surprise]))
    dfresults_meanall = DataFrame()
    dfresults_meanall[!, Symbol("mean_rt_surprise")] = [mean(dfresults[:,:mean_rt_surprise])]
    dfresults_meanall[!, Symbol("std_rt_surprise")] = [std(dfresults[:,:mean_rt_surprise])]
    dfresults_meanall[!, Symbol("se_rt_surprise")] = [std(dfresults[:,:mean_rt_surprise]) / sqrt(length(dfresults[:,:mean_rt_surprise]))]
    dfresults_meanall[!, Symbol("mean_rt_nonsurprise")] = [mean(dfresults[:,:mean_rt_nonsurprise])]
    dfresults_meanall[!, Symbol("std_rt_nonsurprise")] = [std(dfresults[:,:mean_rt_nonsurprise])]
    dfresults_meanall[!, Symbol("se_rt_nonsurprise")] = [std(dfresults[:,:mean_rt_nonsurprise]) / sqrt(length(dfresults[:,:mean_rt_nonsurprise]))]
    #
    dfresults_meanall_after4 = DataFrame()
    dfresults_meanall_after4[!, Symbol("mean_rt_surprise")] = [mean(dfresults_after4[:,:mean_rt_surprise])]
    dfresults_meanall_after4[!, Symbol("std_rt_surprise")] = [std(dfresults_after4[:,:mean_rt_surprise])]
    dfresults_meanall_after4[!, Symbol("se_rt_surprise")] = [std(dfresults_after4[:,:mean_rt_surprise]) / sqrt(length(dfresults_after4[:,:mean_rt_surprise]))]
    dfresults_meanall_after4[!, Symbol("mean_rt_nonsurprise")] = [mean(dfresults_after4[:,:mean_rt_nonsurprise])]
    dfresults_meanall_after4[!, Symbol("std_rt_nonsurprise")] = [std(dfresults_after4[:,:mean_rt_nonsurprise])]
    dfresults_meanall_after4[!, Symbol("se_rt_nonsurprise")] = [std(dfresults_after4[:,:mean_rt_nonsurprise]) / sqrt(length(dfresults_after4[:,:mean_rt_nonsurprise]))]

    CSV.write(joinpath(savepath, "reactiontime_meanall.csv"), dfresults_meanall, delim = " ");

    CSV.write(joinpath(savepath, "reactiontime_meanall_after4.csv"), dfresults_meanall_after4, delim = " ");

    dfresults, pval, dfresults_after4, pval_after4,
    dfresult_scatter, dfresult_scatter_after4,
    dfresults_meanall, dfresults_meanall_after4
end

function rt_removeoutliers(reactiontimes;
                            thrs = 8.,
                            showstr = "surp"
                            )
    k  = findall(reactiontimes .> thrs)
    # @show rt_surp[ksurp]
    if ~isempty(k)
        @show k
        println(showstr)
        @show reactiontimes[k]
    end
    reactiontimes = deleteat!(reactiontimes, k)

    reactiontimes
end


text:
using HorizonSideRobots
include("../basics.jl")

md"""
    ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля

    РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы
"""

"""Поиск в глубину, робот проходится по всем клеткам и маркирует их"""
function dfs!(r, visited, node, side)
    is_not_at_start = side != false
    if is_not_at_start
        cheker = isborder(r, side) == false
    else
        cheker = true
    end
    if (node in visited) == false && cheker
        if is_not_at_start
            move!(r, side)
        end
        putmarker!(r)
        push!(visited, node)

        for tempSide in [Nord, Ost, Sud, West]
            if dfs!(r, visited, side_to_x_y(tempSide, node), tempSide)
                move!(r, opposite_side(tempSide))
            end
        end
        return true
    else
        return false
    end
end

 """Решение задачи 3"""
 function solve_problem_3!(r)
    visited = Set()
    dfs!(r, visited, (0, 0), false)
    show(r)
end


text:
using ROCKS
using Random
using Distributions

Random.seed!(888)
const x = rand(Uniform(-5, 5), 1_000_000)
const logit = -3.0 .+ 0.5 .* x .+ rand(Normal(0, 0.1), length(x))
const prob = @. 1.0 / (1.0 + exp(-logit))
const target = rand(length(x)) .<= prob

kk = kstest(target, prob)

cc = concordance(target, prob)

perf = bcdiag(target, prob)

kp = ksplot(perf)
png(kp, "ksplot.png")

rp = rocplot(perf)
png(rp, "rocplot.png")

bp = biasplot(perf)
png(bp, "biasplot.png")

ap = accuracyplot(perf)
png(ap, "accuracyplot.png")

lp = liftcurve(perf)
png(lp, "liftcurve.png")

clp = cumliftcurve(perf)
png(clp, "cumliftcurve.png")

liftable(perf)

cumliftable(perf)


text:
"""
A general TransientFESolver
"""
struct TransientFESolver
  odes::ODESolver
end

function solve(
  solver::TransientFESolver,op::TransientFEOperator,u0,t0::Real,tf::Real)
  TransientFESolution(solver,op,u0,t0,tf)
end

function solve(
  solver::TransientFESolver,op::TransientFEOperator,u0,v0,a0,t0::Real,tf::Real)
  TransientFESolution(solver,op,u0,v0,a0,t0,tf)
end

function test_transient_fe_solver(solver::TransientFESolver,op::TransientFEOperator,u0,t0,tf)
  solution = solve(solver,op,u0,t0,tf)
  test_transient_fe_solution(solution)
end


text:
include(joinpath(dirname(@__FILE__), "..", "src", "AtanhErf.jl"))

if !success(`which lockfile`)
    warn("could not find the `lockfile` program. Try installing it (it's included with `procmail` – www.procmail.org).")
end


text:
@testset "Metric Wrapper" begin
	X, y = @load_toydata
	model = ConstantClassifier()
	result = evaluate(model, X, y, measure=MetricWrapper(tpr, grp=:Sex))
	@test length(keys(result.per_fold[1][1])) == 3
	result = evaluate(model, X, y, measures=[cross_entropy, MetricWrappers([tpr, fnr], grp=:Sex)...])
	@test length(keys(result.per_fold[2][1])) == 3
	@test length(result.per_fold) == 3
end


text:
using CoordinateTransformations: IdentityTransformation, Transformation
using LinearAlgebra: I

export FractionalFromCartesian,
    CartesianFromFractional,
    FractionalToCartesian,
    CartesianToFractional,
    StandardizedFromPrimitive,
    PrimitiveFromStandardized,
    PrimitiveToStandardized,
    StandardizedToPrimitive

struct CartesianFromFractional{T} <: Transformation
    tf::SMatrix{3,3,T,9}
end
struct FractionalFromCartesian{T} <: Transformation
    tf::SMatrix{3,3,T,9}
end
# This requires the a-vector is parallel to the Cartesian x-axis.
# See https://en.wikipedia.org/wiki/Fractional_coordinates
"""
    CartesianFromFractional(lattice::Union{Lattice,ReciprocalLattice})
    CartesianFromFractional(a, b, c, α, β, γ)

Get the transformation from fractional coordinates to Cartesian coordinates.
"""
CartesianFromFractional(lattice::Lattice) = CartesianFromFractional(lattice.data)
CartesianFromFractional(lattice::ReciprocalLattice) =
    CartesianFromFractional(transpose(lattice.data))
function CartesianFromFractional(a, b, c, α, β, γ)
    Ω = cellvolume(a, b, c, α, β, γ)
    b_sinγ, b_cosγ = b .* (sind(γ), cosd(γ))
    return CartesianFromFractional(
        [
            a b_cosγ c*cosd(β)
            0 b_sinγ c*_auxiliary(α, β, γ)
            0 0 Ω/(a*b_sinγ)
        ],
    )
end
"""
    FractionalFromCartesian(lattice::Union{Lattice,ReciprocalLattice})
    FractionalFromCartesian(a, b, c, α, β, γ)

Get the transformation from Cartesian coordinates to fractional coordinates.
"""
FractionalFromCartesian(lattice::Lattice) = FractionalFromCartesian(inv(lattice.data))
FractionalFromCartesian(lattice::ReciprocalLattice) =
    FractionalFromCartesian(transpose(inv(lattice.data)))
function FractionalFromCartesian(a, b, c, α, β, γ)
    Ω = cellvolume(a, b, c, α, β, γ)
    b_sinγ = b * sind(γ)
    return FractionalFromCartesian(
        [
            1/a -cotd(γ)/a -b*c*_auxiliary(β, α, γ)/Ω
            0 1/b_sinγ -a*c*_auxiliary(α, β, γ)/Ω
            0 0 a*b_sinγ/Ω
        ],
    )
end
const FractionalToCartesian = CartesianFromFractional
const CartesianToFractional = FractionalFromCartesian

# This is a helper function and should not be exported!
_auxiliary(α, β, γ) = (cosd(α) - cosd(β) * cosd(γ)) / sind(γ)

(x::Union{CartesianFromFractional,FractionalFromCartesian})(v) = x.tf * collect(v)
(x::Union{CartesianFromFractional,FractionalFromCartesian})(p::ReciprocalPoint) =
    ReciprocalPoint(x.tf * collect(p.coord), p.weight)

Base.inv(x::FractionalFromCartesian) = CartesianFromFractional(inv(x.tf))
Base.inv(x::CartesianFromFractional) = FractionalFromCartesian(inv(x.tf))
Base.:∘(x::CartesianFromFractional, y::FractionalFromCartesian) = ∘(y, x)
Base.:∘(x::FractionalFromCartesian, y::CartesianFromFractional) =
    x.tf * y.tf ≈ I ? IdentityTransformation() : error("undefined!")

# Idea from https://spglib.github.io/spglib/definition.html#transformation-to-the-primitive-cell
struct StandardizedFromPrimitive{T} <: Transformation
    tf::SMatrix{3,3,T,9}
end
struct PrimitiveFromStandardized{T} <: Transformation
    tf::SMatrix{3,3,T,9}
end
"""
    PrimitiveFromStandardized(tf::AbstractMatrix)

Construct the transformation from a standardized cell to a primitive cell.
"""
PrimitiveFromStandardized(tf::AbstractMatrix) = PrimitiveFromStandardized{eltype(tf)}(tf)
"""
    StandardizedFromPrimitive(tf::AbstractMatrix)

Construct the transformation from a primitive cell to a standardized cell.
"""
StandardizedFromPrimitive(tf::AbstractMatrix) = StandardizedFromPrimitive{eltype(tf)}(tf)
const PrimitiveToStandardized = StandardizedFromPrimitive
const StandardizedToPrimitive = PrimitiveFromStandardized

(x::Union{StandardizedFromPrimitive,PrimitiveFromStandardized})(v) = inv(x.tf) * collect(v)
(x::PrimitiveFromStandardized)(lattice::Lattice) = Lattice(lattice.data * x.tf)
(x::StandardizedFromPrimitive)(lattice::Lattice) = Lattice(lattice.data * inv(x.tf))

PrimitiveFromStandardized(::ACentering) = PrimitiveFromStandardized([
    1 0 0
    0 1//2 -1//2
    0 1//2 1//2
])
PrimitiveFromStandardized(::CCentering) = PrimitiveFromStandardized([
    1//2 1//2 0
    -1//2 1//2 0
    0 0 1
])
PrimitiveFromStandardized(::RhombohedralCentering) = PrimitiveFromStandardized([
    2//3 -1//3 -1//3
    1//3 1//3 -2//3
    1//3 1//3 1//3
])
PrimitiveFromStandardized(::BodyCentering) = PrimitiveFromStandardized([
    -1//2 1//2 1//2
    1//2 -1//2 1//2
    1//2 1//2 -1//2
])
PrimitiveFromStandardized(::FaceCentering) = PrimitiveFromStandardized([
    0 1//2 1//2
    1//2 0 1//2
    1//2 1//2 0
])
StandardizedFromPrimitive(::ACentering) = StandardizedFromPrimitive([
    1 0 0
    0 1 1
    0 -1 1
])
StandardizedFromPrimitive(::CCentering) = StandardizedFromPrimitive([
    1 -1 0
    1 1 0
    0 0 1
])
StandardizedFromPrimitive(::RhombohedralCentering) = StandardizedFromPrimitive([
    1 0 1
    -1 1 1
    0 -1 1
])
StandardizedFromPrimitive(::BodyCentering) = StandardizedFromPrimitive([
    0 1 1
    1 0 1
    1 1 0
])
StandardizedFromPrimitive(::FaceCentering) = StandardizedFromPrimitive([
    -1 1 1
    1 -1 1
    1 1 -1
])

Base.inv(x::StandardizedFromPrimitive) = PrimitiveFromStandardized(inv(x.tf))
Base.inv(x::PrimitiveFromStandardized) = StandardizedFromPrimitive(inv(x.tf))
Base.:∘(x::PrimitiveFromStandardized, y::StandardizedFromPrimitive) = ∘(y, x)
Base.:∘(x::StandardizedFromPrimitive, y::PrimitiveFromStandardized) =
    x.tf * y.tf ≈ I ? IdentityTransformation() : error("undefined!")

const ChangeOfBasis{T} = Union{
    CartesianFromFractional{T},
    FractionalFromCartesian{T},
    StandardizedFromPrimitive{T},
    PrimitiveFromStandardized{T},
}
Base.iterate(x::ChangeOfBasis) = iterate(x.tf)
Base.iterate(x::ChangeOfBasis, state) = iterate(x.tf, state)

Base.eltype(::Type{<:ChangeOfBasis{T}}) where {T} = T

Base.length(::ChangeOfBasis) = 9

Base.size(::ChangeOfBasis) = (3, 3)
Base.size(::ChangeOfBasis, dim::Integer) = dim <= 2 ? 3 : 1

Base.IteratorSize(::Type{<:ChangeOfBasis}) = Base.HasShape{2}()

# Enables `firstindex(x, dim)` and `x[1, 2:end]` or `x[begin:2, 2]`
Base.axes(x::ChangeOfBasis, dim::Integer) = axes(x.tf, dim)

Base.getindex(x::ChangeOfBasis, i) = getindex(x.tf, i)
Base.getindex(x::ChangeOfBasis, I::Vararg) = getindex(x.tf, I...)

Base.firstindex(::ChangeOfBasis) = 1

Base.lastindex(::ChangeOfBasis) = 9

function Base.show(io::IO, x::ChangeOfBasis)
    if get(io, :compact, false) || get(io, :typeinfo, nothing) == typeof(x)
        Base.show_default(IOContext(io, :limit => true), x)  # From https://github.com/mauro3/Parameters.jl/blob/ecbf8df/src/Parameters.jl#L556
    else
        println(io, string(typeof(x)))
        for row in eachrow(x.tf)
            println(io, ' ', join(row, " "))
        end
    end
end


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

include("forwarddiff_utils.jl")
include("forwarddiff.jl")


text:
using LinearAlgebra

A = rand(Int8, 2, 2)
U, S, V = svd(A)

println("A=", A)
println("U=", U)
println("S=", S)
println("V=", V)
println("USV=", U*Diagonal(S)*V')