text:
function sbxsubrm(rm,epochs,rois;fun=nothing)
    nepoch = size(epochs,1)   # number of trial
    minepochlength = floor(Int,minimum(diff(epochs,dims=2)))
    ys=Array{Float64}(undef,length(rois),minepochlength,nepoch)
    for i in 1:nepoch
        y = rm[rois,range(max(1,epochs[i,1]),length=minepochlength)]
        if !isnothing(fun)
            y=fun(y)
        end
        ys[:,:,i] = y
    end
    return nepoch==1 ? dropdims(ys,dims=3) : ys
end

function dFoF(bi)
        return x->begin
        bl = mean(x[:,bi],dims=2)
        x = x./bl .-1
        x
        end
    end

sbxcondin(ctc::Dict)=sbxcondin(DataFrame(ctc))
function sbxcondin(ctc::DataFrame)
    t = [ctc DataFrame(i=1:nrow(ctc))]
    t = by(t, names(ctc),g->DataFrame(n=nrow(g), i=[g[:,:i]]))
    sort!(t.n); sort!(t[1:end-1]); return t
end

"""
# Normalize values of an array to be between -1 and 1
# original sign of the array values is maintained.
"""
function normalized(data;equal=true)

    if isequal(equal,true)
        if abs(min(data...)) > max(data...)
            max_range_value = abs(min(data...))
            min_range_value = min(data...);
        else
            max_range_value = max(data...)
            min_range_value = -max(data...)
        end
        norm_value = 2 .* data / (max_range_value - min_range_value)
    elseif isequal(equal,false)
        max_range_value = max(data...)
        min_range_value = abs(min(data...))
        b=zeros(size(data))
        bw = data.> 0
        idx = findall(x->x==1,bw)
        b[idx] = data[idx] / max_range_value
        bw = data.< 0
        idx = findall(x->x==1,bw)
        b[idx] = data[idx] / min_range_value
        norm_value = b
    end
    return norm_value
end

"""
# Find the extreme value
"""
function sbxexd(stas,stamag,goodTimeidx;sig=true)
    bstdi = argmax(stamag)
    exi = [argmin(stas[:,:,bstdi]),argmax(stas[:,:,bstdi])]
    ex = stas[exi,bstdi]
    absexi = argmax(abs.(ex))
    (ex=sig ? ex[absexi] : 0, mag=stamag[bstdi], pd=goodTimeidx[bstdi])
end

sfdog(x,p...) = dogf(x;aₑ=p[1],μₑ=p[2],σₑ=p[3],aᵢ=p[4],μᵢ=p[5],σᵢ=p[6])

"""
# Peichao:
# Join the results from differernt Hartleys after hartleyFourier analysis
"""
function sbxjoinhartleyFourier(harts)
    cn = length(harts)
    uids = length(harts[1]["result"].cellId)
    kernsize = size(harts[1]["result"].kernnor[1])
    dataset = Dict()
    signifs=Dict();taumaxs=Dict();kstdmaxs=Dict();kdeltas=Dict();slambdas=Dict();
    orimaxs=Dict();sfmaxs=Dict();orimeans=Dict();sfmeans=Dict();
    kerns=Dict();kernraws=Dict();kernests=Dict();oris=Dict();oriLevels=Dict();
    orifits=Dict();sfLevels=Dict(); sffits=Dict();sfs=Dict();

    for u in 1:uids
        kern = Array{Float64}(undef,kernsize...,cn)  # Normlized kernels of all hartley
        kernraw = Array{Float64}(undef,kernsize...,cn)  # Raw kernels of all hartley
        kernest = Array{Float64}(undef,kernsize...,cn)  # Estimated kernels of all hartley
        signif = zeros(1,cn);taumax = zeros(1,cn);kstdmax = zeros(1,cn);
        kdelta = zeros(1,cn);slambda = zeros(1,cn);orimax = zeros(1,cn);
        sfmax = zeros(1,cn);orimean = zeros(1,cn);sfmean = zeros(1,cn);
        ori=[]; oriLevel=[]; sf=[]; sfLevel=[]; orifit=[]; sffit=[];

        for c in 1:cn
            signif[1,c] = harts[c]["result"].signif[u]
            taumax[1,c] = harts[c]["result"].taumax[u]
            kstdmax[1,c] = harts[c]["result"].kstdmax[u]
            kdelta[1,c] = harts[c]["result"].kdelta[u]
            slambda[1,c] = harts[c]["result"].slambda[u]
            orimax[1,c] = harts[c]["result"].orimax[u]
            sfmax[1,c] = harts[c]["result"].sfmax[u]
            orimean[1,c] = harts[c]["result"].orimean[u]
            sfmean[1,c] = harts[c]["result"].sfmean[u]

            kern[:,:,c] = harts[c]["result"].kernnor[u]
            kernraw[:,:,c] = harts[c]["result"].kernraw[u]
            kernest[:,:,c] = harts[c]["result"].kernest[u]

            θ = deg2rad.(harts[c]["result"].oriLevel[u])
            fr = harts[c]["result"].oriResp[u]

            push!(oriLevel,θ);
            push!(ori,fr);
            # fit von Mises for orientation
            fit=()
            try
                # vmfit = curve_fit((x,p)->vmf.(x,p...,n=2),θ,fr,[1.0,0,1])
                # if vmfit.converged
                #     x = 0:0.004:2π
                #     y = vmf.(x,vmfit.param...,n=2)
                #     fit = (fit...,circtuningstats(x,y,od=0.5π,s=:o)...,vm=vmfit)
                # end
                mfit = fitmodel(:vmn2,θ,fr)
                fit = (fit...,circtuningfeature(mfit,od=0.5π,fn=:o)...,vmn2=mfit)
            catch
            end
            push!(orifit,fit)

            fl = harts[c]["result"].sfLevel[u]
            fr = harts[c]["result"].sfResp[u]

            push!(sfLevel,fl);
            push!(sf,fr);
            # fit difference of gaussians
            fit=()
            try
                # dogfit = curve_fit((x,p)->sfdog.(x,p...),fl,fr,[1.0,0,1,1,0,1])
                # if dogfit.converged
                #     fit = (dog=dogfit)
                # end
                mfit = fitmodel(:dog,fl,fr)
                fit = (sftuningfeature(mfit)...,dog=mfit)
            catch
            end
            push!(sffit,fit)

        end
        signifs[u] = signif
        taumaxs[u] = taumax
        kstdmaxs[u] = kstdmax
        kdeltas[u] = kdelta
        slambdas[u] = slambda
        orimaxs[u] = orimax
        sfmaxs[u] = sfmax
        orimeans[u] = orimean
        sfmeans[u] = sfmean

        kerns[u] = kern
        kernraws[u] = kernraw
        kernests[u] = kernest

        oriLevels[u] = oriLevel
        oris[u] = ori
        orifits[u] = orifit
        sfLevels[u] = sfLevel
        sfs[u] = sf
        sffits[u] = sffit
    end
    dataset["signif"] = signifs
    dataset["taumax"] = taumaxs
    dataset["kstdmax"] = kstdmaxs
    dataset["kdelta"] = kdeltas
    dataset["slambda"] = slambdas
    dataset["orimax"] = orimaxs
    dataset["sfmax"] = sfmaxs
    dataset["orimean"] = orimeans
    dataset["sfmean"] = sfmeans

    dataset["kern"] = kerns
    dataset["kernraw"] = kernraws
    dataset["kernest"] = kernests

    dataset["oriraw"] = oris
    dataset["oriLevel"] = oriLevels
    dataset["orifit"] = orifits
    dataset["sfraw"] = sfs
    dataset["sfLevel"] = sfLevels
    dataset["sffit"] = sffits

    return dataset
end

function sbxjoinsta(stas,lbTime,ubTime,blkTime)
    cn = length(stas)
    imagesize = stas[1]["imagesize"]
    maskradius = stas[1]["maskradius"]
    delays = stas[1]["delays"]
    stisize = stas[1]["stisize"]
    ppd = imagesize[1]/stisize[1]
    xi = stas[1]["xi"]     # pixels in stimulated region
    nxi = setdiff(1:prod(imagesize),xi)  # pxiels in background
    cii=CartesianIndices(imagesize)

    goodTimeidx=findall(x -> lbTime .<= x .<= ubTime, delays)
    bdi=findall(x -> (x .<=blkTime), delays)   # select as blanks

    dataset = Dict("stisize"=>stisize,"imagesize"=>imagesize,"maskradius"=>maskradius,"delays"=>delays,"ppd"=>ppd,"bdi"=>bdi,"xi"=>xi)
    # dataset["log"] = [i["log"] for i in stas]
    dataset["color"]=[i["color"] for i in stas]
    # dataset["minmaxcolor"] = [(i["mincolor"],i["maxcolor"]) for i in stas]
    # dataset["minmaxcg"] = map(i->minmaxcolorgradient(i...,n=256),dataset["minmaxcolor"])
    rawsta = Dict();usta = Dict();uzsta = Dict();ustamag = Dict();ustavec = Dict(); ucex=Dict();uresponsive=Dict()
    # uids = mapreduce(c->size(stas[c]["usta"],1),intersect,1:cn)
    uids = size(stas[1]["usta"],1)
    # uids = mapreduce(c->keys(stas[c]["usta"]),intersect,1:cn)
    for u in 1:uids
        rsta = Array{Float64}(undef,imagesize...,length(delays),cn)  # raw sta
        csta = Array{Float64}(undef,imagesize...,length(delays),cn)  # unit sta
        zsta = Array{Float64}(undef,imagesize...,length(delays),cn)  # z-scored sta
        # bsta = Array{Float64}(undef,imagesize...,length(bdi))  # blank sta
        cstamag = Array{Float64}(undef,length(delays),cn)
        cstavec = Array{Float64}(undef,length(xi),length(delays),cn)

        cex = []
        for c in 1:cn
            bsta = stas[c]["usta"][u,bdi,:]
            bm = mean(bsta);bsd = std(bsta)
            zzsta = zscore(stas[c]["usta"][u,:,:],bm,bsd)
            for d in eachindex(delays)
                rawstavec = stas[c]["usta"][u,d,:]  # raw sta of grating stim
                mag = norm(rawstavec)  # magnitude of the raw sta vector
                ucsta = rawstavec/mag   # unit vector sta
                cstamag[d,c] = mag
                csta[cii[nxi],d,c].= 0  # NonS-stimulated region, background(0 0 0)
                csta[cii[xi],d,c] = ucsta   # unit vector STA; Stimulated region,RF  # now STA is 2D matrix, not a vector
                cstavec[:,d,c] = ucsta
                rsta[cii[nxi],d,c].= 0
                rsta[cii[xi],d,c] = rawstavec  # raw vector sta
                zsta[cii[nxi],d,c].= 0
                zsta[cii[xi],d,c] = zzsta[d,:]  # z-scored sta
            end
            push!(cex,sbxexd(csta[:,:,goodTimeidx,c], cstamag[goodTimeidx,c],goodTimeidx))   #PL: To check the extrema of STA, I set the time window constrain here
        end
        rawsta[u] = rsta
        usta[u] = csta
        uzsta[u] = zsta
        ustamag[u] = cstamag
        ustavec[u] = cstavec
        ucex[u] = cex
        uresponsive[u] = false
    end
    dataset["rawsta"] = rawsta
    dataset["usta"] = usta
    dataset["zsta"] = uzsta
    dataset["ustamag"] = ustamag
    dataset["ustavec"] = ustavec
    dataset["ucex"] = ucex
    dataset["uresponsive"] = uresponsive
    return dataset
end

"peak ROI region and its delay"
function sbxpeakroi(clc)
    pi = [Tuple(argmax(clc))...]  # peak index
    plc = clc[:,:,pi[3:end]...]   # peak local contrast at best delay
    segs = seeded_region_growing(plc,[(CartesianIndex(1,1),1),(CartesianIndex(pi[1:2]...),2)])  # find ROI in STA image at best delay
    idx = findall(labels_map(segs).==2)   # find pixels in ROI
    idxlims = dropdims(extrema(mapreduce(i->[Tuple(i)...],hcat,idx),dims=2),dims=2)
    hw = map(i->i[2]-i[1],idxlims)
    center = round.(Int,mean.(idxlims))
    radius = round(Int,maximum(hw)/2)
    return (i=idx,center=center,radius=radius,pd=pi[3])
end
"check local mean or contrast in a region at a delay significently higher than baseline"
# function sbxisresponsive(sta,idx,d,bdi,badTimeidx;msdfactor=3.0,csdfactor=3.0)
#     d in badTimeidx && return false
#     blm = [mean(sta[idx,j]) for j in bdi]
#     bmm = mean(blm);bmsd=std(blm)
#     blc = [std(sta[idx,j]) for j in bdi]
#     bcm = mean(blc);bcsd=std(blc)
#     (std(sta[idx,d]) > bcm+csdfactor*bcsd) || (mean(sta[idx,d]) > bmm+msdfactor*bmsd)
# end

function sbxisresponsive(sta,idx,bdi,badTimeidx;mfactor=3.0,cfactor=3.0)
    lc = [std(sta[idx,j]) for j in 1:size(sta,3)]   # local variance of ROI for each delay
    lm = [mean(sta[idx,j]) for j in 1:size(sta,3)]  # local mean of ROI for each delay
    lcmaxd = argmax(lc); lcmax = lc[lcmaxd]  # maxi local variance
    lmmaxd = argmax(lm); lmmax = lm[lmmaxd]  # maxi local mean
    lmmind = argmin(lm); lmmin = lm[lmmind]  # min local mean
    bmm = mean(lm[bdi]);bmsd=std(lm[bdi])  # mean and std of local mean of blank condition
    bcm = mean(lc[bdi]);bcsd=std(lc[bdi])  # mean and std of local variance of blank condition

    (!(lcmaxd in badTimeidx) && lcmax > bcm+cfactor*bcsd) ||
    (!(lmmaxd in badTimeidx) && lmmax > bmm+mfactor*bmsd) ||
    (!(lmmind in badTimeidx) && lmmin < bmm-mfactor*bmsd)
end

function sbxresponsivesta!(dataset,lbTime,ubTime,respThres;ws=0.5,msdfactor=3.5,csdfactor=3.5)
    imagesize = dataset["imagesize"]
    maskradius = dataset["maskradius"]
    ppd=dataset["ppd"]
    xi = dataset["xi"]     # pixels in stimulated region
    bdi = dataset["bdi"]
    usta = dataset["usta"]
    uzsta = dataset["zsta"]
    ustamag = dataset["ustamag"]
    ucex = dataset["ucex"]
    delays = dataset["delays"]
    uresponsive=dataset["uresponsive"]

    center = [imagesize[1]÷2,imagesize[2]÷2]
    radius = Int64(ceil(max(imagesize[1],imagesize[2])*maskradius))
    goodTimeidx=findall(x -> lbTime .<= x .<= ubTime, delays)
    badTimeidx=findall(x -> (x .< lbTime || x .> ubTime), delays)

    ulsta=Dict();ulroi=Dict();ulcex=Dict();uconeresponsive=Dict();
    p = ProgressMeter.Progress(length(usta),desc="Test STAs ... ")
    for u in keys(usta)
        # u=671
        clc = localcontrast(usta[u],round(Int,ws*ppd))
        # cproi = map(c->sbxpeakroi(clc[:,:,:,c]),1:size(clc,4))
        cproi = map(c->sbxpeakroi(clc[:,:,:,c]),1:size(clc,4))
        # idx,c,cpd = sbxpeakroi(clc)

        # PL: Check responsiveness of each cone type and achromatic
        # PL: based on std and mean of STA image (compare with blank, which are STAs from the stim onset or before that), and delay range
        rs= [sbxisresponsive(uzsta[u][:,:,:,c],cproi[c].i,bdi,badTimeidx,mfactor=msdfactor,cfactor=csdfactor) for c in 1:size(clc,4)]
        # PL: based on threshold and delay range, to filter out low-response and cell 'response' too early or too late
        # rs2 = [((ustamag[u][cproi[i].pd,i]-mean(ustamag[u][bdi,i]))/mean(ustamag[u][bdi,i])>respThres) && (ucex[u][i].pd in goodTimeidx) for i in 1:size(clc,4)]
        rs2 = [(ustamag[u][ucex[u][i].pd,i]-mean(ustamag[u][bdi,i]))/mean(ustamag[u][bdi,i])>respThres  for i in 1:size(clc,4)]
        # rs2 = [(abs(ucex[u][i].ex)>respThres) && (ucex[u][i].pd in goodTimeidx) for i in 1:size(clc,4)]

        uconeresponsive[u] = rs .& rs2  # save cell's responesiveness to each cone & achromatic stimuli

        # PL: if cell response to one of cone stim and pass the threshold, it is responsive.
        if any(rs .& rs2)
            uresponsive[u] = true
            # vr = map(i->intersect(i.+(-radius:radius),1:imagesize[1]),center)
            # radius = (minimum ∘ mapreduce)((r,c)->abs.([r[1],r[end]].-c),append!,vr,center)
            # ulsta[u] = usta[u][xi]
            ulsta[u] = usta[u][map(i->i.+(-radius:radius),center)...,:,:]  # PL: I cut the STA image according to stim radius, not radius detected by localcontrast
            # ulroi[u] = (center=center,stisize=(2*radius+1)/ppd,d=cproi[i].pd,c=c)
            ulcex[u]=map(i->sbxexd(ulsta[u][:,:,goodTimeidx,i],ustamag[u][goodTimeidx,i],goodTimeidx,sig=uconeresponsive[u][i]),1:size(ulsta[u],4))
        else
            uresponsive[u] = false
        end
        next!(p)
    end
    dataset["ulsta"]=ulsta
    # dataset["ulroi"]=ulroi
    dataset["ulcex"]=ulcex
    dataset["uconeresponsive"] = uconeresponsive
    return dataset
end

# rfgabor(x,y,p...) = gaborf(x,y,a=p[1],μ₁=p[2],σ₁=p[3],μ₂=p[4],σ₂=p[3]*p[5],θ=p[6],f=p[7],phase=p[8])
# # rfdog(x,y,p...) = dogf(x,y,aₑ=p[1],μₑ₁=p[2],σₑ₁=p[3],μₑ₂=p[4],σₑ₂=p[3]*p[5],θₑ=p[6],aᵢ=p[7],μᵢ₁=p[2]+p[8],σᵢ₁=p[3]*p[9],μᵢ₂=p[4]+p[10],σᵢ₂=p[3]*p[9]*p[5],θᵢ=p[6])
# rfdog(x,y,p...) = dogf(x,y,aₑ=p[1],μₑ₁=p[2],σₑ₁=p[3],μₑ₂=p[4],σₑ₂=p[3],θₑ=0,aᵢ=p[5],μᵢ₁=p[2],σᵢ₁=p[3]*p[6],μᵢ₂=p[4],σᵢ₂=p[3]*p[6],θᵢ=0)

# function sbxmodelfit(data,ppd;model=:gabor)
#     alb,aub = abs.(extrema(data))
#     ab = max(alb,aub)
#     rpx = (size(data)[1]-1)/2
#     r = rpx/ppd

#     x = (mapreduce(i->[i[2] -i[1]],vcat,CartesianIndices(data)) .+ [-(rpx+1) rpx+1])/ppd
#     y = vec(data)
#     rlt = mfun = missing
#     try
#         if model == :dog
#             if aub >= alb
#                 ai = 3.5alb
#                 ae = aub + ai
#                 es = 0.2r;esl=0.15r;esu=0.3r
#                 ier=2;ierl = 1.1;ieru = 3
#             else
#                 ae = 3.5aub
#                 ai = alb + ae
#                 es = 0.4r;esl=0.16r;esu=0.6r
#                 ier=0.5;ierl = 0.3;ieru = 0.9
#             end
#             # lb=[0,          -0.4sr,    0.1sr,   -0.4sr,    0.5,    0,     0,       -0.1sr,     0.1,    -0.1sr]
#             # ub=[10,         0.4sr,    0.5sr,    0.4sr,    2,      π,     Inf,      0.1sr,     10,       0.1sr]
#             # p0=[0,       0,        0.3sr,    0,        1,      π/4,   aei[2],    0,         0.25,       0]
#             ub=[1.5ae,    0.6r,    esu,    0.6r,     1.5ai,    ieru]
#             lb=[0.2ae,   -0.6r,    esl,   -0.6r,     0.5ai,    ierl]
#             p0=[ae,       0,        es,     0,          ai,    ier]
#             mfun = (x,p) -> rfdog.(x[:,1],x[:,2],p...)
#         elseif model == :gabor
#             ori,sf = freqimagestats(powerspectrum2(data,ppd)...)
#             fub = min(1.5sf,8);flb=max(0.5sf,0.2)
#             ub=[1.5ab,   0.6r,   1.0r,   0.6r,   2.3r,    π,    fub,   1]
#             lb=[0.2ab,  -0.6r,   0.1r,  -0.6r,   0.3r,    0,    flb,   0]
#             p0=[ab,      0,      0.2r,    0,     0.2r,  ori,     sf,   0]
#             mfun = (x,p) -> rfgabor.(x[:,1],x[:,2],p...)
#         end
#         if !ismissing(mfun)
#             mfit = curve_fit(mfun,x,y,p0,lower=lb,upper=ub,
#             maxIter=3000,x_tol=1e-11,g_tol=1e-15,min_step_quality=1e-4,good_step_quality=0.25,lambda_increase=5,lambda_decrease=0.2)
#             rlt = (model=model,radius=r,param=mfit.param,converged=mfit.converged,resid=mfit.resid,r=cor(y,mfun(x,mfit.param)))
#         end
#     catch
#     end
#     return rlt
# end
"fit responsive sta to each type of models"
# function sbxrffit!(dataset;model=[:gabor,:dog])
#     ulsta = dataset["ulsta"]
#     ulroi = dataset["ulroi"]
#     ppd = dataset["ppd"]
#     uconeresponsive = dataset["uconeresponsive"]

#     if !haskey(dataset,"urf")
#         dataset["urf"]=Dict()
#     end
#     urf = dataset["urf"]
#     p = ProgressMeter.Progress(length(ulsta),desc="Fit RFs ... ")
#     for u in keys(ulsta)
#         if !haskey(urf,u)
#             urf[u] = Dict()
#         end
#         exs=map(i->i.ex,dataset["ucex"][u])
#         ds=map(i->i.d,dataset["ucex"][u])
#         for m in model
#             urf[u][m] = map(i->isequal(uconeresponsive[u][i],false) ? nothing : sbxmodelfit(ulsta[u][:,:,ds[i],i],ppd,model=m),1:length(exs))
#         end
#         next!(p)
#     end
#     return dataset
# end

function sbxfitsta!(dataset;model=[:gabor,:dog])
    ulsta = dataset["ulsta"]
    ppd = dataset["ppd"]
    uconeresponsive = dataset["uconeresponsive"]

    if !haskey(dataset,"ulfit")
        dataset["ulfit"]=Dict()
    end
    ulfit = dataset["ulfit"]
    p = ProgressMeter.Progress(length(ulsta),desc="Fit STAs ... ")
    for u in keys(ulsta)
        if !haskey(ulfit,u)
            ulfit[u] = Dict()
        end
        exs=map(i->i.ex,dataset["ulcex"][u])
        ds=map(i->i.pd,dataset["ulcex"][u])
        for m in model
            ulfit[u][m] = map(i->isequal(uconeresponsive[u][i],false) ? missing : fitmodel2(m,ulsta[u][:,:,ds[i],i],ppd), 1:length(exs))
        end
        next!(p)
    end
    return dataset
end

function sbxgetbestconesta(dataset)
    uconeresponsive = dataset["uconeresponsive"]
    delays = dataset["delays"]

    ublsta=Dict(); ubusta=Dict(); ubrawsta=Dict(); ubstamag=Dict(); ubcone=Dict(); ubmag=Dict();ustasign=Dict(); uachro=[];cellId=[];achroResp=DataFrame(); # save Cone Weights
    for u in keys(dataset["ulsta"])
        # check on/off dominance of a cell for each cone (inculde of achromatic), not restrict the delay is the same.
        # if any(uconeresponsive[u][1:4])
        exs=map(i->i.ex,dataset["ulcex"][u][1:4])
            # a = uconeresponsive[u]
            # b = convert(Array{Float64},a)
            # c = sign.(exs)
            # ustasign[u] = b.*c
        ustasign[u] = sign.(exs)
        # end
        # for calculating cone weight, so the time delay is restricted to the same delay.
        if any(uconeresponsive[u][1:3])
            exs=map(i->i.ex,dataset["ulcex"][u][1:3])
            mags=map(i->i.mag,dataset["ulcex"][u][1:3])
            ds=map(i->i.pd,dataset["ulcex"][u][1:3])

            # bestconeidx = findmax(abs.(exs))[2]  # define by max exs(extrem value)
            bestconeidx = findmax(mags)[2]  # define dominant cone by magnitude.
            bestdly = ds[bestconeidx]  # the best deley is just for best cone, not best deley for each cone.

            ulsta = dataset["ulsta"][u]
            usta = dataset["usta"][u]
            urawsta = dataset["rawsta"][u]
            ustamag = dataset["ustamag"][u]
            ublsta[u] = ulsta[:,:,bestdly,:]        # for ploting
            ubusta[u] = usta[:,:,bestdly,:]        # for cone weight calculation
            ubrawsta[u] = urawsta[:,:,bestdly,:]   # for cone weight calculation
            ubstamag[u] = ustamag[bestdly,:]
            ubcone[u] = (resp=sign.(exs), bstidx=bestconeidx, bstdly=delays[bestdly],bstex=exs[bestconeidx],bstmag=mags[bestconeidx])
        end
        if uconeresponsive[u][4]
            exs=map(i->i.ex,dataset["ulcex"][u][4:4])
            push!(uachro, exs[1]);push!(cellId,u)
        end
    end
    cellNum = length(cellId)
    achroResp.py = 0:cellNum-1
    achroResp.cellId=cellId
    achroResp.uachro=uachro
    dataset["ublsta"] = ublsta
    dataset["ubusta"] = ubusta
    dataset["ubrawsta"] = ubrawsta
    dataset["ubstamag"] = ubstamag
    dataset["ubcone"] = ubcone
    dataset["ustasign"] = ustasign
    dataset["achroResp"] = achroResp
    return dataset
end

function coneweight(lcone,mcone,scone)
    tcw = abs(lcone)+abs(mcone)+abs(scone)
    lcw = lcone/tcw
    mcw = mcone/tcw
    scw = scone/tcw
    return lcw,mcw,scw
end

function findextreme(a)
    b=findmax(a)
    c=findmin(a)
    if abs(b[1]) >= abs(c[1])
        d=a[b[2]]
    elseif abs(b[1]) < abs(c[1])
        d=a[c[2]]
    end
    return d
end

function sbxgetconeweight(dataset)
    ubrawsta = dataset["ubrawsta"]
    ubstamag = dataset["ubstamag"]
    ubusta = dataset["ubusta"]
    ubcone = dataset["ubcone"]

    cw = DataFrame()
    cellId=[];Lconemn=[];Mconemn=[];Sconemn=[];Lconemg=[];Mconemg=[];Sconemg=[];Lconemgall=[];Mconemgall=[];Sconemgall=[];

    for u in keys(ubcone)
        ## make STA mask
        if ubcone[u].bstex>=0
            img = ubusta[u][:,:,ubcone[u].bstidx]
        elseif ubcone[u].bstex<0
            img = -ubusta[u][:,:,ubcone[u].bstidx]
        end
        # smoothing
        imgf = imfilter(img, Kernel.gaussian(3))
        mask = Gray.(imgf).> max(imgf...)*0.95
        # calculate relative cone weights
        # coneresp=ubcone[u].resp[1:3]

        lconemn = sum(ubusta[u][:,:,1] .* mask)/sum(mask) #* abs(coneresp[1])
        mconemn = sum(ubusta[u][:,:,2] .* mask)/sum(mask) #* abs(coneresp[2])
        sconemn = sum(ubusta[u][:,:,3] .* mask)/sum(mask) #* abs(coneresp[3])
        lcwmn,mcwmn,scwmn = coneweight(lconemn,mconemn,sconemn)

        idx=findall(mask.==1)
        lconemg = norm(ubrawsta[u][:,:,1][idx]) * sign(findextreme(ubrawsta[u][:,:,1][idx]))
        mconemg = norm(ubrawsta[u][:,:,2][idx]) * sign(findextreme(ubrawsta[u][:,:,2][idx]))
        sconemg = norm(ubrawsta[u][:,:,3][idx]) * sign(findextreme(ubrawsta[u][:,:,3][idx]))
        lcwmg,mcwmg,scwmg = coneweight(lconemg,mconemg,sconemg)

        lconemgall = ubstamag[u][1] * sign(findextreme(ubrawsta[u][:,:,1][idx]))
        mconemgall = ubstamag[u][2] * sign(findextreme(ubrawsta[u][:,:,2][idx]))
        sconemgall = ubstamag[u][3] * sign(findextreme(ubrawsta[u][:,:,3][idx]))
        lcwmgall,mcwmgall,scwmgall = coneweight(lconemgall,mconemgall,sconemgall)

        push!(cellId, u); push!(Lconemn,lcwmn);push!(Mconemn,mcwmn);push!(Sconemn,scwmn);
        push!(Lconemg,lcwmg);push!(Mconemg,mcwmg);push!(Sconemg,scwmg);
        push!(Lconemgall,lcwmgall);push!(Mconemgall,mcwmgall);push!(Sconemgall,scwmgall);
    end
    cellNum = length(cellId)
    cw.py = 0:cellNum-1
    cw.cellId=cellId
    cw.Lconemn=Lconemn
    cw.Mconemn=Mconemn
    cw.Sconemn=Sconemn
    cw.Lconemg=Lconemg
    cw.Mconemg=Mconemg
    cw.Sconemg=Sconemg
    cw.Lconemgall=Lconemgall
    cw.Mconemgall=Mconemgall
    cw.Sconemgall=Sconemgall
    dataset["coneweight"]=cw

    return dataset
end


text:
mutable struct Token
    token::String
    expires::DateTime
    scope::String
    refreshtoken::String
    maxtries::Int
end

"""
    Token(refreshtoken::String)

Returns a `Token` that can be used for requests.

# Arguments:
- `refreshtoken` : Refreshtoken associated with your Chartmetric account.
"""
function Token(refreshtoken::String,  maxtries::Int = 5)
    @assert maxtries > 0 "maxtries have to be positive"
    token_dict, expires = tokenrequest(refreshtoken, maxtries)
    return Token(token_dict["token"],
                expires,
                token_dict["scope"],
                token_dict["refresh_token"],
                maxtries)
end

"""
    newtoken!(token::Token; verbose = false)

Refreshes an existing request `Token`

# Arguments:
- `token` : the request `Token`
- `verbose=false` : print a message when a new token is requested
"""
function newtoken!(token::Token; verbose = false)
    verbose && println("Requesting new Token...")
    token_dict, expires = tokenrequest(token.refreshtoken, token.maxtries)
    token.token = token_dict["token"]
    token.expires = expires
    token.scope = token_dict["scope"]
    return token
end

function tokenprotector!(token::Token)
    token.expires <= (Dates.now() + Dates.Minute(1)) && newtoken!(token)
    return token
end

function tokenrequest(refreshtoken::String, maxtries = 5)
    header = ["Content-Type" => "application/json"]
    body = """{"refreshtoken":"$refreshtoken"}"""
    time = Dates.now()
    req = :(HTTP.request("POST",
            "https://api.chartmetric.com/api/token",
            $header, $body ,
            retry = true,
            status_exception = false))
    resp = eval(req)
    if resp.status ∉ 200:399
        for trial in 1:(maxtries-1)
            resp = eval(req)
            resp.status ∈ 200:399 && break
            trial >= (maxtries - 2) && println(resp); sleep(60)
            trial >= (maxtries - 1) && @error "Could not get token"
        end
    end
    resp.status > 500 && return nothing
    token_dict = JSON.parse(String(resp.body))
    expires = time + Dates.Second(token_dict["expires_in"])
    return (token_dict, expires)
end

ext: jl


text:
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, D},
    ::Type{S},
) where {
    H <: HydroEnergyCascade,
    D <: HydroDispatchReservoirCascade,
    S <: PM.AbstractActivePowerModel,
}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    #Variables
    PSI.add_variables!(optimization_container, PSI.ActivePowerVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.EnergyVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.SpillageVariable, devices, D())

    #Initial Conditions
    PSI.storage_energy_init(optimization_container, devices)

    #Constraints
    energy_balance_cascade_constraint!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    #Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end

function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, D},
    ::Type{S},
) where {
    H <: HydroEnergyCascade,
    D <: HydroDispatchReservoirCascade,
    S <: PM.AbstractPowerModel,
}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    #Variables
    PSI.add_variables!(optimization_container, PSI.ActivePowerVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.EnergyVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.SpillageVariable, devices, D())
    PSI.add_variables!(optimization_container, PSI.ReactivePowerVariable, devices, D())

    #Initial Conditions
    PSI.storage_energy_init(optimization_container, devices)

    #Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ReactivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    energy_balance_cascade_constraint!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    #Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end


text:
#@ Model: Stochastic Unit Commitment
#@ Description: Adapted from model of Patricia Levi, pjlevi@stanford.edu . Find the dispatch schedule, p, for generators g, that minimizes costs across time t, and scenarios o, subject to technical constraints, a sample of which are below. $$\min\sum_{g,t,o}(start\_cost_{g,t,o}+p_{g,t,o} varcost_g)$$$$\sum_g p_{g,t,o}={dem}_t\quad\forall t,o$$$$start\_cost_{g,t,o}\geq v_{g,t,o}startup_g\forall g,t,o$$$$pmin_g u_{g,t,o} \leq p_{g,t,o} \leq pmax_g u_{g,t,o}$$$$v_{g,t+1,o}=u_{g,t+1,o}-u_{g,t,o}$$$$u,v\in(0,1),p\in\mathbb{R}$$

# Original source code available at https://github.com/pajalevi/Stochastic_UC_model_julia

using JuMP, Clp, JSON

# PARAMETERS ###
#@ Helper object: n_t
#@ Description: number of timesteps
n_t = 10

#@ Input Object: n_g
#@ Description: number of generators
n_g = 10

#@ Input Object: n_gsl
#@ Description: number of slow generators
n_gsl = 5

#@ Helper object: n_gf
#@ Description: number of fast generators
n_gf = n_g - n_gsl

#@ Helper object: n_gdr
#@ Description: number of DR generators
n_gdr = 1

#@ Helper Object: n_omega
#@ Description: number of realizations
n_omega=3

#@ Helper object: t_firsts
#@ Description: index of timesteps that are 'first' of a timeblock
t_firsts=1

# Labels
slowgen_timesteps = Dict()
slowgen_tstart = Dict()
dr_timesteps = Dict()
for i in 1:n_gdr
    for j in 1:n_t

    end
end

for i in 1:n_gsl
    for j in 1:n_t
        slowgen_timesteps[(i,j)] = "gen$(i)_t$(j)"
        if j < n_t
            slowgen_tstart[(i,j)] = "gen$(i)_t$(j)"
        end
        if i <= n_gdr
            dr_timesteps[(i,j)] = "gen$(i)_t$(j)"
        end
    end
end

fastgen_tstart_scen = Dict()
for i in (n_gsl+1):n_g
    for j in 1:(n_t-1)
        for k in 1:n_omega
            fastgen_tstart_scen[(i,j,k)] = "gen$(i)_t$(j)_scen$(k)"
        end
    end
end

gen_t_scen = Dict()
for i in 1:n_g
    for j in 1:n_t
        for k in 1:n_omega
            gen_t_scen[(i,j,k)] = "gen$(i)_t$(j)_scen$(k)"
        end
    end
end

gen_nodr_t_scen = Dict()
for i in 1:n_g-n_gdr
    for j in 1:n_t
        for k in 1:n_omega
            gen_nodr_t_scen[(i,j,k)] = "gen$(i+1)_t$(j)_scen$(k)"
        end
    end
end

t_scen = Dict()
for j in 1:n_t
    for k in 1:n_omega
        t_scen[(j,k)] = "t$(j)_scen$(k)"
    end
end

### SETS ###
#@ Helper object: TIME
TIME = 1:n_t

#@ Helper object: SCENARIOS
SCENARIOS = 1:n_omega

#@ Helper object: GENERATORS
#@ Description: all generators
GENERATORS = 1:n_g

#@ Helper object: GEN_NODR
#@ Description: all generators excluding DR
GEN_NODR = (n_gdr + 1):n_g

#@ Helper object: GF
#@ Description: index of fast generators
GF = n_gsl+1:n_g

#@ Helper object: GSL
#@ Description: index of slow generators
GSL = 1:n_gsl

#@ Helper object: GDR
#@ Description: index of DR generator
GDR = 1

### STOCHASTIC PARAMS ###
#@ Helper object: vdr
#@ Description: stochastic parameter vdr
vdr = [0.9,1,1.1]

#@ Helper object: pro
#@ Description: stochastic parameter pro
pro = [0.25,0.5,0.25]
#TODO: insert a check that vdr, pro are same length as n_omega

### INPUT DATA ###
#@ Helper object: dem
#@ Description: demand
dem = [collect(1:(n_t/2)); collect((n_t/2):-1:1)] * 18

#generator min and max
#@ Helper object: pmin
#@ Description: minimum generator
pmin = repeat([1],inner = n_g)

#@ Helper object: pmax
#@ Description: maximum generator
pmax = repeat([10],outer = n_g)

#@ Helper object: startup
#@ Description: startup costs
startup = [0;repeat([2],inner = n_gsl -n_gdr);repeat([1], inner = n_gf)]

#@ Helper object: varcost
#@ Description: variable costs
varcost = [0;collect(1:0.1:(1+0.1*(n_gsl-n_gdr)));collect(2:0.1:(2+0.1*(n_gf)))]

#@ Helper object: pf
#@ Description: generator availability
pf = repeat([1], inner = [n_g, n_t])

### Model ###
#@ Solver: solver
solver = Clp.Optimizer

#@ Problem: m
m = Model(with_optimizer(solver))

### VARIABLES ###
#@ Variable: z
#@ Description: slow generator startup
#@ Labels: slowgen_timesteps
@variable(m, z[1:n_gsl,1:n_t])

#@ Variable: w
#@ Description: slow generator commitment (relaxed binary)
#@ Labels: slowgen_timesteps
@variable(m, 0 <= w[1:n_gsl,1:n_t] <= 1)

#@ Variable: v
#@ Description: real-time commitment, startup
#@ Labels: gen_t_scen
@variable(m, v[1:n_g,1:n_t,1:n_omega])

#@ Variable: u
#@ Description: generator commitment RELAXED BINARY
#@ Labels: gen_t_scen
@variable(m, 0 <= u[1:n_g,1:n_t,1:n_omega] <= 1)

#@ Variable: p
#@ Description: generator production
#@ Labels: gen_t_scen
@variable(m, p[1:n_g,1:n_t,1:n_omega] >= 0)

#@ Variable: p_dr
#@ Description: day-ahead production commitment
#@ Labels: dr_timesteps
@variable(m, p_dr[1:n_gdr, 1:n_t] >= 0)

#@ Variable: start_cost
#@ Description: start_cost
#@ Labels: gen_t_scen
@variable(m, start_cost[1:n_g, 1:n_t, 1:n_omega] >= 0)

### CONSTRAINTS ###

#@ Constraint: SD
#@ Description: Supply from the sum of generators must equal demand $$\sum_g P_{g,t,o}={dem}_t\quad\forall t,o$$
#@ Labels: t_scen
@constraint(m,SD[t=1:n_t, o=1:n_omega],
    sum( p[g,t,o] for g=1:n_g) == dem[t])

#@ Constraint: GENMIN
#@ Description: Minimum generation $$p_{g,t,o}\geq pmin_g u_{g,t,o}$$
#@ Labels: gen_t_scen
@constraint(m, GENMIN[g= 1:n_g, t= 1:n_t, o=1:n_omega ],
    p[g,t,o] >= pmin[g] * u[g,t,o])

#@ Constraint: GENMAX
#@ Description: Maximum generation $$p_{g,t,o)\leq pmax_g u_{g,t,o} pf_{g,t}$$
#@ Labels: gen_nodr_t_scen
@constraint(m,GENMAX[g = GEN_NODR, t = 1:n_t, o=1:n_omega],
    p[g,t,o] <= pmax[g] * u[g,t,o] * pf[g,t] )

#@ Constraint: GENMAXDR
#@ Description: Maximum generation, demand resonse
#@ Labels: dr_timesteps
@constraint(m,GENMAXDR[g = 1:n_gdr, t = 1:n_t],
    p_dr[g,t] <= pmax[g] * w[g,t] * pf[g,t]) #needed for p_da

#@ Constraint: START_S
#@ Description: Generator (slow) start up n
#@ Labels: slowgen_tstart
@constraint(m,START_S[g = GSL, t=1:(n_t-1)],
    z[g,t+1] == w[g,t+1] - w[g,t])

#@ Constraint: START_F
#@ Description: Generator (fast) start up
#@ Labels: fastgen_tstart_scen
@constraint(m,START_F[g=GF ,t=1:(n_t-1), o = SCENARIOS],
    v[g,t+1,o] == u[g,t+1,o] - u[g,t,o])

#@ Constraint: INIT_S
#@ Description: First period commitment for slow generators
#@ Labels: slowgen_init
@constraint(m,INIT_S[g=GSL,t=t_firsts],
    z[g,t] == w[g,t])

slowgen_init = Dict()
for i in 1:n_gsl
    for j in 1:1
        slowgen_init[(i,j)] = "gen$(i)_t$(j)"
    end
end

#@ Constraint: INIT_F
#@ Description: First period commitment for fast generators
#@ Labels: fastgen_init
@constraint(m,INIT_F[g=GF,t=t_firsts, o = SCENARIOS],
    v[g,t,o] == u[g,t,o])

fastgen_init = Dict()
for i in 1:n_gsl
    for j in 1:1
        for k in 1:n_omega
            fastgen_init[(i,j,k)] = "gen$(i)_t$(j)_scen$(k)"
        end
    end
end

#@ Constraint: NAN_ST
@constraint(m,NAN_ST[g = GSL, t = TIME, o = SCENARIOS],
    v[g,t,o] == z[g,t])

#@ Constraint: NAN_CM
@constraint(m,NAN_CM[g = GSL, t = TIME, o = SCENARIOS],
    u[g,t,o] == w[g,t])

#DR_RAND
#think carefully about how to index into p and p_dr
#e.g. need to get from [5,7] for p to [1,2] for p_dr
# one workaround is p[g=GDR,t,o] being its own variable.
# would need to change SUPPLY-DEMAND, GENMAX, and potentially startup/commitment too,
# since these rely on the GENERATORS index...

# another workaround could be to list DR resources first
# so that the indices are the same for both.
# this is kinda hack-y but functional and perhaps cleaner
# for the equations. Lets go with this for now
#@ Constraint: DR_RAND
@constraint(m, DR_RAND[g=GDR,t = TIME, o = SCENARIOS],
     p[g,t,o] == p_dr[g,t] * vdr[o])

#@ Constraint: STARTCOST
#@ Description: Total startup costs are a product of startup costs by generator and the number of startups of that generator $$start\_cost_{g,t,o}\geq v_{g,t,o}startup_g\forall g,t,o$$
#@ Labels: gen_t_scen
@constraint(m, STARTCOST[g=GENERATORS, t = TIME, o = SCENARIOS],
    start_cost[g,t,o] >= v[g,t,o] * startup[g])

#@ Function: objectivefn
#@ Description: minimize costs across time and scenarios $$\min\sum_{g,t,o}(start\_cost_{g,t,o}+p_{g,t,o} varcost_g)$$
@expression(m, objectivefn, sum(pro[o] *
    sum(start_cost[g,t,o] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)   for o = SCENARIOS))

# Solve the model
@objective(m, Min, objectivefn)
println("\nSolving:\n")
JuMP.optimize!(m)
println(string("Status: ",termination_status(m)))

#@ Helper object: startupcost_fraction
#@ Description: Fraction of total costs that are startup costs
startupcost_fraction = JuMP.value.(sum(pro[o] *
    sum(start_cost[g,t,o] for g = GENERATORS, t = TIME)
    for o = SCENARIOS)) / JuMP.value.(objectivefn)

#@ Helper object: varcost_fraction
#@ Description: Fraction of total costs that are var costs
varcost_fraction = JuMP.value.(sum(pro[o] *
    sum(p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS)) / JuMP.value.(objectivefn)

# ok, it solves. Now I just need to figure out how to reasonably
# plot the output so I can check it -- refer to R files to get a
# point of comparison for how I manipulated stuff...

#=
# check production
print("schedule of DR")
x = getvalue(p_dr)
display(x)
print("production of DR:")
y = getvalue(p[1,:,:])
display(y)
print("production of slow generators:")
zs = getvalue(p[2:5,:,:])
display(zs)
print("production of fast generators:")
zf = getvalue(p[6:10,:,:])
display(zf)

# check commitment
print("commitment of slow generators:")
display(getvalue(w))
print("startup of slow generators:")
display(getvalue(z))
print("commitment of all generators")
display(getvalue(u))
print("startup of all generators")
display(getvalue(v))

# check costs
print("total cost")
totcost = getvalue(sum(pro[o] *
    sum(start_cost[g,t,o] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
display(totcost)
print("startup cost")
display("text/plain",getvalue(start_cost))
print("fraction of total costs that are startup costs")
totstartupcost = getvalue(sum(pro[o] *
    sum(start_cost[g,t,o] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
display("text/plain",totstartupcost/totcost)
print("fraction of total costs that are var cost")
totvarcost = getvalue(sum(pro[o] *
    sum(p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
display("text/plain",totvarcost/totcost)
=#

max_stars_repo_name: Fuinn/mos-examples

max_issues_repo_name: Fuinn/mos-examples

max_forks_repo_name: Fuinn/mos-examples


text:
function quantile_ci(x, p=0.05, q=0.5)
    n = length(x)
    zm = quantile(Normal(0.0, 1.0), p/2)
    zp = quantile(Normal(0.0, 1.0), 1-p/2)
    j = max(Int(ceil(n*q + zm*sqrt(n*q*(1-q)))), 1)
    k = min(Int(ceil(n*q + zp*sqrt(n*q*(1-q)))), length(x))
    x_sorted = sort(x)
    return x_sorted[j], quantile(x, q), x_sorted[k]
end

function get_policy(
    policy_type,
    env, num_samples, horizon, λ, α, U₀, cov_mat, pol_log,
    ais_its,
    λ_ais,
    ce_elite_threshold, ce_Σ_est,
    cma_σ, cma_elite_threshold,
)
    if  policy_type == :mppi
        pol = MPPI_Policy(
            env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            log=pol_log,
            rng=MersenneTwister(),
            )
    elseif policy_type == :gmppi
        pol = GMPPI_Policy(
            env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            log=pol_log,
            rng=MersenneTwister(),
        )
    elseif policy_type == :cemppi
        pol = CEMPPI_Policy(
            env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            opt_its=ais_its,
            ce_elite_threshold=ce_elite_threshold,
            Σ_est = ce_Σ_est,
            log=pol_log,
            rng=MersenneTwister(),
        )
    elseif policy_type == :cmamppi
        pol = CMAMPPI_Policy(
            env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            opt_its=ais_its,
            σ=cma_σ,
            elite_perc_threshold=cma_elite_threshold,
            log=pol_log,
            rng=MersenneTwister(),
        )
    elseif policy_type == :μΣaismppi
        pol = μΣAISMPPI_Policy(env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            opt_its=ais_its,
            λ_ais=λ_ais,
            log=pol_log,
            rng=MersenneTwister(),
        )
    elseif policy_type == :μaismppi
        pol = μAISMPPI_Policy(env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            opt_its=ais_its,
            λ_ais=λ_ais,
            log=pol_log,
            rng=MersenneTwister(),
        )
    elseif policy_type == :pmcmppi
        pol = PMCMPPI_Policy(env,
            num_samples=num_samples,
            horizon=horizon,
            λ=λ,
            α=α,
            U₀=U₀,
            cov_mat=cov_mat,
            opt_its=ais_its,
            λ_ais=λ_ais,
            log=pol_log,
            rng=MersenneTwister(),
        )
    else
        error("No policy_type of $policy_type")
    end
    return pol
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
export HopfieldNetwork

type HopfieldNetwork{T<:Number}
	state::Array{T,1}
	weight::Array{T,2}
	n::Integer
	typ::Type
	function HopfieldNetwork(n::Integer)
		state = zeros(T,n)
		Weight = zeros(T,(n,n))
		new(state,Weight,n,T)
	end
end

function Base.show(io::IO, net::HopfieldNetwork)
	@printf io "HopfieldNetwork(type: %s, dim: %d)" string(net.typ) net.n
end


text:
## @primitive

function primitive_sig(sig)
    fn, args = sig.args[1], sig.args[2:end]
    targs = Array{Any}(undef, length(args))
    vnames = Array{Symbol}(undef, length(args))
    for i=1:length(args)
        if args[i] isa Expr
            name, typ_ex = args[i].args
            typ = Core.eval(Base, typ_ex)
            ttyp = typ <: AbstractArray ? TArray : TReal
            targs[i] = :($name :: $ttyp)
            vnames[i] = name
        elseif args[i] isa Symbol
            name = args[i]
            targs[i] = :($name :: $TAny)
            vnames[i] = name
        else
            error("Unexpected signature in @primitive: $(sig)")
        end
    end
    tsig = :($fn($(targs...)))
    return tsig, fn, vnames
end

function _primitive(sig)
    @assert sig.head == :call
    p_sig, fn, vnames = primitive_sig(sig)
    ex = :($p_sig = $record!($(vnames[1]).tape, $Call, $fn, ($(vnames...),)))
    return ex
end

macro primitive(sig)
    return esc(_primitive(sig))
end

## @grad

function grad_sig(sig, idx)
    fn, args = sig.args[1], sig.args[2:end]
    targs = Array{Any}(undef, length(args))
    type_params = Symbol[]
    type_param_idx = 1
    vnames = Array{Symbol}(undef, length(args))
    for i=1:length(args)
        if args[i] isa Expr
            name, typ_ex = args[i].args
            typ = Core.eval(Base, typ_ex)
            if typ  <: AbstractArray
                T, N = Symbol("T$type_param_idx"), Symbol("N$type_param_idx")
                type_param_idx += 1
                push!(type_params, T, N)
                ttyp = Expr(:curly, TArray, T, N)
            else
                ttyp = TReal
            end
            vnames[i] = name
            targs[i] = ttyp
        elseif args[i] isa Symbol
            name = args[i]
            vnames[i] = name
            targs[i] = TAny
        else
            error("Unexpected signature in @primitive: $(sig)")
        end
    end
    grad_fn = Expr(:., Symbol(@__MODULE__), QuoteNode(:grad!))
    if isempty(type_params)
        tsig = :($grad_fn(dy::$TAny, ::Val{$idx}, op::$Call{typeof($fn), Tuple{$(targs...)}}))
    else
        tsig = :($grad_fn(dy::$TAny, ::Val{$idx}, op::$Call{typeof($fn), Tuple{$(targs...)}})
                 where {$(type_params...)})
    end
    return tsig, vnames
end

function grad_body(vnames, body)
    tbody = (body isa Expr && body.head == :block) ? body : :(begin $body end)
    init_vars_ex = :(($(vnames...),) = op.args)
    pushfirst!(tbody.args, init_vars_ex)
    return tbody
end

function _grad(sig, idx, body)
    tsig, vnames = grad_sig(sig, idx)
    tbody = grad_body(vnames, body)
    ex = :($tsig = $tbody)
    return ex
end

macro grad(sig, idx, body)
    esc(_grad(sig, idx, body))
end


text:
module SacredDataStore

using ValueHistories
using JSON
using MacroTools

using DataFrames
using DataFramesMeta

using JSON3
using Measurements
using Dates

include("dataframes_utils.jl")
include("valuehistories_utils.jl")
include("metrics.jl")

include("load.jl")
include("nk_json_log.jl")

end


text:
using SpatialFields
using BenchmarkTools
import StaticArrays: SVector

suite = BenchmarkGroup()

let
    suite["hrbf"] = BenchmarkGroup(["hrbf"])
    points = [rand(SVector{3, Float64}) for i in 1:100]
    normals = [rand(SVector{3, Float64}) for i in 1:100];
    suite["hrbf"]["construct"] = @benchmarkable field = HermiteRadialField($points, $normals)

    field = HermiteRadialField(points, normals)
    x = zeros(SVector{3, Float64})
    suite["hrbf"]["evaluate"] = @benchmarkable ($field)($x)
end

let
    suite["interpolating"] = BenchmarkGroup(["interpolating"])
    edge_points = [rand(SVector{3, Float64}) for i in 1:50]
    interior_points = [rand(SVector{3, Float64}) for i in 1:50]
    points = vcat(edge_points, interior_points)
    values = vcat([0.0 for i in edge_points], [-1.0 for i in interior_points])

    suite["interpolating"]["construct"] = @benchmarkable surface = InterpolatingSurface($points, $values)
    surface = InterpolatingSurface(points, values)
    x = zeros(SVector{3, Float64})
    suite["interpolating"]["evaluate_svector"] = @benchmarkable ($surface)($x)

    x = zeros(3)
    suite["interpolating"]["evaluate_vector"] = @benchmarkable ($surface)($x)
end

tune!(suite)
results = run(suite, verbose=false)
showall(results)
println()


text:
#!/usr/bin/env julia

using Printf
using CUDAnative
using CuArrays

N=10
a_d=CuArrays.fill(2.0::Float64, N)
b_d=CuArrays.fill(3.0::Float64, N)
a  =zeros(Float64,N)
b  =zeros(Float64,N)
r  =zeros(Float64,N)
copyto!(a,a_d)
copyto!(b,b_d)

c = a .+ b
for j in c
  @printf("cpu -> %i\n",j)
end

function add_works!(a,b)
 for i in 1:length(a)
   a[i] += b[i]
 end
 nothing
end

# this works
@cuda  threads=1 add_works!(a_d,b_d)

function add_fails!(a,b)
  a .+= b
  nothing
end

# this fails
@cuda  threads=1 add_fails!(a_d,b_d)

copyto!(r,a_d)

for i in r
 @printf("gpu -> %i\n",i)
end

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
module ReactiveBasics

using DocStringExtensions, DataStructures

export Signal, value, foldp, subscribe!, unsubscribe!, flatmap, flatten, bind!,
       droprepeats, skip, previous, sampleon, preserve, filterwhen, zipmap

# This API mainly follows that of Reactive.jl.

# The algorithms for were derived from the Swift Interstellar package
# https://github.com/JensRavens/Interstellar/blob/master/Sources/Signal.swift
# Copyright (c) 2015 Jens Ravens (http://jensravens.de)
# Offered with the MIT license

"""
A `Signal` is value that will contain a value in the future.
The value of the Signal can change at any time.

Use `map` to derive new Signals, `subscribe!` to subscribe to updates of a Signal,
and `push!` to update the current value of a Signal.
`value` returns the current value of a Signal.
The type of the Signal can optionally be set by the first parameter. Otherwise
it will default to the type of the initial value.

```julia
text = Signal("")

text2 = map(s -> "Bye \$s", text)

subscribe!(text) do s
    println("Hello \$s")
end

push!(text, "world")

value(text)
value(text2)

float_number = Signal(Float64, 1) # Optionally set the type of the Signal
```
"""
mutable struct Signal{T}
   value::T
   callbacks::Vector{Function}   # usually Functions, but could be other callable types
end

Signal(val::T) where {T} = Signal{T}(val, Function[])
Signal(::Type{T}, x) where{T} = Signal{T}(x, Function[])

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

"""
$(SIGNATURES)

The current value of the Signal `u`.
"""
value(u::Signal) = u.value

"""
$(SIGNATURES)

Transform the Signal into another Signal using a function. The initial value of
the output Signal can optionally be set via `init`. Otherwise it defaults to
`f(u.value)`, where f is the passed function and u is the passed Signal. The type
of the output Signal can optionally be set via `typ`. Otherwise it defaults to
the type of the initial value.
"""
function Base.map(f, u::Signal; init = f(u.value), typ = typeof(init))
    signal = Signal(typ, init)
    subscribe!(x -> push!(signal, f(x)), u)
    signal
end
function Base.map(f, u::Signal, v::Signal; init = f(u.value, v.value), typ = typeof(init))
    signal = Signal(typ, init)
    subscribe!(x -> push!(signal, f(x, v.value)), u)
    subscribe!(x -> push!(signal, f(u.value, x)), v)
    signal
end
function Base.map(f, u::Signal, v::Signal, w::Signal; init = f(u.value, v.value, w.value), typ = typeof(init))
    signal = Signal(typ, init)
    subscribe!(x -> push!(signal, f(x, v.value, w.value)), u)
    subscribe!(x -> push!(signal, f(u.value, x, w.value)), v)
    subscribe!(x -> push!(signal, f(u.value, v.value, x)), w)
    signal
end
function Base.map(f, u::Signal, v::Signal, w::Signal, xs::Signal...; init = f((u.value for u in (u,v,w,xs...))...), typ = typeof(init))
    us = (u,v,w,xs...)
    signal = Signal(typ, init)
    for (i,u) in enumerate(us)
        subscribe!(u) do x
            vals = f(((i == j ? x : us[j].value for j in 1:length(us))...)...)
            push!(signal, vals)
        end
    end
    signal
end

"""
$(SIGNATURES)

Transform the Signal into another Signal using a function. It's like `map`,
but it's meant for functions that return `Signal`s. The initial value of
the output Signal can optionally be set via `init`. Otherwise it defaults to
`f(input.value).value`, where f is the passed function and input is the passed
Signal of Signals. The type of the output Signal can optionally be set via `typ`.
Otherwise it defaults to the type of the initial value.
"""
function flatmap(f, input::Signal; init = f(input.value).value, typ = typeof(init))
    signal = Signal(typ, init)
    subscribe!(input) do u
        innersig = f(u)
        push!(signal, innersig.value)
        subscribe!(innersig) do v
            push!(signal, v)
        end
    end
    signal
end

"""
$(SIGNATURES)

Zips given signals first and then applies the `map` function onto the
zipped value. This omits the double calculation when using `map`. The
initial value of the output Signal can optionally be set via `init`. Otherwise
it defaults to `f(zip(u, us...).value...)`, where f is the passed function and
(u, us...) are the passed Signals. The type of the output Signal can optionally
be set via `typ`. Otherwise it defaults to the type of the initial value.

    as = Signal(1)
    bs = map(a -> a * 0.1, as)
    cs = zipmap((a,b) -> a + b, as, bs) # This calculation is done once for
                                        # every change in `as`
"""
function zipmap(f, u::Signal, us::Signal...; init = f(zip(u, us...).value...), typ = typeof(init), max_buffer_size = 0)
    zipped_signal = zip(u, us..., max_buffer_size = max_buffer_size)
    signal = Signal(typ, init)
    subscribe!(x -> push!(signal, f(x...)), zipped_signal)
    signal
end

"""
$(SIGNATURES)

Update the value of a Signal and propagate the change.
"""
function Base.push!(u::Signal, val)
    u.value = val
    foreach(f -> f(val), u.callbacks)
end

"""
$(SIGNATURES)

Subscribe to the changes of this Signal. Every time the Signal is updated, the function `f` runs.
"""
function subscribe!(f, u::Signal)
    push!(u.callbacks, f)
    u
end

"""
$(SIGNATURES)

Unsubscribe to the changes of this Signal.
"""
function unsubscribe!(f, u::Signal)
    u.callbacks = filter(a -> a != f, u.callbacks)
end

"""
$(SIGNATURES)

Zip (combine) Signals into the current Signal. The value of the Signal is a
Tuple of the values of the contained Signals.

    signal = zip(Signal("Hello"), Signal("World"))
    value(signal)    # ("Hello", "World")
"""
function Base.zip(u::Signal, us::Signal...; max_buffer_size = 0)
    signals = (u,us...)
    signal = Signal(Tuple{map(eltype, signals)...}, map(value, signals))
    pasts = map(u -> max_buffer_size != 0 ? CircularDeque{eltype(u)}(max_buffer_size) : Deque{eltype(u)}(), signals)
    for i in 1:length(signals)
        subscribe!(signals[i]) do u
            push!(pasts[i], u)
            if all(map(!isempty, pasts))
                push!(signal, map(popfirst!, pasts))
            end
        end
    end
    signal
end

"""
$(SIGNATURES)

Merge Signals into the current Signal. The value of the Signal is that from
the most recent update.
"""
function Base.merge(u::Signal, us::Signal...)
    signal = Signal(typejoin(map(eltype, (u, us...))...), u.value)
    for v in (u, us...)
        subscribe!(x -> push!(signal, x), v)
    end
    signal
end

"""
$(SIGNATURES)

Fold/map over past values. The first argument to the function `f`
is an accumulated value that the function can operate over, and the
second is the current value coming in. `v0` is the initial value of
the accumulated value.

    a = Signal(2)
    # accumulate sums coming in from a, starting at zero
    b = foldp(+, 0, a) # b == 2
    push!(a, 2)        # b == 4
    push!(a, 3)        # b == 7
"""
function foldp(f, v0, us::Signal...)
    v0r = Base.RefValue(v0)
    map((x...) -> v0r[] = f(v0r[], x...), us...)
end

"""
Return a Signal that updates based on the Signal `u` if `f(value(u))` evaluates to `true`.
"""
function Base.filter(f, default, u::Signal)
    signal = Signal(f(u.value) ? u.value : default)
    subscribe!(result -> f(result) && push!(signal, result), u)
    signal
end

"""
$(SIGNATURES)

Keep updates to `u` only when `predicate` is true.
If `predicate` is false initially, the specified `default` value is used.
"""
function filterwhen(predicate::Signal{Bool}, default, u::Signal)
    signal = Signal(predicate.value ? u.value : default)
    subscribe!(result -> predicate.value && push!(signal, result), u)
    subscribe!(v -> v && push!(signal, u.value), predicate)
    signal
end

"""
$(SIGNATURES)

An asynchronous version of `map` that returns a Signal that is updated after `f`
operates asynchronously. The initial value of the returned Signal
(the `init` arg) must be supplied. The type of the output Signal can optionally
be set via `typ`. Otherwise it defaults to the type of the initial value.
"""
function Base.asyncmap(f, init, input::Signal, inputs::Signal...; typ = typeof(init))
    result = Signal(typ, init)
    map(input, inputs...) do args...
        @async begin
            push!(result, f(args...))
        end
    end
    result
end

"""
$(SIGNATURES)

Flatten a Signal of Signals into a Signal which holds the
value of the current Signal.
"""
function flatten(input::Signal)
    sigref = Base.RefValue(input.value)
    signal = Signal(input.value.value)
    updater = u -> push!(signal, u)
    subscribe!(updater, input.value)
    subscribe!(input) do u
        push!(signal, u.value)
        unsubscribe!(updater, sigref[])
        subscribe!(updater, u)
        sigref[] = u
    end
    push!(input, input.value)
    signal
end

"""
$(SIGNATURES)

For every update to `b` also update `a` with the same value and, if `twoway` is true,
ice-versa.
If `initial` is set to true, `a` is updated immediately and, if `twoway` is true, `b`
is updated immediately as well.
"""
function bind!(a::Signal, b::Signal, twoway = true; initial = true)
    if initial
        push!(a, b.value)
        if twoway
            push!(b, a.value)
        end
    end
    if twoway
        subscribe!(u -> u != value(b) && push!(b, u), a)
    end
    subscribe!(u -> u != value(a) && push!(a, u), b)
end

"""
$(SIGNATURES)

Drop updates to `input` whenever the new value is the same
as the previous value of the Signal.
"""
function droprepeats(input::Signal)
    result = Signal(value(input))
    subscribe!(u -> u != value(result) && push!(result, u), input)
    result
end

"""
$(SIGNATURES)

Continuously skip a predefined number of updates to `input`
"""
function Base.skip(num::Int, input::Signal)
    result = Signal(value(input))
    counter = 1
    subscribe!(input) do u
        mod(counter, num + 1) == 0 && push!(result, u)
        counter = mod(counter, num + 1) + 1
    end
    result
end

"""
$(SIGNATURES)

Create a Signal which holds the previous value of `input`.
You can optionally specify a different initial value.
"""
function previous(input::Signal, default=value(input))
    past = Base.RefValue(default)
    map(input) do u
        res = past[]
        past[] = u
        res
    end
end

"""
$(SIGNATURES)

Sample the value of `b` whenever `a` updates.
"""
function sampleon(a::Signal, b::Signal)
    result = Signal(value(b))
    subscribe!(u -> push!(result, value(b)), a)
    result
end

"""
$(SIGNATURES)

For compatibility with Reactive.
It just returns the original Signal because this isn't needed with direct `push!` updates.
"""
function preserve(x::Signal)
    x
end

end # module


text:
function calc_eta_CD(cur_reactor::AbstractReactor; strong_fail::Bool=false)
  cur_rho_m = _rho_m(cur_reactor)

  isnan(cur_rho_m) && return NaN

  cur_rho_j_list = _rho_j(cur_reactor, rho_sym)

  isempty(cur_rho_j_list) && return NaN

  for (cur_index, cur_rho_j) in enumerate(cur_rho_j_list)
    ( cur_rho_m < cur_rho_j ) && continue

    tmp_rho_j_list = _rho_j(cur_reactor, cur_rho_m)

    if isempty(tmp_rho_j_list)
      cur_rho_j_list[cur_index] = NaN
      continue
    end

    if strong_fail
      @assert length(tmp_rho_j_list) == 1
    else
      if length(tmp_rho_j_list) != 1
        custom_log("Too many rho_j in tmp list")
        custom_log(tmp_rho_j_list)
        custom_log(cur_reactor)
        return NaN
      end
    end

    cur_rho_j_list[cur_index] = first(tmp_rho_j_list)
  end

  filter!(!isnan, cur_rho_j_list)
  isempty(cur_rho_j_list) && return NaN

  cur_rho_j_list = unique(cur_rho_j_list)
  filter_approx!(cur_rho_j_list)

  if strong_fail
    @assert length(cur_rho_j_list) == 1
  else
    if length(cur_rho_j_list) != 1
      custom_log("Too many rho_j in cur list")
      custom_log(cur_rho_j_list)
      custom_log(cur_reactor)
      return NaN
    end
  end

  cur_rho_j = first(cur_rho_j_list)

  cur_eta_twiddle = _wave_eta_twiddle(cur_reactor, cur_rho_j)

  cur_eta_CD = 0.06108 * cur_reactor.eta_LH

  cur_eta_CD /= log( plasma_parameter(cur_reactor) )

  cur_eta_CD *= ( 1 + cur_reactor.nu_T )

  cur_eta_CD *= cur_reactor.T_bar

  cur_eta_CD *= ( 1 - cur_rho_j ^ 2 ) ^ cur_reactor.nu_T

  cur_eta_CD *= cur_eta_twiddle

  cur_eta_CD
end

@symbol_func function _rho_j(cur_reactor::AbstractReactor, nested_rho::AbstractSymbol)
  cur_rho = rho_sym

  cur_nu_T = cur_reactor.nu_T

  cur_equation = _n_para_2(cur_reactor, nested_rho)

  cur_equation *= ( 1 - cur_rho ^ 2 ) ^ cur_nu_T

  cur_equation *= ( 1 + cur_nu_T )

  cur_equation -= ( 28.39 / cur_reactor.T_bar )

  solved_rhos = []

  for cur_attempt in 1:4
    solved_rhos = find_roots(
      cur_equation, min_rho, max_rho,
      no_pts=no_pts_rho, abstol = 1e-4
    )

    isempty(solved_rhos) || break
  end

  filter_approx!(solved_rhos)

  solved_rhos
end

function _rho_m(cur_reactor::AbstractReactor)
  cur_rho = rho_sym

  cur_equation = _n_para_2(cur_reactor, cur_rho)

  is_nan(cur_equation) && return NaN

  cur_equation *= -1

  cur_func = function (tmp_rho)
    cur_value = float(cur_equation(tmp_rho))
    isa(cur_value, Complex) && return Inf
    cur_value
  end

  solved_rho = NaN

  try
    solved_rho = optimize(
      cur_func, min_rho, max_rho
    ).minimizer
  catch cur_error
    isa(cur_error, DomainError) || rethrow(cur_error)
  end

  solved_rho
end

@symbol_func function _n_para_2(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_X = _wave_X(cur_reactor, cur_rho)

  cur_omega_2 = _omega_hat_2(cur_reactor)

  cur_n = sqrt( 1.0 + cur_X - ( cur_X / cur_omega_2 ) )

  cur_n += sqrt( cur_X )

  cur_n ^= 2
end

@symbol_func function _omega_hat_2(cur_reactor::AbstractReactor)
  cur_X = _wave_X(cur_reactor, 0.0)

  cur_omega_2 = 4.0 * wave_gamma ^ 2

  cur_omega_2 *= cur_X * ( 1 + cur_X )

  cur_omega_2 += 1

  cur_omega_2 ^= 1/2

  cur_omega_2 += 1

  cur_omega_2 /= 2

  cur_omega_2 *= cur_X / ( 1 + cur_X )
end

@symbol_func function _wave_X(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_X = omega_pe(cur_reactor, cur_rho) ^ 2

  cur_X /= omega_ce(cur_reactor, cur_rho) ^ 2

  cur_X
end

@symbol_func function _wave_w(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_w = cur_reactor.speed_of_light

  cur_w /= v_te(cur_reactor, cur_rho)

  cur_w /= sqrt(_n_para_2(cur_reactor, cur_rho))

  cur_w
end

@symbol_func function _wave_B_M(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_B = B_rho(cur_reactor, cur_rho, 180.0)
end

@symbol_func function _wave_x_t_2(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_B = B_rho(cur_reactor, cur_rho, cur_reactor.wave_theta)

  cur_x_2 = _wave_w(cur_reactor, cur_rho) ^ 2

  cur_x_2 *= cur_B

  cur_x_2 /= _wave_B_M(cur_reactor, cur_rho) - cur_B

  cur_x_2
end

@symbol_func function _wave_K(cur_reactor::AbstractReactor)
  cur_K = 2.12

  cur_K /= cur_reactor.Z_eff

  cur_K
end

@symbol_func function _wave_D(cur_reactor::AbstractReactor)
  cur_D = 3.83

  cur_D /= cur_reactor.Z_eff ^ sqrt(1/2)

  cur_D
end

@symbol_func function _wave_eta_0(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_w = _wave_w(cur_reactor, cur_rho)

  cur_eta = ( _wave_K(cur_reactor) / cur_w )

  cur_eta += _wave_D(cur_reactor)

  cur_eta += ( 8 * cur_w ^ 2 ) / ( 5 + cur_reactor.Z_eff )

  cur_eta
end

@symbol_func function _wave_C(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_x_2 = _wave_x_t_2(cur_reactor, cur_rho)

  cur_C = 1.0

  cur_C -= exp( -1.0 * ( wave_c * cur_x_2 ) ^ wave_m )

  cur_C
end

@symbol_func function _wave_R(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_w = _wave_w(cur_reactor, cur_rho)

  cur_numerator = sqrt( wave_x_r ^ 2 + cur_w ^ 2 )

  cur_denominator = cur_w

  cur_denominator /= ( cur_reactor.epsilon * cur_rho ) ^ wave_n

  cur_denominator += wave_x_r

  cur_R = 1.0

  cur_R -= ( cur_numerator / cur_denominator )

  cur_R
end

@symbol_func function _wave_eta_twiddle(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_eta = wave_M

  cur_eta *= _wave_C(cur_reactor, cur_rho)

  cur_eta *= _wave_R(cur_reactor, cur_rho)

  cur_eta *= _wave_eta_0(cur_reactor, cur_rho)

  cur_eta
end


text:
using BioMetaheuristics
using BioMetaheuristics.TestFunctions
using Test
using StaticArrays
using Random

@testset "BioMetaheuristics.jl" begin
    include("functests.jl")
    include("populationtests.jl")
    include("solvers-tests/pso.jl")
    include("solvers-tests/results.jl")
    include("solvers-tests/simulated_annealing.jl")
    include("apitests.jl")
end


text:
# Create model
function initialize_model()
    abm = AgentBasedModel(getAgentType(), getSpace())
    agent = SimplePolicyAgent(0,0,(0,0))
    add_agent!(agent, abm)

    abm
end


text:
module Settings

using PyCall: PyObject

using PyQHA: qha_settings, @pyinterface

export from_yaml

mutable struct Config
    o::PyObject
end
Config(args::AbstractDict...) = Config(qha_settings.Settings(args...))

Base.getindex(config::Config, key) = PyObject(config).get(key)
Base.setindex!(config::Config, value, key) = PyObject(config).__setitem__(key, value)

@pyinterface Config

from_yaml(file) = Config(qha_settings.from_yaml(file))

end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function _co_nutate(jd::T, ra::T, dec::T) where {T<:AbstractFloat}
    d_psi, d_eps = nutate(jd)
    eps = mean_obliquity(jd) + sec2rad(d_eps)
    se, ce = sincos(eps)
    sin_ra, cos_ra = sincos(deg2rad(ra))
    sin_dec, cos_dec = sincos(deg2rad(dec))
    x = cos_ra * cos_dec
    y = sin_ra * cos_dec
    z = sin_dec
    x2 = x - sec2rad(y * ce + z * se) * d_psi
    y2 = y + sec2rad(x * ce * d_psi - z * d_eps)
    z2 = z + sec2rad(x * se * d_psi + y * d_eps)
    xyproj = hypot(x2, y2)
    r = hypot(xyproj, z2)
    ra2 = atan(y2, x2)
    dec2 = asin(z2/r)
    ra2 = rad2deg(ra2)

    if ra2 < 0
        ra2 += 360
    end
    d_ra = ra2 - ra
    d_dec = rad2deg(dec2) - dec
    return d_ra, d_dec, eps, d_psi, d_eps
end

"""
    co_nutate(jd, ra, dec) -> d_ra, d_dec, eps, d_psi, d_eps

### Purpose ###

Calculate changes in RA and Dec due to nutation of the
Earth's rotation

### Explanation ###

Calculates necessary changes to ra and dec due to the nutation of the
Earth's rotation axis, as described in Meeus, Chap 23. Uses formulae
from Astronomical Almanac, 1984, and does the calculations in equatorial
rectangular coordinates to avoid singularities at the celestial poles.

### Arguments ###

* `jd`: julian date, scalar or vector
* `ra`: right ascension in degrees, scalar or vector
* `dec`: declination in degrees, scalar or vector

### Output ###

The 5-tuple `(d_ra, d_dec, eps, d_psi, d_eps)`:

* `d_ra`: correction to right ascension due to nutation, in degrees
* `d_dec`: correction to declination due to nutation, in degrees
* `eps`: the true obliquity of the ecliptic
* `d_psi`: nutation in the longitude of the ecliptic
* `d_eps`: nutation in the obliquity of the ecliptic

### Example ###

Example 23a in Meeus: On 2028 Nov 13.19 TD the mean position of Theta
Persei is 2h 46m 11.331s 49d 20' 54.54''. Determine the shift in
position due to the Earth's nutation.

```jldoctest
julia> using AstroLib

julia> jd = jdcnv(2028,11,13,4,56)
2.4620887055555554e6

julia> co_nutate(jd,ten(2,46,11.331) * 15,ten(49,20,54.54))
(0.004400660977140092, 0.00172668646508356, 0.40904016038217555, 14.859389427896476, 2.7038090372350574)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.

The output of `d_ra` and `d_dec` in IDL AstroLib is in arcseconds,
however it is in degrees here.

This function calls [`mean_obliquity`](@ref) and [`nutate`](@ref).
"""
co_nutate(jd::Real, ra::Real, dec::Real) =
    _co_nutate(promote(float(jd), float(ra), float(dec))...)

function co_nutate(jd::AbstractVector{P}, ra::AbstractVector{<:Real},
                   dec::AbstractVector{<:Real}) where {P<:Real}
    @assert length(jd) == length(ra) == length(dec) "jd, ra and dec vectors
                                                     should be of the same length"
    typejd = float(P)
    ra_out  = similar(jd,  typejd)
    dec_out = similar(dec, typejd)
    eps_out = similar(dec, typejd)
    d_psi_out = similar(dec, typejd)
    d_eps_out = similar(dec, typejd)
    for i in eachindex(jd)
        ra_out[i], dec_out[i],eps_out[i], d_psi_out[i], d_eps_out[i]  =
        co_nutate(jd[i], ra[i], dec[i])
    end
    return ra_out, dec_out, eps_out, d_psi_out, d_eps_out
end


text:
mutable struct Fighter <: AbstractAgent
    id::Int
    pos::Dims{3}
    has_prisoner::Bool
    capture_time::Int
    shape::Symbol
end

"""
```julia
battle(; fighters = 50)
```
Same as in [Battle Royale](@ref).
"""
function battle(; fighters = 50)
    model = ABM(
        Fighter,
        GridSpace((100, 100, 10); periodic = false);
        scheduler = Schedulers.randomly,
    )

    n = 0
    while n != fighters
        pos = (rand(model.rng, 1:100, 2)..., 1) # Start at level 1
        if isempty(pos, model)
            add_agent!(pos, model, false, 0, :diamond)
            n += 1
        end
    end

    return model, battle_agent_step!, dummystep
end

space(agent) = agent.pos[1:2]
level(agent) = agent.pos[3]

function closest_target(agent::Fighter, ids::Vector{Int}, model::ABM)
    if length(ids) == 1
        closest = ids[1]
    else
        close_id = argmin(map(id -> edistance(space(agent), space(model[id]), model), ids))
        closest = ids[close_id]
    end
    return model[closest]
end

function battle!(one::Fighter, two::Fighter, model)
    if level(one) == level(two)
        # Odds are equivalent
        one_winner = rand(model.rng) < 0.5
    elseif level(one) > level(two)
        # Odds are in favor of one
        one_winner = 2 * rand(model.rng) > rand(model.rng)
    else
        # Odds are in favor of two
        one_winner = rand(model.rng) > 2 * rand(model.rng)
    end

    one_winner ? (up = one; down = two) : (up = two; down = one)

    new_lvl_up = min(level(up) + 1, 10)
    new_pos_up = clamp.(rand(model.rng, -1:1, 2) .+ space(up), [1, 1], size(model.space)[1:2])
    move_agent!(up, (new_pos_up..., new_lvl_up), model)
    new_lvl_down = level(down) - 1
    if new_lvl_down == 0
        kill_agent!(down, model)
    else
        move_agent!(down, (space(down)..., new_lvl_down), model)
    end
end

function captor_behavior!(agent, model)
    close_ids = collect(nearby_ids(agent, model, (0, 0, 10)))
    if length(close_ids) == 1
        # Taunt prisoner or kill it
        prisoner = model[close_ids[1]]
        if prisoner.capture_time > 10
            agent.shape = :square
            gain = ceil(Int, level(prisoner) / 2)
            new_lvl = min(level(agent) + gain, 10)
            kill_agent!(prisoner, model)
            agent.has_prisoner = false
            move_agent!(agent, (space(agent)..., new_lvl), model)
        end
    else
        # Someone is here to kill the captor. Could be more than one opponent
        prisoner = [model[id] for id in close_ids if model[id].capture_time > 0][1]
        exploiter = rand(model.rng, [
            model[id]
            for
            id in close_ids if
            model[id].capture_time == 0 && model[id].has_prisoner == false
        ])
        exploiter.shape = :square
        gain = ceil(Int, level(agent) / 2)
        new_lvl = min(level(agent) + rand(model.rng, 1:gain), 10)
        kill_agent!(agent, model)
        move_agent!(exploiter, (space(exploiter)..., new_lvl), model)
        # Prisoner runs away in the commotion
        prisoner.shape = :utriangle
        prisoner.capture_time = 0
        walk!(prisoner, (rand(model.rng, -1:1, 2)..., 0), model)
    end
end

function endgame!(agent, model)
    origin = space(agent)
    end_ids = collect(Iterators.filter(
        id -> model[id].shape == :circle && id != agent.id,
        allids(model),
    ))
    agent.shape = :circle
    if !isempty(end_ids)
        opponent = closest_target(agent, end_ids, model)
        target = space(opponent)
        if origin == target
            # Battle
            agent.shape = :square
            opponent.shape = :square
            showdown!(agent, opponent, model)
        else
            walk!(agent, (sign.(target .- origin)..., 0), model)
        end
    end
end

function showdown!(one::Fighter, two::Fighter, model)
    if level(one) == level(two)
        # Odds are equivalent
        one_winner = rand(model.rng) < 0.5
    elseif level(one) > level(two)
        # Odds are in favor of one
        one_winner = level(one) - level(two) * rand(model.rng) > rand(model.rng)
    else
        # Odds are in favor of two
        one_winner = rand(model.rng) > level(two) - level(one) * rand(model.rng)
    end

    one_winner ? kill_agent!(two, model) : kill_agent!(one, model)
end

function battle_agent_step!(agent, model)
    if agent.capture_time > 0
        # Captured agents are powerless, but we need to keep track of how long
        # they have been in this state
        agent.capture_time += 1
    elseif agent.has_prisoner
        captor_behavior!(agent, model)
    else
        origin = space(agent)
        # Find agents that have captives, they are not focused
        occupied_ids = collect(Iterators.filter(
            id -> model[id].has_prisoner,
            nearby_ids(agent, model, (7, 7, 10)),
        ))
        if !isempty(occupied_ids)
            # Sneak up behind them
            target = space(closest_target(agent, occupied_ids, model))
            agent.shape = :pentagon
            walk!(agent, (sign.(target .- origin)..., 0), model)
        else
            # Opponents that are greatly higher in rank that the current agent
            strong_ids = collect(nearby_ids(agent, model, [(1, -5:5), (2, -5:5), (3, 2:4)]))
            if !isempty(strong_ids)
                # Run away from nearest
                target = space(closest_target(agent, strong_ids, model))
                agent.shape = :utriangle
                walk!(agent, (sign.(origin .- target)..., 0), model)
            else
                # There are no distractions. Search for the closest worthy opponent
                worthy_ids = collect(nearby_ids(agent, model, [(3, -1:1)]))
                if !isempty(worthy_ids)
                    opponent = closest_target(agent, worthy_ids, model)
                    target = space(opponent)
                    if origin == target
                        # Battle
                        agent.shape = :square
                        opponent.shape = :square
                        battle!(agent, opponent, model)
                    else
                        # Move towards worthy opponent
                        agent.shape = :diamond
                        walk!(agent, (sign.(target .- origin)..., 0), model)
                    end
                else
                    # Find any weak targets in the vicinity
                    weak_ids = collect(nearby_ids(
                        agent,
                        model,
                        [(1, -10:10), (2, -10:10), (3, -4:-2)],
                    ))
                    if !isempty(weak_ids)
                        prisoner = closest_target(agent, weak_ids, model)
                        target = space(prisoner)
                        if origin == target
                            #Capture and taunt target
                            agent.has_prisoner = true
                            agent.shape = :vline
                            prisoner.capture_time += 1
                            prisoner.shape = :hline
                        else
                            # Chase down nearest (can move 2 steps at a time!)
                            agent.shape = :star4
                            walk!(agent, (2 .* sign.(target .- origin)..., 0), model)
                        end
                    else
                        # Abandon honour. This is the end
                        endgame!(agent, model)
                    end
                end
            end
        end
    end

    return nothing
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export AF_BACKEND_CPU, AF_BACKEND_CUDA, AF_BACKEND_DEFAULT, AF_BACKEND_OPENCL, AF_BINARY_ADD, AF_BINARY_MAX
export AF_BINARY_MIN, AF_BINARY_MUL, AF_CANNY_THRESHOLD_AUTO_OTSU, AF_CANNY_THRESHOLD_MANUAL, AF_COLORMAP_BLUE
export AF_COLORMAP_COLORS, AF_COLORMAP_DEFAULT, AF_COLORMAP_HEAT, AF_COLORMAP_INFERNO, AF_COLORMAP_MAGMA
export AF_COLORMAP_MOOD, AF_COLORMAP_PLASMA, AF_COLORMAP_RED, AF_COLORMAP_SPECTRUM, AF_COLORMAP_VIRIDIS
export AF_CONNECTIVITY_4, AF_CONNECTIVITY_8, AF_CONV_AUTO, AF_CONV_DEFAULT, AF_CONV_EXPAND, AF_CONV_FREQ
export AF_CONV_SPATIAL, AF_DIFFUSION_DEFAULT, AF_DIFFUSION_GRAD, AF_DIFFUSION_MCDE, AF_ERR_ARG, AF_ERR_ARR_BKND_MISMATCH
export AF_ERR_BATCH, AF_ERR_DEVICE, AF_ERR_DIFF_TYPE, AF_ERR_DRIVER, AF_ERR_INTERNAL, AF_ERR_INVALID_ARRAY
export AF_ERR_LOAD_LIB, AF_ERR_LOAD_SYM, AF_ERR_NONFREE, AF_ERR_NOT_CONFIGURED, AF_ERR_NOT_SUPPORTED, AF_ERR_NO_DBL
export AF_ERR_NO_GFX, AF_ERR_NO_MEM, AF_ERR_RUNTIME, AF_ERR_SIZE, AF_ERR_TYPE, AF_ERR_UNKNOWN, AF_FIF_BMP
export AF_FIF_EXR, AF_FIF_HDR, AF_FIF_ICO, AF_FIF_JNG, AF_FIF_JP2, AF_FIF_JPEG, AF_FIF_PNG, AF_FIF_PPM
export AF_FIF_PPMRAW, AF_FIF_PSD, AF_FIF_RAW, AF_FIF_TIFF, AF_FLUX_DEFAULT, AF_FLUX_EXPONENTIAL, AF_FLUX_QUADRATIC
export AF_GRAY, AF_HOMOGRAPHY_LMEDS, AF_HOMOGRAPHY_RANSAC, AF_HSV, AF_ID, AF_INTERP_BICUBIC, AF_INTERP_BICUBIC_SPLINE
export AF_INTERP_BILINEAR, AF_INTERP_BILINEAR_COSINE, AF_INTERP_CUBIC, AF_INTERP_CUBIC_SPLINE, AF_INTERP_LINEAR
export AF_INTERP_LINEAR_COSINE, AF_INTERP_LOWER, AF_INTERP_NEAREST, AF_LSAD, AF_LSSD, AF_MARKER_CIRCLE
export AF_MARKER_CROSS, AF_MARKER_NONE, AF_MARKER_PLUS, AF_MARKER_POINT, AF_MARKER_SQUARE, AF_MARKER_STAR
export AF_MARKER_TRIANGLE, AF_MAT_BLOCK_DIAG, AF_MAT_CONJ, AF_MAT_CTRANS, AF_MAT_DIAG_UNIT, AF_MAT_LOWER
export AF_MAT_NONE, AF_MAT_ORTHOG, AF_MAT_POSDEF, AF_MAT_SYM, AF_MAT_TRANS, AF_MAT_TRI_DIAG, AF_MAT_UPPER
export AF_MOMENT_FIRST_ORDER, AF_MOMENT_M00, AF_MOMENT_M01, AF_MOMENT_M10, AF_MOMENT_M11, AF_NCC, AF_NORM_EUCLID
export AF_NORM_MATRIX_1, AF_NORM_MATRIX_2, AF_NORM_MATRIX_INF, AF_NORM_MATRIX_L_PQ, AF_NORM_VECTOR_1, AF_NORM_VECTOR_2
export AF_NORM_VECTOR_INF, AF_NORM_VECTOR_P, AF_PAD_SYM, AF_PAD_ZERO, AF_RANDOM_ENGINE_DEFAULT, AF_RANDOM_ENGINE_MERSENNE
export AF_RANDOM_ENGINE_MERSENNE_GP11213, AF_RANDOM_ENGINE_PHILOX, AF_RANDOM_ENGINE_PHILOX_4X32_10, AF_RANDOM_ENGINE_THREEFRY
export AF_RANDOM_ENGINE_THREEFRY_2X32_16, AF_RGB, AF_SAD, AF_SHD, AF_SSD, AF_STORAGE_COO, AF_STORAGE_CSC
export AF_STORAGE_CSR, AF_STORAGE_DENSE, AF_SUCCESS, AF_TOPK_DEFAULT, AF_TOPK_MAX, AF_TOPK_MIN, AF_YCC_2020
export AF_YCC_601, AF_YCC_709, AF_YCbCr, AF_ZNCC, AF_ZSAD, AF_ZSSD, afDevice, afHost, af_array, af_backend
export af_binary_op, af_border_type, af_canny_threshold, af_colormap, af_connectivity, af_conv_domain
export af_conv_mode, af_cspace_t, af_diffusion_eq, af_dtype, af_err, af_features, af_flux_function, af_homography_type
export af_image_format, af_interp_type, af_marker_type, af_mat_prop, af_match_type, af_moment_type, af_norm_type
export af_random_engine, af_random_engine_type, af_someenum_t, af_source, af_storage, af_topk_function
export af_window, af_ycc_std, b8, c32, c64, dim_t, f16, f32, f64, intl, s16, s32, s64, u16, u32, u64, u8, uintl

const dim_t = Clonglong
const intl = Clonglong
const uintl = Culonglong

# begin enum af_err
const af_err = UInt32
const AF_SUCCESS = (UInt32)(0)
const AF_ERR_NO_MEM = (UInt32)(101)
const AF_ERR_DRIVER = (UInt32)(102)
const AF_ERR_RUNTIME = (UInt32)(103)
const AF_ERR_INVALID_ARRAY = (UInt32)(201)
const AF_ERR_ARG = (UInt32)(202)
const AF_ERR_SIZE = (UInt32)(203)
const AF_ERR_TYPE = (UInt32)(204)
const AF_ERR_DIFF_TYPE = (UInt32)(205)
const AF_ERR_BATCH = (UInt32)(207)
const AF_ERR_DEVICE = (UInt32)(208)
const AF_ERR_NOT_SUPPORTED = (UInt32)(301)
const AF_ERR_NOT_CONFIGURED = (UInt32)(302)
const AF_ERR_NONFREE = (UInt32)(303)
const AF_ERR_NO_DBL = (UInt32)(401)
const AF_ERR_NO_GFX = (UInt32)(402)
const AF_ERR_LOAD_LIB = (UInt32)(501)
const AF_ERR_LOAD_SYM = (UInt32)(502)
const AF_ERR_ARR_BKND_MISMATCH = (UInt32)(503)
const AF_ERR_INTERNAL = (UInt32)(998)
const AF_ERR_UNKNOWN = (UInt32)(999)
# end enum af_err

# begin enum af_dtype
const af_dtype = UInt32
const f32 = (UInt32)(0)
const c32 = (UInt32)(1)
const f64 = (UInt32)(2)
const c64 = (UInt32)(3)
const b8 = (UInt32)(4)
const s32 = (UInt32)(5)
const u32 = (UInt32)(6)
const u8 = (UInt32)(7)
const s64 = (UInt32)(8)
const u64 = (UInt32)(9)
const s16 = (UInt32)(10)
const u16 = (UInt32)(11)
const f16 = (UInt32)(12)
# end enum af_dtype

# begin enum af_source
const af_source = UInt32
const afDevice = (UInt32)(0)
const afHost = (UInt32)(1)
# end enum af_source

const af_array = Ptr{Nothing}

# begin enum af_interp_type
const af_interp_type = UInt32
const AF_INTERP_NEAREST = (UInt32)(0)
const AF_INTERP_LINEAR = (UInt32)(1)
const AF_INTERP_BILINEAR = (UInt32)(2)
const AF_INTERP_CUBIC = (UInt32)(3)
const AF_INTERP_LOWER = (UInt32)(4)
const AF_INTERP_LINEAR_COSINE = (UInt32)(5)
const AF_INTERP_BILINEAR_COSINE = (UInt32)(6)
const AF_INTERP_BICUBIC = (UInt32)(7)
const AF_INTERP_CUBIC_SPLINE = (UInt32)(8)
const AF_INTERP_BICUBIC_SPLINE = (UInt32)(9)
# end enum af_interp_type

# begin enum af_border_type
const af_border_type = UInt32
const AF_PAD_ZERO = (UInt32)(0)
const AF_PAD_SYM = (UInt32)(1)
# end enum af_border_type

# begin enum af_connectivity
const af_connectivity = UInt32
const AF_CONNECTIVITY_4 = (UInt32)(4)
const AF_CONNECTIVITY_8 = (UInt32)(8)
# end enum af_connectivity

# begin enum af_conv_mode
const af_conv_mode = UInt32
const AF_CONV_DEFAULT = (UInt32)(0)
const AF_CONV_EXPAND = (UInt32)(1)
# end enum af_conv_mode

# begin enum af_conv_domain
const af_conv_domain = UInt32
const AF_CONV_AUTO = (UInt32)(0)
const AF_CONV_SPATIAL = (UInt32)(1)
const AF_CONV_FREQ = (UInt32)(2)
# end enum af_conv_domain

# begin enum af_match_type
const af_match_type = UInt32
const AF_SAD = (UInt32)(0)
const AF_ZSAD = (UInt32)(1)
const AF_LSAD = (UInt32)(2)
const AF_SSD = (UInt32)(3)
const AF_ZSSD = (UInt32)(4)
const AF_LSSD = (UInt32)(5)
const AF_NCC = (UInt32)(6)
const AF_ZNCC = (UInt32)(7)
const AF_SHD = (UInt32)(8)
# end enum af_match_type

# begin enum af_ycc_std
const af_ycc_std = UInt32
const AF_YCC_601 = (UInt32)(601)
const AF_YCC_709 = (UInt32)(709)
const AF_YCC_2020 = (UInt32)(2020)
# end enum af_ycc_std

# begin enum af_cspace_t
const af_cspace_t = UInt32
const AF_GRAY = (UInt32)(0)
const AF_RGB = (UInt32)(1)
const AF_HSV = (UInt32)(2)
const AF_YCbCr = (UInt32)(3)
# end enum af_cspace_t

# begin enum af_mat_prop
const af_mat_prop = UInt32
const AF_MAT_NONE = (UInt32)(0)
const AF_MAT_TRANS = (UInt32)(1)
const AF_MAT_CTRANS = (UInt32)(2)
const AF_MAT_CONJ = (UInt32)(4)
const AF_MAT_UPPER = (UInt32)(32)
const AF_MAT_LOWER = (UInt32)(64)
const AF_MAT_DIAG_UNIT = (UInt32)(128)
const AF_MAT_SYM = (UInt32)(512)
const AF_MAT_POSDEF = (UInt32)(1024)
const AF_MAT_ORTHOG = (UInt32)(2048)
const AF_MAT_TRI_DIAG = (UInt32)(4096)
const AF_MAT_BLOCK_DIAG = (UInt32)(8192)
# end enum af_mat_prop

# begin enum af_norm_type
const af_norm_type = UInt32
const AF_NORM_VECTOR_1 = (UInt32)(0)
const AF_NORM_VECTOR_INF = (UInt32)(1)
const AF_NORM_VECTOR_2 = (UInt32)(2)
const AF_NORM_VECTOR_P = (UInt32)(3)
const AF_NORM_MATRIX_1 = (UInt32)(4)
const AF_NORM_MATRIX_INF = (UInt32)(5)
const AF_NORM_MATRIX_2 = (UInt32)(6)
const AF_NORM_MATRIX_L_PQ = (UInt32)(7)
const AF_NORM_EUCLID = (UInt32)(2)
# end enum af_norm_type

# begin enum af_image_format
const af_image_format = UInt32
const AF_FIF_BMP = (UInt32)(0)
const AF_FIF_ICO = (UInt32)(1)
const AF_FIF_JPEG = (UInt32)(2)
const AF_FIF_JNG = (UInt32)(3)
const AF_FIF_PNG = (UInt32)(13)
const AF_FIF_PPM = (UInt32)(14)
const AF_FIF_PPMRAW = (UInt32)(15)
const AF_FIF_TIFF = (UInt32)(18)
const AF_FIF_PSD = (UInt32)(20)
const AF_FIF_HDR = (UInt32)(26)
const AF_FIF_EXR = (UInt32)(29)
const AF_FIF_JP2 = (UInt32)(31)
const AF_FIF_RAW = (UInt32)(34)
# end enum af_image_format

# begin enum af_moment_type
const af_moment_type = UInt32
const AF_MOMENT_M00 = (UInt32)(1)
const AF_MOMENT_M01 = (UInt32)(2)
const AF_MOMENT_M10 = (UInt32)(4)
const AF_MOMENT_M11 = (UInt32)(8)
const AF_MOMENT_FIRST_ORDER = (UInt32)(15)
# end enum af_moment_type

# begin enum af_homography_type
const af_homography_type = UInt32
const AF_HOMOGRAPHY_RANSAC = (UInt32)(0)
const AF_HOMOGRAPHY_LMEDS = (UInt32)(1)
# end enum af_homography_type

# begin enum af_backend
const af_backend = UInt32
const AF_BACKEND_DEFAULT = (UInt32)(0)
const AF_BACKEND_CPU = (UInt32)(1)
const AF_BACKEND_CUDA = (UInt32)(2)
const AF_BACKEND_OPENCL = (UInt32)(4)
# end enum af_backend

# begin enum af_someenum_t
const af_someenum_t = UInt32
const AF_ID = (UInt32)(0)
# end enum af_someenum_t

# begin enum af_binary_op
const af_binary_op = UInt32
const AF_BINARY_ADD = (UInt32)(0)
const AF_BINARY_MUL = (UInt32)(1)
const AF_BINARY_MIN = (UInt32)(2)
const AF_BINARY_MAX = (UInt32)(3)
# end enum af_binary_op

# begin enum af_random_engine_type
const af_random_engine_type = UInt32
const AF_RANDOM_ENGINE_PHILOX_4X32_10 = (UInt32)(100)
const AF_RANDOM_ENGINE_THREEFRY_2X32_16 = (UInt32)(200)
const AF_RANDOM_ENGINE_MERSENNE_GP11213 = (UInt32)(300)
const AF_RANDOM_ENGINE_PHILOX = (UInt32)(100)
const AF_RANDOM_ENGINE_THREEFRY = (UInt32)(200)
const AF_RANDOM_ENGINE_MERSENNE = (UInt32)(300)
const AF_RANDOM_ENGINE_DEFAULT = (UInt32)(100)
# end enum af_random_engine_type

# begin enum af_colormap
const af_colormap = UInt32
const AF_COLORMAP_DEFAULT = (UInt32)(0)
const AF_COLORMAP_SPECTRUM = (UInt32)(1)
const AF_COLORMAP_COLORS = (UInt32)(2)
const AF_COLORMAP_RED = (UInt32)(3)
const AF_COLORMAP_MOOD = (UInt32)(4)
const AF_COLORMAP_HEAT = (UInt32)(5)
const AF_COLORMAP_BLUE = (UInt32)(6)
const AF_COLORMAP_INFERNO = (UInt32)(7)
const AF_COLORMAP_MAGMA = (UInt32)(8)
const AF_COLORMAP_PLASMA = (UInt32)(9)
const AF_COLORMAP_VIRIDIS = (UInt32)(10)
# end enum af_colormap

# begin enum af_marker_type
const af_marker_type = UInt32
const AF_MARKER_NONE = (UInt32)(0)
const AF_MARKER_POINT = (UInt32)(1)
const AF_MARKER_CIRCLE = (UInt32)(2)
const AF_MARKER_SQUARE = (UInt32)(3)
const AF_MARKER_TRIANGLE = (UInt32)(4)
const AF_MARKER_CROSS = (UInt32)(5)
const AF_MARKER_PLUS = (UInt32)(6)
const AF_MARKER_STAR = (UInt32)(7)
# end enum af_marker_type

# begin enum af_canny_threshold
const af_canny_threshold = UInt32
const AF_CANNY_THRESHOLD_MANUAL = (UInt32)(0)
const AF_CANNY_THRESHOLD_AUTO_OTSU = (UInt32)(1)
# end enum af_canny_threshold

# begin enum af_storage
const af_storage = UInt32
const AF_STORAGE_DENSE = (UInt32)(0)
const AF_STORAGE_CSR = (UInt32)(1)
const AF_STORAGE_CSC = (UInt32)(2)
const AF_STORAGE_COO = (UInt32)(3)
# end enum af_storage

# begin enum af_flux_function
const af_flux_function = UInt32
const AF_FLUX_QUADRATIC = (UInt32)(1)
const AF_FLUX_EXPONENTIAL = (UInt32)(2)
const AF_FLUX_DEFAULT = (UInt32)(0)
# end enum af_flux_function

# begin enum af_diffusion_eq
const af_diffusion_eq = UInt32
const AF_DIFFUSION_GRAD = (UInt32)(1)
const AF_DIFFUSION_MCDE = (UInt32)(2)
const AF_DIFFUSION_DEFAULT = (UInt32)(0)
# end enum af_diffusion_eq

# begin enum af_topk_function
const af_topk_function = UInt32
const AF_TOPK_MIN = (UInt32)(1)
const AF_TOPK_MAX = (UInt32)(2)
const AF_TOPK_DEFAULT = (UInt32)(0)
# end enum af_topk_function

struct af_seq
    _begin::Cdouble
    _end::Cdouble
    step::Cdouble
end

struct af_index_t
    idx::Nothing
    isSeq::Bool
    isBatch::Bool
end

struct af_cfloat
    real::Cfloat
    imag::Cfloat
end

struct af_cdouble
    real::Cdouble
    imag::Cdouble
end

const af_features = Ptr{Nothing}
const af_window = Culonglong

struct af_cell
    row::Cint
    col::Cint
    title::Cstring
    cmap::af_colormap
end

const af_random_engine = Ptr{Nothing}


text:
@testset "Nivre 08" begin

    # See paper:
    # Nivre, 2008: "Algorithms for Deterministic Incremental Dependency Parsing"
    # https://www.aclweb.org/anthology/J08-4003.pdf

    tree_fig1, tree_fig2 = test_treebank("nivre08.conll")

    @testset "Figure 6" begin
        # arc-eager parse for "economic news had little effect on financial markets"
        table = [
            # Transition      Stack        Buffer
            (Shift(),         [0,1],       2:9),
            (LeftArc("NMOD"), [0],         2:9),
            (Shift(),         [0,2],       3:9),
            (LeftArc("SUBJ"), [0],         3:9),
            (RightArc("ROOT"),[0,3],       4:9),
            (Shift(),         [0,3,4],     5:9),
            (LeftArc("NMOD"), [0,3],       5:9),
            (RightArc("OBJ"), [0,3,5],     6:9),
            (RightArc("NMOD"),[0,3,5,6],   7:9),
            (Shift(),         [0,3,5,6,7], 8:9),
            (LeftArc("NMOD"), [0,3,5,6],   8:9),
            (RightArc("PMOD"),[0,3,5,6,8], [9]),
            (Reduce(),        [0,3,5,6],   [9]),
            (Reduce(),        [0,3,5],     [9]),
            (Reduce(),        [0,3],       [9]),
            (RightArc("P"),   [0,3,9],     [])]

        cfg = initconfig(ArcEager(), tree_fig2)
        o = cfg -> static_oracle(cfg, tree_fig2, typed)
        @test stack(cfg) == [0]
        @test buffer(cfg) == 1:9
        for (t, stk, buf) in table
            @test !isfinal(cfg)
            @test t == o(cfg)
            cfg = t(cfg)
            @test stack(cfg) == stk
            @test buffer(cfg) == buf
        end
        @test isfinal(cfg)
        result = DependencyTree(tokens(cfg))
        @test result == tree_fig2
    end

    @testset "Figure 8" begin
        @test !projective_only(ListBasedNonProjective())
        table = [
            # Transition        L1       L2     buffer
            (Shift(),           [0,1],   [],    2:8),
            (RightArc("Atr"),   [0],     [1],   2:8),
            (Shift(),           [0,1,2], [],    3:8),
            (NoArc(),           [0,1],   [2],   3:8),
            (NoArc(),           [0],     [1,2], 3:8),
            (RightArc("Pred"),  [],      0:2,   3:8),
            (Shift(),           0:3,     [],    4:8),
            (Shift(),           0:4,     [],    5:8),
            (LeftArc("AuxZ"),   0:3,     [4],   5:8),
            (RightArc("Sb"),    0:2,     3:4,   5:8),
            (NoArc(),           0:1,     2:4,   5:8),
            (LeftArc("AuxP"),   [0],     1:4,   5:8),
            (Shift(),           0:5,     [],    6:8),
            (NoArc(),           0:4,     [5],   6:8),
            (NoArc(),           0:3,     4:5,   6:8),
            (RightArc("AuxP"),  0:2,     3:5,   6:8),
            (Shift(),           0:6,     [],    7:8),
            (RightArc("Adv"),   0:5,     [6],   7:8),
            (Shift(),           0:7,     [],    [8]),
            (NoArc(),           0:6,     [7],   [8]),
            (NoArc(),           0:5,     6:7,   [8]),
            (NoArc(),           0:4,     5:7,   [8]),
            (NoArc(),           0:3,     4:7,   [8]),
            (NoArc(),           0:2,     3:7,   [8]),
            (NoArc(),           0:1,     2:7,   [8]),
            (NoArc(),           [0],     1:7,   [8]),
            (RightArc("AuxK"),  [],      0:7,   [8]),
            (Shift(),           0:8,     [],    [])]

        o(cfg) = static_oracle(cfg, tree_fig1, typed)
        cfg = initconfig(ListBasedNonProjective(), tree_fig1)
        for (t, l1, l2, buf) in table
            @test !isfinal(cfg)
            @test o(cfg) == t
            cfg = t(cfg)
            @test cfg.λ1 == l1
            @test cfg.λ2 == l2
            @test buffer(cfg) == buf
        end
        @test isfinal(cfg)
        @test replace(showstr(cfg), r"\s+"=>"") ==
            """ListBasedNonProjectiveConfig([0,1,2,3,4,5,6,7,8],[],[])
             1	Z	5
             2	nich	1
             3	je	0
             4	jen	5
             5	jedna	3
             6	na	3
             7	kvalitu	6
             8	.	0""" |> x->replace(x, r"\s"=>"")

        result = DependencyTree(cfg.A, check_single_head=false)
        @test result == tree_fig1

        oracle = StaticOracle(ListBasedNonProjective(), arc=typed)
        pairs = DependencyTrees.xys(oracle, tree_fig1)
        @test last.(pairs) == first.(table)

        cfg1 = initconfig(ListBasedNonProjective(), CoNLLU, [w.form for w in tree_fig1])
        cfg2 = initconfig(ListBasedNonProjective(), tree_fig1)
        @test cfg1 != cfg2 # cfg2 knows about the gold labels, cfg1 doesn't
        @test cfg1.λ1 == cfg2.λ1 && cfg1.λ2 == cfg2.λ2 && cfg1.β == cfg2.β
        @test [token(cfg, 1)] == tokens(cfg, [1])
    end
end


text:
using Revise
using AdFem
using PyPlot

mmesh = Mesh(joinpath(PDATA, "twoholes.stl"))
left = bcnode((x,y)->x<1e-5, mmesh)
right = bcedge((x1,y1,x2,y2)->(x1>0.049-1e-5) && (x2>0.049-1e-5), mmesh)

t1 = eval_f_on_boundary_edge((x,y)->1.0e-4, right, mmesh)
t2 = eval_f_on_boundary_edge((x,y)->0.0, right, mmesh)
rhs = compute_fem_traction_term(t1, t2, right, mmesh)

D = diagm(0=>[1,1,0.5])
K = constant(compute_fem_stiffness_matrix(D, mmesh))

bdval = [eval_f_on_boundary_node((x,y)->0.0, left, mmesh);
        eval_f_on_boundary_node((x,y)->0.0, left, mmesh)]
DOF = [left;left .+ mmesh.ndof]
K, rhs = impose_Dirichlet_boundary_conditions(K, rhs, DOF, bdval)
u = K\rhs
sess = Session(); init(sess)
S = run(sess, u)

# close("all")
# visualize_displacement(S*1e3, mmesh)
# savefig("test.png")

close("all")
figure(figsize=(10,20))
subplot(311)
visualize_scalar_on_fem_points(S[1:mmesh.nnode], mmesh)
subplot(312)
visualize_scalar_on_fem_points(S[1+mmesh.ndof:mmesh.nnode + mmesh.ndof], mmesh)
subplot(313)
visualize_von_mises_stress(D, S, mmesh)
savefig("test.png")


text:
"""
    StepSize

A composite type for the step-size of an [`AbstractRungeKuttaSolver`](@ref).

# Constructors
```julia
StepSize(; h)
```

# Arguments
- `h :: Real` : step-size.

# Functions
- [`show`   ](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
mutable struct StepSize{h_T}
    h::h_T
end

############################################################################################
#                                         FUNCTIONS                                        #
############################################################################################

"""
    show(io::IO, stepsize::StepSize)

prints a full description of `stepsize` and its contents to a stream `io`.
"""
Base.show(io::IO, stepsize::StepSize) = NSDEBase._show(io, stepsize)

"""
    summary(io::IO, stepsize::StepSize)

prints a brief description of `stepsize` to a stream `io`.
"""
Base.summary(io::IO, stepsize::StepSize) = NSDEBase._summary(io, stepsize)