text:
__precompile__()
module BittrexAPI

# module dependencies
import HTTP: get
import JSON: Parser.parse
import Nettle: hexdigest

# exports
export Bittrex
# public functions
export getmarkets,getcurrencies,getticker
export getmarketsummary,getmarketsummaries
export getorderbook,getmarkethistory
# private functions
export buylimit,selllimit,cancel
export getopenorders,getorder
export getbalances,getbalance
export getdepositaddress,withdraw
export getorderhistory,getwithdrawalhistory,getdeposithistory
# additional functions
export gettickshistory

# endpoint formats
# v1.1: "https://bittrex.com/api/$version/$apitype/$method?$parameter=$value
# v2.0: "https://bittrex.com/api/$version/$type/$section/$method?$parameter=$value

URLROOT = "https://bittrex.com/api/"

# class defining the Bittrex API
struct Bittrex
    apikey   ::String
    apisecret::String
    version  ::String
end

# constructor for public API
function Bittrex(version = "v1.1")
    return Bittrex("","",version)
end

# constructor for private API
function Bittrex(apikey, apisecret)
    return Bittrex(apikey,apisecret, "v1.1")
end

# Public API

"""
    getmarkets(api::Bittrex)

Used to get the open and available trading markets at Bittrex along with other meta data.

# Example output

```jldoctest
julia> getmarkets()[1]
Dict{String,Any} with 11 entries:
  "MarketCurrencyLong" => "Litecoin"
  "IsSponsored"        => nothing
  "MarketCurrency"     => "LTC"
  "BaseCurrency"       => "BTC"
  "LogoUrl"            => "https://bittrex.com/..."
  "BaseCurrencyLong"   => "Bitcoin"
  "IsActive"           => true
  "MarketName"         => "BTC-LTC"
  "MinTradeSize"       => 0.0278418
  "Notice"             => nothing
  "Created"            => "2014-02-13T00:00:00"
```
"""
function getmarkets(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getmarkets")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/markets/getmarkets")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = endpoint

    return apirequest(url)
end

"""
    getcurrencies(api::Bittrex)

Used to get all supported currencies at Bittrex along with other meta data.

# Example output

```jldoctest
julia> getcurrencies()[1]
Dict{String,Any} with 8 entries:
  "Currency"        => "BTC"
  "IsActive"        => true
  "BaseAddress"     => "1N52wHoVR79PMDishab2XmRHsbekCdGquK"
  "CurrencyLong"    => "Bitcoin"
  "Notice"          => nothing
  "TxFee"           => 0.001
  "CoinType"        => "BITCOIN"
  "MinConfirmation" => 2
```
"""
function getcurrencies(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getcurrencies")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/currencies/getcurrencies")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = endpoint

    return apirequest(url)
end

"""
    getticker(api::Bittrex, market::String)

Used to get the current tick values for a market.

# Example output

```jldoctest
julia> getticker(api, "BTC-LTC")
Dict{String,Any} with 3 entries:
  "Bid"  => 0.00949003
  "Last" => 0.00950316
  "Ask"  => 0.00950315
```
"""
function getticker(api::Bittrex, market::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getticker")
        query = ["market"=>market]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getlatesttick")
        query = ["marketname"=>market]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = urlquery(endpoint, query)

    return apirequest(url)
end

"""
    getmarketsummaries(api::Bittrex)

Used to get the last 24 hour summary of all active exchanges.

# Example output

```jldoctest
julia> getmarketsummaries(api)[1]
Dict{String,Any} with 13 entries:
  "Bid"            => 4.621e-5
  "TimeStamp"      => "2017-11-26T15:33:37.473"
  "Last"           => 4.645e-5
  "Low"            => 4.3e-5
  "OpenBuyOrders"  => 239
  "OpenSellOrders" => 2870
  "High"           => 4.796e-5
  "MarketName"     => "BTC-1ST"
  "Ask"            => 4.643e-5
  "BaseVolume"     => 41.8878
  "Volume"         => 9.27222e5
  "PrevDay"        => 4.72e-5
  "Created"        => "2017-06-06T01:22:35.727"
```
"""
function getmarketsummaries(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getmarketsummaries")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getmarketsummaries")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = endpoint

    return apirequest(url)
end

"""
    getmarketsummary(api::Bittrex, market::String)

Used to get the last 24 hour summary for a market.

# Example output

```jldoctest
julia> getmarketsummary(api, "BTC-LTC")
Dict{String,Any} with 13 entries:
  "Bid"            => 0.00945801
  "TimeStamp"      => "2017-11-26T15:27:31.04"
  "Last"           => 0.00945801
  "Low"            => 0.00921
  "OpenBuyOrders"  => 4310
  "OpenSellOrders" => 7074
  "High"           => 0.010238
  "MarketName"     => "BTC-LTC"
  "Ask"            => 0.00948002
  "BaseVolume"     => 2285.16
  "Volume"         => 2.33973e5
  "PrevDay"        => 0.00978736
  "Created"        => "2014-02-13T00:00:00"
```
"""
function getmarketsummary(api::Bittrex, market::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getmarketsummary")
        query = ["market"=>market]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getmarketsummary")
        query = ["marketname"=>market]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = urlquery(endpoint, query)

    return apirequest(url)
end

"""
    getorderbook(api::Bittrex, market::String, [booktype::String = "both"])

Used to retrieve the orderbook for a given market.
`booktype` can be any either "buy", "sell" or "both"

# Example output

```jldoctest
julia> getorderbook(api, "BTC-LTC")
Dict{String,Any} with 2 entries:
  "sell" => Any[Dict{String,Any}(Pair{String,Any}("Quantity", 129.906),Pair{String,Anly}("Rate", 0.00945511)), ...]
  "buy"  => Any[Dict{String,Any}(Pair{String,Any}("Quantity", 1.30807),Pair{String,Anly}("Rate", 0.00945508)), ...]
```
"""
function getorderbook(api::Bittrex, market::String, booktype::String = "both")
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getorderbook")
        query = ["market"=>market,"type"=>booktype]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getmarketorderbook")
        query = ["marketname"=>market,"type"=>booktype]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = urlquery(endpoint, query)

    return apirequest(url)
end

"""
    getmarkethistory(api::Bittrex, market::String)

Used to retrieve the latest trades that have occured for a specific market.

# Example output

```jldoctest
julia> getmarkethistory(api, "BTC-LTC")[1]
Dict{String,Any} with 7 entries:
  "Id"        => 105218942
  "Quantity"  => 0.250755
  "FillType"  => "FILL"
  "TimeStamp" => "2017-11-26T15:39:43.637"
  "Total"     => 0.0023681
  "Price"     => 0.00944392
  "OrderType" => "SELL"
```
"""
function getmarkethistory(api::Bittrex, market::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/public/getmarkethistory")
        query = ["market"=>market]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getmarkethistory")
        query = ["marketname"=>market]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = urlquery(endpoint, query)

    return apirequest(url)
end

# End of Public API

# Market API

"""
    buylimit(api::Bittrex, market::String, quantity, rate)

Used to place a buy order in a specific market.
"""
function buylimit(api::Bittrex, market::String, quantity, rate)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/market/buylimit")
        query = ["market"=>market,"quantity"=>quantity,"rate"=>rate]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/buylimit")
        query = ["marketname"=>market,"quantity"=>quantity,"rate"=>rate]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    unshift!(query,"apikey"=>api.apikey,"nonce"=>nonce)
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    selllimit(api::Bittrex, market::String, quantity, rate)

Used to place a sell order in a specific market.
"""
function selllimit(api::Bittrex, market::String, quantity, rate)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/market/selllimit")
        query = ["market"=>market,"quantity"=>quantity,"rate"=>rate]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/selllimit")
        query = ["marketname"=>market,"quantity"=>quantity,"rate"=>rate]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    unshift!(query,"apikey"=>api.apikey,"nonce"=>nonce)
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    cancel(api::Bittrex, uuid::String)

Used to cancel a buy or sell order by specifying the uuid.
"""
function cancel(api::Bittrex, uuid::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/market/cancel")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/cancel")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce,"uuid"=>uuid]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getopenorders(api::Bittrex, [market::String])

Get all orders that you currently have opened. A specific market can be requested.
"""
function getopenorders(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/market/getopenorders")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getopenorders")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

function getopenorders(api::Bittrex, market::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/market/getopenorders")
        query = ["market"=>market]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/market/getopenorders")
        query = ["marketname"=>market]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    unshift!(query,"apikey"=>api.apikey,"nonce"=>nonce)
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

# End of Market API

# Account API

"""
    getbalances(api::Bittrex)

Used to retrieve all balances from your account.
"""
function getbalances(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getbalances")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getbalances")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getbalance(api::Bittrex, currency::String)

Used to retrieve the balance from your account for a specific currency.
"""
function getbalance(api::Bittrex, currency::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getbalance")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getbalance")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce,"currency"=>currency]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getdepositaddress(api::Bittrex, currency::String)

Used to retrieve or generate an address for a specific currency.
If one does not exist, the call will fail and return ADDRESS_GENERATING until one is available.
"""
function getdepositaddress(api::Bittrex, currency::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getdepositaddress")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getdepositaddress")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce,"currency"=>currency]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    withdraw(api::Bittrex, currency::String, quantity, address::String)

Used to withdraw funds from your account.
"""
function withdraw(api::Bittrex, currency::String, quantity, address::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/withdraw")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/withdraw")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["currency"=>currency,"quantity"=>quantity,"address"=>address]
    unshift!(query,"apikey"=>api.apikey,"nonce"=>nonce)
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getorder(api::Bittrex, uuid::String)

Used to retrieve a single order by uuid.
"""
function getorder(api::Bittrex, uuid::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getorder")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getorder")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce,"uuid"=>uuid]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getorderhistory(api::Bittrex, [market::String])

Used to retrieve your order history. A specific market can be requested.
"""
function getorderhistory(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getorderhistory")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getorderhistory")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

function getorderhistory(api::Bittrex, market::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getorderhistory")
        query = ["market"=>market]
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getorderhistory")
        query = ["marketname"=>market]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    unshift!(query,"apikey"=>api.apikey,"nonce"=>nonce)
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getwithdrawalhistory(api::Bittrex, [currency::String])

Used to retrieve your withdrawal history. A specific currency can be requested.
"""
function getwithdrawalhistory(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getwithdrawalhistory")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getwithdrawalhistory")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

function getwithdrawalhistory(api::Bittrex, currency::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getwithdrawalhistory")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getwithdrawalhistory")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce,"currency"=>currency]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

"""
    getdeposithistory(api::Bittrex, [currency::String])

Used to retrieve your deposit history. A specific currency can be requested.
"""
function getdeposithistory(api::Bittrex)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getdeposithistory")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getdeposithistory")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

function getdeposithistory(api::Bittrex, currency::String)
    if api.version == "v1.1"
        endpoint = string(URLROOT,api.version,"/account/getdeposithistory")
    elseif api.version == "v2.0"
        endpoint = string(URLROOT,api.version,"/pub/account/getdeposithistory")
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    nonce   = string(Int(round(time()*1000,0)))
    query   = ["apikey"=>api.apikey,"nonce"=>nonce,"currency"=>currency]
    url     = urlquery(endpoint, query)
    apisign = hexdigest("sha512",api.apisecret,url)

    return apirequest(url,["apisign"=>apisign])
end

# End of Account API

# nice function available in API v2.0 only

"""
    gettickshistory(api::Bittrex, market::String [, interval = "oneMin"])

Used to retrieve the latest ticks that have occured for a specific market and tick interval.
Valid tickintervals = ["oneMin","fiveMin","thirtyMin","hour","day"]

# Example output

```jldoctest
julia> gettickshistory()[1]
Dict{String,Any} with 7 entries:
  "T"  => "2017-10-17T16:00:00"
  "C"  => 0.01076
  "V"  => 2127.41
  "L"  => 0.01074
  "BV" => 22.971
  "H"  => 0.010831
  "O"  => 0.0108139
```
"""
function gettickshistory(api::Bittrex, market::String, interval = "oneMin")
    # API v2.0 only
    if api.version == "v1.1" || api.version == "v2.0"
        endpoint = string(URLROOT,"/v2.0/pub/market/getticks")
        query = ["marketname"=>market,"tickinterval"=>interval]
    else
        error("supported versions are 'v1.1' and 'v2.0'")
    end

    url = urlquery(endpoint, query)

    return apirequest(url)
end

### Helper functions ###

"""Perform a request and parse the JSON resonse into a `Dict` or `Array{Dict}`"""
function apirequest(url::String, headers = [])
    try
        resp = get(url, headers)
        body = parse(String(resp.body))
        if body["success"]
            if length(body["result"]) == 1
                return body["result"][1]
            else
                return body["result"]
            end
        else
            error(body["message"])
        end
    catch exception
        info("caught ",typeof(exception)," exception\n",exception)
    end
end

"""Form a URL `string` by merging an endpoint and query paramenters"""
function urlquery(endpoint::String, query)
    endpoint *= "?"
    for i=1:length(query)
        endpoint *= query[i][1] * "=" * string(query[i][2])
        if i < length(query)
            endpoint *= "&"
        end
    end
    return endpoint
end

end # end of module


text:
module Basic

using Flux
using Requires
using Requires: @init

using ..Transformers: Abstract3DTensor, Container, batchedmul

export PositionEmbedding, Embed, getmask,
    Vocabulary, gather, encode, decode
export OneHotArray, indices2onehot, onehot2indices,
    onehotarray, onehot, tofloat
export Transformer, TransformerDecoder

export NNTopo, @nntopo_str, @nntopo
export Stack, show_stackfunc, stack

export @toNd, Positionwise

export logkldivergence, logcrossentropy

include("./extend3d.jl")
include("./stack/topology.jl")
include("./stack/stack.jl")
include("./embed/Embed.jl")
include("./mh_atten.jl")
include("./transformer.jl")
include("./loss.jl")

end

ext: jl


text:
## A Software Package for Solving Dynamic Rational Inattention Problems (D.R.I.P.)
## Copyright © 2020 Hassan Afrouzi, Choongryul Yang and Miguel Acosta
## When using this package for your work, cite Afrouzi and Yang (2019) for reference.
## Julia code written by Miguel Acosta / Modified by Hassan Afrouzi

using LinearAlgebra

## A General Structure for D.R.I.P.
mutable struct Drip
    ω; β; A; Q; H;                    # primitives
    K; Y; Σ_z; Σ_p; Σ_1; Ω;           # solution
    Drip()          = new();
    Drip(ω,β,A,Q,H) = new(ω,β,A,Q,H);
end

## A General Structure for Transition dynamics of Rational Inattention Problems (T.R.I.P.)
struct Trip
    P::Drip;                # problem and solution in steady state
    T::Int;                 # length of T.R.I.P.
    Σ_1s; Σ_ps; Ωs;         # priors, posteriors and benefit matrices
    Ds;                     # eigenvalues of Σ_t^(0.5)Ω_tΣ_t^(0.5)
end

## A structure for additional signals for treatment
## S₀ = L'x⃗₀+z₀, z₀~N(0,Σ_z)
struct Signal
    L; Σ_z;
end

## A Structure for the impulse responses of D.R.I.P. for the steady state posterior
struct Dripirfs
    T     :: Int            # length of IRFs
    x     :: Array{Float64} # IRFs of the fundamental shocks
    x_hat :: Array{Float64} # IRFs of beliefs
    a     :: Array{Float64} # IRFs of actions
end

## solve_Drip: this function solves for the steady state of of the D.R.I.P.
## Inputs
# ω      : Cost of information
# β      : Discount factor
# A      : Transition matrix: x=Ax+Qu
# Q      : Std. Dev. matrix: x=Ax+Qu
# H      : Mapping of shocks to actions: v=-0.5(a'-x'H)(a-H'x)

## Outputs: a Drip structure with
# Σ_1    : Steady-state Prior uncertainty
# Σ_p    : Steady-state posterior uncertainty
# Λ      : Shadow matrix on the no-forgetting constraint -- Λ*(Σ-Σ_p) = 0
# Ω      : Dynamic benefit matrix
# Y      : Weight vector for evolution of actions
# Σ_z    : Covariance matrix of the rational inattention error
# K      : Kalman gain matrix

function solve_drip(ω,β,A,Q,H;              # primitives of the D.R.I.P.
                    fcap::Bool = false,     # optional: if true then solves the problem with fixed capacity κ = ω.
                    Ω0         = H*H',      # optional: initial guess for steady state information matrix
                    Σ0         = A*A'+Q*Q', # optional: initial guess for steady state prior
                    w          = 1,         # optional: updating weight in iteration
                    tol        = 1e-4,      # optional: tolerance level for convergence
                    maxit      = 10000)     # optional: maximum number of iterations

    ## initialize
    (n,m) = length(size(H)) == 2 ? size(H) : (size(H,1),1);
    # n: dimension of state, m: number of actions
    eye   = Matrix{Float64}(I,n,n);
    err   = 1
    iter  = 0
    SqRΣ  = sqrt(Σ0)
    Ω_c   = H*H';
    Σ1    = Matrix{Float64}(I,n,n)
    Σ_p   = Matrix{Float64}(I,n,n)
    Λ     = Matrix{Float64}(I,n,n)
    κ     = ω
    # iterate
    while (err > tol) & (iter < maxit)
        D, U    = eigen(SqRΣ*Ω0*SqRΣ);
        D       = diagm(getreal(D));
        U       = getreal(U);

        if fcap == true
            ω = (2^(2*κ)/det(max.(ω*eye,D)))^(-1/n);
        end
        Λ       = U*(max.(ω*eye-D,0.0))*U';
        Σ_p     = getreal(ω*SqRΣ*U/(max.(D,ω*eye))*U'*SqRΣ);

        Σ1      = A*Σ_p*A' + Q*Q'
        err     = norm(Σ1 - Σ0,2)/norm(Σ0,2);

        Σ0      = w*Σ1 + (1-w)*Σ0

        SqRΣ    = getreal(sqrt(Σ0));
        invSqRΣ = pinv(SqRΣ);

        Ω0      = w*(Ω_c + β*A'*invSqRΣ*(ω*eye - Λ)*invSqRΣ*A)+(1-w)*Ω0;
        Ω0      = (abs.(Ω0).>1e-10).*Ω0;

        iter   += 1
    end

    if iter == maxit
        print("RI Code hit maxit\n")
    end
    P      = Drip(ω,β,A,Q,H);
    Σ1     = (abs.(Σ1).>1e-10).*Σ1 + diagm(abs.(diag(Σ1)).<= 1e-10)*1e-8
    inv_Σ1 = pinv(Σ1) ;
    P.Σ_1  = Σ1;
    P.Σ_p  = Σ_p;
    P.Y    = (eye - Σ_p*inv_Σ1)'*H;
    P.Σ_z  = H'*(Σ_p - Σ_p*inv_Σ1*Σ_p)*H;
    P.K    = Σ1*P.Y*pinv(P.Y'*Σ1*P.Y + P.Σ_z);
    P.Ω    = Ω0;

    return(P)
end

function solve_drip(P   ::Drip;             # D.R.I.P. to be solved
                    fcap::Bool = false,     # optional: if true then solves the problem with fixed capacity κ = ω.
                    Ω0         = P.H*P.H',  # optional: initial guess for steady state information matrix
                    Σ0         = P.A*P.A'+P.Q*P.Q', # optional: initial guess for steady state prior
                    w          = 1,         # optional: updating weight in iteration
                    tol        = 1e-4,      # optional: tolerance level for convergence
                    maxit      = 10000)     # optional: maximum number of iterations
    P = solve_drip(P.ω,P.β,P.A,P.Q,P.H;
                    fcap  = fcap,
                    Ω0    = Ω0,
                    Σ0    = Σ0,
                    w     = w,
                    tol   = tol,
                    maxit = maxit);
    return(P)
end

## solve_Trip: this function solves for the transition dynamics given an initial prior covariance matrix
function solve_trip(Ss::Drip,             # D.R.I.P. steady state
                    Σ0::Array{Float64,2}; # initial prior matrix
                    T     = 100,          # optional: time until convergence to steady state
                    tol   = 1e-4,         # optional: tolerance for convergence
                    maxit = 1000          # optional: max iterations
                    )
    ## Initialize
    Ωs  = repeat(Ss.Ω, inner = [1,1,T]);
    Ωsp = Ωs;

    Σ_1s  = repeat(Ss.Σ_1, inner = [1,1,T]); Σ_1s[:,:,1] = Σ0;
    Σ_1sp = Σ_1s;
    (n,m) = length(size(Ss.H)) == 2 ? size(Ss.H) : (size(H,1),1);
    # n: dimension of state, m: number of actions

    Σ_ps= repeat(Ss.Σ_p, inner = [1,1,T]); # initialize posteriors
    Ds  = zeros(n,T);                          # initialize eigenvalues

    iter = 0;
    err  = 1;
    eye  = Matrix(I,n,n);
    while (err > tol) & (iter <= maxit)
        ## Given Ωs, find Sigmas using the law of motion for priors
        for i in 1:1:T-1
            SqSigma        = real.(sqrt(Σ_1s[:,:,i]));
            D, U           = eigen(SqSigma*Ωs[:,:,i]*SqSigma);
            Ds[:,i]        = real.(D);
            U              = real.(U);
            D              = diagm(Ds[:,i]);
            Σ_ps[:,:,i]    = Ss.ω*SqSigma*U/(max.(D,Ss.ω*eye))*U'*SqSigma;
            Σ_1sp[:,:,i+1] = Ss.Q*Ss.Q'+Ss.A*Σ_ps[:,:,i]*Ss.A';
        end
        # Given Σ_1s, Find Omegas using the Euler equation
        for i = T-1:-1:1
            SqSigma      = sqrt(Σ_1sp[:,:,i+1]);
            invSqSigma   = real.(pinv(SqSigma));
            D, U         = eigen(SqSigma*Ωs[:,:,i+1]*SqSigma);
            D            = diagm(getreal(D));
            U            = real.(U);
            Ωsp[:,:,i]   = Ss.H*Ss.H'+Ss.β*Ss.A'*invSqSigma*U*min.(D,Ss.ω*eye)*U'*invSqSigma*Ss.A;
        end
        err    = norm(Σ_1sp-Σ_1s)/norm(Σ_1s);
        Σ_1s   = real.(Σ_1sp);
        Ωs     = real.(Ωsp);
        iter  += 1;
    end
    # Throw error if T was too short for convergence to steady state
    con_err = norm(Ss.Q*Ss.Q'+Ss.A*Σ_ps[:,:,end-1]*Ss.A'-Ss.Σ_1)/norm(Ss.Σ_1);
    if (con_err  > tol)
        error("T was too short for convergence of Trip. Try larger T.");
    end
    # Finally, store the eigenvalues of steady state
    SqSigma      = real.(sqrt(Σ_1s[:,:,end]));
    D, U         = eigen(SqSigma*Ωs[:,:,end]*SqSigma);
    Ds[:,end]    = real.(D);
    # return the Trip structure
    return(Trip(Ss,T,Σ_1s,Σ_ps,Ωs,Ds))
end

## solve_Trip: this function solves for the transition dynamics given an initial signal for information treatment
function solve_trip(Ss::Drip,        # D.R.I.P.
                    S::Signal;       # initial Signal for information treatment
                    T     = 100,     # optional: time until convergence to steady state
                    tol   = 1e-4,    # optional: tolerance for convergence
                    maxit = 1000     # optional: max iterations
                    )
    K0   = Ss.Σ_1*S.L/(S.L'*Ss.Σ_1*S.L+S.Σ_z); # Get prior covariance matrix from new signal
    Σ0   = Ss.Σ_1-K0*S.L'*Ss.Σ_1;              # Get Kalman gain from new signal
    trip = solve_trip(Ss,Σ0;T=T,tol=tol,maxit=maxit);
    return(trip)
end

########## Aux. Functions ############

function getreal(M)
    #if maximum(abs.(imag.(M))) > 1e-6
    #    print("Warning: Matrix decomposition returned complex numbers larger than 1e-6")
    #end
    return(real.(M))
end

function infinitesum(func; tol = 1e-6,maxit = 1000,start=0)
    diff  = 1.0
    infsum = func(start)
    it    = start + 1
    while (diff > tol) & (it < maxit)
        func_it = func(it)
        infsum += func_it
        diff = maximum(func_it)
        it += 1
    end
    return(infsum)
end

## calculate the amount of information acquired in bits/nats for a Drip
function capacity(P::Drip;      # Drip structure
                  unit = "bit"  # optional: unit of capacity (bit or nat).
                  )
    if unit == "bit"
        κ = 0.5*log(det(P.Σ_1)/det(P.Σ_p))/log(2); #returns capacity in bits
    elseif unit == "nat"
        κ = 0.5*log(det(P.Σ_1)/det(P.Σ_p));        #returns capacity in nats
    else
        println("Invalid input for unit! Capacity is reported in bits.")
        κ = 0.5*log(det(P.Σ_1)/det(P.Σ_p))/log(2);
    end
    return(κ)
end

## get steady state irfs
function dripirfs(P::Drip,T::Int)
    (n,m) = length(size(P.H)) == 2 ? size(P.H) : (size(P.H,1),1)
    (_,k) = length(size(P.Q)) == 2 ? size(P.Q) : (size(P.Q,1),1)
    x     = zeros(n,k,T);
    x_hat = zeros(n,k,T);
    a     = zeros(m,k,T);
    for kk in 1:k
        e_k = zeros(k,1); e_k[kk] = 1;
        for ii in 1:T
            if ii==1
                x[:,kk,ii]     = P.Q*e_k;
                x_hat[:,kk,ii] = (P.K*P.Y')*(x[:,kk,ii]);
            else
                x[:,kk,ii]     =P.A*x[:,kk,ii-1];
                x_hat[:,kk,ii] =P.A*x_hat[:,kk,ii-1]+(P.K*P.Y')*(x[:,kk,ii]-P.A*x_hat[:,kk,ii-1]);
            end
            a[:,kk,ii]  .= P.H'*x_hat[:,kk,ii];
        end
    end
    return(Dripirfs(T,x,x_hat,a))
end

## get irfs with optimal signals acquired under Pt::Trip
function dripirfs(Pt::Trip,T::Int)
    Ss    = Pt.P; # steady state Drip
    (n,m) = length(size(Ss.H)) == 2 ? size(Ss.H) : (size(Ss.H,1),1) # dimensions of state and actions
    (_,k) = length(size(Ss.Q)) == 2 ? size(Ss.Q) : (size(Ss.Q,1),1) # number of structural shocks
    eye   = Matrix(I,n,n);

    x     = zeros(n,k,T); # initialize IRFs of state
    x_hat = zeros(n,k,T); # initialize IRFs of beliefs
    a     = zeros(m,k,T); # initialize IRFs of actions
    for kk in 1:k
        e_k = zeros(k,1); e_k[kk] = 1;
        for ii in 1:T
            if ii==1
                x[:,kk,ii]     = Ss.Q*e_k;
                x_hat[:,kk,ii] = (eye-Pt.Σ_ps[:,:,ii]*pinv(Pt.Σ_1s[:,:,ii]))*(x[:,kk,ii]);
            else
                x[:,kk,ii]     =Ss.A*x[:,kk,ii-1];
                x_hat[:,kk,ii] =Ss.A*x_hat[:,kk,ii-1]+(eye-Pt.Σ_ps[:,:,ii]*pinv(Pt.Σ_1s[:,:,ii]))*(x[:,kk,ii]-Ss.A*x_hat[:,kk,ii-1]);
            end
            a[:,kk,ii]  .= Ss.H'*x_hat[:,kk,ii];
        end
    end
    return(Dripirfs(T,x,x_hat,a))
end

## get irfs with information treatment S at time zero.
function dripirfs(Ss::Drip,          # Steady state Drip (when treatment happens)
                  T::Int,            # Length of irfs
                  S::Signal;         # Signal for treatment
                  reoptimize = true, # if true gives the irfs with reoptimized signals, if false with steady state signals
                  trip       = false # if false solves the trip, if = P::trip then takes P as the trip
                  )

    (n,m) = length(size(Ss.H)) == 2 ? size(Ss.H) : (size(Ss.H,1),1) # dimensions of state and actions
    (_,k) = length(size(Ss.Q)) == 2 ? size(Ss.Q) : (size(Ss.Q,1),1) # number of structural shocks
    eye   = Matrix(I,n,n);

    if (trip == false) & (reoptimize == true) # get the solution to trip if reoptimize = true
        trip = solve_trip(Ss,S;T = T);
    end
    K0   = Ss.Σ_1*S.L/(S.L'*Ss.Σ_1*S.L+S.Σ_z); # Get prior covariance matrix from new signal
    Σ0   = trip.Σ_1s[:,:,1];

    x     = zeros(n,k,T); # initialize IRFs of state
    x_hat = zeros(n,k,T); # initialize IRFs of beliefs
    a     = zeros(m,k,T); # initialize IRFs of actions
    for kk in 1:k
        e_k = zeros(k,1); e_k[kk] = 1;
        for ii in 1:T
            if ii==1
                x[:,kk,ii]     = Ss.Q*e_k;
                x_hat_p        = K0*S.L'*(x[:,kk,ii]);
                if reoptimize == true
                    x_hat[:,kk,ii] = x_hat_p + (eye-trip.Σ_ps[:,:,ii]*pinv(trip.Σ_1s[:,:,ii]))*(x[:,kk,ii]-x_hat_p);
                else
                    K0 = Σ0*Ss.Y/(Ss.Y'*Σ0*Ss.Y+Ss.Σ_z);
                    Σ0 = Ss.Q*Ss.Q' + Ss.A*(Σ0-K0*Ss.Y'*Σ0)*Ss.A';
                    x_hat[:,kk,ii] = x_hat_p + (K0*Ss.Y')*(x[:,kk,ii]-x_hat_p);
                end
            else
                x[:,kk,ii]     =Ss.A*x[:,kk,ii-1];
                if reoptimize  == true
                    x_hat[:,kk,ii] =Ss.A*x_hat[:,kk,ii-1]+(eye-trip.Σ_ps[:,:,ii]*pinv(trip.Σ_1s[:,:,ii]))*(x[:,kk,ii]-Ss.A*x_hat[:,kk,ii-1]);
                else
                    K0 = Σ0*Ss.Y/(Ss.Y'*Σ0*Ss.Y+Ss.Σ_z);
                    Σ0 = Ss.Q*Ss.Q' + Ss.A*(Σ0-K0*Ss.Y'*Σ0)*Ss.A';
                    x_hat[:,kk,ii] =Ss.A*x_hat[:,kk,ii-1]+(K0*Ss.Y')*(x[:,kk,ii]-Ss.A*x_hat[:,kk,ii-1]);
                end
            end
            a[:,kk,ii]  .= Ss.H'*x_hat[:,kk,ii];
        end
    end
    return(Dripirfs(T,x,x_hat,a))
end


text:
# Julia program to find the number of bits to be flipped to convert one number to another

"""Inorder to compare the individual bits of given the numbers, we can compute the number of
set bits in thier XOR, as we know that XOR operation yields set bits, only when two same
bits are XORed."""
function get_num_flips(n1, n2)
    count = 0
    n = n1 ⊻ n2
    while(n != 0)
        if(n & 1 == 1)
            count = count + 1
        end
        n = n >> 1
    end
    return count
end

print("Enter the first number:")
n1 = readline()
n1 = parse(Int32, n1)
print("Enter the second number:")
n2 = readline()
n2 = parse(Int32, n2)
res = get_num_flips(n1, n2)
print("Number of flips required to inter-convert the given numbers is: $res")

"""
Time Complexity: O(log(n)), where n is the greater of the input number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

Enter the first number:13
Enter the second number:26
Number of flips required to inter-convert the given numbers is: 4
"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
include("calc_n_bar.jl")

include("calc_r_0.jl")
include("calc_b_0.jl")

include("calc_i_p.jl")

include("calc_sigma_v.jl")

include("calc_eta_cd.jl")


text:
using Distributions, Plots; pyplot()

n1, n2 = 10, 15
N = 10^6
mu, sigma = 10, 4
normDist = Normal(mu,sigma)

fValues = Array{Float64}(undef, N)

for i in 1:N
    data1 = rand(normDist,n1)
    data2 = rand(normDist,n2)
    fValues[i] = var(data1)/var(data2)
end

fRange = 0:0.1:5
stephist(fValues, bins=400, c=:blue, label="Simulated", normed=true)
plot!(fRange, pdf.(FDist(n1-1, n2-1), fRange),
	c=:red, label="Analytic", xlims=(0,5), ylims=(0,0.8),
    	xlabel = "F", ylabel = "Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
abstract type LogicalNode
end

Base.show(io::IO, node::LogicalNode) = Base.show(materialize(node))


text:
"""
Check that the dynamic program is solved
"""
function checkdp(model::RustModelCore)
	V1=\(eye(model.pi[1])-model.beta*model.pi[1],model.u[:,1]-log(model.p[:,1]))
	V2=\(eye(model.pi[2])-model.beta*model.pi[2],model.u[:,2]-log(model.p[:,2]))
	norm(V1-V2)
end

#changing autoscale option might solve numerical issues sometimes
function solvedp!(model::DIRMC)
	da=size(model.u)[2]

	function ff!(z,fv)
		fv[:]=-1.0
		for ia=1:da
			fv[:]+=exp(model.u[:,ia]+model.beta*model.pi[ia]*z-z)
		end
	end

	function fj!(z,fm)
		fm[:]=0.0							#inefficient reallocation of temp at each call
		for ia=1:da
			fm[:,:]+=diagm(exp(model.u[:,ia]+model.beta*model.pi[ia]*z-z))*(model.beta*model.pi[ia]-eye(size(fm,1)))
		end
	end

	df = NLsolve.DifferentiableMultivariateFunction(ff!,fj!)

	res=NLsolve.nlsolve(df,model.V,autoscale=false)	#will automatically start at previous value

	# println(res)
	model.V[:]=res.zero
	for ia=1:da
		model.p[:,ia]=exp(model.u[:,ia]+model.beta*model.pi[ia]*model.V-model.V)
	end
	# println(round(model.p,2))
end

function solvedp!(model::SIRMC)
	da=size(model.u)[2]

	function ff!(z,fv)
		fv[:]=-1.0
		for ia=1:da
			# println(sum(model.u[:,ia]+model.beta*model.pi[ia]*z-z))
			fv[:]+=exp(model.u[:,ia]+model.beta*model.pi[ia]*z-z)
		end
	end

	function fj!(z,fm)
		dpj=spzeros(size(fm)...)							#inefficient reallocation of temp at each call
		for ia=1:da
			dpj+=spdiagm(exp(model.u[:,ia]+model.beta*model.pi[ia]*z-z))*(model.beta*model.pi[ia]-speye(size(fm)[1]))
		end
		fm.nzval=dpj.nzval
		fm.rowval=dpj.rowval
		fm.colptr=dpj.colptr
	end

	df = NLsolve.DifferentiableSparseMultivariateFunction(ff!, fj!)
	res= NLsolve.nlsolve(df,model.V,autoscale=false)	#will automatically start at previous value
	# println(res)
	model.V[:]=res.zero
	for ia=1:da
		model.p[:,ia]=exp(model.u[:,ia]+model.beta*model.pi[ia]*model.V-model.V)
	end
	# println(round(model.p,2))
end

### FINITE HORIZON METHODS

#Upon exiting any of the finite-horizon methods, model.v or model.V are NOT calibrated to any meaningful value.

#this works as long as observable states are ranked by chronological order.
function solvedp!(model::DFRMC)
	dk,da=size(model.p)

	model.V[:]=0		#could be slightly more efficient here, see the sparsefinitehorizon case
	for ik in dk:-1:1			#this works as long as observable states are ranked by chronological order.
		for ia=1:da
			model.v[ik,ia]=model.u[ik,ia]+model.beta*dot(full(model.pi[ia])[ik,:],model.V)
			model.p[ik,ia]=exp(model.v[ik,ia]-model.v[ik,1])				#for numerical stability
		end
		model.p[ik,:]/=sum(model.p[ik,:])
		model.V[ik]=dot(model.p[ik,:],model.v[ik,:]-log(model.p[ik,:]))
	end
end

function solvedp!(model::SFRMC)
	dk,da=size(model.p)

	#I could do one loop instead of two
	#but speed matters
	# local pisix::Float64=0
	for ik in model.horizonindex:dk
		model.V[ik]=0
		rho=0
		for ia=1:da
			model.p[ik,ia]=exp(model.u[ik,ia]-model.u[ik,1])				#for numerical stability
			rho+=model.p[ik,ia]
		end
		for ia=1:da
			model.p[ik,ia]/=rho
			model.V[ik]+=model.p[ik,ia]*(model.u[ik,ia]-log(model.p[ik,ia]))
		end
	end
	for ik in model.horizonindex-1:-1:1			#this works as long as observable states are ranked by chronological order.
		model.V[ik]=0
		rho=0
		for ia=1:da
			model.v[ik,ia]=0
			# println(sum(full(model.pi[ia])[ik,model.kindices[ik,ia]]))
			for jk in model.kindices[ik,ia]
				model.v[ik,ia]+=model.pi[ia][ik,jk]*model.V[jk]
			end
			model.v[ik,ia]*=model.beta
			model.v[ik,ia]+=model.u[ik,ia]
			model.p[ik,ia]=exp(model.v[ik,ia]-model.v[ik,1])				#for numerical stability
			rho+=model.p[ik,ia]
		end
		for ia=1:da
			model.p[ik,ia]/=rho
			model.V[ik]+=model.p[ik,ia]*(model.v[ik,ia]-log(model.p[ik,ia]))
		end
	end
	# println(sort(vec(model.p)))
end

## DEPRECATED
## keep around for now for testing purposes

function solvedp_old!(model::RustModelCore)
	da=size(model.u)[2]

	function ff!(z,fv)
		fv[:]=-1
		for ia=1:da
			fv[:]=fv+ exp(model.u[:,ia]+model.beta*model.pi[ia]*z-z)
		end
	end

	function fj!(z,fm)
		dpj=spzeros(size(fm)...)							#inefficient reallocation of temp at each call
		for ia=1:da
			dpj+=spdiagm(exp(model.u[:,ia]+model.beta*model.pi[ia]*z-z))*(model.beta*model.pi[ia]-speye(size(fm)[1]))
		end
		fm.nzval=dpj.nzval
		fm.rowval=dpj.rowval
		fm.colptr=dpj.colptr
	end

	df = NLsolve.DifferentiableSparseMultivariateFunction(ff!, fj!)
	res= NLsolve.nlsolve(df,model.V,autoscale=false)	#will automatically start at previous value
	# println(res)
	model.V[:]=res.zero
	for ia=1:da
		model.p[:,ia]=exp(model.u[:,ia]+model.beta*model.pi[ia]*model.V-model.V)
	end
end


text:
"""
    show(sresult::Simulation)

Print out summary of simulation.
"""
function Base.show(io::IO, sresult::Simulation)
    @printf("===================================================================\n")
    _showtype(sresult)
    @printf("Input parameters: \n")
    @printf("\tMutation rate: %.2f\n", sresult.input.μ)
    _showrates(sresult.input)
    @printf("\tNumber of clonal mutations: %d\n", sresult.input.clonalmutations)
    @printf("\tNumber of subclones: %d\n\n", sresult.input.numclones)
    if sresult.input.numclones > 0
        freqs = subclonefreq(sresult.output)[1]
        for (i, subclone) in enumerate(sresult.output.subclones)
            @printf("Subclone %d \n", i)
            @printf("\tFrequency: %.2f\n", freqs[i])
            @printf("\tNumber of mutations in subclone: %d\n", length(subclone.mutations))
            @printf("\tFitness advantage: %.2f\n", sresult.input.selection[i])
            @printf("\tTime subclone emerges: %.2f\n", subclone.time)
            @printf("\tNumber of divisions: %d\n", subclone.Ndivisions)
            @printf("\tAverage number of divisions per cell: %.2f\n", subclone.avdivisions)
            @printf("\tPopulation size when subclone emerges: %d\n", subclone.N0)
            @printf("\tParent of subclone (0 is host): %d\n\n", subclone.parenttype-1)
        end
    else
        @printf("No clones, tumour growth was neutral\n\n")
    end

end

function Base.show(io::IO, sresult::MultiSimulation)
    @printf("===================================================================\n")
    @printf("Number of simulations = %d\n\n", length(sresult.output))
    @printf("Input parameters: \n")
    _showrates(sresult.input)
    @printf("\tNumber of clonal mutations: %d\n", sresult.input.clonalmutations)
    @printf("\tNumber of subclones: %d\n\n", sresult.input.numclones)

    if sresult.input.numclones > 0
        for i in 1:sresult.input.numclones
            @printf("Subclone %d \n", i)
            @printf("\tFitness advantage: %.2f\n", sresult.input.selection[i])
            @printf("\tEvent time: %.2f\n",
                sresult.input.tevent[i])
            @printf("\tAverage frequency: %.2f (std: %.2f)\n",
                        _get_mean_std(sresult.output, :freq, i)...)
            @printf("\tAverage number of mutations in subclone: %d (std: %.2f)\n",
                        _get_mean_std(sresult.output, :mutations, i)...)
            @printf("\tAverage population size when subclone emerges: %d (std: %.2f)\n",
                        _get_mean_std(sresult.output, :N0, i)...)
            @printf("\tAverage time when subclone emerges: %.2f (std: %.2f)\n",
                        _get_mean_std(sresult.output, :time, i)...)
            @printf("\tDistribution of subclone parents (0 is host): ")
            s = join([@sprintf("%d => %d", parent, freq)
                        for (parent, freq) in parenttype_dist(sresult.output, i)]
                        , ", ")
            @printf("%s\n\n",s)
        end
    else
        @printf("No clones, tumour growth was neutral\n\n")
    end

end

function Base.show(io::IO, population::Population)
    siminput = population.input
    @printf("===================================================================\n")
    @printf("Multilevel branching module simulation\n\n")
    @printf("Input parameters: \n")
    @printf("\tMax age = %.2f\n", siminput.maxtime)
    _showrates(siminput)
    @printf("\tBranching rate = %.4f\n", siminput.branchrate)
    @printf("\tNew module size = %d\n", siminput.branchinitsize)
    @printf("\tMature module size = %d\n", siminput.modulesize)
    @printf("\tNumber of clonal mutations: %d\n", siminput.clonalmutations)
    # @printf("\tNumber of subclones: %d\n\n", siminput.numclones)

    @printf("Population data\n")
    @printf("\tPopulation age = %d\n", age(population))
    @printf("\tNumber of modules = %d\n\n", length(population))

    # @printf("Clonal mutations per module")

    # if siminput.numclones > 0
    #     for i in 1:siminput.numclones
    #         @printf("Subclone %d \n", i)
    #         @printf("\tFitness advantage: %.2f\n", siminput.selection[i])
    #         @printf("\tEvent time: %.2f\n",
    #             siminput.tevent[i])
    #         @printf("\tAverage frequency: %.2f (std: %.2f)\n",
    #                     _get_mean_std(population.output, :freq, i)...)
    #         @printf("\tAverage number of mutations in subclone: %d (std: %.2f)\n",
    #                     _get_mean_std(population.output, :mutations, i)...)
    #         @printf("\tAverage population size when subclone emerges: %d (std: %.2f)\n",
    #                     _get_mean_std(population.output, :N0, i)...)
    #         @printf("\tAverage time when subclone emerges: %.2f (std: %.2f)\n",
    #                     _get_mean_std(population.output, :time, i)...)
    #         @printf("\tDistribution of subclone parents (0 is host): ")
    #         s = join([@sprintf("%d => %d", parent, freq)
    #                     for (parent, freq) in parenttype_dist(population.output, i)]
    #                     , ", ")
    #         @printf("%s\n\n",s)
    #     end
    # else
    #     @printf("No clones, tumour growth was neutral\n\n")
    # end

end

_showtype(sresult::Simulation{BranchingInput}) = @printf("Branching process \n")
_showtype(sresult::Simulation{MoranInput}) = @printf("Moran process \n")
_showtype(sresult::Simulation{BranchingMoranInput}) = @printf("Branching => Moran process\n")
_showtype(sresult::Simulation{MultilevelInput}) = @printf("Branching => Moran proces (multilevel)\n")

function _showrates(siminput::BranchingInput)
    b, d, μ = siminput.b, siminput.d, siminput.μ
    @printf("\tBirth rate of host population: %.2f\n", b)
    @printf("\tDeath rate of host population: %.2f\n", d)
    @printf("\tEffective mutation rate (μ/β): %.2f\n", μ / ((b - d) / b))
end

function _showrates(siminput::MoranInput)
    @printf("\tBirth and death rate of host population: %.2f\n", siminput.bdrate)
end

function _showrates(siminput::Union{BranchingMoranInput, MultilevelInput})
    b, d, μ, bdrate = siminput.b, siminput.d, siminput.μ, siminput.bdrate
    @printf("\tBirth rate of host population (branching): %.2f\n", b)
    @printf("\tDeath rate of host population (branching): %.2f\n", d)
    @printf("\tBirth and death rate of host population (moran): %.2f\n", bdrate)
    @printf("\tMutation rate (μ): %.2f\n", μ)
end

function _get_mean_std(output, subclonefield, i)
    data = [getfield(subclone, subclonefield)[i]
                for simoutput in output
                    for subclone in simoutput.subclones]
    return mean(data), std(data)
end

function parenttype_dist(output, i)
    parenttypes = [subclone.parenttype[i] - 1
                    for simoutput in output
                        for subclone in simoutput.subclones]
    parenttypes = sort(collect(countmap(parenttypes)))
    return parenttypes
end


text:
using Optim: optimize
using LinearAlgebra: diag, diagm, eigvals, ishermitian

export QSpec, TransmonSpec, DuffingSpec, ResonatorSpec, HermitianSpec
export QSystem, label, dimension, spec
export LiteralHermitian, Resonator, DuffingTransmon
export PerturbativeTransmon, ChargeBasisTransmon, DiagonalChargeBasisTransmon, RotatingFrameSystem
export hamiltonian, fit_transmon

######################################################
# Types
######################################################

HermitianMatrix = Matrix{<:Number}
abstract type QSpec end
abstract type QSystem end

# Required QSystem functions with defaults
label(q::QSystem)::AbstractString = q.label
dimension(q::QSystem)::Int64 = q.dimension
spec(q::QSystem)::QSpec = q.spec

################################################################################
# QSpecs
# QSpecs are the physical parameters that define a quantum system
################################################################################

# Structs

struct HermitianSpec <: QSpec
   matrix::HermitianMatrix
   function HermitianSpec(matrix::Matrix{<:Number})
       @assert ishermitian(matrix)
       return new(matrix)
   end
end

"""
Resonator, specified by its level spacing.
"""
struct ResonatorSpec <: QSpec
   frequency::Real
end

"""
Transmon, specified by its charge and junctions energies, EC and EJ1/EJ2, respectively.
"""
struct TransmonSpec <: QSpec
   EC::Real
   EJ1::Real
   EJ2::Real
end

"""
Alternative Transmon specification by its f₀₁ transition frequency and anharmonicity (f₁₂ - f₀₁).
"""
struct DuffingSpec <: QSpec
   frequency::Real
   anharmonicity::Real
end

# Functions

"""
Construct a `TransmonSpec` (EC, EJ1, EJ2) from a `DuffingSpec` (frequency, anharmonicity), and a model
type to perform the conversion.

An optimization routine, fit_transmon, must be run to convert frequency/anharmonicity to EC/EJ.
"""
function TransmonSpec(t::DuffingSpec, model::Type{T}=PerturbativeTransmon,
    num_terms::Int=PERTURBATIVE_NUM_TERMS) where {T<:QSystem}
    return fit_transmon(t.frequency, t.frequency, t.anharmonicity, model, num_terms)
end

"""
Construct a `DuffingSpec` (frequency, anharmonicity) from a TransmonSpec (EC, EJ1, EJ2), at a given
flux quantum ϕ. May optionally supply a custom number of perturbative terms.
"""
function DuffingSpec(t::TransmonSpec, ϕ::Real=0.0, num_terms::Int=PERTURBATIVE_NUM_TERMS)
    freq = perturbative_transmon_freq(t.EC, t.EJ1, t.EJ2, ϕ, num_terms=num_terms)
    anharm = perturbative_transmon_anharm(t.EC, t.EJ1, t.EJ2, ϕ, num_terms=num_terms)
    return DuffingSpec(freq, anharm)
end

################################################################################
# QSystems
# QSystems are Structs comprised of a QSpec, a label and a dimension.
################################################################################

# Structs

struct LiteralHermitian <: QSystem
    label::AbstractString
    spec::HermitianSpec
end

dimension(h::LiteralHermitian)::Int64 = size(h.spec.matrix, 1)
hamiltonian(q::LiteralHermitian)::HermitianMatrix = spec(q).matrix

struct Resonator <: QSystem
    label::AbstractString
    dimension::Int
    spec::ResonatorSpec
end

struct DuffingTransmon <: QSystem
    label::AbstractString
    dimension::Int
    spec::DuffingSpec
end

"""
The `PerturbativeTransmon` type is a `DuffingTransmon` with the f₀₁ and f₁₂ calculated from a
perturbative expansion. For details see Didier, N., Sete, E. A., da Silva, M. P., & Rigetti, C.
(2017). Analytical modeling of parametrically-modulated transmon qubits.
http://arxiv.org/abs/1706.06566 for details.
"""
struct PerturbativeTransmon <: QSystem
    label::AbstractString
    dimension::Int
    spec::TransmonSpec
    num_terms::Int

    PerturbativeTransmon(
        label::AbstractString, dimension::Int, spec::TransmonSpec;
        num_terms::Int=PERTURBATIVE_NUM_TERMS
    ) =
        new(label::AbstractString, dimension::Int, spec::TransmonSpec, num_terms)
end

struct ChargeBasisTransmon <: QSystem
    label::AbstractString
    dimension::Int
    spec::TransmonSpec
    function ChargeBasisTransmon(label::AbstractString, dimension::Int, spec::TransmonSpec)
        @assert mod(dimension, 2) == 1
        return new(label, dimension, spec)
    end
end

const CHARGE_NUM_TERMS = 101

struct DiagonalChargeBasisTransmon <: QSystem
    label::AbstractString
    dimension::Int
    spec::TransmonSpec
    num_terms::Int
    function DiagonalChargeBasisTransmon(label::AbstractString, dimension::Int, spec::TransmonSpec; num_terms::Int=CHARGE_NUM_TERMS)
        @assert mod(num_terms, 2) == 1
        @assert dimension <= num_terms
        return new(label, dimension, spec, num_terms)
    end
end

struct RotatingFrameSystem <: QSystem
    label::AbstractString
    system::QSystem
    rotation_rate::Real
end

# contruct from an existing QSystem without a different label
RotatingFrameSystem(q::QSystem, rotation_rate::Real) = RotatingFrameSystem(label(q), q, rotation_rate)

# defer `dimension` and `spec` to underlying QSystem
dimension(r::RotatingFrameSystem) = dimension(r.system)
spec(r::RotatingFrameSystem) = spec(r.system)

################################################################################
# Hamiltonians
################################################################################

# Resonator

"""
Construct a resonator hamiltonian from energy levels equally spaced by the resonator frequency,
along the diagonal of a square matrix.
"""
hamiltonian(r::Resonator)::HermitianMatrix =
    diagm(0 =>[spec(r).frequency * n for n in 0:dimension(r)-1])

# Duffing Transmon

"""
Construct the hamiltonian of a duffing transmon - constant relative anharmonicity.
```
E_{n+1} - E_{n} = f₀₁ - n*anharmonicity.
```
"""
function duffing_hamiltonian(
    freq::Real, anharm::Real, dimension::Int
)::HermitianMatrix
    diagm(0 => [(freq - 0.5 * anharm) * n + 0.5 * anharm * n^2
                for n in 0:dimension-1])
end

"""
Construct a hamiltonian for a Duffing transmon.
"""
hamiltonian(t::DuffingTransmon)::HermitianMatrix =
    duffing_hamiltonian(spec(t).frequency, spec(t).anharmonicity, dimension(t))

"""
Construct a hamiltonian for a Duffing transmon given a specific f₀₁
"""
hamiltonian(t::DuffingTransmon, frequency::Real)::HermitianMatrix =
    duffing_hamiltonian(frequency, spec(t).anharmonicity, dimension(t))

# Perturbative Transmon

"""
Construct a hamiltonian from a PerturbativeTransmon at a given external flux ϕ in units of Φ₀
"""
function hamiltonian(t::PerturbativeTransmon, ϕ::Real=0.0)::HermitianMatrix
    s = DuffingSpec(spec(t), ϕ, t.num_terms)
    return duffing_hamiltonian(s.frequency, s.anharmonicity, dimension(t))
end

# Charge Basis Transmon

function hamiltonian(t::ChargeBasisTransmon, ϕ::Real=0.0)
    d = dimension(t)
    N = floor(Int, d/2)
    s = spec(t)
    EJ = sqrt(s.EJ1^2 + s.EJ2^2 + 2 * s.EJ1 * s.EJ2 * cos(2π * ϕ))
    EC = s.EC
    charging_term = 4 * EC * diagm(0 => (-N:N).^2)
    tunneling_term = -0.5 * EJ * (diagm(-1 => ones(d-1), 1 => ones(d-1)))
    return charging_term + tunneling_term
end

function hamiltonian(t::DiagonalChargeBasisTransmon, ϕ::Real=0.0)
    h = hamiltonian(ChargeBasisTransmon(label(t), t.num_terms, spec(t)), ϕ)
    # since the Hamiltonian is Hermitian the eigenvalues should already be sorted by the LAPACK
    # solver; however,  since that implementation is not guaranteed by Julia, belts and suspenders
    # style we sort again
    return diagm(0 => sort(real(eigvals(h)))[1:dimension(t)])
end
# TODO: update raising and lowering to be correct for DiagonalChargeBasisTransmon

# Rotating Frame

function hamiltonian(r::RotatingFrameSystem)
    h = hamiltonian(r.system)
    subtract_number!(h, r.rotation_rate)
    return h
end

function hamiltonian(r::RotatingFrameSystem, t::Real)
    h = hamiltonian(r.system, t)
    subtract_number!(h, r.rotation_rate)
    return h
end

# Util

"""
Helper function that subtracts a real number off the diagonal elements of a Real matrix.
"""
function subtract_number!(h::Matrix{<:Real}, r::Real)
    if r != 0
        for i in 1:size(h, 1)
            h[i,i] -= (i-1) * r
        end
    end
end

################################################################################
# Conversion fitting (EC/EJ <> Frequency/Anharmonicity)
################################################################################

"""
    fit_transmon(freq_max::Real, freq_min::Real, anharm_max::Real, model::Type{T}, num_terms::Int) where {T<:QSystem}

Fit a transmon to given frequency and anharmonicity at fmax and the frequency at fmin. The model can
be either PerturbativeTransmon or DiagonalChargeBasisTransmon. If `freq_max == freq_min` then `EJ2
== 0` is enforced.

## args
* `freq_max`: the qubit maximum frequency.
* `freq_min`: the qubit minimum frequency.
* `anharm_max`: the qubit maximum anharmonicity.
* `model`: either PerturbativeTransmon or DiagonalChargeBasisTransmon.
* `num_terms`: the number of terms used in the model.

## returns
A TransmonSpec.
"""
function fit_transmon(
    freq_max::Real, freq_min::Real, anharm_max::Real, model::Type{T}, num_terms::Int
)::TransmonSpec where {T<:Union{PerturbativeTransmon,DiagonalChargeBasisTransmon}}
    function f_fixed(params)
      EC, EJ = params
      t = model("", 3, TransmonSpec(EC, EJ, 0.0), num_terms=num_terms)
      levels = diag(hamiltonian(t))
      test_f_01 = levels[2] - levels[1]
      test_f_12 = levels[3] - levels[2]
      return abs(test_f_01 - freq_max) + abs(test_f_12 - test_f_01 - anharm_max)
    end

    function f_tunable(params)
        t = model("", 3, TransmonSpec(params...), num_terms=num_terms)
        levels_fmax = diag(hamiltonian(t, 0.0))
        levels_fmin = diag(hamiltonian(t, 0.5))
        test_f_01_max = levels_fmax[2] - levels_fmax[1]
        test_f_12_max = levels_fmax[3] - levels_fmax[2]
        test_f_01_min = levels_fmin[2] - levels_fmin[1]
        return abs(test_f_01_max - freq_max) + abs(test_f_12_max - test_f_01_max - anharm_max) + abs(test_f_01_min - freq_min)
    end

    EC_guess = -anharm_max
    EJ_max_guess = -(freq_max - anharm_max)^2 / (8 * anharm_max)
    EJ_min_guess = -(freq_min - anharm_max)^2 / (8 * anharm_max)
    EJ1_guess = 1//2 * (EJ_max_guess + EJ_min_guess)
    EJ2_guess = 1//2 * abs(EJ_max_guess - EJ_min_guess)
    if freq_max == freq_min
        res = optimize(f_fixed, [EC_guess, EJ1_guess])
        EC, EJ = res.minimizer
        return TransmonSpec(EC, EJ, 0.0)
    else
        res = optimize(f_tunable, [EC_guess, EJ1_guess, EJ2_guess])
        return TransmonSpec(res.minimizer...)
    end
end


text:
#Sampling methods
"""
thompson_sampling(f, x::AbstractVector,n_samples, σ²)
Brute force thompson sampling over the inter undiscovered search space.
"""

function thompson_sampling(f, x::AbstractVector,
                n_samples, σ²)
    @assert length(x) ≥ n_samples

    # indices selected
    selected = Int[]
    all_test = copy(x)
    Ndist = f(all_test, σ²)
    while length(selected) < n_samples
        f̂ = rand(Ndist)
        f̂[selected] .= - Inf
       i = argmax(f̂)
       push!(selected, i)
   end
   return (x[selected],selected)
end

"""
    save_thompson_sampling(f, x::AbstractVector,n_samples, σ²)
A brute force thompson sampling algorithm that uses the initer undiscovered search space.
Save version that uses the eigenvalue decomposition of the posterior covariance matrix to assure that the covariance matrix is positive definite
"""

function save_thompson_sampling(f_post, constructs_test::AbstractVector,n_samples,σ²_opt)
    myeigen = (-(eigvals!(cov(f_post(constructs_test)))[1])) + σ²_opt
    if  myeigen > 0
        println("usetrick")
        sigma_add = myeigen+σ²_opt
        nextsample = thompson_sampling(f_post,constructs_test,n_samples,sigma_add)
    else
        nextsample = thompson_sampling(f_post,constructs_test,n_samples,σ²_opt)
    end
    return nextsample
end

#Sampling methods
"""
thompson_sampling(f, x::AbstractVector,n_samples, σ²)
Brute force thompson sampling over the inter undiscovered search space.
"""

function thompson_sampling_fast(f, x::AbstractVector,
                n_samples, σ²)
    @assert length(x) ≥ n_samples

    # indices selected
    selected = Int[]
    all_test = copy(x)
    Ndist = f(all_test, σ²)
    f̂_matrix = rand(Ndist,n_samples)
    for n in eachcol(f̂_matrix)
        n[selected]  .= - Inf
       i = argmax(n)
       push!(selected, i)
   end
   return (x[selected],selected)
end


text:
function main()
    setupgui()
    return
end

function setupgui()
    heading = "ACRLab 2: Inverted Pendulum"

    sub1 = "Pendulum Swing-Up. Set \"F\", \"xinit\", and \"theta1init\"."
    txt1 = textarea("F(theta,dtheta,Vis) = ...")
    btn1 = button("Run Simulation")

    sub2 = "Implementation of a PID Controller. Set \"F\", \"xinit\", and \"theta1init\"."
    txt2 = textarea("function F(e,inte,de,theta,dtheta,Vis)\n\t...\nend")
    btn2 = button("Run Simulation")

    sub3 = "Implementation of an LQ Regulator. Set \"F\", \"xinit\", \"theta1init\", \"Q\", and \"R\"."
    txt3 = textarea("F(z,K,t) = ...")
    btn3 = button("Run Simulation")

    sub4 = "Double Pendulum. Set \"F\", \"xinit\", \"theta1init\", \"theta2init\", \"Q\", and \"R\"."
    txt4 = textarea("F(z,K,t) = ...")
    btn4 = button("Run Simulation")

    widgets = [txt1;btn1;txt2;btn2;txt3;btn3;txt4;btn4]
    widgetnames = ["txt1";"btn1";"txt2";"btn2";"txt3";"btn3";"txt4";"btn4"]
    dict = Dict{String,Widget}()
    for (i,el) in enumerate(widgets)
        dict[widgetnames[i]] = el
    end

    ui = map(update!,widgets...,dict)

    layout = vbox(
        ui,
        heading,
        vskip(1em),
        sub1,
        vskip(1em),
        txt1,
        vskip(1em),
        btn1,
        vskip(1em),
        sub2,
        vskip(1em),
        txt2,
        vskip(1em),
        btn2,
        vskip(1em),
        sub3,
        vskip(1em),
        txt3,
        vskip(1em),
        btn3,
        vskip(1em),
        sub4,
        vskip(1em),
        txt4,
        vskip(1em),
        btn4
    )

    mkpath("Files")

    window = Window()
    body!(window, layout)
end


text:
n = rand(1:10)

A = rand(n,n)

r = polarfact(A, alg = :newton);

# Test unitary matrix U

U = r.U
H = r.H

@test_approx_eq_eps U'*U eye(n) 1e-7

# Test Hermitian positive semifefinite matrix H

@test issym(H)

for i in eigvals(H)
    @test i >= 0.
end

@test_approx_eq A U*H

println("Newton method passed test...")


text:
# Reference: https://neuralpde.sciml.ai/dev/examples/pinns_example/
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
# 3D PDE
@parameters x y t θ
@variables u1(..)
@derivatives Dxx''~x
@derivatives Dyy''~y
@derivatives Dt'~t

# 3D PDE
eq  = Dt(u1(x,y,t,θ)) ~ Dxx(u1(x,y,t,θ)) + Dyy(u1(x,y,t,θ))
# Initial and boundary conditions
bcs = [u1(x,y,0,θ) ~ exp(x+y)*cos(x+y) ,
       u1(0,y,t,θ) ~ exp(y)*cos(y+4t),
       u1(2,y,t,θ) ~ exp(2+y)*cos(2+y+4t) ,
       u1(x,0,t,θ) ~ exp(x)*cos(x+4t),
       u1(x,2,t,θ) ~ exp(x+2)*cos(x+2+4t)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           y ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)]

# Discretization
dx = 0.25; dy= 0.25; dt = 0.25
# Neural network
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN([dx,dy,dt],
                                             chain,
                                             strategy = NeuralPDE.StochasticTraining(include_frac=0.9))
pde_system = PDESystem(eq,bcs,domains,[x,y,t],[u1])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1), progress = false; cb = cb, maxiters=3000)
phi = discretization.phi

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging


text:
"""
    DAsemiinf_DCS

Provides default parameters to use in the semi-infinite space autocorrelation function g2.

# Keyword Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_med`: medium's index of refraction
- `n_ext`: external medium's index of refraction (air or detector)
- `z`: the z-depth orthogonal from the boundary (cm)
- `β`: constant in Siegert relation dependent on collection optics
- `BFi`: Blood flow index ~αDb (cm²/s)
- `λ`: wavelength (nm)

# Examples
```
julia> data = DAsemiinf_DCS() # return default parameters
julia> data = DAsemiinf_DCS(ρ = 1.5) # return ρ = 1.5 with the rest of the parameters given by defaults

# we can now define our correlation times τ
julia> τ = 10 .^(range(-10,stop=0,length=250))
julia> g2_DA_semiinf_CW(τ, data) # can then simulate g2 in semiinf
```
"""
@with_kw struct DAsemiinf_DCS{T <: Real}
    μsp::T = 10.0                                       # reduced scattering coefficient (1/cm)
    μa::T = 0.1                                         # absorption coefficient (1/cm)
    n_ext::T = 1.0                                      # surrounding index of refraction
    n_med::T = 1.0                                      # layers index of refraction

    ρ::Union{T, AbstractVector{T}} = 1.0                # source-detector separation (cm)
    z::T = 0.0                                          # detector depth (cm)

    β::T = 1.0                                          # constant in Siegert relation dependent on collection optics
    BFi::T = 2.0e-8                                     # Blood flow index ~αDb (cm²/s)
    λ::T = 750.0                                        # wavelength (nm)
end

"""
    g2_DA_semiinf_CW(τ::AbstractVector, ρ, μa, μsp; BFi = 2e-8, β = 1.0, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)

Compute the electric field autocorrelation function function g2 in a semi-infinite geometry.

# Arguments
- `τ`: time lags (s)
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `BFi`: Blood flow index ~αDb (cm²/s)
- `β`: constant in Siegert relation dependent on collection optics
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `z`: the z-depth orthogonal from the boundary (cm)
- `λ`: wavelength (nm)

# Examples
```

julia> τ = 10 .^(range(-10,stop=0,length=250)) # need to define τ vector first
julia> g2_DA_semiinf_CW(τ, 1.0,0.1, 10.0) # can simulate with default parameters
julia> g2_DA_semiinf_CW(τ, 1.0,0.1, 10.0, BFi = 2.2e-8, β = 0.8) # manually define Keyword arguments
```
"""
function g2_DA_semiinf_CW(τ::AbstractVector, ρ, μa, μsp; BFi = 2e-8, β = 1.0, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)
    g2 = similar(τ)

    k0 = 2 * π / (λ * 1e-7) # this should be in cm
    tmp = 2 * μsp * BFi * k0^2

    G0 = fluence_DA_semiinf_CW(ρ, μa, μsp, n_ext = n_ext, n_med = n_med, z = z)

    for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        G1 = fluence_DA_semiinf_CW(ρ, μa_dynamic, μsp, n_ext = n_ext, n_med = n_med, z = z)
        g1 = G1 / G0
        g2[ind] = 1 + β * abs(g1)^2
    end

    return g2
end

"""
    g2_DA_semiinf_CW(τ::AbstractVector, data::DAsemiinf_DCS)

Provides a wrapper to `g2_DA_semiinf_CW(τ::AbstractVector, ρ, μa, μsp; BFi = 2e-8, β = 1.0, n_ext = 1.0, n_med = 1.0, z = 0.0, λ = 700.0)`.
Allows for arguments to be defined in a structure with `DAsemiinf_DCS`.

# Examples
```
julia> τ = 10 .^(range(-10,stop=0,length=250)) # need to define τ vector first
julia> data = DAsemiinf_DCS(ρ = 1.5, BFi = 3.1e-8) # return ρ = 1.5 with the rest of the parameters given by defaults
julia> g2_DA_semiinf_CW(τ, data) # can then simulate g2 in semiinf
```
"""
function g2_DA_semiinf_CW(τ::AbstractVector, data::DAsemiinf_DCS)
    return g2_DA_semiinf_CW(τ, data.ρ, data.μa, data.μsp, BFi = data.BFi, n_ext = data.n_ext, n_med = data.n_med, λ = data.λ)
end

"""
    Nlayer_cylinder_DCS

Provides default parameters to use to compute the autocorrelation function g2 in a N-layered cylinder.
The arguments defined with a Vector should be the same length and same type.

# Keyword Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa::Vector{T}`: absorption coefficient (cm⁻¹)
- `μsp::Vector{T}`: reduced scattering coefficient (cm⁻¹)
- `n_med::Vector{T}`: medium's index of refraction
- `n_ext`: external medium's index of refraction (air or detector)
- `l`: the thicknesses of each layer (cm)
- `z`: the z-depth orthogonal from the boundary (cm)
- `a`: the radius of the cylinder (cm)
- `β`: constant in Siegert relation dependent on collection optics
- `BFi::Vector{T}`: Blood flow index ~αDb (cm²/s)
- `λ`: wavelength (nm)
- `bessels::AbstractVector{T}`: roots of J0

# Examples
```
julia> data = Nlayer_cylinder_DCS() # return default parameters
julia> data = Nlayer_cylinder_DCS(ρ = 1.5) # return ρ = 1.5 with the rest of the parameters given by defaults

# we can now define our correlation times τ
julia> τ = 10 .^(range(-10,stop=0,length=250))
julia> g2_DA_Nlay_cylinder_CW(τ, data)
```
"""
@with_kw struct Nlayer_cylinder_DCS{T <: Real}
    μsp::Vector{T} = [10.0, 10.0]                       # reduced scattering coefficient (1/cm)
    μa::Vector{T} = [0.1, 0.1]                          # absorption coefficient (1/cm)
    n_ext::T = 1.0                                      # surrounding index of refraction
    n_med::Vector{T} = [1.0, 1.0]                       # layers index of refraction

    l::Vector{T} = [0.5, 10.0]                          # length of cylinder layers (cm)
    ρ::Union{T, AbstractVector{T}} = 1.0                # source-detector separation (cm)
    a::T = 20.0                                         # radius of cylinder (cm)
    z::T = 0.0                                          # detector depth (cm)

    β::T = 1.0                                          # constant in Siegert relation dependent on collection optics
    BFi::Vector{T} = [2.0e-8, 2.0e-8]                   # Blood flow index ~αDb (cm²/s)
    λ::T = 750.0                                        # wavelength (nm)

    bessels::AbstractVector{T} = besselroots[1:1000]    # roots of J0
end

"""
    g2_DA_Nlay_cylinder_CW(τ::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], β = 1.0, n_ext = 1.0, n_med = 1.0, l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000])

Compute the electric field autocorrelation function function g2 in a N-layered cylinder. μa, μsp, n_med and BFi must be vectors of the same type.

# Arguments
- `τ`: time lags (s)
- `ρ`: the source detector separation (cm⁻¹)
- `μa::Vector{T}`: absorption coefficient (cm⁻¹)
- `μsp::Vector{T}`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `n_med::Vector{T}`: medium's index of refraction
- `n_ext`: external medium's index of refraction (air or detector)
- `z`: the z-depth orthogonal from the boundary (cm)
- `l`: the thicknesses of each layer (cm)
- `a`: the radius of the cylinder (cm)
- `β`: constant in Siegert relation dependent on collection optics
- `BFi::Vector{T}`: Blood flow index ~αDb (cm²/s)
- `λ`: wavelength (nm)
- `bessels::AbstractVector{T}`: roots of J0

# Examples
```
julia> τ = 10 .^(range(-10,stop=0,length=250)) # need to define τ vector first
julia> g2_DA_Nlay_cylinder_CW(τ, 1.0, [0.1, 0.1], [10.0, 10.0], BFi = [2.1e-8, 3.0e-8]) # only define BFi
julia> g2_DA_Nlay_cylinder_CW(τ, 1.0, [0.1, 0.1], [10.0, 10.0], BFi = [2.1e-8, 3.0e-8], l = [0.8, 3.0]) # define BFi with layer thickness
```
"""
function g2_DA_Nlay_cylinder_CW(τ::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], β = 1.0, n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000])
    g2 = similar(τ)
    μa_dynamic = similar(μa)

    k0 = 2 * π / (λ * 1e-7) # this should be in cm
    tmp = @. 2 * μsp * BFi * k0^2

    G0 = fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, bessels)

    Threads.@threads for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        G1 = fluence_DA_Nlay_cylinder_CW(ρ, μa_dynamic, μsp, n_ext, n_med, l, a, z, bessels)
        g1 = G1 / G0
        g2[ind] = 1 + β * abs(g1)^2
    end

    return g2
end

"""
    g2_DA_Nlay_cylinder_CW(τ::AbstractVector, data::Nlayer_cylinder_DCS)

Provides a wrapper to `g2_DA_Nlay_cylinder_CW(τ::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], β = 1.0, n_ext = 1.0, n_med = 1.0, l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000])`.
Allows for arguments to be defined in a structure with `Nlayer_cylinder_DCS`. All vectors must be the same length (you will get a bounds error if not).

# Examples
```
julia> τ = 10 .^(range(-10,stop=0,length=250)) # need to define τ vector first
julia> data = Nlayer_cylinder_DCS(ρ = 1.5, BFi = [3.1e-8, 2.0e-7]) # return ρ = 1.5 and custom BFi
julia> g2_DA_Nlay_cylinder_CW(τ, data) # can then simulate g2 in semiinf
```
"""
function g2_DA_Nlay_cylinder_CW(τ::AbstractVector, data::Nlayer_cylinder_DCS)
    return g2_DA_Nlay_cylinder_CW(τ, data.ρ, data.μa, data.μsp, BFi = data.BFi, n_ext = data.n_ext, n_med = data.n_med, l = data.l, a = data.a, λ = data.λ, bessels = data.bessels)
end

#τ = 10 .^(range(-10,stop=0,length=250))
#BFi = 2e-8

#si = g2_DA_semiinf_CW(τ,1.0,0.1, 10.0, BFi = BFi)

#layered  = g2_DA_Nlay_cylinder_CW(τ, 1.0, [0.1, 0.1], [10.0, 10.0]; n_ext = 1.0, n_med = [1.0, 1.0], l = [0.5, 10.0], a = 20.0, z = 0.0, bessels = besselroots[1:1000])
#layered2  = g2_DA_Nlay_cylinder_CW(τ, 1.0, [0.1, 0.1], [10.0, 10.0]; BFi = [2e-8, 5e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [0.5, 10.0], a = 20.0, z = 0.0, bessels = besselroots[1:1000])
#layered3  = g2_DA_Nlay_cylinder_CW(τ, 1.0, [0.1, 0.1], [10.0, 10.0]; BFi = [5e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [0.5, 10.0], a = 20.0, z = 0.0, bessels = besselroots[1:1000])

#scatter(log10.(τ), si, label = "Semi-inf - BFi = 2e-8 cm²/s", ylabel = "g2(τ)", xlabel = "log(τ (s))")
#plot!(log10.(τ), layered, label = "Layered - BFi = [2e-8, 2e-8] cm²/s", lw = 1.5, linecolor =:black)
#plot!(log10.(τ), layered2, label = "Layered - BFi = [2e-8, 5e-8] cm²/s", linestyle=:dash, lw = 1.5, linecolor =:red)
#plot!(log10.(τ), layered3, label = "Layered - BFi = [5e-8, 2e-8] cm²/s", linestyle=:dot, lw = 1.5, linecolor =:blue4)

max_stars_repo_name: heltonmc/LightPropagation

max_issues_repo_name: heltonmc/LightPropagation

max_forks_repo_name: heltonmc/LightPropagation


text:
using SweepContractor

L=16; d=2

LTN = LabelledTensorNetwork{Tuple{Int,Int}}()
for i∈1:L, j∈1:L
    adj=Tuple{Int,Int}[];
    i>1 && push!(adj,(i-1,j))
    j>1 && push!(adj,(i,j-1))
    i<L && push!(adj,(i+1,j))
    j<L && push!(adj,(i,j+1))
    LTN[i,j] = Tensor(adj, randn(d*ones(Int,length(adj))...), i, j)
end

println("example_2d_open")
for χ ∈ d.^(4:9)
    sweep = sweep_contract(LTN, χ, 2*χ; fast=true)
    println("χ=$χ:\t", ldexp(sweep...))
end


text:
""" Dias de vida
# Leia a idade de uma pessoa expressa em anos, meses e dias
# e escreva a idade dessa pessoa expressa apenas em dias.
# Considerar ano com 365 dias e mês com 30 dias.
"""

println("informe sua idade")
days = ( parse(UInt8, readline()) * 365 )

println("informar quantos meses de vida ")
days += ( parse(UInt8, readline()) * 30 )

println("informar quantos dias de vida")
days += parse(UInt8, readline())

println("Total de ", days, " de vida")

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
# /!\ LazyBridgeOptimizer is still experimental and is meant to be used by JuMP. Use the SingleBridgeOptimizer in the solver tests as explained in the MOI documentation.
"""
    LazyBridgeOptimizer{OT<:MOI.ModelLike} <: AbstractBridgeOptimizer

The `LazyBridgeOptimizer` combines several bridges, which are added using the [`add_bridge`](@ref) function.
Whenever a constraint is added, it only attempts to bridge it if it is not supported by the internal model (hence its name `Lazy`).
When bridging a constraint, it selects the minimal number of bridges needed.
For instance, a constraint `F`-in-`S` can be bridged into a constraint `F1`-in-`S1` (supported by the internal model) using bridge 1 or
bridged into a constraint `F2`-in-`S2` (unsupported by the internal model) using bridge 2 which can then be
bridged into a constraint `F3`-in-`S3` (supported by the internal model) using bridge 3,
it will choose bridge 1 as it allows to bridge `F`-in-`S` using only one bridge instead of two if it uses bridge 2 and 3.
"""
mutable struct LazyBridgeOptimizer{OT<:MOI.ModelLike} <: AbstractBridgeOptimizer
    model::OT   # Internal model
    con_to_name::Dict{CI, String}
    name_to_con::Union{Dict{String, MOI.ConstraintIndex}, Nothing}
    # Constraint Index of bridged constraint -> Bridge.
    # It is set to `nothing` when the constraint is deleted.
    bridges::Vector{Union{Nothing, AbstractBridge}}
    # Constraint Index of bridged constraint -> Constraint type.
    constraint_types::Vector{Tuple{DataType, DataType}}
    # For `SingleVariable` constraints: (variable, set type) -> bridge
    single_variable_constraints::Dict{Tuple{Int64, DataType}, AbstractBridge}
    bridgetypes::Vector{Any} # List of types of available bridges
    dist::Dict{Tuple{DataType, DataType}, Int}      # (F, S) -> Number of bridges that need to be used for an `F`-in-`S` constraint
    best::Dict{Tuple{DataType, DataType}, DataType} # (F, S) -> Bridge to be used for an `F`-in-`S` constraint
end
function LazyBridgeOptimizer(model::MOI.ModelLike)
    return LazyBridgeOptimizer{typeof(model)}(
        model, Dict{CI, String}(), nothing,
        Union{Nothing, AbstractBridge}[], Tuple{DataType, DataType}[],
        Dict{Tuple{Int64, DataType}, AbstractBridge}(),
        Any[], Dict{Tuple{DataType, DataType}, Int}(),
        Dict{Tuple{DataType, DataType}, DataType}())
end

function _dist(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    if MOI.supports_constraint(b.model, F, S)
        return 0
    else
        return get(b.dist, (F, S), typemax(Int))
    end
end

# Update `b.dist` and `b.dest` for constraint types in `constraints`
function update_dist!(b::LazyBridgeOptimizer, constraints)
    # Bellman-Ford algorithm
    changed = true # Has b.dist changed in the last iteration ?
    while changed
        changed = false
        for BT in b.bridgetypes
            for (F, S) in constraints
                if MOI.supports_constraint(BT, F, S) && all(C -> supports_constraint_no_update(b, C[1], C[2]), added_constraint_types(BT, F, S))
                    # Number of bridges needed using BT
                    dist = 1 + mapreduce(
                        C -> _dist(b, C[1], C[2]), +,
                        added_constraint_types(BT, F, S), init = 0)
                    # Is it better that what can currently be done ?
                    if dist < _dist(b, F, S)
                        b.dist[(F, S)] = dist
                        b.best[(F, S)] = concrete_bridge_type(BT, F, S)
                        changed = true
                    end
                end
            end
        end
    end
end

function fill_required_constraints!(required::Set{Tuple{DataType, DataType}}, b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    if supports_constraint_no_update(b, F, S)
        return # The constraint is supported
    end
    if (F, S) in required
        return # The requirements for this constraint have already been added or are being added
    end
    # The constraint is not supported yet, add in `required` the required constraint types to bridge it
    push!(required, (F, S))
    for BT in b.bridgetypes
        if MOI.supports_constraint(BT, F, S)
            for C in added_constraint_types(BT, F, S)
                fill_required_constraints!(required, b, C[1], C[2])
            end
        end
    end
end

# Compute dist[(F, S)] and best[(F, S)]
function update_constraint!(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    required = Set{Tuple{DataType, DataType}}()
    fill_required_constraints!(required, b, F, S)
    update_dist!(b, required)
end

"""
    add_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})

Enable the use of the bridges of type `BT` by `b`.
"""
function add_bridge(b::LazyBridgeOptimizer, BT::Type{<:AbstractBridge})
    push!(b.bridgetypes, BT)
    # Some constraints (F, S) in keys(b.best) may now be bridged
    # with a less briges than `b.dist[(F, S)] using `BT`
    update_dist!(b, keys(b.best))
end

# It only bridges when the constraint is not supporting, hence the name "Lazy"
function is_bridged(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    !MOI.supports_constraint(b.model, F, S)
end
# Same as supports_constraint but do not trigger `update_constraint!`. This is
# used inside `update_constraint!`.
function supports_constraint_no_update(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    return MOI.supports_constraint(b.model, F, S) || (F, S) in keys(b.best)
end
function supports_bridging_constraint(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    update_constraint!(b, F, S)
    return (F, S) in keys(b.best)
end
function bridge_type(b::LazyBridgeOptimizer{BT}, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet}) where BT
    update_constraint!(b, F, S)
    result = get(b.best, (F, S), nothing)
    if result === nothing
        throw(MOI.UnsupportedConstraint{F, S}())
    end
    return result
end


text:
#===============================================================================

    mc_analysis.jl

    perform melting curve analysis

    Author: Tom Price
    Date:   July 2019

===============================================================================#

import DataStructures.OrderedDict
import DataFrames.DataFrame
import Memento: debug

#===============================================================================
    field names >>
===============================================================================#

# const MC_TF_KEYS = [:temperature, :fluorescence]
const MC_OUTPUT_FIELDS = OrderedDict(
    :observed_data          => :melt_curve_data,
    :peaks_filtered         => :melt_curve_analysis)

#===============================================================================
    function definitions >>
===============================================================================#

"Analyse a melting curve experiment via calls to `calibrate` and `mc_peak_analysis`."
function mc_analysis(i ::McInput)

    # ## function: format fluorescence data for calibration
    # ##
    # ## >> PROBLEM >>
    # ## this function generates a lot of intermediate representations
    # ## and might be sped up by creating an appropriate container
    # ## and mutating it in place
    # function get_mc_data(channel_num ::Integer)
    #
    #     ## subset melting curve data by channel (curried)
    #     select_mcdata_by_channel(channel_num ::Integer) =
    #         mc_data ::DataFrame ->
    #             Dict(
    #                 map([:temperature, :fluorescence, :well]) do f
    #                     f => mc_data[f][mc_data[:channel] .== channel_num]
    #                 end)
    #
    #     ## split temperature and fluorescence data by well
    #     ## return vector of TF Dicts
    #     split_tf_by_well(fluo_sel ::Associative) =
    #         map(i.wells) do well
    #             Dict(
    #                 map(MC_TF_KEYS) do key
    #                     key => fluo_sel[key][fluo_sel[:well] .== well]
    #                 end)
    #         end
    #
    # "Curried function that extends a vector with NaN values to a specified length."
    # extend_NaN(len ::Integer) =
    #     vec ::AbstractVector ->
    #         len - length(vec) |>
    #             m ->
    #                 m >= 0 ?
    #                     vcat(vec, fill(NaN_T, m)) :
    #                     error(logger, "vector is too long")
    #
    #     ## extend array elements with NaNs to length of longest element
    #     @inline extend(x ::AbstractArray) =
    #         map(extend_NaN(y |> mold(length) |> maximum), x)
    #
    #     ## extend data vectors with NaN values where necessary to make them equal in length
    #     ## this is performed to convert the fluorescence data to a 3D array
    #     ##
    #     ## >> PROBLEM >>
    #     ## although it is unlikely that the maximum vector length varies by channel
    #     ## if it does the arrays (which are generated separately by channel)
    #     ## will not be conformable and the data transformation will throw an error
    #     extend_tf_vecs(tf_dict_vec ::AbstractArray) =
    #         map(tf_dict_vec) do tf_dict
    #             Dict(
    #                 map(MC_TF_KEYS) do key
    #                     key => extend_NaN(
    #                                 maximum(
    #                                     map(length ∘ index(:temperature),
    #                                         tf_dict_vec)))(tf_dict[key])
    #                 end)
    #         end
    #
    #     ## convert to MeltCurveTF object
    #     toMeltCurveTF(tf_nv_adj ::AbstractArray) =
    #         MeltCurveTF(
    #             map(MC_TF_KEYS) do key
    #                 mapreduce(index(key), hcat, tf_nv_adj)
    #             end...)
    #
    # ## << end of function definitions nested in get_mc_data()
    #
    # ## calculate
    #     return i.raw_df |>
    #             select_mcdata_by_channel(channel_num) |>
    #             split_tf_by_well |>
    #             extend_tf_vecs |>
    #             toMeltCurveTF
    # end ## get_mc_data

    "Convert raw melting curve data from a DataFrame to 3D array of fluorescences
    as required by `calibrate`."
    function transform_3d(df ::DataFrame)
        ## split-apply-combine style
        # extended = copy(df)
        # extended[:fluorescence] = Vector{Float_T}(extended[:fluorescence]) ## to allow NaN
        # const lengths = by(extended, [:channel, :well],
        #     df -> DataFrame(len=length(df[:temperature])))
        # const longest = maximum(lengths[:len])
        # for j in 1:size(lengths, 1),
        #     k in range(1, longest - lengths[:len][j])
        #     push!(extended, Dict(
        #         :temperature    => NaN_T,
        #         :fluorescence   => NaN_T,
        #         :channel        => lengths[:channel][j],
        #         :well           => lengths[:well][k]))
        # end ## next j, k
        # const f =
        #     reshape(
        #         hcat(by(extended, [:channel, :well], df -> df[:fluorescence])[:x1]...)',
        #         i.num_channels, i.num_wells, longest) |>
        #     RawData{Float_T}
        #
        ## devectorized style
        cs = [df[:channel] .== c for c in i.channels]
        ws = [df[:well   ] .== w for w in i.wells   ]
        selection =
            [   find(cs[ci] .& ws[wi])
                for ci in eachindex(i.channels), wi in eachindex(i.wells)   ]
        # @inline selector(c ::Int, w ::Int) = (df[:channel] .== c) .& (df[:well] .== w)
        # const selection =
        #     [   find(selector(c, w))
        #         for c in i.channels, w in i.wells   ]
        longest = selection |> mold(length) |> maximum
        dims = (longest, i.num_wells, i.num_channels)
        t = Array{Float_T}(dims)
        f = Array{Float_T}(dims)
        for ci in eachindex(i.channels)
            # c = i.channels[ci]
            for wi in eachindex(i.wells)
                # w = i.wells[wi]
                location = selection[ci, wi]
                let ## to future-proof access to `ti` outside the for-loop:
                    ## scoping rules change after v0.6
                    ti = 0
                    for ti in eachindex(location)
                        t[ti, wi, ci] = df[location[ti], :temperature ]
                        f[ti, wi, ci] = df[location[ti], :fluorescence]
                    end ## next ti
                    for ti in (ti + 1):longest
                        t[ti, wi, ci] = NaN_T
                        f[ti, wi, ci] = NaN_T
                    end ## next ti
                end ## let
            end ## next wi
        end ## next ci
        return RawData(t), RawData(f)
    end ## transform_3d()

    normalize_tf(ci ::Int_T, wi ::Int_T) =
        normalize_fluos!(
            remove_when_temperature_NaN(
                DataFrame(
                    :temperature    => raw_temps.data[ :, wi, ci],
                    :fluorescence   => calibrated_data[:, wi, ci])))

    "Take as input a vector of temperatures and a vector of fluorescences, and set
    the fluorescence to NaN wherever the corresponding temperature value is NaN."
    remove_when_temperature_NaN(df ::DataFrame) =
        df[.!isnan.(df[:temperature]), :]

    "Normalize fluorescences values by subtracting the lowest value. Note that
    if any value is NaN, the result will be a vector of NaNs."
    function normalize_fluos!(df ::DataFrame)
        df[:fluorescence] = sweep(minimum)(-)(df[:fluorescence])
        return df
    end

    ## << end of function definitions nested in mc_analysis()

    debug(logger, "at mc_analysis()")
    # const (channels, wells) =
    #     map((:channel, :well)) do fieldname
    #         i.raw_df[fieldname] |> unique |> sort
    #     end ## do fieldname
    # const num_channels      = length(channels)
    # const num_wells         = length(wells)
    #
    ## get data arrays by channel
    ## output is Vector{MeltCurveTF}
    # const mc_data_bychannel = map(get_mc_data, i.channels)
    #
    ## reshape raw fluorescence data to 3-dimensional array
    ## dimensions 1,2,3 = temperature,well,channel
    # const mc_data_array     = cat(3, map(field(:fluorescence), mc_data_bychannel)...)
    # const raw_fluos         = RawData(mc_data_array)

    ## alternative method, manipulating DataFrame directly
    (raw_temps, raw_fluos) = transform_3d(i.raw_df)
    #
    ## deconvolute and normalize
    peak_format = peak_output_format(i.out_format)
    if peak_format == McPeakLongOutput
        (   background_subtracted_data,
            k_deconv,
            deconvoluted_data,
            norm_data,
            norm_wells,
            calibrated_data ) =
                calibrate(i, i.calibration_data, i.calibration_args, raw_fluos, array)
    else
        ## McPeakShortOutput
        (_, _, _, _, norm_wells, calibrated_data) = ## discard other output fields
            calibrate(i, i.calibration_data, i.calibration_args, raw_fluos, array)
    end ## if peak_format
    #
    ## subset temperature/fluorescence data by channel / well
    ## then smooth the fluorescence/temperature data and calculate Tm peak, area
    ## result is mc_matrix: dim1 = well, dim2 = channel
    mcpa_matrix =
        eachindex(i.channels) |> ## do for each channel
        moose(hcat) do ci
            map(eachindex(i.wells)) do wi
                if i.wells[wi] in norm_wells
                    mc_peak_analysis(peak_format, i, normalize_tf(ci, wi))
                else
                    McPeakShortOutput()
                end ## if
            end ## next wi
        end #= next ci =# |>
        morph(length(norm_wells), i.num_channels) ## coerce to 2d array
    #
    # if (i.out_format == full_output)
    if peak_format == McPeakLongOutput
        return McLongOutput(
            i.wells,
            i.channels,
            raw_fluos.data,
            background_subtracted_data,
            k_deconv,
            deconvoluted_data,
            norm_data,
            i.wells,
            calibrated_data,
            mcpa_matrix)
    else
        ## json_output, pre_json_output
        ## McPeakShortOutput
        output_dict = OrderedDict{Symbol,Any}(
            map(keys(MC_OUTPUT_FIELDS)) do f
                MC_OUTPUT_FIELDS[f] =>
                    [   i.reporting(getfield(mcpa_matrix[wi, ci], f))
                        for wi in 1:i.num_wells, ci in 1:i.num_channels ]
            end) ## do f
        output_dict[:valid] = true
        return output_dict
    end ## if out_format
end ## mc_analysis()

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
function RIC2ECI3x3(RIC, r, v)
    y = LinearAlgebra.cross(r, v)
    r̂ = normalize(r)
    ĉ = normalize(y)
    î = LinearAlgebra.cross(ĉ, r̂)
    rotate(RIC, r̂, î, ĉ)
end

@inline function combinecov(Σ₁, Σ₂, y, z)
    Σ = Σ₁ + Σ₂
    x = LinearAlgebra.cross(y, z)
    a = Σ * x
    SymmetricM2(
        a' * x,
        a' * z,
        quadform(z, Σ)
    )
end

### Is this actually faster than just multiplying matrices???
@inline function rotate(RIC, r̂, î, ĉ)

    @fastmath @inbounds begin
        a = r̂[1]*RIC[1,1] + î[1]*RIC[1,2] + ĉ[1]*RIC[1,3]
        b = r̂[1]*RIC[2,1] + î[1]*RIC[2,2] + ĉ[1]*RIC[2,3]
        c = r̂[1]*RIC[3,1] + î[1]*RIC[3,2] + ĉ[1]*RIC[3,3]

        ECI11 = a * r̂[1] + b * î[1] + c * ĉ[1]
        ECI12 = a * r̂[2] + b * î[2] + c * ĉ[2]
        ECI13 = a * r̂[3] + b * î[3] + c * ĉ[3]

        a = r̂[2]*RIC[1,1] + î[2]*RIC[1,2] + ĉ[2]*RIC[1,3]
        b = r̂[2]*RIC[2,1] + î[2]*RIC[2,2] + ĉ[2]*RIC[2,3]
        c = r̂[2]*RIC[3,1] + î[2]*RIC[3,2] + ĉ[2]*RIC[3,3]

        ECI22 = a * r̂[2] + b * î[2] + c * ĉ[2]
        ECI23 = a * r̂[3] + b * î[3] + c * ĉ[3]

        a = r̂[3]^2*RIC[1,1] + 2r̂[3]*î[3]*RIC[1,2] + 2r̂[3]*ĉ[3]*RIC[1,3]
        b = î[3]^2*RIC[2,2] + 2î[3]*ĉ[3]*RIC[2,3] + ĉ[3]^2*RIC[3,3]

        ECI33 = a + b
    end

    SymmetricM3( ECI11, ECI12, ECI22, ECI13, ECI23, ECI33 )
end

@inline function RIC_to_2D(r₁::SVector{3,T}, v₁::SVector{3,T}, Σ₁::SymmetricM3{T},
                            r₂::SVector{3,T}, v₂::SVector{3,T}, Σ₂::SymmetricM3{T}) where {T}
    δr = r₁ - r₂
    δr₀ = norm(δr)
    δv = v₁ - v₂
    z = LinearAlgebra.cross(δr, δv)

    # ! Relative encounter frame
    nδv = normalize(δv)
    nz  = normalize(z)

    Cp = combinecov(RIC2ECI3x3(Σ₁, r₁, v₁), RIC2ECI3x3(Σ₂, r₂, v₂), nδv, nz)
    Cp, δr₀
end


text:
function ackley(x; a=20, b=0.2, c=2π)
    d = length(x)
    return -a*exp(-b*sqrt(sum(x.^2)/d)) -
        exp(sum(cos.(c*x))/d) + a + exp(1)
end

booth(x) = (x[1] + 2*x[2] - 7)^2 + (2*x[1] + x[2] -5)^2

function branin(x; a=1, b=5.1/(4π^2), c=5/π, r=6, s=10, t=1/(8π))
    return a*(x[2] - b*x[1]^2 + c*x[1] - r)^2 + s*(1-t)*cos(x[1]) + s
end

function flower(x; a=1, b=1, c=4)
    return a*norm(x) + b*sin(c*atan(x[2], x[1]))
end

function michalewicz(x; m=10)
    return -sum(sin(xᵢ)*sin(i*xᵢ^2/π)^(2m) for (i, xᵢ) in enumerate(x))
end

rosenbrock(x; a=1, b=100) = (a - x[1])^2 + b*(x[2] - x[1]^2)^2

wheeler(x, a=1.5) = -exp(-(x[1]*x[2] - a)^2 - (x[2] - a)^2)

function circle(x)
    θ = x[1]
    r = 0.5 + 0.5 * (2*x[2]/(1 + x[2]^2))
    return [1-r*cos(θ); 1-r*sin(θ)]
end