text:
############################################################
# joPDEfull ################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export joPDEfull, joPDEfullException, subsample

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

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections

struct joPDEfull{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
	info::Info
	model::ModelAll
	srcGeometry::Geometry
	recGeometry::Geometry
	options::Options
	fop::Function              # forward
	fop_T::Nullable{Function}  # transpose
end

mutable struct joPDEfullException <: Exception
	msg :: String
end

############################################################
## Constructor

function joModeling(info::Info,model::ModelAll, srcGeometry::Geometry, recGeometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
	(DDT == Float32 && RDT == Float32) || throw(joPDEfullException("Domain and range types not supported"))

	# Determine dimensions
	if typeof(recGeometry) == GeometryOOC
		m = sum(recGeometry.nsamples)
	else
		m = 0
		for j=1:info.nsrc m += length(recGeometry.xloc[j])*recGeometry.nt[j] end

	end
	if typeof(srcGeometry) == GeometryOOC
		n = sum(srcGeometry.nsamples)
	else
		n = 0
		for j=1:info.nsrc n += length(srcGeometry.xloc[j])*srcGeometry.nt[j] end
	end

	if info.nsrc > 1
		srcnum = 1:info.nsrc
	else
		srcnum = 1
	end

	return F = joPDEfull{Float32,Float32}("Proj*F*Proj'", m, n, info, model, srcGeometry, recGeometry, options,
							  src -> time_modeling(model, srcGeometry, src.data, recGeometry, [], [], srcnum, 'F', 1, options),
							  rec -> time_modeling(model, srcGeometry, [], recGeometry, rec.data, [], srcnum, 'F', -1, options),
							  )
end

############################################################
## overloaded Base functions

# conj(joPDEfull)
conj{DDT,RDT}(A::joPDEfull{DDT,RDT}) =
	joPDEfull{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		get(A.fop),
		A.fop_T
		)

# transpose(joPDEfull)
transpose{DDT,RDT}(A::joPDEfull{DDT,RDT}) =
	joPDEfull{DDT,RDT}("Proj*F'*Proj'",A.n,A.m,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		get(A.fop_T),
		A.fop
		)

# ctranspose(joPDEfull)
ctranspose{DDT,RDT}(A::joPDEfull{DDT,RDT}) =
	joPDEfull{DDT,RDT}("Proj*F'*Proj'",A.n,A.m,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		get(A.fop_T),
		A.fop
		)

############################################################
## overloaded Base *(...joPDEfull...)

# *(joPDEfull,joData)
function *{ADDT,ARDT,vDT}(A::joPDEfull{ADDT,ARDT},v::joData{vDT})
	A.n == size(v,1) || throw(joPDEfullException("shape mismatch"))
	if compareGeometry(A.srcGeometry,v.geometry) == false && compareGeometry(A.recGeometry,v.geometry) == false
		throw(joPDEfullException("geometry mismatch"))
	end
	jo_check_type_match(ADDT,vDT,join(["DDT for *(joPDEfull,joData):",A.name,typeof(A),vDT]," / "))
	V = A.fop(v)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joPDEfull,joData):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(num,joPDEfull)
function *{ADDT,ARDT}(a::Number,A::joPDEfull{ADDT,ARDT})
	return joPDEfull{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		v1 -> jo_convert(ARDT,a*A.fop(v1),false),
		v2 -> jo_convert(ADDT,a*A.fop_T(v2),false)
		)
end

############################################################
## overloaded Bases +(...joPDEfull...), -(...joPDEfull...)

# +(joPDEfull,num)
function +{ADDT,ARDT}(A::joPDEfull{ADDT,ARDT},b::Number)
	return joPDE{ADDT,ARDT}("("*A.name*"+N)",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		v1 -> A.fop(v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
		v2 -> get(A.fop_T)(v2)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
		)
end

# -(joPDEfull,num)
function -{ADDT,ARDT}(A::joPDEfull{ADDT,ARDT},b::Number)
	return joPDE{ADDT,ARDT}("("*A.name*"-N)",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		v1 -> A.fop(v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
		v2 -> get(A.fop_T)(v2)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
		)
end

# -(joPDEfull)
-{DDT,RDT}(A::joPDEfull{DDT,RDT}) =
	joPDEfull{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
		v1 -> -A.fop(v1),
		v2 -> -get(A.fop_T)(v2)
		)

############################################################
## Additional overloaded functions

# Subsample Modeling operator
function subsample{ADDT,ARDT}(F::joPDEfull{ADDT,ARDT}, srcnum)

	srcGeometry = subsample(F.srcGeometry,srcnum)		# Geometry of subsampled data container
	recGeometry = subsample(F.recGeometry,srcnum)

	info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
	return joModeling(info, F.model, srcGeometry, recGeometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

getindex(F::joPDEfull,a) = subsample(F,a)

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
using Documenter, Nabla

makedocs(
    modules=[Nabla],
    format=:html,
    pages=[
        "Home" => "index.md",
        "API" => "pages/api.md",
        "Custom Sensitivities" => "pages/custom.md",
        "Details" => "pages/autodiff.md",
    ],
    sitename="Nabla.jl",
    authors="Invenia Labs",
    assets=[
        "assets/invenia.css",
    ],
)

deploydocs(
    repo = "github.com/invenia/Nabla.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing,
)

ext: jl


text:
@testset "GammaLikelihood" begin
    lik = GammaLikelihood(1.)
    test_interface(lik, SqExponentialKernel(), rand(10); functor_args=(:α, :invlink))
end


text:
"""
Save file .JSON of the boundingbox in path.
"""
function save_AABB2json(folder::String, aabb::AABB)
	@assert isdir(folder) "save_AABB2json: $path not a valid directory"
	name = splitdir(folder)[2]
	filename = name*".json"
	scale = DataStructures.OrderedDict{String,Any}("x"=>aabb.x_max-aabb.x_min, "y"=>aabb.y_max-aabb.y_min, "z"=>aabb.z_max-aabb.z_min)
	position = DataStructures.OrderedDict{String,Any}("x"=>(aabb.x_max+aabb.x_min)/2, "y"=>(aabb.y_max+aabb.y_min)/2, "z"=>(aabb.z_max+aabb.z_min)/2)
	rotation = DataStructures.OrderedDict{String,Any}("x"=>0., "y"=>0., "z"=>0.)
	data = DataStructures.OrderedDict{String,Any}("clip"=>true, "name"=>name,
			"scale"=>scale,"position"=>position,"rotation"=>rotation,
			"permitExtraction"=>true)
	open(joinpath(folder,filename),"w") do f
  		JSON.print(f, data,4)
	end
end

# """
# {
#    "object":"plane",
#    "normal":{
# 	  "x":0.000,
# 	  "y":0.000,
# 	  "z":1.000
#    },
#    "position":{
#       "x":0.0,
#       "y":0.0,
#       "z":0.0
#    },
# }
# """
# function plane2json(plane::Plane, filename::String)
# 	pos = plane.centroid
# 	dir = plane.normal
# 	position = DataStructures.OrderedDict{String,Any}("x"=>pos[1], "y"=>pos[2], "z"=>pos[3])
# 	normal = DataStructures.OrderedDict{String,Any}("x"=>dir[1], "y"=>dir[2], "z"=>dir[3])
# 	data = DataStructures.OrderedDict{String,Any}("object"=> "plane", "position"=>position,"normal"=>normal)
# 	open(filename,"w") do f
#   		JSON.print(f, data,4)
# 	end
# end

# """
# ============================================  CAMERA
# """
# function cameraparameters(path::String)
# 	dict=Dict{String,Any}[]
# 	open(path, "r") do f
# 	    dict = JSON.parse(f)  # parse and transform data
# 	end
# 	position = dict["position"]
# 	target = dict["target"]
# 	return position, target
# end

#function viewcoordinatesystem(camera)
#     position,target = camera
# 	up = [0,0,1.]
# 	dir = target-position
# 	x = -dir/Lar.norm(dir)
# 	if x != [0,0,1] && x != [0,0,-1]
# 		y = Lar.cross(x,up)
# 		z = Lar.cross(y,x)
# 	end
#     return [y';z';x']
# end
#
#
# function coordsystemcamera(file::String)
#     mat = PointClouds.cameramatrix(file)
#     return Matrix(mat[1:3,1:3]')
# end

# """
# camera parameters from JSON.
# """
# function cameramatrix(path::String)
# 	dict = Dict{String,Any}[]
# 	open(path, "r") do f
# 	    dict = JSON.parse(f)  # parse and transform data
# 	end
# 	mat = dict["object"]["matrix"]
# 	return [mat[1] mat[5] mat[9] mat[13];
# 			mat[2] mat[6] mat[10] mat[14];
# 			mat[3] mat[7] mat[11] mat[15];
# 			mat[4] mat[8] mat[12] mat[16]]
# end
#

# ====================================  read points for polygon
# """
# extract verteces from area tools.
# """
# function vertspolygonfromareaJSON(file::String)
# 	dict = Dict{String,Any}[]
# 	open(file, "r") do f
# 	    dict = JSON.parse(f)  # parse and transform data
# 	end
# 	features = dict["features"]
# 	for feature in features
# 		type = feature["geometry"]["type"]
# 		if type == "Polygon"
# 			points = feature["geometry"]["coordinates"]
# 			V = hcat(points[1][1:end-1]...)
# 			return V
# 		end
# 	end
# end
#
#
# """
# create polygon model.
# """
# function polygon(file::String)
# 	verts = vertspolygonfromareaJSON(file)
# 	EV = [[i,i+1] for i in 1:size(verts,2)-1]
# 	push!(EV,[size(verts,2),1])
# 	axis,centroid = PointClouds.planefit(verts)
# 	if Lar.dot(axis,Lar.cross(verts[:,1]-centroid,verts[:,2]-centroid))<0
# 		axis = -axis
# 	end
# 	PointClouds.projectpointson(verts,(axis,centroid),"plane")
# 	return verts,EV
# end


text:
"""
    legendre(x, b)

Compute Legendre polynomials.

Legendre polynomials are Jacobi polynomials with parameters a and b equals to zero.

Same interface as `jacobi`. `dlegendre` computes the derivative of the Legendre
polynomial.

# Example
```julia-repl
julia> x = legendre(0.4, 5)
0.27063999999999994

julia> x = legendre(2//5, 5)
3383//12500

julia> x = dlegendre(0.4, 5)
-1.3170000000000002

julia> x = dlegendre(2//5, 5)
-1317//1000
```
"""
function legendre(x, n)

    if n==0
        return one(x)
    elseif n==1
        return(x)
    end

    p0 = one(x)
    p1 = x

    for i = 2:n
        p2 = ( (2i-1)*x*p1 - (i-1)*p0 ) / i
        p0 = p1
        p1 = p2
    end

    return p1
end

dlegendre(x, n) = djacobi(x, n)

legendre_zeros(m) = jacobi_zeros(m)

@doc (@doc legendre) dlegendre
@doc (@doc jacobi_zeros!) legendre_zeros


text:
#module snake_forward

#export parallel_forward

using Dierckx

function grouper(arr)
    arr = append!(arr,arr[1])
    result = Array[]
    group = []
    i = 1
    while i < length(arr)
        group = [i]
        i = i+1
        if  i <= length(arr) && arr[i-1] < arr[i]
            while i <= length(arr) && arr[i-1] < arr[i]
                if i == length(arr)
                    push!(group,1)
                else
                    push!(group,i)
                end
                i = i+1
            end
            push!(result, group)
            i = i-1
        elseif i <= length(arr) && arr[i-1] > arr[i]
            while i <= length(arr) && arr[i-1] > arr[i]
                if i == length(arr)
                    push!(group,1)
                else
                    push!(group,i)
                end
                i = i+1
            end
            push!(result, group)
            i = i-1
        end
    end
    if result[end][end] == result[1][1] && sign(result[end][end-1]-result[end][end]) == sign(result[1][1]-result[1][2])
        e = pop!(result)
        pop!(e)
        prepend!(result[1],e)
    end

    return result
end

function piece_wise_spline(splines, bins)
    values = zeros(Float64, length(bins))
    for i = 1:length(bins)
        b = bins[i]
        values[i] = sum(map(spl -> spl(b), splines))
    end
    return values
end

function project_curve(vertices, angle, bins)
    projection = [cos(angle) sin(angle)]'; # projection direction [along with detector]
    normal = [-projection[2] projection[1]]'; # normal direction [along with rays]

    # preparing for work
    vertex_coordinates = (vertices*projection)[:,1]; # vertex coordinates in projection
    distance = vertices*normal;# distances from vertex to projection
    groups = grouper(vertex_coordinates)
    projection = zeros(Float64, length(bins))
    splines = []
    for i = 1:length(groups)
        coordinates_ordered = vertex_coordinates[groups[i]]
        heights = distance[groups[i]]
        if coordinates_ordered[1] > coordinates_ordered[2]
            reverse!(coordinates_ordered)
            heights = -reverse(heights)
        end
        spl= Spline1D(coordinates_ordered,heights, k=1, bc="zero")
        push!(splines, spl)
    end
    return piece_wise_spline(splines, bins)
end

#forward model of list of vertices
function parallel_forward(vertices,angles,bins)
    N = size(vertices,1);
    indices = 1:N
    sinogram = zeros(Float64, length(bins),length(angles));
    vertex_coordinates = zeros(Float64,N,length(angles));
    position = zeros(Float64,N,1); # pre-allocating

    for k = 1:length(angles)

        angle = angles[k]
        projection = project_curve(vertices, angle, bins)

        sinogram[:,k] = projection
    end

    return sinogram
end

#end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
# mostly copied from Kafka.jl
# primitives, strings and arrays

writeobj(io::IO, n::Integer) = write(io, hton(n))
readobj{T<:Integer}(io::IO, ::Type{T}) = ntoh(read(io, T))

function writeobj(io::IO, s::ASCIIString)
    len = Int16(length(s))
    writeobj(io, len > 0 ? len : -1)
    write(io, s)
end
function readobj(io::IO, ::Type{ASCIIString})
    len = readobj(io, Int16)
    return len > 0 ? bytestring(readbytes(io, len)) : ""
end

function writeobj{T}(io::IO, arr::Vector{T})
    len = Int32(length(arr))
    writeobj(io, len)
    for x in arr
        writeobj(io, x)
    end
end
function readobj{T}(io::IO, ::Type{Vector{T}})
    len = readobj(io, Int32)
    arr = Array(T, len)
    for i=1:len
        arr[i] = readobj(io, T)
    end
    return arr
end

# composite types
function writeobj(io::IO, o)
    for f in fieldnames(o)
        writeobj(io, getfield(o, f))
    end
end
function readobj{T}(io::IO, ::Type{T})
    vals = Array(Any, length(T.types))
    for (i, t) in enumerate(T.types)
        vals[i] = readobj(io, t)
    end
    return T(vals...)
end


text:
"""
This module provides:
- An extension of the `AbstractArray` interface in order to properly deal with mutable caches.
- A mechanism to generate lazy arrays resulting from operations between arrays.
- A collection of concrete implementations of `AbstractArray`.

The exported names in this module are:

$(EXPORTS)
"""
module Arrays

using Gridap.Helpers
using Gridap.Inference
using Gridap.Algebra

using DocStringExtensions
using Test
using FillArrays
using Base: @propagate_inbounds
using LinearAlgebra
using ForwardDiff
using BlockArrays

export BlockArrayCoo
export BlockVectorCoo
export BlockMatrixCoo
export is_zero_block
export is_nonzero_block
export enumerateblocks
export eachblockindex
export VectorOfBlockArrayCoo
export VectorOfBlockVectorCoo
export VectorOfBlockMatrixCoo
export zeros_like
export TwoLevelBlockedUnitRange

export array_cache
export getindex!
export getitems!
export getitems
export testitem
export uses_hash
export test_array
export testitems
export array_caches
export get_array
export get_arrays
export add_to_array!

export CachedArray
export CachedMatrix
export CachedVector
export setsize!
export setaxes!

export CompressedArray
export LocalToGlobalArray
export LocalToGlobalPosNegArray
export FilteredCellArray
export FilterKernel

export kernel_cache
export kernel_caches
export apply_kernels!
export apply_kernel!
export apply_kernel
export test_kernel
export bcast
export elem
export contract
export MulKernel
export MulAddKernel
export kernel_return_type
export kernel_return_types
export kernel_testitem
export Kernel

export apply
export apply_all

export Table
export identity_table
export empty_table
export rewind_ptrs!
export length_to_ptrs!
export append_ptrs
export append_ptrs!
export get_ptrs_eltype
export get_data_eltype
export generate_data_and_ptrs
export find_inverse_index_map
export find_inverse_index_map!
export append_tables_globally
export append_tables_locally
export flatten_partition
export collect1d
export UNSET
export get_local_item
export find_local_index

export reindex
export identity_vector

export SubVector
export pair_arrays
export unpair_arrays

export lazy_append
export lazy_split
export AppendedArray

export autodiff_array_gradient
export autodiff_array_jacobian
export autodiff_array_hessian

export VectorWithEntryRemoved
export VectorWithEntryInserted

import Base: size
import Base: getindex, setindex!
import Base: similar
import Base: IndexStyle

import Gridap.Io: to_dict
import Gridap.Io: from_dict

import Gridap.Algebra: scale_entries!

include("Interface.jl")

include("BlockArraysCoo.jl")

include("CachedArrays.jl")

include("Kernels.jl")

include("Apply.jl")

include("CompressedArrays.jl")

include("Tables.jl")

include("LocalToGlobalArrays.jl")

include("LocalToGlobalPosNegArrays.jl")

include("FilteredArrays.jl")

include("Reindex.jl")

include("IdentityVectors.jl")

include("SubVectors.jl")

include("ArrayPairs.jl")

include("AppendedArrays.jl")

include("VectorsOfBlockArrayCoo.jl")

include("Autodiff.jl")

include("VectorsWithEntryRemoved.jl")

include("VectorsWithEntryInserted.jl")

end # module


text:
@doc raw"""
    thermometer(x::Float64, n::Int)
    thermometer(x::Vector{Float64}, n::Int)
"""
function thermometer(x::Float64, n::Int)
    y = zeros(Bool, n)
    k = trunc(Int, clamp(x, 0.0, 1.0) * n)
    y[1:k] .= true
    return y
end

function thermometer(x::Vector{Float64}, n::Int)
    m = length(x)
    y = zeros(Bool, m * n)
    for i = 1:m
        j = (i - 1) * n + 1
        k = trunc(Int, clamp(x[i], 0.0, 1.0) * n) - 1
        y[j:j+k] .= true
    end
    return y
end


text:
"""
# CInclude.jl

Include C language header files in Julia source code
(using [Clang.jl](https://github.com/JuliaInterop/Clang.jl)).

    julia> using CInclude
    julia> @cinclude "termios.h" quiet

    [ Info: @cinclude "/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/termios.h"
    ...

    julia> term = termios()
    julia> settings.c_lflag = ICANON
    julia> settings.c_cflag = CS8
    julia> cfsetspeed(Ref(settings), B38400)
    julia> tcsetattr(fd, TCSANOW, Ref(settings))
"""
module CInclude

export @cinclude

using ReadmeDocs
using Clang

function system_include_path()
    path = ["/usr/include"]
    if Sys.isapple()
        sdk = chomp(read(`xcrun --show-sdk-path`, String))
        push!(path, joinpath(sdk, "usr/include"))
    elseif Sys.islinux()
        try
            x = eachline(`sh -c "gcc -xc -E -v /dev/null 2>&1"`)
            line, state = iterate(x)
            while line != nothing &&
                  line != "#include <...> search starts here:"
                line, state = iterate(x, state)
            end
            line, state = iterate(x, state)
            while line != nothing &&
                  line != "End of search list."
                push!(path, strip(line))
                line, state = iterate(x, state)
            end
        catch err
            @warn err
        end
    end
    path
end

function find_header(header)
    header = replace(header, r"[<>]" => "")
    if !isfile(header)
        for d in system_include_path()
            h = joinpath(d, header)
            if isfile(h)
                return h
            end
        end
    end
    header
end

function macro_values(headers, names)
    headers = [startswith(h, "<") ? h : "\"$h\"" for h in headers]
    mktempdir() do d
        cfile = joinpath(d, "tmp.cpp")
        delim = "aARf6F3fWe6"
        write(cfile, """
            #include <iostream>
            #include <string>
            #include <iomanip>
            #include <typeinfo>
            $(join(["#include $h" for h in headers], "\n"))

            #define T(x) typeid(x).name()[0]
            #define isstr(x) (T(x) == 'A')
            #define ischr(x) (T(x) == 'c')
            #define isint8(x) (T(x) == 'h'|| T(x) == 'a')
            #define dump(n,x) std::cout << "const " << n << " = " << x
            #define jlquote(x) std::quoted((const char*)(size_t)x)
            #define jlchar(x) std::string("Char(") +                         \\
                              std::to_string((size_t)(x)) + ")"

            #define wrap(x) ({                                               \\
                std::cout << "$delim";                                       \\
                     if(isstr(x))  dump(#x, jlquote(x));                     \\
                else if(ischr(x))  dump(#x,  jlchar(x));                     \\
                else if(isint8(x)) dump(#x,  size_t(x));                     \\
                else               dump(#x,        (x));                     \\
            })

            int main() {
                $(join(["wrap($n);" for n in names], "\n"))
            }
            """)
        #println(read(cfile, String))
        write("cinclude_tmp.c", read(cfile, String))
        binfile = joinpath(d, "tmp.bin")
        try
            run(`g++ -w -o $binfile $cfile`)
            output = read(`$binfile`, String)
            map(Meta.parse, split(output, delim; keepempty=false))
        catch err
            @error err
            return []
        end
    end
end

function wrap_headers(headers; lib="libc", include="", exclude=r"^!", args=[])

    long_headers = map(find_header, headers)
    @info "@cinclude $long_headers"

    ctx = DefaultContext()
    ctx.libname = lib
    ctx.options["is_function_strictly_typed"] = false
    ctx.options["is_struct_mutable"] = false

    cargs::Vector{String} = vcat((["-I", d] for d in system_include_path())...)
    append!(cargs, args)
    @info cargs

    parse_headers!(ctx, long_headers, args=cargs)

    macro_names = []

    for unit in ctx.trans_units
        ctx.children = children(getcursor(unit))
        for (i, child) in enumerate(ctx.children)
            ctx.force_name = ""
            ctx.children_index = i
            child_name = name(child)

            # choose which cursor to wrap
            if haskey(ctx.common_buffer, child_name) || # skip already wrapped
               (exclude != "" && occursin(exclude, child_name) &&
               (include == "" || !occursin(include, child_name)) &&
               !(child isa Clang.CLEnumDecl))

                continue
            end

            if child isa Clang.CLEnumDecl && child_name == ""
                ctx.anonymous_counter += 1
                ctx.force_name = "ANONYMOUS_ENUM_$(ctx.anonymous_counter)"
            end

            if child isa Clang.CLMacroDefinition
                if startswith(child_name, "_")
                    continue
                end
                tokens = try tokenize(child) catch end
                if tokens != nothing &&
                   tokens.size > 1 &&
                   tokens[2].text != "(" &&
                   tokens[2].text != "{" &&
                   (tokens.size < 3 || tokens[3].text != ".")
                    push!(macro_names, child_name)
                    continue
                end
            end

            try
                wrap!(ctx, child)
            catch err
                if !startswith(child_name, "_")
                    @info "Can't wrap $child_name ($err)"
                end
            end
        end
    end

    template = [
        :(import Base),
        :(using Base:|),
        :(using CEnum),
        :(const Ctm = Base.Libc.TmStruct),
        :(const Ctime_t = UInt),
        :(const Cclock_t = UInt),
        :(const Cstring = Base.Cstring),
        :(const Clong = Base.Clong)]

    for t in values(Clang.CLANG_JULIA_TYPEMAP)
        if isdefined(Base, t)
            push!(template, :(const $t = Base.$t))
        end
    end

    api = [template;
           dump_to_buffer(ctx.common_buffer);
           ctx.api_buffer;
           macro_values(headers, unique(macro_names))]

    constructors = []

    for e in api
        if !(e isa Expr)
            continue
        end

        # Remove library name for built-in functions.
        # `ccall((:cfunc, libc), ...` => `ccall(:func, ...`)
        if e.head == :function
           args = e.args[2].args[1].args
           if args[2].args[2] == :libc
               args[2] = args[2].args[1]
           end
        end

        # Generate zero-value default constructors.
        if e.head == :struct
            for T in [e.args[2], Symbol("$(e.args[2])_m")]
                push!(constructors, :(
                    function $T()
                        $T(CInclude.czeros($T)...)
                    end))
            end
        end
    end

    [api; constructors]
end

function symbol_name(e::Expr)
    if e.head in (:using, :const, :function)
        return e.args[1].args[1]
    elseif e.head == :struct
        return e.args[2]
    elseif e.head == :macrocall
        return e.args[3].args[1]
    else
        dump(e)
    end
    return nothing
end

czeros(T) =
    (hasmethod(zero, (t,)) ? zero(t) :
                t <: Tuple ? (czeros(t)...,) :
                            t(czeros(t)...)
     for t in fieldtypes(T))

README"## Interface"

eval_int_constant(mod, name) = Base.eval(mod, :(try Int($name) catch end))
eval_enum_constant(mod, name) = Base.eval(mod, :(CEnum.namemap($name)))

README"""
    @cinclude "header.h" [quiet] [exclude=""] [include=""]
    @cinclude ["foo.h" "bar.h"] [quiet] [exclude=""] [include=""]

Import symbols from C language `header.h` into the current module.

Note: C language headers may define a large number of symbols. To avoid
clashes with local names `@cinclude` can be used inside a sub-module. e.g.

    module SOCKET
        using CInclude
        @cinclude "sys/socket.h"
    end
    SOCKET.send(fd, buf, length(buf), SOCKET.MSG_OOB)
"""
macro cinclude(h, options...)
    if :quiet in options
        options=filter(x->x!=:quiet, options)
        logger=:NullLogger
    else
        logger=:current_logger
    end
    static = :static in options
    options=filter(x->x!=:static, options)
    verbose = :verbose in options
    options=filter(x->x!=:verbose, options)
    for o in options
        o.args[1] = esc(o.args[1])
    end
    if h isa String
        h = [h]
    end
    quote

        static_code = String[]
        function add_static_code(e)
            $static || return
            push!(static_code, string(e))
        end

        constants = Dict{Int,Vector{Symbol}}()
        function add_constant(n::Symbol, v::Integer)
            v = convert(Int, v)
            names = get(constants, v, nothing)
            if names == nothing
                constants[v] = [n]
            else
                push!(names, n)
            end
        end

        function doceval(e)
            if VERSION >= v"1.6"
                if e.head in (:function, :const)
                    doc = string("```\n", e, "\n```\n")
                    e = :(Core.@doc($doc, $e))
                end
            end
            Base.eval(@__MODULE__, e)
        end

        Base.CoreLogging.with_logger(Base.CoreLogging.$logger()) do
            for e in CInclude.wrap_headers($h; $(options...))
                if e isa String
                    @info e
                else
                    $verbose && @info e
                    # Check for duplicate methods.
                    if e.head == :function
                        args = e.args[1].args
                        f = args[1]
                        if Base.isdefined(@__MODULE__, f)
                            args = args[2:end]
                            if all(a->a isa Symbol, args)
                                types = ((Any for a in args)...,)
                                f = Base.eval(@__MODULE__, f)
                                if !isempty(methods(f, types))
                                    @info "Skipping duplicate $f$types"
                                    continue
                                end
                            end
                        end
                    end
                    try
                        doceval(e)
                        add_static_code(e)
                        # Create `_m` mutable struct variant.
                        if e.head == :struct
                            e = copy(e)
                            e.args[2] = Symbol("$(e.args[2])_m")
                            e.args[1] = true
                            Base.eval(@__MODULE__, e)
                            add_static_code(e)
                        end
                    catch err
                        #exception=(err, Base.catch_backtrace())
                        @warn "$err in $e" #exception
                    end
                    if e.head == :macrocall && e.args[1] == Symbol("@cenum")
                        name = e.args[3].args[1]
                        for (v, n) in eval_enum_constant(@__MODULE__, name)
                            add_constant(n, v)
                        end
                    elseif e.head == :const
                        name = e.args[1].args[1]
                        v = eval_int_constant(@__MODULE__, name)
                        if v != nothing
                            add_constant(name, v)
                        end
                    end
                end
            end
            e = Expr(:const, Expr(Symbol("="), :constants, constants))
            Base.eval(@__MODULE__, e)

            add_static_code(e)
            @info "Saving output in $(pwd())"
            write("c_include_static_output.jl", join(static_code, "\n\n"))
            nothing
        end
    end
end

end # module


text:
#= If running in console:
using Test
=#
using Cytof5, Flux, Distributions, RCall
import Random
using BSON
@rimport graphics
@rimport grDevices

@testset "Variational Inference" begin
  println("Test Variational Inference...")
  # Set random seed for generating data
  Random.seed!(2)

  # Generate Data
  L = Dict(0=>3, 1=>3)
  J = 20
  K = 5

  println("Simulate Data...")
  N = [8, 1, 2] * 5000

  I = length(N)
  mus=Dict(0 => -[1.0, 2.3, 3.5],
           1 => +[1.0, 2.0, 3.0])
  a_W=rand(K)*10
  a_eta=Dict(z => rand(L[z])*10 for z in 0:1)
  Z = Cytof5.Model.genSimpleZ(J, K)
  @time dat = Cytof5.Model.genData(J=J, N=N, K=K, L=L, Z=Z,
                                   beta=[-9.2, -2.3],
                                   sig2=[0.2, 0.1, 0.3],
                                   mus=mus,
                                   a_W=a_W,
                                   a_eta=a_eta,
                                   sortLambda=false,
                                   propMissingScale=0.7,
                                   eps=fill(.005, I))

  tau = .0001 # FIXME: why can't I do .001?
  use_stickbreak = false
  noisy_var = 10.0
  K_MCMC = 10
  L_MCMC = Dict(false=>5, true=>3)
  priors = Cytof5.VB.Priors(K=K_MCMC, L=L_MCMC, use_stickbreak=use_stickbreak)
  mc = Cytof5.Model.defaultConstants(Cytof5.Model.Data(dat[:y]),
                                     K_MCMC, Dict{Int64,Int64}(L_MCMC),
                                     yQuantiles=[0.0, 0.25, 0.5],
                                     pBounds=[.05, .8, .05])

  beta = [mc.beta[:, i] for i in 1:I]
  c = Cytof5.VB.Constants(I, N, J, K_MCMC, L_MCMC, tau, beta, use_stickbreak,
                          noisy_var, priors)
  y = dat[:y]

  println("fit model...")
  # niters = 20000
  # batchsize = 2000
  # to test compilation
  niters = 100
  batchsize = 200
  opt = ADAM(1e-2)
  out = Cytof5.VB.fit(y=y, niters=niters, batchsize=batchsize, c=c, opt=opt,
                      seed=0, nsave=30)

  # Save results
  mkpath("result/vb")
  @time BSON.bson("result/vb/vb-out.bson", out)

  #= Load results
  using Cytof5, Flux, Distributions, BSON
  =#
  out = BSON.load("result/vb/vb-out.bson")
  samples = [Cytof5.VB.rsample(out[:state])[2] for b in 1:100]
  Z = [reshape(s.v, 1, length(s.v)) .> s.H for s in samples]

  grDevices.pdf("result/vb/results.pdf")
    graphics.plot(out[:metrics][:elbo] / sum(out[:c].N),
                  typ="l", xlab="iterations", ylab="elbo");
    graphics.image(mean(Z)' .> .5)
  grDevices.dev_off()

  @test true
end


text:
abstract type SlidingVector{T,P <: AbstractArray{T}} <: AbstractVector{T} end
mutable struct SlidingVec{T,P} <: SlidingVector{T, P}
    parent::P
    offset::Int
    l::Int
end
struct ChildVec{T,P<: SlidingVector{T}} <: SlidingVector{T,P}
    parent::P
    offset::Int
    l::Int
end
struct SlidingVecFun{V, F}
    f::F
    v::V
end

Base.IndexStyle(::Type{<: SlidingVector}) = IndexLinear()
function Base.getindex(V::SlidingVector, i::Int)
    @inbounds r = V.parent[V.offset + i]
    r
end
function Base.setindex!(V::SlidingVector, x, i::Int)
    @inbounds V.parent[V.offset + i] = x
    V
end

function Base.getindex(V::ChildVec, i::Int)
    @inbounds r = V.parent.parent[V.offset + V.parent.offset + i]
    r
end
function Base.setindex!(V::ChildVec, x, i::Int)
    @inbounds V.parent.parent[V.offset + V.parent.offset + i] = x
    V
end
Base.size(A::SlidingVector) = (A.l, )
Base.length(A::SlidingVector) = A.l
function SlidingVec(n::P, i::Int = 0, l::Int = 1) where {T, P <: AbstractArray{T}}
    SlidingVec{T, P}(n, i, l)
end
function SlidingVec(n::P, i::Int = 0, l::Int = size(n,1)) where {T, P <: AbstractMatrix{T}}
    SlidingVec{T, P}(n, i, l)
end
function SlidingVec(n::P, i::UnitRange{Int}) where {T, P <: AbstractArray{T,2}}
    SlidingVec{T, P}(n, i[1], length(i))
end
function view(A::SlidingVec{T, P}, i::UnitRange{Int}) where {T, P}
    ChildVec{T, P}(A, i[1], length(i))
end
function slide!(A::SlidingVec)
    A.offset += A.l
end
function reset!(A::SlidingVec)
    A.offset = 0
end
function set!(A::SlidingVecFun, M::AbstractMatrix)
    A.v.parent = M
    A.v.offset = - A.v.l
end
function eval!(A::SlidingVecFun)
    A.v.offset += A.v.l
    A.f()
end
function (A::SlidingVecFun)(v::AbstractVecOrMat)
    A.v.parent[1+A.v.offset:length(v)+A.v.offset] .= v
    A.f()
end
function (A::SlidingVecFun)(v::SVector{p}) where p
    @inbounds for (i, v_i) ∈ enumerate(v)
        A.v.parent[i + A.v.offset] = v_i
    end
    A.v.parent[1+A.v.offset:A.v.l+A.v.offset] .= v
    A.f()
end
@inline (A::SlidingVecFun)() = eval!(A)
@inline eval!(A::SlidingVecFun, v::AbstractVector) = A(v)
#@inline (A::SlidingVecFun)(i::Int) = getindex(A, i)
function getindex(A::SlidingVecFun, i::Int)
    A.v.offset = A.v.l*(i-1)
    A.f()
end
function setindex!(A::SlidingVecFun, v::AbstractVector, i::Int)
    A.v.offset = A.v.l*(i-1)
    A.v.parent[1+A.v.offset:A.v.l+A.v.offset] .= v
    A.f()
end
function SlidingVecFun(f::Function, M::AbstractMatrix)
    l = size(M, 1)
    v = SlidingVec(M, -l, l)
    g() = f(v)
    SlidingVecFun(g, v)
end

struct randIter
    n::Int
end
function Base.start(A::SlidingVecFun)
    A.v.offset = 0
    A.f()
end
function Base.next(A::SlidingVecFun, state)
    eval!(A), state + 1
end
function Base.done(A::SlidingVecFun, state)
    state >= size(A.v.parent, 2)
end
function Base.eltype(A::SlidingVecFun)
    Core.Inference.return_type(A.f, ())
end
function Base.endof(A::SlidingVecFun)
    eval!(A, size(A.v.parent, 2))
end
Base.length(A::SlidingVecFun) = size(A.v.parent, 2)
Base.size(A::SlidingVecFun) = ( size(A.v.parent, 2), )


text:
struct Amphipod
    kind::Symbol
    position::Tuple{Int, Int}
    energy_used::Int
end

struct Burrow
    amphipods::Vector{Amphipod}
    spaces::Vector{Tuple{Int, Int}}
end

function cost(amphipod::Amphipod)
    amphipod.kind == :A && return 1
    amphipod.kind == :B && return 10
    amphipod.kind == :C && return 100
    amphipod.kind == :D && return 1000
end

function homecol(amphipod::Amphipod)
    amphipod.kind == :A && return 3
    amphipod.kind == :B && return 5
    amphipod.kind == :C && return 7
    amphipod.kind == :D && return 9
end

function ishome(amphipod::Amphipod, burrow::Burrow)
    row, col = amphipod.position
    (col != homecol(amphipod) || row == 1) && return false
    others = [get(burrow, (r, col)) for (r, _) in burrow.spaces if r > row]
    for other in others
        !isnothing(other) && other.kind != amphipod.kind && return false
    end
    true
end

function isvacant(burrow::Burrow, position::Tuple{Int, Int})
    position in burrow.spaces || return false
    for amphipod in burrow.amphipods
        amphipod.position == position && return false
    end
    true
end

function Base.get(burrow::Burrow, position::Tuple{Int, Int})
    for amphipod in burrow.amphipods
        amphipod.position == position && return amphipod
    end
    nothing
end

function Base.parse(::Type{Burrow}, lines::Vector{String})
    amphipods = Amphipod[]
    spaces = Tuple{Int, Int}[]
    for (row, line) in enumerate(lines[2:end-1])
        chars = collect(line[2:end])  # ignore first char so col begins counting at start of hallway
        for (col, char) in enumerate(chars)
            (char == '#' || char == ' ') && continue
            char == '.' && push!(spaces, (row, col))
            if char in ['A', 'B', 'C', 'D']
                push!(spaces, (row, col))
                push!(amphipods, Amphipod(Symbol(char), (row, col), 0))
            end
        end
    end
    Burrow(amphipods, spaces)
end

function isinhallway(amphipod::Amphipod)
    amphipod.position[1] == 1
end

function canenterhallway(amphipod::Amphipod, burrow::Burrow)
    isinhallway(amphipod) && return false
    row, col = amphipod.position
    all(isvacant(burrow, (r, col)) for (r, _) in burrow.spaces if r < row)
end

function canwalkhallway(amphipod::Amphipod, burrow::Burrow)
    _, col = amphipod.position
    home_col = homecol(amphipod)
    home_col == col && return true
    step = sign(home_col - col)
    all(isvacant(burrow, (1, c)) for c in col+step:step:home_col)
end

function canenterhome(amphipod::Amphipod, burrow::Burrow)
    others = [get(burrow, (r, homecol(amphipod))) for (r, _) in burrow.spaces if r > 1]
    all(isnothing(other) || other.kind == amphipod.kind for other in others)
end

function cangohome(amphipod::Amphipod, burrow::Burrow)
    (canenterhallway(amphipod, burrow) || isinhallway(amphipod)) &&
        canwalkhallway(amphipod, burrow) &&
        canenterhome(amphipod, burrow)
end

function moves(amphipod::Amphipod, burrow::Burrow)
    row, col = amphipod.position
    moves = Amphipod[]
    ishome(amphipod, burrow) && return moves
    if cangohome(amphipod, burrow)
        home_col = homecol(amphipod)
        home_row = sum(isvacant(burrow, (r, c)) for (r, c) in burrow.spaces if c == home_col)
        energy_used = amphipod.energy_used + cost(amphipod) * (row - 1 + abs(home_col - col) + home_row - 1)
        push!(moves, Amphipod(amphipod.kind, (home_row, home_col), energy_used))
    elseif canenterhallway(amphipod, burrow)
        energy_used = amphipod.energy_used + cost(amphipod) * (row - 1)
        left_energy_used = energy_used
        for left_col in col-1:-1:1
            !isvacant(burrow, (1, left_col)) && break
            left_energy_used += cost(amphipod)
            left_col in [3, 5, 7, 9] && continue
            push!(moves, Amphipod(amphipod.kind, (1, left_col), left_energy_used))
        end
        right_energy_used = energy_used
        for right_col in col+1:11
            !isvacant(burrow, (1, right_col)) && break
            right_energy_used += cost(amphipod)
            right_col in [3, 5, 7, 9] && continue
            push!(moves, Amphipod(amphipod.kind, (1, right_col), right_energy_used))
        end
    end
    moves
end

function totalenergy(burrow::Burrow)
    sum(amphipod.energy_used for amphipod in burrow.amphipods)
end

function isorganized(burrow::Burrow)
    all(ishome(amphipod, burrow) for amphipod in burrow.amphipods)
end

function moves(burrow::Burrow)
    burrows = Burrow[]
    for (i, amphipod) in enumerate(burrow.amphipods)
        for move in moves(amphipod, burrow)
            amphipods_copy = copy(burrow.amphipods)
            splice!(amphipods_copy, i, [move])
            push!(burrows, Burrow(amphipods_copy, burrow.spaces))
        end
    end
    burrows
end

function organize!(burrow::Burrow, energy::Vector{Int})
    burrows = moves(burrow)
    isempty(burrows) && return
    for b in burrows
        if isorganized(b)
            push!(energy, totalenergy(b))
        else
            organize!(b, energy)
        end
    end
end

lines = readlines("day23/burrow.txt")

burrow = parse(Burrow, lines)

energy = Int[]
@time organize!(burrow, energy)

p1 = minimum(energy)
@show p1

insert!(lines, 4, "  #D#C#B#A#")
insert!(lines, 5, "  #D#B#A#C#")

burrow = parse(Burrow, lines)

energy = Int[]
@time organize!(burrow, energy)

p2 = minimum(energy)
@show p2


text:
#=
run_spineopt_test:
- Julia version:
- Author: jhjiang
- Date: 2020-10-01
=#
# syntax to run this script

using SpineOpt
using IJulia

if isdefined(Main, :IJulia) && Main.IJulia.inited
    IJulia.set_max_stdio(1 << 25)
end

input_db_url = ARGS[1]
output_db_url = ARGS[2]
@time begin
    m = run_spineopt(input_db_url, output_db_url; cleanup=true, optimize=true)
end

# Show active variables and constraints
println("*** Active constraints: ***")
for key in keys(m.ext[:constraints])
    !isempty(m.ext[:constraints][key]) && println(key)
end
println("*** Active variables: ***")
for key in keys(m.ext[:variables])
    !isempty(m.ext[:variables][key]) && println(key)
end


text:
# This is just some experimental code
# to implement EntityEmbeddings for purely
# categorical features

using Flux

mutable struct EmbeddingMatrix
    e
    levels

    function EmbeddingMatrix(levels; dim=4)
        if dim <= 0
            dimension = div(length(levels), 2)
        else
            dimension = min(length(levels), dim)      # Dummy function for now
        end
        return new(Dense(length(levels), dimension), levels), dimension
    end

end

Flux.@treelike EmbeddingMatrix

function (embed::EmbeddingMatrix)(ip)
    return embed.e(Flux.onehot(ip, embed.levels))
end

mutable struct EntityEmbedding
    embeddingmatrix

    function EntityEmbedding(a...)
        return new(a)
    end
end

Flux.@treelike EntityEmbedding

# ip is an array of tuples
function (embed::EntityEmbedding)(ip)
    return hcat((vcat((embed.embeddingmatrix[i](ip[idx][i]) for i=1:length(ip[idx]))...) for idx =1:length(ip))...)
end

#   Q1. How should this be called in the API?
#   nn = NeuralNetworkClassifier(builder=builder, optimiser = .., embeddingdimension = 5)
#
#
#


text:
function update_manifests(path::AbstractString = pwd();
                          registries::Vector{Pkg.Types.RegistrySpec} = default_registries,
                          delete_old_manifest::Bool = false)
    environments_to_update = Vector{String}(undef, 0)
    for (root, dirs, files) in walkdir(path)
        for file in files
            if (file == "Manifest.toml") | (file == "JuliaManifest.toml")
                environment = joinpath(path, root)
                push!(environments_to_update, environment)
            end
        end
    end
    unique!(environments_to_update)
    _update_manifests(environments_to_update;
                      registries = registries,
                      delete_old_manifest = delete_old_manifest)
    return nothing
end

function _update_manifests(environments::AbstractVector{<:AbstractString};
                           registries::Vector{Pkg.Types.RegistrySpec},
                           delete_old_manifest::Bool)
    for environment in environments
        _update_manifest(environment;
                         registries = registries,
                         delete_old_manifest = delete_old_manifest)
    end
    return nothing
end

function _has_project(environment::AbstractString)
    return (ispath(joinpath(environment, "Project.toml"))) |
           (ispath(joinpath(environment, "JuliaProject.toml")))
end

function _update_manifest(environment::AbstractString;
                          registries::Vector{Pkg.Types.RegistrySpec},
                          delete_old_manifest::Bool)
    @info("Updating environment: $(environment)")
    if _has_project(environment) && delete_old_manifest
        @debug("Removing Manifest.toml if it exists")
        rm(joinpath(environment, "Manifest.toml");
           force = true, recursive = true)
        @debug("Removing JuliaManifest.toml if it exists")
        rm(joinpath(environment, "JuliaManifest.toml");
           force = true, recursive = true)
    end
    with_tmp_dir() do tmp_dir
        code = """
            using Pkg;
            using UUIDs;
            Pkg.Types.RegistrySpec(name::Union{Nothing, String},
                                   uuid::Union{Nothing, UUID},
                                   url::Union{Nothing, String},
                                   path::Union{Nothing, String}) = Pkg.Types.RegistrySpec(;
                                                                                name = name,
                                                                                uuid = uuid,
                                                                                url = url,
                                                                                path = path)
            registries = $(registries);
            for registry in registries;
                Pkg.Registry.add(registry);
            end;
            Pkg.activate("$(environment)"; shared = false);
            Pkg.instantiate();
            Pkg.update();
            """
        cmd = Cmd(`$(Base.julia_cmd()) -e $(code)`;
                  env = Dict("PATH" => ENV["PATH"],
                             "JULIA_DEPOT_PATH" => tmp_dir))
        run(pipeline(cmd, stdout=stdout, stderr=stderr))
        return nothing
    end
    return nothing
end


text:
__precompile__()

module Bifurcations

    using FileIO, Images, Formatting

    include("utils.jl")
    include("populations.jl")
    include("gradients.jl")

    function plot(pop::Population,
            colmap::Gradient;
            background::Colorant=RGB(1,1,1),
            name::String="test",
            iters::Int64=600,
            osa::Bool=true,
            sx::Int64=1080, sy::Int64=1920,
            xmin::Float64=0.0, xmax::Float64=1.0,
            ymin::Float64=0.0, ymax::Float64=4.0)
        # Image size
        if osa
            sx, sy = 2sx, 2sy
        end
        sx, sy = sx-1, sy-1

        # Preallocates image array
        img::Array = fill(background, sx, sy)

        # We precompute colors for each iteration for performance
        colors::Array = [grad(colmap, n, 1, iters/2) for n = 1:iters]

        # Pixels populate
        for pix_y = 1:sy
            p::Float64 = pop.initial
            y::Float64 = interp(pix_y, 1, sy, ymin, ymax)

            old_values = Set{Float64}()
            for n = 1:iters
                # Detects cyclical behaviour
                approx_p::Float64 = round(p, digits=10)
                if in(approx_p, old_values)
                    break
                end
                push!(old_values, approx_p)

                # n > card(f) to avoid initial horizontal lines on the plot
                if xmin < p < xmax && n > length(pop.fertilities)
                    pix_x = trunc(Int64, interp(p, xmin, xmax, sx, 1))
                    img[pix_x, pix_y] = colors[n]
                end
                f = pop.fertilities[n % end + 1](y)
                p = pop.law(f, p)
            end
        end
        if osa
            img = restrict(img)
        end
        save("$name.png", img)
    end

    # Obsolete exemple for plotting animated zooms

    #= function __init__()
        function compute(y1::Real, y2::Real, x1::Real, x2::Real, i::Int)
            println("=== STARTING COMPUTATION ===")
            colormap = Gradient(RGB(1,0,0), RGB(0,0,0))
            pop = Population(logistic, [id, cst(3.8)], .5)
            plot(pop, colormap, name=format("u/{:.5f}-{:.5f}__{:.5f}-{:.5f}", y1, y2, x1, x2), iters=trunc(i), ymin=y1, ymax=y2, xmin=x1, xmax=x2)
            println("Done.")
        end

        # function compute(x::Real)
        #     red2black = Gradient(RGB(1,0,0), RGB(0,0,0))
        #     pop = Population(logistic, [id, cst(x)], .5)
        #     plot(pop, red2black, format("i/{:.5f}", x))
        # end

        function expspace(start::Real, stop::Real, steps::Int)
            logStart = log(start)
            logStop  = log(stop)
            space = linspace(logStart, logStop, steps)
            exp.(space)
        end

        function zoom(lowInit, upInit, fact, frames)
            lows, ups = [], []
            low,  up  = lowInit, upInit
            for n = 1:frames
                push!(lows, low)
                push!(ups,  up)
                range = up - low
                low += range * fact
                up  -= range * fact
            end
            lows, ups
        end

        const FRAMES = 2000
        const Y_LIMITS = zoom(0.99680, 4.00000, 0.025, FRAMES)
        const X_LIMITS = zoom(0.00000, 0.82896, 0.02, FRAMES)

        # crops a window [lower, upper] to make <center> its midpoint
        function limits(lower, center, upper)
            range = upper -lower
            if center < range/2
                return lower, 2 * center
            elseif center > range/2
                return 2 * center - range, upper
            else
                lower, upper
            end

        # pmap(compute, Y_LIMITS[1], Y_LIMITS[2], fill(0.4, FRAMES), fill(0.6, FRAMES), Int64.(trunc.(linspace(600, 6000, FRAMES))))
        pmap(compute, Y_LIMITS[1], Y_LIMITS[2], X_LIMITS[1], X_LIMITS[2], Int64.(trunc.(linspace(600, 100000, FRAMES))))
    end =#

    function __init__()
        # pop = Population(logistic, [id, cst(3.9), cst(3.9), cst(3.9)], .5) # pas mal
        # pop = Population(logistic, [id, id, id, cst(3.75), id, cst(3.75), cst(3.75)], .5) # trop stylé
        # pop = Population(logistic, [id, cst(3.85), cst(3.75)], .5) #vraiment cool
        # pop = Population(logistic, [id, id, id, cst(3.9), id, cst(3.9), cst(3.9)], .5) #celle affichée dans ma chambre
        # pop = Population(logistic, [id, cst(3.8), cst(3.7)], .5)
        pop = Population(logistic, [id], .5)

        black = RGB(0,0,0)
        red   = RGB(1,0,0)
        white = RGB(1,1,1)
        alpha = RGBA(1,1,1,0)

        black_to_red = Gradient(black, red)
        red_to_black = Gradient(red,   black)
        red_to_white = Gradient(red,   white)
        black     = Gradient(black, black)

        # Precomp func for benchmarks

        # plot(pop, colormap; name="precomp", sx=128, sy=128)

        # Usual demo func

        plot(pop, black;
            name = string(rand(0000:9999)),
            iters = 1000,
            osa = true,
            sx = 1080*2, sy = 1920*2,
            xmin = 0.0, xmax = 1.0,
            ymin = 3.6, ymax = 4.0
        )

        # Batch multiple views

        #= n = 22
        for i = 1:n
            println("Computing frame $i/$n...")

            plot(pop, red2black;
                name = format("{:02d}", i),
                iters = 1000,
                osa = true,
                # sx = 297, sy = 420,
                sx = 2970, sy = 4200,
                xmin = 0.0, xmax = 1.0,
                ymin = 4(i-1)/n, ymax = 4i/n
            )
        end =#
    end
end


text:
"""
    This implements the backends
"""

import Base: show

export  backend, Backend,
        # Possibiities
        BackendPETSc,  BackendParallelStencil, BackendNone, @init_backend,

        # MPI, flobal variables
        comm, mpirank, mpisize

"""
    AbstractBackend{type, FT}

Abstract supertype that specifies the backend of the simulation (PETSc.jl, ParallelStencil, None) as well as how
the simulation is being setup (MPI or not)
"""
abstract type AbstractBackend{type, FT} end

"""
    AbstractTypeBackend

Abstract supertype that determines the types of the backends that are implemented
"""
abstract type AbstractTypeBackend end

"""
    BackendPETSc

Use PETSc as a backend
"""
abstract type BackendPETSc <: AbstractTypeBackend end

"""
    BackendParallelStencil

Use ParallelStencil as a backend
"""
abstract type BackendParallelStencil <: AbstractTypeBackend end

"""
    BackendNone

No backend; julia native only. Note that this is not compatible with MPI
"""
abstract type BackendNone <: AbstractTypeBackend end

"""
    Backend(type::Symbol, arch::Symbol, mpi::Bool)

Structure that defines the backend used for the computional grid

- `type`: `:ParallelStencil`, `:PETSc` (requires loading the corresponding packages), or `:Julia` (julia native grid)
- `arch`: `:CPU`,`:GPU`,`:CUDA` (requires `CUDA` to be loaded first)
- `mpi`: Boolean that indicates if we employ MPI

"""
struct Backend{T,FT} <: AbstractBackend{T,FT}
    type   :: Symbol
    arch   :: Symbol
    mpi    :: Bool
    Scalar :: DataType
end

"""
    Backend(; type::Symbol=:PETSc, arch::Symbol=:CPU, mpi::Bool=true)

Specify the backend we are using as well as whether the simulation is done on an MPI parallel machine or not

- `type`: `:ParallelStencil`, `:PETSc` (requires loading the corresponding packages), or `:Julia` (julia native)
- `arch`: `:CPU`,`:GPU`,`:CUDA`
- `mpi`: Boolean that indicates if we employ MPI or not. Note that PETSc always requires MPI
- `dim`: Number of dimensions (only relevant for `ParallelStencil`)

"""
function Backend(; arch::Symbol=:CPU, type::Symbol=:PETSc, mpi::Bool=true) where FT
    global backend

    @show type, typeof(type)

    # Determine the main backend we are using (PETSc/ParallelStencil or None)
    if type == :PETSc
        backend_type = BackendPETSc;
    elseif type == :ParallelStencil
        backend_type = BackendParallelStencil;
    else
        backend_type = BackendNone;
    end

    return Backend{backend_type, FT}(type, arch, mpi, Scalar)
end

# Printing
function show(io::IO, b::Backend{type,FT}) where {type,FT}
    if b.mpi
        println(io, "Backend{$(FT)}: $(b.type) ($(b.arch) | MPI)")
    else
        println(io, "Backend{$(FT)}: $(b.type) ($(b.arch) | no MPI)")
    end
end

"""
    @init_backend(type, arch, mpi=false, Scalar=Float64)

Specify the backend we are using as well as whether the simulation is done on an MPI parallel machine or not

- `type`: `:ParallelStencil`, `:PETSc` (requires loading the corresponding packages), or `:Julia` (julia native)
- `arch`: `:CPU`,`:GPU`,`:CUDA`
- `mpi`: Boolean that indicates if we employ MPI or not. Note that PETSc always requires MPI (and currently only works with CPU)
- `Scalar`: Type of Scalar

"""
macro init_backend(type, arch, mpi::Bool=false, ndims::Int=3, FT=Float64)
    global backend, comm, mpirank, mpisize

    mpirank = 0
    mpisize = 1
    comm    = 0

    if mpi
        # Initialize MPI & make
        if isdefined(Main, :MPI)
            comm = MPI.COMM_WORLD
            if !MPI.Initialized()
                MPI.Init()
            end
            mpirank = MPI.Comm_rank(MPI.COMM_WORLD)
            mpisize = MPI.Comm_size(MPI.COMM_WORLD)
        else
            error("Load MPI package first with using MPI")
        end
    end

    #
    if type == :PETSc
        backend_type = BackendPETSc;
        type_sym = :PETSc
        if !mpi
            error("PETSc backend requires MPI")
        end
    elseif type == :ParallelStencil
        backend_type = BackendParallelStencil;
        type_sym = :ParallelStencil
    else
        backend_type = BackendNone;
        type_sym = :Julia
    end

    if arch == :CUDA
        arch_sym = :CUDA
    elseif arch == :Threads
        arch_sym = :Threads
    elseif arch == :Base
        arch_sym = :CPU
    else
        arch_sym = :CPU
    end

    # Initialize the global backend
    backend = Backend{backend_type, eval(FT)}(type_sym, arch_sym, mpi, eval(FT))

    if type == :PETSc
        return nothing
    elseif type == :ParallelStencil
        return esc(:(ParallelStencil.@init_parallel_stencil($arch, $FT, $ndims);))
    else
        return nothing
    end
end

"""
    localfromglobalsize(mpi,N, opts)

Computes the local grid size given a global grid and an MPI communicator (not taking halo into account).

"""
function localfromglobalsize(mpi,N, opts)
    dim = length(N);
    dims = zeros(Int64,dim)

    if !mpi
        N_l = N
    else
        dims[dim+1:end] .= 1;       # fix
        if typeof(opts) <: Dict # process command-line options
            dims[1] = get(opts, :dimx, dims[1]);
            if dim>1;   dims[2] = get(opts, :dimy, dims[2]);    end
            if dim>2;   dims[3] = get(opts, :dimz, dims[3]);    end
        end

        # Use MPI to compute processor partitioning (dims contain the # of processors in each direction)
        MPI.Dims_create!(mpisize,dims)

        # Compute local from global size
        N_l = N./dims
        if any(abs.(mod.(N_l,1.0)) .> 0)
            error("Grid size cannot be distributed over the given number of MPI ranks; N_l=$N_l; dims=$dims ")
        end

        N_l = (Int64.(N_l)...,)

    end

    return N_l, dims
end


text:
# ---
# title: 688. Knight Probability in Chessboard
# id: problem688
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/knight-probability-in-chessboard/description/>
# hidden: true
# ---
#
# On an `N`x`N` chessboard, a knight starts at the `r`-th row and `c`-th column
# and attempts to make exactly `K` moves. The rows and columns are 0 indexed, so
# the top-left square is `(0, 0)`, and the bottom-right square is `(N-1, N-1)`.
#
# A chess knight has 8 possible moves it can make, as illustrated below. Each
# move is two squares in a cardinal direction, then one square in an orthogonal
# direction.
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/knight.png)
#
#
#
# Each time the knight is to move, it chooses one of eight possible moves
# uniformly at random (even if the piece would go off the chessboard) and moves
# there.
#
# The knight continues moving until it has made exactly `K` moves or has moved
# off the chessboard. Return the probability that the knight remains on the
# board after it has stopped moving.
#
#
#
# **Example:**
#
#
#
#     Input: 3, 2, 0, 0
#     Output: 0.0625
#     Explanation: There are two moves (to (1,2), (2,1)) that will keep the knight on the board.
#     From each of those positions, there are also two moves that will keep the knight on the board.
#     The total probability the knight stays on the board is 0.0625.
#
#
#
#
# **Note:**
#
#   * `N` will be between 1 and 25.
#   * `K` will be between 0 and 100.
#   * The knight always initially starts on the board.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using FastFibration
using Test

@testset "FastFibration.jl" begin
    # Testing vertex obj
    @test FastFibration.Vertex(0).index == 0
    @test length(FastFibration.Vertex(0).edges_source) == 0
    @test length(FastFibration.Vertex(0).edges_target) == 0

    # Testing initialization algorithm for fibration
    # Open files
    graph_test_1 = FastFibration.graph_from_csv("../data/test_1.csv", true)
    graph_test3scc = FastFibration.graph_from_csv("../data/test_3SCC.csv", true)
    graph_test3scc2 = FastFibration.graph_from_csv("../data/test_3SCC_2.csv", true)
    graph_test6scc = FastFibration.graph_from_csv("../data/test_6SCC.csv", true)

    part_test1, pq_test1 = FastFibration.initialize(graph_test_1)
    part_test3scc, pq_test3scc = FastFibration.initialize(graph_test3scc)
    part_test3scc2, pq_test3scc2 = FastFibration.initialize(graph_test3scc2)
    part_test6scc, pq_test6scc = FastFibration.initialize(graph_test6scc)

end


text:
module dlgo

export Player, black, white, void
export other, isvoid, nbrs, zobristhash
export Point, Move, GoString, GoStringBuffer, IdType, Board, GameState,
       Territory, TerritoryStatus, GameResult
export getindex, setindex!, ==, is_on_grid
export new_game, apply_move, is_valid_move, is_over, resign, play, pass_turn, legal_moves
export grid, stringdict, stringbuffer, getpointid, setpointid!, getstring,
       setstring!, place_stone!, remove_string!, remove_point!
export tosymbol, print_board, point_from_coords, print_move
export evaluate_territory, score, winnercolor, compute_game_result
export features

include("gotypes.jl")
include("zobrist.jl")
include("goboard.jl")
include("stringrep.jl")
include("scoring.jl")
include("features.jl")

include("agent/agent.jl")

end


text:
@testset "setup_constraints" begin

FL=64
TF=Float64

comp_grid = compgrid((1.0, 1.0),(100, 201))
comp_grid3D = compgrid((1.0, 1.0, 1.0),(50, 60, 30))

options    = PARSDMM_options()
options.FL = TF

#test project_bounds! with scalar bounds
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = -0.11
  m_max     = 0.01
  set_type  = "bounds"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)

  x = randn(100)
  P_sub[1](x)
  @test maximum(x) <= m_max
  @test minimum(x) >= m_min

#test project_bounds! with vector bounds
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = randn(100) .- 10.0
  m_max     = randn(100) .+ 10.0
  set_type  = "bounds"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)

  x = 100.0 .* randn(100)
  P_sub[1](x)
  [@test x[i] <= m_max[i] for i=1:100]
  [@test x[i] >= m_min[i] for i=1:100]

#test project l1_Duchi!
  #see if projector returns untouched vector
  x=randn(100)
  y=deepcopy(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = norm(x,1)*2
  set_type  = "l1"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test x == y

  #choose smaller constraint than vector l1 norm
  x=randn(100)
  y=deepcopy(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = norm(x,1)*0.234
  set_type  = "l1"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test isapprox(norm(x,1),m_max,rtol=10*eps())

#test project_cardinality!

  #test vector mode
  x=randn(100)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 5
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test count(!iszero, x)==m_max

  #test on closed form solution (vector)
  x = [0;0;-1;2;-3]
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 2
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test x==[0;0;0;2;-3]

  #test mode that projects each column of a matrix separately
  X = randn(comp_grid.n)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 7
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("fiber","x")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](vec(X))
  X = reshape(X,comp_grid.n)
  [@test count(!iszero, X[:,i]) == 7 for i=1:size(X,2)]

  #test mode that projects each row of a matrix separately
  X = randn(comp_grid.n)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 11
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("fiber","z")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](vec(X))
  X = reshape(X,comp_grid.n)
  [@test count(!iszero, X[i,:]) == 11 for i=1:size(X,1)]

  #test mode that projects each fiber of a tensor separately (either in x,y or z direction)
  X = randn(comp_grid3D.n)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 7
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("fiber","x")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](vec(X))
  X = reshape(X,comp_grid3D.n)
  [@test count(!iszero, X[:,i,j]) == 7 for i=1:size(X,2) for j=1:size(X,3)]

  X = randn(comp_grid3D.n)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 6
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("fiber","y")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](vec(X))
  X = reshape(X,comp_grid3D.n)
  [@test count(!iszero, X[i,:,j]) == 6 for i=1:size(X,1) for j=1:size(X,3)]

  X = randn(comp_grid3D.n)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 4
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("fiber","z")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](vec(X))
  X = reshape(X,comp_grid3D.n)
  [@test count(!iszero, X[i,j,:]) == 4 for i=1:size(X,1) for j=1:size(X,2)]

  #test mode that projects onto each slice of a tensor separately (either in x,y or z direction)
  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 7
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("slice","x")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  X = P_sub[1](X)
  X = reshape(X,comp_grid3D.n)
  [@test count(!iszero, X[i,:,:]) == 7 for i=1:comp_grid3D.n[1]]

  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 6
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("slice","y")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  X = P_sub[1](X)
  X = reshape(X,comp_grid3D.n)
  [@test count(!iszero, X[:,i,:]) == 6 for i=1:comp_grid3D.n[2]]

  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 5
  set_type  = "cardinality"
  TD_OP     = "identity"
  app_mode  = ("slice","z")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  X = P_sub[1](X)
  X = reshape(X,comp_grid3D.n)
  [@test count(!iszero, X[:,:,i]) == 5 for i=1:comp_grid3D.n[3]]

# test project_l2!
  #test for sigma smaller than 2-norm of vector
  x=randn(100)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 0.123
  set_type  = "l2"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](x)
  @test isapprox(norm(x,2),0.123,rtol=eps()*10)

  #test for sigma larger than 2-norm of vector
  x=randn(100)
  y=deepcopy(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 1.234*norm(x,2)
  set_type  = "l2"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](x)
  @test x==y

#test project_rank! (codes assumes input and output are vectors)

  # test if input is returned untouched if input matrix rank is less than the constraint
  X = vec(randn(comp_grid.n))
  XX = deepcopy(X)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = minimum(comp_grid.n)
  set_type  = "rank"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](X)
  @test isapprox(XX,X,rtol=20*eps())

  # test if output matrix has correct rank when constraint is smaller than matrix rank
  X = vec(randn(comp_grid.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 12
  set_type  = "rank"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid.n)
  r_X=rank(X)
  @test r_X==12

  #test rank projection on slices of a 3D tensor
  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 7
  set_type  = "rank"
  TD_OP     = "identity"
  app_mode  = ("slice","x")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid3D.n)
  [@test rank(X[i,:,:]) == 7 for i=1:size(X,1)]

  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 6
  set_type  = "rank"
  TD_OP     = "identity"
  app_mode  = ("slice","y")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid3D.n)
  [@test rank(X[:,i,:]) == 6 for i=1:size(X,2)]

  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 5
  set_type  = "rank"
  TD_OP     = "identity"
  app_mode  = ("slice","z")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid3D.n)
  [@test rank(X[:,:,i]) == 5 for i=1:size(X,3)]

#test project nuclear! (codes assumes input and output are vectors)
  # test if input is returned untouched if nuclear norm is less than constraint
  X    = randn(comp_grid.n)
  nn_X = norm(svdvals(X),1)
  X    = vec(X)
  XX   = deepcopy(X)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = nn_X*1.1
  set_type  = "nuclear"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](X)
  @test isapprox(XX,X,rtol=20*eps())

  # test if output matrix has correct nuclear norm when constraint is smaller than matrix nuclear norm
  X    = randn(comp_grid.n)
  nn_X = norm(svdvals(X),1)
  X    = vec(X)
  XX   = deepcopy(X)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = nn_X*0.567
  set_type  = "nuclear"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](X)
  nn_Xp = norm(svdvals(reshape(X,comp_grid.n)),1)
  @test isapprox(nn_Xp,nn_X*0.567,rtol=20*eps())

  #test nuclear norm projection onto slices of a 3D tensor
  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 1.234
  set_type  = "nuclear"
  TD_OP     = "identity"
  app_mode  = ("slice","x")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid3D.n)
  [@test isapprox(norm(svdvals(X[i,:,:]),1), 1.234,rtol=100*eps()) for i=1:size(X,1)]

  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 1.234
  set_type  = "nuclear"
  TD_OP     = "identity"
  app_mode  = ("slice","y")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid3D.n)
  [@test isapprox(norm(svdvals(X[:,i,:]),1), 1.234,rtol=100*eps()) for i=1:size(X,2)]

  X = vec(randn(comp_grid3D.n))
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 1.234
  set_type  = "nuclear"
  TD_OP     = "identity"
  app_mode  = ("slice","z")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid3D,options.FL)
  P_sub[1](X)
  X=reshape(X,comp_grid3D.n)
  [@test isapprox(norm(svdvals(X[:,:,i]),1), 1.234,rtol=100*eps()) for i=1:size(X,3)]

#test project_subspace!
  #on a vector inputs
  #orthogonal subspace
  M=randn(100,50)
  F=svd(M)
  x=randn(100)
  y=deepcopy(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 0.0
  set_type  = "subspace"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = (F.U,true)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test isapprox(x,F.U*(F.U'*y),rtol=eps()*100)

  #non-orthogonal subspace
  M=randn(100,50)
  x=randn(100)
  y=deepcopy(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 0.0
  set_type  = "subspace"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = (M,false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test isapprox(x,M*((M'*M)\(M'*y)),rtol=eps()*10)

  #on matrix input: project every column onto the subspace
  M=randn(prod(comp_grid.n[1]),50)
  x=randn(prod(comp_grid.n))
  y=deepcopy(x)
  x=vec(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 0.0
  set_type  = "subspace"
  TD_OP     = "identity"
  app_mode  = ("fiber","x")
  custom_TD_OP = (M,false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test isapprox(x,vec(M*((M'*M)\(M'*reshape(y,comp_grid.n)))),rtol=eps()*10)

  #on matrix input: project every row onto the subspace
  M=randn(prod(comp_grid.n[2]),50)
  x=randn(prod(comp_grid.n))
  y=deepcopy(x)
  x=vec(x)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = 0
  m_max     = 0.0
  set_type  = "subspace"
  TD_OP     = "identity"
  app_mode  = ("fiber","z")
  custom_TD_OP = (M,false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test isapprox(x,vec((M*((M'*M)\(M'*reshape(y,comp_grid.n)')))'),rtol=eps()*10)

#test projection onto relaxed histogram
  #first test exact histogram projection
  ref = sort(randn(100))
  x   = randn(100)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = ref
  m_max     = ref
  set_type  = "histogram"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  @test ref == sort(x)

  #relaxed histogram projection
  x   = randn(100)
  constraint = Vector{SetIntersectionProjection.set_definitions}()
  m_min     = sort(randn(100))
  m_max     = m_min.+0.7
  set_type  = "histogram"
  TD_OP     = "identity"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  P_sub[1](x)
  x = sort(x)
  [@test x[i] <= m_max[i] for i=1:100]
  [@test x[i] >= m_min[i] for i=1:100]

# #test projection onto l1-ball with DFT operator
# constraint=Dict()
# constraint["use_TD_l1_2"]=true
# constraint["TD_l1_operator_2"]="DFT"
# constraint["TD_l1_sigma_2"] = 1.234
# (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,FL)
#
# (TD_OP, dummy1, dummy2, dummy3)=get_TD_operator(comp_grid,"DFT",TF)
#
# #project vector with ||TD_OP*x||_1 > sigma norm
# x=randn(TF,100,201)
# x=vec(x)
# P_sub[1](x)
# @test isapprox(norm(TD_OP*x,1),1.234,rtol=100*eps(TF))
#
# #project vector with ||TD_OP*x||_1 < sigma norm
# constraint=Dict()
# constraint["use_TD_l1_2"]=true
# constraint["TD_l1_operator_2"]="DFT"
# constraint["TD_l1_sigma_2"] = 10000000.0
# (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,FL)
# x=randn(TF,100,201)
# x=vec(x)
# y=deepcopy(x)
# P_sub[1](x)
# @test isapprox(x,y,rtol=10*eps(TF))

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, CMBMissionSim

makedocs(modules = [CMBMissionSim],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CMBMissionSim.jl",
    pages = [
        "Introduction" => "index.md",
        "Pointing generation" => "genpointings.md"
    ])

deploydocs(repo = "github.com/lspestrip/CMBMissionSim.jl.git",
)


text:
# Autogenerated wrapper script for NOVAS_jll for i686-w64-mingw32
export libnovas

JLLWrappers.@generate_wrapper_header("NOVAS")
JLLWrappers.@declare_library_product(libnovas, "libnovas-3.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnovas,
        "bin\\libnovas-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Base.Filesystem
# using MLJ
# using MLJScikitLearnInterface,
#     MLJDecisionTreeInterface,
#     MLJClusteringInterface,
#     MLJGLMInterface,
#     MLJLIBSVMInterface,
#     MLJMultivariateStatsInterface,
#     MLJXGBoostInterface
using Pandas
using DataFrames
using CSV

const REGRESSION_PRED_COLUMN_NAME = "Predictions"
const RUNNING_LANG_MSG = "Running environment language: Julia"
const CUSTOM_MODEL_FILE_EXTENSION = ".jlso"

global defined_hooks = Dict(
    "init" => false,
    "read_input_data" => false,
    "load_model" => false,
    "transform" => false,
    "score" => false,
    "unstructured_predict" => false,
    "post_process" => false,
)

## initialize

function init(code_dir, target_type)

    ## include custom.jl
    custom_path = Filesystem.joinpath(code_dir, "custom.jl")
    if Filesystem.isfile(custom_path)
        custom_loaded = (
            try
                include(custom_path)
                @eval using Main.Custom
                @info "successfully included $custom_path"
                true
            catch e
                bt = backtrace()
                msg = sprint(showerror, e, bt)
                @error """\nload custom.jl failed"""
                @error msg
                throw(Exception)
            end
        )
    else
        @info "no custom.jl detected"
        custom_loaded = false
    end

    if custom_loaded
        hooks = names(Custom)
        for hook in hooks
            hook_str = string(hook)
            if haskey(defined_hooks, hook_str)
                @info "$hook_str is defined"
                global defined_hooks[hook_str] = true
            end
        end
    end
    if defined_hooks["init"]
        ## no idea why this work :(
        @info "running init from defined hooks with code directory $code_dir"
        # @eval Custom.init($code_dir)
        Base.invokelatest(Custom.init, code_dir)
    end

end

function load_serialized_model(model_dir)

    model = nothing

    files = Filesystem.readdir(model_dir)
    artifacts = [x for x in files if endswith(x, "jlso")]
    if length(artifacts) == 0
        @error """\nCould not find a serialized model artifact with $CUSTOM_MODEL_FILE_EXTENSION extension, supported by default Julia predictor.
        If your artifact is not supported by default predictor, implement Custom.load_model hook."""
        throw(Exception)
    elseif length(artifacts) > 1
        @error """Multiple serialized model artifacts found: $(join(artifacts, ", "))"""
        throw(Exception)
    else
        artifact = artifacts[1]
        artifact_path = Filesystem.joinpath(model_dir, artifact)
        mlj_libs = Filesystem.joinpath(@__DIR__, "mlj_libs.jl")
        @info include(mlj_libs)
        ## thinking about loading necessary MLJ modules instead
        ## of top leve.  should save some time and lighten lib
        ## requirement on user if there is no intention to use MLJ
    end
    if defined_hooks["load_model"]
        @info "Loading model via hook"
        model = Custom.load_model(model_dir)
    else
        @info "Loading serialized model artifact found at $artifact_path"
        model = Main.MLJLibs.MLJ.machine(artifact_path)
    end

    return model
end

function predict_regression(data, model; kwargs...)
    if defined_hooks["score"]
        predictions = Custom.score(data, model; kwargs...)
    else
        predictions = Main.MLJLibs.MLJ.predict(model, data)
    end
    return Pandas.DataFrame(predictions, columns=[REGRESSION_PRED_COLUMN_NAME])
end

function predict_classification(
    data,
    model;
    positive_class_label=nothing,
    negative_class_label=nothing,
    class_labels=nothing,
    target_type=nothing,
)
    @debug "postive_class_label = $positive_class_label"
    @debug "negative_class_label = $negative_class_label"
    @debug "class_labels = $class_labels"
    if !isnothing(positive_class_label) & !isnothing(negative_class_label)
        labels = [positive_class_label, negative_class_label]
    elseif !isnothing(class_labels)
        labels = class_labels
    else
        throw("something is not right")
    end
    yhat = Main.MLJLibs.MLJ.predict(model, data)
    probs = Main.MLJLibs.MLJ.pdf(yhat, labels)
    return Pandas.DataFrame(probs, columns=labels)
end

predictors = Dict(
    "binary" => predict_classification,
    "multiclass" => predict_classification,
    "regression" => predict_regression,
    "anomaly" => predict_regression,
)

function model_predict(data, model; kwargs)
    target_type = kwargs.target_type
    if !haskey(predictors, target_type)
        throw("Target type $target_type is not supported by Julia predictor")
    else
        predictors[target_type](data, model; kwargs...)
    end
end

function outer_predict(
    target_type;
    binary_data=nothing,
    mimetype=nothing,
    model=nothing,
    positive_class_label=nothing,
    negative_class_label=nothing,
    class_labels=nothing,
)

    if defined_hooks["read_input_data"]
        data = Custom.read_input_data(binary_data)
    elseif !isnothing(mimetype) && mimetype == "text/mtx"
        @error "MTX is not yet supported"
        throw(Exception)
    else
        data = CSV.read(
            IOBuffer(binary_data),
            DataFrames.DataFrame,
            delim=","
        )
    end

    if isnothing(model)
        model = load_serialized_model()
    end

    if defined_hooks["transform"]
        data = Custom.transform(data, model)
    end

    kwargs = (
        positive_class_label = positive_class_label,
        negative_class_label = negative_class_label,
        class_labels = class_labels,
    )

    if defined_hooks["score"]
        predictions = Custom.score(data, model; kwargs)
    else
        kwargs = merge(kwargs, (target_type = target_type,))
        predictions = model_predict(data, model; kwargs)
    end

    if defined_hooks["post_process"]
        predictions = Custom.post_process(predictions, model)
    end

    return Pandas.DataFrame(predictions)
end

function predict_unstructured(data; model=nothing, mimetype=nothing, charset=nothing, query=nothing)
    kwargs = (mimetype = mimetype, charset = charset, query = query)
    predictions = Custom.score_unstructured(data; model, kwargs)
    return predictions
end

max_stars_repo_name: andreakropp/datarobot-user-models

max_issues_repo_name: andreakropp/datarobot-user-models

max_forks_repo_name: andreakropp/datarobot-user-models


text:
module MCMCChainsStorage

import Base: read, write

using HDF5
using MCMCChains

"""
    write(f::Union{HDF5.File, HDF5.Group}, chains::Chains)

Write MCMCChains object to the HDF5 file or group.
"""
function Base.write(f::Union{HDF5.File, HDF5.Group}, c::Chains)
    for s in sections(c)
        g = create_group(f, string(s))
        for n in names(c, s)
            g[string(n), shuffle=(), compress=3] = Array(c[n])
        end
    end
end

"""
    read(f::Union{HDF5.File, HDF5.Group}, ::Type{Chains})

Read a chain object from the given HDF5 file our group.

"""
function Base.read(f::Union{HDF5.File, HDF5.Group}, ::Type{Chains})
    secs = keys(f)
    pns = []
    datas = []
    name_map = Dict()
    for s in secs
        ns = keys(f[s])
        name_map[Symbol(s)] = ns
        for n in ns
            push!(pns, n)
            push!(datas, read(f[s], n))
        end
    end

    nc, ns = size(datas[1])
    np = size(datas,1)

    a = zeros(nc, np, ns)

    for i in 1:np
        a[:,i,:] = datas[i]
    end

    Chains(a, pns, name_map)
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.06"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data1 = CSV.read(datadir("purchaseData.csv"), copycols=true)
	Text("Original value: $(data1.Name[1])")
end

# ╔═╡ e303495a-243f-11eb-1440-a1d6ef299e6a
begin
	data2 = data1
	data2.Name[1] = "EMILY"
	data1.Name[1]
end

# ╔═╡ e30380ae-243f-11eb-351f-058cb130306d
begin
	data3 = CSV.read(datadir("purchaseData.csv"), copycols=true)
	data4 = copy(data3)
	data4.Name[1] = "EMILY"
	data3.Name[1]
end

# ╔═╡ e30411fc-243f-11eb-25b0-a33492c9b3fb
begin
	data5 = DataFrame()
	data5.X = [[0,1],[100,101]]
	data6 = copy(data5)
	data6.X[1][1] = -1
	data5.X[1][1]
end

# ╔═╡ e30b3752-243f-11eb-2e8d-bb6dc37d5c08
begin
	data7 = DataFrame(X = [[0,1],[100,101]])
	data8 = deepcopy(data7)
	data8.X[1][1] = -1
	data7.X[1][1]
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.06"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═e303495a-243f-11eb-1440-a1d6ef299e6a
# ╠═e30380ae-243f-11eb-351f-058cb130306d
# ╠═e30411fc-243f-11eb-25b0-a33492c9b3fb
# ╠═e30b3752-243f-11eb-2e8d-bb6dc37d5c08
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
function main()
    ands::Vector{Bool} = []
    ors::Vector{Bool} = []
    xors::Vector{Bool} = []
    for a in [false, true]
        for b in [false, true]
            push!(ands, a & b)
            push!(ors, a | b)
            push!(xors, a ⊻ b)
        end
    end
    @assert(ands == [false, false, false, true])
    @assert(ors == [false, true, true, true])
    @assert(xors == [false, true, true, false])
    println("OK")
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
# This file was generated, do not modify it. # hide
mat = Matrix(boston)
mat[1:3, 1:3]

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
An OMOP Common Data Model (CDM) type.
"""
abstract type CDMType end


text:
# # Solutions to exercises

# ## Setup

# Inspect Julia version:

VERSION

# The following instantiates a package environment.

# The package environment has been created using **Julia 1.6** and may not
# instantiate properly for other Julia versions.

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# ## Exercise 1

# Write a function named `total` that adds the elements of its vector input.

# ### Solution

function total(v)
    sum = 0.0 # better is zero(v)
    for i in 1:length(v)
        sum = sum + v[i]
    end
    return sum
end

#-

total(1:10)

# The built-in function is called `sum`.

# ## Exercise 2

# Generate a 1000 random samples from the standard normal
# distribution. Create a second such sample, and add the two samples
# point-wise.  Compute the (sample) mean and variance of the combined
# samples. In the same plot, show a frequency-normalized histogram of
# the combined samples and a plot of the pdf for normal distribution
# with zero mean and variance `2`.

# ### Solution

using Distributions, Plots, Statistics

samples1 = randn(1000); # or rand(Normal(), 1000)
samples2 = randn(1000);

samples = samples1 .+ samples2;

mu = mean(samples)
var = std(samples)^2

@show mu var

#-

d = Normal(0, sqrt(2))
f(x) = pdf(d, x)

xs = -5:(0.1):5
ys = f.(xs);

using Plots

plt = histogram(samples, normalize=true)
plot!(xs, ys)
plt

# ## Exercise 3

# The following shows that named tuples share some behaviour with dictionaries:

t = (x = 1, y = "cat", z = 4.5)
keys(t)

#-

t[:y]

#-

# Write a function called `dict` that converts a named tuple to an
# actual dictionary. You can create an empty dictionary using `Dict()`.

# ### Solution

function dict(t)
    d = Dict()
    for k in keys(t)
        d[k] = t[k]
    end
    return d
end

#-

dict(t)

# Or alternatively:

dict(t) = Dict(k => t[k] for k in keys(t))