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

"""
    Filter(func)

Filters the table returning only the rows where the `func` returns true.

# Examples

```julia
Filter(row -> sum(row) > 10)
Filter(row -> row.a == true && row.b < 30)
```

## Notes

* The schema of the table is preserved by the transform.
"""
struct Filter{F} <: Stateless
  func::F
end

isrevertible(::Type{<:Filter}) = true

function apply(transform::Filter, table)
  rows = Tables.rowtable(table)

  # selected and rejected rows/inds
  sinds = findall(transform.func, rows)
  rinds = setdiff(1:length(rows), sinds)
  srows = rows[sinds]
  rrows = rows[rinds]

  newtable = srows |> Tables.materializer(table)
  return newtable, zip(rinds, rrows)
end

function revert(::Filter, newtable, cache)
  rows = Tables.rowtable(newtable)

  for (i, row) in cache
    insert!(rows, i, row)
  end

  rows |> Tables.materializer(newtable)
end

"""
    DropMissing()
    DropMissing(:)

Drop all rows with missing values in table.

    DropMissing(col₁, col₂, ..., colₙ)
    DropMissing([col₁, col₂, ..., colₙ])
    DropMissing((col₁, col₂, ..., colₙ))

Drop all rows with missing values in selected columns `col₁`, `col₂`, ..., `colₙ`.

    DropMissing(regex)

Drop all rows with missing values in columns that match with `regex`.

# Examples

```julia
DropMissing()
DropMissing("b", "c", "e")
DropMissing([2, 3, 5])
DropMissing((:b, :c, :e))
DropMissing(r"[bce]")
```

## Notes

* The transform can alter the element type of columns from `Union{Missing,T}` to `T`.
"""
struct DropMissing{S<:ColSpec} <: Stateless
  colspec::S
end

DropMissing(::Tuple{}) = throw(ArgumentError("Cannot create a DropMissing object with empty tuple."))

DropMissing() = DropMissing(:)

DropMissing(cols::T...) where {T<:ColSelector} =
  DropMissing(cols)

isrevertible(::Type{<:DropMissing}) = true

_ftrans(::DropMissing{Colon}, cols) =
  Filter(row -> all(!ismissing, row))

_ftrans(::DropMissing, cols) =
  Filter(row -> all(!ismissing, getindex.(Ref(row), cols)))

# nonmissing
_nonmissing(::Type{T}, x) where {T} = x
_nonmissing(::Type{Union{Missing,T}}, x) where {T} = collect(T, x)
_nonmissing(x) = _nonmissing(eltype(x), x)

function apply(transform::DropMissing, table)
  names = Tables.columnnames(table)
  types = Tables.schema(table).types
  snames = choose(transform.colspec, names)
  ftrans = _ftrans(transform, snames)
  newtable, fcache = apply(ftrans, table)

  # post-processing
  cols = Tables.columns(newtable)
  pcols = map(names) do n
    x = Tables.getcolumn(cols, n)
    n ∈ snames ? _nonmissing(x) : x
  end
  𝒯 = (; zip(names, pcols)...)
  ptable = 𝒯 |> Tables.materializer(newtable)

  ptable, (ftrans, fcache, types)
end

function revert(::DropMissing, newtable, cache)
  ftrans, fcache, types = cache

  # pre-processing
  cols = Tables.columns(newtable)
  names = Tables.columnnames(newtable)
  pcols = map(zip(types, names)) do (T, n)
    x = Tables.getcolumn(cols, n)
    collect(T, x)
  end
  𝒯 = (; zip(names, pcols)...)
  ptable = 𝒯 |> Tables.materializer(newtable)

  revert(ftrans, ptable, fcache)
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

struct UnitCell
    coordinateframe #this should be computed once when the UnitCell is created. It incorporates the parent coordinateframe transformation and the transformation from the parent to the unit cell.
    emitter
    # lens::LensAssembly
    detectors #have multiple eye positions that get computed all at once. Maybe. Lots of space so might be inefficient.
end

""" coordinateframe is a coordinate frame defined at a point in the repeating structure which all other positions will be measured with respect to. Used to define a globally consistent coordinate frame that pixels in each unit call map to."""
struct RepeatingStructure{N}
    coordinateframe #need to think about how to define this. Maybe using a Transform?
    cells::SVector{N,UnitCell} #may want to access cells not just by vector index, i.e., x,y, or something more complicated for hexagonal systems.
end

""" Maps a pixel position from an emitter in a unit cell into a globally consistent angular coordinate frame with coordinates (θ,ϕ). This allows the contributions from multiple cells to be summed correctly."""
function globalcoordinates(px,py,cell::UnitCell,repeat::RepeatingStructure)::NamedTuple{(:θ,:ϕ),Tuple{Real,Real}}
    #result will be like this (θ = val1, ϕ = val2)
end

function rectangulararray(nrows,ncols)
end

function hexagonalarray()
end

ext: jl


text:
# Equity Option example
# various pricing methods

using QuantLib

function main()
  cal = QuantLib.Time.TargetCalendar()
  todaysDate = Date(1998, 5, 15)
  settlementDate = Date(1998, 5, 17)

  set_eval_date!(settings, todaysDate)

  optionType = Put()
  underlying = 36.0
  strike = 40.0
  dividendYield = 0.00
  riskFreeRate = 0.06
  vol = 0.20
  mat = Date(1999, 5, 17)
  dc = QuantLib.Time.Actual365()

  exerciseDates = Date[settlementDate + Dates.Month(3 * i) for i = 1:4]

  europeanExercise = EuropeanExercise(mat)
  bermudanExercise = BermudanExercise(exerciseDates)
  americanExercise = AmericanExercise(settlementDate, mat)

  underlyingH = Quote(underlying)

  flatTermStructure = FlatForwardTermStructure(settlementDate, cal, Quote(riskFreeRate), dc)
  flatDividendTS = FlatForwardTermStructure(settlementDate, cal, Quote(dividendYield), dc)
  flatVolTS = BlackConstantVol(settlementDate, cal, vol, dc)

  payoff = PlainVanillaPayoff(optionType, strike)

  bsmProcess = BlackScholesMertonProcess(underlyingH, flatTermStructure, flatDividendTS, flatVolTS)

  # options
  bsPE = AnalyticEuropeanEngine(bsmProcess) # black scholes PE
  europeanOption = VanillaOption(payoff, europeanExercise, bsPE)

  # black scholes for European
  println("BlackScholes (Euro):                   ", npv(europeanOption))

  # heston process
  hestonProcess = HestonProcess(flatTermStructure, flatDividendTS, underlyingH, vol * vol, 1.0, vol * vol, 0.001, 0.0)
  hestonModel = HestonModel(hestonProcess)
  hestonPE = AnalyticHestonEngine(hestonModel)

  europeanOption = update_pricing_engine(europeanOption, hestonPE)
  println("Heston semi-analytic (Euro):           ", npv(europeanOption))

  # bates process
  batesProcess = BatesProcess(flatTermStructure, flatDividendTS, underlyingH, vol * vol, 1.0, vol * vol,
                              0.001, 0.0, 1e-14, 1e-14, 1e-14)
  batesModel = BatesModel(batesProcess)
  batesPE = BatesEngine(batesModel)

  europeanOption = update_pricing_engine(europeanOption, batesPE)
  println("Bates semi-analytic (Euro):            ", npv(europeanOption))

  # American Option
  baroneAdesiPE = BaroneAdesiWhaleyApproximationEngine(bsmProcess)
  americanOption = VanillaOption(payoff, americanExercise, baroneAdesiPE)
  println("Barone-Adesi/Whaley (Amer):            ", npv(americanOption))

  # Bjerksund Stensland engine
  bjerksundStenslandPE = BjerksundStenslandApproximationEngine(bsmProcess)
  americanOption = update_pricing_engine(americanOption, bjerksundStenslandPE)
  println("Bjerksund Stensland (Amer):            ", npv(americanOption))

  # Integral Engine
  integralPE = IntegralEngine(bsmProcess)
  europeanOption = update_pricing_engine(europeanOption, integralPE)
  println("Integral (Euro):                       ", npv(europeanOption))

  # FD Engines
  timeSteps = 801
  fdEuroPE = FDEuropeanEngine(bsmProcess, CrankNelson, timeSteps, timeSteps - 1)
  europeanOption = update_pricing_engine(europeanOption, fdEuroPE)
  println("FD (Euro):                             ", npv(europeanOption))

  # Bermudan Option
  fdBermudanPE = FDBermudanEngine(bsmProcess, CrankNelson, timeSteps, timeSteps - 1)
  bermudanOption = VanillaOption(payoff, bermudanExercise, fdBermudanPE)
  println("FD (Berm):                             ", npv(bermudanOption))

  fdAmericanPE = FDAmericanEngine(bsmProcess, CrankNelson, timeSteps, timeSteps - 1)
  americanOption = update_pricing_engine(americanOption, fdAmericanPE)
  println("FD (Amer):                             ", npv(americanOption))

  # Binomial Tree engines
  # Jarrow Rudd
  btJarrowRuddEuro = BinomialVanillaEngine(bsmProcess, timeSteps, JarrowRudd)
  europeanOption = update_pricing_engine(europeanOption, btJarrowRuddEuro)
  println("Binomial Jarrow-Rudd (Euro):           ", npv(europeanOption))

  btJarrowRuddBermudan = BinomialVanillaEngine(bsmProcess, timeSteps, JarrowRudd)
  bermudanOption = update_pricing_engine(bermudanOption, btJarrowRuddBermudan)
  println("Binomial Jarrow-Rudd (Berm):           ", npv(bermudanOption))

  btJarrowRuddAmerican = BinomialVanillaEngine(bsmProcess, timeSteps, JarrowRudd)
  americanOption = update_pricing_engine(americanOption, btJarrowRuddAmerican)
  println("Binomial Jarrow-Rudd (Amer):           ", npv(americanOption))

  # CoxRossRubinstein
  btCoxRubinsteinEuro = BinomialVanillaEngine(bsmProcess, timeSteps, CoxRossRubinstein)
  europeanOption = update_pricing_engine(europeanOption, btCoxRubinsteinEuro)
  println("Binomial Cox-Ross-Rubinstein (Euro):   ", npv(europeanOption))

  btCoxRubinsteinBermudan = BinomialVanillaEngine(bsmProcess, timeSteps, CoxRossRubinstein)
  bermudanOption = update_pricing_engine(bermudanOption, btCoxRubinsteinBermudan)
  println("Binomial Cox-Ross-Rubinstein (Berm):   ", npv(bermudanOption))

  btCoxRubinsteinAmerican = BinomialVanillaEngine(bsmProcess, timeSteps, CoxRossRubinstein)
  americanOption = update_pricing_engine(americanOption, btCoxRubinsteinAmerican)
  println("Binomial Cox-Ross-Rubinstein (Amer):   ", npv(americanOption))

  # AdditiveEQP
  btAdditiveEQPEuro = BinomialVanillaEngine(bsmProcess, timeSteps, AdditiveEQP)
  europeanOption = update_pricing_engine(europeanOption, btAdditiveEQPEuro)
  println("Additive Equiprobs (Euro):             ", npv(europeanOption))

  btAdditiveEQPBermudan = BinomialVanillaEngine(bsmProcess, timeSteps, AdditiveEQP)
  bermudanOption = update_pricing_engine(bermudanOption, btAdditiveEQPBermudan)
  println("Additive Equiprobs (Berm):             ", npv(bermudanOption))

  btAdditiveEQPAmerican = BinomialVanillaEngine(bsmProcess, timeSteps, AdditiveEQP)
  americanOption = update_pricing_engine(americanOption, btAdditiveEQPAmerican)
  println("Additive Equiprobs (Amer):             ", npv(americanOption))

  # Trigeorgis
  btTrigeorgisEuro = BinomialVanillaEngine(bsmProcess, timeSteps, Trigeorgis)
  europeanOption = update_pricing_engine(europeanOption, btTrigeorgisEuro)
  println("Binomial Trigeorgis (Euro):            ", npv(europeanOption))

  btTrigeorgisBermudan = BinomialVanillaEngine(bsmProcess, timeSteps, Trigeorgis)
  bermudanOption = update_pricing_engine(bermudanOption, btTrigeorgisBermudan)
  println("Binomial Trigeorgis (Berm):            ", npv(bermudanOption))

  btTrigeorgisAmerican = BinomialVanillaEngine(bsmProcess, timeSteps, Trigeorgis)
  americanOption = update_pricing_engine(americanOption, btTrigeorgisAmerican)
  println("Binomial Trigeorgis (Amer):            ", npv(americanOption))

  # Tian
  btTianEuro = BinomialVanillaEngine(bsmProcess, timeSteps, Tian)
  europeanOption = update_pricing_engine(europeanOption, btTianEuro)
  println("Binomial Tian (Euro):                  ", npv(europeanOption))

  btTianBermudan = BinomialVanillaEngine(bsmProcess, timeSteps, Tian)
  bermudanOption = update_pricing_engine(bermudanOption, btTianBermudan)
  println("Binomial Tian (Berm):                  ", npv(bermudanOption))

  btTianAmerican = BinomialVanillaEngine(bsmProcess, timeSteps, Tian)
  americanOption = update_pricing_engine(americanOption, btTianAmerican)
  println("Binomial Tian (Amer):                  ", npv(americanOption))

  # LeisenReimer
  btLeisenReimerEuro = BinomialVanillaEngine(bsmProcess, timeSteps, LeisenReimer)
  europeanOption = update_pricing_engine(europeanOption, btLeisenReimerEuro)
  println("Binomial Leisen-Reimer (Euro):         ", npv(europeanOption))

  btLeisenReimerBermudan = BinomialVanillaEngine(bsmProcess, timeSteps, LeisenReimer)
  bermudanOption = update_pricing_engine(bermudanOption, btLeisenReimerBermudan)
  println("Binomial Leisen-Reimer (Berm):         ", npv(bermudanOption))

  btLeisenReimerAmerican = BinomialVanillaEngine(bsmProcess, timeSteps, LeisenReimer)
  americanOption = update_pricing_engine(americanOption, btLeisenReimerAmerican)
  println("Binomial Leisen-Reimer (Amer):         ", npv(americanOption))

  # Joshi4
  btJoshi4Euro = BinomialVanillaEngine(bsmProcess, timeSteps, Joshi4)
  europeanOption = update_pricing_engine(europeanOption, btJoshi4Euro)
  println("Binomial Joshi (Euro):                 ", npv(europeanOption))

  btJoshi4Bermudan = BinomialVanillaEngine(bsmProcess, timeSteps, Joshi4)
  bermudanOption = update_pricing_engine(bermudanOption, btJoshi4Bermudan)
  println("Binomial Joshi (Berm):                 ", npv(bermudanOption))

  btJoshi4American = BinomialVanillaEngine(bsmProcess, timeSteps, Joshi4)
  americanOption = update_pricing_engine(americanOption, btJoshi4American)
  println("Binomial Joshi (Amer):                 ", npv(americanOption))

  ## MONTE CARLO ##
  timeSteps = 1
  mcSeed = 42
  mcengine1 = MCEuropeanEngine(bsmProcess; timeSteps = timeSteps, requiredTolerance = 0.02, seed = mcSeed)
  europeanOption = update_pricing_engine(europeanOption, mcengine1)
  println("MC-crude (Euro):                       ", npv(europeanOption))

  nSamples = 32768
  mcengine2 = MCEuropeanEngine(bsmProcess; timeSteps = timeSteps, requiredSamples = nSamples)
  europeanOption = update_pricing_engine(europeanOption, mcengine2)
  println("MC-Sobol (Euro):                       ", npv(europeanOption))

  mcengine3 = MCAmericanEngine(bsmProcess; timeSteps = 100, antitheticVariate=true, requiredTolerance=0.02, seed=mcSeed, nCalibrationSamples=4096)
  americanOption = update_pricing_engine(americanOption, mcengine3)
  println("MC-Longstaff Schwartz (Amer):          ", npv(americanOption))
  # americanOption
end


text:
module TestScitypes

using MLJ
import MLJBase
using Test

S = scitype(ConstantClassifier())
@test S().prediction_type == :probabilistic
U = scitype(FeatureSelector())
@test U().input_scitype == MLJ.Table(Scientific)

# XXX fix this in the next release cycle (Nov 5, 2019)
# M = scitype(rms)
# @test_broken M().prediction_type == :deterministic

for handle in localmodels()
    name = Symbol(handle.name)
    eval(quote
         scitype(($name)())
         end)
end

end
true


text:
#forward evaluation for a scalar function
function forward_pass_0ord{I,V}(ttstart::I, ttend::I, tt::Vector{I}, vvals::Vector{V}, pvals::Vector{V}, stk::Vector{Float64})
    idx = ttstart
    stklen = zero(I)

    @inbounds while(idx < ttend)
        ntype = tt[idx]
        # @show ntype
        # @show stklen, stk

        # if ntype == TYPE_P
        #     @assert false

        # else
        if ntype == TYPE_V
            stklen += 1
            @inbounds stk[stklen] = vvals[tt[idx+1]]
            idx += 3

        elseif ntype == TYPE_O
            @inbounds oc = tt[idx+2]
            @inbounds n = tt[idx+3]

            if n == 1 # 1-argument functions
                @inbounds stk[stklen] = eval_0ord(OP[oc],stk[stklen])
            elseif n == 2
                # @show OP[oc],stk
                @inbounds val = eval_0ord(OP[oc],stk[stklen-1],stk[stklen])
                stklen -= 1
                @inbounds stk[stklen] = val
            else
                @inbounds val = eval_0ord(OP[oc],stk,stklen-n+1,stklen)
                stklen -= n-1
                @inbounds stk[stklen] = val
            end
            idx += 5

        elseif ntype == TYPE_O1
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds stk[stklen] = eval_0ord(OP[oc],pval,stk[stklen])
            idx += 5

        elseif ntype == TYPE_O2
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds stk[stklen] = eval_0ord(OP[oc],stk[stklen],pval)
            idx += 5

        elseif ntype == TYPE_O3
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds vval = vvals[tt[idx+4]]
            stklen += 1
            @inbounds stk[stklen] = eval_0ord(OP[oc],pval,vval)
            idx += 6

        elseif ntype == TYPE_O4
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds vval = vvals[tt[idx+4]]
            stklen += 1
            @inbounds stk[stklen] = eval_0ord(OP[oc],vval,pval)
            idx += 6

        elseif ntype == TYPE_O5
            @inbounds oc = tt[idx+2]
            @inbounds vval = vvals[tt[idx+3]]
            @inbounds stk[stklen] = eval_0ord(OP[oc],vval,stk[stklen])
            idx += 5

        elseif ntype == TYPE_O6
            @inbounds oc = tt[idx+2]
            @inbounds vval = vvals[tt[idx+3]]
            @inbounds stk[stklen] = eval_0ord(OP[oc],stk[stklen],vval)
            idx += 5

        elseif ntype == TYPE_O7
            @inbounds oc = tt[idx+2]
            @inbounds vval = vvals[tt[idx+3]]
            stklen += 1
            @inbounds stk[stklen] = eval_0ord(OP[oc],vval)
            idx += 5

        # else
        #     @assert false
        end
    end
    @assert idx == ttend + 1
    return stk[1]
end

## Interface method
function feval{I,V}(ts::I, te::I, tt::Vector{I}, vvals::Vector{V}, pvals::Vector{V}, stk::Vector{V})
    return forward_pass_0ord(ts, te, tt, vvals,pvals, stk)
end


text:
module Runtests

using Modia
using Test

@time begin
    usePlotPackage("SilentNoPlot")
    include("$(Modia.ModiaLang.path)/test/runtests.jl")
    include("../examples/runexamples.jl")
    usePreviousPlotPackage()
end

end


text:
"""
Computes line profiles and associated quantities.
"""

const invSqrtPi = 1. / sqrt(π)

"""
Compute scaled complex complementary error function using
[Humlicek (1982, JQSRT 27, 437)](https://doi.org/10.1016/0022-4073(82)90078-4)
W4 rational approximations.
Here, z is defined as z = v + i * a, and returns w(z) = H(a,v) + i * L(a, v).
"""
function humlicek(z::Complex)
    s = abs(real(z)) + imag(z)
    if s > 15.0
        # region I
        w = im * invSqrtPi * z / (z * z - 0.5)
    elseif s > 5.5
        # region II
        zz = z * z
        w = im * (z * (zz * invSqrtPi - 1.4104739589)) / (0.75 + zz * (zz - 3.0))
    else
        x, y = real(z), imag(z)
        t = y - im * x
        if y >= 0.195 * abs(x) - 0.176
            # region III
            w = ((16.4955 + t * (20.20933 + t * (11.96482 + t * (3.778987 + 0.5642236 * t))))
               / (16.4955 + t * (38.82363 + t * (39.27121 + t * (21.69274 + t * (6.699398 + t))))))
        else
            # region IV
            u = t * t
            nom = t * (36183.31 - u * (3321.99 - u * (1540.787 -  u *
                   (219.031 - u * (35.7668 - u * (1.320522 - u * .56419))))))
            den = 32066.6 - u * (24322.8 - u * (9022.23 - u * (2186.18 -
                    u * (364.219 - u * (61.5704 - u * (1.84144 - u))))))
            w = exp(u) - nom / den
        end
    end
    return w
end

"""
    voigt_profile(a::AbstractFloat, v::AbstractFloat, ΔD::Unitful.Quantity)

Compute the normalised Voigt profile, given a damping constant `a`, dimensionless
wavelength or frequency `v`, and Doppler width `ΔD` (wavelength or frequency).
In the case of wavelength, v = (λ - λ0) / ΔD.
Uses Humlicek's W4 approximation. Returns in inverse units of ΔD.
"""
function voigt_profile(a::AbstractFloat, v::AbstractFloat, ΔD::Unitful.Quantity)
    z = v + a * im
    profile = real(humlicek(z))
    return profile * invSqrtPi / ΔD
end

"""
    dispersion_profile(a::AbstractFloat, v::AbstractFloat, ΔD::Unitful.Quantity)

Compute the normalised dispersion (or Faraday) profile, given a damping constant `a`,
dimensionless wavelength or frequency `v`, and Doppler width `ΔD` (wavelength or frequency).
In the case of wavelength, v = (λ - λ0) / ΔD.
Uses Humlicek's W4 approximation. Returns in inverse units of ΔD.
"""
function dispersion_profile(a::AbstractFloat, v::AbstractFloat, ΔλD::Unitful.Length)
    z = v + a * im
    profile = imag(humlicek(z))
    return profile * invSqrtPi / ΔλD
end

"""
    doppler_width(
        λ0::Unitful.Length,
        mass::Unitful.Mass,
        temperature::Unitful.Temperature
    )
    doppler_width(
        ν0::Unitful.Frequency,
        mass::Unitful.Mass,
        temperature::Unitful.Temperature
    )

Compute Doppler width in wavelength or frequency units, given a rest wavelength/frequency,
mass of atom/ion, and temperature. Returns in typical units for UV/optical/infrared lines.
"""
function doppler_width(
    λ0::Unitful.Length,
    mass::Unitful.Mass,
    temperature::Unitful.Temperature
)
    return (λ0 / c_0 * sqrt(2 * k_B * temperature / mass)) |> u"nm"
end

function doppler_width(
    ν0::Unitful.Frequency,
    mass::Unitful.Mass,
    temperature::Unitful.Temperature
)
    return (ν0 / c_0 * sqrt(2 * k_B * temperature / mass)) |> u"THz"
end


text:
"""
$(SIGNATURES)

Track the signal `signal` based on the current tracking `state`, the sample frequency
`sample_frequency` and PRN `prn`. Optional configurations are:
- Post correlation filter `post_corr_filter` defaults to `get_default_post_corr_filter(...)`
- Intermediate frequency `intermediate_frequency` defaults to 0Hz
- Maximal integration time `max_integration_time` defaults to 1ms. The actual integration
  time can be less in order to find the secondary code or the bit shift. Once this is found
  the integration time is identical to the maximal integration time.
- Minimal integration time `min_integration_time` defaults to 0.75ms. It's the minimal
  integration time, that leads to a valid correlation result. It is only used for the first
  integration period.
- Sample shift between early and late `early_late_sample_shift` defaults to
  `get_early_late_sample_shift(...)`
- Bandwidth of the carrier loop `carrier_loop_filter_bandwidth` defaults to 18Hz
- Bandwidth of the code loop `code_loop_filter_bandwidth` defaults to 1Hz
- Velocity aiding `velocity_aiding` defaults to 0Hz
"""
function track(
        signal,
        state::TrackingState{S, C, CALF, COLF, CN},
        prn::Integer,
        sample_frequency;
        post_corr_filter = get_default_post_corr_filter(get_correlator(state)),
        intermediate_frequency = 0.0Hz,
        max_integration_time::typeof(1ms) = 1ms,
        min_integration_time::typeof(1.0ms) = 0.75ms,
        early_late_sample_shift = get_early_late_sample_shift(S,
            get_correlator(state), sample_frequency, 0.5),
        carrier_loop_filter_bandwidth = 18Hz,
        code_loop_filter_bandwidth = 1Hz,
        velocity_aiding = 0Hz
) where {
    S <: AbstractGNSSSystem,
    C <: AbstractCorrelator,
    CALF <: AbstractLoopFilter,
    COLF <: AbstractLoopFilter,
    CN <: AbstractCN0Estimator
}
    if get_data_frequency(S) != 0Hz
        @assert rem(1 / get_data_frequency(S), max_integration_time) == 0ms
    end
    init_carrier_doppler = get_init_carrier_doppler(state)
    init_code_doppler = get_init_code_doppler(state)
    carrier_doppler = get_carrier_doppler(state)
    code_doppler = get_code_doppler(state)
    carrier_phase = get_carrier_phase(state)
    code_phase = get_code_phase(state)
    sc_bit_detector = get_sc_bit_detector(state)
    correlator = get_correlator(state)
    carrier_loop_filter = get_carrier_loop_filter(state)
    code_loop_filter = get_code_loop_filter(state)
    prompt_accumulator = get_prompt_accumulator(state)
    integrated_samples = get_integrated_samples(state)
    cn0_estimator = get_cn0_estimator(state)
    signal_start_sample = 1
    bit_buffer = BitBuffer()
    valid_correlator = zero(correlator)
    got_correlator = false
    while true
        num_samples_left_to_integrate = get_num_samples_left_to_integrate(
            S,
            max_integration_time,
            sample_frequency,
            code_doppler,
            code_phase,
            found(sc_bit_detector)
        )
        signal_samples_left = get_num_samples(signal) - signal_start_sample + 1
        num_samples_left = min(num_samples_left_to_integrate, signal_samples_left)
        carrier_frequency = get_current_carrier_frequency(
            intermediate_frequency,
            carrier_doppler
        )
        code_frequency = get_current_code_frequency(S, code_doppler)
        correlator = correlate(
            S,
            correlator,
            signal,
            prn,
            early_late_sample_shift,
            signal_start_sample,
            num_samples_left,
            carrier_frequency,
            code_frequency,
            sample_frequency,
            carrier_phase,
            code_phase
        )
        integrated_samples += num_samples_left
        carrier_phase = update_carrier_phase(
            num_samples_left,
            carrier_frequency,
            sample_frequency,
            carrier_phase
        )
        prev_code_phase = code_phase
        code_phase = update_code_phase(
            S,
            num_samples_left,
            code_frequency,
            sample_frequency,
            code_phase,
            found(sc_bit_detector)
        )
        integration_time = integrated_samples / sample_frequency
        if num_samples_left == num_samples_left_to_integrate &&
                integration_time >= min_integration_time
            got_correlator = true

            correlator = normalize(correlator, integrated_samples)
            valid_correlator = correlator
            filtered_correlator = filter(post_corr_filter, correlator)
            pll_discriminator = pll_disc(S, filtered_correlator)
            dll_discriminator = dll_disc(
                S,
                filtered_correlator,
                early_late_sample_shift,
                code_frequency / sample_frequency
            )
            carrier_freq_update, carrier_loop_filter = filter_loop(
                carrier_loop_filter,
                pll_discriminator,
                integration_time,
                carrier_loop_filter_bandwidth
            )
            code_freq_update, code_loop_filter = filter_loop(
                code_loop_filter,
                dll_discriminator,
                integration_time,
                code_loop_filter_bandwidth
            )
            carrier_doppler, code_doppler = aid_dopplers(
                S,
                init_carrier_doppler,
                init_code_doppler,
                carrier_freq_update,
                code_freq_update,
                velocity_aiding
            )
            cn0_estimator = update(cn0_estimator, get_prompt(filtered_correlator))
            bit_buffer, prompt_accumulator = buffer(
                S,
                bit_buffer,
                prompt_accumulator,
                found(sc_bit_detector),
                prev_code_phase,
                code_phase,
                max_integration_time,
                get_prompt(filtered_correlator)
            )
            sc_bit_detector = find(S, sc_bit_detector, get_prompt(filtered_correlator))
            correlator = zero(correlator)
            integrated_samples = 0
        end

        num_samples_left == signal_samples_left && break
        signal_start_sample += num_samples_left
    end
    next_state = TrackingState{S, C, CALF, COLF, CN}(
        init_carrier_doppler,
        init_code_doppler,
        carrier_doppler,
        code_doppler,
        carrier_phase,
        code_phase,
        correlator,
        carrier_loop_filter,
        code_loop_filter,
        sc_bit_detector,
        integrated_samples,
        prompt_accumulator,
        cn0_estimator
    )
    estimated_cn0 = estimate_cn0(cn0_estimator, max_integration_time)
    TrackingResults(next_state, valid_correlator, got_correlator, bit_buffer, estimated_cn0)
end

"""
$(SIGNATURES)

Default post corr filter
"""
function get_default_post_corr_filter(correlator::AbstractCorrelator{T}) where T <: SVector
    x -> x[1]
end

function get_default_post_corr_filter(correlator)
    x -> x
end

"""
$(SIGNATURES)

Returns the appropiate integration time. It will be the maximum integration time once the
secondary code or the bit shift has been found.
"""
function get_integration_time(
    ::Type{S},
    max_integration_time,
    secondary_code_or_bit_found::Bool
) where S <: AbstractGNSSSystem
    ifelse(
        secondary_code_or_bit_found,
        max_integration_time,
        min(
            convert(typeof(1ms), get_code_length(S) / get_code_frequency(S)),
            max_integration_time
        )
    )
end

"""
$(SIGNATURES)

Calculates the number of chips to integrate.
"""
function get_num_chips_to_integrate(
    ::Type{S},
    max_integration_time,
    current_code_phase,
    secondary_code_or_bit_found
) where S <: AbstractGNSSSystem
    max_phase = Int(upreferred(get_code_frequency(S) *
        get_integration_time(S, max_integration_time, secondary_code_or_bit_found)))
    current_phase_mod_max_phase = mod(current_code_phase, max_phase)
    max_phase - current_phase_mod_max_phase
end

"""
$(SIGNATURES)

Calculates the number of samples to integrate.
"""
function get_num_samples_left_to_integrate(
    ::Type{S},
    max_integration_time,
    sample_frequency,
    current_code_doppler,
    current_code_phase,
    secondary_code_or_bit_found
) where S <: AbstractGNSSSystem
    phase_to_integrate = get_num_chips_to_integrate(
        S,
        max_integration_time,
        current_code_phase,
        secondary_code_or_bit_found
    )
    code_frequency = get_code_frequency(S) + current_code_doppler
    ceil(Int, phase_to_integrate * sample_frequency / code_frequency)
end

"""
$(SIGNATURES)

Calculates the current carrier frequency.
"""
function get_current_carrier_frequency(intermediate_frequency, carrier_doppler)
    carrier_doppler + intermediate_frequency
end

"""
$(SIGNATURES)

Calculates the current code frequency.
"""
function get_current_code_frequency(::Type{S}, code_doppler) where S <: AbstractGNSSSystem
    code_doppler + get_code_frequency(S)
end

"""
$(SIGNATURES)

Updates the carrier phase.
"""
function update_carrier_phase(
    num_samples,
    carrier_frequency,
    sample_frequency,
    start_carrier_phase
)
    mod(
        carrier_frequency * num_samples / sample_frequency + start_carrier_phase + 0.5,
        1
    ) - 0.5
end

"""
$(SIGNATURES)

Updates the code phase.
"""
function update_code_phase(
    ::Type{S},
    num_samples,
    code_frequency,
    sample_frequency,
    start_code_phase,
    secondary_code_or_bit_found
) where S <: AbstractGNSSSystem
    secondary_code_or_bit_length =
        get_data_frequency(S) == 0Hz ?
        get_secondary_code_length(S) :
        Int(get_code_frequency(S) / (get_data_frequency(S) * get_code_length(S)))
    code_length = get_code_length(S) *
        (secondary_code_or_bit_found ? secondary_code_or_bit_length : 1)
    mod(code_frequency * num_samples / sample_frequency + start_code_phase, code_length)
end

"""
$(SIGNATURES)

Correlate the signal the the replicas.
"""
function correlate(
    ::Type{S},
    correlator,
    signal,
    prn::Integer,
    early_late_sample_shift::Integer,
    start_sample::Integer,
    num_samples::Integer,
    carrier_frequency,
    code_frequency,
    sample_frequency,
    start_carrier_phase,
    start_code_phase
) where S <: AbstractGNSSSystem
    total_code_length = get_code_length(S) * get_secondary_code_length(S)
    start_code_phase = mod(start_code_phase, total_code_length)
    carrier_phase_wrap = 0
    code_phase_wrap = 0
    code_register = 0
    max_early_late_sample_shift = maximum(early_late_sample_shift)
    for i = -max_early_late_sample_shift:max_early_late_sample_shift-1
        code = get_code(S, code_frequency * i / sample_frequency + start_code_phase, prn)
        code_register = code_register << 1 + (code + 1) >> 1
    end
    @inbounds @fastmath for i = 0:num_samples - 1
        carrier_phase = carrier_frequency * i / sample_frequency + start_carrier_phase
        carrier_phase_wrap += (carrier_phase - carrier_phase_wrap >= 0.5) * 1
        carrier = cis(-2π * (carrier_phase - carrier_phase_wrap))
        current_signal = get_signal(correlator, signal, start_sample + i)
        earliest_code_phase = code_frequency * (i + max_early_late_sample_shift) /
            sample_frequency + start_code_phase
        code_phase_wrap += (earliest_code_phase - code_phase_wrap >= total_code_length) *
            total_code_length
        earliest_code = get_code_unsafe(S, earliest_code_phase - code_phase_wrap, prn)
        code_register = code_register << 1 + (earliest_code + 1) >> 1
        correlator = correlate_iteration(
            S,
            correlator,
            current_signal,
            carrier,
            code_register,
            early_late_sample_shift,
            earliest_code
        )
    end
    correlator
end

"""
$(SIGNATURES)

Aid dopplers. That is velocity aiding for the carrier doppler and carrier aiding
for the code doppler.
"""
function aid_dopplers(
    ::Type{S},
    init_carrier_doppler,
    init_code_doppler,
    carrier_freq_update,
    code_freq_update,
    velocity_aiding
) where S <: AbstractGNSSSystem
    carrier_doppler = carrier_freq_update + velocity_aiding
    code_doppler = code_freq_update + carrier_doppler * get_code_center_frequency_ratio(S)
    init_carrier_doppler + carrier_doppler, init_code_doppler + code_doppler
end

"""
$(SIGNATURES)

Get the signal at sample `sample`.
"""
Base.@propagate_inbounds @inline function get_signal(
    correlator::AbstractCorrelator{T},
    signal,
    sample::Integer
) where T <: SVector
    T(view(signal, :, sample))
end

Base.@propagate_inbounds @inline function get_signal(correlator, signal, sample::Integer)
    signal[sample]
end

"""
$(SIGNATURES)

Get the number of samples in the signal.
"""
@inline function get_num_samples(signal)
    length(signal)
end

@inline function get_num_samples(signal::AbstractMatrix)
    size(signal, 2)
end


text:
module RulesDotProductOutTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_rules
import LinearAlgebra: dot

@testset "rules:typeof(dot):out" begin

     @testset "Belief Propagation: (m_in1::PointMass, m_in2::NormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0), meta = NoCorrection()), output = NormalMeanVariance(2.0, 2.0)),
            (input = (m_in1 = PointMass(-1.0), m_in2 = NormalMeanPrecision(3.0, 1.0), meta = NoCorrection()), output = NormalMeanVariance(-3.0, 1.0)),
            (input = (m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(2.0, 0.5), meta = NoCorrection()), output = NormalMeanVariance(8.0, 8.0)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0), meta = TinyCorrection()), output = NormalMeanVariance(2.0, 2.0)),
            (input = (m_in1 = PointMass(-1.0), m_in2 = NormalMeanPrecision(3.0, 1.0), meta = TinyCorrection()), output = NormalMeanVariance(-3.0, 1.0)),
            (input = (m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(2.0, 0.5), meta = TinyCorrection()), output = NormalMeanVariance(8.0, 8.0)),
        ]

    end

    @testset "Belief Propagation: (m_in1::NormalDistributionsFamily, m_in2::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = NormalMeanVariance(2.0, 2.0), m_in2 = PointMass(4.0), meta = NoCorrection()), output = NormalMeanVariance(8.0, 32.0)),
            (input = (m_in1 = NormalMeanPrecision(2.0, inv(3.0)), m_in2 = PointMass(2.0), meta = NoCorrection()), output = NormalMeanVariance(4.0, 12.0)),
            (input = (m_in1 = NormalWeightedMeanPrecision(2.0, 0.5), m_in2 = PointMass(1.0), meta = NoCorrection()), output = NormalMeanVariance(4.0, 2.0)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = NormalMeanVariance(2.0, 2.0), m_in2 = PointMass(4.0), meta = TinyCorrection()), output = NormalMeanVariance(8.0, 32.0)),
            (input = (m_in1 = NormalMeanPrecision(2.0, inv(3.0)), m_in2 = PointMass(2.0), meta = TinyCorrection()), output = NormalMeanVariance(4.0, 12.0)),
            (input = (m_in1 = NormalWeightedMeanPrecision(2.0, 0.5), m_in2 = PointMass(1.0), meta = TinyCorrection()), output = NormalMeanVariance(4.0, 2.0)),
        ]

    end

    @testset "Belief Propagation: (m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = MvNormalMeanCovariance([-1.0, 1.0], [2.0 -1.0; -1.0 4.0]), m_in2 = PointMass([4.0, 1.0]), meta = TinyCorrection()), output = NormalMeanVariance(-3, 28)),
            (input = (m_in1 = MvNormalMeanPrecision([2.0, 1.0], [2.0 -0.5; -0.5 5.0]), m_in2 = PointMass([2.0, 2.0]), meta = TinyCorrection()), output = NormalMeanVariance(6.0, 128/39)),
            (input = (m_in1 = MvNormalWeightedMeanPrecision([3.0, 2.0], [10.0 1.0; 1.0 20.0]), m_in2 = PointMass([-1.0, 3.0]), meta = TinyCorrection()), output = NormalMeanVariance(-7/199, 116/199)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = MvNormalMeanCovariance([-1.0, 1.0], [2.0 -1.0; -1.0 4.0]), m_in2 = PointMass([4.0, 1.0]), meta = NoCorrection()), output = NormalMeanVariance(-3, 28)),
            (input = (m_in1 = MvNormalMeanPrecision([2.0, 1.0], [2.0 -0.5; -0.5 5.0]), m_in2 = PointMass([2.0, 2.0]), meta = NoCorrection()), output = NormalMeanVariance(6.0, 128/39)),
            (input = (m_in1 = MvNormalWeightedMeanPrecision([3.0, 2.0], [10.0 1.0; 1.0 20.0]), m_in2 = PointMass([-1.0, 3.0]), meta = NoCorrection()), output = NormalMeanVariance(-7/199, 116/199)),
        ]

    end

    @testset "Belief Propagation: (m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = PointMass([4.0, 1.0]), m_in2 = MvNormalMeanCovariance([-1.0, 1.0], [2.0 -1.0; -1.0 4.0]), meta = TinyCorrection()), output = NormalMeanVariance(-3, 28)),
            (input = (m_in1 = PointMass([2.0, 2.0]), m_in2 = MvNormalMeanPrecision([2.0, 1.0], [2.0 -0.5; -0.5 5.0]), meta = TinyCorrection()), output = NormalMeanVariance(6.0, 128/39)),
            (input = (m_in1 = PointMass([-1.0, 3.0]), m_in2 = MvNormalWeightedMeanPrecision([3.0, 2.0], [10.0 1.0; 1.0 20.0]), meta = TinyCorrection()), output = NormalMeanVariance(-7/199, 116/199)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:out, Marginalisation) [
            (input = (m_in1 = PointMass([4.0, 1.0]), m_in2 = MvNormalMeanCovariance([-1.0, 1.0], [2.0 -1.0; -1.0 4.0]), meta = NoCorrection()), output = NormalMeanVariance(-3, 28)),
            (input = (m_in1 = PointMass([2.0, 2.0]), m_in2 = MvNormalMeanPrecision([2.0, 1.0], [2.0 -0.5; -0.5 5.0]), meta = NoCorrection()), output = NormalMeanVariance(6.0, 128/39)),
            (input = (m_in1 = PointMass([-1.0, 3.0]), m_in2 = MvNormalWeightedMeanPrecision([3.0, 2.0], [10.0 1.0; 1.0 20.0]), meta = NoCorrection()), output = NormalMeanVariance(-7/199, 116/199)),
        ]

    end

end

end


text:
using Stats, Distributions, Turing
using Gadfly

# include("ASCIIPlot.jl")
import Gadfly.ElementOrFunction

# First add a method to the basic Gadfly.plot function for QQPair types (generated by Distributions.qqbuild())
Gadfly.plot(qq::QQPair, elements::ElementOrFunction...) = Gadfly.plot(x=qq.qx, y=qq.qy, Geom.point, Theme(highlight_width=0px), elements...)

# Now some shorthand functions
qqplot(x, y, elements::ElementOrFunction...) = Gadfly.plot(qqbuild(x, y), elements...)
qqnorm(x, elements::ElementOrFunction...) = qqplot(Normal(), x, Guide.xlabel("Theoretical Normal quantiles"), Guide.ylabel("Observed quantiles"), elements...)

NSamples = 5000

@model gdemo_fw() = begin
  s ~ InverseGamma(2,3)
  # s = 1
  m ~ Normal(0,sqrt.(s))
  y ~ MvNormal([m; m; m], [sqrt.(s) 0 0; 0 sqrt.(s) 0; 0 0 sqrt.(s)])
end

@model gdemo_bk(x) = begin
  # Backward Step 1: theta ~ theta | x
  s ~ InverseGamma(2,3)
  # s = 1
  m ~ Normal(0,sqrt.(s))
  x ~ MvNormal([m; m; m], [sqrt.(s) 0 0; 0 sqrt.(s) 0; 0 0 sqrt.(s)])
  # Backward Step 2: x ~ x | theta
  y ~ MvNormal([m; m; m], [sqrt.(s) 0 0; 0 sqrt.(s) 0; 0 0 sqrt.(s)])
end

fw = IS(NSamples)
# bk = Gibbs(10, PG(10,10, :s, :y), HMC(1, 0.25, 5, :m));
bk_sample_n = 20
bk = NUTS(bk_sample_n, 100, 0.65);

s = sample(gdemo_fw(), fw);
# describe(s)

N = div(NSamples, bk_sample_n)

x = [s[:y][1]...]
s_bk = Array{Turing.Chain}(undef, N)

simple_logger = Base.CoreLogging.SimpleLogger(stderr, Base.CoreLogging.Debug)
with_logger(simple_logger) do
  i = 1
  while i <= N
    s_bk[i] = sample(gdemo_bk(x), bk);
    x = [s_bk[i][:y][end]...];
    i += 1
  end
end

s2 = vcat(s_bk...);
# describe(s2)

# qqplot(s[:m], s2[:m])
# qqplot(s[:s], s2[:s])

qqm = qqbuild(s[:m], s2[:m])

using UnicodePlots

qqm = qqbuild(s[:m], s2[:m])
show(scatterplot(qqm.qx, qqm.qy, title = "QQ plot for m", canvas = DotCanvas))
show(scatterplot(qqm.qx[51:end-50], qqm.qy[51:end-50], title = "QQ plot for m (removing first and last 50 quantiles):", canvas = DotCanvas))
show(scatterplot(qqm.qx, qqm.qy, title = "QQ plot for m"))
show(scatterplot(qqm.qx[51:end-50], qqm.qy[51:end-50], title = "QQ plot for m (removing first and last 50 quantiles):"))

qqs = qqbuild(s[:s], s2[:s])
show(scatterplot(qqs.qx, qqs.qy, title = "QQ plot for s"))
show(scatterplot(qqs.qx[51:end-50], qqs.qy[51:end-50], title = "QQ plot for s (removing first and last 50 quantiles):"))

X = qqm.qx
y = qqm.qy
slope = (1 / (transpose(X) * X)[1] * transpose(X) * y)[1]

print("  slopeₛ = $slope ≈ 1 (ϵ = 0.1)")
ans1 = abs(slope - 1.0) <= 0.1
if ans1
  printstyled(" ✓\n", color=:green)
else
  printstyled(" X\n", color=:red)
  printstyled("    slope = $slope, diff = $(slope - 1.0)\n", color=:red)
end

# qqs = qqbuild(s[:s], s2[:s])


text:
using SymbolicMath
using Test

@time @testset "symbol creation and equality" begin
    @test @symbol(x) == :x
    @test @symbols(y,z) == (:y, :z)
    @test x == :x
    @test :x == x
    @test x == x
end

@time @testset "binary operators" begin
    @symbols x y z
    for op in SymbolicMath.binary_opeators
        @test eval(op)(x, y) == :($op(x, y))
        @test eval(op)(x, :y) == :($op(x, y))
        @test eval(op)(:x, y) == :($op(x, y))
        @test eval(op)(x + y, z) == :($op(x + y, z))
        @test eval(op)(x, y + z) == :($op(x, y + z))
    end
end

@time @testset "unary operators" begin
    @symbols x y
    for op in SymbolicMath.unary_operators
        @test eval(op)(x) == :($op(x))
        @test eval(op)(x+y) == :($op(x+y))
    end
end

max_stars_repo_name: timkittel/SymbolicMath


text:
@testset "Schottky genus 4" begin
    τ = [0.95870734+0.73587725im 0.22092477+0.76863646im 0.53877459+0.87577267im 0.68177023+0.867436im; 0.22092477+0.76863646im 0.98812562+1.79674905im 0.54859032+1.10626215im 0.63310305+1.30158981im; 0.53877459+0.87577267im 0.54859032+1.10626215im 0.50173043+1.27729044im 0.49163557+1.33147334im; 0.68177023+0.867436im 0.63310305+1.30158981im 0.49163557+1.33147334im 0.35312207+1.60745975im];
    @test schottky_genus_4(τ) ≈ 2.16230-0.94677im atol=1e-3
end


text:
"""
`QueryFlags` correspond to the MongoDB wire protocol.
They may be bitwise or’d together.
They may modify how a query is performed in the MongoDB server.

These flags are passed as optional argument
for the aggregation function `Mongoc.aggregate`.

This data type mirrors C struct `mongoc_query_flags_t`.
See [libmongoc docs](http://mongoc.org/libmongoc/current/mongoc_query_flags_t.html)
for more information.

# Constants

`Mongoc.QUERY_FLAG_NONE`:
Specify no query flags.

`Mongoc.QUERY_FLAG_TAILABLE_CURSOR`:
Cursor will not be closed when the last data is retrieved. You can resume this cursor later.

`Mongoc.QUERY_FLAG_SLAVE_OK`:
Allow query of replica set secondaries.

`Mongoc.QUERY_FLAG_OPLOG_REPLAY`:
Used internally by MongoDB.

`Mongoc.QUERY_FLAG_NO_CURSOR_TIMEOUT`:
The server normally times out an idle cursor after an inactivity period (10 minutes).
This prevents that.

`Mongoc.QUERY_FLAG_AWAIT_DATA`:
Use with `Mongoc.MONGOC_QUERY_TAILABLE_CURSOR`. Block rather than returning no data.
After a period, time out.

`Mongoc.QUERY_FLAG_EXHAUST`:
Stream the data down full blast in multiple "reply" packets.
Faster when you are pulling down a lot of data and you know you want to retrieve it all.

`Mongoc.QUERY_FLAG_PARTIAL`:
Get partial results from mongos if some shards are down (instead of throwing an error).
"""
primitive type QueryFlags sizeof(Cint) * 8 end

"""
Adds one or more flags to the `FindAndModifyOptsBuilder`.
These flags can be *ORed* together,
as in `flags = Mongoc.FIND_AND_MODIFY_FLAG_UPSERT | Mongoc.FIND_AND_MODIFY_FLAG_RETURN_NEW`.

* `FIND_AND_MODIFY_FLAG_NONE`: Default. Doesn’t add anything to the builder.

* `FIND_AND_MODIFY_FLAG_REMOVE`: Will instruct find_and_modify to remove the matching document.

* `FIND_AND_MODIFY_FLAG_UPSERT`: Update the matching document or, if no document matches, insert the document.

* `FIND_AND_MODIFY_FLAG_RETURN_NEW`: Return the resulting document.
"""
primitive type FindAndModifyFlags sizeof(Cint) * 8 end

const QueryOrFindAndModifyFlags = Union{QueryFlags, FindAndModifyFlags}

Base.convert(::Type{T}, t::QueryOrFindAndModifyFlags) where {T<:Number} = reinterpret(Cint, t)
Base.convert(::Type{Q}, n::T) where {Q<:QueryOrFindAndModifyFlags, T<:Number} = reinterpret(Q, n)
Cint(flags::QueryOrFindAndModifyFlags) = convert(Cint, flags)
Base.:(|)(flag1::T, flag2::T) where {T<:QueryOrFindAndModifyFlags} = T( Cint(flag1) | Cint(flag2) )
Base.:(&)(flag1::T, flag2::T) where {T<:QueryOrFindAndModifyFlags} = T( Cint(flag1) & Cint(flag2) )

QueryFlags(u::Cint) = convert(QueryFlags, u)
QueryFlags(i::Number) = QueryFlags(Cint(i))
Base.show(io::IO, flags::QueryFlags) = print(io, "QueryFlags($(Cint(flags)))")

FindAndModifyFlags(u::Cint) = convert(FindAndModifyFlags, u)
FindAndModifyFlags(i::Number) = FindAndModifyFlags(Cint(i))
Base.show(io::IO, flags::FindAndModifyFlags) = print(io, "FindAndModifyFlags($(Cint(flags)))")

const QUERY_FLAG_NONE              = QueryFlags(0)
const QUERY_FLAG_TAILABLE_CURSOR   = QueryFlags(1 << 1)
const QUERY_FLAG_SLAVE_OK          = QueryFlags(1 << 2)
const QUERY_FLAG_OPLOG_REPLAY      = QueryFlags(1 << 3)
const QUERY_FLAG_NO_CURSOR_TIMEOUT = QueryFlags(1 << 4)
const QUERY_FLAG_AWAIT_DATA        = QueryFlags(1 << 5)
const QUERY_FLAG_EXHAUST           = QueryFlags(1 << 6)
const QUERY_FLAG_PARTIAL           = QueryFlags(1 << 7)

const FIND_AND_MODIFY_FLAG_NONE = FindAndModifyFlags(0)
const FIND_AND_MODIFY_FLAG_REMOVE = FindAndModifyFlags(1 << 0)
const FIND_AND_MODIFY_FLAG_UPSERT = FindAndModifyFlags(1 << 1)
const FIND_AND_MODIFY_FLAG_RETURN_NEW = FindAndModifyFlags(1 << 2)

# `URI` is a wrapper for C struct `mongoc_uri_t`."
mutable struct URI
    uri::String
    handle::Ptr{Cvoid}

    function URI(uri_string::String)
        err_ref = Ref{BSONError}()
        handle = mongoc_uri_new_with_error(uri_string, err_ref)
        if handle == C_NULL
            throw(err_ref[])
        end
        new_uri = new(uri_string, handle)
        finalizer(destroy!, new_uri)
        return new_uri
    end
end

mutable struct ClientPool
    uri::String
    handle::Ptr{Cvoid}

    function ClientPool(uri::URI; max_size::Union{Nothing, Integer}=nothing)
        client_pool_handle = mongoc_client_pool_new(uri.handle)
        @assert client_pool_handle != C_NULL "Failed to create client pool from URI $(uri.uri)."
        client_pool = new(uri.uri, client_pool_handle)
        finalizer(destroy!, client_pool)

        if max_size != nothing
            set_max_size(client_pool, max_size)
        end

        return client_pool
    end
end

# `Client` is a wrapper for C struct `mongoc_client_t`.
mutable struct Client{T<:Union{Nothing, ClientPool}}
    uri::String
    handle::Ptr{Cvoid}
    pool::T
end

function Client(uri::URI)
    client_handle = mongoc_client_new_from_uri(uri.handle)
    @assert client_handle != C_NULL "Failed to create client handle from URI $(uri.uri)."
    client = Client(uri.uri, client_handle, nothing)
    finalizer(destroy!, client)
    return client
end

function Client(pool::ClientPool; try_pop::Bool=false)

    local client_handle::Ptr{Cvoid}

    if try_pop
        client_handle = mongoc_client_pool_try_pop(pool.handle)
    else
        client_handle = mongoc_client_pool_pop(pool.handle)
    end

    @assert client_handle != C_NULL "Failed to create client handle from Pool."
    client = Client(pool.uri, client_handle, pool)
    finalizer(destroy!, client)
    return client
end

abstract type AbstractDatabase end
abstract type AbstractCollection end

# `Database` is a wrapper for C struct `mongoc_database_t`.
mutable struct Database <: AbstractDatabase
    client::Client
    name::String
    handle::Ptr{Cvoid}

    function Database(client::Client, db_name::String)
        db = new(client, db_name, mongoc_client_get_database(client.handle, db_name))
        finalizer(destroy!, db)
        return db
    end
end

mutable struct Bucket
    database::Database
    handle::Ptr{Cvoid}

    function Bucket(database::Database; options::Union{Nothing, BSON}=nothing)
        options_handle = options == nothing ? C_NULL : options.handle
        err_ref = Ref{BSONError}()

        gridfs_handle = mongoc_gridfs_bucket_new(database.handle,
                                                 options_handle,
                                                 C_NULL, err_ref)

        if gridfs_handle == C_NULL
            throw(err_ref[])
        end

        gridfs = new(database, gridfs_handle)
        finalizer(destroy!, gridfs)
        return gridfs
    end
end

# `Collection` is a wrapper for C struct `mongoc_collection_t`.
mutable struct Collection <: AbstractCollection
    database::Database
    name::String
    handle::Ptr{Cvoid}

    function Collection(database::Database, collection_name::String)
        collection_handle = mongoc_database_get_collection(database.handle, collection_name)
        @assert collection_handle != C_NULL "Failed to create a collection handle to $collection_name on db $(database.name)."
        collection = new(database, collection_name, collection_handle)
        finalizer(destroy!, collection)
        return collection
    end
end

const CursorSource = Union{Client, AbstractDatabase, AbstractCollection, Bucket}

# `Cursor` is a wrapper for C struct `mongoc_cursor_t`.
mutable struct Cursor{T<:CursorSource}
    source::T
    handle::Ptr{Cvoid}

    function Cursor(source::T, handle::Ptr{Cvoid}) where {T<:CursorSource}
        cursor = new{T}(source, handle)
        finalizer(destroy!, cursor)
        return cursor
    end
end

mutable struct BulkOperation
    collection::Collection
    handle::Ptr{Cvoid}
    executed::Bool

    function BulkOperation(collection::Collection; options::Union{Nothing, BSON}=nothing)
        options_handle = options == nothing ? C_NULL : options.handle
        handle = mongoc_collection_create_bulk_operation_with_opts(collection.handle, options_handle)
        @assert handle != C_NULL "Failed to create a bulk operation handle."
        bulk_operation = new(collection, handle, false)
        finalizer(destroy!, bulk_operation)
        return bulk_operation
    end
end

struct InsertOneResult{T}
    reply::BSON
    inserted_oid::T
end

struct BulkOperationResult{T}
    reply::BSON
    server_id::UInt32
    inserted_oids::Vector{T}
end

#
# Basic functions for types
#

function destroy!(uri::URI)
    if uri.handle != C_NULL
        mongoc_uri_destroy(uri.handle)
        uri.handle = C_NULL
    end
    nothing
end

function destroy!(client::Client{Nothing})
    if client.handle != C_NULL
        mongoc_client_destroy(client.handle)
        client.handle = C_NULL
    end
    nothing
end

function destroy!(client::Client{ClientPool})
    if client.handle != C_NULL
        mongoc_client_pool_push(client.pool.handle, client.handle)
        client.handle = C_NULL
    end
    nothing
end

function destroy!(client_pool::ClientPool)
    if client_pool.handle != C_NULL
        mongoc_client_pool_destroy(client_pool.handle)
        client_pool.handle = C_NULL
    end
    nothing
end

function destroy!(database::Database)
    if database.handle != C_NULL
        mongoc_database_destroy(database.handle)
        database.handle = C_NULL
    end
    nothing
end

function destroy!(gridfs::Bucket)
    if gridfs.handle != C_NULL
        mongoc_gridfs_bucket_destroy(gridfs.handle)
        gridfs.handle = C_NULL
    end
end

function destroy!(collection::Collection)
    if collection.handle != C_NULL
        mongoc_collection_destroy(collection.handle)
        collection.handle = C_NULL
    end
    nothing
end

function destroy!(cursor::Cursor)
    if cursor.handle != C_NULL
        mongoc_cursor_destroy(cursor.handle)
        cursor.handle = C_NULL
    end
    nothing
end

function destroy!(bulk_operation::BulkOperation)
    if bulk_operation.handle != C_NULL
        mongoc_bulk_operation_destroy(bulk_operation.handle)
        bulk_operation.handle = C_NULL
        bulk_operation.executed = true
    end
    nothing
end

#
# Session Types
#

mutable struct SessionOptions
    handle::Ptr{Cvoid}

    function SessionOptions(; casual_consistency::Bool=true)
        session_options_handle = mongoc_session_opts_new()
        @assert session_options_handle != C_NULL "Failed to create session options handle."
        session_options = new(session_options_handle)
        finalizer(destroy!, session_options)
        set_casual_consistency!(session_options, casual_consistency)
        return session_options
    end
end

mutable struct Session
    client::Client
    options::SessionOptions
    handle::Ptr{Cvoid}

    function Session(client::Client; options::SessionOptions=SessionOptions())
        err_ref = Ref{BSONError}()
        session_handle = mongoc_client_start_session(client.handle, options.handle, err_ref)
        if session_handle == C_NULL
            throw(err_ref[])
        end
        session = new(client, options, session_handle)
        finalizer(destroy!, session)
        return session
    end
end

struct DatabaseSession <: AbstractDatabase
    database::Database
    session::Session
end

struct CollectionSession <: AbstractCollection
    database_session::DatabaseSession
    collection::Collection
end

#=
struct BulkOperationSession
    collection_session::CollectionSession
    bulk_operation::BulkOperation

    function BulkOperationSession(collection_session::CollectionSession, options::Union{Nothing, BSON})
        options_with_session = _join(options, get_session(collection_session))
        bulk_operation = BulkOperation(collection_session.collection, options=options_with_session)
        return BulkOperationSession(collection_session, bulk_operation)
    end
end
=#

function destroy!(session_options::SessionOptions)
    if session_options.handle != C_NULL
        mongoc_session_opts_destroy(session_options.handle)
        session_options.handle = C_NULL
    end
    nothing
end

function destroy!(session::Session)
    if session.handle != C_NULL
        mongoc_client_session_destroy(session.handle)
        session.handle = C_NULL
    end
    nothing
end

mutable struct FindAndModifyOptsBuilder
    handle::Ptr{Cvoid}

    function FindAndModifyOptsBuilder(;
        update::Union{Nothing, BSON}=nothing,
        sort::Union{Nothing, BSON}=nothing,
        fields::Union{Nothing, BSON}=nothing,
        flags::Union{Nothing, FindAndModifyFlags}=nothing,
        bypass_document_validation::Bool=false,
    )

        opts = new(mongoc_find_and_modify_opts_new())
        finalizer(destroy!, opts)

        if update != nothing
            opts.update = update
        end

        if sort != nothing
            opts.sort = sort
        end

        if fields != nothing
            opts.fields = fields
        end

        if flags != nothing
            opts.flags = flags
        end

        opts.bypass_document_validation = bypass_document_validation

        return opts
    end
end

function destroy!(opts::FindAndModifyOptsBuilder)
    if opts.handle != C_NULL
        mongoc_find_and_modify_opts_destroy(opts.handle)
        opts.handle = C_NULL
    end
    nothing
end

#
# Streams
#

"""
# Interface

## Mutable fields

* `handle::Ptr{Cvoid}`

* `isopen::Bool`

* `timeout_msec::Int`

* `chunk_size::UInt`
"""
abstract type AbstractMongoStream <: IO end

function destroy!(stream::AbstractMongoStream)
    if stream.handle != C_NULL
        mongoc_stream_destroy(stream.handle)
        stream.handle = C_NULL
    end
end

# 30sec
const DEFAULT_TIMEOUT_MSEC = 30000

# 56kb of chunk size when reading from streams
const DEFAULT_CHUNK_SIZE = 56 * 1024

mutable struct MongoStreamFile{T} <: AbstractMongoStream
    owner::T
    handle::Ptr{Cvoid}
    isopen::Bool
    timeout_msec::Int
    chunk_size::UInt

    function MongoStreamFile(owner::T, handle::Ptr{Cvoid};
            timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
            chunk_size::Integer=DEFAULT_CHUNK_SIZE
        ) where {T}

        @assert handle != C_NULL
        stream = new{T}(owner, handle, true, Int(timeout_msec), UInt(chunk_size))
        finalizer(destroy!, stream)
        return stream
    end
end

mutable struct MongoIOStream{T} <: AbstractMongoStream
    owner::T
    handle::Ptr{Cvoid}
    isopen::Bool
    timeout_msec::Int
    chunk_size::UInt

    function MongoIOStream(owner::T, handle::Ptr{Cvoid};
            timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
            chunk_size::Integer=DEFAULT_CHUNK_SIZE
        ) where {T}

        @assert handle != C_NULL
        stream = new{T}(owner, handle, true, Int(timeout_msec), UInt(chunk_size))
        finalizer(destroy!, stream)
        return stream
    end
end


text:
# This code generates the 2l+l dimensional real irreducible representations of O(3).
#
# For a function defined on the sphere (f: S² → R) and represented in the basis of the
# real spherical harmonics, these matrix representations rotate/reflect the function
# about an axis in 3D space.
#
# The routines in this file are based upon:
#
#   Rotation Matrices for Real Spherical Harmonics, Direct Determination by Recursion
#   Joseph Ivanic and Klaus Ruedenberg
#   J. Phys. Chem., Vol. 100, No. 15, 1996
#
# The article contains numerous errors, and the "Additions and Corrections"
# — J. Phys. Chem. A, Vol. 102, No. 45, 1998 — also contains at least one error,
# for the matrices generated using Table 2 do not constitute an O(3) group
# homomorphism. A few empirical tests reveal that the term "sqrt(1 - δ(m,-1))" should
# probably be changed to "sqrt(1 + δ(m,-1))" in Table 2 for "V_mm'" and "m < 0".
#
# This file contains an optimized implementation of the (corrected) methods in
# the article. A reference implementation is provided in reference_irreps_O3.jl
#
# These methods are dependently typed on (and compile for) distinct values of l.
# Because of this, the first function call for a new value of l may be slow, but
# subsequent calls should be very fast.
#
# Compilation becomes rather slow beyond about l = 12, but if the goal is to rotate
# high bandwidth-limited functions on the sphere, then there exist better methods that
# do not involve explicitly constructing the entire rotation matrix.

δ(i::Int, j::Int) = i == j ? 1 : 0

function P(l::Int, i::Int, m::Int, n::Int)
    r = 2
    q = l

    if abs(n) < l
        :(R[$(r + i), $(r + 0)] * M[$(q + m), $(q + n)])
    elseif n == l
        :(R[$(r + i), $(r + 1)] * M[$(q + m), $(q + l - 1)] - R[$(r + i), $(r - 1)] * M[$(q + m), $(q - l + 1)])
    elseif n == -l
        :(R[$(r + i), $(r + 1)] * M[$(q + m), $(q - l + 1)] + R[$(r + i), $(r - 1)] * M[$(q + m), $(q + l - 1)])
    end
end

function irrep_O3_coeff(l::Int, m::Int, n::Int)
    nu = (l + m) * (l - m)
    nv = (1 + δ(m, 0)) * (l + abs(m) - 1) * (l + abs(m))
    nw = (l - abs(m) - 1) * (l - abs(m))
    d = abs(n) < l ? (l + n) * (l - n) : (2l) * (2l - 1)
    u = sqrt(nu / d)
    v = (1/2) * sqrt(nv / d) * (1 - 2δ(m, 0))
    w = (-1/2) * sqrt(nw / d) * (1 - δ(m, 0))

    U = P(l, 0, m, n)

    if m > 0
        V = :($(sqrt(1 + δ(m, 1))) * $(P(l, 1, m - 1, n)) - $(1 - δ(m, 1)) * $(P(l, -1, -m + 1, n)))
        W = :($(P(l, 1, m + 1, n)) + $(P(l, -1, -m - 1, n)))
    elseif m < 0
        V = :($(1 - δ(m, -1)) * $(P(l, 1, m + 1, n)) + $(sqrt(1 + δ(m, -1))) * $(P(l, -1, -m - 1, n)))
        W = :($(P(l, 1, m - 1, n)) - $(P(l, -1, -m + 1, n)))
    else
        V = :($(P(l, 1, 1, n)) + $(P(l, -1, -1, n)))
        W = :(0)
    end

    ex = :($(u == 0 ? 0 : :($u * ($U))) + $(v == 0 ? 0 : :($v * ($V))) + $(w == 0 ? 0 : :($w * ($W))))
    :(N[$(m + l + 1), $(n + l + 1)] = $(ex))
end

#@generated irrep_O3!{L}(N::AbstractMatrix, M::AbstractMatrix, R::AbstractMatrix, ::Type{Val{L}}) = Expr(:block, (irrep_O3_coeff(L, i, j) for i in -L:L, j in -L:L)..., :N)
@generated irrep_O3!{L}(N::AbstractMatrix, M::AbstractMatrix, R::AbstractMatrix, ::Type{Val{L}}) = Expr(:block, [irrep_O3_coeff(L, i, j) for i in -L:L, j in -L:L]..., :N)

"""
    irrep_O3!(N, M, R)

Generate an irrep of O(3) of dimension 2l+1 from irrep M of dimension 2l-1 and irrep R of dimension 3, and store it in matrix N.
"""
irrep_O3!{T}(N::AbstractMatrix{T}, M::AbstractMatrix{T}, R::AbstractMatrix{T}) = irrep_O3!(N, M, R, Val{fld(size(N, 1), 2)})

"""
    irrep_O3(M, R)

Generate an irrep of O(3) of dimension 2l+1 from irrep M of dimension 2l-1 and irrep R of dimension 3.
"""
irrep_O3{T}(M::AbstractMatrix{T}, R::AbstractMatrix{T}) = irrep_O3!(Matrix{T}(size(M, 1) + 2, size(M, 2) + 2), M, R)

"""
    irrep_O3(R, l)

Generate an irrep of O(3) of dimension 2l+1 from irrep R of dimension 3.
"""
irrep_O3{T}(R::AbstractMatrix{T}, l::Int) = l == 1 ? R : irrep_O3(irrep_O3(R, l - 1), R)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# benchmark/benchmarks.jl
# Runs benchmarks on all supported McCormick operators.
# Allocations are expected for the following operations:
# - MC(Float64, Interval, 1): Nonstatic definition of static type.
# - exp2, exp10, tanh, atanh, acosh: Underlying interval operation may allocate.
#############################################################################

using BenchmarkTools, IntervalArithmetic

# don't really want to measure interval performance, so we're using the
# fastest interval rounding mode
setrounding(Interval, :none)
using McCormick, ReverseMcCormick

const SUITE = BenchmarkGroup()

for T in (NS, Diff, MV)
    begin
    end
end


text:
const Ensemble = Vector{Path}

function Base.show(io::IO, ens::Ensemble)
    printlog(length(ens), " Trajectories")
end

function addEnsemble(sys::System; nSamples::Int = sys.routine.nSamples, method::Symbol = :Sobol, gn = 0)
    ens = Ensemble(undef, nSamples*sys.routine.nRepeat)
    if method == :Sobol
        if sys.routine.sampling == :Normal
            if sys.state.i == 1
                sob = SobolSeq(length(sys.model.bounds[1]), sys.model.bounds[1],  sys.model.bounds[2])
            else
                low = sys.state.θ .- 2 .* sqrt.(diag(sys.state.σ2))
                up = sys.state.θ .+ 2 .* sqrt.(diag(sys.state.σ2))
                sob = SobolSeq(length(sys.model.bounds[1]), low, up)
            end
            setEnsemble!(sys, ens, sob, gn)
        elseif sys.routine.sampling == :log
            if sys.state.i == 1
                sob = SobolSeq(length(sys.model.bounds[1]), log.(sys.model.bounds[1]),  log.(sys.model.bounds[2]))
            else
                low = sys.state.θ .- 2 .* sqrt.(diag(sys.state.σ2))
                up = sys.state.θ .+ 2 .* sqrt.(diag(sys.state.σ2))
                sob = SobolSeq(length(sys.model.bounds[1]), low, up)
            end
            setEnsembleLog!(sys, ens, sob, gn)
        end
    elseif method == :Dist
        if sys.routine.sampling == :Normal
            dist = MvNormal(sys.state.θ, sys.state.σ2)
            setEnsemble!(sys, ens, dist, gn)
        elseif sys.routine.sampling == :log
            dist = MvLogNormal(sys.state.θ, sys.state.σ2)
            setEnsemble!(sys, ens, dist, gn)
        end
    end
    return ens
end

function setEnsemble!(sys::System, ens::Ensemble, sob::Sobol.ScaledSobolSeq, gn)
    l = round(Int, length(ens)/sys.routine.nRepeat,RoundUp)
    j = 0
    for i in gn+1:gn+l
        ps = next!(sob)
        for r in 1:sys.routine.nRepeat
            j += 1
            if j > length(ens)
                return
            end
            ens[j] = Path(i, ps, sys.model.x0, sys.times[1], zeros(Int, sys.model.nr), length(sys.data))
        end
    end
end

function setEnsemble!(sys::System, ens::Ensemble, dist, gn)
    l = round(Int, length(ens)/sys.routine.nRepeat,RoundUp)
    j = 0
    for i in gn+1:gn+l
        ps = rand(dist)
        for r in 1:sys.routine.nRepeat
            j += 1
            if j > length(ens)
                return
            end
            ens[j] = Path(i, ps, sys.model.x0, sys.times[1], zeros(Int, sys.model.nr), length(sys.data))
        end
    end
end

function setEnsembleLog!(sys::System, ens::Ensemble, sob::Sobol.ScaledSobolSeq, gn)
    l = round(Int, length(ens)/sys.routine.nRepeat,RoundUp)
    j = 0
    for i in gn+1:gn+l
        ps = exp.(next!(sob))
        for r in 1:sys.routine.nRepeat
            j += 1
            if j > length(ens)
                return
            end
            ens[j] = Path(i, ps, sys.model.x0, sys.times[1], zeros(Int, sys.model.nr), length(sys.data))
        end
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
function ksdist{T<:Real, S<:Real}(x::AbstractVector{T}, y::AbstractVector{S})

  #adapted from HypothesisTest.jl
  n_x, n_y = length(x), length(y)
  sort_idx = sortperm([x; y])
  pdf_diffs = [ones(n_x)/n_x; -ones(n_y)/n_y][sort_idx]
  cdf_diffs = cumsum(pdf_diffs)
  δp = maximum(cdf_diffs)
  δn = -minimum(cdf_diffs)
  δ = max(δp, δn)

  return δ
end

function setupSMCparticles(ABCrejresults::ABCrejectionresults, ABCsetup)

  weights = ones(ABCsetup.nparticles)./ABCsetup.nparticles
  scales = (maximum(ABCrejresults.parameters, 1) -
                  minimum(ABCrejresults.parameters, 1) ./2)[:]

  particles = Array{ParticleSMC}(ABCsetup.nparticles)

  for i in 1:length(particles)
    particles[i] = ParticleSMC(ABCrejresults.particles[i].params, weights[1],
    scales, ABCrejresults.particles[i].distance,
    ABCrejresults.particles[i].other)
  end

  return particles, weights
end

function setupSMCparticles(ABCrejresults::ABCrejectionmodelresults, ABCsetup)

  weights = ones(ABCsetup.Models[1].nparticles)./ABCsetup.Models[1].nparticles
  scales = map(x -> collect((maximum(x, 1) -
                  minimum(x, 1) ./2)[:]), ABCrejresults.parameters)

  particles = Array{ParticleSMCModel}(ABCsetup.Models[1].nparticles)

  for i in 1:length(particles)
    particles[i] = ParticleSMCModel(ABCrejresults.particles[i].params,
    weights[1], scales[ABCrejresults.particles[i].model],
    ABCrejresults.particles[i].model, ABCrejresults.particles[i].distance,
    ABCrejresults.particles[i].other)
  end

  return particles, weights
end

function getscales(particles, ABCsetup::ABCSMC)

  parameters = hcat(map(x -> x.params, particles)...)'
  scales = ((maximum(parameters, 1) -
                  minimum(parameters, 1)) ./ABCsetup.scalefactor)[:]

  for i in 1:length(particles)
    particles[i].scales = scales
  end

  return particles
end

function getscales(particles, ABCsetup::ABCSMCModel)

  modelindex = trues(ABCsetup.nparticles, ABCsetup.nmodels)
  for i in 1:ABCsetup.nmodels
      modelindex[:, i] = map(x -> x.model, particles) .== i
  end

  modelfreq = sum(modelindex, 1)
  scales =  Array{Float64,1}[]

  for i in 1:ABCsetup.nmodels
    if modelfreq[i] == 0
      push!(scales, [0.0])
    elseif modelfreq[i] == 1
      push!(scales, particles[modelindex[:, i]][1].scales)
    else
      parameters = hcat(map(x -> x.params, particles[modelindex[:, i]])...)'
      push!(scales, ((maximum(parameters, 1) -
                      minimum(parameters, 1)) ./ABCsetup.scalefactor)[:])
    end
  end

  for i in 1:length(particles)
    particles[i].scales = scales[particles[i].model]
  end

  return particles
end

function show(io::IO, ABCresults::ABCrejectionresults)

  upperci = zeros(Float64, size(ABCresults.parameters, 2))
  lowerci = zeros(Float64, size(ABCresults.parameters, 2))
  parametermeans = zeros(Float64, size(ABCresults.parameters, 2))
  parametermedians = zeros(Float64, size(ABCresults.parameters, 2))

  for i in 1:size(ABCresults.parameters, 2)
    parametermeans[i] = mean(ABCresults.parameters[:, i])
    parametermedians[i] = median(ABCresults.parameters[:, i])
    (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[:, i], [0.025,0.975])
  end

  @printf("Number of simulations: %.2e\n", ABCresults.numsims)
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)

  print("Median (95% intervals):\n")
  for i in 1:length(parametermeans)
      @printf("Parameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
  end
end

function show(io::IO, ABCresults::ABCSMCresults)

  upperci = zeros(Float64, size(ABCresults.parameters, 2))
  lowerci = zeros(Float64, size(ABCresults.parameters, 2))
  parametermeans = zeros(Float64, size(ABCresults.parameters, 2))
  parametermedians = zeros(Float64, size(ABCresults.parameters, 2))

  for i in 1:size(ABCresults.parameters, 2)
    parametermeans[i] = mean(ABCresults.parameters[:, i],
    weights(ABCresults.weights))
    parametermedians[i] = median(ABCresults.parameters[:, i],
    weights(ABCresults.weights))
    (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[:, i],
    weights(ABCresults.weights),
    [0.025,0.975])
  end

  @printf("Total number of simulations: %.2e\n", sum(ABCresults.numsims))
  println("Cumulative number of simulations = $(cumsum(ABCresults.numsims))")
  @printf("Acceptance ratio: %.2e\n", ABCresults.accratio)
  println("Tolerance schedule = $(round.(ABCresults.ϵ, 2))\n")

  print("Median (95% intervals):\n")
  for i in 1:length(parametermeans)
      @printf("Parameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
  end
end

function show(io::IO, ABCresults::ABCrejectionmodelresults)

  @printf("Number of simulations: %.2e\n", ABCresults.numsims)
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)
  print("Model frequencies:\n")
  for j in 1:length(ABCresults.modelfreq)
    @printf("\tModel %d: %.2f\n", j, ABCresults.modelfreq[j])
  end

  print("\nParameters:\n\n")

  for j in 1:length(ABCresults.parameters)
    print("Model $j\n")

    upperci = zeros(Float64, size(ABCresults.parameters[j], 2))
    lowerci = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermeans = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermedians = zeros(Float64, size(ABCresults.parameters[j], 2))

    for i in 1:size(ABCresults.parameters[j], 2)
      parametermeans[i] = mean(ABCresults.parameters[j][:, i])
      parametermedians[i] = median(ABCresults.parameters[j][:, i])
      (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[j][:, i], [0.025,0.975])
    end

    print("\tMedian (95% intervals):\n")
    for i in 1:length(parametermeans)
        @printf("\tParameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
    end
  end
end

function show(io::IO, ABCresults::ABCSMCmodelresults)

  @printf("Total number of simulations: %.2e\n", sum(ABCresults.numsims))
  println("Cumulative number of simulations = $(cumsum(ABCresults.numsims))")
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)
  println("Tolerance schedule = $(round.(ABCresults.ϵ, 2))\n")

  print("Model probabilities:\n")
  for j in 1:length(ABCresults.modelprob)
    @printf("\tModel %d: %.3f\n", j, ABCresults.modelprob[j])
  end

  print("\nParameters:\n\n")

  for j in 1:length(ABCresults.parameters)
    print("Model $j\n")

    upperci = zeros(Float64, size(ABCresults.parameters[j], 2))
    lowerci = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermeans = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermedians = zeros(Float64, size(ABCresults.parameters[j], 2))

    for i in 1:size(ABCresults.parameters[j], 2)
      parametermeans[i] = mean(ABCresults.parameters[j][:, i],
      weights(ABCresults.weights[j]))
      parametermedians[i] = median(ABCresults.parameters[j][:, i],
      weights(ABCresults.weights[j]))
      (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[j][:, i],
      weights(ABCresults.weights[j]), [0.025,0.975])
    end

    print("\tMedian (95% intervals):\n")
    for i in 1:length(parametermeans)
        @printf("\tParameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
    end
  end
end

function getparticleweights(particles, ABCsetup)

  w = zeros(Float64, ABCsetup.nmodels, ABCsetup.nparticles)
  for i in 1:ABCsetup.nparticles
    w[particles[i].model, i] = particles[i].weight
  end
  weights = w ./ sum(w, 2)

  return weights, sum(w, 2)[:]
end

function priorprob(parameters::Array{Float64, 1}, prior::Prior)
  pprob = 1
  for i in 1:length(parameters)
      pprob = pprob * pdf(prior.distribution[i], parameters[i])
  end

  return pprob
end

"""
    writeoutput(results; <keyword arguments>)

Write the results of an ABC inference to a text file. For model selection algorithms a text file with the parameters of each model will be written and a text file with model probabilities.
...
## Arguments
- `dir = ""`: Directory where the text file will be written to.
- `file= ""`: Filename to write to, default depends on the type of inference.
...
"""
function writeoutput(results::ABCSMCresults; dir = "", file = "SMC-output.txt")
  distance = map(x -> x.distance, results.particles)
  wts = map(x -> x.weight, results.particles)
  nparams =  size(results.parameters)[2]

  head = map(x -> "parameter$x\t", 1:nparams)
  append!(head, ["distance\t", "weights\n"])

  out = hcat(results.parameters, distance, wts)

  f = open(joinpath(dir, file), "w")
  write(f, "## ABC SMC algorithm\n")
  write(f, "## Number of simulations: $(sum(results.numsims))\n")
  write(f, "## Acceptance ratio: $(round(results.accratio, 4))\n")
  write(f, "## Tolerance schedule : $(results.ϵ)\n")
  write(f, head...)
  writedlm(f, out)
  close(f)

end

function writeoutput(results::ABCSMCmodelresults; dir = "", file = "SMCModel-output")

  for i in 1:length(results.modelprob)
      prtcles = results.particles[map(x -> x.model, results.particles).==i]
      distance = map(x -> x.distance, prtcles)
      wts = map(x -> x.weight, prtcles)
      nparams =  size(results.parameters[i])[2]

      head = map(x -> "parameter$x\t", 1:nparams)
      append!(head, ["distance\t", "weights\n"])

      out = hcat(results.parameters[i], distance, wts)

      #write parameters to file
      f = open(joinpath(dir, "$(file)model$i.txt"), "w")
      write(f, "## ABC SMC algorithm\n")
      write(f, "## Model: $i\n")
      write(f, "## Number of simulations: $(sum(results.numsims))\n")
      write(f, "## Acceptance ratio: $(round(results.accratio, 4))\n")
      write(f, "## Tolerance schedule : $(results.ϵ)\n")
      write(f, head...)
      writedlm(f, out)
      close(f)
  end

  #write model probabilities to file
  head = ["Model\t", "Probability\n"]
  f = open(joinpath(dir, "$(file)modelprobabilities.txt"), "w")
  write(f, head...)
  writedlm(f, hcat(collect(1:length(results.modelprob)), results.modelprob))
  close(f)

end

function writeoutput(results::ABCrejectionresults; dir = "", file = "Rejection-output.txt")
  distance = map(x -> x.distance, results.particles)
  nparams =  size(results.parameters)[2]

  head = map(x -> "parameter$x\t", 1:nparams)
  append!(head, ["distance\n"])

  out = hcat(results.parameters, distance)

  f = open(joinpath(dir, file), "w")
  write(f, "## ABC Rejection algorithm\n")
  write(f, "## Number of simulations: $(results.numsims)\n")
  write(f, "## Acceptance ratio: $(round(results.accratio, 4))\n")
  write(f, head...)
  writedlm(f, out)
  close(f)

end

function writeoutput(results::ABCrejectionmodelresults; dir = "", file = "RejectionModel-output")

    for i in 1:length(results.modelfreq)
      prtcles = results.particles[map(x -> x.model, results.particles).==i]
      distance = map(x -> x.distance, prtcles)
      nparams =  size(results.parameters[i])[2]

      head = map(x -> "parameter$x\t", 1:nparams)
      append!(head, ["distance\n"])
      out = hcat(results.parameters[i], distance)

      f = open(joinpath(dir, "$(file)model$i.txt"), "w")
      write(f, "## ABC Rejection algorithm\n")
      write(f, "## Number of simulations: $(results.numsims)\n")
      write(f, "## Acceptance ratio: $(round(results.accratio, 4))\n")
      write(f, head...)
      writedlm(f, out)
      close(f)
  end

  #write model probabilities to file
  head = ["Model\t", "Probability\n"]
  f = open(joinpath(dir, "$(file)modelprobabilities.txt"), "w")
  write(f, head...)
  writedlm(f, hcat(collect(1:length(results.modelfreq)), results.modelfreq))
  close(f)

end


text:
"""
    metropolisStep!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T}, β::Float64)

Performs one metropolis step on Society soc using MoralVector x and temperature-like parameter β
"""
function metropolisStep!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K,T}, β::Float64)
    i = rand(1:N)
    j = sample(weights(soc[i, :]))

    # Sample a 'proposed' MoralVector using a MultivariateGaussian centered at the old Agent
    proposed = MoralVector(rand(MvNormal(soc[i][:], ones(5))))
    newcost  = cogcost(proposed, soc, i, j, x)

    ΔV = newcost - cogcost(soc, i, j, x)

    # transistion probability
    p_trans  =  (ΔV < 0 ? 1 : exp(-β*soc[i, j]*ΔV ) )

    # if (p = rand()) < p_trans
    if rand() < p_trans
        insertagent!(soc, proposed, i)
        return proposed
    end

    return i, p_trans
end

"""
    metropolis!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K, T}; β::Float64 = βDEF)

Performs metropolis on Society soc using MoralVector x and optional temperature-like parameter β
"""
function metropolis!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K, T}; β::Float64 = βDEF)
    # Maybe change this to a more fundamented convergence test
    # We could also argue that people have a maximum number of interactions in society
    iter  = 150*length(soc)

    for i in 1:iter
        # ignoring the MoralVector 'metropolisStep()' outputs
        metropolisStep!(soc, x, β);
    end

    return iter
end

"""
    metropolis!{N,K,T}(soc::Society{N,K,T}; β::Float64 = βDEF)

Performs metropolis on Society soc using random MoralVectors and optional temperature-like parameter β
"""
function metropolis!{N,K,T}(soc::Society{N,K,T}; β::Float64 = βDEF)
    # Maybe change this to a more fundamented convergence test
    # We could also argue that people have a maximum number of interactions in society
    iter  = 150*length(soc)

    # The society discusses a single issue, some kind of Zeitgeist
    x = MoralVector()

    for i in 1:iter
        # ignoring the MoralVector 'metropolisStep()' outputs
        metropolisStep!(soc, x, β);
    end

    return iter, x
end

"""
    discreteStep!{N,K,T}(soc::DistrustSociety{N,K,T}, x::MoralVector{K,T}; freeze = :none, verbose = false)

Performs one discrete update step on Society soc using MoralVector x

* freeze: freezes the dynamics of one of more variables (:none, :variances, :distrust, :slowtrust, :slowopinion)
"""
function discreteStep!{N,K,T}(soc::DistrustSociety{N,K,T}, x::MoralVector{K,T}; freeze = :none, verbose = false)
    i = rand(1:N)
    j = sample(weights(soc[i, :]))

    # i learns from j
    deltas = computeDeltas(soc, i, j, x)

    # update the parameters according to the 'freeze' keyword argument
    soc[i] = MoralVector(soc[i][:] + deltas[1])

    if freeze ∉ [:variances, :slowopinion]
        soc.C[i] += deltas[2]
    end

    if freeze ∉ [:distrust]
        soc.mu[i, j] += deltas[3]
    end

    if freeze ∉ [:variances, :distrust, :slowtrust]
        soc.s2[i, j] += deltas[4]
    end

    if verbose
        return (i, j, deltas...)
    end
    return i, j
end

"""
    discreteStep!{N,K,T}(soc::BasicSociety{N,K,T}, x::MoralVector{K,T}; freeze = :none, verbose = false)

Performs one discrete update step on Society soc using MoralVector x
"""
function discreteStep!{N,K,T}(soc::BasicSociety{N,K,T}, x::MoralVector{K,T}; verbose = false, freeze = :none)
    i = rand(1:N)
    j = sample(weights(soc[i, :]))

    # i learns from j
    delta = computeDeltas(soc, i, j, x)

    # update the parameters
    soc[i] = MoralVector(soc[i][:] + delta)

    if verbose
        return (i, j, delta)
    end
    return i, j
end

"""
    discreteEvol!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K, T}; freeze = :none)

Performs the discrete evolution on Society soc using MoralVector x

* freeze: freezes the dynamics of one of more variables (:none, :variances, :distrust, :slowtrust, :slowopinion)
"""
function discreteEvol!{N,K,T}(soc::Society{N,K,T}, x::MoralVector{K, T}; freeze = :none)
    # Maybe change this to a more fundamented convergence test
    # We could also argue that people have a maximum number of interactions in society
    iter  = 150*length(soc)

    for i in 1:iter
        # ignoring the indices 'discreteStep!' outputs
        discreteStep!(soc, x, freeze = freeze);
    end

    return iter
end

"""
    discreteEvol!{N,K,T}(soc::Society{N,K,T}; freeze = :none)

Performs one discrete update step on Society soc using several random MoralVectors

* freeze: freezes the dynamics of one of more variables (:none, :variances, :distrust, :slowtrust, :slowopinion)
"""
function discreteEvol!{N,K,T}(soc::Society{N,K,T}; freeze = :none)
    # Maybe change this to a more fundamented convergence test
    # We could also argue that people have a maximum number of interactions in society
    iter  = 150*length(soc)

    for i in 1:iter
        # ignoring the indices 'discreteStep!' outputs
        discreteStep!(soc, MoralVector(), freeze = freeze);
    end

    return iter
end

"""
    societyHistory(soc::DistrustSociety, N::Int, P::Int; freeze = :none, sequential::Bool = false, verbose::Bool = false)

Performs the discrete update on Society soc using a number of P MoralVectors in N*P iterations

* freeze: freezes the dynamics of one of more variables (:none, :variances, :distrust, :slowtrust, :slowopinion)
* sequential: whether the random vectors are presented sequentially, one after the other, or in random order

Returns the history of the evolution and the deltas used in the evolution. If verbose is set to true, also returns the order of the MoralVectors presented
"""
function societyHistory!(soc::DistrustSociety, N::Int, P::Int; freeze = :none, sequential::Bool = false, verbose::Bool = false)
    xs   = [MoralVector() for i in 1:P];

    if sequential && N > 1
         xi = repeat(xs, inner = N)
    else xi = rand(xs, N*P)
    end

    history = [deepcopy(soc)]
    deltas  = zeros(1, size(fieldnames(soc), 1) - 1)

    for i in 1:P*N
        deltasi = discreteStep!(soc, xi[i], freeze = freeze, verbose = true)[3:end]
        append!(history, [deepcopy(soc)])
        deltas = vcat(deltas, [deltasi...]')
    end

    if verbose
        return history, deltas, xi
    else return history, deltas
    end
end

"""
    societyHistory(soc::BasicSociety, N::Int, P::Int; issues = :none, freeze = :none, sequential::Bool = false, verbose::Bool = false)

Performs the discrete update on Society soc using a number of P MoralVectors in N*P iterations

* freeze: freezes the dynamics of one of more variables (:none, :variances, :distrust, :slowtrust, :slowopinion)
* sequential: whether the random vectors are presented sequentially, one after the other, or in random order

Returns the history of the evolution and the deltas used in the evolution. If verbose is set to true, also returns the order of the MoralVectors presented
"""
function societyHistory!(soc::BasicSociety, N::Int, P::Int; issues = :none, sequential::Bool = false, verbose::Bool = false)
    if typeof(issues) == Matrix{Float64} && size(issues, 2) == size(soc, 2)
         xs = [MoralVector(issues[i, :]) for i in rand(1:size(issues, 1), P)]
    else xs = [MoralVector() for i in 1:P];
    end

    if sequential && N > 1
         xi = repeat(xs, inner = N)
    else xi = rand(xs, N*P)
    end

    history = [deepcopy(soc)]
    deltas  = zeros(1, size(fieldnames(soc), 1) - 3)

    for i in 1:P*N
        dW = discreteStep!(soc, xi[i], freeze = :none, verbose = true)[3]
        append!(history, [deepcopy(soc)])
        deltas = vcat(deltas, [dW]')
    end

    if verbose
        return history, deltas, xi
    else return history, deltas
    end
end

"""
    computeDeltas{K, T}(agi::MoralVector{K, T}, agj::MoralVector{K, T}, x::MoralVector{K,T}, C::Matrix{Float64}, μ::Float64, s2::Float64)

Compute the evolution delta for each of the variables of an agent i listening agent j about issue x given all other parameters
"""
function computeDeltas{K, T}(agi::MoralVector{K, T}, agj::MoralVector{K, T}, x::MoralVector{K,T}, C::Matrix{Float64}, μ::Float64, s2::Float64)
    σ  = sign(agj ⋅ x)
    h  = agi ⋅ x
    Cx = C * x[:]
    γ  = sqrt(x[:]' * C * x[:])

    Fw, Fϵ = modfunc(h*σ/γ, μ/sqrt(1 + s2))

    return  (Fw * σ * Cx / γ,
            -Fw * (Fw + h*σ/γ) * Cx * Cx' / γ^2,
             Fϵ * s2 / sqrt(1 + s2),
            -Fϵ * (Fϵ + μ/sqrt(1 + s2)) * s2^2 / (1 + s2) )

end

"""
    computeDeltas{N, K, T}(soc::DistrustSociety{N,K,T}, i::Int, j::Int, x::MoralVector{K,T})

Compute the evolution delta for each of the variables of an agent i from society soc listening agent j about issue x
"""
function computeDeltas{N, K, T}(soc::DistrustSociety{N,K,T}, i::Int, j::Int, x::MoralVector{K,T})
    σ  = sign(soc[j] ⋅ x)
    h  = soc[i] ⋅ x
    C  = soc.C[i]
    Cx = C * x[:]
    γ  = gammasoc(soc, i, x)
    μ  = soc.mu[i, j]
    s2 = soc.s2[i, j]

    hσ1γ   = h*σ/γ
    μ1sqs2 = μ/sqrt(1 + s2)

    Fw, Fϵ = modfunc(hσ1γ, μ1sqs2)

    return  (Fw * σ * Cx / γ,
            -Fw * (Fw + hσ1γ) * Cx * Cx' / γ^2,
             Fϵ * s2 / sqrt(1 + s2),
            -Fϵ * (Fϵ + μ1sqs2) * s2^2 / (1 + s2) )

end

"""
    computeDeltas{N, K, T}(soc::BasicSociety{N,K,T}, i::Int, j::Int, x::MoralVector{K,T})

Compute the evolution delta for each of the variables of an agent i from society soc listening agent j about issue x
"""
function computeDeltas{N, K, T}(soc::BasicSociety{N,K,T}, i::Int, j::Int, x::MoralVector{K,T})
    σj  = sign(soc[j] ⋅ x)
    hi  = soc[i] ⋅ x
    γi  = gammasoc(soc, i) #, x)
    ϵij = epssoc(soc, i, j)

    hσ1γ = hi*σj/γi
    phiw = phi(hσ1γ)

    Z    = ϵij + phiw - 2*ϵij*phiw

    ret = (1 - 2ϵij) * G(hσ1γ) * σj * γi * x[:] / Z
    if any(isnan.(ret))
        return 0.
    end
    return ret

end

"""
    pickIssue{N, K, T}(soc::DistrustSociety{N, K, T}, i::Int; err::Float64 = 1e-4, niter::Int = 500, α::Float64 = 0.1, verbose = false)

Pick a specific issue x trying to maximize the change in trust of agent-listener i
"""
function pickIssue{N, K, T}(soc::DistrustSociety{N, K, T}, i::Int;
                  err::Float64 = 1e-4, niter::Int = 500, α::Float64 = 0.1, verbose = false)

    w = soc[i]
    C = soc.C[i]

    x = - sqrt(K)*w

    t = 1.
    n = 0.

    if verbose xs = zeros(K+1, niter) end

    @inbounds while t > err && n < niter

        h = w ⋅ x
        Γ = x ⋅ (C * x)/K

        xnew = α .* x + (1 - α) .* (2 .* (C * x) ./ Γ .- K .* w ./ h)

        t = norm(xnew - x)
        n += 1.

        if verbose
            xs[1:K, Int(n)] = xnew
            xs[K+1, Int(n)] = t
        end

        x = sqrt(K) * normalize(xnew[:])
    end

    if verbose
        return xs[:, 1:Int(n)], n
    end

    return x
end

"""
    pickIssue(w::Vector{Float64}, C::Matrix{Float64}; err::Float64 = 1e-4, niter::Int = 500, α::Float64 = 0.1, verbose = false)

Pick a specific issue x trying to maximize the change in trust of agent-listener w which has covariance-doubts C
"""
function pickIssue(w::Vector{Float64}, C::Matrix{Float64};
                  err::Float64 = 1e-4, niter::Int = 500, α::Float64 = 0.1, verbose = false)
    t = 1.
    n = 0.

    K = length(w)
    x = - w

    if ! issymmetric(C) error("Fix entries of C. It is not symmetric") end

    xnew = zeros(x)

    if verbose xs = zeros(K+1, niter) end

    @inbounds while t > err && n < niter

        h = w ⋅ x
        Γ = x ⋅ (C * x)/K

        xnew = α .* x + (1 - α) .* (2 .* (C * x) ./ Γ .- K .* w ./ h)

        t = norm(xnew - x)
        n += 1.

        if verbose
            xs[1:K, Int(n)] = xnew
            xs[K+1, Int(n)] = t
        end

        x = sqrt(K) * normalize(xnew[:])
    end

    if verbose
        return xs, n
    end

    return x, n
end


text:
export
   execute,
   execute_single,
   execute_single_measured,
   execute_all,
   execute_all_measured

"""
    execute(qwd, initstate, runtime[, all=false, measure=false])

Run proper execution function depending on given keywords. `all` and `measure`
keywords defaults to `false`. In case of `all` being true, all intermidiate states
are returned. Note that for `all` equal to `true`, the model
specified by `qwd` needs to be disrete. For `measure` equal to true, measurement
distributions are returned.
"""
function execute(qwd::QWDynamics,
                 initstate,
                 runtime::Real;
                 all::Bool = false,
                 measure::Bool = false)
   if !all && !measure
      execute_single(qwd, initstate, runtime)
   elseif !all && measure
      execute_single_measured(qwd, initstate, runtime)
   elseif all && !measure
      execute_all(qwd, initstate, runtime)
   else
      execute_all_measured(qwd, initstate, runtime)
   end
end

"""
    execute_single(qwd, initstate, runtime)

Evolve `initstate` according to QWDynamics `qwd` for time `runtime`.
`runtime` needs to be nonnegative.
"""
function execute_single(qwd::QWDynamics{<:QWModelDiscr},
                        initstate,
                        runtime::Int)
   @assert runtime>=0 "Argument 'runtime' needs to be nonnegative"

   state = initstate
   for t=1:runtime
      state = evolve(qwd, state)
   end

   state
end

function execute_single(qwd::QWDynamics{<:QWModelCont},
                        initstate,
                        runtime::Real)
   @assert runtime>=0 "Argument 'runtime' needs to be nonnegative"

   evolve(qwd, initstate, runtime)
end

"""
    execute_single_measured(qwd, initstate, runtime)

Evolve `initstate` acording to `qwd` for time `runtime` and measure
it in the end. `runtime` needs to be nonnegative.
"""
function execute_single_measured(qwd::QWDynamics,
                                 initstate,
                                 runtime::Real)
   measure(qwd, execute_single(qwd, initstate, runtime))
end

"""
    execute_all(qwd, initstate, runtime)

Returns list of all states including the `initstate` according to `qwd` for
times from 0 to `runtime`. `runtime` needs to be nonnegative. Wroks only for
discrete-time evolution.
"""
function execute_all(qwd::QWDynamics{<:QWModelDiscr},
                     initstate::S,
                     runtime::Int) where S
   @assert runtime>=0 "Argument 'runtime' needs to be nonnegative"

   result = Vector{S}([initstate])
   state = initstate
   for t=1:runtime
      state = evolve(qwd, state)
      push!(result, state)
   end

   result
end

"""
    execute_all_measured(qwd, initstate, runtime)

Evolve `initstate` acording to `qwd` for time `runtime`. Returns matrix of type
`Matrix{Float64}` for which `i`-th column is the probability distribution obtained
 from the measurement in (`i-1`)-th step. `runtime` needs to be nonnegative.
 Wroks only for discrete-time evolution.
"""
function execute_all_measured(qwd::QWDynamics{<:QWModelDiscr},
                              initstate,
                              runtime::Int)
   @assert runtime>=0 "Argument 'runtime' needs to be nonnegative"

   result = zeros(Float64, (nv(graph(qwd)), runtime+1)) # +1 to include 0

   state = initstate
   result[:,1] = measure(qwd, state)
   for t=1:runtime
      state = evolve(qwd, state)
      result[:,t+1] = measure(qwd, state) #evolution starts with t=0
   end

   result
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s = get_simulation("quadruped", "flat_2D_lc", "flat")
model = s.model
env = s.env

# ## Run policy
function run_policy(s::Simulation; H_sim::Int = 2000, verbose = false,
		qinit::AbstractVector = zeros(s.model.nq))
	model = s.model
	env = s.env

	ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(model, env,
	    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
	    load_type = :split_traj_alt))

	H = ref_traj.H
	h = ref_traj.h
	N_sample = 5
	H_mpc = 10
	h_sim = h / N_sample
	H_sim = H_sim
	κ_mpc = 1.0e-4

	obj = TrackingObjective(model, env, H_mpc,
	    q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.75 * ones(model.nq-3)]) for t = 1:H_mpc],
	    u = [Diagonal(3e-2 * ones(model.nu)) for t = 1:H_mpc],
		γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
	    b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc])

	p = ci_mpc_policy(ref_traj, s, obj,
	    H_mpc = H_mpc,
	    N_sample = N_sample,
	    κ_mpc = κ_mpc,
		mode = :configuration,
	    n_opts = NewtonOptions(
			r_tol = 3e-4,
			max_iter = 5,
			max_time = ref_traj.h, # HARD REAL TIME
			),
	    mpc_opts = CIMPCOptions(
	        ),
		ip_opts = InteriorPointOptions(
			max_iter = 100,
			verbose = false,
			r_tol = 1.0e-4,
			κ_tol = 1.0e-4,
			diff_sol = true,
			solver = :empty_solver,
			),
	    )

	q1_ref = copy(qinit)
	q0_ref = copy(qinit)
	q1_sim = ContactImplicitMPC.SVector{model.nq}(q1_ref)
	q0_sim = ContactImplicitMPC.SVector{model.nq}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
	@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

	v1_sim = (ref_traj.q[2] - ref_traj.q[1]) ./ h

	sim = simulator(s, H_sim, h=h_sim, policy=p);

	status = simulate!(sim, q1_sim, v1_sim);
	return deepcopy(sim.traj)
end

# ## Collect runs
function collect_runs(s::Simulation; n::Int = 1, H_sim::Int = 2000, verbose::Bool = false)
	nq = s.model.nq
	trajs = []
	ContactImplicitMPC.Random.seed!(100)
	conf_max = [0.05, 0.8, 0.8, 0.8, +0.2,  0.10]
	conf_min = [0.00, 0.6, 0.6, 0.6, -0.2, -0.30]
	conf_Δ = conf_max .- conf_min
	for i = 1:n
		verbose && println("sample = $i/$n")
		conf = conf_min .+ conf_Δ .* ContactImplicitMPC.rand(length(conf_min))
		conf[end] = max(conf[end], 0.0)
		qinit = initial_configuration(s.model, conf...)
		traj = run_policy(s, H_sim = H_sim, qinit = qinit, verbose = verbose)
		push!(trajs, traj)
	end
	return trajs
end

# ## Initial configurations
function initial_configuration(model::Quadruped, θ0, θ1, θ2, θ3, x, Δz)
    q = zeros(model.nq)
	q[1] = x
    q[3] = pi / 2.0
    q[4] = -θ1
    q[5] = θ2

    q[8] = -θ1
    q[9] = θ2

    q[2] = model.l_thigh1 * cos(q[4]) + model.l_calf1 * cos(q[5])

    q[10] = -θ3
    q[11] = acos((q[2] - model.l_thigh2 * cos(q[10])) / model.l_calf2)

    q[6] = -θ3
    q[7] = acos((q[2] - model.l_thigh2 * cos(q[6])) / model.l_calf2)

	q[2] += Δz
	q[3] += θ0
    return q
end

# ## Visualize runs
function visualize_runs!(vis::ContactImplicitMPC.Visualizer, model::Model, trajs::AbstractVector;
		sample=max(1, Int(floor(trajs[1].H / 100))), h=trajs[1].h*sample,  α=1.0,
		anim::ContactImplicitMPC.MeshCat.Animation=ContactImplicitMPC.MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=ContactImplicitMPC.model_name(model))

	for (i, traj) in enumerate(trajs)
		name_i = Symbol(string(name) * string(i))
		anim = visualize_meshrobot!(vis, model, traj, sample = sample, h = h, α = α, anim = anim, name = name_i)
	end
	return anim
end

# ## Experiments
H_sim = 1000
trajs = collect_runs(s, n = 100, H_sim = H_sim, verbose = true)

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
anim = visualize_runs!(vis, model, trajs, α=0.3, sample=5)

# ## High drop
conf = [-0.05, 0.6, 0.6, 0.6, 0.0, 0.30]
qinit = initial_configuration(s.model, conf...)
traj = run_policy(s, H_sim = 300, qinit=qinit)
anim = visualize_meshrobot!(vis, model, traj, name=:drop, anim=anim, sample=5, α=1.0)