text:
import LibGEOS
import GeoJSON
import Base.convert
using DataFrames
using Distributions
import CSV
import JSON

import GeoInterface

const SchDistr = Int

const SALE_PRICE = Symbol("SALE PRICE")
const SQFT = Symbol("GROSS SQUARE FEET")
const BUILDING_CLASS_AT_TIME_OF_SALE = Symbol("BUILDING CLASS AT TIME OF SALE")
const BUILDING_CLASS_CATEGORY = Symbol("BUILDING CLASS CATEGORY")
const TAX_CLASS_AT_TIME_OF_SALE = Symbol("TAX CLASS AT TIME OF SALE")

function read_distr_shapes(; data_dir="NYC_data")
    json=GeoJSON.parsefile(joinpath(data_dir, "nysd_16c", "nysd.json"))
    schdistr_shape_dict = Dict{SchDistr, GeoRDD.RegionType}()
    for feature in json.features
        schdistr = convert(SchDistr, feature.properties["SchoolDist"])
        shape = convert(GeoRDD.RegionType, feature.geometry)
        schdistr_shape_dict[schdistr] = shape
    end
    return schdistr_shape_dict
end

function read_processed_sales(; data_dir="NYC_data")
    NYC_sales=CSV.read(joinpath(data_dir, "processed", "NYC_sales.csv"),
                       DataFrame,
                       types=Dict("TAX CLASS AT PRESENT" => Union{Missings.Missing, String},
                                  "TAX CLASS AT TIME OF SALE" => Union{Missings.Missing, String}),
                       strings=:raw,
                       allowmissing=:all,
                      )
    nyc_schdistrs = CategoricalVector(NYC_sales[:SchDistr])
    # schd_strings = [ismissing(sd)?missing:dec(sd,2) for sd in  nyc_schdistrs]
    # str_schdistrs = CategoricalVector(schd_strings)
    NYC_sales[:SchDistr] = nyc_schdistrs
    # categorical variables
    categorical!(NYC_sales, Symbol[
        "BOROUGH",
        "BUILDING CLASS CATEGORY",
        "BUILDING CLASS AT TIME OF SALE",
        "ZIP CODE",
        "TAX CLASS AT PRESENT",
        "TAX CLASS AT TIME OF SALE",
        "NEIGHBORHOOD",
        ])
    sort!(NYC_sales, :SchDistr)
    return NYC_sales
end

const DWELLINGS_DICT = Dict(
    "01  ONE FAMILY DWELLINGS"=>true,
    "02  TWO FAMILY DWELLINGS"=>true,
    "03  THREE FAMILY DWELLINGS"=>true,
    "04  TAX CLASS 1 CONDOS"=>false,
    "05  TAX CLASS 1 VACANT LAND"=>false,
    "06  TAX CLASS 1 - OTHER"=>false,
    "07  RENTALS - WALKUP APARTMENTS"=>false,
    "08  RENTALS - ELEVATOR APARTMENTS"=>false,
    "09  COOPS - WALKUP APARTMENTS"=>false,
    "10  COOPS - ELEVATOR APARTMENTS"=>false,
    "11  SPECIAL CONDO BILLING LOTS"=>false,
    "11A CONDO-RENTALS"=>false,
    "12  CONDOS - WALKUP APARTMENTS"=>false, # why are these duplicated?
    "13  CONDOS - ELEVATOR APARTMENTS"=>false,
    "14  RENTALS - 4-10 UNIT"=>false,
    "15  CONDOS - 2-10 UNIT RESIDENTIAL"=>false,
    "16  CONDOS - 2-10 UNIT WITH COMMERCIAL UNIT"=>false,
    "17  CONDO COOPS"=>false,
    "18  TAX CLASS 3 - UNTILITY PROPERTIES"=>false,
    "21  OFFICE BUILDINGS"=>false,
    "22  STORE BUILDINGS"=>false,
    "23  LOFT BUILDINGS"=>false,
    "25  LUXURY HOTELS"=>false,
    "26  OTHER HOTELS"=>false,
    "27  FACTORIES"=>false,
    "28  COMMERCIAL CONDOS"=>false,
    "29  COMMERCIAL GARAGES"=>false,
    "30  WAREHOUSES"=>false,
    "31  COMMERCIAL VACANT LAND"=>false,
    "32  HOSPITAL AND HEALTH FACILITIES"=>false,
    "33  EDUCATIONAL FACILITIES"=>false,
    "34  THEATRES"=>false,
    "35  INDOOR PUBLIC AND CULTURAL FACILITIES"=>false,
    "36  OUTDOOR RECREATIONAL FACILITIES"=>false,
    "37  RELIGIOUS FACILITIES"=>false,
    "38  ASYLUMS AND HOMES"=>false,
    "41  TAX CLASS 4 - OTHER"=>false,
    "42  CONDO CULTURAL/MEDICAL/EDUCATIONAL/ETC"=>false,
    "43  CONDO OFFICE BUILDINGS"=>false,
    "44  CONDO PARKING"=>false,
    "45  CONDO HOTELS"=>false,
    "46  CONDO STORE BUILDINGS"=>false,
    "47  CONDO NON-BUSINESS STORAGE"=>false,
    "48  CONDO TERRACES/GARDENS/CABANAS"=>false,
    "49  CONDO WAREHOUSES/FACTORY/INDUS"=>false,
    )
const RESIDENTIAL_DICT = Dict(
    "01  ONE FAMILY DWELLINGS"=>true,
    "02  TWO FAMILY DWELLINGS"=>true,
    "03  THREE FAMILY DWELLINGS"=>true,
    "04  TAX CLASS 1 CONDOS"=>true,
    "05  TAX CLASS 1 VACANT LAND"=>false,
    "06  TAX CLASS 1 - OTHER"=>false,
    "07  RENTALS - WALKUP APARTMENTS"=>false,
    "08  RENTALS - ELEVATOR APARTMENTS"=>false,
    "09  COOPS - WALKUP APARTMENTS"=>true,
    "10  COOPS - ELEVATOR APARTMENTS"=>true,
    "11  SPECIAL CONDO BILLING LOTS"=>false,
    "11A CONDO-RENTALS"=>false,
    "12  CONDOS - WALKUP APARTMENTS"=>true, # why are these duplicated?
    "13  CONDOS - ELEVATOR APARTMENTS"=>true,
    "14  RENTALS - 4-10 UNIT"=>false,
    "15  CONDOS - 2-10 UNIT RESIDENTIAL"=>true,
    "16  CONDOS - 2-10 UNIT WITH COMMERCIAL UNIT"=>false,
    "17  CONDO COOPS"=>true,
    "18  TAX CLASS 3 - UNTILITY PROPERTIES"=>false,
    "21  OFFICE BUILDINGS"=>false,
    "22  STORE BUILDINGS"=>false,
    "23  LOFT BUILDINGS"=>false,
    "25  LUXURY HOTELS"=>false,
    "26  OTHER HOTELS"=>false,
    "27  FACTORIES"=>false,
    "28  COMMERCIAL CONDOS"=>false,
    "29  COMMERCIAL GARAGES"=>false,
    "30  WAREHOUSES"=>false,
    "31  COMMERCIAL VACANT LAND"=>false,
    "32  HOSPITAL AND HEALTH FACILITIES"=>false,
    "33  EDUCATIONAL FACILITIES"=>false,
    "34  THEATRES"=>false,
    "35  INDOOR PUBLIC AND CULTURAL FACILITIES"=>false,
    "36  OUTDOOR RECREATIONAL FACILITIES"=>false,
    "37  RELIGIOUS FACILITIES"=>false,
    "38  ASYLUMS AND HOMES"=>false,
    "41  TAX CLASS 4 - OTHER"=>false,
    "42  CONDO CULTURAL/MEDICAL/EDUCATIONAL/ETC"=>false,
    "43  CONDO OFFICE BUILDINGS"=>false,
    "44  CONDO PARKING"=>false,
    "45  CONDO HOTELS"=>false,
    "46  CONDO STORE BUILDINGS"=>false,
    "47  CONDO NON-BUSINESS STORAGE"=>false,
    "48  CONDO TERRACES/GARDENS/CABANAS"=>false,
    "49  CONDO WAREHOUSES/FACTORY/INDUS"=>false,
    )

function filter_sales(NYC_sales::DataFrame)

    logSalePricePerSQFT = log.(NYC_sales[SALE_PRICE]) .- log.(NYC_sales[SQFT])
    believable = zeros(Bool, size(NYC_sales,1))
    removed_reason = Dict(
        :residential => 0,
        :noprice => 0,
        :nosqft => 0,
        :geocode => 0,
        :toosmall => 0,
        :outlier => 0,
       )
    for i in 1:size(NYC_sales,1)
        if ismissing(NYC_sales[i,BUILDING_CLASS_AT_TIME_OF_SALE])
            # remove sales with missing covariates
            removed_reason[:residential] += 1
        elseif ismissing(NYC_sales[i,BUILDING_CLASS_CATEGORY])
            # remove sales with missing covariates
            removed_reason[:residential] += 1
        elseif !DWELLINGS_DICT[NYC_sales[i,BUILDING_CLASS_CATEGORY]]
            # remove non-residential properties
            # in fact, remove things that aren't dwellings (houses?)
            removed_reason[:residential] += 1
        elseif ismissing(NYC_sales[i,SALE_PRICE])
            # remove sales without a sale price
            removed_reason[:noprice] += 1
        elseif ismissing(NYC_sales[i,SQFT])
            # remove sales with missing GROSS SQUARE FEET information
            removed_reason[:nosqft] += 1
        elseif ismissing(NYC_sales[i,:XCoord])
            # remove properties with failed geocoding
            removed_reason[:geocode] += 1
        elseif ismissing(NYC_sales[i,:YCoord])
            # remove properties with failed geocoding
            removed_reason[:geocode] += 1
        elseif ismissing(NYC_sales[i,:SchDistr])
            # remove properties without a school district
            removed_reason[:geocode] += 1
        # elseif ismissing(NYC_sales[i,TAX_CLASS_AT_TIME_OF_SALE])
            # # remove sales with missing covariates
            # removed_reason[7] += 1
        elseif NYC_sales[i,SQFT] < 100.0
            # remove properties smaller to 100sqft (seems unlikely to be real)
            removed_reason[:toosmall] += 1
        elseif logSalePricePerSQFT[i] < 3.0
            # that's too cheap (remove outliers)
            removed_reason[:outlier] += 1
        elseif logSalePricePerSQFT[i] > 8.0
            # that's too expensive (remove outliers)
            removed_reason[:outlier] += 1
        else
            # otherwise keep
            believable[i] = true
        end
    end

    filtered = copy(NYC_sales[believable,:])
    filtered[:logSalePricePerSQFT] = convert.(Float64, logSalePricePerSQFT[believable])
    for col in (:XCoord, :YCoord)
        filtered[col] = convert.(Float64, filtered[col])
    end
    return Dict(
        :filtered => filtered,
        :believable => believable,
        :removed_reason => removed_reason,
        )
end

# function sales_dicts(NYC_sales::DataFrame)
    # schdistrs = sort(NYC_sales[:SchDistr].pool.levels)
    # schdistr_indices = Dict{SchDistr,Vector{Int}}()
    # schdistrs_col = NYC_sales[:SchDistr]
    # for distr in schdistrs
        # indices = find(schdistrs_col.refs .== find(schdistrs_col.pool.index .== distr)[1])
        # nobsv_schdistr = length(indices)
        # @printf("District %d has %d sales\n", distr, nobsv_schdistr)
        # schdistr_indices[distr] = indices
    # end
    # Y_dict=Dict{SchDistr, Vector{Float64}}()
    # X_dict=Dict{SchDistr, Array{Float64,2}}()
    # for distr in schdistrs
        # Y_dict[distr] = NYC_sales[schdistr_indices[distr], :logSalePricePerSQFT]
        # X_dict[distr] = NYC_sales[schdistr_indices[distr],[:XCoord, :YCoord]]
    # end
    # return Dict(
            # :schdistr_indices => schdistr_indices,
            # :schdistrs => schdistrs,
            # :X_dict => X_dict,
            # :Y_dict => Y_dict
            # )
# end

function table_τpair(τpost_pairs::Dict{Tuple{NYC.SchDistr,NYC.SchDistr},Normal}, include_distr::Vector{Int})
    τ_pair_nested = Dict{NYC.SchDistr, Dict{NYC.SchDistr, Normal}}()
    for (distrA, distrB) in keys(τpost_pairs)
        τpost = τpost_pairs[distrA, distrB]
        if distrA ∉ include_distr
            continue
        end
        if distrB ∉ include_distr
            continue
        end
        if distrA ∉ keys(τ_pair_nested)
            τ_pair_nested[distrA] = Dict{NYC.SchDistr, Normal}()
        end
        τ_pair_nested[distrA][distrB] = τpost
    end
    ;

    for distrA in sort(collect(keys(τ_pair_nested)))
        if isempty(τ_pair_nested[distrA])
            continue
        end
        @printf("\\( \\mathbf{%2d} \\)", distrA)
        for distrB in sort(collect(keys(τ_pair_nested[distrA])))
            τ = τ_pair_nested[distrA][distrB]
            @printf("& \\( \\mathbf{%2d:}~%+.2f \\pm %.2f \\)", distrB, -mean(τ), std(τ))
        end
        print("\\\\ \n")
    end
end


text:
export getHybridKaczmarzPrecond,getHybridKaczmarz, getHybridKaczmarzParam
const parRelax_lib  = abspath(joinpath(splitdir(Base.source_path())[1],"../..","deps","builds","parRelax"))

const ArrIdxsType = UInt32;

mutable struct hybridKaczmarz{VAL,IND}
	numDomains 	:: Array{Int64,1}
	invDiag 	:: Array
	numCores	:: Int64
	omega_damp	:: Float64
	ArrIdxs 	:: Array{ArrIdxsType,2}
	precond 	:: Function
	numit 		:: Int64
end

mutable struct hybridKaczmarzParams
	numDomains 	:: Array{Int64,1}
	numCores	:: Int64
	omega_damp	:: Float64
	numit 		:: Int64
end

import Base.copy
function copy(p::hybridKaczmarzParams)
	return hybridKaczmarzParams(p.numDomains,p.numCores,p.omega_damp,p.numit);
end

function getHybridKaczmarzParam(numDomains::Array{Int64,1},omega_damp::Float64,numCores::Int64,numit::Int64)
	return hybridKaczmarzParams(numDomains,numCores,omega_damp,numit);
end

function getHybridKaczmarz(VAL::Type,IND::Type,AT::SparseMatrixCSC,Mesh::RegularMesh, numDomains::Array{Int64,1},
						getIndicesOfCell::Function,omega_damp::Float64,numCores::Int64,numit::Int64)
if prod(numDomains) < numCores
	println("*** WARNING: getHybridKaczmarz: numDomains < numCores. ***");
end
invDiag = convert(Array{VAL,1},omega_damp./(vec(sum(conj(AT).*AT,dims=1))));
ArrIdxs = getIndicesOfCellsArray(Mesh, zeros(Int64,size(numDomains)),numDomains,getIndicesOfCell);
return hybridKaczmarz{VAL,IND}(numDomains,invDiag,numCores,omega_damp,ArrIdxs,identity,numit);
end

function getHybridKaczmarzPrecond(param::hybridKaczmarz{VAL,IND},AT::SparseMatrixCSC{VAL,IND},nrhs::Int64) where {VAL,IND}
# void applyHybridKaczmarz(spIndType *rowptr , spValType *valA ,spIndType *colA, long long numDomains, long long domainLength ,
						# unsigned int *ArrIdxs, spValType *x, spValType *b, spValType *invD,long long numit, long long numCores){
x = zeros(VAL,size(AT,2),nrhs);
numDomains = prod(param.numDomains);
param.precond = (r)->(x[:].=0.0;applyHybridKaczmarz(param::hybridKaczmarz{VAL,IND},AT::SparseMatrixCSC{VAL,IND},r::Array{VAL},x::Array{VAL},numDomains);return x;);
return param.precond;
end

function applyHybridKaczmarz(param::hybridKaczmarz{Float64,Int64},AT::SparseMatrixCSC{Float64,Int64},r::Array{Float64},x::Array{Float64},numDomains)
	ccall((:applyHybridKaczmarz_FP64_INT64,parRelax_lib),Nothing,(Ptr{Int64},Ptr{Float64},Ptr{Int64},Int64,Int64,Ptr{ArrIdxsType},Ptr{Float64},Ptr{Float64},Int64,Int64,Ptr{Float64},Int64,Int64,),
		AT.colptr,AT.nzval,AT.rowval,numDomains,size(param.ArrIdxs,1),param.ArrIdxs,x,r,size(r,2),size(r,1),param.invDiag,param.numit,param.numCores);
end

function applyHybridKaczmarz(param::hybridKaczmarz{Float32,Int64},AT::SparseMatrixCSC{Float32,Int64},r::Array{Float32},x::Array{Float32},numDomains)
	ccall((:applyHybridKaczmarz_FP32_INT64,parRelax_lib),Nothing,(Ptr{Int64},Ptr{Float32},Ptr{Int64},Int64,Int64,Ptr{ArrIdxsType},Ptr{Float32},Ptr{Float32},Int64,Int64,Ptr{Float32},Int64,Int64,),
		AT.colptr,AT.nzval,AT.rowval,numDomains,size(param.ArrIdxs,1),param.ArrIdxs,x,r,size(r,2),size(r,1),param.invDiag,param.numit,param.numCores);
end

function applyHybridKaczmarz(param::hybridKaczmarz{ComplexF64,Int64},AT::SparseMatrixCSC{ComplexF64,Int64},r::Array{ComplexF64},x::Array{ComplexF64},numDomains)
	ccall((:applyHybridKaczmarz_CFP64_INT64,parRelax_lib),Nothing,(Ptr{Int64},Ptr{ComplexF64},Ptr{Int64},Int64,Int64,Ptr{ArrIdxsType},Ptr{ComplexF64},Ptr{ComplexF64},Int64,Int64,Ptr{ComplexF64},Int64,Int64,),
		AT.colptr,AT.nzval,AT.rowval,numDomains,size(param.ArrIdxs,1),param.ArrIdxs,x,r,size(r,2),size(r,1),param.invDiag,param.numit,param.numCores);
end

function applyHybridKaczmarz(param::hybridKaczmarz{ComplexF32,Int64},AT::SparseMatrixCSC{ComplexF32,Int64},r::Array{ComplexF32},x::Array{ComplexF32},numDomains)
	ccall((:applyHybridKaczmarz_CFP32_INT64,parRelax_lib),Nothing,(Ptr{Int64},Ptr{ComplexF32},Ptr{Int64},Int64,Int64,Ptr{ArrIdxsType},Ptr{ComplexF32},Ptr{ComplexF32},Int64,Int64,Ptr{ComplexF32},Int64,Int64,),
		AT.colptr,AT.nzval,AT.rowval,numDomains,size(param.ArrIdxs,1),param.ArrIdxs,x,r,size(r,2),size(r,1),param.invDiag,param.numit,param.numCores);
end

# function applyHybridKaczmarz(param::hybridKaczmarz{VAL,IND},AT::SparseMatrixCSC{VAL,IND},b::Array{VAL},x::Array{VAL},numDomains)
	# ccall((:applyHybridKaczmarz,parRelax_lib),Nothing,(Ptr{Int64},Ptr{spValType},Ptr{Int64},Int64,Int64,Ptr{ArrIdxsType},Ptr{spValType},Ptr{spValType},Ptr{spValType},Int64,Int64,),
		# AT.colptr,AT.nzval,AT.rowval,numDomains,size(ArrIdxs,1),param.ArrIdxs,x,r,param.invDiag,param.numit,param.numCores);
# end

ext: jl


text:
# This file was generated, do not modify it. # hide
using DelimitedFiles

function cnt(z,c)
    k::UInt8 = 0
    while (k < 255)
        z = z^2 + c
        if abs2(z)>4.
            return k
        end
        k += 1
    end
    return k
end

function mand(M, N)
    grid = Array{UInt8}(undef, M, N)
    init_z = complex(0.0, 0.0)
    xs = range(-2,2,length=N)
    ys = range(-2,2,length=M)
    grid = cnt.(init_z,complex.(xs',ys))
    return grid
end

function main()
    s = time()
    M = N = 5000
    grid = mand(M, N)
    t = time()
    @show(t - s)
    writedlm("/tmp/resultJ2.txt",grid,",")
end

main()

using BenchmarkTools
@btime mand(5000, 5000);


text:
module JuliaWeBWorK

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

## USeful  docs
## https://www.uleth.ca/sites/default/files/2017/09/problemgenerationwithwebwork.pdf
## https://webwork.maa.org/mediawiki_new/images/a/ab/WeBWorK_Problem_Authoring_Tutorial.pdf
## https://webwork.maa.org/wiki/SampleProblem4
## https://webwork.maa.org/wiki/Display_Macros

#using Markdown
using CommonMark
using Mustache
using Base64
using Random

export Page, QUESTIONS
export randomizer,  randomq, numericq, stringq, yesnoq, plotq, essayq
export radioq, multiplechoiceq
export iframe, label,  hint
export List,  Formula,  Interval
export Plot, File
export @L_str, @L_mstr  # deprecate this
export @q_str           # deprecate this
export @MT_str          # deprecate this
export @mt_str
export @jmt_str         # main macro for authoring allows $interpolation, single \
export numbers_only

#include("markdown-additions.jl")
include("commonmark-additions.jl")
include("answers.jl")
include("page.jl")

end


text:
"""
    UncertainIndexValueDataset{
        IDXTYP<:AbstractUncertainIndexDataset,
        VALSTYP<:AbstractUncertainValueDataset}

A generic dataset type consisting of a set of uncertain `indices` (e.g. time,
depth, order, etc...) and a set of uncertain `values`.

The i-th index is assumed to correspond to the i-th value. For example, if
`data` is an instance of a `UncertainIndexValueDataset`, then

- `data.indices[2]` is the index for the value `data.values[2]`
- `data.values[7]` is the value for the index `data.indices[7]`.
- `data[3]` is an index-value tuple `(data.indices[3], data.values[3])`.

## Fields

- **`indices::T where {T <: AbstractUncertainIndexDataset}`**: The uncertain indices,
    represented by some type of uncertain index dataset.
- **`values::T  where {T <: AbstractUncertainValueDataset}`**: The uncertain values,
    represented by some type of uncertain index dataset.

## Example

```julia
# Simulate some data values measured a specific times.
times = 1:100
values = sin.(0.0:0.1:100.0)

# Assume the data were measured by a device with normally distributed
# measurement uncertainties with fluctuating standard deviations
σ_range = (0.1, 0.7)

uncertain_values = [UncertainValue(Normal, val, rand(Uniform(σ_range...)))
    for val in values]

# Assume the clock used to record the times is uncertain, but with uniformly
# distributed noise that doesn't change through time.
uncertain_times = [UncertainValue(Uniform, t-0.1, t+0.1) for t in times]

# Pair the time-value data. If vectors are provided to the constructor,
# the first will be interpreted as the indices and the second as the values.
data = UncertainIndexValueDataset(uncertain_times, uncertain_values)

# A safer option is to first convert to UncertainIndexDataset and
# UncertainValueDataset, so you don't accidentally mix the indices
# and the values.
uidxs = UncertainIndexDataset(uncertain_times)
uvals = UncertainValueDataset(uncertain_values)

data = UncertainIndexValueDataset(uidxs, uvals)
```
"""
struct UncertainIndexValueDataset{IDXTYP <: AbstractUncertainIndexDataset, VALSTYP <: AbstractUncertainValueDataset} <: AbstractUncertainIndexValueDataset

    """ The indices of the uncertain index-value dataset. """
    indices::IDXTYP

    """ The values of the uncertain index-value dataset. """
    values::VALSTYP

    function UncertainIndexValueDataset(
        indices::UncertainIndexDataset,
        values::UncertainValueDataset)

        IDXTYP = UncertainIndexDataset
        VALSTYP = UncertainValueDataset

        new{IDXTYP, VALSTYP}(
            UncertainIndexDataset(indices),
            UncertainValueDataset(values)
        )
    end

    function UncertainIndexValueDataset(
        indices::ConstrainedUncertainIndexDataset,
        values::UncertainValueDataset)

        IDXTYP = ConstrainedUncertainIndexDataset
        VALSTYP = UncertainValueDataset

        new{IDXTYP, VALSTYP}(indices, values)
    end

    function UncertainIndexValueDataset(
        indices::UncertainIndexDataset,
        values::ConstrainedUncertainValueDataset)

        IDXTYP = UncertainIndexDataset
        VALSTYP = ConstrainedUncertainValueDataset

        new{IDXTYP, VALSTYP}(indices, values)
    end

    function UncertainIndexValueDataset(
        indices::ConstrainedUncertainIndexDataset,
        values::ConstrainedUncertainValueDataset)

        IDXTYP = ConstrainedUncertainIndexDataset
        VALSTYP = ConstrainedUncertainValueDataset

        new{IDXTYP, VALSTYP}(indices, values)
    end

    function UncertainIndexValueDataset(
            indices::Vector{<:AbstractUncertainValue},
            values::Vector{<:AbstractUncertainValue})

        IDXTYP = UncertainIndexDataset
        VALSTYP = UncertainValueDataset

        new{IDXTYP, VALSTYP}(
            UncertainIndexDataset(indices),
            UncertainValueDataset(values)
        )
    end

    function UncertainIndexValueDataset(
            indices::DT,
            values::Vector{<:AbstractUncertainValue}) where {DT <: AbstractUncertainIndexDataset}

        IDXTYP = DT
        VALSTYP = UncertainValueDataset

        new{IDXTYP, VALSTYP}(indices, UncertainValueDataset(values))
    end

    function UncertainIndexValueDataset(
            indices::Vector{<:AbstractUncertainValue},
            values::DT) where {DT <: AbstractUncertainValueDataset}

        IDXTYP = UncertainIndexDataset
        VALSTYP = DT

        new{IDXTYP, VALSTYP}(UncertainIndexDataset(indices), values)
    end
end

function UncertainIndexValueDataset(x::AbstractArray{T, 1}, y::AbstractArray{T, 1}) where T
    idxs = UncertainIndexDataset(x)
    vals = UncertainValueDataset(y)
    UncertainIndexValueDataset(idxs, vals)
end

Base.length(u::UncertainIndexValueDataset) = length(u.values)
Base.size(u::UncertainIndexValueDataset) = length(u.values)

Base.getindex(u::UncertainIndexValueDataset, i) = (u.indices[i], u.values[i])
Base.getindex(u::UncertainIndexValueDataset, i::AbstractVector) =
    [(u.indices[i], u.values[i]) for i in 1:length(u)]

Base.getindex(u::UncertainIndexValueDataset, i::Colon) =
        [(u.indices[i], u.values[i]) for i in 1:length(u)]

Base.firstindex(u::UncertainIndexValueDataset) = 1
Base.lastindex(u::UncertainIndexValueDataset) = length(u.values)

Base.eachindex(u::UncertainIndexValueDataset) = Base.OneTo(length(u))
Base.iterate(u::UncertainIndexValueDataset, state = 1) = iterate((u.indices, u.values), state)

index(u::UncertainIndexValueDataset, i) = u.indices[i]
value(u::UncertainIndexValueDataset, i) = u.values[i]

export
UncertainIndexValueDataset,
index, value


text:
module AssetRegistry
using SHA
using JSON
using Pidfile

const baseurl = Ref{String}("")
const registry = Dict{String,String}()

withlock(f,lockf) = (lock=mkpidlock(lockf, stale_age=5); f(); close(lock))

# Need to remove the basepath for file database
filekey(key) = key[length(baseurl[])+1:end]

"""

    register("path/to/asset")

    Register an asset. Returns a unique key which is
the unique URL where the asset can be accessed.
For example:

```julia
julia> key = AssetRegistry.register("/Users/ranjan/.julia/v0.6/Tachyons/assets/tachyons.min.css")
"/assetserver/97a47bdda5bd9274ad1a9cd10a0337f3b033a790-tachyons.min.css"
```
"""
function register(path; registry_file = joinpath(homedir(), ".jlassetregistry.json"))
    target = normpath(abspath(expanduser(path)))
    (isfile(target) || isdir(target)) ||
                    error("Asset not found")

    key = getkey(target)
    if haskey(registry, key)
        # We have already registered this in our session
        return key
    end

    # update global registry file

    # touch(registry_file) -- this doesn't work on Azure fs
    if !isfile(registry_file)
        # WARN: may need a lock here
        open(_->nothing, registry_file, "w")
    end

    pidlock = joinpath(homedir(), ".jlassetregistry.lock")

    withlock(pidlock) do
        fkey = filekey(key)
        io = open(registry_file, "r+") # open in read-and-write mode
        # first get existing entries
        prev_registry =  filesize(io) > 0 ? JSON.parse(io) : Dict{String,Tuple{String,Int}}()

        if haskey(prev_registry, fkey)
            prev_registry[fkey] = (target, prev_registry[fkey][2]+1) # increment ref count
        else
            prev_registry[fkey] = (target, 1) # init refcount to 1
        end

        # write the updated registry to file
        seekstart(io)
        JSON.print(io, prev_registry)

        # truncate it to current length
        truncate(io, position(io))
        close(io)
    end

    registry[key] = target # keep this in current process memory

    key
end

function deregister(path; registry_file = joinpath(homedir(), ".jlassetregistry.json"))
    target = normpath(abspath(expanduser(path)))

    key = getkey(target)
    if !haskey(registry, key)
        return
    end

    pop!(registry, key)

    touch(registry_file)

    pidlock = joinpath(homedir(), ".jlassetregistry.lock")
    withlock(pidlock) do
        fkey = filekey(key)
        io = open(registry_file, "r+")

        # get existing information
        prev_registry =  filesize(io) > 0 ? JSON.parse(io) : Dict{String,Tuple{String, Int}}()

        if haskey(prev_registry, fkey)
            val, count = prev_registry[fkey]
            if count == 1
                pop!(prev_registry, fkey)
            else
                prev_registry[fkey] = (val, count-1) # increment ref count
            end
        end

        seekstart(io)
        JSON.print(io, prev_registry)

        # truncate it to current length
        truncate(io, position(io))
        close(io)
    end

    key
end

generatekey(path) =  baseurl[] * "/assetserver/" * bytes2hex(sha1(path)) * "-" * basename(path)

isregistered(path) = normpath(abspath(expanduser(path))) in values(registry)

function getkey(path)
    target = normpath(abspath(expanduser(path)))
    key = findfirst(x -> x == target, registry)
    key === nothing ? generatekey(path) : key
end

function __init__()
    baseurl[] = get(ENV, "JULIA_ASSETREGISTRY_BASEURL", get(ENV, "JUPYTERHUB_SERVICE_PREFIX", ""))
    atexit() do
        for (key, path) in AssetRegistry.registry
            AssetRegistry.deregister(path)
        end
    end
end

end # module


text:
import Dates: now

"""
    RealTimeClock(functor::TimeFunc)

A real time clock (system time) which use a `TimeFunc` functor

`real_time_clock` is default system clock. It's using UTC DateTime
"""
struct RealTimeClock <: AbstractClock
    func::TimeFunc
end

real_time_clock = RealTimeClock(UTCDateTimeFunc)

"""
    now(clock::RealTimeClock)

Return current instant from a (system) clock
"""
function now(clock::RealTimeClock)
    clock.func()
end

"""
    set(clock::RealTimeClock, value)

Setting clock to a system clock is not an allowed operation
"""
function set(clock::RealTimeClock, value)
    throw(ClockException("RealTimeClock can't be set"))
end


text:
# Literals
0
'a'
"test"

# Math
1+1
1+1/(7+5)

# Logic expressions
1 == 1
true > false

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
include("test_examples_3d_part1.jl")
include("test_examples_3d_part2.jl")


text:
using GitDataManager
using Test

@testset "GitDataManager.jl" begin
    # Write your tests here.
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the acoustic perturbation equations

equations = AcousticPerturbationEquations2D(v_mean_global=(0.0, 0.0), c_mean_global=0.0,
                                            rho_mean_global=0.0)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-20.6, 0.0) # minimum coordinates (min(x), min(y))
coordinates_max = (30.6, 51.2) # maximum coordinates (max(x), max(y))

"""
  initial_condition_monopole(x, t, equations::AcousticPerturbationEquations2D)

Initial condition for the monopole in a boundary layer setup, used in combination with
[`boundary_condition_monopole`](@ref).
"""
function initial_condition_monopole(x, t, equations::AcousticPerturbationEquations2D)
  m = 0.3 # Mach number

  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = 0.0

  v1_mean = x[2] > 1 ? m : m * (2*x[2] - 2*x[2]^2 + x[2]^4)
  v2_mean = 0.0
  c_mean = 1.0
  rho_mean = 1.0

  prim = SVector(v1_prime, v2_prime, p_prime, v1_mean, v2_mean, c_mean, rho_mean)

  return prim2cons(prim, equations)
end
initial_condition = initial_condition_monopole # does not use the global mean values given above

"""
  boundary_condition_monopole(u_inner, orientation, direction, x, t, surface_flux_function,
                              equations::AcousticPerturbationEquations2D)

Boundary condition for a monopole in a boundary layer at the -y boundary, i.e. `direction = 3`.
This will return an error for any other direction. This boundary condition is used in combination
with [`initial_condition_monopole`](@ref).
"""
function boundary_condition_monopole(u_inner, orientation, direction, x, t, surface_flux_function,
                                     equations::AcousticPerturbationEquations2D)
  if direction != 3
    error("expected direction = 3, got $direction instead")
  end

  # Wall at the boundary in -y direction with a monopole at -0.05 <= x <= 0.05. In the monopole area
  # we use a sinusoidal boundary state for the perturbed variables. For the rest of the -y boundary
  # we set the boundary state to the inner state and multiply the perturbed velocity in the
  # y-direction by -1.
  if -0.05 <= x[1] <= 0.05 # Monopole
    v1_prime = 0.0
    v2_prime = p_prime = sin(2 * pi * t)

    prim_boundary = SVector(v1_prime, v2_prime, p_prime, u_inner[4], u_inner[5], u_inner[6], u_inner[7])

    u_boundary = prim2cons(prim_boundary, equations)
  else # Wall
    u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4], u_inner[5], u_inner[6],
                         u_inner[7])
  end

  # Calculate boundary flux
  flux = surface_flux_function(u_boundary, u_inner, orientation, equations)

  return flux
end

"""
    boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::AcousticPerturbationEquations2D)

Boundary condition that uses a boundary state where the state variables are zero and the mean
variables are the same as in `u_inner`.
"""
function boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                                 equations::AcousticPerturbationEquations2D)
  value = zero(eltype(u_inner))
  u_boundary = SVector(value, value, value, cons2mean(u_inner, equations)...)

  # Calculate boundary flux
  if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

boundary_conditions = (x_neg=boundary_condition_zero,
                       x_pos=boundary_condition_zero,
                       y_neg=boundary_condition_monopole,
                       y_pos=boundary_condition_zero)

# Create a uniformly refined mesh with periodic boundaries
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=6,
                n_cells_max=100_000,
                periodicity=false)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 24.0
tspan =  (0.0, 24.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=0.8)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks)

# Print the timer summary
summary_callback()


text:
import DielectricUpdateTechique as DUTCH
using Test

nCells = [25,54,37];
p,q,r = nCells
resolution = [1f-3,1f-3,1f-3];
m = DUTCH.get_constants(300f6);
k_b = m.kb

AToE = DUTCH._electric_vector_potential_to_electric_field_operator(nCells,resolution,k_b);

xl = (p+2)*(q+1)*(r+1)
yl = (p+1)*(q+2)*(r+1)
zl = (p+1)*(q+1)*(r+2)

@test size(AToE) == (xl+yl+zl, xl+yl+zl)
@test real(AToE[2,2]) == -2/resolution[1]^2 + k_b^2
@test real(AToE[3,2]) == 1/resolution[1]*1/resolution[2]

Ig = DUTCH._create_Greens_functions_restriction_operators(nCells);
dx,dy,dz = 2 .^(ceil.(Int,log2.(nCells.+2).+1))

@test size(Ig) == (xl+yl+zl, dx*dy*dz*3)


text:
@testset "Array variables" begin
    lfp = LinearFractionalModel(with_optimizer(Clp.Optimizer))
    a = [-2, 0]
    x = @variable(lfp, [i=1:2], base_name="x", lower_bound=0)
    @constraint(lfp, x[2] <= 6)
    @constraint(lfp, -x[1] + x[2] <= 4)
    @constraint(lfp, 2x[1] + x[2] <= 14)
    numer = @expression(lfp, sum(a[i] * x[i] for i in 1:2) + 2)
    denom = @expression(lfp, x[1] + 3x[2] + 4)
    set_objective(lfp, MOI.MIN_SENSE, numer, denom)
    optimize!(lfp)
    @test termination_status(lfp) == MOI.OPTIMAL
    @test all(value.(x) .≈ [7.0, 0.0])
end

@testset "Array constraints" begin
    lfp = LinearFractionalModel(with_optimizer(Clp.Optimizer))
    x = @variable(lfp, [i=1:2], base_name="x")
    a = [2.0, 1.0]
    upbs = [4.0, 20.0]
    lbs = [-4.0, -20.0]
    @constraint(lfp, [i=1:2], x[i] <= upbs[i])
    @constraint(lfp, [i=1:2], x[i] >= lbs[i])
    numer = @expression(lfp,  sum(a[i] * x[i] for i in 1:2) + 2)
    denom = @expression(lfp,  -x[1] + 20*x[2] + 4)
    set_objective(lfp, MOI.MIN_SENSE, numer, denom)
    optimize!(lfp)
    @test termination_status(lfp) === MOI.OPTIMAL
    @test value.(x) ≈ [-4.0, -0.2]
end

@testset "Match LP" begin
    lfp = LinearFractionalModel(with_optimizer(Clp.Optimizer))
    a = [4, 2]
    upbs = [4.0, 20.0]
    lbs = [-1.0, -10.0]
    x = @variable(lfp, [i=1:2], base_name="x", lower_bound=lbs[i], upper_bound=upbs[i])
    @constraint(lfp, x[1] + x[2] <= 5.0)
    @constraint(lfp, x[1] - 2*x[2] >= 10.0)
    numer = @expression(lfp,  dot(a, x))
    denom = @expression(lfp,  sum(x))
    set_objective(lfp, MOI.MIN_SENSE, numer, denom)
    optimize!(lfp)
    m = Model(with_optimizer(Clp.Optimizer))
    xm = @variable(m, [i=1:2], base_name="x")
    t = @variable(m, lower_bound=0.0)
    a = [4, 2]
    upbs = [4, 20]
    lbs = [-1,-10]
    @constraint(m, xm[1] + xm[2] <= 5.0t)
    @constraint(m, xm[1] - 2*xm[2] >= 10.0t)
    @constraint(m, [i=1:2], xm[i] <= upbs[i] * t)
    @constraint(m, [i=1:2], xm[i] >= lbs[i] * t)
    @constraint(m, sum(xm) == 1)
    @objective(m,  Min, sum(a[i] * xm[i] for i in 1:2))
    optimize!(m)
    @test value.(xm)/value(t) ≈ value.(x)
    @test objective_value(lfp) ≈ objective_value(m)
end

@testset "Constant denominator" begin
    lfp = LinearFractionalModel(with_optimizer(Clp.Optimizer))
    a = [4, 2]
    upbs = [4, 20]
    lbs = [-1, -10]
    x = @variable(lfp, [i=1:2], base_name="x", lower_bound=lbs[i], upper_bound=upbs[i])
    @constraint(lfp, x[1] + x[2] <= 5.0)
    @constraint(lfp, x[1] - 2*x[2] >= 10.0)
    @constraint(lfp, [i=1:2], x[i] <= upbs[i])
    @constraint(lfp, [i=1:2], x[i] >= lbs[i])
    numer = @expression(lfp, sum(a[i] * x[i] for i in 1:2))
    denom = @expression(lfp, 5.0)
    set_objective(lfp, MOI.MIN_SENSE, numer, denom)
    optimize!(lfp)

    m = Model(with_optimizer(Clp.Optimizer))
    a = [4, 2]
    upbs = [4, 20]
    lbs = [-1,-10]
    xm = @variable(m, [i=1:2], base_name="x", lower_bound=lbs[i], upper_bound=upbs[i])
    @constraint(m, xm[1] + xm[2] <= 5.0)
    @constraint(m, xm[1] - 2*xm[2] >= 10.0)
    @objective(m,  Min, sum(a[i] * xm[i] for i in 1:2)/5.0)
    optimize!(m)
    @test value.(xm) ≈ value.(x)
    @test objective_value(lfp) ≈ objective_value(m)
end


text:
#######################################################################
# randomkets
#######################################################################
"""
    `randomkets(h, A, n)`

Creates `n` random vectors with length `D = dim(h)` and stacks them into a `D×n` matrix.
The resulting ket amplitudes are normalized by `1/√n`.

This type of ket structure is useful e.g. in stochastic trace evaluation of KPM methods,
where the amplitude, `a = cis(2pi*rand())` is chosen as a random function s.t. `⟨a⟩ = 0`,
`⟨aa⟩ = 0` and `⟨a'a⟩ = 1`. Then,
the normalized trace of an operator `A` can be estimated with `Tr[A]/N₀ ≈ ∑⟨ket|A|ket⟩`,
where the sum is taken over the `n` random kets `|ket⟩` of norm `1/√n` produced by
`randomkets`, and `N₀` is the total number of orbitals in the full unit cell.

By default when applied to a multiorbital system with a maximum of `N` orbitals per site,
the generated kets have independent, complex, normally-distributed random components.

    `randomkets(h, n) = randomkets(h, I, n)`

useful for dosKPM
"""
randomkets(h, n) = randomkets(h, I, n)

function randomkets(h, A, n)
    dim = size(h,1)
    rkets = zeros(promote_type(eltype(h), eltype(A)), dim, n)
    [rkets[i,j] = cis(2pi*rand())/√n for i in 1:dim for j in 1:n]
    return rkets
end

#######################################################################
# Kernel Polynomial Method : momenta
#######################################################################
using Base.Threads

struct MomentaKPM{T, B<:Tuple}
    mulist::Vector{T}
    bandbracket::B
end

struct KPMBuilder{AM,HM<:AbstractMatrix,T,K<:AbstractMatrix,B}
    hmat::HM
    Amat::AM
    ket0::K
    ket1::K
    ket2::K
    bandbracket::Tuple{B,B}
    order::Int
    mulist::Vector{T}
end

function KPMBuilder(h, A, ket_or_randkets, order, bandrange) where {L}
    eh = eltype(eltype(h))
    eA = eltype(eltype(A))
    mulist = zeros(promote_type(eh, eA), order + 1)
    bandbracket = bandbracketKPM(h, bandrange)
    zeroket = ketseed_KPM(ket_or_randkets)
    builder = KPMBuilder(h, A, zeroket, similar(zeroket), similar(zeroket), bandbracket,
        order, mulist)
    return builder
end

ketseed_KPM(k::AbstractVector) =
    size(k,1) > size(k,2) ? reshape(k, :, 1) : reshape(Matrix(k'), :, 1)

ketseed_KPM(k::AbstractMatrix) = reshape(k[:,1],:,1)

"""
    `momentaKPM(h, A = I; kets = randomkets(h, A, n), order = 10, bandrange = missing)`

Compute the Kernel Polynomial Method (KPM) momenta `μₙ = ⟨k|Tₙ(h) A|k⟩`, where `|k⟩` may
either refer to a single ket (kets::AbstractVector), or a set of `kets::AbstractMatrix`
stored as columns, or to a collection of `n` random kets generated by `randomkets(h, A, n)`.
`A` is an observable (`AbstractMatrix`) and `Tₙ(h)` is the order-`n` Chebyshev polynomial of
the Hamiltonian `h`.

For the randomkets case, `μₙ` is summed over all kets, `μₙ= ∑_k ⟨k|Tₙ(h) A|k⟩`.
If `kets` are a multi-column ket, a sum over columns will similarly be performed.

`kets` generated with `randomkets()` are useful to estimate momenta of normalized traces
using the stochastic trace approach, whereby `μ_n = Tr[A T_n(h)]/N₀ ≈ ∑ₖ⟨k|A T_n(h)|k⟩`.
Here the `|k⟩`s are `n` random kets of norm `1/√n` and `N₀` is the total number of orbitals
per unit cell of `h` (see `randomkets`).

The order of the Chebyshev expansion is `order`. The `bandbrange = (ϵmin, ϵmax)` should
completely encompass the full bandwidth of `hamiltonian`. If `missing` it is computed
automatically using `ArnoldiMethod` (must be loaded `using ArnoldiMethod`).

# Examples

```
julia> momentaKPM(h, I, kets = randomkets(h, I, 5), order = 20)

julia> momentaKPM(h, bandrange = (-6,6)).mulist |> length

```
"""
function momentaKPM(h, A = I; kets = randomkets(h, A, 1), order = 10, bandrange = missing)
    builder = KPMBuilder(h, A, kets, order, bandrange)
    momentaKPM!(builder, kets)
    jackson!(builder.mulist)
    return MomentaKPM(builder.mulist, builder.bandbracket)
end

function momentaKPM!(b::KPMBuilder, kets::AbstractMatrix)
    aux_ket = b.ket0 # pointer to b.ket0
    for i in 1:size(kets, 2) #number of vectors
        aux_ket .= kets[:, i] #updating builder content
        pmeter = Progress(b.order, "Computing moments: ")
        addmomentaKPM!(b, pmeter)
    end
    return nothing
end

function momentaKPM!(b::KPMBuilder, kets::AbstractVector)
    momentaKPM!(b, reshape(kets,:,1))
end

# This iterates bras <psi_n| = <psi_0|AT_n(h) instead of kets (faster CSC multiplication)
# In practice we iterate their conjugate |psi_n> = T_n(h') A'|psi_0>, and do the projection
# onto the start ket, |psi_0>
function addmomentaKPM!(b::KPMBuilder{<:AbstractMatrix,<:AbstractSparseMatrix}, pmeter)
    mulist = b.mulist
    kmat, kmat0, kmat1 = parent.((b.ket0, b.ket1, b.ket2))
    h, A, bandbracket = b.hmat, b.Amat, b.bandbracket
    order = length(mulist) - 1
    mul!(kmat0, A', kmat)
    mulscaled!(kmat1, h', kmat0, bandbracket)
    mulist[1] += proj(kmat0, kmat)
    mulist[2] += proj(kmat1, kmat)
    for n in 3:(order+1)
        ProgressMeter.next!(pmeter; showvalues = ())
        iterateKPM!(kmat0, h', kmat1, bandbracket)
        mulist[n] += proj(kmat0, kmat)
        kmat0, kmat1 = kmat1, kmat0
    end
    return mulist
end

function addmomentaKPM!(b::KPMBuilder{<:UniformScaling, <:AbstractSparseMatrix,T}, pmeter) where {T}
    mulist = b.mulist
    kmat, kmat0, kmat1 = parent.((b.ket0, b.ket1, b.ket2))
    h, A, bandbracket = b.hmat, b.Amat, b.bandbracket
    order = length(mulist) - 1
    kmat0 .= kmat
    mulscaled!(kmat1, h', kmat0, bandbracket)
    mulist[1] += μ0 = proj(kmat0, kmat0)
    mulist[2] += μ1 = proj(kmat1, kmat0)
    # This is not used in the currently activated codepath (BLAS mul!), but is needed in the
    # commented out @threads codepath
    thread_buffers = (zeros(T, Threads.nthreads()), zeros(T, Threads.nthreads()))
    for n in 3:2:(order+1)
        ProgressMeter.next!(pmeter; showvalues = ())
        ProgressMeter.next!(pmeter; showvalues = ()) # twice because of 2-step
        proj11, proj10 = iterateKPM!(kmat0, h', kmat1, bandbracket, thread_buffers)
        mulist[n] += 2 * proj11 - μ0
        n + 1 > order + 1 && break
        mulist[n + 1] += 2 * proj10 - μ1
        kmat0, kmat1 = kmat1, kmat0
    end
    A.λ ≈ 1 || (mulist .*= A.λ)
    return mulist
end

function mulscaled!(y, h´, x, (center, halfwidth))
    mul!(y, h´, x)
    invhalfwidth = 1/halfwidth
    @. y = (y - center * x) * invhalfwidth
    return y
end

function iterateKPM!(kmat0, h´, kmat1, (center, halfwidth), buff = ())
    α = 2/halfwidth
    β = 2center/halfwidth
    mul!(kmat0, h´, kmat1, α, -1)
    @. kmat0 = kmat0 - β * kmat1
    return proj_or_nothing(buff, kmat0, kmat1)
end

proj_or_nothing(::Tuple{}, kmat0, kmat1) = nothing
proj_or_nothing(buff, kmat0, kmat1) = (proj(kmat1, kmat1), proj(kmat1, kmat0))

# This is equivalent to tr(kmat1'*kmat2) for matrices
proj(kmat1, kmat2) = dot(kmat1, kmat2)

# function iterateKPM!(ket0, h´, ket1, (center, halfwidth), thread_buffers = ())
#     h = parent(h´)
#     nz = nonzeros(h)
#     rv = rowvals(h)
#     α = -2 * center / halfwidth
#     β = 2 / halfwidth
#     reset_buffers!(thread_buffers)
#     @threads for row in 1:size(ket0, 1)
#         ptrs = nzrange(h, row)
#         @inbounds for col in 1:size(ket0, 2)
#             k1 = ket1[row, col]
#             tmp = α * k1 - ket0[row, col]
#             for ptr in ptrs
#                 tmp += β * adjoint(nz[ptr]) * ket1[rv[ptr], col]
#             end
#             # |k0⟩ → (⟨k1|2h - ⟨k0|)' = 2h'|k1⟩ - |k0⟩
#             ket0[row, col] = tmp
#             update_buffers!(thread_buffers, k1, tmp)
#         end
#     end
#     return sum_buffers(thread_buffers)
# end

# reset_buffers!(::Tuple{}) = nothing
# function reset_buffers!((q, q´))
#     fill!(q, zero(eltype(q)))
#     fill!(q´, zero(eltype(q´)))
#     return nothing
# end

# @inline update_buffers!(::Tuple{}, k1, tmp) = nothing
# @inline function update_buffers!((q, q´), k1, tmp)
#     q[threadid()]  += dot(k1, k1)
#     q´[threadid()] += dot(tmp, k1)
#     return nothing
# end

# @inline sum_buffers(::Tuple{}) = nothing
# @inline sum_buffers((q, q´)) = (sum(q), sum(q´))

function jackson!(μ::AbstractVector)
    order = length(μ) - 1
    @inbounds for n in eachindex(μ)
        μ[n] *= ((order - n + 1) * cos(π * n / (order + 1)) +
                sin(π * n / (order + 1)) * cot(π / (order + 1))) / (order + 1)
    end
    return μ
end

function bandbracketKPM(h, ::Missing)
    bandbracketKPM(h, bandrangeKPM(h))
end

bandbracketKPM(h, (ϵmin, ϵmax)::Tuple{T,T}, pad = float(T)(0.01)) where {T} =
    ((ϵmax + ϵmin) / 2, (ϵmax - ϵmin) / (2 - pad))

function bandrangeKPM(h::AbstractMatrix{T}; quiet = false) where {T}
    if quiet == true
        nothing
    else @warn "Computing spectrum bounds... Consider using the `bandrange`
            option for faster performance."
    end
    R = real(T)
    decompl, _ = Main.ArnoldiMethod.partialschur(h, nev=1, tol=1e-4,
        which = Main.ArnoldiMethod.LR());
    decomps, _ = Main.ArnoldiMethod.partialschur(h, nev=1, tol=1e-4,
        which = Main.ArnoldiMethod.SR());
    ϵmax = R(real(decompl.eigenvalues[1]))
    ϵmin = R(real(decomps.eigenvalues[1]))
    if quiet == true nothing
    else @warn("Computed bandrange = ($ϵmin, $ϵmax)") end
    return (ϵmin, ϵmax)
end

#######################################################################
# Kernel Polynomial Method : observables
#######################################################################
"""
    dosKPM(h::AbstractMatrix; resolution = 2, kets = randomkets(h, 1), kw...)

Compute, using the Kernel Polynomial Method (KPM), the local density of states `ρₖ(ϵ) =
⟨k|δ(ϵ-h)|k⟩` for a single ket `|k⟩`, a collection of `{|k_i⟩}`, or a set of randomly-
generated kets for a stochastic trace calculation(see `momentaKPM`).

The result is a tuple of energy points `ϵᵢ::Vector` spanning the band range,
and real `ρₖ(ϵᵢ)::Vector` values (any residual imaginary part in `ρₖ` is dropped). The
number of energy points `ϵᵢ` is `order * resolution`, rounded to the closest integer.

If `kets` contains more than a single ket, the sum `∑ₖρₖ(ε)` over all models will be
computed. In the case of the default `kets = randomkets(h, n)`, this results in an
estimate of the total density of states per orbital, computed through an stochastic trace,
`ρ(ϵ) = ∑ₖ⟨k|δ(ϵ-h)|k⟩/n ≈ Tr[δ(ϵ-h)]/N₀`, where `N₀` is the total number of orbitals in the
unit cell.

`dosKPM` is a particular case of `densityKPM` for an operator `A = I` and with any residual
imaginary parts dropped

    `dosKPM(μ::MomentaKPM; resolution = 2)`

Same as above with KPM momenta `μ` as input.

# See also
    `momentaKPM`, `densityKPM`, `averageKPM`
"""
dosKPM(h; resolution = 2, kets = randomkets(h, 1), kw...) =
    dosKPM(momentaKPM(h, I; kets = kets, kw...), resolution = resolution)

dosKPM(μ::MomentaKPM; resolution = 2) = real.(densityKPM(μ; resolution = resolution))

"""
    densityKPM(h::AbstractMatrix, A; resolution = 2, kets = randomkets(1), kw...)

Compute, using the Kernel Polynomial Method (KPM), the spectral density of `A`, `ρᴬₖ(ϵ) =
⟨k|A δ(ϵ-h)|k⟩` for a single ket `|k⟩`, a collection of `{|k_i⟩}`, or a set of randomly-
generated kets for a stochastic trace calculation(see `momentaKPM`).

The result is a tuple of energy points `ϵᵢ::Vector` spanning the band  range, and real
`ρᴬₖ(ϵᵢ)::Vector` values. The number of energy points `ϵᵢ` is `order *resolution`,
rounded to the closest integer.

If `kets` contains more than a single ket, the sum `∑ₖρₖ(ε)` over all models will be
computed. In the case of the default `kets = randomkets(h, I, n)`, this results in an
estimate of the total density of states per orbital, computed through an stochastic trace,
`ρ(ϵ) = ∑ₖ⟨k|δ(ϵ-h)|k⟩/n ≈ Tr[δ(ϵ-h)]/N₀`, where `N₀` is the total number of orbitals in the
unit cell.

    `densityKPM(μ::MomentaKPM; resolution = 2)`

Same as above with KPM momenta `μ` as input.

# See also
    `dosKPM`, `momentaKPM`, `averageKPM`
"""
densityKPM(h, A; resolution = 2, kw...) =
    densityKPM(momentaKPM(h, A; kw...); resolution = resolution)

function densityKPM(momenta::MomentaKPM{T}; resolution = 2) where {T}
    (center, halfwidth) = momenta.bandbracket
    numpoints = round(Int, length(momenta.mulist) * resolution)
    ρlist = zeros(T, numpoints)
    copyto!(ρlist, momenta.mulist)
    Main.FFTW.r2r!(ρlist, Main.FFTW.REDFT01, 1)  # DCT-III in FFTW
    xk = [cos(π * (k + 0.5) / numpoints) for k in 0:numpoints - 1]
    @. ρlist = center + halfwidth * ρlist / (π * sqrt(1.0 - xk^2))
    @. xk = center + halfwidth * xk
    return xk, ρlist
end

"""
    `averageKPM(h::AbstractMatrix, A; kBT = 0, Ef = 0, kw...)`

Compute, using the Kernel Polynomial Method (KPM), the thermal expectation value `⟨A⟩ = Σ_k
f(E_k) ⟨k|A|k⟩ = ∫dE f(E) Tr [A δ(E-H)]/N₀ = Tr [A f(H)]/N₀` for a given hermitian operator
`A` and a zero-dimensional hamiltonian `h` with a total of `N₀` orbitals (see `momentaKPM`
and its options `kw` for further details). `f(E)` is the Fermi-Dirac distribution function,
`|k⟩` are `h` eigenstates with energy `E_k`, kBT` is the temperature in energy units and
`Ef` the Fermi energy.

    `averageKPM(μ::MomentaKPM, A; kBT = 0, Ef = 0)``

Same as above with the KPM momenta as input (see `momentaKPM`).

# See also
    `dosKPM`, `momentaKPM`, `averageKPM`
"""
averageKPM(h, A; kBT = 0, Ef = 0, kw...) = averageKPM(momentaKPM(h, A; kw...); kBT = kBT,
    Ef = Ef)

function averageKPM(momenta::MomentaKPM; kBT = 0.0, Ef = 0.0)
    (center, halfwidth) = momenta.bandbracket
    order = length(momenta.mulist) - 1
    # if !iszero(kBT)
    #     @warn "Finite temperature requires numerical evaluation of the integrals"
    #     checkloaded(:QuadGK)
    # end
    average = sum(n -> momenta.mulist[n + 1] * fermicheby(n, Ef, kBT, center, halfwidth),
         0:order)
    return average
end

# Pending issue: Unexpected behaviour with center != 0.
function fermicheby(n, Ef, kBT, center, halfwidth)
    kBT´ = kBT / halfwidth;
    Ef´ = (Ef - center) / halfwidth;
    T = typeof(Ef´)
    if kBT´ == 0
        int = n == 0 ? 0.5+asin(Ef´)/π : -2.0*sin(n*acos(Ef´))/(n*π)
    else
        throw(error("Finite temperature not yet implemented"))
        # η = 1e-10
        # int = Main.QuadGK.quadgk(E´ -> _intfermi(n, E´, Ef´, kBT´), -1.0+η, 1.0-η,
        #   atol= 1e-10, rtol=1e-10)[1]
    end
    return T(int)
end

_intfermi(n, E´, Ef´, kBT´) = fermifun(E´, Ef´, kBT´) * 2/(π*(1-E´^2)^(1/2)) *
    chebypol(n, E´) / (1+(n == 0 ? 1 : 0))

fermifun(E´, Ef´, kBT´) = kBT´ == 0 ? (E´<Ef´ ? 1.0 : 0.0) : (1/(1+exp((E´-Ef´)/(kBT´))))

function chebypol(m::Int, x::T) where {T<:Number}
    cheby0 = one(T)
    cheby1 = x
    if m == 0
        chebym = cheby0
    elseif m == 1
        chebym = cheby1
    else
        for i in 2:m
            chebym = 2x * cheby1 - cheby0
            cheby0, cheby1 = cheby1, chebym
        end
    end
    return chebym
end

max_stars_repo_name: fernandopenaranda/ChebyshevExpansions

max_issues_repo_name: fernandopenaranda/ChebyshevExpansions

max_forks_repo_name: fernandopenaranda/ChebyshevExpansions


text:
using ErrorStateKalmanFilter
using Test
using LinearAlgebra
using Quaternions

@testset "types" begin
    x = collect(range(0, length=19, step=0.5))
    @test vec(NominalState(x)) == x
    x = collect(range(0, length=18, step=0.1))
    @test vec(ErrorState(x)) == x
end

@testset "ESKF" begin
    @test skew_matrix([1, 2, 3]) == [0 -3 2; 3 0 -1; -2 1 0]

    dt = 1e-3
    times = 0:dt:5

    nominal = NominalState(
        [0.0, 0, 0],
        [0.0, 0, 0],
        quat(1.0),
        zeros(3),
        zeros(3),
        [0.0, 0.0, -1.0]
    )

    error = ErrorState(zeros(18))

    error_cov = Matrix{Float64}(I, 18, 18)

    measurement = Measurement([0.0, 0, 0], [0.0, 0, 0])

    noise = Noise(
        2e-7  * (1.0/512),
        1e-10 * (1.0/512),
        1e-9,
        1e-9
    )

    new_nominal, new_error, new_error_cov = predict(nominal, error, error_cov, measurement, noise, dt)

    # error state is unchanged
    @test new_error == error
    # error covariance increases
    @test all(abs.(new_error_cov) .>= abs.(error_cov))
end


text:
using AbstractOperators
using ProximalOperators

A1 = [ 1.0  -2.0   3.0  -4.0  5.0;
       2.0  -1.0   0.0  -1.0  3.0]
A2 = [-1.0   0.0   4.0  -3.0  2.0;
      -1.0  -1.0  -1.0   1.0  3.0]
A = [A1; A2]
opA = VCAT(MatrixOp(A1), MatrixOp(A2))

b1 = [1.0, 2.0]
b2 = [3.0, 4.0]
b = [b1; b2]

f = SeparableSum(Translate(SqrNormL2(), -b1), Translate(SqrNormL2(), -b2))
lam = 0.1*norm(A'*b, Inf)
g = NormL1(lam)

x_star = [-3.877278911564627e-01, 0, 0, 2.174149659863943e-02, 6.168435374149660e-01]

# Nonfast/Nonadaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x .- x_star, Inf) <= 1e-4
@test it < 150
println(sol)

# Nonfast/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true)
@test norm(x .- x_star, Inf) <= 1e-4
@test it < 300
println(sol)

# Fast/Nonadaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2, fast=true)
@test norm(x .- x_star, Inf) <= 1e-4
@test it < 100
println(sol)

# Fast/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true, fast=true)
@test norm(x .- x_star, Inf) <= 1e-4
@test it < 200
println(sol)

# ZeroFPR/Nonadaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 15
println(sol)

# ZeroFPR/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=opA, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20

# PANOC/Nonadaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)

# PANOC/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=opA, g=g, adaptive=true)
@test norm(x - x_star, Inf) <= 1e-4
@test it < 20
println(sol)


text:
include("base.jl")
include("../providers/data.jl")

function analyze(attr::DataAttribute, ctx::PlotContext)
    # TODO: We really should pass the context here to get dates and titles specific to this context, but analyze can't handle contexts that good so we can't do that. But we definetly need to implement that!
    #=local srcCtx
    if "source context" in keys(ctx.styles)
        if typeof(ctx.styles["source context"]) <: Array
            srcCtx = Dict("$x"=>x for x in ctx.styles["source context"])
        else
            @assert typeof(ctx.styles["source context"]) <: DataAttributeContext
            srcCtx = Dict("$(ctx.styles["source context"])"=>ctx.styles["source context"])
        end
    else
        srcCtx = fetchContext!(ctx.store, attr, ctx.minDate, ctx.maxDate)
    end=#

    local srcCtx = fetchContext!(ctx.store, attr, ctx.minDate, ctx.maxDate)

    @assert length(srcCtx) > 0 && typeof(srcCtx) <: Dict{String, DataAttributeContext}
    for (k,lctx) in srcCtx
        local dates, values, title = fetchData!(ctx.store, attr, lctx, ctx.minDate, ctx.maxDate)

        # TODO: mini und max
        return dates, srcCtx, minimum(values), maximum(values), [title]
    end
end

import Base.+
+(arr::Array{<:Number,1}, b::Number) = begin
    arr = copy(arr)
    if b == 0
        return arr
    end
    for i in 1:length(arr)
        arr[i] += b
    end
    return arr
end

function renderWidget!(attr::DataAttribute, ctx::PlotContext)

    local dates, values, _ = fetchData!(ctx.store, attr, ctx.styles["source context"], ctx.minDate, ctx.maxDate)

    local label = attr.property
    if "stacked" in keys(ctx.styles)
        label *= " (stacked)"
    end
    if "append to label" in keys(ctx.styles)
        label *= " " * ctx.styles["append to label"]
    end

    local color = ctx.styles["color"]
    local fcolor = ctx.styles["fillcolor"]
    if ctx.styles["color"] === nothing
        color = nextColor!(ctx.palette)
    end
    if fcolor === nothing
        fcolor = color
    end

    # For example in stacked plots we need some offset
    values += ctx.styles["offset"]

    local fillrange = ctx.styles["fillrange"]
    # Filling a log plot is tricky! Don't go below the minimum!
    if "log" in keys(ctx.styles)
        local l = ctx.styles["log"]
        if typeof(fillrange) <: Array
            fillrange = map(x -> max(x, l), fillrange)
        elseif typeof(fillrange) <: Number
            fillrange = max(fillrange, l)
        end
    end

    local ylims = nothing
    if "ylims" in keys(ctx.styles)
        ylims = ctx.styles["ylims"]
    end

    # add to graph
    plot!(
        ctx.plot,
        dates,
        values,
        label = label,
        color = color,
        ylims = ylims,
        fillrange = fillrange,
        fillcolor = fcolor,
        fillalpha = ctx.styles["fillalpha"],
    )

    return dates, values
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
using PaddleBattle

PaddleBattle.julia_main([""])

max_stars_repo_name: NHDaly/pongClone

max_issues_repo_name: NHDaly/pongClone

max_forks_repo_name: NHDaly/pongClone


text:
using Quandl, DataFrames, Redis, Winston

qapi_key = "aaaaaaaaaaaaaaaaaaaa"	# Your Quandl key
set_auth_token(qapi_key)

qdf1 = quandl("WIKI/AAPL", format="DataFrame", order="asc");
aapl = convert(Array, qdf1[:Close]);
scf1 = 1.0/aapl[1];  aapl = scf1 .* aapl;

qdf2 = quandl("WIKI/MSFT", format="DataFrame", order="asc");
msft = convert(Array, qdf2[:Close]);
scf2 = 1.0/msft[1];  msft = scf2 .* msft;

n = [1:length(aapl)];
conn = RedisConnection()
for i = 1:n
  rpush(conn,'APPL~Close',aapl[i])
  rpush(conn,'MSFT~Close',msft[i])
end

t = [1:n];		# Just plot the data not the dates

aapl-data = float32(lrange(conn,'AAPL~Close',0,-1);
msft-data = float32(lrange(conn,'MSFT~Close',0,-1);
plot(t, aapl-data,"b--", t, msft-data, "r.");


text:
set_threadid(task::Task, i) = set_raw_tid(task, i - 1)
unset_threadid(task::Task) = set_raw_tid(task, -1)

function set_raw_tid(task::Task, raw_tid)
    raw_tid = convert(Int16, raw_tid)
    ptr = Ptr{Int16}(pointer_from_objref(task)) + TASK_TID_OFFSET
    GC.@preserve task begin
        return unsafe_store!(ptr, raw_tid)
    end
end

function check_set_raw_tid()
    task = Task(nothing)
    for raw_tid in Int16(0):Int16(20)
        set_raw_tid(task, raw_tid)
        actual = Threads.threadid(task) - 1
        if raw_tid != actual
            error("expected tid ($raw_tid) does not match with the actual tid ($actual)")
        end
    end
end


text:
#=

  From Yi Wu, Lei Li, Stuart Russell
  "BFiT: From Possible-World Semantics to Random-Evaluation Semantics in Open Universe"
  Page 3 (BLOG model)

  Distributions of variable numPeople (num:0)
  4.00000 =>    1763  (0.176300)
  5.00000 =>    1757  (0.175700)
  6.00000 =>    1519  (0.151900)
  3.00000 =>    1406  (0.140600)
  7.00000 =>    1067  (0.106700)
  2.00000 =>     845  (0.084500)
  8.00000 =>     632  (0.063200)
  9.00000 =>     372  (0.037200)
  1.00000 =>     326  (0.032600)
  10.00000 =>     180  (0.018000)
  11.00000 =>      83  (0.008300)
  12.00000 =>      28  (0.002800)
  13.00000 =>      12  (0.001200)
  14.00000 =>       5  (0.000500)
  15.00000 =>       4  (0.000400)
  16.00000 =>       1  (0.000100)

  Distributions of variable sumHonest (num:0)
  4.00000 =>    1893  (0.189300)
  5.00000 =>    1743  (0.174300)
  3.00000 =>    1727  (0.172700)
  6.00000 =>    1314  (0.131400)
  2.00000 =>    1124  (0.112400)
  7.00000 =>     840  (0.084000)
  1.00000 =>     472  (0.047200)
  8.00000 =>     461  (0.046100)
  9.00000 =>     221  (0.022100)
  10.00000 =>     106  (0.010600)
  11.00000 =>      40  (0.004000)
  0.00000 =>      36  (0.003600)
  12.00000 =>      15  (0.001500)
  14.00000 =>       5  (0.000500)
  13.00000 =>       3  (0.000300)

  Distributions of variable honest1 (num:0)
  1.00000 =>    9061  (0.906100)
  0.00000 =>     939  (0.093900)

  Distributions of variable login1 (num:0)
  1.00000 =>    9208  (0.920800)
  0.00000 =>     750  (0.075000)
  2.00000 =>      30  (0.003000)
  3.00000 =>       9  (0.000900)
  4.00000 =>       3  (0.000300)

  Note: Currently this model got problems when using SMC() etc.
  For a discussion about this, see https://discourse.julialang.org/t/turing-jl-sometimes-this-error-is-thrown-cannot-convert-an-object-of-type-geometric-float64-to-an-object-of-type-dirac-float64/65202/3

  Since this model don't have any observable data, the use of an inference sampler might
  give problems (and it does now). However, using Prior() seems to work without any problem.

   Cf ~/webppl/person_login.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function person_login()

    # Number of people (exclude 0 people)
    numPeople ~ truncated(Poisson(5),0,Inf)

    # Is this person honest?
    honest ~ filldist(Bernoulli(0.9),numPeople)
    sumHonest ~ Dirac(sum(honest))

    # An honest person have just one login
    login = tzeros(numPeople)
    for p in 1:numPeople
        if honest[p]
            login[p] ~ Dirac(1.0)
        else
            login[p] ~ Geometric(0.8)
        end
    end

    login1 ~ Dirac(login[1])
    honest1 ~ Dirac(honest[1])

end

model = person_login()

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

# display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:numPeople)
show_var_dist_pct(chns,:sumHonest)
show_var_dist_pct(chns,:honest1)
show_var_dist_pct(chns,:login1)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@testset "calc_position func" begin
    time_vec = [0.0, 2.0]
    u_vec = [1.0, 1.0]
    v_vec = [0.0, 0.0]
    r_vec = [0.0, 0.0]
    x, y, ψ = calc_position(time_vec, u_vec, v_vec, r_vec, x0 = 1.0, y0 = 1.0, ψ0 = 0.0)
    @test x[2] ≈ 3.0
    @test y[2] ≈ 1.0
    @test ψ[2] ≈ 0.0
end

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text:
using GLMakie
using GraphMakie
using LightGraphs
using GLMakie.Colors

using GraphMakie: plot_controlpoints!, SquareGrid
g = complete_graph(3)
tangents = Dict(1 => ((1,1),(0,-1)),
                2 => ((0,1),(0,-1)),
                3 => ((0,-1),(1,0)))
tfactor = [0.5, 0.75, (0.5, 0.25)]
v_pos = [[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]]
f, ax, p = graphplot(g; node_pos = v_pos, tangents = tangents, tfactor = tfactor,
                     edge_width = [2.0 for _ = 1:ne(g)],
                     edge_color=[colorant"red", colorant"green", colorant"blue"],
                     node_size = [10 for _ = 1:nv(g)],
                     node_color = [colorant"black" for _ = 1:nv(g)],
                     elabels="Edge ".*repr.(1:ne(g)), elabels_distance=10)
hidedecorations!(ax); hidespines!(ax);
ax.aspect = DataAspect()
deregister_interaction!(ax, :rectanglezoom)
plot_controlpoints!(ax, p) # show control points for demonstration

function node_hover_action(state, idx, event, axis)
    p.node_size[][idx] = state ? 20 : 10
    p.node_size[] = p.node_size[] # trigger observable
end
nhover = NodeHoverHandler(node_hover_action)
register_interaction!(ax, :nhover, nhover)

function edge_hover_action(state, idx, event, axis)
    p.edge_width[][idx]= state ? 5.0 : 2.0
    p.edge_width[] = p.edge_width[] # trigger observable
end
ehover = EdgeHoverHandler(edge_hover_action)
register_interaction!(ax, :ehover, ehover)

function node_click_action(idx, args...)
    println(idx)
    p.node_color[][idx] = rand(RGB)
    p.node_color[] = p.node_color[]
end
nclick = NodeClickHandler(node_click_action)
register_interaction!(ax, :nclick, nclick)

function edge_click_action(idx, args...)
    println(idx)
    p.edge_color[][idx] = rand(RGB)
    p.edge_color[] = p.edge_color[]
end
eclick = EdgeClickHandler(edge_click_action)
register_interaction!(ax, :eclick, eclick)

function node_drag_action(state, idx, event, axis)
    p[:node_pos][][idx] = event.data
    p[:node_pos][] = p[:node_pos][]
end
ndrag = NodeDragHandler(node_drag_action)
register_interaction!(ax, :ndrag, ndrag)


text:
function check(fields, required_fields)
    return all(f in fields for f in required_fields)
end

function run(s)
    # Your code here
    passports = split(s, "\n\n")
    counter = 0
    required_fields = Set(
        ("byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid")
    )
    for passport_str in passports
        fields = split(replace(passport_str, '\n' => ' '), ' ')
        fields = map(x -> split(x, ":")[1], fields)
        if check(fields, required_fields)
            counter += 1
        end
    end
    return counter
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
module Runtests

using Test

@testset "Test ModiaBase" begin
    include("TestSymbolic.jl")
    include("TestBLTandPantelides.jl")
    include("TestDifferentiate.jl")
    include("TestTearing.jl")
    include("TestLinearIntegerEquations.jl")
end

end

max_stars_repo_name: ModiaSim/ModiaBase

max_issues_repo_name: ModiaSim/ModiaBase

max_forks_repo_name: ModiaSim/ModiaBase


text:
VERSION >= v"0.4.0-dev+6641" && __precompile__()

module Codecs

using Compat

export encode, decode, Base64, Zlib, BCD

abstract Codec

function encode{T <: Codec}(codec::Type{T}, s::AbstractString)
    encode(codec, convert(Vector{UInt8}, s))
end

function decode{T <: Codec}(codec::Type{T}, s::AbstractString)
    decode(codec, convert(Vector{UInt8}, s))
end

# RFC3548/RFC4648 base64 codec

abstract Base64 <: Codec

const b64enc_tbl = UInt8[
    'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
    'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P',
    'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X',
    'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f',
    'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n',
    'o', 'p', 'q', 'r', 's', 't', 'u', 'v',
    'w', 'x', 'y', 'z', '0', '1', '2', '3',
    '4', '5', '6', '7', '8', '9', '+', '/']

const base64_pad = @compat UInt8('=')
const sentinel = typemax(UInt8)

const b64dec_tbl = fill(sentinel, 256)
for (val, ch) in enumerate(b64enc_tbl)
    b64dec_tbl[ch] = val - 1
end

function base64dec(c::UInt8)
    v = b64dec_tbl[c]
    if v == sentinel
        error("Invalid base64 symbol: $(char(c))")
    end
    v
end

function encode(::Type{Base64}, input::Vector{UInt8})
    n = length(input)
    if n == 0
        return Array(UInt8, 0)
    end

    m = @compat Int(4 * ceil(n / 3))
    output = Array(UInt8, m)

    k = 1
    for ii = 1:3:length(input)-2
        u, v, w = input[ii], input[ii+1], input[ii+2]
        output[k]     = b64enc_tbl[  (u >> 2)                     + 1]
        output[k + 1] = b64enc_tbl[(((u << 4) | (v >> 4)) & 0x3f) + 1]
        output[k + 2] = b64enc_tbl[(((v << 2) | (w >> 6)) & 0x3f) + 1]
        output[k + 3] = b64enc_tbl[  (w                   & 0x3f) + 1]
        k += 4
    end

    if n % 3 == 1
        output[end - 3] = b64enc_tbl[ (input[end] >> 2)         + 1]
        output[end - 2] = b64enc_tbl[((input[end] << 4) & 0x3f) + 1]
        output[end - 1] = base64_pad
        output[end]     = base64_pad
    elseif n % 3 == 2
        output[end - 3] = b64enc_tbl[  (input[end - 1] >> 2)      + 1]
        output[end - 2] = b64enc_tbl[(((input[end - 1] << 4) |
                                       (input[end] >> 4)) & 0x3f) + 1]
        output[end - 1] = b64enc_tbl[ ((input[end] << 2) & 0x3f)  + 1]
        output[end] = base64_pad
    end

    output
end

function decode(::Type{Base64}, input::Vector{UInt8})
    n = length(input)
    if n == 0
        return Array(UInt8, 0)
    end

    if n % 4 != 0
        warn("Length of Base64 input is not a multiple of four.")
    end

    m = 3 * div(n,  4)
    if input[end] == base64_pad
        m -= 1
    end
    if input[end - 1] == base64_pad
        m -= 1
    end
    output = Array(UInt8, m)

    k = 1
    @inbounds for ii = 1:4:length(input)-4
        # This loop is performance-critical, so inline base64dec and consolidate into a single branch point for error-checking
        u = b64dec_tbl[input[ii]]
        v = b64dec_tbl[input[ii + 1]]
        w = b64dec_tbl[input[ii + 2]]
        z = b64dec_tbl[input[ii + 3]]
        if u | v | w | z == sentinel
            if u == sentinel
                error("Invalid base64 symbol: $(char(input[ii]))")
            elseif v == sentinel
                error("Invalid base64 symbol: $(char(input[ii + 1]))")
            elseif w == sentinel
                error("Invalid base64 symbol: $(char(input[ii + 2]))")
            else
                error("Invalid base64 symbol: $(char(input[ii + 3]))")
            end
        end
        output[k]     = (u << 2) | (v >> 4)
        output[k + 1] = (v << 4) | (w >> 2)
        output[k + 2] = (w << 6) | z
        k += 3
    end

    u = input[end-3]
    v = input[end-2]
    w = input[end-1]
    z = input[end]
    u = base64dec(u)
    v = base64dec(v)

    if w == base64_pad && z == base64_pad
        output[end] = (u << 2) | (v >> 4)
    elseif z == base64_pad
        w = base64dec(w)
        output[end - 1] = (u << 2) | (v >> 4)
        output[end]     = (v << 4) | (w >> 2)
    else
        w = base64dec(w)
        z = base64dec(z)
        output[end - 2] = (u << 2) | (v >> 4)
        output[end - 1] = (v << 4) | (w >> 2)
        output[end]     = (w << 6) | z
    end

    output
end

# Zlib/Gzip

abstract Zlib <: Codec

const Z_NO_FLUSH      = 0
const Z_PARTIAL_FLUSH = 1
const Z_SYNC_FLUSH    = 2
const Z_FULL_FLUSH    = 3
const Z_FINISH        = 4
const Z_BLOCK         = 5
const Z_TREES         = 6

const Z_OK            = 0
const Z_STREAM_END    = 1
const Z_NEED_DICT     = 2
const ZERRNO          = -1
const Z_STREAM_ERROR  = -2
const Z_DATA_ERROR    = -3
const Z_MEM_ERROR     = -4
const Z_BUF_ERROR     = -5
const Z_VERSION_ERROR = -6

if is_unix()
    const libz = "libz"
elseif is_windows()
    const libz = "zlib1"
end

# The zlib z_stream structure.
type z_stream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in::Culong

    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out::Culong

    msg::Ptr{UInt8}
    state::Ptr{Void}

    zalloc::Ptr{Void}
    zfree::Ptr{Void}
    opaque::Ptr{Void}

    data_type::Cint
    adler::Culong
    reserved::Culong

    function z_stream()
        strm = new()
        strm.next_in   = C_NULL
        strm.avail_in  = 0
        strm.total_in  = 0
        strm.next_out  = C_NULL
        strm.avail_out = 0
        strm.total_out = 0
        strm.msg       = C_NULL
        strm.state     = C_NULL
        strm.zalloc    = C_NULL
        strm.zfree     = C_NULL
        strm.opaque    = C_NULL
        strm.data_type = 0
        strm.adler     = 0
        strm.reserved  = 0
        strm
    end
end

function zlib_version()
    ccall((:zlibVersion, libz), Ptr{UInt8}, ())
end

function encode(::Type{Zlib}, input::Vector{UInt8}, level::Integer)
    if !(1 <= level <= 9)
        error("Invalid zlib compression level.")
    end

    strm = z_stream()
    ret = ccall((:deflateInit_, libz),
                Int32, (Ptr{z_stream}, Int32, Ptr{UInt8}, Int32),
                &strm, level, zlib_version(), sizeof(z_stream))

    if ret != Z_OK
        error("Error initializing zlib deflate stream.")
    end

    strm.next_in = input
    strm.avail_in = length(input)
    strm.total_in = length(input)
    output = Array(UInt8, 0)
    outbuf = Array(UInt8, 1024)
    ret = Z_OK

    while ret != Z_STREAM_END
        strm.avail_out = length(outbuf)
        strm.next_out = outbuf
        flush = strm.avail_in == 0 ? Z_FINISH : Z_NO_FLUSH
        ret = ccall((:deflate, libz),
                    Int32, (Ptr{z_stream}, Int32),
                    &strm, flush)
        if ret != Z_OK && ret != Z_STREAM_END
            error("Error in zlib deflate stream ($(ret)).")
        end

        if length(outbuf) - strm.avail_out > 0
            append!(output, outbuf[1:(length(outbuf) - strm.avail_out)])
        end
    end

    ret = ccall((:deflateEnd, libz), Int32, (Ptr{z_stream},), &strm)
    if ret == Z_STREAM_ERROR
        error("Error: zlib deflate stream was prematurely freed.")
    end

    output
end

function encode(::Type{Zlib}, input::AbstractString, level::Integer)
    encode(Zlib, convert(Vector{UInt8}, input), level)
end

encode(::Type{Zlib}, input::Vector{UInt8}) = encode(Zlib, input, 9)

function decode(::Type{Zlib}, input::Vector{UInt8})
    strm = z_stream()
    ret = ccall((:inflateInit_, libz),
                Int32, (Ptr{z_stream}, Ptr{UInt8}, Int32),
                &strm, zlib_version(), sizeof(z_stream))

    if ret != Z_OK
        error("Error initializing zlib inflate stream.")
    end

    strm.next_in = input
    strm.avail_in = length(input)
    strm.total_in = length(input)
    output = Array(UInt8, 0)
    outbuf = Array(UInt8, 1024)
    ret = Z_OK

    while ret != Z_STREAM_END
        strm.next_out = outbuf
        strm.avail_out = length(outbuf)
        ret = ccall((:inflate, libz),
                    Int32, (Ptr{z_stream}, Int32),
                    &strm, Z_NO_FLUSH)
        if ret == Z_DATA_ERROR
            error("Error: input is not zlib compressed data.")
        elseif ret != Z_OK && ret != Z_STREAM_END
            error("Error in zlib inflate stream ($(ret)).")
        end

        if length(outbuf) - strm.avail_out > 0
            append!(output, outbuf[1:(length(outbuf) - strm.avail_out)])
        end
    end

    ret = ccall((:inflateEnd, libz), Int32, (Ptr{z_stream},), &strm)
    if ret == Z_STREAM_ERROR
        error("Error: zlib inflate stream was prematurely freed.")
    end

    output
end

# Packed binary-coded decimal (BCD) integers
# Two decimal digits per byte, each represented with 4 bits

abstract BCD <: Codec

function encode(::Type{BCD}, i::Integer, bigendian::Bool = false)
    if i < 0
        error("Integer must not be negative")
    end
    ndig = ndigits(i)
    ndig += isodd(ndig)
    v = digits(i, 10, ndig)
    nbytes = @compat trunc(Integer, ndig/2)
    out = Array(UInt8, nbytes)
    dj = 1-2*bigendian
    for i = 1:nbytes
        j = bigendian*length(v) + dj*2*(i-1) + 1-bigendian
        out[i] = (v[j]<<4) | v[j+dj]
    end
    return out
end

function decode(::Type{BCD}, v::Vector{UInt8}, bigendian::Bool = false)
    ret = 0
    if bigendian
        for i = 1:length(v)
            ret = 100*ret + 10 * (@compat Int(v[i]>>>4)) + (@compat Int(v[i]&0x0f))
        end
    else
        for i = length(v):-1:1
            ret = 100*ret + (@compat Int(v[i]>>>4)) + 10 * (@compat Int(v[i]&0x0f))
        end
    end
    return ret
end

end # module Codecs


text:
module ITensorGPU

using CUDA
using CUDA.CUTENSOR
using CUDA.CUBLAS
using CUDA.CUSOLVER
using LinearAlgebra
using Random, Strided
using TimerOutputs
using StaticArrays
using ITensors
using ITensors.NDTensors
using Strided
import CUDA: CuArray, CuMatrix, CuVector, cu
import CUDA.CUTENSOR: cutensorContractionPlan_t, cutensorAlgo_t

import CUDA.Mem: pin
#=
const devs = Ref{Vector{CUDAdrv.CuDevice}}()
const dev_rows = Ref{Int}(0)
const dev_cols = Ref{Int}(0)
function __init__()
  voltas    = filter(dev->occursin("V100", CUDAdrv.name(dev)), collect(CUDAdrv.devices()))
  pascals    = filter(dev->occursin("P100", CUDAdrv.name(dev)), collect(CUDAdrv.devices()))
  devs[] = voltas[1:1]
  #devs[] = pascals[1:2]
  CUBLASMG.cublasMgDeviceSelect(CUBLASMG.mg_handle(), length(devs[]), devs[])
  dev_rows[] = 1
  dev_cols[] = 1
end
=#
import ITensors:
  randn!,
  compute_contraction_labels,
  eigen,
  tensor,
  scale!,
  unioninds,
  array,
  matrix,
  vector,
  polar,
  tensors,
  truncate!,
  leftlim,
  rightlim,
  permute,
  BroadcastStyle,
  Indices
import ITensors.NDTensors:
  similartype,
  ContractionProperties,
  contract!!,
  _contract!!,
  _contract!,
  contract!,
  contract,
  contraction_output,
  UniformDiagTensor,
  CombinerTensor,
  contraction_output_type,
  UniformDiag,
  Diag,
  DiagTensor,
  NonuniformDiag,
  NonuniformDiagTensor,
  zero_contraction_output,
  outer!,
  outer!!,
  is_trivial_permutation,
  ind,
  permutedims!!,
  Dense,
  DenseTensor,
  Combiner,
  Tensor,
  data,
  getperm,
  compute_contraction_properties!,
  Atrans,
  Btrans,
  Ctrans,
  _contract_scalar!,
  _contract_scalar_noperm!
import Base.*, Base.permutedims!
import Base: similar
include("tensor/cudense.jl")
include("tensor/dense.jl")
include("tensor/culinearalgebra.jl")
include("tensor/cutruncate.jl")
include("tensor/cucombiner.jl")
include("tensor/cudiag.jl")
include("cuitensor.jl")
include("mps/cumps.jl")
include("mps/cumpo.jl")

#const ContractionPlans = Dict{String, Tuple{cutensorAlgo_t, cutensorContractionPlan_t}}()
const ContractionPlans = Dict{String,cutensorAlgo_t}()

export cu,
  cpu, cuITensor, randomCuITensor, cuMPS, randomCuMPS, productCuMPS, randomCuMPO, cuMPO

end #module


text:
#########################################
###Remove duplicates from a sorted vector
#########################################
function uniquesorted(in::AbstractVector)
  out = similar(in)
  c::Int = 0
  if length(in) > 0
    out[c+=1] = in[1]
    @inbounds for i=2:length(in)
      if in[i]!=in[i-1]
        out[c+=1]=in[i]
      end
    end
  end
  out[1:c]
end

############################################################
###Calculate a log-normal value for given location and scale
############################################################
@inline function randlognormal!{T<:AbstractFloat}(v::Vector{T},μ::T,σ::T)
    #tested different variations and the way this is implemented runs fastest
    #there is about 30% benefit in speed of using Float32 over Float64 for large vectors
    l = length(v)
    @simd for i=1:l
        @inbounds v[i] = T(randn())
    end
    @simd for i=1:l
        @inbounds v[i] = μ + σ*v[i]
    end
    @simd for i=1:l
        @inbounds v[i] = exp(v[i])
    end
    v
end


text:
using ConstraintSolver, JuMP, MathOptInterface, JSON

if !@isdefined CS
    const CS = ConstraintSolver
end
const MOI = MathOptInterface
const MOIU = MOI.Utilities
include("../../test/sudoku_fcts.jl")
# include("../../visualizations/plot_search_space.jl")

function parseJSON(json_sums)
    sums = []
    for s in json_sums
        indices = Tuple[]
        for ind in s["indices"]
            push!(indices, tuple(ind...))
        end

        push!(sums, (result = s["result"], indices = indices, color = s["color"]))
    end
    return sums
end

function solve_all(filenames; benchmark = false, single_times = true)
    ct = time()
    for (i, filename) in enumerate(filenames)
        sums = parseJSON(JSON.parsefile("data/$(filename)"))

        # plot_killer(zeros(Int, (9,9)), sums, filename; fill=false)
        # continue

        m = CS.Optimizer(logging = [], time_limit = 20)

        x = [[MOI.add_constrained_variable(m, MOI.Integer()) for i in 1:9] for j in 1:9]
        for r in 1:9, c in 1:9
            MOI.add_constraint(m, x[r][c][1], MOI.GreaterThan(1.0))
            MOI.add_constraint(m, x[r][c][1], MOI.LessThan(9.0))
        end

        for s in sums
            saf = MOI.ScalarAffineFunction{Float64}(
                [MOI.ScalarAffineTerm(1.0, x[ind[1]][ind[2]][1]) for ind in s.indices],
                0.0,
            )
            MOI.add_constraint(m, saf, MOI.EqualTo(convert(Float64, s.result)))
            # MOI.add_constraint(m, [x[ind[1]][ind[2]][1] for ind in s.indices], CS.CPE.AllDifferent(length(s.indices)))
        end

        # sudoku constraints
        moi_add_sudoku_constr!(m, x)

        if single_times
            GC.enable(false)
            MOI.optimize!(m)
            status = MOI.get(m, MOI.TerminationStatus())
            GC.enable(true)
            println(i - 1, ", ", MOI.get(m, MOI.SolveTime()))
        else
            GC.enable(false)
            MOI.optimize!(m)
            status = MOI.get(m, MOI.TerminationStatus())
            GC.enable(true)
        end
        if !benchmark
            println("Status: ", status)
            var_x = fill(MOI.VariableIndex(0), (9, 9))
            for r in 1:9
                var_x[r, :] = [x[r][c][1] for c in 1:9]
            end
            if status == MOI.OPTIMAL
                solution = zeros(Int, 9, 9)
                for r in 1:9
                    solution[r, :] =
                        [MOI.get(m, MOI.VariablePrimal(), x[r][c][1]) for c in 1:9]
                end
                @assert jump_fulfills_sudoku_constr(solution)
            else
                println("NOT SOLVED TO OPTIMALITY")
            end
        end
        com = nothing
        GC.gc()
    end
    println("")
    tt = time() - ct
    println("total time: ", tt)
    println("avg: ", tt / length(filenames))
end

function main(; benchmark = false, single_times = true)
    solve_all(
        [
            "niallsudoku_5500",
            # "niallsudoku_5501",
            # "niallsudoku_5502",
            # "niallsudoku_5503",
            # "niallsudoku_6417",
            # "niallsudoku_6249",
        ];
        benchmark = benchmark,
        single_times = single_times,
    )
    # solve_all(from_file("hardest.txt"), "hardest")
end


text:
quoted(x) = typof(x) === Quote || typof(x) === Quotenode
unquoted(x) = is_unary_call(x) && isoperator(x[1]) && kindof(x[1]) == CSTParser.Tokens.EX_OR

function remove_ref(x::EXPR)
    if hasref(x) && refof(x) isa Binding && refof(x).refs isa Vector
        for ia in enumerate(refof(x).refs)
            if ia[2] == x
                deleteat!(refof(x).refs, ia[1])
                setref!(x, nothing)
                return
            end
        end
        error()
    end
end

function clear_binding(x::EXPR)
    if bindingof(x) isa Binding
        for r in bindingof(x).refs
            if r isa EXPR
                setref!(r, nothing)
            elseif r isa Binding
                clear_binding(r)
            end
        end
        x.meta.binding = nothing
    end
end
function clear_scope(x::EXPR)
    if hasmeta(x) && scopeof(x) isa Scope
        setparent!(scopeof(x), nothing)
        empty!(scopeof(x).names)
        if typof(x) === CSTParser.FileH && scopeof(x).modules isa Dict && scopehasmodule(scopeof(x), :Base) && scopehasmodule(scopeof(x), :Core)
            m1, m2 = getscopemodule(scopeof(x), :Base), getscopemodule(scopeof(x), :Core)
            empty!(scopeof(x).modules)
            addmoduletoscope!(scopeof(x), m1)
            addmoduletoscope!(scopeof(x), m2)
        else
            scopeof(x).modules = nothing
        end
        if scopeof(x).overloaded !== nothing
            empty!(scopeof(x).overloaded)
        end
    end
end

function clear_ref(x::EXPR)
    if refof(x) isa Binding
        if refof(x).refs isa Vector
            for i in 1:length(refof(x).refs)
                if refof(x).refs[i] == x
                    deleteat!(refof(x).refs, i)
                    break
                end
            end
        end
        setref!(x, nothing)
    elseif refof(x) !== nothing
        setref!(x, nothing)
    end
end
function clear_error(x::EXPR)
    if hasmeta(x) && x.meta.error !== nothing
        x.meta.error = nothing
    end
end
function clear_meta(x::EXPR)
    clear_binding(x)
    clear_ref(x)
    clear_scope(x)
    clear_error(x)
    for a in x
        clear_meta(a)
    end
end

function get_root_method(b, server)
    return b
end

function get_root_method(b::Binding, server, b1=nothing, visited_bindings=Binding[])
    if b.prev === nothing || b == b.prev || !(b.prev isa Binding) || b in visited_bindings
        return b
    end
    push!(visited_bindings, b)
    if b.type == b.prev.type == CoreTypes.Function
        return get_root_method(b.prev, server, b, visited_bindings)
    elseif b.type == CoreTypes.Function && b.prev.type == CoreTypes.DataType
        return b.prev
    else
        return b
    end
end

function retrieve_delayed_scope(x)
    if (CSTParser.defines_function(x) || CSTParser.defines_macro(x)) && scopeof(x) !== nothing
        if parentof(scopeof(x)) !== nothing
            return parentof(scopeof(x))
        else
            return scopeof(x)
        end
    else
        return retrieve_scope(x)
    end
    return nothing
end

function retrieve_scope(x)
    if scopeof(x) !== nothing
        return scopeof(x)
    elseif parentof(x) isa EXPR
        return retrieve_scope(parentof(x))
    end
    return
end

function find_return_statements(x::EXPR)
    rets = EXPR[]
    if CSTParser.defines_function(x)
        find_return_statements(x[3], true, rets)
    end
    return rets
end

function find_return_statements(x::EXPR, last_stmt, rets)
    if last_stmt && !(typof(x) === CSTParser.Block || typof(x) === CSTParser.If || iskw(x))
        push!(rets, x)
        return rets, false
    end

    if typof(x) === CSTParser.Return
        push!(rets, x)
        return rets, true
    end

    for i = 1:length(x)
        _, stop_iter = find_return_statements(x[i], last_stmt && (i == length(x) || (typof(x) === CSTParser.If && typof(x[i]) === CSTParser.Block)), rets)
        stop_iter && break
    end
    return rets, false
end

last_method(b::Binding, visited=Binding[]) = b.next isa Binding && b.next.type === CoreTypes.Function && !(b in visited) ? (push!(visited, b);last_method(b.next, visited)) : b

function find_exported_names(x::EXPR)
    exported_vars = EXPR[]
    for i in 1:length(x[3])
        expr = x[3][i]
        if typof(expr) == CSTParser.Export &&
            for j = 2:length(expr)
                if isidentifier(expr[j]) && hasref(expr[j])
                    push!(exported_vars, expr[j])
                end
            end
        end
    end
    return exported_vars
end

hasreadperm(p::String) = (uperm(p) & 0x04) == 0x04

# check whether a path is in (including subfolders) the julia base dir. Returns "" if not, and the path to the base dir if so.
function _is_in_basedir(path::String)
    i = findfirst(r".*base", path)
    i === nothing && return ""
    path1 = path[i]::String
    !hasreadperm(path1) && return ""
    !isdir(path1) && return ""
    files = readdir(path1)
    if all(f -> f in files, ["Base.jl", "coreio.jl", "essentials.jl", "exports.jl"])
        return path1
    end
    return ""
end

isexportedby(k::Symbol, m::SymbolServer.ModuleStore) = haskey(m, k) && k in m.exportednames
isexportedby(k::String, m::SymbolServer.ModuleStore) = isexportedby(Symbol(k), m)
isexportedby(x::EXPR, m::SymbolServer.ModuleStore) = isexportedby(valof(x), m)
isexportedby(k, m::SymbolServer.ModuleStore) = false

function retrieve_toplevel_scope(x::EXPR)
    if scopeof(x) !== nothing && (CSTParser.defines_module(x) || typof(x) === CSTParser.FileH)
        return scopeof(x)
    elseif parentof(x) isa EXPR
        return retrieve_toplevel_scope(parentof(x))
    else
        @info "Tried to reach toplevel scope, no scope found. Final expression $(typof(x))"
        return nothing
    end
end
retrieve_toplevel_scope(s::Scope) = (CSTParser.defines_module(s.expr) || typof(s.expr) === CSTParser.FileH || !(parentof(s) isa Scope)) ? s : retrieve_toplevel_scope(parentof(s))

# b::SymbolServer.FunctionStore or DataTypeStore
# tls is a top-level Scope (expected to contain loaded modules)
# for a FunctionStore b, checks whether additional methods are provided by other packages
# f is a function that returns `true` if we want to break early from the loop

iterate_over_ss_methods(b, tls, server, f) = false
function iterate_over_ss_methods(b::SymbolServer.FunctionStore, tls::Scope, server, f)
    for m in b.methods
        ret = f(m)
        ret && return true
    end
    if b.extends in keys(getsymbolextendeds(server)) && tls.modules !== nothing
        # above should be modified,
        rootmod = SymbolServer._lookup(b.extends.parent, getsymbolserver(server)) # points to the module containing the initial function declaration
        if rootmod !== nothing && haskey(rootmod, b.extends.name) # check rootmod exists, and that it has the variable
            rootfunc = rootmod[b.extends.name]
            # find extensoions
            if haskey(getsymbolextendeds(server), b.extends) # method extensions listed
                for vr in getsymbolextendeds(server)[b.extends] # iterate over packages with extensions
                    !(SymbolServer.get_top_module(vr) in keys(tls.modules)) && continue
                    rootmod = SymbolServer._lookup(vr, getsymbolserver(server))
                    !(rootmod isa SymbolServer.ModuleStore) && continue
                    if haskey(rootmod.vals, b.extends.name) && (rootmod.vals[b.extends.name] isa SymbolServer.FunctionStore || rootmod.vals[b.extends.name] isa SymbolServer.DataTypeStore)# check package is available and has ref
                        for m in rootmod.vals[b.extends.name].methods #
                            ret = f(m)
                            ret && return true
                        end
                    end
                end
            end
        end
    end
    return false
end

function iterate_over_ss_methods(b::SymbolServer.DataTypeStore, tls::Scope, server, f)
    if b.name isa SymbolServer.VarRef
        bname = b.name
    elseif b.name isa SymbolServer.FakeTypeName
        bname = b.name.name
    end
    for m in b.methods
        ret = f(m)
        ret && return true
    end
    if (bname in keys(getsymbolextendeds(server))) && tls.modules !== nothing
        # above should be modified,
        rootmod = SymbolServer._lookup(bname.parent, getsymbolserver(server), true) # points to the module containing the initial function declaration
        if rootmod !== nothing && haskey(rootmod, bname.name) # check rootmod exists, and that it has the variable
            rootfunc = rootmod[bname.name]
            # find extensoions
            if haskey(getsymbolextendeds(server), bname) # method extensions listed
                for vr in getsymbolextendeds(server)[bname] # iterate over packages with extensions
                    !(SymbolServer.get_top_module(vr) in keys(tls.modules)) && continue
                    rootmod = SymbolServer._lookup(vr, getsymbolserver(server))
                    !(rootmod isa SymbolServer.ModuleStore) && continue
                    if haskey(rootmod.vals, bname.name) && (rootmod.vals[bname.name] isa SymbolServer.FunctionStore || rootmod.vals[bname.name] isa SymbolServer.DataTypeStore)# check package is available and has ref
                        for m in rootmod.vals[bname.name].methods #
                            ret = f(m)
                            ret && return true
                        end
                    end
                end
            end
        end
    end
    return false
end

# CST utils
fcall_name(x::EXPR) = is_call(x) && length(x) > 0 && valof(x[1])
is_call(x::EXPR) = typof(x) === CSTParser.Call
is_unary_call(x::EXPR) = typof(x) === CSTParser.UnaryOpCall && length(x) == 2
is_binary_call(x::EXPR) = typof(x) === CSTParser.BinaryOpCall && length(x) == 3
is_binary_call(x::EXPR, opkind) = is_binary_call(x) && kindof(x[2]) === opkind
is_macro_call(x::EXPR) = typof(x) === CSTParser.MacroCall
is_macroname(x::EXPR) = (typof(x) === CSTParser.MacroName && length(x) == 2) || (is_getfield_w_quotenode(x) && is_macroname(x[3][1]))
is_id_or_macroname(x::EXPR) = isidentifier(x) || is_macroname(x)
is_getfield(x) = x isa EXPR && is_binary_call(x) && kindof(x[2]) == CSTParser.Tokens.DOT
is_getfield_w_quotenode(x) = x isa EXPR && is_binary_call(x) && kindof(x[2]) == CSTParser.Tokens.DOT && typof(x[3]) === CSTParser.Quotenode && length(x[3]) > 0
is_declaration(x::EXPR) = is_binary_call(x) && kindof(x[2]) === CSTParser.Tokens.DECLARATION
is_where(x::EXPR) = typof(x) === CSTParser.WhereOpCall
isnonstdid(x::EXPR) = typof(x) === CSTParser.NONSTDIDENTIFIER
is_kwarg(x::EXPR) = typof(x) === CSTParser.Kw
is_parameters(x::EXPR) = typof(x) === CSTParser.Parameters
is_tuple(x::EXPR) = typof(x) === CSTParser.TupleH
is_curly(x::EXPR) = typof(x) === CSTParser.Curly
is_invis_brackets(x::EXPR) = typof(x) === CSTParser.InvisBrackets
iswherecall(x::EXPR) = typof(x) === CSTParser.WhereOpCall
rem_wheres(x::EXPR) = iswherecall(x) ? rem_wheres(x[1]) : x
rem_wheres_decls(x::EXPR) = iswherecall(x) || CSTParser.isdeclaration(x) ? rem_wheres_decls(x[1]) : x
hasparent(x::EXPR) = parentof(x) isa EXPR

"""
    is_in_fexpr(x::EXPR, f)
Check whether `x` isa the child of an expression for which `f(parent) == true`.
"""
is_in_fexpr(x::EXPR, f) = f(x) || (parentof(x) isa EXPR && is_in_fexpr(parentof(x), f))

"""
    get_in_fexpr(x::EXPR, f)
Get the `parent` of `x` for which `f(parent) == true`. (is_in_fexpr should be called first.)
"""
get_parent_fexpr(x::EXPR, f) = f(x) ? x : get_parent_fexpr(parentof(x), f)

issigoffuncdecl(x::EXPR) = parentof(x) isa EXPR ? issigoffuncdecl(x, parentof(x)) : false
function issigoffuncdecl(x::EXPR, p::EXPR)
    if CSTParser.is_where(p) || CSTParser.isdeclaration(p)
        return issigoffuncdecl(parentof(p))
    elseif CSTParser.defines_function(p)
        return true
    else
        return false
    end
end
issigoffuncdecl(x::EXPR, p) = false


text:
module ObservableCollector

import Base: push!

export  timeseries,
        @at,
        @every,
        @condition,
        @observations

function Base.push!(D::Dict{Symbol, T}, m::NamedTuple{(:name, :val), Tuple{Symbol, V}}) where {T,V}
    if !haskey(D, m.name)
        push!(D, m.name => Union{V, Missing}[m.val])
    else
        push!(D[m.name], m.val)
    end
    D
end

"""
    timeseries(x)

Turn a series of named tuples `(name=x, val=y)` into a dictionary with one series for each `name`.
Types are `Union{Missing, T}` where `T` is the type of the first value encountered
for a given name.
"""
function timeseries(x::AbstractArray)
    D = Dict{Symbol,Any}()
    for m in x
        if !haskey(D, m.name)
            push!(D, m.name => Union{typeof(m.val), Missing}[m.val])
        else
            push!(D[m.name], m.val)
        end
    end
    # Flatten all arrays of length 1
    for d in D
        if length(d[2])==1
            D[d[1]] = d[2][1]
        end
    end
    D
end

"""
    @condition(cond, block)

Define multiple observables to be evaluated under the same condition. Use `-->`
to define pairs of `name --> map`.

# Example
```julia
@condition (s,t)->t>=3 begin
    "S" --> (s,t)->s^2
    "T" --> (s,t)->t
end
```
"""
macro condition(condition, block)
    typeof(block)==Expr && block.head == :block || throw(ArgumentError("b in @condition(a, b) must be a block."))
    subexpr = filter(x->(typeof(x)==Expr && x.head==:-->), block.args)

    esc(quote
        local observables = []
        for ob in $subexpr
            name, func = eval.(ob.args)
            typeof(name) == String && typeof(func) <: Function || throw(ArgumentError("String and function expected."))
            name = Symbol(name)
            push!(observables, (name = name, map=func))
        end
        (condition=$condition, observables=observables)
    end)
end

"""
    @condition(cond, name, map)

Creates a clause to collect an observable `map(state,time)` called `name` under
condition `cond(state,time)`.

See also: [`@at`](@ref), [`@every`](@ref)
"""
macro condition(cond, name, map)
    esc(quote
        @condition $cond begin $name --> $map end
    end)
end

"""
    @at(N::Integer, name, map)

Creates a clause to collect an observable `map(state,time)` called `name` at
timestep `N`.

See also: [`@condition`](@ref)
"""
macro at(args...)
    esc(quote
        @condition (s,t)->t==$(args[1]) $(args[2]) $(args[3])
    end)
end

"""
    @every(N::Integer, name, map)

Creates a clause to collect an observable `map(state,time)` called `name` every
`N` timesteps, i.e `t%N==0`.

See also: [`@condition`](@ref)
"""
macro every(args...)
    esc(quote
        @condition (s,t)->t % $(args[1])==0 $(args[2]) $(args[3])
    end)
end

"""
    @observations(block)

Creates an anonymous function with arguments `(output::Vector{Any}, state, time)`
    that collects observables under conditions defined in `block`.

# Example:
```julia
output = []
obs = @observations begin
    @condition (s,t)->t>=3 begin
        "A" --> (s,t)->s
        "T" --> (s,t)->t
    end
end
for t in 1:4
    obs(output,1,t)
end

julia> @show output
output = Any[(name = :A, val = 1), (name = :T, val = 3), (name = :A, val = 1), (name = :T, val = 4)]
4-element Array{Any,1}:
 (name = :A, val = 1)
 (name = :T, val = 3)
 (name = :A, val = 1)
 (name = :T, val = 4)
```
"""
macro observations(ex::Expr)
    bodyparts = filter(x->typeof(x)==Expr && x.head==:macrocall,ex.args)

    esc(quote
        local clauses = eval.($bodyparts)
        local func = :(;)
        local inner = :(;)
        for clause in clauses
            if haskey(clause, :observables)
                inner = :(;)
                for m in clause.observables
                    inner = quote
                        $inner;
                        push!(output, (name=($m).name,val=(($m).map)(state,t)))
                    end
                end
            elseif haskey(clause, :map)
                inner = :( push!(output, (name=($clause).name,val=(($clause).map)(state,t))) )
            end
            func = quote
                $func;
                if (($clause).condition)(state,t)
                    $inner;
                end
            end
        end
        eval( :( (output, state, t) -> begin $func; nothing end ) )
    end)

end

end