text:
using Soss: predictive

function posterior_predictive(model, posterior_draws; kwargs...)
    pred = predictive(model, keys(first(posterior_draws))...)(; kwargs...)
    [rand(pred(draw)) for draw in posterior_draws]
end

function posterior_predictive(model::DynamicPPL.Model, posterior_draws; variables=(:complaints, ), kwargs...)
    predictions = predict(model, posterior_draws)
    groups = (eachrow(Array(MCMCChains.group(predictions, k))) for k in variables)
    return map(zip(groups...)) do vals
        NamedTuple{variables}(vals)
    end
end


text:
export ft_clrrts

function ft_clrrts(ft_handle::Culong)
  ft_status = ccall((:FT_ClrRts, d2xx),Cuint,(Culong,),ft_handle)
  checkstatus(ft_status)
  return nothing
end

ext: jl


text:
################################################################################
#
#   CONSTRAINT RELATED
#
################################################################################

# function to compute the hard spin constraint
function constraintFunction(
    k :: Vector{Float64},
    eigenvectors :: Vector{Vector{Float64}},
    alpha :: Vector{Float64},
    sites :: Vector{S},
    d_spin :: Int64
    ) :: Float64 where {L<:Integer,D,S<:AbstractSite{L, D}}

    d_eig   = length(eigenvectors)
    N_sites = length(sites)

    sum = 0
    spin_vectors = Vector{Float64}[]

    # iterate over all sites in the lattice
    for r in 1 : N_sites

        site  = sites[r]
        index = label(site)
        pos   = point(site)

        spin_vector = zeros(Float64, d_spin)
        product     = dot(pos, k)

        # iterate over the (degenerate) subset of eigenvectors with minimal energy for vector k
        for s in 1 : d_eig
            spin_vector .+= (alpha[s] * cos(product) + alpha[s + d_eig] * sin(product)) * eigenvectors[s][d_spin * (index - 1) + 1 : d_spin * index]
        end

        # compute deviation from unit length
        sum += abs(norm(spin_vector) - 1.0)

    end

    # average over all sites
    return sum / N_sites

end

# function to compute the constraint for a given k on a finite lattice
function computeConstraint(
    k :: Vector{Float64},
    hamiltonian :: LTHamiltonian{L,U,HB},
    sites :: Vector{S},
    d_spin :: Int64
    ) :: Float64 where {L<:Integer,NS,U,HB<:AbstractBondHamiltonian{L,NS},D,S<:AbstractSite{L, D}}

    # compute the interaction matrix
    M_k = getMatrixAtK(hamiltonian, k)

    # diagonalize it
    eigenfactorization = eigen(M_k)
    eigenvalues        = real(eigenfactorization.values)
    eigenvectors       = real(eigenfactorization.vectors)

    # get (degenerate) subset of minimal energy eigenvectors
    E_min      = minimum(eigenvalues)
    deviations = [abs(E_min - E) for E in eigenvalues]
    degenerate = collect(1 : length(eigenvalues))[deviations .<= 1e-10]

    # minimize the constraint on the lattice
    alpha            = ones(2 * length(degenerate))
    constraint_value = Optim.minimum(Optim.optimize(x->constraintFunction(k, [eigenvectors[:, i] for i in degenerate], x, sites, d_spin), alpha, NelderMead(), Optim.Options(iterations = 10000)))

    return constraint_value

end

################################################################################
#
#   GENERAL FUNCTION
#
################################################################################

# function to compute the groundstate given a mesh of random k_vectors and connections between high symmetry points
function getLTGroundstateKSpace(
        N_random :: Int64,
        symmetric_k_vectors :: Vector{Vector{Float64}},
        hamiltonian :: LTHamiltonian{L,U,HB},
        sites :: Vector{S},
        ruc :: RU
        ;
        initial_tries :: Int64 = 10,
        d_spin :: Int64 = 3,
        epsilon :: Float64 = 1e-6,
        epsilon_k :: Float64 = 1e-10,
        groundstate_energy :: Float64 = Inf
    ) :: LTGroundstate where {L<:Integer,NS,U,HB<:AbstractBondHamiltonian{L,NS},D,S<:AbstractSite{L, D}, P,B,RU<:AbstractReciprocalUnitcell{P, B}}

    d_spatial    = length(point(sites[1]))
    bounds_lower = -2 * pi .* ones(Float64, d_spatial)
    bounds_upper = 2 * pi .* ones(Float64, d_spatial)

    # size of the symmetric mesh
    N_symmetric = length(symmetric_k_vectors)

    # check if the GS energy has to be calculated
    if groundstate_energy == Inf
        # search for 10 random points and compare the found minimum
        for tries in 1 : initial_tries
            # find a suitable starting point for the Newton algorithm
            k_start = zeros(Float64, d_spatial)
            for j in 1 : d_spatial
                comp = rand()
                k_start[j] = comp * bounds_lower[j] + (1 - comp) * bounds_upper[j]
            end
            # optimize the energy
            groundstate_energy = min(Optim.minimum(Optim.optimize(x->energy(hamiltonian, x), k_start)), groundstate_energy)
        end
        # print the groundstate_energy
        println("Groundstate energy is E_min = $(groundstate_energy)")
    end

    # initialize random points, minimize their energies and fold back to first BZ
    function dE(k :: Vector{Float64}) :: Float64
        M_k = getMatrixAtK(hamiltonian, k)
        eigenvalues = eigvals(M_k) .- groundstate_energy
        return minimum(eigenvalues .* eigenvalues)
    end

    random_k_vectors = zeros(Float64, N_random, d_spatial)
    index = 1
    while index <= N_random

        # initialize a random starting point
        k = zeros(Float64, d_spatial)
        for j in 1 : d_spatial
            comp = rand()
            k[j] = comp * bounds_lower[j] + (1 - comp) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = dE(k)
         # iterate i over 100 newton steps (maximum)
        for i in 1 : 100
            # check if the energy is already converged
            if e0 < epsilon
                # fold back to 1.BZ and save the k vector
                random_k_vectors[index, :] = shiftToFirstBZ(ruc, k)
                # increment the index
                index += 1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps = zeros(Float64, d_spatial)
            for j in 1 : d_spatial
                shift = zeros(Float64, d_spatial)
                shift[j] = epsilon_k
                H_eps[j] = dE(k .+ shift)
            end
            dH = (H_eps .- H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = dE(k)
        end
    end

    # check which points in the symmetric mesh also have minimal energy
    symmetric_k_vectors_min = Vector{Float64}[]
    for i in 1 : N_symmetric
        k = symmetric_k_vectors[i]
        if abs(energy(hamiltonian, k) - groundstate_energy) <= epsilon
            push!(symmetric_k_vectors_min, k)
        end
    end

    # iterate over all points and compute the contraint
    N_symmetric        = length(symmetric_k_vectors_min)
    N_total            = N_random + N_symmetric
    k_vectors          = zeros(Float64, N_total, d_spatial)
    constraint_values  = zeros(Float64, N_total)
    for i in 1 : N_total
        if i <= N_random
            k                    = random_k_vectors[i, :]
            constraint_values[i] = computeConstraint(k, hamiltonian, sites, d_spin)
            k_vectors[i, :]      = k
        else
            k                    = symmetric_k_vectors_min[i - N_random]
            constraint_values[i] = computeConstraint(k, hamiltonian, sites, d_spin)
            k_vectors[i, :]      = k
        end
    end

    # return LTGroundstate
    return LTGroundstate(groundstate_energy, k_vectors, constraint_values)

end


text:
# ------------------------------------------------------------------
# Code to simulate logit data and estimate nonparametrically
# Simplest reasonable example of nonparametric demand estimation
# Written by James Brand
# ------------------------------------------------------------------
using Statistics, NPDemand
using RCall, DataFrames
@rlibrary ggplot2

J = 2; # number of products
T =2000; # number of markets
beta = -0.4; # price coefficient
sdxi = 0.15; # standard deviation of xi

# Returns market shares, prices, instruments, and the market demand shock, respectively
s, p, z  = simulate_logit(J, T, beta, sdxi);
df = toDataFrame(s,p,z);

# Estimate demand nonparametrically
inv_sigma, designs = inverse_demand(df);

# Calculate price elasticities at realized prices and market shares
elast, jacobians = price_elasticity(inv_sigma, df, p);
true_elast = beta.*p.*(1 .- s[:,1]) # equation for own-price elasticities in logit model

# Plot kernel densities of estimated and true own-price elasticities
df2 = DataFrame(Estimate = elast, True = true_elast[:,1])
ggplot(df2, aes(x=:True))+
 geom_density(aes(x=:Estimate), color = "gray", linetype = "dashed") +
 geom_density(aes(x=:True), color = "black") + xlab("Elasticity") +
 ylab("Density")


text:
using CSV, DataFrames, Plots, StatsPlots, Missings
using MLDataUtils
using PyCall
using JSON, HTTP
using Dates

function getCoordinates(addressList)
    """
    getCoordinates(addressList::Array{String,1})

    returns a Vector of coordinates in (lat,lon) format as well as a data Array that contains all the data returned by nominatim
    """
    #nominatim requires as user agent and will block you otherwise
    HTTP.setuseragent!("Mozilla/5.0 (Windows NT 10.0; rv:68.0) Gecko/20100101 Firefox/68.0")
    coordinatesList = Array{Tuple{Float64,Float64},1}(undef, length(addressList)) # will contain (lat,lon) of addresses in addressList
    data = Array{Any,1}(undef, length(addressList)) # will contain all data nominatim gives for each request
    for (i, address) in enumerate(addressList)
        addressFormatted = join(split(address), "+")
        rawdata = HTTP.get(string(
            "https://nominatim.openstreetmap.org/search?q=",
            addressFormatted,
            "&format=json&limit=1",
        ))
        # exception if we aren't able to parse lat/long for any reason
        # we simply set the values to zero
        try
            data[i] = JSON.parse(String(rawdata.body))[1]
            coordinatesList[i] = parse(Float64, data[i]["lat"]), parse(Float64, data[i]["lon"])
            sleep(1) #nominatim will block you if you have more than 1 request per second
        catch e
            coordinatesList[i] = (0., 0.)
        end
    end
    return coordinatesList, data
end

P = download("https://docs.google.com/spreadsheets/d/1uv0IOtBCdpql3y926MtoT0wrHU3T7KGtLKIJ1FZMkjo/export?format=csv&gid=176397632", "responses.csv")

data = CSV.read("responses.csv");
# remove the last two columns where Marie wrote some stuff
data = data[:,1:7];

# rename the columns
rename!(data, [:Date, :Email, :Name, :TimeZone, :Org, :Position, :Hotdog]);

# remove entries with missing responses
dropmissing!(data);

# encode hotdog responses as bool for logistic regression later
hotdog_bool = Dict("Yes" => Float32(1.), "No" => Float32(0.));

# Create a new column
data.HotDogEncoding = map(val -> hotdog_bool[val], data.Hotdog);

# Generate encodings for each postiion type
data[!, :PosEncoding] = convertlabel(LabelEnc.Indices{Float32}, data.Position);

# Instead of running it for every row in the dataframe, do it on only
# the unique organizations to cut down time in half!
orgs = unique(data[:,:Org]);
coordinates = getCoordinates(orgs);
# Generate a mapping dictionary that will fill in blanks
coord_mapping = Dict(org => coord for (org, coord) in zip(orgs, coordinates[1]))
coord_mapping["Center for Astrophysics | Harvard & Smithsonian"] = (42.381405, -71.128036);
coord_mapping["Fritz-Haber Institute"] = (52.448343, 13.282765);
for key in ["U Paris Saclay/CNRS/ISMO", "CNRS/Université Paris Saclay", "University of Paris-Saclay"]
    coord_mapping[key] = (48.706617, 2.178913);
end
coord_mapping["Institut de Physique de Rennes"] = (48.116743, -1.640843);
coord_mapping["University of Reims Champagne-Ardenne"] = (49.255836, 4.040983);
coord_mapping["New England Biolab"] = (42.650274, -70.842718);
coord_mapping["Laboratoire de Physique des Lasers, CNRS-Université Sorbone Paris Nord"] = (48.956389, 2.340612);
coord_mapping["Leiden Observatory"] = (52.154681, 4.483809);
coord_mapping["Imperial College London (Earth Science & Engineering)"] = (51.498866, -0.175457);
coord_mapping["JILA, University of Colorado Boulder"] = coord_mapping["University of Colorado Boulder"];
coord_mapping["Laboratoire Interdiscipliniare de Physique (LIPhy)"] = coord_mapping["Laboratoire interdisciplinaire de Physique"];
coord_mapping["LMD  Ecole Molytechnique, France"] = (48.714226, 2.210970);
# Get lat, long for all the institutions
# coords = getCoordinates(data[:,:Org]);

data[!,:temp] = map(loc->coord_mapping[loc], data.Org);
lat_array = [];
long_array = [];
for coord in data.temp
    lat, long = coord
    push!(lat_array, lat)
    push!(long_array, long)
end
data[!,:Lat] = lat_array;
data[!,:Long] = long_array;

# Generate a relative time from first registration
reg_times = DateTime.(data.Date, dateformat"d/m/y H:M:S");
reg_times .-= reg_times[1]
relative_time = Dates.value.(reg_times) ./ (1000. * 60. * 60. * 24.);
data[!,:RelativeTime] = relative_time;

CSV.write("cleaned.csv", data);


text:
using ApproxFun, RatFun, Test

x = Fun()
f = cos(x) // exp(x)

@test f(0.1) ≈ cos(0.1)/exp(0.1)
@test 2f isa RationalFun
@test (2f)(0.1) ≈ 2f(0.1)
@test inv(f)(0.1) ≈ exp(0.1)/cos(0.1)
@test (f*f)(0.1) ≈ f(0.1)^2
@test (f/f)(0.1) ≈ 1


text:
using JLD2
using Base.Test

include("lookup3.jl")
include("internal.jl")
include("rw.jl")
include("recon.jl")
include("append.jl")
#include("customserialization.jl") currently broken due to #265