text:
module PooksoftAssetModelingKit

# include -
include("Include.jl")

# export methods -
export price
export evaluate
export sample_underlying_gbm_model
export estimate_single_index_model
export estimate_random_walk_model
export compute_linear_return_array
export compute_log_return_array
export compute_return_volatility

# lattice -
export compute_ternary_price_tree
export compute_binary_price_tree

# technical indicators methods -
export compute_fibonacci_retracement_levels
export compute_exponential_moving_average

# factory methods -
export build_geometric_brownian_motion_model

# export types -
export PSGeometricBrownianMotionModelParameters
export PSHestonAssetPricingModelParameters
export PSSingleIndexModelParameters
export PSRandomWalkModelParameters

end # module


text:
ABCS  = zeros(4)
AA    = zeros(2,2)
BB    = zeros(2,2)
CC    = zeros(2,2)
SS    = zeros(2)

svars    = tuple(ones(11)...)
fvars    = tuple(ones(37)...)
sdvars   = tuple([ZeroPotential(); ones(70)]...)
bdgvars  = tuple(ones(71)...)
phidvars = tuple([ZeroPotential(); ones(72)]...)
avars    = tuple([ZeroPotential(); ones(76)]...)

@testset "S equation inner grid coefficients:" begin
    AdS5_3_1.S_eq_coeff!(ABCS, svars, AdS5_3_1.Inner())

    @test all(ABCS .≈ [6.0, 48.0, 133.42988060987634, 130.85976121975267])
end

@testset "Fxy equations inner grid coefficients:" begin
    AdS5_3_1.Fxy_eq_coeff!(AA, BB, CC, SS, fvars, AdS5_3_1.Inner())

    @test all( AA .≈ [48.92907291226281 0.0; 0.0 18.0] )
    @test all( BB .≈ [811.6226411252798 -151.92523101186953; 119.40115727043062 41.42071634074196] )
    @test all( CC .≈ [9249.882276773546 -1318.4601133442727; 2930.339069731014 256.1897346300231] )
    @test all( SS .≈ [9577.64621761626, 3981.7746343640238])
end

@testset "Sd equation inner grid coefficients:" begin
    AdS5_3_1.Sd_eq_coeff!(ABCS, sdvars, AdS5_3_1.Inner())

    @test all( ABCS .≈ [0.0, -880.7233124207306, -2544.311791437666, 273918.94409251545])
end

@testset "B2d equation inner grid coefficients:" begin
    AdS5_3_1.B2d_eq_coeff!(ABCS, bdgvars, AdS5_3_1.Inner())

    @test all( ABCS .≈ [0.0, -2642.169937262192, -9687.956436628037, -23245.124131500335] )
end

@testset "B1dGd equations inner grid coefficients:" begin
    AdS5_3_1.B1dGd_eq_coeff!(AA, BB, CC, SS, bdgvars, AdS5_3_1.Inner())

    @test all( AA .≈ [0.0 0.0; 0.0 0.0] )
    @test all( BB .≈ [-2642.169937262192 0.0; 0.0 -2642.169937262192] )
    @test all( CC .≈ [-15724.739986410723 -6036.783549782685; 14374.17230438983 -9687.956436628037] )
    @test all( SS .≈ [-143418.3032738444, -58153.507064913196] )
end

@testset "phid equations inner grid coefficients:" begin
    AdS5_3_1.phid_eq_coeff!(ABCS, phidvars, AdS5_3_1.Inner())

    @test all( ABCS .≈ [0.0, -587.1488749471532, -1565.7303331924084, 12974.144452661105] )
end

@testset "A equations inner grid coefficients:" begin
    AdS5_3_1.A_eq_coeff!(ABCS, avars, AdS5_3_1.Inner())

    @test all( ABCS .≈ [1321.084968631096, 7926.5098117865755, 7926.5098117865755, -827587.7882250317 ])
end

ext: jl


text:
#Comparison tools for NetwAnalysis
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const SPARAM_THRESH = 1e-12
const TPARAM_THRESH = 1e-12
const ABCD_THRESH = 1e-12

#==Error of NetworkParameters
===============================================================================#
_assert_z0match(n1::NetworkParameters, n2::NetworkParameters) = nothing
function _assert_z0match(n1::NetworkParametersRef, n2::NetworkParametersRef)
	if n1.z0 != n2.z0
		throw(ArgumentError("Type mismatch: NetworkParameters"))
	end
	nothing
end

function maxerror(v1::DataF1, v2::DataF1)
	@assert(v1.x == v2.x, "x-values do not match")
	Δ = v2-v1
	ϵ = maximum(abs(Δ))
end

function maxerror(n1::NetworkParameters, n2::NetworkParameters)
	throw(ArgumentError("Type mismatch: NetworkParameters"))
end

#==Error of NetworkParameters
===============================================================================#
#For matrix of scalars
#-------------------------------------------------------------------------------
function _maxerror(::Type{Number}, n1::NetworkParameters, n2::NetworkParameters)
	_assert_z0match(n1, n2)
	Δ = n2-n1
	ϵ = maximum(abs.(Δ.m))
	return ϵ
end

function maxerror(n1::NetworkParametersRef{TP, NP, T},
	n2::NetworkParametersRef{TP, NP, T2}) where {TP, NP, T<:Number, T2<:Number}
	return _maxerror(Number, n1, n2)
end
function maxerror(n1::NetworkParametersNoRef{TP, NP, T},
	n2::NetworkParametersNoRef{TP, NP, T2}) where {TP, NP, T<:Number, T2<:Number}
	return _maxerror(Number, n1, n2)
end

#For matrix of Vectors
#-------------------------------------------------------------------------------
function _maxerror(::Type{Vector}, n1::NetworkParameters, n2::NetworkParameters)
	_assert_z0match(n1, n2)
	Δv = n2-n1
	Δ = [maximum(abs2.(elem)) for elem in Δv.m]
	ϵ = maximum(Δ)
	return ϵ
end
function maxerror(n1::NetworkParametersRef{TP, NP, T},
	n2::NetworkParametersRef{TP, NP, T}) where {TP, NP, T<:Vector}
	return _maxerror(Vector, n1, n2)
end
function maxerror(n1::NetworkParametersNoRef{TP, NP, T},
	n2::NetworkParametersNoRef{TP, NP, T}) where {TP, NP, T<:Vector}
	return _maxerror(Vector, n1, n2)
end

#For matrix of DataF1
#-------------------------------------------------------------------------------
function _maxerror(::Type{DataF1}, n1::NetworkParameters, n2::NetworkParameters)
	_assert_z0match(n1, n2)
	Δv = n2-n1
	Δ = [maximum(abs2(elem)) for elem in Δv.m]
	ϵ = maximum(Δ)
	return ϵ
end
function maxerror(n1::NetworkParametersRef{TP, NP, T},
	n2::NetworkParametersRef{TP, NP, T2}) where {TP, NP, T<:DataF1, T2<:DataF1}
	return _maxerror(DataF1, n1, n2)
end
function maxerror(n1::NetworkParametersNoRef{TP, NP, T},
	n2::NetworkParametersNoRef{TP, NP, T2}) where {TP, NP, T<:DataF1, T2<:DataF1}
	return _maxerror(DataF1, n1, n2)
end

#Last line


text:
# # [Grover Search](@id Grover)
using Yao
using YaoExtensions: variational_circuit
using LinearAlgebra

# ## Grover Step
# A single grover step is consist of applying oracle circuit and reflection circuit.
# The `reflection_circuit` function takes the wave function generator `U` as the input and returns `U|0><0|U'`.
function grover_step!(reg::AbstractRegister, oracle, U::AbstractBlock)
    apply!(reg |> oracle, reflect_circuit(U))
end

function reflect_circuit(gen::AbstractBlock{N}) where N
    reflect0 = control(N, -collect(1:N-1), N=>-Z)
    chain(gen', reflect0, gen)
end

# Compute the propotion of target states to estimate the number of iterations,
# which requires computing the output state.
function solution_state(oracle, gen::AbstractBlock{N}) where N
    reg= zero_state(N) |> gen
    reg.state[real.(statevec(ArrayReg(ones(ComplexF64, 1<<N)) |> oracle)) .> 0] .= 0
    normalize!(reg)
end

function num_grover_step(oracle, gen::AbstractBlock{N}) where N
    reg = zero_state(N) |> gen
    ratio = abs2(solution_state(oracle, gen)'*reg)
    Int(round(pi/4/sqrt(ratio)))-1
end

# #### Run
# First, we define the problem by an oracle, it finds bit string `bit"000001100100"`.
num_bit = 12
oracle = matblock(Diagonal((v = ones(ComplexF64, 1<<num_bit); v[Int(bit"000001100100")+1]*=-1; v)))

# then solve the above problem
gen = repeat(num_bit, H, 1:num_bit)
reg = zero_state(num_bit) |> gen

target_state = solution_state(oracle, gen)

for i = 1:num_grover_step(oracle, gen)
    grover_step!(reg, oracle, gen)
    overlap = abs(reg'*target_state)
    println("step $(i-1), overlap = $overlap")
end

# ## Rejection Sampling
# In practise, it is often not possible to determine the number of iterations before actual running.
# we can use rejection sampling technique to avoid estimating the number of grover steps.

using Random; Random.seed!(2)  #src

# In a single try, we `apply` the grover algorithm for `nstep` times.
function single_try(oracle, gen::AbstractBlock{N}, nstep::Int; nbatch::Int) where N
    reg = zero_state(N+1; nbatch=nshot)
    focus!(reg, 1:N) do r
        r |> gen
        for i = 1:nstep
            grover_step!(r, oracle, gen)
        end
        return r
    end
    reg |> checker
    res = measure!(RemoveMeasured(), reg, (N+1))
    return res, reg
end

# After running the grover search, we have a checker program that flips the ancilla qubit
# if the output is the desired value, we assume the checker program can be implemented in polynomial time.
# to gaurante the output is correct.
# We contruct a checker "program", if the result is correct, flip the ancilla qubit
ctrl = -collect(1:num_bit); ctrl[[3,6,7]] *= -1
checker = control(num_bit+1,ctrl, num_bit+1=>X)

# The register is batched, with batch dimension `nshot`.
# [`focus!`](@ref Yao.focus!) views the first 1-N qubts as system.
# For a batched register, [`measure!`](@ref Yao.measure!)
# returns a vector of bitstring as output.

# #### Run
maxtry = 100
nshot = 3

for nstep = 0:maxtry
    println("number of iter = $nstep")
    res, regi = single_try(oracle, gen, nstep; nbatch=3)

    ## success!
    if any(==(1), res)
        overlap_final = viewbatch(regi, findfirst(==(1), res))'*target_state
        println("success, overlap = $(overlap_final)")
        break
    end
end

# The final state has an overlap of `1` with the target state.

# ## Amplitude Amplification
# Given a circuit to generate a state,
# now we want to project out the subspace with [1,3,5,8,9,11,12] fixed to 1 and [4,6] fixed to 0.
# We can construct an oracle
evidense = [1, 3, -4, 5, -6, 8, 9, 11, 12]
function inference_oracle(nbit::Int, locs::Vector{Int})
    control(nbit, locs[1:end-1], abs(locs[end]) => (locs[end]>0 ? Z : -Z))
end
oracle = inference_oracle(nqubits(reg), evidense)

# We use a variational circuit generator defined in `YaoExtensions`
gen = dispatch!(variational_circuit(num_bit), :random)
reg = zero_state(num_bit) |> gen

# #### Run
solution = solution_state(oracle, gen)
for i = 1:num_grover_step(oracle, gen)
    grover_step!(reg, oracle, gen)
    println("step $(i-1), overlap = $(abs(reg'*solution))")
end

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
# basic test that parsing works correctly

testdir = dirname(@__FILE__)

@test_throws Gumbo.InvalidHTMLException parsehtml("", strict=true)

let
    page = open("$testdir/example.html") do example
        example |> readall |> parsehtml
    end
    @test page.doctype == "html"
    root = page.root
    @test tag(root[1][1]) == :meta
    @test root[2][1][1].text == "A simple test page."
    @test is(root[2][1][1].parent, root[2][1])
end


text:
using Mustache
using Test

mutable struct CallableFunction <: Function
  calls::Int
end
(F::CallableFunction)() = (F.calls += 1; string(F.calls))

@testset " lambdas " begin
# overview: |
#   Lambdas are a special-cased data type for use in interpolations and
#   sections.

#   When used as the data value for an Interpolation tag, the lambda MUST be
#   treatable as an arity 0 function, and invoked as such.  The returned value
#   MUST be rendered against the default delimiters, then interpolated in place
#   of the lambda.

#   When used as the data value for a Section tag, the lambda MUST be treatable
#   as an arity 1 function, and invoked as such (passing a String containing the
#   unprocessed section contents).  The returned value MUST be rendered against
#   the current delimiters, then interpolated in place of the section.
# tests:
#   - name: Interpolation
#     desc: A lambda's return value should be interpolated.
#     data:
#       lambda: !code
#         ruby:    'proc { "world" }'
#         perl:    'sub { "world" }'
#         js:      'function() { return "world" }'
#         php:     'return "world";'
#         python:  'lambda: "world"'
#         clojure: '(fn [] "world")'
#     template: "Hello, {{lambda}}!"
#     expected: "Hello, world!"

    template = "Hello, {{lambda}}!"
    expected = "Hello, world!"
    data = Dict("lambda"=> () -> "world")
    @test Mustache.render(template, data) == expected

#   - name: Interpolation - Expansion
#     desc: A lambda's return value should be parsed.
#     data:
#       planet: "world"
#       lambda: !code
#         ruby:    'proc { "{{planet}}" }'
#         perl:    'sub { "{{planet}}" }'
#         js:      'function() { return "{{planet}}" }'
#         php:     'return "{{planet}}";'
#         python:  'lambda: "{{planet}}"'
#         clojure: '(fn [] "{{planet}}")'
#     template: "Hello, {{lambda}}!"
#     expected: "Hello, world!"

    template = "Hello, {{lambda}}!"
    expected = "Hello, world!"
    data = Dict("lambda" => () -> "{{planet}}", "planet" => "world")
    @test Mustache.render(template, data) == expected

#   - name: Interpolation - Alternate Delimiters
#     desc: A lambda's return value should parse with the default delimiters.
#     data:
#       planet: "world"
#       lambda: !code
#         ruby:    'proc { "|planet| => {{planet}}" }'
#         perl:    'sub { "|planet| => {{planet}}" }'
#         js:      'function() { return "|planet| => {{planet}}" }'
#         php:     'return "|planet| => {{planet}}";'
#         python:  'lambda: "|planet| => {{planet}}"'
#         clojure: '(fn [] "|planet| => {{planet}}")'
#     template: "{{= | | =}}\nHello, (|&lambda|)!"
#     expected: "Hello, (|planet| => world)!"

    template = "{{= | | =}}\nHello, (|&lambda|)!"
    expected = "Hello, (|planet| => world)!"
    data = Dict("planet"=>"world", "lambda"=> () -> "|planet| => {{planet}}")
    @test Mustache.render(template, data) == expected

#   - name: Interpolation - Multiple Calls
#     desc: Interpolated lambdas should not be cached.
#     data:
#       lambda: !code
#         ruby:    'proc { $calls ||= 0; $calls += 1 }'
#         perl:    'sub { no strict; $calls += 1 }'
#         js:      'function() { return (g=(function(){return this})()).calls=(g.calls||0)+1 }'
#         php:     'global $calls; return ++$calls;'
#         python:  'lambda: globals().update(calls=globals().get("calls",0)+1) or calls'
#         clojure: '(def g (atom 0)) (fn [] (swap! g inc))'
#     template: '{{lambda}} == {{{lambda}}} == {{lambda}}'
#     expected: '1 == 2 == 3'

    template = "{{lambda}} == {{{lambda}}} == {{lambda}}"
    expected =  "1 == 2 == 3"
    data = Dict("lambda" => CallableFunction(0))
    @test Mustache.render(template, data) == expected

#   - name: Escaping
#     desc: Lambda results should be appropriately escaped.
#     data:
#       lambda: !code
#         ruby:    'proc { ">" }'
#         perl:    'sub { ">" }'
#         js:      'function() { return ">" }'
#         php:     'return ">";'
#         python:  'lambda: ">"'
#         clojure: '(fn [] ">")'
#     template: "<{{lambda}}{{{lambda}}}"
#     expected: "<>>"

    template = "<{{lambda}}{{{lambda}}}"
    expected = "<>>"
    data = Dict("lambda" => () -> ">")
    @test Mustache.render(template, data) == expected

#   - name: Section
#     desc: Lambdas used for sections should receive the raw section string.
#     data:
#       x: 'Error!'
#       lambda: !code
#         ruby:    'proc { |text| text == "{{x}}" ? "yes" : "no" }'
#         perl:    'sub { $_[0] eq "{{x}}" ? "yes" : "no" }'
#         js:      'function(txt) { return (txt == "{{x}}" ? "yes" : "no") }'
#         php:     'return ($text == "{{x}}") ? "yes" : "no";'
#         python:  'lambda text: text == "{{x}}" and "yes" or "no"'
#         clojure: '(fn [text] (if (= text "{{x}}") "yes" "no"))'
#     template: "<{{#lambda}}{{x}}{{/lambda}}>"
#     expected: "<yes>"

    template = "<{{#lambda}}{{x}}{{/lambda}}>"
    expected = "<yes>"
    data = Dict("x" => "Error!", "lambda" => (txt) ->  txt == "{{x}}" ? "yes" : "no")
    @test Mustache.render(template, data) == expected

#   - name: Section - Expansion
#     desc: Lambdas used for sections should have their results parsed.
#     data:
#       planet: "Earth"
#       lambda: !code
#         ruby:    'proc { |text| "#{text}{{planet}}#{text}" }'
#         perl:    'sub { $_[0] . "{{planet}}" . $_[0] }'
#         js:      'function(txt) { return txt + "{{planet}}" + txt }'
#         php:     'return $text . "{{planet}}" . $text;'
#         python:  'lambda text: "%s{{planet}}%s" % (text, text)'
#         clojure: '(fn [text] (str text "{{planet}}" text))'
#     template: "<{{#lambda}}-{{/lambda}}>"
#     expected: "<-Earth->"

    template = "<{{#lambda}}-{{/lambda}}>"
    expected = "<-Earth->"
    data = Dict("planet"=> "Earth", "lambda" => (txt) -> txt * "{{planet}}" * txt)
    @test Mustache.render(template, data) == expected

#   - name: Section - Alternate Delimiters
#     desc: Lambdas used for sections should parse with the current delimiters.
#     data:
#       planet: "Earth"
#       lambda: !code
#         ruby:    'proc { |text| "#{text}{{planet}} => |planet|#{text}" }'
#         perl:    'sub { $_[0] . "{{planet}} => |planet|" . $_[0] }'
#         js:      'function(txt) { return txt + "{{planet}} => |planet|" + txt }'
#         php:     'return $text . "{{planet}} => |planet|" . $text;'
#         python:  'lambda text: "%s{{planet}} => |planet|%s" % (text, text)'
#         clojure: '(fn [text] (str text "{{planet}} => |planet|" text))'
#     template: "{{= | | =}}<|#lambda|-|/lambda|>"
#     expected: "<-{{planet}} => Earth->"

template = "{{= | | =}}<|#lambda|-|/lambda|>"
expected =  "<-{{planet}} => Earth->"
data = Dict("planet"=>"Earth", "lambda"=> (txt) -> txt * "{{planet}} => |planet|" * txt)
@test Mustache.render(template, data) == expected

#   - name: Section - Multiple Calls
#     desc: Lambdas used for sections should not be cached.
#     data:
#       lambda: !code
#         ruby:    'proc { |text| "__#{text}__" }'
#         perl:    'sub { "__" . $_[0] . "__" }'
#         js:      'function(txt) { return "__" + txt + "__" }'
#         php:     'return "__" . $text . "__";'
#         python:  'lambda text: "__%s__" % (text)'
#         clojure: '(fn [text] (str "__" text "__"))'
#     template: '{{#lambda}}FILE{{/lambda}} != {{#lambda}}LINE{{/lambda}}'
#     expected: '__FILE__ != __LINE__'

 template = "{{#lambda}}FILE{{/lambda}} != {{#lambda}}LINE{{/lambda}}"
expected = "__FILE__ != __LINE__"
data = Dict("lambda" => (txt) ->  "__" * txt * "__")
@test Mustache.render(template, data) == expected

#   - name: Inverted Section
#     desc: Lambdas used for inverted sections should be considered truthy.
#     data:
#       static: 'static'
#       lambda: !code
#         ruby:    'proc { |text| false }'
#         perl:    'sub { 0 }'
#         js:      'function(txt) { return false }'
#         php:     'return false;'
#         python:  'lambda text: 0'
#         clojure: '(fn [text] false)'
#     template: "<{{^lambda}}{{static}}{{/lambda}}>"
#     expected: "<>"

 template = "<{{^lambda}}{{static}}{{/lambda}}>"
expected = "<>"
data = Dict("static" => "static", "lambda" => (txt) -> false)
@test Mustache.render(template, data) == expected

end


text:
"""
The OutputChannel.jl module defines the OutputChannel type and methods

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Exports

# ------ Types

export OutputChannel

# ------ Functions

export publish, get_last_value

# --- Type definitions

using ZMQ

"""
    struct OutputChannel

Component that encodes and publishes output data. The `data` field stores the
data most recently sent by the channel.
"""
struct OutputChannel
    #=
      Fields
      ------
      * `data`: most recently sent data

      * `url`: URL that OutputChannel publishes output data on

      * `socket`: ZMQ Socket for sending data
    =#
    data::AbstractChannelData
    url::String
    socket::Socket
end

"""
    OutputChannel(data_type::Type{<:AbstractChannelData}, url::String;
                 use_bind=false)

Construct an output channel that publishes data of type `data_type` to `url`.
When `use_bind` is true, the ZMQ Socket that messages are published to is
connected to `url` using the `bind()` method. Otherwise, the ZMQ Socket is
connected to `url` using the `connect()` method.
"""
function OutputChannel(data_type::Type{<:AbstractChannelData},
                       url::String;
                       use_bind=true)

    # Create data storage
    data = data_type()

    # Create Socket to publish output data
    socket = Socket(PUB)

    # Connect socket to URL
    if use_bind
        bind(socket, url)
    else
        connect(socket, url)
    end

    # Return new OutputChannel
    OutputChannel(data, url, socket)
end

# --- Method definitions

"""
    publish(channel::OutputChannel, value)

Publish `value` to `channel`.
"""
function publish(channel::OutputChannel, value)
    # Encode value and send message
    message = encode_value(typeof(channel.data), value)
    send(channel.socket, message)

    # Update channel data
    set_data!(channel.data, value, decode=false)
end

"""
    get_last_value(channel::OutputChannel)

Return the most recent data value sent by `channel`.
"""
get_last_value(channel::OutputChannel) = get_data(channel.data)


text:
module Dec5Part1
function row_to_int(row::String)::Int64
    binary_array = map(x -> x == 'F' ? '0' : '1', row)
    return binary_to_int(binary_array)
end

function col_to_int(col::String)::Int64
    binary_array = map(x -> x == 'L' ? '0' : '1', col)
    return binary_to_int(binary_array)
end

function binary_to_int(x::String)::Int64
    return parse(Int64, x, base=2)
end

function get_seat(row::Int64, col::Int64)::Int64
    return row * 8 + col
end

function get_codes(x::String)::Tuple{String,String}
    return x[begin:7], x[8:end]
end
function get_seat_number(code::String)::Int64
    row, col = get_codes(code)
    row_id = row_to_int(row)
    col_id = col_to_int(col)
    return get_seat(row_id, col_id)
end

function main(x::String)::Int64
    boardingpasses = readlines(x)
    seat_numbers = get_seat_number.(boardingpasses)
    return max(seat_numbers...)
end

end


text:
struct FieldNull <: Field end
struct Field1{E,I}  <: Field
    inner::I
    es::Vector{E}
    Field1(es::Vector{E}) where E = new{E,Nothing}(nothing,es)
    Field1(inner::F,es::Vector{E}) where {E,F} = new{E,typeof(inner)}(inner,es)
    Field1(::Nothing,es::Vector{E}) where E = new{E,Nothing}(nothing,es)
end
struct IndexedExpression{E <: Expression} <: Entry
    index::Int
    e::E
end

FIELD_NULL = FieldNull()

@inline (F::FieldNull)(y,x,p=nothing) = nothing
@inline function (F::Field1{E,I})(y,x,p=nothing) where {E,I}
    inner(F)(y,x,p)
    @simd for i in eachindex(F.es)
        F.es[i](y,x,p)
    end
end
@inline function (F::Field1{E,Nothing})(y,x,p=nothing) where E
    @simd for i in eachindex(F.es)
        F.es[i](y,x,p)
    end
end
@inline (e::IndexedExpression{E})(y,x,p=nothing) where {E} = (y[e.index] = e.e(x,p))

Field() = Sink{Field}(FIELD_NULL)
# Field(e::E) where E <: IndexedExpression= Field1([e])
Field(f::FieldNull,e::E) where E <: IndexedExpression = Field1([e])
Field(f::F,e::E) where {F <: Field, E <: IndexedExpression} = _field(f,e) ? f : Field1(f,[e])

function _field(F::Field1{E,I},e) where {E,I}
    if e isa eltype(F.es)
        push!(F.es,e)
        return true
    else
        return _field(inner(F),e)
    end
end
function _field(F::Field1{E,Nothing},e) where E
    if e isa eltype(F.es)
        push!(F.es,e)
        return true
    else
        return false
    end
end


text:
function get_data(dataobject::ClumpDataType,
                vars::Array{Symbol,1},
                units::Array{Symbol,1},
                direction::Symbol,
                center::Array{<:Any,1},
                mask::MaskType,
                ref_time::Real)

    vars_dict = Dict()
    #vars = unique(vars)
    boxlen = dataobject.boxlen
    descriptor =[:x, :y, :z, :mass]

    if direction == :z
        apos = :peak_x
        bpos = :peak_y
        cpos = :peak_z

        avel = :vx
        bvel = :vy
        cvel = :vz

    elseif direction == :y
        apos = :peak_z
        bpos = :peak_x
        cpos = :peak_y

        avel = :vz
        bvel = :vx
        cvel = :vy
    elseif direction == :x
        apos = :peak_z
        bpos = :peak_y
        cpos = :peak_x

        avel = :vz
        bvel = :vy
        cvel = :vx
    end

    column_names = propertynames(dataobject.data.columns)

    for i in vars

        # quantities that are in the datatable
        if in(i, column_names) || in(i, descriptor)#|| occursin("var", string(i))
            selected_unit = getunit(dataobject, i, vars, units)

            if i == :peak_x || i == :x
                vars_dict[i] = ( select(dataobject.data, apos) .-  boxlen * center[1]) .* selected_unit
            elseif i == :peak_y || i == :y
                vars_dict[i] = ( select(dataobject.data, bpos) .-  boxlen * center[2]) .* selected_unit
            elseif i == :peak_z || i == :z
                vars_dict[i] = ( select(dataobject.data, cpos) .-  boxlen * center[3]) .* selected_unit

            elseif i == :vx
                vars_dict[i] =  select(dataobject.data, avel) .* selected_unit
            elseif i == :vy
                vars_dict[i] =  select(dataobject.data, bvel) .* selected_unit
            elseif i == :vz
                vars_dict[i] =  select(dataobject.data, cvel) .* selected_unit
            elseif i == :mass
                vars_dict[i] =  select(dataobject.data, :mass_cl) .* selected_unit
            else
                vars_dict[i] =  select(dataobject.data, i) .* selected_unit
            end

        # quantities that are derived from the variables in the data table
        elseif i == :v
            selected_unit = getunit(dataobject, :v, vars, units)
            vars_dict[:v] =  sqrt.(select(dataobject.data, :vx).^2 .+
                                   select(dataobject.data, :vy).^2 .+
                                   select(dataobject.data, :vz).^2 ) .* selected_unit

        elseif i == :ekin
            selected_unit = getunit(dataobject, :ekin, vars, units)
            vars_dict[:ekin] =   0.5 .* getvar(dataobject, vars=[:mass_cl])  .*
                                (select(dataobject.data, :vx).^2 .+
                                select(dataobject.data, :vy).^2 .+
                                select(dataobject.data, :vz).^2 ) .* selected_unit
        end

    end

    if length(mask) > 1
        for i in keys(vars_dict)
            vars_dict[i]=vars_dict[i][mask]
        end
    end

    if length(vars)==1
            return vars_dict[vars[1]]
    else
            return vars_dict
    end

end


text:
# given a simplex and a face returns:
# +v if face is the v-th face of the simplex oriented according to the simplex
# -v if face is the v-th face of the simplex oriented oppositely to the simplex
# 0 is face is not a face of the simplex
function relorientation(face, simplex)

    v = setdiff(simplex, face)
    length(v) == 1 || return 0

    # find the position of the missing vertex
    v = v[1]
    i = Base.findfirst(simplex, v)
    s = (-1)^(i-1)

    # remove that vertex from the simplex
    face2 = Array{Int}(length(simplex)-1)
    for j in 1 : i-1
        face2[j] = simplex[j]
    end
    for j in i : length(simplex)-1
        face2[j] = simplex[j+1]
    end

    # get the permutation that maps face to face2
    #p = indexin(face, face2)
    p = [ findfirst(face2,v) for v in face ]

    return s * levicivita(p) * i
end

"""
    getcommonedge(cell1, cell2) -> e1, e2, edge

Returns in edge the common vertices of cell1 and cell2. e1 contains the index
of the vertex of cell1 opposite to this common edge, and with a plus or minus
sign depending on whether the orientation of the common edge is along or
against the internal orientation of cell1. Similar for e2.
"""
function getcommonedge(cell1, cell2)
    isct = intersect(cell1, cell2)
    relorientation(isct, cell1), relorientation(isct, cell2), isct
end


text:
export flush_print

function flush_print_string(s::String, width::Int, right::Bool)::String
    n = length(s)
    @assert n >= 0 "Width must be nonnegative, got $n"

    if width < n
        @warn "Trunctated to fit width"
        return s[1:width]
    end

    spacer = " "^(width - n)
    if right
        return spacer * s
    end
    return s * spacer
end

"""
`flush_print(x,width)` returns a `String` version of `x` right justified
in a string of length `width`.

Use `flush_print(x,width,false)` for left-justified.
"""
function flush_print(x, width::Int, right::Bool = true)::String
    return flush_print_string(string(x), width, right)
end


text:
using Test, AdmittanceModels
using LinearAlgebra: eigen, rank, nullspace
using SparseArrays: sparse, SparseMatrixCSC

@testset "closest_permutation" begin
    mat = [[1, 0, 2, 0, 1.5] [1, 0, 3, 0, 2] [1, 0, 1.3, 0, 1.4]]
    u = AdmittanceModels.closest_permutation(mat)
    @test vcat(sum(u, dims=1)...) == [1,1,1]
    _, indices = findmax(u, dims=1)
    @test [index[1] for index in vcat(indices...)] == [5,3,1]
end

@testset "inv_power_eigen" begin
    M = [1im 2.2 4; 3.1 0.2 3; 4 1 2im]
    F = eigen(M)
    for i in 1:size(M, 1)
        v0, λ0 = F.vectors[:,i], F.values[i]
        λ, v = AdmittanceModels.inv_power_eigen(M, v0=v0)
        @test λ ≈ λ0
        λ, v = AdmittanceModels.inv_power_eigen(M, λ0=λ0)
        @test λ ≈ λ0
        λ, v = AdmittanceModels.inv_power_eigen(M, v0=v0, λ0=λ0)
        @test λ ≈ λ0
    end
end

@testset "row_column_graph" begin
    rows = [[1.0, 2.0, -3.0,    0, 0, 0,  0],
        [    0,     0,  1.0, -4.0, 0, 0,  0],
        [    0.1,   0,    0, -0.3, 0, 0,  0],
        [    0,     0,    0,    0, 1, 1,  0],
        [    0,     0,    0,    0, 0, 1, -1]]
    mat = sparse(transpose(hcat(rows...)))
    correct_adj = zeros(Int, size(rows, 1)+1, size(rows, 1)+1)
    correct_adj[1, 2] = correct_adj[2, 1] = 3
    correct_adj[1, 3] = correct_adj[3, 1] = 1
    correct_adj[1, 6] = correct_adj[6, 1] = 2
    correct_adj[2, 3] = correct_adj[3, 2] = 4
    correct_adj[4, 5] = correct_adj[5, 4] = 6
    correct_adj[4, 6] = correct_adj[6, 4] = 5
    correct_adj[5, 6] = correct_adj[6, 5] = 7
    adj = AdmittanceModels.row_column_graph(mat)
    @test adj == correct_adj
end

@testset "sparse_nullbasis and nullbasis" begin
    # example with no redundancies or inconsistentices
    rows = [[1.0, 2.0, -3.0,    0, 0, 0,  0],
            [  0,   0,  1.0, -4.0, 0, 0,  0],
            [0.1,   0,    0, -0.3, 0, 0,  0],
            [  0,   0,    0,    0, 1, 1,  0],
            [  0,   0,    0,    0, 0, 1, -1]]
    mat = sparse(transpose(hcat(rows...)))
    null_basis = AdmittanceModels.sparse_nullbasis(mat)
    @test count(!iszero, mat * null_basis) == 0
    @test rank(collect(mat)) + rank(collect(null_basis)) == size(mat, 2) # rank nullity theorem
    @test count(!iszero, nullspace(collect(mat))) == 14 # dense
    @test count(!iszero, null_basis) == 7 # sparse
    @test AdmittanceModels.nullbasis(mat) == null_basis
    # example with redundancies and inconsistentices
    rows = [[1, -1,  0, 0,  0,  0],
            [1,  0, -1, 0,  0,  0],
            [0,  1,  1, 0,  0,  0],
            [0,  0,  0, 1, -1,  0],
            [0,  0,  0, 1,  0, -1],
            [0,  0,  0, 0,  1, -1]]
    mat = sparse(transpose(hcat(rows...)))
    null_basis = AdmittanceModels.sparse_nullbasis(mat)
    @test count(!iszero, mat * null_basis) == 0
    @test rank(collect(mat)) + rank(collect(null_basis)) == size(mat, 2) # rank nullity theorem
    @test AdmittanceModels.sparse_nullbasis(sparse(mat)) == null_basis
    @test AdmittanceModels.nullbasis(mat) == null_basis
    # example where sparse nullbasis algorithm does not apply
    mat = reshape(collect(1:9), 3,3)
    @test_throws AssertionError AdmittanceModels.sparse_nullbasis(sparse(mat))
    @test AdmittanceModels.nullbasis(mat, warn=false) == nullspace(mat)
end


text:
function bc_rect!(nxe,nye,nf,dir)
  #
  # This subroutine generates the nf array for a rectangular mesh
  # of 8-node quadrilaterals fully fixed on the base and with
  # vertical rollers on the left and right sides. Nodes numbered
  # in the x- or y-direction
  #
 if dir == :y
   nm = 0
   ic = 0
   for j in 1:2nye
     nm = nm+1
     nf[1, nm] = 0
     ic = ic + 1
     nf[2, nm]=ic
   end
   nm += 1
   nf[1, nm] = 0
   nf[2, nm] = 0

   for i in 1:nxe-1
     for j in 1:nye
       nm = nm + 1
       ic = ic + 1
       nf[1, nm] = ic
       ic += 1
       nf[2, nm] = ic
     end
     nm += 1
     nf[1, nm] = 0
     nf[2, nm] = 0

     for j in 1:2nye
       nm += 1
       ic += 1
       nf[1, nm] = ic
       ic += 1
       nf[2, nm] = ic
     end
     nm += 1
     nf[1, nm] = 0
     nf[2, nm] = 0
   end

   for j in 1:nye
     nm += 1
     ic += 1
     nf[1, nm] = ic
     ic += 1
     nf[2, nm] = ic
   end
   nm += 1
   nf[1, nm] = 0
   nf[2, nm] = 0

   for j in 1:2nye
     nm += 1
     nf[1, nm] = 0
     ic += 1
     nf[2, nm] = ic
   end
   nm += 1
   nf[1, nm] = 0
   nf[2, nm] = 0
 else
   nm=0
   ic=0
   for j in 1:nye
     nm += 1
     nf[1, nm] = 0
     ic += 1
     nf[2, nm] = ic
     for i in 1:2nxe-1
       nm += 1
       ic += 1
       nf[1, nm] = ic
       ic += 1
       nf[2, nm] = ic
     end
     nm += 1
     nf[1, nm] = 0
     ic += 1
     nf[2, nm] = ic

     nm += 1
     nf[1, nm] = 0
     ic += 1
     nf[2, nm] = ic
     for i in 1:nxe-1
       nm += 1
       ic += 1
       nf[1, nm] = ic
       ic += 1
       nf[2, nm] = ic
     end
     nm += 1
     nf[1, nm] = 0
     ic += 1
     nf[2, nm] = ic
   end

   for i in 1:2*nxe+1
     nm += 1
     nf[1, nm] = 0
     nf[2, nm] = 0
   end
 end
end


text:
module TestBasic

using Test
using TreeParzen

# Quadratic TPE
objective(params) = (params[:x] - 3)^2
best = fmin(objective, Dict(:x => HP.Uniform(:x, -5.0, 5.0)), 50)
@test abs(best[:x] - 3.0) < .25

# Generate trials to calculate
objective(params) = params[:x]^2 + params[:y]^2
space = Dict(
    :x => HP.Uniform(:x, -10.0, 10.0),
    :y => HP.Uniform(:y, -10.0, 10.0),
)
points = [Dict(:x => 0.0, :y => 0.0), Dict(:x => 1.0, :y => 1.0)]
best = fmin(objective, space, 10, points)
@test best == Dict(:x => 0.0, :y => 0.0)

end
true


text:
##########################################################################################
#
#    Definition of distributions
#
#    links logpdfXXXX() functions to logpdf(XXXX(), x) in Distributions.jl
#
##########################################################################################

########## locally defined distributions #############
# defined here because there is a speed gain over using Distributions.jl
# TODO : file an issue in Distributions.jl

# function logpdfNormal(mu::Real, sigma::Real, x::Real)
# 	local const fac = -log(sqrt(2pi))
# 	@assert sigma > 0. "calling logpdfNormal with negative or null stdev"
# 	local r = (x-mu)/sigma
# 	return -r*r*0.5-log(sigma)+fac
# end

# function logpdfUniform(a::Real, b::Real, x::Real)
# 	@assert a < b "calling logpdfUniform with lower limit >= upper limit"
# 	return (a <= x <= b) ? -log(b-a) : throw("give up eval")
# end

########### distributions using libRmath #########

#          Name           Arity
dists = { (:Weibull, 	  3),
		  (:Binomial, 	  3),
		  (:Gamma,  	  3),
		  (:Cauchy,  	  3),
		  (:LogNormal,    3),
		  (:Beta, 	      3),
		  (:Poisson,  	  2),
		  (:TDist,  	  2),
		  (:Exponential,  2),
		  (:Uniform, 	  3),
	      (:Normal,  	  3),
		  (:Bernoulli,    2)}

for d in dists  # d = dists[11]
	fsym = symbol("logpdf$(d[1])")

	npars = d[2]
	args = [symbol("x$i") for i in [npars, 1:(npars-1)...]] # in correct order for libRmath

	fex = quote
		function ($fsym)($([Expr(:(::), symbol("x$i"), :Real) for i in 1:npars]...))
			local res = 0.

	        res = logpdf( ($(d[1]))($([symbol("x$i") for i in 1:(npars-1)]...)), $(symbol("x$npars")) )

			if res == -Inf
				throw("give up eval")
			elseif isnan(res)
				local ar = $(Expr(:tuple, [args..., 1]))
				error(string("calling ", $(string(fsym)), ar, " returned an error"))
			end
			res
		end
	end

	eval(fex)
end

####### all possible vectorized versions  ##########

for d in dists # d = dists[1]
	fsym = symbol("logpdf$(d[1])")

	arity = d[2]
	ps = [ ifloor((i-1) / 2^(j-1)) % 2 for i=2:2^arity, j=1:arity]

	for l in 1:size(ps,1) # l = 4
		pars = Symbol[ps[l,j]==0 ? :Real : :AbstractArray for j in 1:arity]

		rv = symbol("x$(findfirst(pars .== :AbstractArray))")
		npars = length(pars)

		fex = quote
			function ($fsym)($([Expr(:(::), symbol("x$i"), pars[i]) for i in 1:npars]...))
		        local res = 0.
		        for i in 1:length($rv)
		        	res += ($fsym)($([pars[i]==:Real ? symbol("x$i") : Expr(:ref, symbol("x$i"), :i) for i in 1:npars]...))
		        end
		        res
		    end
		end

		eval(fex)
	end
end


text:
function happy(x)
	happy_ints = ref(Int)
	int_try = 1
	while length(happy_ints) < x
		n = int_try
		past = ref(Int)
		while n != 1
	        	n = sum([y^2 for y in digits(n)])
	        	contains(past,n) ? break : push!(past,n)
	    	end
		n == 1 && push!(happy_ints,int_try)
		int_try += 1
	end
	return happy_ints
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ArgParse

push!(LOAD_PATH, pwd())

using Generator

# Main entry
function main(sequence, time, silence_time, noise, out)

	generate_signal(sequence,
			time=time,
			silence_time=silence_time,
			noise=noise,
			out_file=out)
end

if PROGRAM_FILE == "SignalGenerator.jl"

	# Parse the arguments
	s = ArgParseSettings("Audio dialing generator")

	@add_arg_table s begin
		"--sequence"
			help = "The sequence to be generated (ie: *1_5_0_#)" *
			       "The '_' represent the silences"
			required = true
		"--time"
			help = "The tone duration time in [ms]"
			default = 90
			arg_type = Int
		"--silence-time"
			help = "The silence duration time in [ms]"
			default = 70
			arg_type = Int
		"--noise"
			help = "Add white noise to the signal"
			default = false
			arg_type = Bool
		"--out"
			help = "Output file name with the file extension"
			default = "tone.wav"
	end

	parsed_args = parse_args(ARGS, s)

	main(parsed_args["sequence"],
	     parsed_args["time"],
	     parsed_args["silence-time"],
	     parsed_args["noise"],
	     parsed_args["out"])
end

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
module DiffusionProfile
# This module calculates the corrrelation function <δz(t)δ(0)> and its itegration in time
# for local diffusion coefficient from time series of single observable.

#include
include("my_functions.jl")

#using
using DelimitedFiles
using Statistics
using Plots
using BenchmarkTools

Nwins = 101
z0 = -15.0
dz = 0.3
#Cycle through all umbrella windows to extract D(z)
for win in 1:1
    #import umbrella window time series, calculate <δz(t)δ(0)>.
    obj = readdlm("/Users/seb/Desktop/CO2O2HCO3_graphene_project/wham/HCO3/O2_charged_window$win.dat")
    (corr,corr_cut,var0,tCorr) = corr_i(obj,-15.0+(win-1)*dz)

    #Integrate correlation function from 0 → "∞"
    int_corr = Nint_corr(corr_cut)
    println("Window $win:")
    println("< δz(0)δz(0) > = ",var0," Å^2")
    println("∫< δz(t)δz(0) >dt = ",int_corr," Å^2 ⋅ ps")

    D_i = (var0^2/int_corr)*1e-4
    println("Local diffusion coefficient = ",D_i," cm^2/s")
    println(" ")

    #visualize <δz(t)δ(0)>
    l = length(corr_cut[:,1])
    dash0 = zeros(l,2)
    display(plot(corr_cut[:,1],[corr_cut[:,2],dash0[:,2]],color = [:gray :red],linewidth = [4 1.2],xlabel = "t (ps)",ylabel = "< δz(t)δz(0) >",label = ["< δz(t)δz(0) >"  " 0 "], ylims = (-0.017,0.03),linestyle = [:solid :dashdot]))

    savefig("Corr_O2_charged_win1.png")
end

#@btime corr_i(obj,-15.0)

end


text:
""" 2-category of categories, functors, and natural transformations.

Categories in mathematics appear in the large, often as categories of sets with
extra structure, and in the small, as algebraic structures that generalize
groups, monoids, preorders, and graphs. This division manifests in Catlab as
well. Large categories (in spirit, if not in the [technical
sense](https://ncatlab.org/nlab/show/large+category)) occur throughout Catlab as
`@instance`s of the theory of categories. For computational reasons, small
categories are usually presented by generators and relations.

This module provides a minimal interface to accomodate both situations. Category
instances are supported through the wrapper type [`TypeCat`](@ref). Finitely
presented categories are provided by another module, [`FinCats`](@ref).
"""
module Categories
export Cat, TypeCat, Functor, Transformation, dom, codom, compose, id,
  ob, hom, is_hom_equal, ob_map, hom_map, dom_ob, codom_ob, component,
  OppositeCat, op, co

using AutoHashEquals

using ...GAT
import ...Theories: Category2, ob, hom, dom, codom, compose, ⋅, ∘, id,
  composeH, *

# Categories
############

""" Size of a category, used for dispatch and subtyping purposes.

A [`Cat`](@ref) type having a particular `CatSize` means that categories of that
type are *at most* that large.
"""
abstract type CatSize end

""" Size of a large category, such as Set.

To the extent that they form a category, we regard Julia types and functions
([`TypeCat`](@ref)) as forming a large category.
"""
struct LargeCatSize <: CatSize end

""" Abstract base type for a category.

The objects and morphisms in the category have Julia types `Ob` and `Hom`,
respectively. Note that these types do *not* necessarily form an `@instance` of
the theory of categories, as they may not meaningfully form a category outside
the context of this object. For example, a finite category regarded as a
reflexive graph with a composition operation might have type `Cat{Int,Int}`,
where the objects and morphisms are numerical identifiers for vertices and edges
in the graph.

The basic operations available in any category are: [`dom`](@ref),
[`codom`](@ref), [`id`](@ref), [`compose`](@ref).
"""
abstract type Cat{Ob,Hom,Size<:CatSize} end

""" Coerce or look up object in category.

Converts the input to an object in the category, which should be of type `Ob` in
a category of type `Cat{Ob,Hom}`. How this works depends on the category, but a
common case is to look up objects, which might be integers or GAT expressions,
by their human-readable name, usually a symbol.

See also: [`hom`](@ref).
"""
function ob end

""" Coerce or look up morphism in category.

See also: [`ob`](@ref).
"""
function hom end

""" Domain of morphism in category.
"""
dom(C::Cat, f) = dom(f)

""" Codomain of morphism in category.
"""
codom(C::Cat, f) = codom(f)

""" Identity morphism on object in category.
"""
id(C::Cat, x) = id(x)

""" Compose morphisms in a category.
"""
compose(C::Cat, fs...) = compose(fs...)

""" Are two morphisms in a category equal?

By default, just checks for equality of Julia objects using ``==``. In some
categories, checking equality of morphisms may involve nontrivial reasoning.
"""
is_hom_equal(C::Cat, f, g) = is_hom_equal(f, g)
is_hom_equal(f, g) = f == g

# Instances
#----------

""" Pair of Julia types regarded as a category.

The Julia types should form an `@instance` of the theory of categories
(`Theories.Category`).
"""
struct TypeCat{Ob,Hom} <: Cat{Ob,Hom,LargeCatSize} end

TypeCat(Ob::Type, Hom::Type) = TypeCat{Ob,Hom}()

# FIXME: Type conversion isn't practical because types are often too tight.
#ob(::TypeCat{Ob,Hom}, x) where {Ob,Hom} = convert(Ob, x)
#hom(::TypeCat{Ob,Hom}, f) where {Ob,Hom} = convert(Hom, f)
ob(::TypeCat, x) = x
hom(::TypeCat, f) = f

Base.show(io::IO, ::TypeCat{Ob,Hom}) where {Ob,Hom} =
  print(io, "TypeCat($Ob, $Hom)")

# Functors
##########

""" Abstract base type for a functor between categories.

A functor has a domain and a codomain ([`dom`](@ref) and [`codom`](@ref)), which
are categories, and object and morphism maps, which can be evaluated using
[`ob_map`](@ref) and [`hom_map`](@ref). The functor object can also be called
directly when the objects and morphisms have distinct Julia types. This is
sometimes but not always the case (see [`Cat`](@ref)), so when writing generic
code one should prefer the `ob_map` and `hom_map` functions.
"""
abstract type Functor{Dom<:Cat,Codom<:Cat} end

""" Evaluate functor on object.
"""
@inline ob_map(F::Functor, x) = do_ob_map(F, x)

""" Evaluate functor on morphism.
"""
@inline hom_map(F::Functor, f) = do_hom_map(F, f)

""" Identity functor on a category.
"""
@auto_hash_equals struct IdentityFunctor{Dom<:Cat} <: Functor{Dom,Dom}
  dom::Dom
end

codom(F::IdentityFunctor) = F.dom

do_ob_map(F::IdentityFunctor, x) = ob(F.dom, x)
do_hom_map(F::IdentityFunctor, f) = hom(F.dom, f)

function Base.show(io::IO, F::IdentityFunctor)
  print(io, "id(")
  show_domains(io, F, codomain=false)
  print(io, ")")
end

""" Composite of functors.
"""
@auto_hash_equals struct CompositeFunctor{Dom,Codom,
    F<:Functor{Dom,<:Cat},G<:Functor{<:Cat,Codom}} <: Functor{Dom,Codom}
  fst::F
  snd::G
end
Base.first(F::CompositeFunctor) = F.fst
Base.last(F::CompositeFunctor) = F.snd

dom(F::CompositeFunctor) = dom(first(F))
codom(F::CompositeFunctor) = codom(last(F))

do_ob_map(F::CompositeFunctor, x) = ob_map(F.snd, ob_map(F.fst, x))
do_hom_map(F::CompositeFunctor, f) = hom_map(F.snd, hom_map(F.fst, f))

function Base.show(io::IO, F::CompositeFunctor)
  print(io, "compose(")
  show(io, first(F))
  print(io, ", ")
  show(io, last(F))
  print(io, ")")
end

show_type_constructor(io::IO, ::Type{<:Functor}) = print(io, "Functor")

function show_domains(io::IO, f; domain::Bool=true, codomain::Bool=true,
                      recurse::Bool=true)
  if get(io, :hide_domains, false)
    print(io, "…")
  else
    if domain
      show(IOContext(io, :compact=>true, :hide_domains=>!recurse), dom(f))
    end
    if codomain
      domain && print(io, ", ")
      show(IOContext(io, :compact=>true, :hide_domains=>!recurse), codom(f))
    end
  end
end

# Callables
#----------

Functor(f::Function, g::Function, C::Cat, D::Cat) = FunctorCallable(f, g, C, D)

""" Functor defined by two Julia callables, an object map and a morphism map.
"""
@auto_hash_equals struct FunctorCallable{Dom,Codom} <: Functor{Dom,Codom}
  ob_map::Any
  hom_map::Any
  dom::Dom
  codom::Codom
end

dom(F::FunctorCallable) = F.dom
codom(F::FunctorCallable) = F.codom
do_ob_map(F::FunctorCallable, x) = F.ob_map(x)
do_hom_map(F::FunctorCallable, f) = F.hom_map(f)

# Instances
#----------

(F::Functor{TypeCat{Ob,Hom}})(x::Ob) where {Ob,Hom} = ob_map(F, x)
(F::Functor{TypeCat{Ob,Hom}})(f::Hom) where {Ob,Hom} = hom_map(F, f)

# Natural transformations
#########################

""" Abstract base type for a natural transformation between functors.

A natural transformation ``α: F ⇒ G`` has a domain ``F`` and codomain ``G``
([`dom`](@ref) and [`codom`](@ref)), which are functors ``F,G: C → D`` having
the same domain ``C`` and codomain ``D``. The transformation consists of a
component ``αₓ: Fx → Gx`` in ``D`` for each object ``x ∈ C``, accessible using
[`component`](@ref) or indexing notation (`Base.getindex`).
"""
abstract type Transformation{C<:Cat,D<:Cat,Dom<:Functor,Codom<:Functor} end

""" Component of natural transformation.
"""
function component end

@inline Base.getindex(α::Transformation, c) = component(α, c)

""" Domain object of natural transformation.

Given ``α: F ⇒ G: C → D``, this function returns ``C``.
"""
dom_ob(α::Transformation) = dom(dom(α)) # == dom(codom(α))

""" Codomain object of natural transformation.

Given ``α: F ⇒ G: C → D``, this function returns ``D``.
"""
codom_ob(α::Transformation) = codom(dom(α)) # == codom(codom(α))

@auto_hash_equals struct IdentityTransformation{C<:Cat,D<:Cat,Dom<:Functor{C,D}} <:
     Transformation{C,D,Dom,Dom}
  dom::Dom
end

codom(α::IdentityTransformation) = α.dom

function component(α::IdentityTransformation, x)
  F = dom(α)
  id(codom(F), ob_map(F, x))
end

const IdIdTransformation{C<:Cat} = IdentityTransformation{C,C,IdentityFunctor{C}}

# 2-category of categories
##########################

@instance Category2{Cat,Functor,Transformation} begin
  dom(F::Functor) = F.dom
  codom(F::Functor) = F.codom
  id(C::Cat) = IdentityFunctor(C)

  function compose(F::Functor, G::Functor; strict::Bool=true)
    !strict || codom(F) == dom(G) ||
      error("Domain mismatch in composition $F ⋅ $G")
    compose_id(F, G)
  end

  dom(α::Transformation) = α.dom
  codom(α::Transformation) = α.codom
  id(F::Functor) = IdentityTransformation(F)

  function compose(α::Transformation, β::Transformation; strict::Bool=true)
    !strict || codom(α) == dom(β) ||
      error("Domain mismatch in vertical composition $α ⋅ $β")
    compose_id(α, β)
  end
  function composeH(α::Transformation, β::Transformation; strict::Bool=true)
    !strict || codom_ob(α) == dom_ob(β) ||
      error("Domain mismatch in horizontal composition $α * $β")
    composeH_id(α, β)
  end
  function composeH(α::Transformation, H::Functor; strict::Bool=true)
    !strict || codom_ob(α) == dom(H) ||
      error("Domain mismatch in whiskering $α * $H")
    composeH_id(α, H)
  end
  function composeH(F::Functor, β::Transformation; strict::Bool=true)
    !strict || codom(F) == dom_ob(β) ||
      error("Domain mismatch in whiskering $F * $β")
    composeH_id(F, β)
  end
end

# XXX: Is this normalization of identities using multiple dispatch a good idea?
# Unlike in `Sets`, it doesn't feel great since it requires so much boilerplate.

@inline compose_id(F::Functor, G::Functor) = do_compose(F, G)
@inline compose_id(F::Functor, ::IdentityFunctor) = F
@inline compose_id(::IdentityFunctor, G::Functor) = G
@inline compose_id(F::IdentityFunctor, ::IdentityFunctor) = F

@inline compose_id(α::Transformation, β::Transformation) = do_compose(α, β)
@inline compose_id(α::Transformation, ::IdentityTransformation) = α
@inline compose_id(::IdentityTransformation, β::Transformation) = β
@inline compose_id(α::IdentityTransformation, ::IdentityTransformation) = α

@inline composeH_id(α::Transformation, β::Transformation) = do_composeH(α, β)
@inline composeH_id(α::Transformation, ::IdIdTransformation) = α
@inline composeH_id(::IdIdTransformation, β::Transformation) = β
@inline composeH_id(α::IdIdTransformation, ::IdIdTransformation) = α

@inline composeH_id(α::Transformation, H::Functor) = do_composeH(α, H)
@inline composeH_id(α::Transformation, ::IdentityFunctor) = α
@inline composeH_id(α::IdentityTransformation, H::Functor) =
  id(compose_id(dom(α), H))
@inline composeH_id(α::IdentityTransformation, ::IdentityFunctor) = α

@inline composeH_id(F::Functor, β::Transformation) = do_composeH(F, β)
@inline composeH_id(::IdentityFunctor, β::Transformation) = β
@inline composeH_id(F::Functor, β::IdentityTransformation) =
  id(compose_id(F, dom(β)))
@inline composeH_id(::IdentityFunctor, β::IdentityTransformation) = β

do_compose(F::Functor, G::Functor) = CompositeFunctor(F, G)

@inline function do_composeH(α::Transformation, β::Transformation)
  do_composeH(α, β, Val{:covariant})
end
function do_composeH(α::Transformation, β::Transformation,
                     ::Type{Val{:covariant}})
  G, H = codom(α), dom(β)
  compose_id(composeH_id(α, H), composeH_id(G, β))
end
function do_composeH(α::Transformation, β::Transformation,
                     ::Type{Val{:contravariant}})
  F, K = dom(α), codom(β)
  compose_id(composeH_id(F, β), composeH_id(α, K))
end

# Oppositization 2-functor
#-------------------------

""" Opposite category, where morphism are reversed.

Call `op(::Cat)` instead of directly instantiating this type.
"""
@auto_hash_equals struct OppositeCat{Ob,Hom,Size<:CatSize,C<:Cat{Ob,Hom,Size}} <:
    Cat{Ob,Hom,Size}
  cat::C
end

ob(C::OppositeCat, x) = ob(C.cat, x)
hom(C::OppositeCat, f) = hom(C.cat, f)

dom(C::OppositeCat, f) = codom(C.cat, f)
codom(C::OppositeCat, f) = dom(C.cat, f)
id(C::OppositeCat, x) = id(C.cat, x)
compose(C::OppositeCat, f, g) = compose(C.cat, g, f)

""" Opposite functor, given by the same mapping between opposite categories.

Call `op(::Functor)` instead of directly instantiating this type.
"""
@auto_hash_equals struct OppositeFunctor{C,D,F<:Functor{C,D}} <: Functor{C,D}
    # XXX: Requires more type parameters: ObC, HomC, ObD, HomD.
    #Functor{OppositeCat{C},OppositeCat{D}}
  func::F
end

dom(F::OppositeFunctor) = op(dom(F.func))
codom(F::OppositeFunctor) = op(codom(F.func))

do_ob_map(F::OppositeFunctor, x) = ob_map(F.func, x)
do_hom_map(F::OppositeFunctor, f) = hom_map(F.func, f)

do_compose(F::OppositeFunctor, G::OppositeFunctor) =
  OppositeFunctor(do_compose(F.func, G.func))

#= Not yet needed because the only natural transformations we currently support
#are `FinTransformationMap`, for which can just implement `op` directly.

""" Opposite natural transformation between opposite functors.

Call `op(::Transformation)` instead of directly instantiating this type.
"""
@auto_hash_equals struct OppositeTransformation{C,D,F,G,T<:Transformation{C,D,F,G}} <: Transformation{C,D,F,G}
    # XXX: Requires more type parameters: ObC, HomC, ObD, HomD.
    #Transformation{OppositeCat{C},OppositeCat{D},OppositeFunctor{C,D,G},OppositeFunctor{C,D,F}}
  trans::T
end

dom(α::OppositeTransformation) = op(codom(α.trans))
codom(α::OppositeTransformation) = op(dom(α.trans))

component(α::OppositeTransformation, x) = component(α.trans, x)

do_compose(α::OppositeTransformation, β::OppositeTransformation) =
  OppositeTransformation(do_compose(β.trans, α.trans))
do_composeH(α::OppositeTransformation, β::OppositeTransformation) =
  OppositeTransformation(do_composeH(α.trans, β.trans))
do_composeH(F::OppositeFunctor, β::OppositeTransformation) =
  OppositeTransformation(do_composeH(F.func, β.trans))
do_composeH(α::OppositeTransformation, H::OppositeFunctor) =
  OppositeTransformation(do_composeH(α.trans, H.func))
=#

""" Oppositization 2-functor.

The oppositization endo-2-functor on Cat, sending a category to its opposite, is
covariant on objects and morphisms and contravariant on 2-morphisms, i.e., is a
2-functor ``op: Catᶜᵒ → Cat``. For more explanation, see the
[nLab](https://ncatlab.org/nlab/show/opposite+category).
"""
op(C::Cat) = OppositeCat(C)
op(F::Functor) = OppositeFunctor(F)
#op(α::Transformation) = OppositeTransformation(α)
op(C::OppositeCat) = C.cat
op(F::OppositeFunctor) = F.func
#op(α::OppositeTransformation) = α.trans

""" 2-cell dual of a 2-category.
"""
function co end

end

max_forks_repo_name: epatters/CompCat


text:
using DrWatson
@quickactivate "INN_Velocity-Migration"

using JLD, PyPlot, LinearAlgebra, Distributed, SharedArrays
addprocs(20)
@everywhere include("../src/RTM_64x200.jl")
@everywhere using .RTM_64x200, JUDI.TimeModeling, Dates

m_all = load(datadir("overthrust_4k_models_200x64.jld"), "m_all1")
n2, n1, num = size(m_all)

# build a migration velocity m0 whose value increases linearly with depth
idx_wb = 9 # index of the water bottom
v0 = ones(1, n1)
v0[idx_wb+1:end] = collect(range(3.5, 4.3, length = n1-idx_wb))
v0 = (1f0 ./ v0).^2
v0[1:idx_wb] .= m_all[1,1,1]
m0 = repeat(v0, n2, 1)

# Right-hand preconditioners (model topmute)
# assign a mask where depth 0 to idx_wb-2 should be exactly 0,
# idx_wb-2 to idx_wb go from 0 to 1 smoothly,
# after idx_wb is exactly 1
n = size(m0)
d = (25., 25.)
o = (0., 0.)
model0 = Model(n, d, o, m0)
# dm = m - m0
# idx_wb = find_water_bottom(dm) # find the index of the water bottom (you could also manually add this)
Tm = judiTopmute(model0.n, idx_wb, 2)  # Mute water column
S = judiDepthScaling(model0) # design a depth scaling operator
Mr = S * Tm

## Computing RTM images
rtm1_all = SharedArray{Float32,3}(n1, n2, num)
rtm2_all = SharedArray{Float32,3}(n1, n2, num)

@sync @distributed for i = 1:num

    println(string("for the ", i, "th model at ", now()))

    m = m_all[:,:,i]
    # m0 = m0_all[:,:,i] # Migration velocity
    # m0 = ones(Float32, n) .* m[1,1] # constant velocity

    # Computing RTM images
    rtm = rtm_isic(adjoint(m), adjoint(m0) )
    rtm1_all[:,:,i] = adjoint(reshape(rtm, n))

    rtm = Mr' * rtm # mute water layers
    rtm2_all[:,:,i] = adjoint(reshape(rtm, n))

end

rtm1_all = Array(rtm1_all)
rtm2_all = Array(rtm2_all)

println("Parallel computation finished")
save(datadir("overthrust_4k_rtm_64x200_lin_vel.jld"), "rtm1_all", rtm1_all, "rtm2_all", rtm2_all)

# figures of rtm examples
figure(figsize=[20,8])
ax1 = subplot(2, 3, 1); imshow(adjoint(m_all[:,:,1]), cmap="jet"); title("True velocity model1")
ax2 = subplot(2, 3, 2); imshow(adjoint(m_all[:,:,2]), cmap="jet"); title("True velocity model2")
ax3 = subplot(2, 3, 3); imshow(adjoint(m_all[:,:,3]), cmap="jet"); title("True velocity model3")
ax4 = subplot(2, 3, 4); plot_rtm(adjoint(rtm1_all[:,:,1]), d; new_fig=false); title("rtm1")
ax5 = subplot(2, 3, 5); plot_rtm(adjoint(rtm1_all[:,:,2]), d; new_fig=false); title("rtm2")
ax6 = subplot(2, 3, 6); plot_rtm(adjoint(rtm1_all[:,:,3]), d; new_fig=false); title("rtm3")
savefig(datadir("figs", "examples_rtm_64x200_1.png"))

figure(figsize=[20,8])
ax1 = subplot(2, 3, 1); imshow(adjoint(m_all[:,:,4]), cmap="jet"); title("True velocity model4")
ax2 = subplot(2, 3, 2); imshow(adjoint(m_all[:,:,5]), cmap="jet"); title("True velocity model5")
ax3 = subplot(2, 3, 3); imshow(adjoint(m_all[:,:,6]), cmap="jet"); title("True velocity model6")
ax4 = subplot(2, 3, 4); plot_rtm(adjoint(rtm1_all[:,:,4]), d; new_fig=false); title("rtm4")
ax5 = subplot(2, 3, 5); plot_rtm(adjoint(rtm1_all[:,:,5]), d; new_fig=false); title("rtm5")
ax6 = subplot(2, 3, 6); plot_rtm(adjoint(rtm1_all[:,:,6]), d; new_fig=false); title("rtm6")
savefig(datadir("figs", "examples_rtm_64x200_2.png"))

####################################################################################################
# generate INN dataset of size n1 x n2 x nc x num
m_all = permutedims(m_all, [2,1,3])
m_all = reshape(m_all, (n1, n2, 1, num))
rtm1_all = reshape(rtm1_all, (n1, n2, 1, num))
rtm2_all = reshape(rtm2_all, (n1, n2, 1, num))

save(datadir("vel_4k_samples_64x200_lin_vel.jld"), "m_all", m_all)
save(datadir("rtm_4k_samples_64x200_lin_vel.jld"), "rtm1_all", rtm1_all, "rtm2_all", rtm2_all)


text:
using Test, GpABC

@testset "RBF Kernel Tests" begin
    @test 2 == get_hyperparameters_size(SquaredExponentialIsoKernel(), zeros(0, 0))
    @test 5 == get_hyperparameters_size(SquaredExponentialArdKernel(), zeros(0, 4))

    # lifted from Matlab
    D2 = [
        1.67814152778737 1.20785927753982 0.624461685654236 1.46408315059547 8.48310688500901 1.45871370121154 6.42543323212624;
        1.14246534460365 1.85035457131647 1.19053740612926 0.611497365423743 6.99118646721199 2.43597535692027 5.07544002153986;
        7.23733046161335 3.16542101637104 2.72647762438244 8.92586928526950 17.1642767758545 1.17728913861911 14.7362524561371;
        19.0995257794008 4.56588070446465 5.99211031953106 11.1739924041597 34.8720573587766 15.7837441497418 30.3384444420246;
        0.466988806886986 3.56126964624601 2.47101254290117 2.18123029882355 4.99339934819349 1.06659412805776 3.56150053043329
    ]
    x1 = [
        0.459530675532078 0.485182659460782;
        0.332063667630806 0.464957557369470;
        0.994061870963342 0.541666708673195;
        0.197865280880530 0.792372670393611;
        0.521775759068213 0.407471305414788
    ]
    x2 = [
        0.370230106343622 0.358754987995492;
        0.457397317261154 0.595083240231072;
        0.476711208222765 0.564018563630068;
        0.0852766000131502 0.459998930552974;
        0.243816421015103 0.202025803831418;
        0.822588000829605 0.447684994641545;
        0.236529943177723 0.241704940442790
    ]
    log_theta = [-log(10)/2, -log(100)/2]
    D2D = [
        0 0.203383856493429 3.17628077043479 10.1212580692727 0.642649959010387;
        0.203383856493429 0 4.97084560152883 10.9001576937249 0.690373692755269;
        3.17628077043479 4.97084560152883 0 12.6246380248108 4.03138234045883;
        10.1212580692727 10.9001576937249 12.6246380248108 0 15.8640860550538;
        0.642649959010387 0.690373692755269 4.03138234045883 15.8640860550539 0
    ]

    @test 4 * ones(5, 1) ≈ covariance_diagonal(SquaredExponentialArdKernel(), [log(2), 0.0, 0.0], x1)
    @test 9 * ones(5, 1) ≈ covariance_diagonal(SquaredExponentialIsoKernel(), [log(3), 0], x1)

    log_theta1 = [log(2); log_theta]
    @test 4 * exp.(-D2D/2) ≈ covariance_training(SquaredExponentialArdKernel(), log_theta1, x1)
    @test 4 * exp.(-D2/2) ≈ covariance(SquaredExponentialArdKernel(), log_theta1, x1, x2)
    @test 4 * exp.(-D2D/2) ≈ covariance_training(SquaredExponentialArdKernel(), log_theta1, x1)
    @test 4 * exp.(-D2/2) ≈ covariance(SquaredExponentialArdKernel(), log_theta1, x1, x2)
    @test 4 * exp.(-D2D/2) ≈ covariance_training(SquaredExponentialArdKernel(), log_theta1, x1)
    @test 4 * exp.(-D2/2) ≈ covariance(SquaredExponentialArdKernel(), log_theta1, x1, x2)

    x1 = [
        0.853623959633347	0.0784714031398629;
        0.857889557656781	0.868872571675513;
        0.549795537511723	0.141312150350452;
        0.413598959526792	0.0333109400390233;
        0.357712766787835	0.856969235863988
    ]
    K = [
        403.428793492735	401.127219311054	403.073312456853	402.706563960842	400.293352208840;
        401.127219311054	403.428793492735	401.129005018406	400.133662383833	402.505043923819;
        403.073312456853	401.129005018406	403.428793492735	403.317183316184	401.405364130706;
        402.706563960842	400.133662383833	403.317183316184	403.428793492735	400.918680124180;
        400.293352208840	402.505043923819	401.405364130706	400.918680124180	403.428793492735
    ]
    @test K ≈ covariance_training(SquaredExponentialIsoKernel(), [3.0, 2.0], x1)
    @test K ≈ covariance(SquaredExponentialIsoKernel(), [3.0, 2.0], x1, x1)
    K = [
        401.977811520080	401.541155407376	402.523543965735	400.671895243459	402.732109940195	402.482565433579	401.263637928358;
        403.301972026813	399.081819073422	403.069411359787	400.973386644735	402.486964554685	400.132636562887	403.418742502755;
        402.078736824856	402.775457037114	402.568058976288	402.175937911770	403.104367431456	403.259534426607	401.152026633490;
        401.293095434847	403.147715014431	401.890001738463	402.191713363450	402.707451571754	403.407040991986	400.131332629428;
        402.684468987700	400.872263943178	402.472762274165	402.816140740862	402.544213669122	401.157890532104	402.329417475725
    ]
    x2 = [
        0.778909885221585	0.701258592634117;
        0.138440919173155	0.0531664015972830;
        0.776528543222942	0.567711645535990;
        0.115225947981340	0.529626523391170;
        0.602507722378945	0.432978635389814;
        0.347290022926074	0.0719249965957948;
        0.903415318533506	0.843417785523188;
    ]
    @test K ≈ covariance(SquaredExponentialIsoKernel(), [3.0, 2.0], x1, x2)

    K = [
        2980.95798704173	2978.65028672416	2978.42445044689	2975.66944742414	2972.01871140493;
        2978.65028672416	2980.95798704173	2976.41449546484	2973.00056583059	2974.13568651361;
        2978.42445044689	2976.41449546484	2980.95798704173	2980.40855933366	2978.05996749357;
        2975.66944742414	2973.00056583059	2980.40855933366	2980.95798704173	2978.36743523035;
        2972.01871140493	2974.13568651361	2978.05996749357	2978.36743523035	2980.95798704173
    ]
    @test K ≈ covariance_training(SquaredExponentialArdKernel(), [4.0, 2.0, 3.0], x1)
    K = [
        2979.37304583894	2967.02516765859	2979.91160562529	2965.36261803110	2978.77301466686	2973.96726404099	2978.72931416283;
        2980.68391806226	2964.41563041355	2980.44223635749	2965.51614192689	2978.47660458701	2971.50931134187	2980.89901384615;
        2978.36770986479	2976.31355555597	2978.88359595894	2975.25091626029	2980.56786900222	2979.82092270343	2975.72768863743;
        2975.67123379203	2978.89038176959	2976.31073840956	2977.61944811412	2979.39404570740	2980.83245035378	2971.99726674807;
        2976.02944167625	2977.26070750597	2975.86478266464	2978.95759612153	2978.65883077957	2978.67897590742	2972.83895895038
    ]
    @test K ≈ covariance(SquaredExponentialArdKernel(), [4.0, 2.0, 3.0], x1, x2)

    R = [
        0.0819471660451718	0.460473784459975	0.831553137317092	0.174000919877552	0.260121663222146;
        0.286019721201630	0.182895418807350	0.913578888361146	0.433460382015193	0.864793774608892;
        0.133454235798425	0.923303895128567	0.520376685783442	0.168640529349745	0.588969851652612;
        0.965341064131512	0.968565425820912	0.970957962316907	0.341376781156928	0.816296725448359;
        0.487264858320013	0.679111772067431	0.890472660255839	0.519560904484038	0.499064610572096;
    ]
    @test [5610.54080009540; 21.9062443257787] * 2 ≈
        covariance_grad(SquaredExponentialIsoKernel(), [3.0, 2.0], x1, R)
    @test [41563.6677709956; 38.4313787854587; 16.7982248189535] * 2 ≈
        covariance_grad(SquaredExponentialArdKernel(), [4.0, 2.0, 3.0], x1, R)

    # R2 = [
    #     0.322721288224087	0.758281509279405	0.220444597406520	0.463054206437922	0.659690255606403	0.473267965239132	0.313020285830377;
    #     0.976691182226265	0.437326458952076	0.0123486113792046	0.552047602158062	0.414125976811222	0.214691774088280	0.931943951142926;
    #     0.248699676826922	0.145795607541593	0.543111450291756	0.264049470547610	0.215292266204624	0.324647101724822	0.181986859153019;
    #     0.895468430079272	0.386775951546978	0.0741450535478746	0.879451493676496	0.417118101109340	0.343838689077570	0.653141784785847;
    #     0.137904342328154	0.543417650062320	0.155914979128127	0.941254220875564	0.503925344675745	0.896847326726116	0.353976895731863
    # ]
    # @test [6374.93996516720; 21.9355980134406] * 2 ≈
    #     covariance_grad(SquaredExponentialIsoKernel(), [3.0, 2.0], x1, x2, mode_regression, R2)
    # @test [4.719203177727137e+04;61.483491155964536;13.680112702542862] * 2 ≈
    #     covariance_grad(SquaredExponentialArdKernel(), [4.0, 2.0, 3.0], x1, x2, mode_regression, R2)
end;


text:
#= = # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Dst Index
#
#   World Data Center for Geomagnetism
#
#   For more information, see:
#
#       http://dx.doi.org/10.17593/14515-74000
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # = =#

################################################################################
#                       Private Structures and Variables
################################################################################

"""
_Dst_Structure

Structure to store the interpolations of the data in Dst files.

# Fields

* `Dst`: Dst index.

"""
struct _Dst_Structure
    Dst::_space_indices_itp_constant{SVector{24,Float64}}
end

# Remote files: *.html (prefered) or for.request (some months are missed)
#
# This set will be configured in the function `init_space_indices()`.
_dstfiles = RemoteFileSet(".dst files",
                    Dict{Symbol,RemoteFile}(),
                )

# Optional variable that will store the WDC data.
@OptionalData(_dst_data, _Dst_Structure,
          "Run `init_space_indices()` with `:dstfiles` in `enabled_files` array to initialize required data.")

################################################################################
#                               Public Functions
################################################################################

export get_Dst

#                                   Getters
# ==============================================================================

"""
get_Dst(DT::DateTime)

Return the Kp index at DateTime `DT`.

"""
function get_Dst(DT::DateTime)
    # @show DT
    Dst_day = get(_dst_data).Dst(datetime2julian(DT))
    # @show Kp_day
    # Get the hour of the day and return the appropriate Kp.
    # y, m, d, h, min, sec = JDtoDate(DT)
    # @show  floor(Int, h / 3)
    # @show Kp_day[ floor(Int, hour(DT) / 3) + 1 ]
    return Dst_day[  hour(DT) + 1 ]
end

################################################################################
#                              Private Functions
################################################################################

"""
_init_dstfiles(;force_download = false, local_dir = nothing, dstfiles_oldest_year = year(now())-3)

Initialize the data in the WDC files by creating `_dstfiles_data`. The
initialization process is composed of:

1. Download the files, if it is necessary;
2. Parse the files;
3. Create the interpolations and the structures.

If the keyword `force_download` is `true`, then the files will always be
downloaded.

The user can also specify a location for the directory with the WDC files using
the keyword `local_dir`. If it is `nothing`, which is the default, then the
file will be downloaded.

The user can select what is the oldest year in which the data will be downloaded
by the keyword `dstfiles_oldest_year`. By default, it will download the data
from 3 previous years.

The user can select what is the newest year in which the data will be downloaded
by the keyword `dstfiles_newest_year`. It it is `nothing`, which is the default,
then it is set to the current year.

"""
function _init_dstfiles(;force_download=false, local_dir=nothing,
                    dstfiles_oldest_year=year(now()) - 3,
                    dstfiles_newest_year=nothing)

    yms     = Int[]
    filepaths = String[]

    if local_dir == nothing
        local_dir = joinpath(path_si_root, "dst")
        isdir(local_dir) || mkpath(local_dir)
        (dstfiles_newest_year == nothing) && (dstfiles_newest_year = year(now()))
        _prepare_dst_remote_files(dstfiles_oldest_year, dstfiles_newest_year, local_dir)
        download(_dstfiles; force=force_download,)

        # Get the files available and sort them by the year.
        for (sym, dstfile) in _dstfiles.files
            #
            # The year must not be obtained by the data inside the file,
            # because it contains only 2 digits and will break in 2032.
            # We will obtain the year by the symbol of the remote file. The
            # symbol name is:
            #
            #       kpYYYY
            #
            # where `YYYY` is the year.
            push!(yms, parse(Int, String(sym)[4:9]))
            push!(filepaths, path(dstfile))
        end
    else
    # If the user provided a directory, check what files are available.
    # Notice that the name must be the same as the ones online.
        for (root, dirs, files) in walkdir(local_dir)
            for file in files
                if occursin(r"^dst[1-2][0-9][0-9][0-9].for.request$", file) ||
                    occursin(r"^dst[1-2][0-9][0-9][0-9].html$", file)
                    year = parse(Int, file[4:7])
                    # Check if the year is not older than the oldest year.
                    if year >= dstfiles_oldest_year
                        @info "Found Dst file `$file` related to the year `$year`."
                        push!(filepaths, joinpath(root, file))
                        push!(yms,     file[4:9])
                    end
                end
            end
        end
    end

    p = sortperm(yms)

    push!(_dst_data,       _parse_dstfiles(filepaths[p]))

    nothing
end

"""
_parse_dstfiles(filepaths::Vector{String}, years::Vector{Int})

Parse the Dst files with paths in `filepaths` related to the years in `years`.

**Notice that the files must be sorted by the year!**

"""
function _parse_dstfile_html_line(ln::SubString)
    Dst_k = zeros(Float64, 24)
    for i in 1:3
        for j in 1:8
            k = i * 8 - 8 + j
            Dst_k[k] = parse(Int, ln[ i - 1 + k * 4:i - 1 + 3 + k * 4])
        end
    end
    Dst_k
end

function _parse_dstfile_html(filepath::String,
                            DT::Vector{DateTime},
                            Dst::Vector{SVector{24,Float64}})

    bn = basename(filepath)
    y = parse(Int, bn[4:7])
    m = parse(Int, bn[8:9])

    open(filepath, "r") do file
        text = read(file, String)
        a = r"""
\n((.|\n)+?)\n
"""
        b = match(a, text)
        if b == nothing
            @warn "no match " * filepath
            return
        end

        alllines = split(b.captures[1], "\n")
        # utline = alllines[6]
        lines = alllines[8:end - 1]
        filter!(!isempty, lines)

        for ln in lines
            d   = parse(Int, ln[1:2])
            DT_k  = DateTime(y, m, d, 12, 0, 0)

            # Get the vector of Dst
            Dst_k = _parse_dstfile_html_line(ln)
            # Add data to the vector.
            push!(DT, DT_k)
            push!(Dst, Dst_k)
        end
    end

end

function _parse_dstfile_forrequest(filepath::String,
                                DT::Vector{DateTime},
                                Dst::Vector{SVector{24,Float64}})

    bn = basename(filepath)
    y  = parse(Int, bn[4:7])
    m  = parse(Int, bn[8:9])

    open(filepath, "r") do file

        for ln in eachline(file)
            # Get the Julian Day.
            d   = parse(Int, ln[9:10])

            # The DT of the data will be computed at noon. Hence, we will be
            # able to use the nearest-neighbor algorithm in the
            # interpolations.
            DT_k  = DateTime(y, m, d, 12, 0, 0)

            # Get the vector of Kps and Aps.
            Dst_k = zeros(Float64, 24)

            for i = 1:24
                Dst_k[i] = parse(Int, ln[17 + i * 4:20 + i * 4])
            end

            # Add data to the vector.
            push!(DT, DT_k)
            push!(Dst, Dst_k)
        end
    end
end

function _parse_dstfiles(filepaths::Vector{String}, )
    # Allocate the raw data.
    DT = DateTime[]
    Dst = SVector{24,Float64}[]

    for filepath in filepaths
        if endswith(filepath, "html")

            _parse_dstfile_html(filepath, DT, Dst)
        elseif endswith(filepath, "for.request")
            _parse_dstfile_forrequest(filepath, DT, Dst)
        else
            @warn "can't recognise filetype " * filepath
        end
    end

    # Create the interpolations for each parameter.
    knots    = (datetime2julian.(DT),)
    # Create the interpolations.
    itp_Dst = interpolate(knots, Dst, Gridded(Constant()))
    _Dst_Structure(itp_Dst)
end

"""
_prepare_dst_remote_files(oldest_year::Number, newest_year::Number)

Configure all the WDC remote files between `newest_year` and `oldest_year`.
Notice that previous years will never be updated whereas the current year will
be updated daily.

If `oldest_year` is greater than current year, then only the files from the
current year will be downloaded.

If `newest_year` is smaller than `oldest_year`, then only the files from the
`oldest_year` will be downloaded.

This function modifies the global variable `_wdcfiles`.

"""

function _prepare_dst_remote_files(oldest_year::Number, newest_year::Number, local_dir::String)
    # Get the current year.
    current_year = year(now())

    # If `oldest_year` is greater than current year, then consider only the
    # current year.
    (oldest_year > current_year) && (oldest_year = current_year)
    (newest_year < oldest_year)  && (newest_year = oldest_year)
    (newest_year > current_year) && (newest_year = current_year)

    # For the current year, we must update the remote file every day. Otherwise,
    # we do not need to update at all.

    # http://wdc.kugi.kyoto-u.ac.jp/dst_final/199412/index.html
    baseurl = "http://wdc.kugi.kyoto-u.ac.jp/"
    # end_time => dst type
    url_prefix = [
            DateTime(2014, 12, 31, 23, 59, 59) => "dst_final",
            DateTime(2017, 12, 31, 23, 59, 59) => "dst_provisional",
            now() => "dst_realtime",
        ]

    yms = Date(oldest_year):Month(1):Date(newest_year + 1) - Day(1)
    start_index = 1
    for ym in yms
        ymstr = @sprintf("%04d%02d",year(ym),month(ym))
        forrequest = string("dst",
                        @sprintf("%02d%02d",
                                    year(ym) % 100,
                                    month(ym)),
                        ".for.request",
                )
        httpindex = "index.html"
        # use html (for.request files are missed in some months)
        fileremote = httpindex

        filename = string("dst", ymstr, ".html")
        sym = Symbol(filename)
        for (idt, isp) in url_prefix[start_index:end]
            if ym > idt
                start_index += 1
                continue
            end
            file_y = @RemoteFile(
                    string(baseurl, isp, "/", ymstr, "/", fileremote),
                    file = filename, dir = joinpath(local_dir, isp),
                    updates = (isp == "dst_final") ? :never : :yearly,
                )
            merge!(_dstfiles.files, Dict(sym => file_y))
            break
        end
    end
    nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1RollingUpdateDaemonSet } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy

const _property_map_IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiExtensionsV1beta1RollingUpdateDaemonSet", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1DaemonSetUpdateStrategy }, name::Symbol, val)
end


text:
import Base: .^

cuda10 = [
# ("add",".+","s+xi"),
# ("sub",".-","s-xi"),
# ("mul",".*","s*xi"),
# ("div","./","s/xi"),
("pow",".^","pow(xi,s)"),
# "hypot",
# "rhypot",
# "atan2",
# "frexp",
# "ldexp",
# "scalbn",
# "scalbln",
# "jn",
# "yn",
# "fmod",
# "remainder",
# "mod",
# "fdim",
]

(.^){T}(a::KnetArray{T},s::Number)=(.^)(a,T(s))

function cuda10def(f, j=f, o...)
    J=Symbol(j)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$(S)_10"
        @eval begin
            function $J(x::KnetArray{$T},s::$T)
                y = similar(x)
                ccall(($F,$libknet8),Void,(Cint,Ptr{$T},$T,Ptr{$T}),length(y),x,s,y)
                return y
            end
        end
    end
end

#if isdefined(:libknet8)
    for f in cuda10
        isa(f,Tuple) || (f=(f,))
        cuda10def(f...)
    end
#end


text:
using ReservoirComputing, GaussianProcesses, MLJLinearModels, Random, Statistics, LIBSVM

const res_size=20
const ts = 0.:0.1:50.0
const data = sin.(ts)
const train_len = 400
const predict_len = 100
const input_data = reduce(hcat, data[1:train_len-1])
const target_data = reduce(hcat, data[2:train_len])
const test = reduce(hcat, data[train_len+1:train_len+predict_len])
const reg = 10e-6

Random.seed!(77)
esn = ESN(input_data;
    reservoir=RandSparseReservoir(res_size, 1.2, 0.1))

training_methods = [StandardRidge(regularization_coeff=reg), LinearModel(RidgeRegression, regression_kwargs=(;lambda=reg)), LinearModel(regression=RidgeRegression, regression_kwargs=(;lambda=reg)), GaussianProcess(MeanZero(), Poly(1.0, 1.0, 2)), EpsilonSVR()]

for t in training_methods
    output_layer = train(esn, target_data, t)
    output = esn(Predictive(input_data), output_layer)
    if t isa GaussianProcess
        output = output[1]
    else
        output = output
    end
    @test mean(abs.(target_data .- output)) ./ mean(abs.(target_data)) < 0.21
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
module KernelDensityDistributionEsimates

using   KernelDensity,
        Interpolations,
        Distributions,
        StatsBase, Statistics,
        Gadfly

export KDE, kdemean

struct KDE{T,ITP <: ScaledInterpolation} <: ContinuousUnivariateDistribution
    x::StepRangeLen{T,Base.TwicePrecision{T},Base.TwicePrecision{T}}
    density::Vector{T}
    cumulative_density::Vector{T}
    pdf::ITP
    cdf::ITP
end
function KDE(kde::UnivariateKDE)
    x = kde.x
    # if length(x) == 1
    #     @show x
    #     return KDE(x, kde.density, [1.0], [1.0], [1.0])
    # end
    density = kde.density
    cumulative_density = cumsum(density)
    cumulative_density ./= cumulative_density[end]

    pdf = Interpolations.scale(interpolate(density, BSpline(Quadratic(Line(OnGrid())))), x)
    cdf = Interpolations.scale(interpolate(cumulative_density, BSpline(Quadratic(Line(OnGrid())))), x)

    KDE(x, density, cumulative_density, pdf, cdf)
end
KDE(distances::AbstractVector) = KDE(kde(distances))
KDE(distances::AbstractVector, weights::AbstractVector) = KDE(kde(distances, weights = weights))

kdemean(kde::KDE) = mean(kde)
function kdemean(kde::UnivariateKDE)
    if length(kde.x) == 1
        return kde.x[1]
    else
        return mean(KDE(kde))
    end
end

# Parameters are summarized by the full (x, density) set
StatsBase.params(kde::KDE) = (kde.x, kde.density)
Statistics.mean(kde::KDE{T}) where T = kde.x' * kde.density * T(kde.x.step)
function Statistics.var(kde::KDE{T}) where T
    μ = zero(T)
    σ² = zero(T)
    x = kde.x; density = kde.density
    @inbounds @fastmath for i ∈ eachindex(x)
        xd = x[i] * density[i]
        μ += xd
        σ² += x[i] * xd
    end
    σ² * T(kde.x.step) - (μ * T(kde.x.step))^2
    # σ², μ, T(kde.x.step)
end
Statistics.std(kde::KDE) = sqrt(var(kde))
function Statistics.median(kde::KDE)
    kde.x[findfirst(p -> p > 0.5, kde.cumulative_density)]
end
function StatsBase.mode(kde::KDE)
    kde.x[argmax(kde.density)]
end
function StatsBase.entropy(kde::KDE{T}) where T
    out = zero(T)
    density = kde.density
    @inbounds for i ∈ eachindex(density)
        out -= density[i] * log(density[i])
    end
    out * T(kde.x.step)
end
function Distributions.pdf(kde::KDE, x::T) where T <: Real
    if (x < minimum(kde.x)) || (x > maximum(kde.x))
        zero(T)
    else
        max(zero(T), kde.pdf(x))
    end
end
Distributions.logpdf(kde::KDE, x::Real) = log(pdf(kde, x))
function Distributions.cdf(kde::KDE{T}, x::Real) where T
    isnan(x) && return T(NaN)
    if x < minimum(kde.x)
        p = zero(T)
    elseif x > maximum(kde.x)
        p = one(T)
    else
        # the min and max should ideally be unnecessary, but
        # special cases have come up in practice.
        # until we have a better way
        p = min(max(zero(T), kde.cdf(x)),one(T))
    end
    p
end

function Gadfly.layer(kde::KDE, args::Vararg{Union{Function, Gadfly.Element, Theme, Type},N} where N)
    layer(x = kde.x, y = kde.density, Geom.line, args...)
end
Gadfly.plot(kde::KDE, args::Vararg{Union{Function, Gadfly.Element, Theme, Type},N} where N) = plot(layer(kde, args...))

end # module


text:
using Documenter
using PWF

makedocs(
    modules = [PWF],
    format = Documenter.HTML(analytics = "UA-367975-10", mathengine = Documenter.MathJax()),
    sitename = "PWF",
    authors = "Iago Chávarry and Pedro Hamacher",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting Started" => "quickguide.md",
        ],
        "PWF File" => [
            "Overview" => "pwf_overview.md",
            "DBAR" => "dbar.md",
            "DLIN" => "dlin.md",
            "DGER" => "dger.md",
            "DGBT" => "dgbt.md",
            "DGLT" => "dglt.md",
            "DSHL" => "dshl.md",
            "DELO" => "delo.md",
            "DCBA" => "dcba.md",
            "DCLI" => "dcli.md",
            "DCNV" => "dcnv.md",
            "DCCV" => "dccv.md",
            "DBSH" => "dbsh.md",
            "DCER" => "dcer.md",
            "DCSC" => "dcsc.md",
            "DOPC" => "dopc.md",
            "DCTE" => "dcte.md",
        ],
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
        repo="github.com/LAMPSPUC/PWF.jl.git",
        devbranch = "main",
        push_preview = true
    )