text:
using DrWatson
@quickactivate "StatReth"

# %%
using DataFrames
using CSV
using StatsBase
using Distributions
using StatsPlots

# %% 4.7 - 4.11
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))
# precis(d)
d.height

d2 = filter(row -> row.age >= 18, d)    # either
d2 = d[d.age .>= 18, :]                 # or

# %% 4.12, 4.13
plot(p -> pdf(Normal(178, 20), p), 100:250)
plot(p -> pdf(Uniform(0, 50), p), -10:60)

# %% 4.14
N = 10_000
d_μ = Normal(178, 20)
sample_mu = rand(d_μ, N)
d_σ = Uniform(0, 50)
sample_sigma = rand(d_σ, N)

prior_h = rand.(Normal.(sample_mu, sample_sigma))
density(prior_h)

# %% 4.15
d_μ_100 = Normal(178, 100)
sample_μ_100 = rand(d_μ_100, N)

sample_height_100 = rand.(Normal.(sample_μ_100, sample_sigma))
# or is
# mix = MixtureModel(Normal.(sample_μ_100, sample_sigma))
# sample_height_100 = rand(mix, N)
# more correct?
density(sample_height_100)

# %% 4.16, 4.17
μ_grid = range(150, 160, length = 100)
σ_grid = range(7, 9, length = 100)
post = [(; :μ => μ, :σ => σ) for μ in μ_grid, σ in σ_grid]  # this is `post`

ll = [loglikelihood(Normal(μ, σ), d2.height) for (μ, σ) in post]
prior = [logpdf(d_μ, μ) + logpdf(d_σ, σ) for (μ, σ) in post]
lposterior = ll .+ prior
lposterior .-= maximum(lposterior)
posterior = exp.(lposterior)

wireframe(μ_grid, σ_grid, posterior)

# %% 4.19 - 4.22
samples = sample(post, pweights(posterior), N) |> DataFrame

scatter(samples.μ, samples.σ, alpha = 0.05, ms = 5, xlabel = "μ", ylabel = "σ")

density(samples.μ)
density(samples.σ)

quantile(samples.μ, (0.1, 0.9))
quantile(samples.σ, (0.1, 0.9))

# %% 4.23 - 4.25
d3_height = d2.height[1:20]
μ_grid = range(140, 160, length = 200)
σ_grid = range(4, 20, length = 200)
post = [(μ = μ, σ = σ) for μ in μ_grid, σ in σ_grid]

ll = [loglikelihood(Normal(μ, σ), d3_height) for (μ, σ) in post]
prior = [logpdf(d_μ, μ) + logpdf(d_σ, σ) for (μ, σ) in post]
lposterior = ll .+ prior
lposterior .-= maximum(lposterior)
posterior = exp.(lposterior)

wireframe(μ_grid, σ_grid, posterior)

samples2 = sample(post, pweights(posterior), N) |> DataFrame

scatter(samples2.μ, samples2.σ, alpha = 0.05, ms = 5)
density(samples2.μ)
density(samples2.σ)


text:
using DifferentialEquations
using Distributions
using Distances
using Gadfly
using ApproxBayes
using DataFrames
using Plots

function getsolution(sol, times)
    x1 = map(x -> sol(x)[1], times)
    y1 = map(x -> sol(x)[2], times)
  return x1, y1
end

# define Lotka Voltera ODE where a and b can be modified
f = @ode_def LV begin
  dx = a*x - b*x*y
  dy = b*x*y - y
end a b

x0 = [1.0; 0.5]
tspan = (0.0, 15.0)
p = [2.0, 0.8]
prob = ODEProblem(f, x0, tspan, p)
sol = solve(prob)
Plots.plot(sol)

#generate target data by sampling 15 points and then adding Gaussian noise
times = 1.0:2.0:15.0
x, y = getsolution(sol, times)
x .+= rand(Normal(0.0, 1), length(times))
y .+= rand(Normal(0.0, 1), length(times))
targetdata = [x, y]

#simulations function for ABC. return distance (sum of squared distances) and solution
function simLV(params, constants, targetdata)
  a = params[1]
  b = params[2]
  x0 = [1.0; 0.5]
  tspan = (0.0, 15.0)
  times = 1.0:2.0:15.0
  prob = ODEProblem(f, x0, tspan, [a, b])
  sol = solve(prob)
  x1, y1 = getsolution(sol, times)
  d = sum((x1 .- targetdata[1]).^2 + (y1 .- targetdata[2]).^2)
  return d, sol
end

#define ABC setup type
setup = ABCSMC(simLV,
  2,
  0.1,
  Prior([Uniform(0.0, 5.0), Uniform(0.0, 5.0)]),
  maxiterations = 10^6,
  convergence = 0.001,
  nparticles = 1000
  )
  #run ABC SMC algorithm
@time ressmc = runabc(setup, targetdata, verbose = true, progress = true, parallel=true)

#show results
show(ressmc)

#plot posterior parameters
plotparameterposterior(ressmc)

#plots results using Gadfly
#Create data frame for target data
DFt = DataFrame(x = targetdata[1], y = targetdata[2], time = times)
d1 = stack(DFt, [:x, :y])

#get solution for
res = mean(ressmc.parameters, weights(ressmc.weights), 1)
x0 = [1.0; 0.5]
tspan = (0.0, 15.0)
prob = ODEProblem(f, x0, tspan, [res[1], res[2]])
sol = solve(prob)
times2 = 1.0:0.001:15.0
x2, y2 = getsolution(sol, times2)

DF1 = DataFrame(x = x2, y = y2, time = times2)
d2 = stack(DF1, [:x, :y])

l1 = layer(d1,x = :time, y = :value, color = :variable, Geom.point)
l2 = layer(d2,x = :time, y = :value, color = :variable, Geom.line)
Gadfly.plot(l1, l2)

ext: jl


text:
__precompile__(true)
module QuartzImageIO
#import Base: error, size
using Images, Colors, ColorVectorSpace, FixedPointNumbers, Compat
import FileIO: @format_str, File, Stream, filename, stream

# We need to export writemime_, since that's how ImageMagick does it.
export writemime_

typealias CFURLRef Ptr{Void}
typealias CFStringRef Ptr{UInt8}
typealias CFDictionaryRef Ptr{Void}
typealias CGImageDestinationRef Ptr{Void}
typealias CGImageRef Ptr{Void}
typealias CGColorSpaceRef Ptr{Void}
typealias CGContextRef Ptr{Void}

image_formats = [
    format"BMP",
    format"GIF",
    format"JPEG",
    format"PNG",
    format"TIFF",
    format"TGA",
]

# There's a way to get the mapping through
# UTTypeCreatePreferredIdentifierForTag, but a dict is less trouble for now
const apple_format_names = Dict(format"BMP" => "com.microsoft.bmp",
                                format"GIF" => "com.compuserve.gif",
                                format"JPEG" => "public.jpeg",
                                format"PNG" => "public.png",
                                format"TIFF" => "public.tiff",
                                format"TGA" => "com.truevision.tga-image")

# The rehash! is necessary because of a precompilation issue
function __init__() Base.rehash!(apple_format_names) end

get_apple_format_name(format) = apple_format_names[format]

for format in image_formats
    eval(quote
        load(image::File{$format}, args...; key_args...) = load_(filename(image), args...; key_args...)
        load(io::Stream{$format}, args...; key_args...) = load_(read(io), args...; key_args...)
        save(fname::File{$format}, img::Image, args...; key_args...) =
            save_(filename(fname), img, get_apple_format_name($format), args...;
                  key_args...)
        save(io::Stream{$format}, img::Image, args...; key_args...) =
            save_(stream(io), img, get_apple_format_name($format), args...;
                  key_args...)
    end)
end

function load_(b::Array{UInt8, 1})
    data = CFDataCreate(b)
    imgsrc = CGImageSourceCreateWithData(data)
    CFRelease(data)
    read_and_release_imgsrc(imgsrc)
end

function load_(filename)
    myURL = CFURLCreateWithFileSystemPath(abspath(filename))
    imgsrc = CGImageSourceCreateWithURL(myURL)
    CFRelease(myURL)
    read_and_release_imgsrc(imgsrc)
end

## core, internal function
function read_and_release_imgsrc(imgsrc)
    if imgsrc == C_NULL
        warn("OSX reader created no image source")
        return nothing
    end
    # Get image information
    imframes = convert(Int, CGImageSourceGetCount(imgsrc))
    if imframes == 0
        # Bail out to ImageMagick
        warn("OSX reader found no frames")
        CFRelease(imgsrc)
        return nothing
    end
    dict = CGImageSourceCopyPropertiesAtIndex(imgsrc, 0)
    imheight = CFNumberGetValue(CFDictionaryGetValue(dict, "PixelHeight"), Int16)
    imwidth = CFNumberGetValue(CFDictionaryGetValue(dict, "PixelWidth"), Int16)
    isindexed = CFBooleanGetValue(CFDictionaryGetValue(dict, "IsIndexed"))
    if isindexed
        # Bail out to ImageMagick
        warn("OSX reader: indexed color images not implemented")
        CFRelease(imgsrc)
        return nothing
    end
    hasalpha = CFBooleanGetValue(CFDictionaryGetValue(dict, "HasAlpha"))

    pixeldepth = CFNumberGetValue(CFDictionaryGetValue(dict, "Depth"), Int16)
    # Colormodel is one of: "RGB", "Gray", "CMYK", "Lab"
    colormodel = CFStringGetCString(CFDictionaryGetValue(dict, "ColorModel"))
    if colormodel == ""
        # Bail out to ImageMagick
        warn("OSX reader found empty colormodel string")
        CFRelease(imgsrc)
        return nothing
    end
    imtype = CFStringGetCString(CGImageSourceGetType(imgsrc))
    alphacode, storagedepth = alpha_and_depth(imgsrc)

    # Get image description string
    imagedescription = ""
    if imtype == "public.tiff"
        tiffdict = CFDictionaryGetValue(dict, "{TIFF}")
        imagedescription = tiffdict != C_NULL ?
            CFStringGetCString(CFDictionaryGetValue(tiffdict, "ImageDescription")) : nothing
    end
    CFRelease(dict)

    # Allocate the buffer and get the pixel data
    sz = imframes > 1 ? (convert(Int, imwidth), convert(Int, imheight), convert(Int, imframes)) : (convert(Int, imwidth), convert(Int, imheight))
    const ufixedtype = Dict(10=>UFixed10, 12=>UFixed12, 14=>UFixed14, 16=>UFixed16)
    T = pixeldepth <= 8 ? UFixed8 : ufixedtype[pixeldepth]
    if colormodel == "Gray" && alphacode == 0 && storagedepth == 1
        buf = Array(Gray{T}, sz)
        fillgray!(reinterpret(T, buf, tuple(sz...)), imgsrc)
    elseif colormodel == "Gray" && in(alphacode, [1, 3])
        buf = Array(GrayA{T}, sz)
        fillgrayalpha!(reinterpret(T, buf, tuple(2, sz...)), imgsrc)
    elseif colormodel == "Gray" && in(alphacode, [2, 4])
        buf = Array(AGray{T}, sz)
        fillgrayalpha!(reinterpret(T, buf, tuple(2, sz...)), imgsrc)
    elseif colormodel == "RGB" && in(alphacode, [1, 3])
        buf = Array(RGBA{T}, sz)
        fillcolor!(reinterpret(T, buf, tuple(4, sz...)), imgsrc, storagedepth)
    elseif colormodel == "RGB" && in(alphacode, [2, 4])
        buf = Array(ARGB{T}, sz)
        fillcolor!(reinterpret(T, buf, tuple(4, sz...)), imgsrc, storagedepth)
    elseif colormodel == "RGB" && alphacode == 0
        buf = Array(RGB{T}, sz)
        fillcolor!(reinterpret(T, buf, tuple(3, sz...)), imgsrc, storagedepth)
    elseif colormodel == "RGB" && in(alphacode, [5, 6])
        buf = alphacode == 5 ? Array(RGB4{T}, sz) : Array(RGB1{T}, sz)
        fillcolor!(reinterpret(T, buf, tuple(4, sz...)), imgsrc, storagedepth)
    else
        warn("Unknown colormodel ($colormodel) and alphacode ($alphacode) found by OSX reader")
        CFRelease(imgsrc)
        return nothing
    end
    CFRelease(imgsrc)

    # Set the image properties
    prop = Dict(
        "spatialorder" => ["x", "y"],
        "pixelspacing" => [1, 1],
        "colorspace" => colormodel,
        "imagedescription" => imagedescription,
        "suppress" => Set(Any["imagedescription"])
    )
    if imframes > 1
        prop["timedim"] = ndims(buf)
    end
    Image(buf, prop)
end

function alpha_and_depth(imgsrc)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)  # Check only first frame
    alphacode = CGImageGetAlphaInfo(CGimg)
    bitspercomponent = CGImageGetBitsPerComponent(CGimg)
    bitsperpixel = CGImageGetBitsPerPixel(CGimg)
    CGImageRelease(CGimg)
    # Alpha codes documented here:
    # https://developer.apple.com/library/mac/documentation/graphicsimaging/reference/CGImage/Reference/reference.html#//apple_ref/doc/uid/TP30000956-CH3g-459700
    # Dividing bits per pixel by bits per component tells us how many
    # color + alpha slices we have in the file.
    alphacode, convert(Int, div(bitsperpixel, bitspercomponent))
end

function fillgray!{T}(buffer::AbstractArray{T, 2}, imgsrc)
    imwidth, imheight = size(buffer, 1), size(buffer, 2)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)
    imagepixels = CopyImagePixels(CGimg)
    pixelptr = CFDataGetBytePtr(imagepixels, eltype(buffer))
    imbuffer = unsafe_wrap(Array, pixelptr, (imwidth, imheight), false)
    buffer[:, :] = imbuffer
    CFRelease(imagepixels)
    CGImageRelease(CGimg)
end

# Image stack
function fillgray!{T}(buffer::AbstractArray{T, 3}, imgsrc)
    imwidth, imheight, nimages = size(buffer, 1), size(buffer, 2), size(buffer, 3)
    for i in 1:nimages
        CGimg = CGImageSourceCreateImageAtIndex(imgsrc, i - 1)
        imagepixels = CopyImagePixels(CGimg)
        pixelptr = CFDataGetBytePtr(imagepixels, T)
        imbuffer = unsafe_wrap(Array, pixelptr, (imwidth, imheight), false)
        buffer[:, :, i] = imbuffer
        CFRelease(imagepixels)
        CGImageRelease(CGimg)
    end
end

function fillgrayalpha!(buffer::AbstractArray{UInt8, 3}, imgsrc)
    imwidth, imheight = size(buffer, 2), size(buffer, 3)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)
    imagepixels = CopyImagePixels(CGimg)
    pixelptr = CFDataGetBytePtr(imagepixels, UInt16)
    imbuffer = unsafe_wrap(Array, pixelptr, (imwidth, imheight), false)
    buffer[1, :, :] = imbuffer & 0xff
    buffer[2, :, :] = div(imbuffer & 0xff00, 256)
    CFRelease(imagepixels)
    CGImageRelease(CGimg)
end
fillgrayalpha!(buffer::AbstractArray{UFixed8, 3}, imgsrc) = fillgrayalpha!(reinterpret(UInt8, buffer), imgsrc)

function fillcolor!{T}(buffer::AbstractArray{T, 3}, imgsrc, nc)
    imwidth, imheight = size(buffer, 2), size(buffer, 3)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)
    imagepixels = CopyImagePixels(CGimg)
    pixelptr = CFDataGetBytePtr(imagepixels, T)
    imbuffer = unsafe_wrap(Array, pixelptr, (nc, imwidth, imheight), false)
    buffer[:, :, :] = imbuffer
    CFRelease(imagepixels)
    CGImageRelease(CGimg)
end

function fillcolor!{T}(buffer::AbstractArray{T, 4}, imgsrc, nc)
    imwidth, imheight, nimages = size(buffer, 2), size(buffer, 3), size(buffer, 4)
    for i in 1:nimages
        CGimg = CGImageSourceCreateImageAtIndex(imgsrc, i - 1)
        imagepixels = CopyImagePixels(CGimg)
        pixelptr = CFDataGetBytePtr(imagepixels, T)
        imbuffer = unsafe_wrap(Array, pixelptr, (nc, imwidth, imheight), false)
        buffer[:, :, :, i] = imbuffer
        CFRelease(imagepixels)
        CGImageRelease(CGimg)
    end
end

## Saving Images ###############################################################

function save_and_release(cg_img::Ptr{Void}, # CGImageRef
                          fname, image_type::AbstractString)
    out_url = CFURLCreateWithFileSystemPath(fname)
    out_dest = CGImageDestinationCreateWithURL(out_url, image_type, 1)
    CGImageDestinationAddImage(out_dest, cg_img)
    CGImageDestinationFinalize(out_dest)
    CFRelease(out_dest)
    CFRelease(out_url)
    CGImageRelease(cg_img)
    nothing
end

""" `save_(fname, img::Image, image_type)`

- fname is the name of the file to save to
- image_type should be one of Apple's image types (eg. "public.jpeg")
"""
function save_(fname, img::AbstractImage, image_type)
    # TODO:
    # - avoid this convert call where possible
    # - support writing greyscale images
    # - spatialorder? It seems to work already, maybe because of convert.
    img2 = convert(Image{RGBA{UFixed8}}, img)
    buf = reinterpret(FixedPointNumbers.UInt8, Images.data(img2))
    nx, ny = size(img2)
    colspace = CGColorSpaceCreateDeviceRGB()
    bmp_context = CGBitmapContextCreate(buf, nx, ny, 8, nx*4, colspace,
                                        kCGImageAlphaPremultipliedLast)
    CFRelease(colspace)
    cgImage = CGBitmapContextCreateImage(bmp_context)
    CFRelease(bmp_context)

    save_and_release(cgImage, fname, image_type)
end

function save_(io::IO, img::AbstractImage, image_type)
    write(io, getblob(img, image_type))
end

function getblob(img::AbstractImage, format)
    # In theory we could save the image directly to a buffer via
    # CGImageDestinationCreateWithData - TODO. But I couldn't figure out how
    # to get the length of the CFMutableData object. So I take the inefficient
    # route of saving the image to a temporary file for now.
    @assert format == "png" || format == "public.png" # others not supported for now
    temp_file = "/tmp/QuartzImageIO_temp.png"
    save_(temp_file, img, "public.png")
    read(open(temp_file))
end

@deprecate writemime_(io::IO, ::MIME"image/png", img::AbstractImage) save(Stream(format"PNG", io), img)

## OSX Framework Wrappers ######################################################

# Commented out functions remain here because they might be useful for future
# debugging.

const foundation = Libdl.find_library(["/System/Library/Frameworks/Foundation.framework/Resources/BridgeSupport/Foundation"])
const imageio = Libdl.find_library(["/System/Library/Frameworks/ImageIO.framework/ImageIO"])

const kCFNumberSInt8Type = 1
const kCFNumberSInt16Type = 2
const kCFNumberSInt32Type = 3
const kCFNumberSInt64Type = 4
const kCFNumberFloat32Type = 5
const kCFNumberFloat64Type = 6
const kCFNumberCharType = 7
const kCFNumberShortType = 8
const kCFNumberIntType = 9
const kCFNumberLongType = 10
const kCFNumberLongLongType = 11
const kCFNumberFloatType = 12
const kCFNumberDoubleType = 13
const kCFNumberCFIndexType = 14
const kCFNumberNSIntegerType = 15
const kCFNumberCGFloatType = 16
const kCFNumberMaxType = 16

# enum defined at https://developer.apple.com/library/mac/documentation/GraphicsImaging/Reference/CGImage/index.html#//apple_ref/c/tdef/CGImageAlphaInfo
const kCGImageAlphaNone = 0
const kCGImageAlphaPremultipliedLast = 1
const kCGImageAlphaPremultipliedFirst = 2
const kCGImageAlphaLast = 3
const kCGImageAlphaFirst = 4
const kCGImageAlphaNoneSkipLast = 5
const kCGImageAlphaNoneSkipFirst = 6
const kCGImageAlphaOnly = 7

# Objective-C and NS wrappers
oms{T}(id, uid, ::Type{T}=Ptr{Void}) =
    ccall(:objc_msgSend, T, (Ptr{Void}, Ptr{Void}), id, selector(uid))

ogc{T}(id, ::Type{T}=Ptr{Void}) =
    ccall((:objc_getClass, "Cocoa.framework/Cocoa"), Ptr{Void}, (Ptr{UInt8}, ), id)

selector(sel::AbstractString) = ccall(:sel_getUid, Ptr{Void}, (Ptr{UInt8}, ), sel)

NSString(init::AbstractString) = ccall(:objc_msgSend, Ptr{Void},
                               (Ptr{Void}, Ptr{Void}, Ptr{UInt8}, UInt64),
                               oms(ogc("NSString"), "alloc"),
                               selector("initWithCString:encoding:"), init, 4)

# NSLog(str::AbstractString, obj) = ccall((:NSLog, foundation), Ptr{Void},
#                                 (Ptr{Void}, Ptr{Void}), NSString(str), obj)

# NSLog(str::AbstractString) = ccall((:NSLog, foundation), Ptr{Void},
#                            (Ptr{Void}, ), NSString(str))

# NSLog(obj::Ptr) = ccall((:NSLog, foundation), Ptr{Void}, (Ptr{Void}, ), obj)

# Core Foundation
# General

# CFRetain(CFTypeRef::Ptr{Void}) = CFTypeRef != C_NULL &&
#     ccall(:CFRetain, Void, (Ptr{Void}, ), CFTypeRef)

CFRelease(CFTypeRef::Ptr{Void}) = CFTypeRef != C_NULL &&
    ccall(:CFRelease, Void, (Ptr{Void}, ), CFTypeRef)

# function CFGetRetainCount(CFTypeRef::Ptr{Void})
#     CFTypeRef == C_NULL && return 0
#     ccall(:CFGetRetainCount, Clonglong, (Ptr{Void}, ), CFTypeRef)
# end

CFShow(CFTypeRef::Ptr{Void}) = CFTypeRef != C_NULL &&
    ccall(:CFShow, Void, (Ptr{Void}, ), CFTypeRef)

# function CFCopyDescription(CFTypeRef::Ptr{Void})
#     CFTypeRef == C_NULL && return C_NULL
#     ccall(:CFCopyDescription, Ptr{Void}, (Ptr{Void}, ), CFTypeRef)
# end

# CFCopyTypeIDDescription(CFTypeID::Cint) = CFTypeRef != C_NULL &&
#     ccall(:CFCopyTypeIDDescription, Ptr{Void}, (Cint, ), CFTypeID)

# function CFGetTypeID(CFTypeRef::Ptr{Void})
#     CFTypeRef == C_NULL && return nothing
#     ccall(:CFGetTypeID, Culonglong, (Ptr{Void}, ), CFTypeRef)
# end

# CFURLCreateWithString(filename) =
#     ccall(:CFURLCreateWithString, Ptr{Void},
#           (Ptr{Void}, Ptr{Void}, Ptr{Void}), C_NULL, NSString(filename), C_NULL)

CFURLCreateWithFileSystemPath(filename::AbstractString) =
    ccall(:CFURLCreateWithFileSystemPath, Ptr{Void},
          (Ptr{Void}, Ptr{Void}, Cint, Bool), C_NULL, NSString(filename), 0, false)

# CFDictionary

# CFDictionaryGetKeysAndValues(CFDictionaryRef::Ptr{Void}, keys, values) =
#     CFDictionaryRef != C_NULL &&
#     ccall(:CFDictionaryGetKeysAndValues, Void,
#           (Ptr{Void}, Ptr{Ptr{Void}}, Ptr{Ptr{Void}}), CFDictionaryRef, keys, values)

function CFDictionaryGetValue(CFDictionaryRef::Ptr{Void}, key)
    CFDictionaryRef == C_NULL && return C_NULL
    ccall(:CFDictionaryGetValue, Ptr{Void},
          (Ptr{Void}, Ptr{Void}), CFDictionaryRef, key)
end

CFDictionaryGetValue(CFDictionaryRef::Ptr{Void}, key::AbstractString) =
    CFDictionaryGetValue(CFDictionaryRef::Ptr{Void}, NSString(key))

# CFNumber
function CFNumberGetValue(CFNum::Ptr{Void}, numtype)
    CFNum == C_NULL && return nothing
    out = Cint[0]
    ccall(:CFNumberGetValue, Bool, (Ptr{Void}, Cint, Ptr{Cint}), CFNum, numtype, out)
    out[1]
end

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{Int8}) =
    CFNumberGetValue(CFNum, kCFNumberSInt8Type)

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{Int16}) =
    CFNumberGetValue(CFNum, kCFNumberSInt16Type)

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{Int32}) =
    CFNumberGetValue(CFNum, kCFNumberSInt32Type)

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{Int64}) =
    CFNumberGetValue(CFNum, kCFNumberSInt64Type)

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{Float32}) =
    CFNumberGetValue(CFNum, kCFNumberFloat32Type)

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{Float64}) =
    CFNumberGetValue(CFNum, kCFNumberFloat64Type)

CFNumberGetValue(CFNum::Ptr{Void}, ::Type{UInt8}) =
    CFNumberGetValue(CFNum, kCFNumberCharType)

#CFBoolean
CFBooleanGetValue(CFBoolean::Ptr{Void}) =
    CFBoolean != C_NULL &&
    ccall(:CFBooleanGetValue, Bool, (Ptr{Void}, ), CFBoolean)

# CFString
function CFStringGetCString(CFStringRef::Ptr{Void})
    CFStringRef == C_NULL && return ""
    buffer = Array(UInt8, 1024)  # does this need to be bigger for Open Microscopy TIFFs?
    res = ccall(:CFStringGetCString, Bool, (Ptr{Void}, Ptr{UInt8}, UInt, UInt16),
                CFStringRef, buffer, length(buffer), 0x0600)
    res == C_NULL && return ""
    return unsafe_string(pointer(buffer))
end

# These were unsafe, can return null pointers at random times.
# See Apple Developer Docs
# CFStringGetCStringPtr(CFStringRef::Ptr{Void}) =
#     ccall(:CFStringGetCStringPtr, Ptr{UInt8}, (Ptr{Void}, UInt16), CFStringRef, 0x0600)
#
# getCFString(CFStr::Ptr{Void}) = CFStringGetCStringPtr(CFStr) != C_NULL ?
#     bytestring(CFStringGetCStringPtr(CFStr)) : ""

# Core Graphics
# CGImageSource

CGImageSourceCreateWithURL(myURL::Ptr{Void}) =
    ccall((:CGImageSourceCreateWithURL, imageio), Ptr{Void}, (Ptr{Void}, Ptr{Void}), myURL, C_NULL)

CGImageSourceCreateWithData(data::Ptr{Void}) =
    ccall((:CGImageSourceCreateWithData, imageio), Ptr{Void}, (Ptr{Void}, Ptr{Void}), data, C_NULL)

CGImageSourceGetType(CGImageSourceRef::Ptr{Void}) =
    ccall(:CGImageSourceGetType, Ptr{Void}, (Ptr{Void}, ), CGImageSourceRef)

CGImageSourceGetStatus(CGImageSourceRef::Ptr{Void}) =
    ccall(:CGImageSourceGetStatus, UInt32, (Ptr{Void}, ), CGImageSourceRef)

CGImageSourceGetStatusAtIndex(CGImageSourceRef::Ptr{Void}, n) =
    ccall(:CGImageSourceGetStatusAtIndex, Int32,
          (Ptr{Void}, Csize_t), CGImageSourceRef, n) #Int32?

# CGImageSourceCopyProperties(CGImageSourceRef::Ptr{Void}) =
#     ccall(:CGImageSourceCopyProperties, Ptr{Void},
#           (Ptr{Void}, Ptr{Void}), CGImageSourceRef, C_NULL)

CGImageSourceCopyPropertiesAtIndex(CGImageSourceRef::Ptr{Void}, n) =
    ccall(:CGImageSourceCopyPropertiesAtIndex, Ptr{Void},
          (Ptr{Void}, Csize_t, Ptr{Void}), CGImageSourceRef, n, C_NULL)

CGImageSourceGetCount(CGImageSourceRef::Ptr{Void}) =
    ccall(:CGImageSourceGetCount, Csize_t, (Ptr{Void}, ), CGImageSourceRef)

CGImageSourceCreateImageAtIndex(CGImageSourceRef::Ptr{Void}, i) =
    ccall(:CGImageSourceCreateImageAtIndex, Ptr{Void},
          (Ptr{Void}, UInt64, Ptr{Void}), CGImageSourceRef, i, C_NULL)

# CGImageGet

CGImageGetAlphaInfo(CGImageRef::Ptr{Void}) =
    ccall(:CGImageGetAlphaInfo, UInt32, (Ptr{Void}, ), CGImageRef)

# Use this to detect if image contains floating point values
# CGImageGetBitmapInfo(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetBitmapInfo, UInt32, (Ptr{Void}, ), CGImageRef)

CGImageGetBitsPerComponent(CGImageRef::Ptr{Void}) =
    ccall(:CGImageGetBitsPerComponent, Csize_t, (Ptr{Void}, ), CGImageRef)

CGImageGetBitsPerPixel(CGImageRef::Ptr{Void}) =
    ccall(:CGImageGetBitsPerPixel, Csize_t, (Ptr{Void}, ), CGImageRef)

# CGImageGetBytesPerRow(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetBytesPerRow, Csize_t, (Ptr{Void}, ), CGImageRef)

CGImageGetColorSpace(CGImageRef::Ptr{Void}) =
    ccall(:CGImageGetColorSpace, UInt32, (Ptr{Void}, ), CGImageRef)

# CGImageGetDecode(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetDecode, Ptr{Float64}, (Ptr{Void}, ), CGImageRef)

# CGImageGetHeight(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetHeight, Csize_t, (Ptr{Void}, ), CGImageRef)

# CGImageGetRenderingIntent(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetRenderingIntent, UInt32, (Ptr{Void}, ), CGImageRef)

# CGImageGetShouldInterpolate(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetShouldInterpolate, Bool, (Ptr{Void}, ), CGImageRef)

# CGImageGetTypeID() =
#     ccall(:CGImageGetTypeID, Culonglong, (),)

# CGImageGetWidth(CGImageRef::Ptr{Void}) =
#     ccall(:CGImageGetWidth, Csize_t, (Ptr{Void}, ), CGImageRef)

CGImageRelease(CGImageRef::Ptr{Void}) =
    ccall(:CGImageRelease, Void, (Ptr{Void}, ), CGImageRef)

# Get pixel data
# See: https://developer.apple.com/library/mac/qa/qa1509/_index.html
CGImageGetDataProvider(CGImageRef::Ptr{Void}) =
    ccall(:CGImageGetDataProvider, Ptr{Void}, (Ptr{Void}, ), CGImageRef)

CGDataProviderCopyData(CGDataProviderRef::Ptr{Void}) =
    ccall(:CGDataProviderCopyData, Ptr{Void}, (Ptr{Void}, ), CGDataProviderRef)

CopyImagePixels(inImage::Ptr{Void}) =
    CGDataProviderCopyData(CGImageGetDataProvider(inImage))

CFDataGetBytePtr{T}(CFDataRef::Ptr{Void}, ::Type{T}) =
    ccall(:CFDataGetBytePtr, Ptr{T}, (Ptr{Void}, ), CFDataRef)

# CFDataGetLength(CFDataRef::Ptr{Void}) =
#     ccall(:CFDataGetLength, Ptr{Int64}, (Ptr{Void}, ), CFDataRef)

CFDataCreate(bytes::Array{UInt8,1}) =
    ccall(:CFDataCreate,Ptr{Void},(Ptr{Void},Ptr{UInt8},Csize_t),C_NULL,bytes,length(bytes))

### For output #################################################################

""" `check_null(x)`

Triggers an error if `x` is `NULL`, else returns `x` """
function check_null(x)
    # Poor-man's error handling. TODO: raise more specific exceptions.
    if x == C_NULL
        error("C call returned NULL")
    else
        x
    end
end

function CGImageDestinationCreateWithURL(url::CFURLRef,
                                         filetype::AbstractString,
                                         count::Integer,
                                         options::CFDictionaryRef=C_NULL)
    check_null(ccall((:CGImageDestinationCreateWithURL, imageio),
                     CGImageDestinationRef,
                     (CFURLRef, CFStringRef, Csize_t, CFDictionaryRef),
                     url, NSString(filetype), count, options))
end

function CGImageDestinationAddImage(dest::CGImageDestinationRef,
                                    image::CGImageRef,
                                    properties::CFDictionaryRef=C_NULL)
    # Returns NULL.
    # From the Apple docs: "The function logs an error if you add more images
    # than what you specified when you created the image destination. "
    # Maybe we should catch that somehow?
    ccall((:CGImageDestinationAddImage, imageio),
          Ptr{Void},
          (CGImageDestinationRef, CGImageRef, CFDictionaryRef),
          dest, image, properties)
end

type WritingImageFailed <: Exception end
function CGImageDestinationFinalize(dest::CGImageDestinationRef)
    rval = ccall((:CGImageDestinationFinalize, imageio),
                 Bool, (CGImageDestinationRef,), dest)
    # See https://developer.apple.com/library/mac/documentation/GraphicsImaging/Reference/CGImageDestination/index.html#//apple_ref/c/func/CGImageDestinationFinalize
    if !rval throw(WritingImageFailed()) end
end

CGColorSpaceCreateDeviceRGB() =
    check_null(ccall((:CGColorSpaceCreateDeviceRGB, imageio),
                     CGColorSpaceRef, ()))

function CGBitmapContextCreate(data, # void*
                               width, height, # size_t
                               bitsPerComponent, bytesPerRow, # size_t
                               space::CGColorSpaceRef,
                               bitmapInfo) # uint32_t
    check_null(ccall((:CGBitmapContextCreate, imageio),
                     CGContextRef,
                     (Ptr{Void}, Csize_t, Csize_t, Csize_t, Csize_t,
                      CGColorSpaceRef, UInt32),
                     data, width, height, bitsPerComponent, bytesPerRow, space,
                     bitmapInfo))
end

function CGBitmapContextCreateImage(context_ref::CGContextRef)
    check_null(ccall((:CGBitmapContextCreateImage, imageio),
                     CGImageRef,
                     (CGContextRef,), context_ref))
end

end # Module


text:
module FactorGraphTest

using Test
using ForneyLab: FactorGraph, FactorNode, Interface, Edge, Variable, generateId, currentGraph, setCurrentGraph,
generateId, addNode!, hasNode, addVariable!, hasVariable, Clamp

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function MockNode(; id=generateId(MockNode))
        self = new(id, Interface[], Dict{Symbol,Interface}())

        return self
    end
end

@testset "FactorGraph" begin
    # FactorGraph should initiate a new factor graph
    g = FactorGraph()

    @test isa(g.nodes, Dict{Symbol, FactorNode})
    @test isa(g.edges, Vector{Edge})
    @test isa(g.variables, Dict{Symbol, Variable})
    @test isa(g.counters, Dict{String, Int})
    @test isa(g.placeholders, Dict{Clamp, Tuple{Symbol, Int}})

    # currentGraph() should point to the current graph
    @test ===(currentGraph(), g)
    f = FactorGraph()
    @test ===(currentGraph(), f)
    setCurrentGraph(g)
    @test ===(currentGraph(), g)
end

@testset "generateId" begin
    g = FactorGraph()

    # generateId should generate a unique id based on type
    @test generateId(MockNode) == :mocknode_1
    @test generateId(MockNode) == :mocknode_2
    @test g.counters["mocknode"] == 2
end

@testset "addNode!" begin
    g = FactorGraph()
    nd = MockNode(id=:nd)

    # addNode! should register a node with a graph
    @test isempty(g.nodes)
    addNode!(g, nd)
    @test ===(g.nodes[:nd], nd)
    @test hasNode(g, nd)

    # addNode should not add the same node twice
    @test_throws Exception addNode!(g, nd)
end

end


text:
using TypeClasses
using Test
using Distributed

wait_a_little(f::Function, seconds=0.3) = @spawnat :any begin
  sleep(seconds)
  f()
end

wait_a_little(x, seconds=0.3) = wait_a_little(() -> x, seconds)

# Functor, Applicative, Monad
# ===========================

squared = map(wait_a_little(4)) do x
  x*x
end;  # returns a Task

@test fetch(squared) == 16

@test typeof(mapn(+, wait_a_little(11), wait_a_little(12))) === Future
@test fetch(mapn(+, wait_a_little(11), wait_a_little(12))) == 23

monadic = @syntax_flatmap begin
  a = wait_a_little(5)
  b = wait_a_little(a + 3)
  @pure a, b
end;  # returns a Task

@test fetch(monadic) == (5, 8)

@test fetch(pure(Future, "a")) == "a"

# Monoid, Alternative
# ===================

@test fetch(wait_a_little("hello.") ⊕ wait_a_little("world.")) == "hello.world."
@test fetch(wait_a_little(:a, 1.0) ⊘ wait_a_little(:b, 2.0)) == :a
@test fetch(wait_a_little(:a, 3.0) ⊘ wait_a_little(:b, 2.0)) == :b

@test fetch(wait_a_little(() -> error("fails"), 0.1) ⊘ wait_a_little(:succeeds, 0.3)) == :succeeds

@test fetch(wait_a_little(:succeeds, 0.3) ⊘ wait_a_little(() -> error("fails"), 0.1)) == :succeeds

@test isa(fetch(wait_a_little(() -> error("fails1")) ⊘ wait_a_little(() -> error("fails2")) ⊘ wait_a_little(() -> error("fails3")) ⊘ wait_a_little(() -> error("fails4"))),
   RemoteException)


text:
using ExamplePackage
using Test
using LinearAlgebra

@testset "ExamplePackage.jl" begin
    @test printDocstrDict(LinearAlgebra, :eigen) == IdDict{Any, Any}(Union{Tuple{TB}, Tuple{TA}, Tuple{AbstractMatrix{TA}, AbstractMatrix{TB}}} where {TA, TB} => Base.Docs.DocStr(svec("    eigen(A, B) -> GeneralizedEigen\n\nComputes the generalized eigenvalue decomposition of `A` and `B`, returning a\n[`GeneralizedEigen`](@ref) factorization object `F` which contains the generalized eigenvalues in\n`F.values` and the generalized eigenvectors in the columns of the matrix `F.vectors`.\n(The `k`th generalized eigenvector can be obtained from the slice `F.vectors[:, k]`.)\n\nIterating the decomposition produces the components `F.values` and `F.vectors`.\n\nAny keyword arguments passed to `eigen` are passed through to the lower-level\n[`eigen!`](@ref) function.\n\n# Examples\n```jldoctest\njulia> A = [1 0; 0 -1]\n2×2 Matrix{Int64}:\n 1   0\n 0  -1\n\njulia> B = [0 1; 1 0]\n2×2 Matrix{Int64}:\n 0  1\n 1  0\n\njulia> F = eigen(A, B);\n\njulia> F.values\n2-element Vector{ComplexF64}:\n 0.0 - 1.0im\n 0.0 + 1.0im\n\njulia> F.vectors\n2×2 Matrix{ComplexF64}:\n  0.0+1.0im   0.0-1.0im\n -1.0+0.0im  -1.0-0.0im\n\njulia> vals, vecs = F; # destructuring via iteration\n\njulia> vals == F.values && vecs == F.vectors\ntrue\n```\n"), nothing, Dict{Symbol, Any}(:typesig => Union{Tuple{TB}, Tuple{TA}, Tuple{AbstractMatrix{TA}, AbstractMatrix{TB}}} where {TA, TB}, :module => LinearAlgebra, :linenumber => 456, :binding => LinearAlgebra.eigen, :path => "/buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl")), Tuple{Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}} where {T<:Real, S}, UnitRange} => Base.Docs.DocStr(svec("    eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> Eigen\n\nComputes the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`\nwhich contains the eigenvalues in `F.values` and the eigenvectors in the columns of the\nmatrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)\n\nIterating the decomposition produces the components `F.values` and `F.vectors`.\n\nThe following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).\n\nThe [`UnitRange`](@ref) `irange` specifies indices of the sorted eigenvalues to search for.\n\n!!! note\n    If `irange` is not `1:n`, where `n` is the dimension of `A`, then the returned factorization\n    will be a *truncated* factorization.\n"), nothing, Dict{Symbol, Any}(:typesig => Tuple{Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}} where {T<:Real, S}, UnitRange}, :module => LinearAlgebra, :linenumber => 685, :binding => LinearAlgebra.eigen, :path => "/buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl")), Union{Tuple{AbstractMatrix{T}}, Tuple{T}} where T => Base.Docs.DocStr(svec("    eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> Eigen\n\nComputes the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`\nwhich contains the eigenvalues in `F.values` and the eigenvectors in the columns of the\nmatrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)\n\nIterating the decomposition produces the components `F.values` and `F.vectors`.\n\nThe following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).\n\nFor general nonsymmetric matrices it is possible to specify how the matrix is balanced\nbefore the eigenvector calculation. The option `permute=true` permutes the matrix to become\ncloser to upper triangular, and `scale=true` scales the matrix by its diagonal elements to\nmake rows and columns more equal in norm. The default is `true` for both options.\n\nBy default, the eigenvalues and vectors are sorted lexicographically by `(real(λ),imag(λ))`.\nA different comparison function `by(λ)` can be passed to `sortby`, or you can pass\n`sortby=nothing` to leave the eigenvalues in an arbitrary order.   Some special matrix types\n(e.g. [`Diagonal`](@ref) or [`SymTridiagonal`](@ref)) may implement their own sorting convention and not\naccept a `sortby` keyword.\n\n# Examples\n```jldoctest\njulia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])\nEigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\nvalues:\n3-element Vector{Float64}:\n  1.0\n  3.0\n 18.0\nvectors:\n3×3 Matrix{Float64}:\n 1.0  0.0  0.0\n 0.0  1.0  0.0\n 0.0  0.0  1.0\n\njulia> F.values\n3-element Vector{Float64}:\n  1.0\n  3.0\n 18.0\n\njulia> F.vectors\n3×3 Matrix{Float64}:\n 1.0  0.0  0.0\n 0.0  1.0  0.0\n 0.0  0.0  1.0\n\njulia> vals, vecs = F; # destructuring via iteration\n\njulia> vals == F.values && vecs == F.vectors\ntrue\n```\n"), nothing, Dict{Symbol, Any}(:typesig => Union{Tuple{AbstractMatrix{T}}, Tuple{T}} where T, :module => LinearAlgebra, :linenumber => 179, :binding => LinearAlgebra.eigen, :path => "/buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl")), Tuple{Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}} where {T<:Real, S}, Real, Real} => Base.Docs.DocStr(svec("    eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> Eigen\n\nComputes the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`\nwhich contains the eigenvalues in `F.values` and the eigenvectors in the columns of the\nmatrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)\n\nIterating the decomposition produces the components `F.values` and `F.vectors`.\n\nThe following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).\n\n`vl` is the lower bound of the window of eigenvalues to search for, and `vu` is the upper bound.\n\n!!! note\n    If [`vl`, `vu`] does not contain all eigenvalues of `A`, then the returned factorization\n    will be a *truncated* factorization.\n"), nothing, Dict{Symbol, Any}(:typesig => Tuple{Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}} where {T<:Real, S}, Real, Real}, :module => LinearAlgebra, :linenumber => 711, :binding => LinearAlgebra.eigen, :path => "/buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl")))
    @test printDocstrDict(LinearAlgebra, :triu) == IdDict{Any, Any}(Tuple{AbstractMatrix{T} where T, Integer} => Base.Docs.DocStr(svec("    triu(M, k::Integer)\n\nReturns the upper triangle of `M` starting from the `k`th superdiagonal.\n\n# Examples\n```jldoctest\njulia> a = fill(1.0, (4,4))\n4×4 Matrix{Float64}:\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n\njulia> triu(a,3)\n4×4 Matrix{Float64}:\n 0.0  0.0  0.0  1.0\n 0.0  0.0  0.0  0.0\n 0.0  0.0  0.0  0.0\n 0.0  0.0  0.0  0.0\n\njulia> triu(a,-3)\n4×4 Matrix{Float64}:\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n```\n"), nothing, Dict{Symbol, Any}(:typesig => Tuple{AbstractMatrix{T} where T, Integer}, :module => LinearAlgebra, :linenumber => 367, :binding => LinearAlgebra.triu, :path => "/buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl")), Tuple{AbstractMatrix{T} where T} => Base.Docs.DocStr(svec("    triu(M)\n\nUpper triangle of a matrix.\n\n# Examples\n```jldoctest\njulia> a = fill(1.0, (4,4))\n4×4 Matrix{Float64}:\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n 1.0  1.0  1.0  1.0\n\njulia> triu(a)\n4×4 Matrix{Float64}:\n 1.0  1.0  1.0  1.0\n 0.0  1.0  1.0  1.0\n 0.0  0.0  1.0  1.0\n 0.0  0.0  0.0  1.0\n```\n"), nothing, Dict{Symbol, Any}(:typesig => Tuple{AbstractMatrix{T} where T}, :module => LinearAlgebra, :linenumber => 319, :binding => LinearAlgebra.triu, :path => "/buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl")))
end

max_stars_repo_name: AnirudhPraveen/ExamplePackage

max_issues_repo_name: AnirudhPraveen/ExamplePackage

max_forks_repo_name: AnirudhPraveen/ExamplePackage


text:
"Goal specification for backward search."
struct BackwardSearchGoal{G <: Goal,C} <: Goal
    goal::G # Original goal specification
    start::State # Start state to be reached via backward search
    constraint_diff::C # State constraints as a diff, if any
end

BackwardSearchGoal(goal::Goal, start::State) =
    BackwardSearchGoal(goal, start, nothing)
BackwardSearchGoal(goal::StateConstrainedGoal, start::State) =
    BackwardSearchGoal(goal, start, precond_diff(constraints))

Base.hash(spec::BackwardSearchGoal, h::UInt) =
    hash(spec.start, hash(spec.goal, h))
Base.isequal(s1::BackwardSearchGoal, s2::BackwardSearchGoal) =
    s1.state == s2.state && s1.goal == s2.goal

is_goal(spec::BackwardSearchGoal, domain::Domain, state::State) =
    issubset(state, spec.start)
is_violated(spec::BackwardSearchGoal, domain::Domain, state::State) =
    is_violated(spec.goal, domain, state)
get_cost(spec::BackwardSearchGoal, domain::Domain, s1::State, a::Term, s2::State) =
    get_cost(spec.goal, domain, s2, a, s1)
get_reward(spec::BackwardSearchGoal, domain::Domain, s1::State, a::Term, s2::State) =
    get_reward(spec.goal, domain, s2, a, s1)
get_goal_terms(spec::BackwardSearchGoal) =
    get_goal_terms(spec.goal)

add_constraints!(spec::BackwardSearchGoal, state::State) =
    nothing
add_constraints!(spec::BackwardSearchGoal{G,PDDL.Diff}, state::State) where {G} =
    update!(state, spec.constraint_diff)


text:
function foo()
    lines = readlines(open("data/2.txt", "r"))
    for i in 1:26
        x = Set{String}()
        for line in lines
            if i == 1
                str = line[2:end]
            elseif i == 26
                str = line[1:end-1]
            else
                str = line[1:i-1] * line[i+1:end]
            end
            if str in x
                return str
            end
            push!(x, str)
        end
    end
end

println(foo())


text:
using Test, LinearAlgebra, DynamicPolynomials
using LearningAlgebraicVarieties
const LAV = LearningAlgebraicVarieties

### auxiliary_functions
@testset "Auxiliary Functions" begin
    data = rand(2,10)

    D = LAV.ScaledEuclidean(data)
    @test maximum(D) <= 1

    D = LAV.FubiniStudyDistances(data)
    @test size(D) == (10,10)

    D = zeros(10,10)
    LAV.FubiniStudyDistances!(D, data)
    @test maximum(D) > 0

    D = LAV.ScaledFubiniStudy(data)
    @test maximum(D) <= 1

    @polyvar x y
    f = x^2 + 2y^2 - 1

    D = LAV.EllipsoidDistances(data, f, 0.5)
    @test size(D) == (10,10)
    D = LAV.EllipsoidDistances(data, [f], 0.5)
    @test size(D) == (10,10)
end

### estimate_dimension
@testset "PCA" begin
    data_affine = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)
    data_linear = hcat(map(i -> [i[1];-i[1]], rand(10))...)

    d = LAV.EstimateDimensionPCA(data_affine, false)
    @test d == 1
    d = LAV.EstimateDimensionPCA(data_linear, false)
    @test d == 1

    d = LAV.EstimateDimensionPCA(data_affine, true)
    @test d == 1
    d = LAV.EstimateDimensionPCA(data_linear, true)
    @test d == 0

    ϵ_array = Array{Float64}(range(0.1, length = 25, stop = 0.9))

    d = LAV.EstimateDimensionNPCA(data_affine, ϵ_array, false)
    @test d[end] == 1
    d = LAV.EstimateDimensionNPCA(data_affine, false)
    @test d[end] == 1
    d = LAV.EstimateDimensionNPCA(data_linear, ϵ_array, false)
    @test d[end] == 1

    d = LAV.EstimateDimensionNPCA(data_affine, ϵ_array, true)
    @test d[end] == 1
    d = LAV.EstimateDimensionNPCA(data_affine, true)
    @test d[end] == 1
    d = LAV.EstimateDimensionNPCA(data_linear, ϵ_array, true)
    @test d[end] == 0
end

@testset "CorrSum" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)
    ϵ_array = Array{Float64}(range(0.1, length = 25, stop = 0.9))

    d = LAV.EstimateDimensionCorrSum(data, ϵ_array, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionCorrSum(data, ϵ_array, true)
    @test length(d) == 25

    d = LAV.EstimateDimensionCorrSum(data, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionCorrSum(data, true)
    @test length(d) == 25
end

@testset "BoxCounting" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)
    ϵ_array = Array{Float64}(range(0.1, length = 25, stop = 0.9))

    d = LAV.EstimateDimensionBoxCounting(data, ϵ_array, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionBoxCounting(data, ϵ_array, true)
    @test length(d) == 25

    d = LAV.EstimateDimensionBoxCounting(data, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionBoxCounting(data, true)
    @test length(d) == 25
end

@testset "Bickel-Levina" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)
    ϵ_array = Array{Float64}(range(0.1, length = 25, stop = 0.9))

    d = LAV.EstimateDimensionMLE(data, ϵ_array, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionMLE(data, ϵ_array, true)
    @test length(d) == 25

    d = LAV.EstimateDimensionMLE(data, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionMLE(data, true)
    @test length(d) == 25
end

@testset "ANOVA" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)
    ϵ_array = Array{Float64}(range(0.1, length = 25, stop = 0.9))

    d = LAV.EstimateDimensionANOVA(data, ϵ_array, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionANOVA(data, ϵ_array, true)
    @test length(d) == 25

    d = LAV.EstimateDimensionANOVA(data, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionANOVA(data, true)
    @test length(d) == 25
end

@testset "PHCurve" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)
    ϵ_array = Array{Float64}(range(0.1, length = 25, stop = 0.9))

    d = LAV.EstimateDimensionPHCurve(data, ϵ_array, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionPHCurve(data, ϵ_array, true)
    @test length(d) == 25

    d = LAV.EstimateDimensionPHCurve(data, false)
    @test length(d) == 25
    d = LAV.EstimateDimensionPHCurve(data, true)
    @test length(d) == 25
end

###estimate_equations
@testset "FindEquations" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)

    f = LAV.FindEquations(data, :with_svd, 1, false)
    @test abs(f[1]([1;0])) < 1e-15

    f = LAV.FindEquations(data, :with_svd, 2, true)
    @test f[1] == 0
end

@testset "Vandermonde" begin
    data = hcat(map(i -> [i[1];-i[1]+1], rand(10))...)

    M = LAV.MultivariateVandermondeMatrix(data, [[1,0], [1,1]])
    @test size(M.Vandermonde) == (10,2)

    M = LAV.MultivariateVandermondeMatrix(data, 1, false)
    @test size(M.Vandermonde) == (10,3)

    f = LAV.FindEquations(M, :with_svd)
    @test abs(f[1]([1;0])) < 1e-15

    f = LAV.FindEquations(M, :with_svd, 1e-12)
    @test abs(f[1]([1;0])) < 1e-15

    f = LAV.FindEquations(M, :with_qr)
    @test abs(f[1]([1;0])) < 1e-15

    f = LAV.FindEquations(M, :with_qr, 1e-12)
    @test abs(f[1]([1;0])) < 1e-15

    f = LAV.FindEquations(M, :with_rref)
    @test abs(f[1]([1;0])) < 1e-15

    f = LAV.FindEquations(M, :with_rref, 1e-12)
    @test abs(f[1]([1;0])) < 1e-15
end


text:
function formdataseries(;
  refreshcomp::Bool = true,
  refreshtvolmcapcrsp::Bool = true,
  datapath::String=DATA_PATH,
  levname::String = LEV_NAME,
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM,
  binextension::String = BIN_EXTENSION,
  outpath::String = OUT_PATH,
  refreshdataseries::Bool = true,
  trialrunonly::Bool = false,
  retaineddataseriescolumns = RETAINED_COLS_DATA_SERIES)::DataFrame

  local lev::DataFrame
  local fyear::Int

  local dailyunlesstrial::Symbol = trialrunonly ? :monthly : :daily
  refreshdataseries = (refreshdataseries || refreshcomp)

  levname = "$levname-$dailyunlesstrial"

  if refreshdataseries
    lev = prepcomp(refreshcomp=refreshcomp)

    #supplements the market value and constructs the trailing daily vol
    lev = maketvolmcapseries!(
      lev, refreshtvolmcapcrsp=refreshtvolmcapcrsp, crsptype=dailyunlesstrial)

    #lev = marketretseries!(lev, refreshtretcrsp = refreshtretcrsp, crsptype=:monthly)

    issorted(lev, [:lpermno, :fdate]) || error("lev should be sorted here")
    CSV.write("output\\$(levname)_ts.csv", lev[1:20_000, :])

    winsorizekeydata!(lev)
    constructportfolios!(lev)
    doubleincreasegroups!(lev)

    issorted(lev, [:lpermno, :fdate]) || sort!(lev, [:lpermno, :fdate])
    if TEST_OUTPUT
      fyear = 2015
      CSV.write("output\\$(levname)_$fyear.csv", lev[(y->y==fyear).(lev.fyear), :])

      fyear = 2014
      CSV.write("output\\$(levname)_$fyear.csv", lev[(y->y==fyear).(lev.fyear), :])
    end

    select!(lev, RETAINED_COLS_DATA_SERIES)

    outbinstream("$datapath\\$levname.$binextension", lev)

  else

    lev = inbinstream("$datapath\\$levname.$binextension")

  end

  return lev
end

function winsorizekeydata!(lev::DataFrame;
    fields2winsorize::AbstractArray{Symbol}=[CQUART_FVALS; CQUART_FVALS_D; CQUART_FVALS_BRS_D],
    winsorizeprop=WINSORIZE_PROP)

  #dups will lead to race errors
  @assert allunique(fields2winsorize)

  @mpar for s ∈ fields2winsorize
    local nbefore::Int = length(unique(lev[!, s]))
    winsorize!(lev[!, s], prop = winsorizeprop)
    (length(unique(lev[!, s])) == nbefore) && (@warn "winsorization didn't do anything for field $s.
      Something is probably wrong.")
  end

  return nothing
end

function maketvolmcapseries!(lev::DataFrame;
  refreshtvolmcapcrsp::Bool = true,
  crsptype::Symbol = :daily,
  months2stale::Int = MONTHS_2_STALE_SHORT,
  requiredannualdatacols::Vector{Symbol} = REQUIRED_ANNUAL_DATA_COLS,
  lagmonthrange::UnitRange = LAG_MONTH_RANGE)

  crspcolumns::Vector{Symbol} = [:date, :permno, :mktcap,
    :vol252dnet, :vol252dret, :ret252d, :lret, :ret]
  periods2stale = Month(months2stale)

  validpermnos::Vector{Int} = unique(lev.lpermno)

  print("time to load vol: ")
  @time crsp::DataFrame = prepcrsp(
    crspcolumns = crspcolumns,
    crsptype = crsptype,
    refreshcrsp=refreshtvolmcapcrsp,
    validpermnos = validpermnos,
    prefix = "tvolmcap") do df::DataFrame
    computecrspmktcap!(df)
    computenetreturns!(df)
    #winsorize!(df.net, prop=0.999)

    #it will be useful to flag end of months (or end of periods if cut short)
    issorted(df, [:permno, :date]) || error("crsp must be sorted")
    df.eom = Vector{Bool}(undef, size(df,1))
    sdfs::GroupedDataFrame = groupby(df, :permno)
    @mpar for i ∈ 1:length(sdfs)
      sdf::SubDataFrame = sdfs[i]
      months::Vector{Int} = (month).(sdf.date)
      sdf.eom .=  [months[1:(end-1)] .≠ months[2:end]; true]
    end

    #compute the vols
    computetrailingvol!(df, crsptype==:daily ? 252 : 12, [:net, :ret],
      Fvols = [:vol252dnet, :vol252dret],
      minpointsperwindow = crsptype==:daily ? REQUIRED_DAILY_TRAILING : REQUIRED_MONTHLY_TRAILING,
      eomonly = true)

    #compute the trailing rets
    computetrailingreturns!(df, crsptype==:daily ? 252 : 12,
      Flret = :lret, Fnewret = :ret252d,
      minpointsperwindow = crsptype==:daily ? REQUIRED_DAILY_TRAILING : REQUIRED_MONTHLY_TRAILING,
      eomonly = true)

    #winsorize per Ivo
    df.vol252dnet .= (f::MFloat64->winsorizelevel(f, 0., 0.1)).(df.vol252dnet)
    df.vol252dret .= (f::MFloat64->winsorizelevel(f, 0., 0.1)).(df.vol252dret)

  end

  lev = mergecrsp(crsp, lev,
    refreshmerge=true,
    cachedata = false,
    crsptype=crsptype,
    prefix="tvolmcap")

  obliterate!(crsp)

  #supplement the market value
  sort!(lev, [:permno, :date])
  lev.last = falses(size(lev,1))

  #gets the row in a dataframe given a date or the msot recent past row
  @inline function rowfromdatesorted(df::AbstractDataFrame, dt::Date)::Union{Nothing, DataFrameRow}
    local r::Union{DataFrameRow,Nothing}

    i::Int = searchsortedfirst(df.date, dt)

    #check if the index is valid for our purposes
    r = ((i>size(df,1)) #this means the date or a larger date wasn't found
      ? nothing : df[i,:])

    return r
  end

  slevs::GroupedDataFrame = groupby(lev, :permno)
  mkequityadded::Int = sum((ismissing).(lev.mkequity))
  @mpar for i ∈ 1:length(slevs)
    slev::SubDataFrame = slevs[i]

    for sslev ∈ groupby(slev, :compid)
      sslev.last[end] = true

      #now supplement the market equity
      local Nmissing::Int = sum((ismissing).(sslev.mkequity))
      if Nmissing > 1
        (Nmissing < size(sslev,1)) && error("Inconsistent presence of mkequity")

        dt::Date = eom(sslev.fyrmonth[1])

        #find the first date on or before the end of the fiscal year if available
        local r::Union{Nothing, DataFrameRow} = rowfromdatesorted(slev, dt)

        #note- need to adjust for different scaling between crsp and compustat
        #NOTE: In PreprocessCRSP we scaled mktcap to dollars, now we scale it down to millions
        (!isnothing(r)) && ( dt - periods2stale ≤ r.date) && (sslev.mkequity .= r.mktcap/1_000_000)
      end
    end
  end

  mkequityadded -= sum((ismissing).(lev.mkequity))
  @info "Supplemented $mkequityadded mkequity records using $crsptype data"

  #compute the returns and standard deviation over the fiscal year
  lev.fret = Vector{MFloat64}(missing, size(lev,1))
  lev.fvol = Vector{MFloat64}(missing, size(lev,1))
  periods2stalelong = Month(MONTHS_2_STALE_LONG)
  slevs = groupby(lev, [:permno]) #refresh the grouped df, not sure if this is necessary
  #=@mpar =#for i ∈ 1:length(slevs)
    slev::SubDataFrame = slevs[i]
    sslevs::GroupedDataFrame = groupby(slev, :compid)
    for j ∈ 1:length(sslevs)
      sslev = sslevs[j]
      fdate::Date = sslev.fdate[end]
      Lfdate::MDate = sslev.Lfdate[end]
      if ismissing(Lfdate)
        #as a check, if the lag is missing look for fdates within the acceptable window for lags
        if (j≠1) && (sslev[1,:fdate] - periods2stalelong ≤ sslevs[j-1][1, :fdate])
          println(view(slev, slev.last, [:permno, :date, :fdate, :Lfdate, :fyear]))
          error("unexpected missing lagged fdate")
        end
        Lfdate = sslev[1,:date] - Day(1)
      end

      #make sure the fdate interval is acceptable
      ((Lfdate) ≥ firstdayofmonth(fdate) - Month(18)) || continue
      ((Lfdate) ≤ lastdayofmonth(fdate) - Month(6)) || continue

      #ismissing(Lfdate) && continue
      s3lev::SubDataFrame = view(slev, (r::DataFrameRow->
        (r.date ≤ fdate) && (r.date > Lfdate)).(eachrow(slev)), :)
      sslev.fret .= exp(sum(skipmissing(s3lev.lret))) - 1.0
      sslev.fvol .= std(skipmissing(s3lev.ret))
    end
  end

  lev[1:20_000,:] |> CSV.write("output\\fretvollevcheck.csv")

  #now that we have patched the market value, drop all other rows
  #can either free the memory or do this in place
  @info "rows before dropping all but last in retvol: $(size(lev,1))"
  filter!(r->r.last, lev)
  @info "rows after dropping all but last in retvol: $(size(lev,1))"

  #improve mkequity quality, and make sure we have the field
  lev.mkequity .= (floororx).(lev.mkequity)
  lev = lev[(m::MFloat64->(!ismissing(m)) && (m>MIN_MKEQUITY)).(lev.mkequity),:]

  #make the cash field
  lev.mkat = lev.at .- lev.bkequity .+ lev.mkequity #old: lev.lt .+ lev.mkequity
  lev.mkat .= (floororx).(lev.mkat)
  lev.mknegcash = 1.0 .- (lev.cash ./ lev.mkat)
  lev.mknegcash .= (boundunit).(lev.mknegcash)

  #compute liab
  lev.mkliab = lev.netliab ./ lev.mkat #netliab is compa.lt .- compa.cash from preprocesscomp
  lev.mkliab .= (boundunit).(lev.mkliab)

  #compute #flev field
  lev.mkflev = lev.netdebt ./ (lev.netdebt .+ lev.mkequity)
  lev.mkflev .= (boundunit).(lev.mkflev)

  #make finite (note we winsorize later on)
  @mpar for s ∈ [:mknegcash, :mkflev, :mkliab]
    lev[!,s] .= (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(lev[!,s])
  end

  #sort!(lev, [:permno, :fdate]) #not completely sure we need this sort

  lev.bm = lev.bkequity ./ lev.mkequity
  finiteormissing!(lev, :bm)
  lev.bm .= ((x::MFloat64) -> winsorizelevel(x, 0.0, 5.0)).(lev.bm)

  #make lagged change in assets as a control variable
  #println(describe(lev))
  #println("Busy work: $(mean(rand(5000,5000)))")
  lev.lat = (log).(lev.at)
  finiteormissing!(lev, :lat)

  #maybe I should lag the book stuff earlier? maybe not
  issorted(lev, [:permno, :fdate]) || error("lev should be sorted on permno, fdate")
  differencewithin2!(lev, [:mknegcash, :bknegcash, :mkflev,
    :bkflev, :mkliab, :bkliab, :lat, :vol252dnet, :vol252dret,:mkequity], :permno, date=:fdate, sorted=true)

  lev.Dat = (exp).(lev.Dlat)
  lev.lmkequity = (log).(lev.mkequity*1_000_000) #rescale

  finiteormissing!(lev, :lmkequity)
  lagwithin2!(lev, [:vol252dnet, :Dvol252dnet,:vol252dret, :Dvol252dret,
    :mkat, :mkequity, :lmkequity, :bm, :Dlat, :nop,
    :mknegcash, :bknegcash, :mkflev, :bkflev, :mkliab, :bkliab,
    :Dmknegcash, :Dbknegcash, :Dmkflev, :Dbkflev, :Dmkliab, :Dbkliab,
    :netdebt, :lt, :at, :netliab],
    :permno, date=:fdate, sorted=true)
  lagwithin2!(lev, [:Lvol252dnet, :Lvol252dret], :permno, date=:fdate, sorted=true)

  ############################ return-related procedures
  #this is remaining content from the old return buildup procedure
  issorted(lev, [:permno, :fdate]) || error("lev should be sorted on permno, fdate")
  lagwithin2!(lev, [:ret252d], :permno, date=:fdate, sorted=true)
  lagwithin2!(lev, [:Lret252d], :permno, date=:fdate, sorted=true)

  #println(describe(lev))
  ##now make the table 5 leverage versions net of stock returns

  #this is the stock return constribution of leverage
  lev.RLmkequity = lev.Lmkequity .* (1 .+ lev.fret)
  lev.RDmkflev = lev.Lnetdebt ./ (lev.Lnetdebt .+ lev.RLmkequity) .- lev.Lmkflev
  lev.RDmkliab = lev.Lnetliab ./ (lev.Llt .+ lev.RLmkequity) .- lev.Lmkliab

  #the residual is the managerial component
  lev.MDmkflev = lev.Dmkflev .- lev.RDmkflev
  lev.MDmkliab = lev.Dmkliab .- lev.RDmkliab

  issorted(lev, [:permno, :fdate]) || error("lev should be sorted on permno, fdate")
  #WARNING: FILTERS
  #now we only allow records in Ivo's range
  #WARNING this drops a lot of records. Implemented 2152020
  #unlike Ivo I drop both next and lagged (C-L and L-LL)
  filter!(lev) do r
    (!ismissing(r[:fyrmonth])) &&
    (!ismissing(r[:Lfyrmonth])) &&
    (!ismissing(r[:LLfyrmonth])) &&
    (r[:fyrmonth] - r[:Lfyrmonth] ∈ lagmonthrange) &&
    (r[:Lfyrmonth] - r[:LLfyrmonth] ∈ lagmonthrange)
  end

  dropmissing!(lev, requiredannualdatacols)

  ######END IVO FILTERS
  rename!(lev, [:date=>:dailyenddate, :compid=>:compidtvolmcap])
  select!(lev, Not([:permno, :last])) #dup of lpermno

  if TEST_OUTPUT
    d = Date(2015,12,31)
    CSV.write("output\\levuniv_$d.csv", lev[(cd->cd==d).(lev.dailyenddate), :])
  end

  if TEST_OUTPUT
    TEST_OUTPUT && CSV.write(
      "output\\volretmcap_$(crsptype).csv",
        lev[(end-20000):end, :])
  end

  return lev
end

#creates the control variables with appropriate lagging
function makepanelmonthly(panel::DataFrame;
  crsptype::Symbol = :daily,
  refreshtretcrsp::Bool = true)

  local crsptvol::DataFrame
  local crsptret::DataFrame
  local crsp::DataFrame

  #load the previous crsp fields
  crsptvol = prepcrsp(
    crspcolumns = nothing,
    crsptype = crsptype,
    refreshcrsp=false,
    validpermnos = nothing,
    prefix = "tvolmcap") do df::DataFrame
      @assert false #do not refresh the crsp file here
    end
  select!(crsptvol, [:permno, :date, :vol252dnet, :vol252dret,:ret252d])

  crsptret = prepcrsp(
    crspcolumns = nothing,
    crsptype = :monthly,
    refreshcrsp=false,
    validpermnos = nothing,
    prefix = "tret") do df::DataFrame
    computenetreturns!(df)
    end
  select!(crsptret, [:permno, :date, :ret, :lret])

  #use the dialy data to supplement the existing monthly data (new dataset is monthly)
  crsp = join(crsptvol, crsptret, on=[:permno, :date])
  obliterate!(crsptvol)
  obliterate!(crsptret)
  crsp[1:20_000,:] |> CSV.write("$OUT_PATH\\crspret_$(crsptype)_ts.csv")

  #this will make the lookups easier
  crsp.dyrmonth = (d->YearMonth(Year(d).value, Month(d).value)).(crsp.date)

  select!(panel, Not([:vol252dnet, :ret252d])) #drop these since we will work with the monthly versions
  panel = mergecrsp(crsp, panel, refreshmerge=true, cachedata = false,
    crsptype=crsptype, prefix="panel")

  #preallocate
  N::Int = size(panel,1)
  panel.L12ret252d = Vector{MFloat64}(undef, N)
  panel.L24ret252d = Vector{MFloat64}(undef, N)
  #panel.ayret12m = Vector{MFloat64}(undef, N)

  panel.L12vol252dret = Vector{MFloat64}(undef, N)
  panel.L24vol252dret= Vector{MFloat64}(undef, N)
  #panel.ayvol252dnet= Vector{MFloat64}(undef, N)

  panel.ayoffset = panel.dyrmonth .- panel.ayrmonth

  sort!(panel, [:permno, :dyrmonth])
  @info "beginning panel creation"
  spanels::GroupedDataFrame = groupby(panel, [:permno])
  #=@mpar =#for i ∈ 1:length(spanels)
    local spanel::SubDataFrame = spanels[i]
    local Nspanel::Int = size(spanel,1)
    local ptr::Int
    local target::YearMonth

    #this should be true always
    if (!allunique(spanel.dyrmonth))
      println(spanel[!, [:permno, :gvkey, :adate, :begindate, :enddate, :dyrmonth]])
      error("Dates for dyrmonth are not unique w/in permno")
    end

    for (j, r) ∈ enumerate(eachrow(spanel))

      if j > 12 #compute teh 1yr lagged controls
        target = r.dyrmonth - 12
        ptr = searchsortedfirst(spanel.dyrmonth, target)
        if spanel.dyrmonth[ptr] == target
          r.L12ret252d = spanel[ptr, :ret252d]
          r.L12vol252dret = spanel[ptr, :vol252dret]
        end
      end

      if j > 24 #compute teh 2yr lagged controls
        target = r.dyrmonth - 24
        ptr = searchsortedfirst(spanel.dyrmonth, target)
        if spanel.dyrmonth[ptr] == target
          r.L24ret252d = spanel[ptr, :ret252d]
          r.L24vol252dret = spanel[ptr, :vol252dret]
        end
      end

    end
  end

  #println(unique(panel.fyoffset))
  #error("stop here")

  #println(describe(panel))
  @info "panel creation complete"

  return panel
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "test_analysis" begin
## TZERO ##
# Examples from the Emami-Naeini & Van Dooren Paper
# Example 3
A = [0 1 0 0 0 0;
     0 0 1 0 0 0;
     0 0 0 0 0 0;
     0 0 0 0 1 0;
     0 0 0 0 0 1;
     0 0 0 0 0 0]
B = Matrix([0 0 1 0 0 0;
     0 0 0 0 0 1]')
C = [1 1 0 0 0 0;
     0 0 0 1 -1 0]
D = [1 0;
     1 0]

ex_3 = ss(A, B, C, D)
@test tzero(ex_3) ≈ [0.3411639019140099 + 1.161541399997252im,
                             0.3411639019140099 - 1.161541399997252im,
                             0.9999999999999999 + 0.0im,
                             -0.6823278038280199 + 0.0im]
# Example 4
A = [-0.129    0.0   0.396e-1  0.25e-1    0.191e-1;
     0.329e-2  0.0  -0.779e-4  0.122e-3  -0.621;
     0.718e-1  0.0  -0.1       0.887e-3  -0.385e1;
     0.411e-1  0.0   0.0      -0.822e-1   0.0;
     0.361e-3  0.0   0.35e-4   0.426e-4  -0.743e-1]
B = [0.0  0.139e-2;
     0.0  0.359e-4;
     0.0  -0.989e-2;
     0.249e-4  0.0;
     0.0  -0.534e-5]
C = [1 0 0 0 0;
     0 1 0 0 0]
D = zeros(2, 2)
ex_4 = ss(A, B, C, D)
@test tzero(ex_4) ≈ [-0.06467751189940692,-0.3680512036036696]

# Example 5
s = tf("s")
ex_5 = 1/s^15
@test tzero(ex_5) == Float64[]
@test tzero(ss(ex_5)) == Float64[]

# Example 6
A = [2 -1 0;
     0 0 0;
     -1 0 0]
B = [0; 0; 1]
C = [0 -1 0]
D = [0]
ex_6 = ss(A, B, C, D)
@test tzero(ex_6) == Float64[]

@test ss(A, [0 0 1]', C, D) == ex_6

# Example 7
ex_7 = ss(zeros(2, 2), [0;1], [-1 0], [0])
@test tzero(ex_7) == Float64[]

# Example 8
A = [-2 1 0 0 0 0;
     1 -2 1 0 1 -1;
     0 1 -2 1 0 0;
     0 0 1 -1 0 1;
     0 -1 0 0 0 0;
     0 1 0 -1 0 0]
B = [1; 0; 0; 0; 1; 0]
C = [0 0 0 1 0 0]
D = [0]
ex_8 = ss(A, B, C, D)
# TODO : there may be a way to improve the precision of this example.
@test tzero(ex_8) ≈ [-1.0, -1.0] atol=1e-7

# Example 9
ex_9 = (s - 20)/s^15
@test tzero(ex_9) ≈ [20.0]
@test tzero(ss(ex_9)) ≈ [20.0]

# Example 11
A = [-2 -6 3 -7 6;
     0 -5 4 -4 8;
     0 2 0 2 -2;
     0 6 -3 5 -6;
     0 -2 2 -2 5]
B = Matrix([-2 -8 -3 1 -8;
     7 -5 0 5 0]')
C = [0 -1 2 -1 -1;
     1 1 1 0 -1;
     0 3 -2 3 -1]
D = [0 0;
     0 0;
     0 0]
ex_11 = ss(A, B, C, D)
@test tzero(ex_11) ≈ [4.0, -3.0]

# Test for multiple zeros, siso tf
sys = s*(s + 1)*(s^2 + 1)*(s - 3)/((s + 1)*(s + 4)*(s - 4))
@test tzero(sys) ≈ [3.0, -1.0, im, -im, 0.0]

## POLE ##
@test pole(sys) ≈ [4.0, -4.0, -1.0]
@test pole([sys sys]) ≈ [4.0, -4.0, -1.0] # Issue #81
@test pole(ex_11) ≈ eigvals(ex_11.A)
@test pole([2/(s+1) 3/(s+2); 1/(s+1) 1/(s+1)]) ≈ [-1, -1, -2]

poles = [-3.383889568918823 + 0.000000000000000im
                            -2.199935841931115 + 0.000000000000000im
                            -0.624778101910111 + 1.343371895589931im
                            -0.624778101910111 - 1.343371895589931im
                            -0.083309192664918 + 0.487701968391972im
                            -0.083309192664918 - 0.487701968391972im]
approxin2(el,col) = any(el.≈col)
# Compares the computed poles with the expected poles
# TODO: Improve the test for testing equalifity of sets of complex numbers
# i.e. simplify and handle doubles.
@test all(approxin(p,poles) for p in pole(ex_8)) && all(approxin2(p,pole(ex_8)) for p in poles)

ex_12 = ss(-3, 2, 1, 2)
@test pole(ex_12) ≈ [-3]

ex_13 = ss([-1 1; 0 -1], [0; 1], [1 0], 0)
@test pole(ex_13) ≈ [-1, -1]

## ZPKDATA ##

H = [tf(0) tf([3, 0],[1, 1, 10]) ; tf([1, 1],[1, 5]) tf([2],[1, 6])]
G = ss(H)
sol_z = vecarray(ComplexF64, 2, 2, ComplexF64[], ComplexF64[0.0 + 0.0im],
        ComplexF64[-1.0 + 0.0im], ComplexF64[])
sol_p = vecarray(ComplexF64, 2, 2, ComplexF64[], ComplexF64[-0.5 - 3.1224989991991996im,
        -0.5 + 3.1224989991991996im],
        ComplexF64[-5.0 + 0.0im], ComplexF64[-6.0 + 0.0im])
sol_k = [0.0 3.0; 1.0 2.0]
z, p, k = zpkdata(H)

@test_array_vecs_eps z sol_z 2*eps(Float64)
@test_broken true == false # order of poles changed below, should probably be consistent
#@test_array_vecs_eps p sol_p 2*eps(Float64)
@test k == sol_k
z, p, k = zpkdata(G)
@test_array_vecs_eps z sol_z 10*eps(Float64)
@test_array_vecs_eps p sol_p 10*eps(Float64)
@test k == sol_k

## GAIN ## #Gain is confusing when referring to zpkdata. Test dcgain instead
@test [dcgain(H[1, 1]) dcgain(H[1, 2]); dcgain(H[2, 1]) dcgain(H[2, 2])] ≈ [0 0; 0.2 1/3]
@test [dcgain(G[1, 1]) dcgain(G[1, 2]); dcgain(G[2, 1]) dcgain(G[2, 2])] ≈ [0 0; 0.2 1/3]

## MARKOVPARAM ##
@test markovparam(G, 0) == [0.0 0.0; 1.0 0.0]
@test markovparam(G, 1) == [0.0 3.0; -4.0 2.0]
@test markovparam(G, 2) == [0.0 -3.0; 20.0 -12.0]

## DAMP ##
@test damp(sys)[1] ≈ [1.0, 4.0, 4.0]
@test damp(sys)[2] ≈ [1.0, -1.0, 1.0]

damp_output = damp(ex_11)
@test damp_output[1] ≈ [1.0, 1.0, 2.0, 2.0, 3.0]
# THe order of the poles in ±1 and ±2 may come out in different order
@test damp_output[2][1:2] ≈ [-1.0, 1.0] || damp_output[2][1:2] ≈ [1.0, -1.0]
@test damp_output[2][3:4] ≈ [-1.0, 1.0] || damp_output[2][3:4] ≈ [1.0, -1.0]
@test damp_output[2][5] ≈ -1.0

## DAMPREPORT ##
@test sprint(dampreport, sys) == (
     "|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
     "|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
     "+---------------+---------------+---------------+---------------+\n"*
     "|  -1.000e+00   |  1.000e+00    |  1.000e+00    |  1.000e+00    |\n"*
     "|  4.000e+00    |  -1.000e+00   |  4.000e+00    |  -2.500e-01   |\n"*
     "|  -4.000e+00   |  1.000e+00    |  4.000e+00    |  2.500e-01    |\n")

@test sprint(dampreport, 1/(s+1+2im)/(s+2+3im)) == (
     "|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
     "|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
     "+---------------+---------------+---------------+---------------+\n"*
     "|  -1.000e+00   |  4.472e-01    |  2.236e+00    |  1.000e+00    |\n"*
     "|  -2.000e+00 im|               |               |               |\n"*
     "|  -2.000e+00   |  5.547e-01    |  3.606e+00    |  5.000e-01    |\n"*
     "|  -3.000e+00 im|               |               |               |\n")

# Example 5.5 from http://www.control.lth.se/media/Education/EngineeringProgram/FRTN10/2017/e05_both.pdf
G = [1/(s+2) -1/(s+2); 1/(s+2) (s+1)/(s+2)]
@test_broken length(pole(G)) == 1
@test length(tzero(G)) == 1
@test_broken size(minreal(ss(G)).A) == (2,2)

end


text:
using FixYourWorkaround
using Test

package_name = "package"
dir = mktempdir()
toml_path = joinpath(dir, "Project.toml")

@testset "Exceptions" begin
    @testset "CompatNotFound" begin
        write(toml_path, "")

        @test_throws CompatNotFound package_compatible(package_name, "0.0"; toml_path=toml_path)
    end

    @testset "PackageNotInCompat" begin
        write(toml_path, "[compat]")

        @test_throws PackageNotInCompat package_compatible(package_name, "0.0"; toml_path=toml_path)
    end
end

@testset "Package outside of Versions" begin
    write(
        toml_path,
        """
        [compat]
        $package_name = "1"
        """
    )

    @test_throws VersionNotCompatible package_compatible(package_name, "0.0"; toml_path=toml_path)
end

@testset "Package in Versions" begin
    write(
        toml_path,
        """
        [compat]
        $package_name = "1"
        """
    )

    @test package_compatible(package_name, "1.1"; toml_path=toml_path)
end

@testset "Version::Int64" begin
    write(
        toml_path,
        """
        [compat]
        $package_name = "1"
        """
    )

    @test package_compatible(package_name, 1; toml_path=toml_path)
end

@testset "Version::Float64" begin
    write(
        toml_path,
        """
        [compat]
        $package_name = "1"
        """
    )

    @test package_compatible(package_name, 1.0; toml_path=toml_path)
end

@testset "Version::VersionNumber" begin
    write(
        toml_path,
        """
        [compat]
        $package_name = "1"
        """
    )

    @test package_compatible(package_name, v"1"; toml_path=toml_path)
end


text:
abstract type AbstractKalmanFilter <: AbstractFilter end

@with_kw struct KalmanFilter{AT,BT,CT,DT,R1T,R2T,R2DT,D0T,XT,RT} <: AbstractKalmanFilter
    A::AT
    B::BT
    C::CT
    D::DT
    R1::R1T
    R2::R2T
    R2d::R2DT
    d0::D0T
    x::XT
    R::RT
    t::Ref{Int} = Ref(1)
end

"""
KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1))
"""
function KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1))
    all(iszero, D) || throw(ArgumentError("Nonzero D matrix not supported yet"))
    KalmanFilter(A,B,C,D,R1,R2,MvNormal(R2), d0, Vector(d0.μ), Matrix(d0.Σ), Ref(1))
end

sample_state(kf::AbstractKalmanFilter) = rand(kf.d0)
sample_state(kf::AbstractKalmanFilter, x, u, t) = kf.A*x .+ kf.B*u .+ rand(MvNormal(kf.R1))
sample_measurement(kf::AbstractKalmanFilter, x, t) = kf.C*x .+ rand(MvNormal(kf.R2))
particletype(kf::AbstractKalmanFilter) = typeof(kf.x)
covtype(kf::AbstractKalmanFilter)      = typeof(kf.R)
state(kf::AbstractKalmanFilter)        = kf.x
covariance(kf::AbstractKalmanFilter)   = kf.R

function reset!(kf::AbstractKalmanFilter)
    kf.x .= Vector(kf.d0.μ)
    kf.R .= copy(Matrix(kf.d0.Σ))
    kf.t[] = 1
end

# UKF ==========================================================================

@with_kw struct UnscentedKalmanFilter{DT,MT,R1T,R2T,R2DT,D0T,VT,XT,RT} <: AbstractKalmanFilter
    dynamics::DT
    measurement::MT
    R1::R1T
    R2::R2T
    R2d::R2DT
    d0::D0T
    xs::Vector{VT}
    x::XT
    R::RT
    t::Ref{Int} = Ref(1)
end

"""
    UnscentedKalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1))
"""
function UnscentedKalmanFilter(dynamics,measurement,R1,R2,d0=MvNormal(R1))
    n = size(R1,1)
    p = size(R2,1)
    R1 = SMatrix{n,n}(R1)
    R2 = SMatrix{n,n}(R2)
    xs = sigmapoints(mean(d0), cov(d0))
    UnscentedKalmanFilter(dynamics,measurement,R1,R2,MvNormal(Matrix(R2)), d0, xs, Vector(d0.μ), Matrix(d0.Σ), Ref(1))
end

sample_state(kf::UnscentedKalmanFilter) = rand(kf.d0)
sample_state(kf::UnscentedKalmanFilter, x, u, t) = kf.dynamics(x,u) .+ rand(MvNormal(Matrix(kf.R1)))
sample_measurement(kf::UnscentedKalmanFilter, x, t) = kf.measurement(x) .+ rand(MvNormal(Matrix(kf.R2)))

# function transform_moments!(S,X,m,L)
#     X .-= mean(X) # Normalize the sample
#     for i in eachindex(X)
#         S[i] = m .+ L*X[i]
#     end
#     S
# end

function sigmapoints(m,Σ)
    n = max(length(m), size(Σ,1))
    xs = [@SVector zeros(n) for _ in 1:(2n+1)]
    sigmapoints!(xs,m,Σ)
end

function sigmapoints!(xs, m, Σ::AbstractMatrix)
    n = length(xs[1])
    X = sqrt(Symmetric(n*Σ))
    for i in 1:n
        xs[i] = X[:,i]
        xs[i+n] = -xs[i] .+ m
        xs[i] = xs[i] .+ m
    end
    xs[end] = m
    xs
end

function predict!(ukf::UnscentedKalmanFilter, u, t = index(ukf))
    @unpack dynamics,measurement,x,xs,R,R1 = ukf
    ns = length(xs)
    sigmapoints!(xs,x,R) # TODO: these are calculated in the update step
    r = copy(R)
    for i in eachindex(xs)
        xs[i] = dynamics(xs[i],u)
    end
    x .= mean(xs)
    # for i in eachindex(xs)
    #     d = xs[i]-x
    #     r .+= Symmetric(d*d')
    # end
    # R .= Symmetric(r)./n + R1
    R .= cov(xs) + R1
    ukf.t[] += 1
end

function correct!(ukf::UnscentedKalmanFilter, y, t = index(ukf))
    @unpack measurement,x,xs,R,R1,R2,R2d = ukf
    n = size(R1,1)
    p = size(R2,1)
    ns = length(xs)
    sigmapoints!(xs,x,R) # Update sigmapoints here since untransformed points required
    # S = @SMatrix zeros(p,p)
    C = @SMatrix zeros(n,p)
    # ym = @SVector zeros(p) # y mean
    ys = map(measurement, xs)
    # for i in eachindex(xs)
    #     ys  = measurement(xs[i])
    #     ym += ys
    # end
    # ym  = ym ./ ns
    ym = mean(ys)
    for i in eachindex(ys)
        d   = ys[i]-ym
        # S .+= Symmetric(d*d')
        C  += Symmetric((xs[i]-x)*d')
    end
    e   = y .- ym
    S = cov(ys) + R2
    # S   = S./ns + R2
    K   = (C./ns)/S
    x .+= K*e
    R  .= R - K*S*K'
    0. #logpdf(MvNormal(Matrix(S)), e) - 1/2*logdet(S) # TODO: this is not correct
end

# SigmaFilter ==========================================================================

struct SigmaFilter{DT,MT,MLT,D0T,DFT,VT} <: AbstractParticleFilter
    dynamics::DT
    measurement::MT
    measurement_likelihood::MLT
    w::Vector{Float64}
    we::Vector{Float64}
    initial_density::D0T
    dynamics_density::DFT
    x::Vector{VT}
    xprev::Vector{VT}
    xm::Vector{Float64}
    R::Matrix{Float64}
    t::Ref{Int}
end

"""
    SigmaFilter(dynamics,measurement,d0)
"""
function SigmaFilter(N,dynamics,measurement,measurement_likelihood,df,d0)
    @show n = length(d0)
    xs = rand(d0,N)
    xs = reinterpret(SVector{n,Float64}, xs)[:]
    SigmaFilter(dynamics,measurement,measurement_likelihood,
        zeros(length(xs)), # w
        zeros(length(xs)), # we
        d0, df,
        xs,                # x
        copy(xs),          # xprev
        Vector(d0.μ),      # xm
        Matrix(d0.Σ),      # R
        Ref(1))
end

sample_state(sf::SigmaFilter) = rand(sf.initial_density)
sample_state(sf::SigmaFilter, x, u, t) = sf.dynamics(x,u,t,true)
sample_measurement(sf::SigmaFilter, x, t) = sf.measurement(x,t,true)
num_particles(sf::SigmaFilter) = length(sf.x)
particles(sf::SigmaFilter) = sf.x
weights(sf::SigmaFilter) = sf.w
expweights(sf::SigmaFilter) = sf.we
state(sf::SigmaFilter) = sf
rng(sf::SigmaFilter) = Random.GLOBAL_RNG

function predict!(sf::SigmaFilter, u, t = index(sf))
    @unpack dynamics,measurement,x,xprev,xm,R,w,we = sf
    N = length(x)
    n = length(x[1])
    xm .= sum(x->x[1]*x[2], zip(x,we)) # Must update based on correction
    R .= 0
    for i in eachindex(x)
        δ = x[i]-xm
        R .+= (δ*δ')*we[i]
    end
    R ./= (1-sum(abs2,we))
    R .= (R + R')/2 + 1e-16I
    d = MvNormal(xm,R)
    noisevec = zeros(length(d))
    for i in eachindex(x)
        xi = SVector{n,Float64}(rand!(d, noisevec))
        x[i] = dynamics(xi, u, t, true)
    end
    w  .= -log(N)
    we .= 1/N
    xm .= mean(x)
    R  .= cov(x)
    sf.t[] += 1
end

Base.@propagate_inbounds propagate_particles!(sf::SigmaFilter, u, j::Vector{Int}, t::Int, noise=true) = propagate_particles!(sf, u, t, noise)

# Base.@propagate_inbounds function measurement_equation!(sf::SigmaFilter, y, t, w = weights(sf))
#     g = measurement_likelihood(sf)
#     any(ismissing.(y)) && return w
#     x = particles(sf)
#     @inbounds for i = 1:num_particles(sf)
#         w[i] += g(x[i],y,t)
#     end
#     w
# end


text:
include("_common.jl")
include("_calib_funcs.jl")

if nworkers() < 2
    addprocs(min(length(APPROACHES), CPU_CORES), exeflags="--project=..")
end

include_everywhere("_common.jl")
include_everywhere("_calib_funcs.jl")

@everywhere baseline_calib = joinpath(DATA_PATH, "baselines") * "/"
@everywhere calib_param_path = joinpath(DATA_PATH, "calib_params") * "/"
@everywhere calib_fig_path = joinpath(FIG_PATH, "calib_state") * "/"
mkpath(calib_fig_path)
mkpath(calib_param_path)

@everywhere function run_state_calibration(metric)
    approach, objfunc = metric

    wid = Distributed.myid()
    @info "Worker $wid : Calibrating $approach"

    calib_start = (1825*2)+1  # 10-year burn in
    sn, n_id = setup_network("$(baseline_calib)cotter_baseline_IHACRES_$approach.yml")
    func = nothing
    if !occursin("RMSE", approach)
        func = (obs, sim) -> abs(1.0 - objfunc(obs, sim))
    else
        func = (obs, sim) -> objfunc(obs, sim)
    end

    # Use baseline calibration parameters as the base
    target_idx = (1,2,4,5,6,7,8)
    thresholds = [0.0, 0.1, 0.9]
    # thresholds = [0.0, 0.16, 0.84]
    popsize = 64 * length(target_idx)^2
    state = :storage
    state_str = String(state)

    # Calibrate network using the BlackBoxOptim package
    # keyword arguments will be passed to the `bboptimize()` function
    # This will set node parameters to the optimal values found
    metric = (active_param_set, param_idxs, obs, sim) -> bin_metric(active_param_set[calib_start:CALIB_LN], param_idxs, obs[calib_start:CALIB_LN], sim[calib_start:CALIB_LN], func, nothing)
    result, optobj = burn_state_based_calibrate(sn, n_id, state, CALIB_CLIMATE, CALIB_FLOW, metric,
                                                target_idx, thresholds;
                                                burn_in=calib_start,
                                                log_transform=true,
                                                Method=:borg_moea,
                                                # ϵ=0.025,
                                                FitnessScheme=ParetoFitnessScheme{length(thresholds)}(is_minimizing=true),
                                                PopulationSize=popsize, MaxTime=TRAINING_TIME, TraceInterval=300)

    # Save params
    params = best_candidate(result)
    outfile = "$(calib_param_path)log_$(state_str)_$(approach)_10_year.txt"
    open(outfile, "w") do f
        print(f, join(params, ","))
    end

    Streamfall.save_calibration!(result, optobj, outfile*".opt")

    best_score = best_fitness(result)

    reset!(sn)
    node = sn[n_id]
    _, x0, __ = param_info(node; with_level=false)
    mod_params = update_partial(x0, target_idx, params)
    active_param_set, param_idxs = online_burn_state_node!(sn, n_id, FULL_CLIMATE, state, mod_params, thresholds; burn_in=calib_start, log_transform=true, releases=nothing)

    sim_flow = node.outflow
    @assert best_score == metric(active_param_set[1:CALIB_LN], param_idxs, CALIB_OBS[1:CALIB_LN], sim_flow[1:CALIB_LN])

    calib_data = FULL_DATASET[calib_start:CALIB_LN, "410730_Q"]
    calib_sim = sim_flow[calib_start:CALIB_LN]
    rmse = round(Streamfall.RMSE(calib_data, calib_sim), digits=2)
    nse = round(Streamfall.NSE(calib_data, calib_sim), digits=2)
    mkge = round(Streamfall.mKGE(calib_data, calib_sim), digits=2)
    plot(FULL_DATASET.Date[calib_start:CALIB_LN], calib_data,
         title="$approach\n(KGE': $mkge; NSE: $nse; RMSE: $rmse)",
         legend=:topright,
         label="Historic",
         xlabel="Date",
         ylabel="Streamflow [ML]")
    plot!(CALIB_DATES[calib_start:end], calib_sim, label="Calibration", alpha=0.7)
    savefig(joinpath(calib_fig_path, "log_$(state_str)_$(approach)_10_year_calib.png"))

    qqplot(calib_data, calib_sim, title="QQ Calibration\n($(state_str); $approach)")
    xlabel!("Historic")
    ylabel!("Simulated")
    savefig(joinpath(calib_fig_path, "log_$(state_str)_$(approach)_10_year_qq_calib.png"))

    valid_data = FULL_DATASET[CALIB_LN+1:end, "410730_Q"]
    valid_sim = sim_flow[CALIB_LN+1:end]
    rmse = round(Streamfall.RMSE(valid_data, valid_sim), digits=2)
    nse = round(Streamfall.NSE(valid_data, valid_sim), digits=2)
    mkge = round(Streamfall.mKGE(valid_data, valid_sim), digits=2)
    plot(FULL_DATASET.Date[CALIB_LN+1:end], valid_data,
         title="$approach\n(KGE': $mkge; NSE: $nse; RMSE: $rmse)",
         legend=:topleft,
         label="Historic",
         xlabel="Date",
         ylabel="Streamflow [ML]")
    plot!(VALID_DATES, valid_sim, label="Validation", alpha=0.7)
    savefig(joinpath(calib_fig_path, "log_$(state_str)_$(approach)_10_year_valid.png"))

    qqplot(valid_data, valid_sim, title="QQ Validation\n($(state_str); $approach)")
    xlabel!("Historic")
    ylabel!("Simulated")
    savefig(joinpath(calib_fig_path, "log_$(state_str)_$(approach)_10_year_qq_valid.png"))
end

pmap(run_state_calibration, zip(APPROACHES, OBJFUNCS))

@info "Finished"


text:
using Markdown

"""
    InteractiveLogger(stream::IO; min_level=Info))

Logger for logging to an interactive terminal, formatting logs in a readable
way.
"""
mutable struct InteractiveLogger <: AbstractLogger
    stream::IO
    default_min_level::LogLevel
    catch_exceptions::Bool
    prev_progress_key
    message_counts::Dict{Any,Int}
    module_limits::Dict{Module,LogLevel}
    blacklisted_ids::Set{Any}
end

function InteractiveLogger(stream::IO; min_level=Info, catch_exceptions=true)
    Base.depwarn("InteractiveLogger is deprecated - use ConsoleLogger instead.", :InteractiveLogger)
    InteractiveLogger(stream, min_level, catch_exceptions, nothing,
             Dict{Any,Int}(), Dict{Module,LogLevel}(), Set{Any}())
end

function configure_logging(logger::InteractiveLogger, _module=nothing;
                           min_level=Info)
    min_level = parse_level(min_level)
    if _module == nothing
        empty!(logger.module_limits)
        logger.default_min_level = min_level
    else
        # Per-module log limiting
        logger.module_limits[_module] = min_level
    end
    logger
end

function shouldlog(logger::InteractiveLogger, level, _module, group, id)
    if level < get(logger.module_limits, _module, logger.default_min_level)
        return false
    end
    if id in logger.blacklisted_ids
        return false
    end
    return true
end

function min_enabled_level(logger::InteractiveLogger)
    min_level = logger.default_min_level
    for (_,level) ∈ logger.module_limits
        if level < min_level
            min_level = level
        end
    end
    return min_level
end

catch_exceptions(logger::InteractiveLogger) = logger.catch_exceptions

formatmsg(msg) = string(msg)

# Formatting of exceptions
formatmsg(e::Tuple{Exception,Any}) = formatmsg(e[1], e[2])
function formatmsg(ex::Exception)
    bt = catch_backtrace()
    formatmsg(ex, isempty(bt) ? nothing : bt)
end
function formatmsg(ex_msg::Exception, bt)
    io = IOBuffer()
    showerror(io, ex_msg, bt; backtrace = bt!=nothing)
    String(take!(io))
end

function levelstyle(level::LogLevel)
    if     level < Info  return ((:blue,   true), "D-")
    elseif level < Warn  return ((:cyan,   true), "I-")
    elseif level < Error return ((:yellow, true), "W-")
    else                 return ((:red,    true), "E-")
    end
end

function handle_message(logger::InteractiveLogger, level, msg, _module, group,
                        id, file, line; kwargs...)
    # TODO: filter maxlog here...
    handle_message(logger, level, formatmsg(msg), _module, group, id,
                   file, line; kwargs...)
end

function handle_message(logger::InteractiveLogger, level, msg::AbstractString,
                        _module, group, id, filepath, line;
                        max_log=nothing, maxlog=nothing, kwargs...)
    if maxlog == nothing && max_log != nothing
        Base.depwarn("the max_log logging keyword has become maxlog", :handle_message)
        maxlog = max_log
    end
    if maxlog !== nothing
        count = get!(logger.message_counts, id, 0)
        count += 1
        logger.message_counts[id] = count
        if count > maxlog
            push!(logger.blacklisted_ids, id)
            return
        end
    end
    display_message(logger, level, msg, _module, group, id, filepath, line; kwargs...)
    nothing
end

# Print a string with prefixes on each line, and a suffix on the last line
function print_with_decorations(io, prefixes, decoration_color, str, suffix)
    @assert !isempty(prefixes)
    p = 1
    splitter = r"(\n|\e\[[0-9;]*m)"
    i = firstindex(str)
    n = lastindex(str)
    colorcode = ""
    stylecode = ""
    printstyled(io, prefixes[p], bold=decoration_color[2], color=decoration_color[1])
    linelen = length(prefixes[p])
    p == length(prefixes) || (p += 1)
    print(io, colorcode, stylecode)
    while i <= n
        r = findnext(splitter, str, i)
        @static if VERSION < v"0.7"
            if isempty(r)
                r = nothing
            end
        end
        j = (r === nothing) ? n : last(r)
        linelen += length(SubString(str, i, ((r === nothing) ? n : prevind(str,first(r)))))
        print(io, SubString(str, i, j))
        i = nextind(str,j)
        if r !== nothing
            if str[r[1]] == '\n'
                !Base.have_color || print(io, "\e[0m")
                printstyled(io, prefixes[p], bold=decoration_color[2], color=decoration_color[1])
                linelen = length(prefixes[p])
                p == length(prefixes) || (p += 1)
                !Base.have_color || print(io, "\e[0m")
                print(io, colorcode, stylecode)
            else
                code = str[r]
                if code == "\e[m"
                elseif code == "\e[22m"
                    stylecode = ""
                elseif code == "\e[39m"
                    colorcode = ""
                else
                    icode = parse(Int, code[3:end-1])
                    if 1 <= icode <= 29
                        stylecode = code
                    elseif 30 <= icode <= 38
                        colorcode = code
                    end
                end
            end
        end
    end
    if !isempty(suffix)
        width = displaysize(io)[2]
        if length(suffix) > width - linelen
            printstyled(io, "\n", prefixes[p], bold=decoration_color[2], color=decoration_color[1])
            linelen = length(prefixes[p])
        end
        print(io, " "^max(0, (width - linelen - length(suffix))))
        printstyled(io, suffix, bold=decoration_color[2], color=decoration_color[1])
    end
end

function display_message(logger::InteractiveLogger, level, msg::AbstractString,
                         _module, group, id, filepath, line;
                         progress=nothing, once=nothing, kwargs...)
    # TODO: progress throttling?
    # Log printing
    filename = basename(String(filepath))
    color, prefix = levelstyle(convert(LogLevel, level))
    # Attempt at avoiding the problem of distracting metadata in info log
    # messages - print metadata to the right hand side.
    metastr = " $level $filename:$line"
    msg = rstrip(msg, '\n')

    # Hack: render markdown to temporary buffer, and munge it to remove
    # trailing newlines, even with confounding color codes.
    buf = IOBuffer()
    dsize = displaysize(logger.stream)
    width = dsize[2]
    Markdown.term(IOContext(buf, :displaysize=>(dsize[1],width-2)), Markdown.parse(msg))
    msg = String(take!(buf))
    msg = replace(msg, r"\n(\e\[[0-9]+m)$"=>s"\1")
    msg = rstrip(msg)

    if progress === nothing
        if logger.prev_progress_key !== nothing
            print(logger.stream, "\n")
        end
        logger.prev_progress_key = nothing

        print_with_decorations(logger.stream, [prefix, "| "], color, msg, metastr)
        print(logger.stream, "\n")

        for (k,v) in kwargs
            valstr = formatmsg(v)
            kvstr = '\n' in valstr ? string(k, " =\n", formatmsg(v)) :
                                     string(k, " = ", formatmsg(v))
            print_with_decorations(logger.stream, ["|   ", "|    "], color, kvstr, "")
            print(logger.stream, "\n")
        end
    else
        progress_key = msg
        if logger.prev_progress_key !== nothing && logger.prev_progress_key != progress_key
            print(logger.stream, "\n")
        end
        nbar = max(1, width - (termlength(msg) + length(prefix) + length(metastr)) - 3)
        nfilledbar = round(Int, clamp(progress, 0, 1)*nbar)
        msgbar = string(msg, " [", "-"^nfilledbar, " "^(nbar - nfilledbar), "]")
        print(logger.stream, "\r")
        print_with_decorations(logger.stream, [prefix], color, msgbar, metastr)
        logger.prev_progress_key = progress_key
    end
end


text:
## Classification of MNIST dataset
## with the convolutional neural network know as LeNet5.
## This script also combines various
## packages from the Julia ecosystem  with Flux.
using Flux
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, WeightDecay
using Flux: onehotbatch, onecold, logitcrossentropy
using Statistics, Random
using Parameters: @with_kw
using Logging: with_logger, global_logger
using TensorBoardLogger: TBLogger, tb_overwrite, set_step!, set_step_increment!
import ProgressMeter
import MLDatasets
import DrWatson: savename, struct2dict
import BSON
using CUDAapi

# LeNet5 "constructor".
# The model can be adapted to any image size
# and number of output classes.
function LeNet5(; imgsize=(28,28,1), nclasses=10)
    out_conv_size = (imgsize[1]÷4 - 3, imgsize[2]÷4 - 3, 16)

    return Chain(
            x -> reshape(x, imgsize..., :),
            Conv((5, 5), imgsize[end]=>6, relu),
            MaxPool((2, 2)),
            Conv((5, 5), 6=>16, relu),
            MaxPool((2, 2)),
            x -> reshape(x, :, size(x, 4)),
            Dense(prod(out_conv_size), 120, relu),
            Dense(120, 84, relu),
            Dense(84, nclasses)
          )
end

function get_data(args)
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32, dir=args.datapath)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32, dir=args.datapath)

    # MLDatasets uses HWCN format, Flux works with WHCN
    xtrain = permutedims(reshape(xtrain, 28, 28, 1, :), (2, 1, 3, 4))
    xtest = permutedims(reshape(xtest, 28, 28, 1, :), (2, 1, 3, 4))

    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    train_loader = DataLoader(xtrain, ytrain, batchsize=args.batchsize, shuffle=true)
    test_loader = DataLoader(xtest, ytest,  batchsize=args.batchsize)

    return train_loader, test_loader
end

loss(ŷ, y) = logitcrossentropy(ŷ, y)

function eval_loss_accuracy(loader, model, device)
    l = 0f0
    acc = 0
    ntot = 0
    for (x, y) in loader
        x, y = x |> device, y |> device
        ŷ = model(x)
        l += loss(ŷ, y) * size(x)[end]
        acc += sum(onecold(ŷ |> cpu) .== onecold(y |> cpu))
        ntot += size(x)[end]
    end
    return (loss = l/ntot |> round4, acc = acc/ntot*100 |> round4)
end

## utility functions

num_params(model) = sum(length, Flux.params(model))

round4(x) = round(x, digits=4)

# arguments for the `train` function
@with_kw mutable struct Args
    η = 3e-4             # learning rate
    λ = 0                # L2 regularizer param, implemented as weight decay
    batchsize = 128      # batch size
    epochs = 20          # number of epochs
    seed = 0             # set seed > 0 for reproducibility
    cuda = true          # if true use cuda (if available)
    infotime = 1 	     # report every `infotime` epochs
    checktime = 5        # Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    tblogger = false      # log training with tensorboard
    savepath = nothing    # results path. If nothing, construct a default path from Args. If existing, may overwrite
    datapath = joinpath(homedir(), "Datasets", "MNIST") # data path: change to your data directory
end

function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.cuda && CUDAapi.has_cuda_gpu()
    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ## DATA
    train_loader, test_loader = get_data(args)
    @info "Dataset MNIST: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    model = LeNet5() |> device
    @info "LeNet5 model: $(num_params(model)) trainable params"

    ps = Flux.params(model)

    opt = ADAM(args.η)
    if args.λ > 0
        opt = Optimiser(opt, WeightDecay(args.λ))
    end

    ## LOGGING UTILITIES
    if args.savepath == nothing
        experiment_folder = savename("lenet", args, scientific=4,
                    accesses=[:batchsize, :η, :seed, :λ]) # construct path from these fields
        args.savepath = joinpath("runs", experiment_folder)
    end
    if args.tblogger
        tblogger = TBLogger(args.savepath, tb_overwrite)
        set_step_increment!(tblogger, 0) # 0 auto increment since we manually set_step!
        @info "TensorBoard logging at \"$(args.savepath)\""
    end

    function report(epoch)
        train = eval_loss_accuracy(train_loader, model, device)
        test = eval_loss_accuracy(test_loader, model, device)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
        if args.tblogger
            set_step!(tblogger, epoch)
            with_logger(tblogger) do
                @info "train" loss=train.loss  acc=train.acc
                @info "test"  loss=test.loss   acc=test.acc
            end
        end
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch in 1:args.epochs
        p = ProgressMeter.Progress(length(train_loader))

        for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                ŷ = model(x)
                loss(ŷ, y)
            end
            Flux.Optimise.update!(opt, ps, gs)
            ProgressMeter.next!(p)   # comment out for no progress bar
        end

        epoch % args.infotime == 0 && report(epoch)
        if args.checktime > 0 && epoch % args.checktime == 0
            !ispath(args.savepath) && mkpath(args.savepath)
            modelpath = joinpath(args.savepath, "model.bson")
            let model=cpu(model), args=struct2dict(args)
                BSON.@save modelpath model epoch args
            end
            @info "Model saved in \"$(modelpath)\""
        end
    end
end

## Execution as a script
if abspath(PROGRAM_FILE) == @__FILE__
    train()
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
using AnimatedPlots

log_graph = StaticGraph(log, color=SFML.magenta)
exponential_graph = StaticGraph(x -> 2^x, color=SFML.cyan)

plot(log_graph)
plot(exponential_graph)

waitfor(current_window())


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Polytope{N,Dim,T}

We say that a geometry is a N-polytope when it is a collection of "flat" sides
that constitue a `N`-dimensional subspace. They are called polygon and polyhedron
respectively for 2D (`N=2`) and 3D (`N=3`) subspaces, embedded in a `Dim`-dimensional
space. The parameter `N` is also known as the rank or parametric dimension of the
polytope: https://en.wikipedia.org/wiki/Abstract_polytope.

The term polytope expresses a particular combinatorial structure. A polyhedron,
for example, can be decomposed into faces. Each face can then be decomposed into
edges, and edges into vertices. Some conventions act as a mapping between vertices
and higher dimensional features (edges, faces, cells...), removing the need to
store all features. We follow the ordering conventions of the GMSH project:
https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering

Additionally, the following property must hold in order for a geometry to be considered
a polytope: the boundary of a (N+1)-polytope is a collection of N-polytopes, which may
have (N-1)-polytopes in common. See https://en.wikipedia.org/wiki/Polytope.
"""
abstract type Polytope{N,Dim,T} <: Geometry{Dim,T} end

(::Type{PL})(vertices::Vararg{P}) where {PL<:Polytope,P<:Point} = PL(SVector(vertices))
(::Type{PL})(vertices::AbstractVector{TP}) where {PL<:Polytope,TP<:Tuple} = PL(Point.(vertices))
(::Type{PL})(vertices::Vararg{TP}) where {PL<:Polytope,TP<:Tuple} = PL(collect(vertices))

# type aliases for convenience
const Polygon = Polytope{2}
const Polyhedron = Polytope{3}
area(p::Polygon) = measure(p)
volume(p::Polyhedron) = measure(p)

"""
    paramdim(polytope)

Return the parametric dimension or rank of the polytope.
"""
paramdim(::Type{<:Polytope{N}}) where {N} = N

"""
    vertices(polytope)

Return the vertices of a `polytope`.
"""
vertices(p::Polytope) = p.vertices

"""
    nvertices(polytope)

Return the number of vertices in the `polytope`.
"""
nvertices(p::Polytope) = length(p.vertices)

"""
    facets(polytope, convention=GMSH)

Return the facets of a `polytope` according to some
ordering `convention`. Default convention is [`GMSH`](@ref).
See https://en.wikipedia.org/wiki/Facet_(geometry)
"""
function facets(p::Polytope{N}, convention=GMSH) where {N}
  faces(p, N-1, convention)
end

"""
    faces(polytope, rank, convention=GMSH)

Return the faces of the `polytope` of given `rank` according
to some ordering `convention`. Default convention is [`GMSH`](@ref).
"""
function faces(p::Polytope, rank, convention=GMSH)
  (materialize(ord, p.vertices) for ord in connectivities(typeof(p), rank, convention))
end

"""
    p1 == p2

Tells whether or not polytopes `p1` and `p2` are equal.
"""
==(p1::Polytope, p2::Polytope) = p1.vertices == p2.vertices

"""
    centroid(polytope)

Return the centroid of the `polytope`.
"""
centroid(p::Polytope) = Point(sum(coordinates.(p.vertices)) / length(p.vertices))

function Base.show(io::IO, p::Polytope)
  kind = nameof(typeof(p))
  vert = join(p.vertices, ", ")
  print(io, "$kind($vert)")
end

function Base.show(io::IO, ::MIME"text/plain", p::Polytope{N,Dim,T}) where {N,Dim,T}
  kind = nameof(typeof(p))
  lines = ["  └─$v" for v in p.vertices]
  println(io, "$kind{$Dim,$T}")
  print(io, join(lines, "\n"))
end

include("polytopes/segment.jl")
include("polytopes/triangle.jl")
include("polytopes/quadrangle.jl")
include("polytopes/pyramid.jl")
include("polytopes/tetrahedron.jl")
include("polytopes/hexahedron.jl")
include("polytopes/chain.jl")
include("polytopes/polyarea.jl")