text:
# This file is part of IntegerTriangles.jl.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module TrianglesUtils

using Nemo, TrianglesBase, HTTP

export Show, GetSeqnum, GetSeqnumUri, MakeSeqUri, SeqToString, SeqToFixlenString
export profilespath, datapath, docsrcpath, csv_files
export oeis_notinstalled, oeis_search, search_failed
export devdoctrait

"""
* Search the OEIS for a sequence.
* Save a file in the 'data' directory in json format.

* Triangle in standard format (without typeinfo):
```
[1]
[1, 1]
[2, 4, 1]
[6, 18, 9, 1]
[24, 96, 72, 16, 1]
[120, 600, 600, 200, 25, 1]
[720, 4320, 5400, 2400, 450, 36, 1]
```
* Triangle as a nested array:
```
[[1], [1, 1], [2, 4, 1], [6, 18, 9, 1], [24, 96, 72, 16, 1], [120, 600, 600, 200, 25, 1], [720, 4320, 5400, 2400, 450, 36, 1]]
```
* Triangle in flattened format:
```
1, 1, 1, 2, 4, 1, 6, 18, 9, 1, 24, 96, 72, 16, 1, 120, 600, 600, 200, 25, 1, 720, 4320, 5400, 2400, 450, 36, 1,
```
* Triangle in info format with sequence ID:
```
[A021009] 1, 1, 1, 2, 4, 1, 6, 18, 9, 1, 24, 96, 72, 16, 1, ...
```
* Triangle in mapped format:
```
0 ↦ 1
1 ↦ 1
2 ↦ 1
3 ↦ 2
4 ↦ 4
5 ↦ 1
6 ↦ 6
7 ↦ 18
8 ↦ 9
9 ↦ 1
```
"""
const ModuleTrianglesUtils = ""

const srcdir = realpath(joinpath(dirname(@__FILE__)))
const ROOTDIR = dirname(srcdir)
const datadir = joinpath(ROOTDIR, "data")
const profilesdir = joinpath(ROOTDIR, "profiles")
const docsdir = joinpath(ROOTDIR, "docs")
const docsrcdir = joinpath(docsdir, "src")

datapath(name) = joinpath(datadir, name)
profilespath(name) = joinpath(profilesdir, name)
docsrcpath(name) = joinpath(docsrcdir, name)

oeis_file() = joinpath(datadir, "stripped")
is_oeis_installed() = isfile(oeis_file())

csv_files() = filter!(s -> occursin(r"\.csv$", s), readdir(profilesdir))

function oeis_notinstalled()
    if !is_oeis_installed()
        @warn("OEIS data not installed! Download stripped.gz from oeis.org,")
        @warn("expand it and put it in the directory ../data.")
        return true
    end
    return false
end

function search_failed(name)
    for ln ∈ eachline(name)
        occursin("greeting", ln) && continue
        occursin("{", ln) && continue
        occursin("start", ln) && continue
        if occursin("\"results\": null", ln)
            return true
        end
    end
    false
end

"""
Search the OEIS for a sequence. The file is saved in the 'data' directory in json format.
"""
function oeis_search(seq)

    seqstr = AbsSeqToString(seq[1:min(end,12)])
    filename = joinpath(datadir, seqstr[1:min(end,12)] * ".json")
    uristr = "https://oeis.org/search?q=" * seqstr * "&go=Search&fmt=json"

    url = HTTP.URI(uristr)
    tries = 3
    r = nothing
    for i = 1:tries
        try
            r = HTTP.get(url; readtimeout = 2)
            getfield(r, :status) == 200 && break
            getfield(r, :status) == 302 && break
        catch e
            @warn(e)
        end
        sleep(2)
    end
    if !(r === nothing) && getfield(r, :status) == 200
        open(filename, "w") do f
            write(f, getfield(r, :body))
        end
        @info("Dowloaded " * basename(filename) * " to " * datadir)
    else
        if r === nothing
            @warn("Search did not succeed, connection timed out.\n")
        else
            @warn("Search did not succeed, $url\nStatus: $(getfield(r, :status))")
        end
    end
    filename
end

function AbsSeqToString(seq::ℤSeq, max=100)
    separator = ","
    str = ""
    c = 1
    for trm in seq
        str *= string(abs(trm)) * separator
        c += 1
        c > max && break
    end
    str
end

function SeqToString(seq::ℤSeq, max::Int=100)
    separator = " "
    str = "["
    c = 1
    for trm in seq
        str *= string(trm) * separator
        c += 1
        c > max && break
    end
    str * "]"
end

function SeqToFixlenString(seq::ℤSeq, maxlen::Int=72)
    #maxlen = 72
    separator = " "
    str = "["
    len = 3
    for trm in seq
        s = string(trm) * separator
        len += length(s)
        len > maxlen && break
        str *= s
    end
    str * "]"
end

# increases accuracy and prevents premature matches
const minlen = 30  # fragil! do not reduce!

function GetSeqnum(seq::ℤSeq)
    str = AbsSeqToString(seq)
    soff = 1; loff = 10
    #println(seq)
    #println(str)
    for ln ∈ eachline(oeis_file())
        ln[1] == '#' && continue
        l = replace(ln, "-" => "")
        m = min(length(str), length(l))
        len = min(length(str), length(l)-loff) - 1
        len < minlen && continue
        cmp(str[soff:soff+len], l[loff:loff+len]) == 0 && return ln[1:7]
    end
    toff = findfirst(",", str)[1] + 1
    for ln ∈ eachline(oeis_file())
        ln[1] == '#' && continue
        l = replace(ln, "-" => "")
        len = min(length(str)-toff, length(l)-loff) - 1
        len < minlen && continue
        cmp(str[toff:toff+len], l[loff:loff+len]) == 0 && return "A" * ln[2:7]
    end
    for ln ∈ eachline(oeis_file())
        ln[1] == '#' && continue
        l = replace(ln, "-" => "")
        ff = findnext(",", l, 10)
        ff === nothing && continue
        loff = ff[1] + 1; soff = 1
        len = min(length(str)-soff, length(l)-loff) - 1
        len < minlen && continue
        cmp(str[soff:soff+len], l[loff:loff+len]) == 0 && return "A" * ln[2:7]
        len = min(length(str)-toff, length(l)-loff) - 1
        len < minlen && continue
        cmp(str[toff:toff+len], l[loff:loff+len]) == 0 && return "A" * ln[2:7]
    end
    nothing
end

GetSeqnum(seq::ℤTri) = GetSeqnum(Flat(seq))

function GetSeqnum(seq, search=false)
    if seq === nothing || seq == []
        return nothing
    end
    anum = GetSeqnum(seq)
    if search
        if anum === nothing
            filename = oeis_search(seq)
            if ! search_failed(filename)
                println()
                Println(seq[1:min(18, end)])
                @warn "Not found in the local base but possibly in the OEIS!"
                println()
            end
        end
    end
    anum
end

function GetSeqnumUri(seq::ℤSeq, len=10)
    anum = GetSeqnum(seq)
    if anum === nothing
        s = SeqToFixlenString(seq, 72)
        return "" * "nomatch"
    end
    uri = joinpath("https://oeis.org/", anum)
    return "" * anum * ""
end

function MakeSeqUri(anum, seq)
    if anum == "nothing"
        return "" * "nomatch"
    end
    uri = joinpath("https://oeis.org/", anum)
    return "" * anum * ""
end

function print_without_type(io, v::AbstractVector)
    print(io, "[")
    for (i, el) ∈ enumerate(v)
        i > 1 && print(io, ", ")
        print(io, el)
    end
    println(io, "]")
end

"""
Print the array ``A`` in the format ``n ↦ A[n]`` for n in the given range.
"""
function MappedShow(A::Array, R::AbstractRange, offset=0)
    for k ∈ R
        if isassigned(A, k)
            println(k + offset, " ↦ ", A[k])
        else
            println(k + offset, " ↦ ", "undef")
        end
    end
end

"""
Print an integer triangle without typeinfo.
"""
function Show(T::ℤTri, format="std")

    if format == "std"  # default
        for t in T[1:min(10, end)]
            print_without_type(IOContext(stdout), t)
        end

    elseif format == "nest"
        T |> println

    elseif format == "flat"
        for n in Flat(T)
            print(n, ", ")
        end

    elseif format == "map"
        offset = -1 # Julia is unfortunately 1-based.
        for n in T
            MappedShow(n, 1:length(n), offset)
            offset += length(n)
        end

    elseif format == "info"
        S = Flat(T)
        print("[", GetSeqnum(S), "] ")
        for t in T[1:min(5, end)], k in t
            print(k, ", ")
        end
        println("... ")
    end
end

# https://openlibmathseq.github.io/IntegerTriangles.jl/dev/#IntegerTriangles.AltSum
const devdoc = "https://openlibmathseq.github.io/IntegerTriangles.jl/dev/#IntegerTriangles."
devdoctrait(trait) =  "" * trait * ""

# START-TEST-########################################################

using TrianglesExamples

function test()

    T = LaguerreTriangle(7)

    println("\nThe $Laguerre triangle in different formats:")

    println("\nTriangle in standard format:\n")
    Show(T)

    println("\nTriangle as a nested array:\n")
    Show(T, "nest")

    println("\n\nTriangle in flattened format:\n")
    Show(T, "flat")

    println("\n\nTriangle in info format with sequence ID:\n")
    Show(T, "info")

    println("\nTriangle in mapped format:\n")
    Show(T, "map")
end

function demo()
    N = 7

    println("\nTriangle in standard format:\n")
    Show(LaguerreTriangle(N))

    println("\nTriangle as a nested array:\n")
    Show(LaguerreTriangle(N), "nest")

    println("\n\nTriangle in flattened format:\n")
    Show(LaguerreTriangle(N), "flat")

    println("\n\nTriangle in mapped format:\n")
    Show(LaguerreTriangle(N), "map")

    println("\nTriangle transform of the squares:\n")
    LaguerreTransform([ZZ(k^2) for k in 0:N - 1]) |> println

    println("\nTriangle transform of the IN:\n")
    LaguerreTransform([ZZ(k) for k in 0:N])   |> println
    LaguerreTransform([ZZ(k) for k in 1:N + 1]) |> println

    println()
    T = LaguerreTriangle(9)
    PosHalf(T) |> println
    NegHalf(T) |> println

    println()
    P = Polynomial(T)
    println.(P)
    E = Evaluate(P, ZZ(3))
    println(E)

    println()
    Println.(PolyTriangle(T))
    println()
    Println.(PolyArray(T))
    println()

    A = ℤInt[0, 1, 2, 7, 44, 361, 3654, 44207, 622552, 10005041,
    180713290, 3624270839, 79914671748, 1921576392793, 50040900884366,
    1403066801155039, 42142044935535536]
    SeqToString(A, 12) |> println
    SeqToFixlenString(A, 12) |> println
end

function perf()
    GetSeqnum(ℤInt[1, 1, -2, 3, -3, 3, -5, 7, -6, 6, -10,
    12, -11, 13, -17, 20, -21, 21, -27, 34, -33, 36, -46, 51,
    -53, 58, -68, 78, -82, 89, -104]) |> println
    GetSeqnum(ℤInt[0, 1, 1, 1, 2, 1, 2, 1, 5, 2, 2, 1, 5, 1,
    2, 1, 14, 1, 5, 1, 5, 2, 2, 1, 15, 2, 2, 5, 4, 1, 4, 1, 51,
    1, 2, 1, 14, 1, 2, 2, 14, 1, 6, 1, 4, 2, 2, 1, 52]) |> println
    GetSeqnum(ℤInt[1, 1, 7, 37, 241, 2101, 18271, 201097, 2270017,
    29668681, 410815351, 6238931821, 101560835377, 1765092183037,
    32838929702671, 644215775792401]) |> println
    GetSeqnum(ℤInt[1, 1, 1, 7, 37, 241, 2101, 18271, 201097, 2270017,
    29668681, 410815351, 6238931821, 101560835377, 1765092183037,
    32838929702671, 644215775792401]) |> println
    GetSeqnum(ℤInt[0, 1, 2, 7, 44, 361, 3654, 44207, 622552, 10005041,
    180713290, 3624270839, 79914671748, 1921576392793, 50040900884366,
    1403066801155039, 42142044935535536]) |> println
    GetSeqnum(ℤInt[0, 70, 3783, 338475, 40565585, 6061961733,
    1083852977811, 225615988054171, 53595807366038234, 14308700593468127485,
    4241390625289880226714]) |> println
    GetSeqnumUri(ℤInt[1, 1, 7, 37, 241, 2101, 18271, 201097, 2270017,
    29668681, 410815351, 6238931821, 101560835377, 1765092183037,
    32838929702671, 644215775792401]) |> println
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
#=
vector regression problem with options for loss, regularization, and prior information
expressed through convex constraints

sets:
i ∈ 1..n   observations
j ∈ 1..m   independent variables
l ∈ 1..L   prior constraints

parameters:
Xᵢⱼ ∈ ℝ   jth element of independent variables for observation i
Yᵢ  ∈ ℝ   dependent variable for observation i
λ   ∈ ℝ₊  regularization parameter
f   loss function (domain ℝ)
hₗ   lth prior function (domain ℝ₊)
aₗ  ∈ ℝ   lth prior constant

variables:
βⱼ ∈ ℝ₊   jth parameter (assumed nonnegative)

objective:
min f(Y - X β) + λ L₀(β)   minimize loss plus L₀ regularization

constraints:
hₗ(β) ≤ aₗ, ∀l   lth prior constraint
=#

struct VectorRegression <: ExampleInstance
    n::Int
    m::Int
    f_l2::Bool # if true, f is L2 norm, else L1 norm
    hs::Vector{<:VecSpecExt}
end

function build(inst::VectorRegression)
    (n, m) = (inst.n, inst.m)
    hs = inst.hs

    # generate parameters
    βmax = 3
    β0 = βmax * rand(m) # random solution used to construct data
    X = randn(n, m)
    Y = X * β0 + 0.1 * randn(n)
    λ = 0.1
    a = [get_val(1 .+ β0, h_l) for h_l in hs]

    K_f = if inst.f_l2
        Hypatia.EpiNormEuclCone{Float64}(1 + n)
    else
        Hypatia.EpiNormInfCone{Float64, Float64}(1 + n, true)
    end

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, β[1:m] >= 0)
    JuMP.@variable(model, z[1:m], Bin)
    JuMP.@constraint(model, β .<= βmax * z)

    JuMP.@variable(model, f_epi)
    f_aff = JuMP.@expression(model, vcat(f_epi, Y - X * β))
    JuMP.@constraint(model, f_aff in K_f)

    JuMP.@objective(model, Min, f_epi + λ * sum(z))

    for (h_l, a_l) in zip(hs, a)
        add_homog_spectral(h_l, m, vcat(a_l, 1 .+ β), model)
    end

    # save for use in tests
    model.ext[:z] = z
    model.ext[:β] = β
    model.ext[:f_aff] = f_aff
    model.ext[:a] = a

    return model
end

function test_extra(inst::VectorRegression, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    z = JuMP.value.(model.ext[:z])
    @test z ≈ round.(Int, z) atol = tol rtol = tol
    @test all(-tol .<= z .<= 1 + tol)
    # check conic feasibility
    f_aff = JuMP.value.(model.ext[:f_aff])
    f_val = norm(f_aff[2:end], (inst.f_l2 ? 2 : 1))
    @test f_aff[1] ≈ f_val atol = tol rtol = tol
    β = JuMP.value.(model.ext[:β])
    a = model.ext[:a]
    @test all(a_l >= get_val(1 .+ β, h_l) - tol for (h_l, a_l) in zip(inst.hs, a))
    return
end

ext: jl


text:
"""
    Set a mapping properties for the vertices.

    Depending on the type of the elements of 'arr' the
    properties are set in different variables.

    If there is already a property with the given 'name',
    the function will replace the old array with the new
    one, without any warning.

    'name':
        name of the property.
    'arr':
        array of typed elements.
    'graph':
        graph to be assigned with the property.
"""
function set_vertices_properties(name::String, arr::Array{Int}, graph::Graph)
    if length(graph.vertices) == length(arr)
        graph.int_vproperties[name] = arr
    else
        print("size array does not match number of vertices\n")
    end
end

function set_vertices_properties(name::String, arr::Array{Float64}, graph::Graph)
    if length(graph.vertices) == length(arr)
        graph.float_vproperties[name] = arr
    else
        print("size array does not match number of vertices\n")
    end
end

function set_vertices_properties(name::String, arr::Array{String}, graph::Graph)
    if length(graph.vertices) == length(arr)
        graph.string_vproperties[name] = arr
    else
        print("size array does not match number of vertices\n")
    end
end

"""
    Set a mapping properties for the edges.

    Depending on the type of the elements of 'arr' the
    properties are set in different variables.

    If there is already a property with the given 'name',
    the function will replace the old array with the new
    one, without any warning.

    'name':
        name of the property.
    'arr':
        array of typed elements.
    'graph':
        graph to be assigned with the property.
"""
function set_edges_properties(name::String, arr::Array{Int}, graph::Graph)
    if length(graph.edges) == length(arr)
        graph.int_eproperties[name] = arr
    else
        print("size array does not match number of edges\n")
    end
end

function set_edges_properties(name::String, arr::Array{Float64}, graph::Graph)
    if length(graph.edges) == length(arr)
        graph.float_eproperties[name] = arr
    else
        print("size array does not match number of edges\n")
    end
end

function set_edges_properties(name::String, arr::Array{String}, graph::Graph)
    if length(graph.edges) == length(arr)
        graph.string_eproperties[name] = arr
    else
        print("size array does not match number of edges\n")
    end
end

"""
    Print all the properties associated with 'graph'.

    graph:
        Graph object.
"""
function list_properties(graph::Graph)
    print("$(keys(graph.int_vproperties)) - Int - Vertex\n")
    print("$(keys(graph.float_vproperties)) - Float - Vertex\n")
    print("$(keys(graph.string_vproperties)) - String - Vertex\n")

    print("$(keys(graph.int_eproperties)) - Int - Edge\n")
    print("$(keys(graph.float_eproperties)) - Float - Edge\n")
    print("$(keys(graph.string_eproperties)) - String - Edge\n")
end


text:
using JLD2
using PyPlot
using Statistics
using Printf
nfile = 8
include("../../../src/loglinspace.jl")
include("../../../src/histogram_code.jl")

# Load in likelihood profile:
#@load "T1_likelihood_profile_student_all_3.0sig.jld2"
@load "../../../data/T1_likelihood_profile_student_all_3.0sig_v02.jld2"

# Load in the Markov chain dataset:
@load "../../../data/T1_hmc_total_02212020.jld2"
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v15/Hyak/T1_hmc_total_02212020.jld2"

## Now, make plots of each parameter versus the others:
#planet =["b","c","d","e","f","g","h"]
#var = ["m","P","t0","ecos","esin","log(nu)","V1 e^{1/(2nu)}"]
#
#for i=1:36
#  if i < 36
#    ip1 = ceil(Int64,i/5); ip2 = mod(i-1,5)+1
#    x = vec(elements_grid_all[:,ip1+1,ip2,:])
#    xname = string(planet[ip1]," ",var[ip2])
#  else
#    x = log.(vec(ndof_grid_all))
#    xname = var[6]
#  end
#  for j=i+1:37
##  for j=36:37
#    if j < 36
#      jp1 = ceil(Int64,j/5); jp2 = mod(j-1,5)+1
#      y = vec(elements_grid_all[:,jp1+1,jp2,:])
#      yname = string(planet[jp1]," ",var[jp2])
#    elseif j < 37
#      y = log.(vec(ndof_grid_all))
#      yname = var[6]
#    else
#      y = exp.(vec(lnV1_grid_all) .+ 1 ./(2vec(ndof_grid_all)))
#      yname = var[7]
#    end
#    if(abs(cor(x,y)) > 0.5)
#      clf()
#      cmap = exp.(-0.5 .*(vec(chi_grid_all) .- minimum(chi_grid_all)))
#      scatter(x[cmap .> 0.6],y[cmap .> 0.6],c=cmap[cmap .> 0.6])
#      println(i," ",xname," ",j," ",yname," ",cor(x,y))
#      read(stdin,Char)
#    end
#  end
#end

planet =["b","c","d","e","f","g","h"]
cp = ["C0","C1","C2","C3","C4","C5","C6","C7","C8","C9"]
x2=collect(linearspace(-0.015,0.015,1000))
fig,axes = subplots(4,2,figsize=(10,8),sharex="col",sharey="row")
for i=1:7
  ax = axes[i]
  x= elements_grid_all[(i-1)*5+4,i+1,4,:]; ecos0 = elements_grid_all[(i-1)*5+4,i+1,4,16]
  ecc = sqrt.(x.^2 .+ elements_grid_all[(i-1)*5+5,i+1,5,:].^2)
  prob = exp.(-0.5*(chi_grid_all[(i-1)*5+4,:] .-chi_grid_all[(i-1)*5+4,16]))
  ax.plot(x,prob./maximum(prob),color=cp[i],linewidth=1)
  prob2 = exp.(-0.5 .*(x2 .-ecos0).^2 ./cov_save[(i-1)*5+4,(i-1)*5+4])
  ax.plot(x2,prob2,color=cp[i],alpha=0.3,linewidth=1)
  # Plot histogram:
  ecos_bin,ecos_hist,ecos_bin_square,ecos_hist_square = histogram(state_total[(i-1)*5+4,:],50)
#  ax.plot(ecos_bin_square,ecos_hist_square ./maximum(ecos_hist_square),color=cp[i],linewidth=3,label=L"$e\cos{\omega}$)
  ax.plot(ecos_bin_square,ecos_hist_square ./maximum(ecos_hist_square),color=cp[i],linewidth=3,label=L"$e\cos \omega$")
  x = elements_grid_all[(i-1)*5+5,i+1,5,:]; esin0 = elements_grid_all[i*5,i+1,5,16]
  ecc = sqrt.(x.^2 + elements_grid_all[5i-4,i+1,4,:].^2)
  prob = exp.(-0.5*(chi_grid_all[(i-1)*5+5,:] .-chi_grid_all[(i-1)*5+5,16]))
  ax.plot(x,prob./maximum(prob),linestyle=":",color=cp[i],linewidth=1)
  prob2 = exp.(-0.5 .*(x2 .-esin0).^2 ./cov_save[i*5,i*5])
  ax.plot(x2,prob2,color=cp[i],linestyle=":",alpha=0.3,linewidth=1)
  # Plot histogram:
  esin_bin,esin_hist,esin_bin_square,esin_hist_square = histogram(state_total[(i-1)*5+5,:],50)
  ax.plot(esin_bin_square,esin_hist_square ./maximum(esin_hist_square),color=cp[i],linestyle=":",linewidth=3,label=L"$e\sin{\omega}$")
#  ax.plot(esin_bin_square,esin_hist_square ./maximum(esin_hist_square),color=cp[i],linestyle=":",linewidth=3)
  ax.plot([0,0],[0,1.05],linestyle="--",color=cp[i],linewidth=2)
  ax.legend();
  ax.axis([-0.0175,0.0175,0,1.05]); ax.annotate(string("(",planet[i],")"),xy=[-0.014;0.8])
  println(planet[i]," ",@sprintf("%6.4f",ecos0),"+-",@sprintf("%6.4f",sqrt(cov_save[5i-1,5i-1])),
     " ",@sprintf("%6.4f",esin0),"+-",@sprintf("%6.4f",sqrt(cov_save[5i,5i])))
  ax.grid(linestyle=":")
end
ax = axes[8]
#ax.axis("off")
# Plot the prior:
include("compute_ecc_prior.jl")
ax.grid(linestyle=":")
tight_layout()
subplots_adjust(hspace = 0,wspace=0)
savefig("../T1_eccentricity_vectors_likelihood_profile_hmc.pdf",bbox_inches="tight")
#read(stdin,Char)

# Make a plot of eccentricity histogram:

clf()
for i=1:7
  ecc_bin,ecc_hist,ecc_bin_square,ecc_hist_square = histogram(sqrt.(state_total[(i-1)*5+4,:].^2+state_total[(i-1)*5+5,:].^2),50)
  plot(ecc_bin_square,ecc_hist_square./maximum(ecc_hist_square),color=cp[i],linewidth=3,label=planet[i])
end
xlabel("Eccentricity",fontsize=15)
ylabel("Probability",fontsize=15)
legend(fontsize=15)
axis([0,0.015,0,1.05])
xticks(fontsize=15)
yticks(fontsize=15)
#read(stdin,Char)

savefig("../eccentricity_posterior.pdf",bbox_inches="tight")

#read(stdin,Char)

fig,axes = subplots(4,2)
for i=1:7
  ax = axes[i]
  # Plot period:
  P0 = elements_grid_all[(i-1)*5+2,i+1,2,16]; x= elements_grid_all[(i-1)*5+2,i+1,2,:] .- P0
  prob = exp.(-0.5*(chi_grid_all[(i-1)*5+4,:] .-chi_grid_all[(i-1)*5+4,16]))
  sigP = sqrt(cov_save[(i-1)*5+2,(i-1)*5+2])
  x2=collect(linearspace(-4sigP,4sigP,1000))
  ax.plot(x,prob,label="Period",color=cp[i])
  prob2 = exp.(-0.5 .*x2.^2 ./cov_save[(i-1)*5+2,(i-1)*5+2])
  ax.plot(x2,prob2,color=cp[i],alpha=0.3)
  t0 = elements_grid_all[(i-1)*5+3,i+1,3,16]
  prob = exp.(-0.5*(chi_grid_all[(i-1)*5+3,:] .-chi_grid_all[(i-1)*5+3,16]))
  x = elements_grid_all[(i-1)*5+3,i+1,3,:] .- t0
  ax.plot(x,prob,label=L"$t_0$",linestyle="--",color=cp[i])
  sigt0 = sqrt(cov_save[(i-1)*5+3,(i-1)*5+3])
  x2=collect(linearspace(-4sigt0,4sigt0,1000))
  prob2 = exp.(-0.5 .*x2.^2 ./cov_save[(i-1)*5+3,(i-1)*5+3])
  ax.plot(x2,prob2,color=cp[i],linestyle="--",alpha=0.3)
  ax.plot([0,0],[0,1],linestyle=":",color=cp[i])
  ax.legend();
  ax.axis([-0.015,0.015,0,1]); ax.annotate(string("(",planet[i],")"),xy=[-0.014;0.8])
  println(planet[i]," ",@sprintf("%6.4f",P0),"+-",@sprintf("%6.4f",sigP),
     " ",@sprintf("%6.4f",t0),"+-",@sprintf("%6.4f",sigt0))
end

# Finally plot log(ndof) & V1e^{1/(2nu)}:

fig,axes = subplots(1,2,sharey="row")

ax = axes[1]

#ndof0 = ndof_grid_all[nparam-1,16]; x= ndof_grid_all[nparam-1,:]
#prob = exp.(-0.5*(chi_grid_all[nparam-1,:] .-chi_grid_all[nparam-1,16]))
#signdof = sqrt(cov_save[nparam-1,nparam-1])
#x2=collect(linearspace(-4signdof,4signdof,1000))
#ax.plot(x,prob,label=L"$\nu$",color=cp[n])
#prob2 = exp.(-0.5 .*x2.^2 ./cov_save[nparam-1,nparam-1])
#x2 .+= ndof0
#ax.plot(x2,prob2,color=cp[n],alpha=0.3)
#lgndof0 = log(ndof_grid_all[nparam-1,16]); x= log.(ndof_grid_all[nparam-1,:])
lgndof0 = lgndof_grid_all[nparam-1,16]; x= lgndof_grid_all[nparam-1,:]
prob = exp.(-0.5*(chi_grid_all[nparam-1,:] .-chi_grid_all[nparam-1,16]))
signdof = sqrt(cov_save[nparam-1,nparam-1])/exp(lgndof0)
x2=collect(linearspace(-4signdof,4signdof,1000))
ax.plot(x,prob,label=L"$\nu$",color=cp[n],linewidth=3)
prob2 = exp.(-0.5 .*x2.^2 ./(cov_save[nparam-1,nparam-1]/exp(2lgndof0)))
x2 .+= lgndof0
ax.plot(x2,prob2,color=cp[n],alpha=0.3,linewidth=3)
# Plot histogram of log(ndof) parameter:
lndof_bin,lndof_hist,lndof_bin_square,lndof_hist_square = histogram(state_total[36,:],50)
ax.plot(lndof_bin_square,lndof_hist_square./maximum(lndof_hist_square))
ax.set_xlabel(L"$\log{\nu}$",fontsize=15)
ax.set_ylabel("Probability",fontsize=15)

ax = axes[2]
V1exp2nuinv0 = V1exp2nuinv_grid_all[nparam,16]; x= V1exp2nuinv_grid_all[nparam,:]
prob = exp.(-0.5*(chi_grid_all[nparam,:] .-chi_grid_all[nparam,16]))
#sigV1exp2nuinv = sqrt(cov_save[nparam,nparam])
sigV1exp2nuinv =  0.09
x2=collect(linearspace(-4sigV1exp2nuinv,4sigV1exp2nuinv,1000))
ax.plot(x,prob,label=L"$\ln(V_1)$",color=cp[n+1],linewidth=3)
#prob2 = exp.(-0.5 .*x2.^2 ./cov_save[nparam,nparam])
prob2 = exp.(-0.5 .*x2.^2 ./sigV1exp2nuinv^2)
x2 .+= V1exp2nuinv0
ax.plot(x2,prob2,color=cp[n+1],alpha=0.3,linewidth=3)
# Plot histogram of this parameter:
V1expinv2nu_bin,V1expinv2nu_hist,V1expinv2nu_bin_square,V1expinv2nu_hist_square = histogram(state_total[37,:],50)
ax.plot(V1expinv2nu_bin_square,V1expinv2nu_hist_square./maximum(V1expinv2nu_hist_square))
ax.set_xlabel(L"$V_1 e^{1/(2\nu)}$",fontsize=15)
subplots_adjust(wspace=0)
savefig("../T1_students_params_transformed.pdf",bbox_inches="tight")


text:
sp = SplineParams(linspace(0, 5, 10), 0, 3)
cp = ChebParams(6, 2.0, 5.0)
lp = LinParams(linspace(0, 5, 10), 0)

@testset "Test subtype structure" begin
    @testset "Basis families" begin
        for T in [Cheb, Lin, Spline]
            @test T <: BasisFamily
        end
    end

    @testset "params" begin
        for T in [ChebParams, SplineParams, LinParams]
            @test T <: BasisParams
        end
    end

    @testset "basis structure representations" begin
        for T in [Tensor, Direct, Expanded]
            @test T <: AbstractBasisMatrixRep
        end
    end
end

@testset "test _param method" begin
    for (T, TP) in [(Cheb, ChebParams), (Lin, LinParams), (Spline, SplineParams)]
        @test BasisMatrices._param(T) == TP
        @test BasisMatrices._param(T()) == TP
    end
end

@testset "test more Param methods" begin

    @testset "Test extra outer constructors" begin
        @test sp == SplineParams(10, 0, 5, 3)
        @test lp == LinParams(10, 0, 5)
    end

    @testset "test equality of params" begin
        @test (  ==(sp, lp) ) == false
        @test (  ==(sp, cp) ) == false
        @test (  ==(cp, lp) ) == false

        @test   ==(sp, sp)
        @test   ==(cp, cp)
        @test   ==(lp, lp)
    end
end


text:
module GlobalTest

    #
    # Performance experiments with global variables.
    # see: https://docs.julialang.org/en/v1/manual/performance-tips/
    # This is in a module so I can declare a global const
    #

    using BenchmarkTools
    using InteractiveUtils

    const N = 10_000

    GLOB_VAR = rand(N)
    const GLOB_CONST = rand(N)
    # const GLOB_CONST_TYPED::Vector{Float64} = rand(N)
    # syntax: type declarations on global variables are not yet supported

    const T = eltype(GLOB_VAR)

    function fv(x::Vector)::Real
        s = 0.0
        for i in x
          s += i
        end
        s
    end

    function f_glob_const()::Real
        s = 0.0
        for i in GLOB_CONST
             s += i
        end
        s
    end

    function f_glob_var()::Real
        s = 0.0
        for i in GLOB_VAR
            s += i
        end
        s
    end

    function f_glob_var_typed()::Real
        s = 0.0
        for i in GLOB_VAR::Vector{T}
            s += i
        end
        s
    end

    function f_glob_var_semi_typed()::Real
        s = 0.0
        for i in GLOB_VAR::Vector
            s += i
        end
        s
    end

    function f_glob_var_demi_typed()::Real
        s = 0.0
        for i in GLOB_VAR::Vector{<:Real}
            s += i
        end
        s
    end

    function do_all()::BenchmarkGroup
        t = BenchmarkGroup()
        t[:fv_glob_const] = @benchmark fv(GLOB_CONST)
        t[:fv_glob_var] = @benchmark fv(GLOB_VAR)
        t[:f_glob_const] = @benchmark f_glob_const()
        t[:f_glob_var] = @benchmark f_glob_var()
        t[:f_glob_var_typed] = @benchmark f_glob_var_typed()
        t[:f_glob_var_semi_typed] = @benchmark f_glob_var_semi_typed()
        t[:f_glob_var_demi_typed] = @benchmark f_glob_var_demi_typed()
        t
    end
end

GlobalTest.do_all()

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
include("libavcodec_h.jl")

include("avcodec.jl")

const AV_NOPTS_VALUE  =  reinterpret(Int64, 0x8000000000000000)

function AVFrame()
    ns = fieldnames(AVFrame)
    ts = AVFrame.types
    parms = [T<:Ptr ? C_NULL : zero(T) for T in ts]

    fmt_pos = findfirst(ns, :format)
    parms[fmt_pos] = -one(ts[fmt_pos])

    pts_pos = findfirst(ns, :pts)
    parms[pts_pos] = AV_NOPTS_VALUE

    AVFrame(parms...)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using Random
using Test
using GeoPhyInv
using LinearAlgebra

const JF=GeoPhyInv.Fdtd

function Δxtest(nx)
	(nx<6) && error("test requires nx>6")

	p=Random.randn(nx)
	p[1:2] .= 0.0  # it is necessary that starting models have these boundary conditions
	p[nx-1:nx] .= 0.0

	pout=JF.Δxsample(p)

	pout2=randn(nx)
	pout2=pout

	p2=JF.Δxsample(pout2)

	@test dot(pout2,pout) ≈ dot(p,p2)
end

Δxtest(11)


text:
nmrpipe_header_fields = [
  # name               type               word offset (4-byte words)
    ("FDMAGIC",        Float32,           0),
    ("FDFLTFORMAT",    Float32,           1),
    ("FDFLTORDER",     Float32,           2),
    ("FDDIMCOUNT",     Float32,           9),
    ("FDF3OBS",        Float32,           10),
    ("FDF3SW",         Float32,           11),
    ("FDF3ORIG",       Float32,           12),
    ("FDF3FTFLAG",     Float32,           13),
    ("FDPLANELOC",     Float32,           14),
    ("FDF3SIZE",       Float32,           15),
    ("FDF2LABEL",      FixedString{8},    16),
    ("FDF1LABEL",      FixedString{8},    18),
    ("FDF3LABEL",      FixedString{8},    20),
    ("FDF4LABEL",      FixedString{8},    22),
    ("FDDIMORDER1",    Float32,           24),
    ("FDDIMORDER2",    Float32,           25),
    ("FDDIMORDER3",    Float32,           26),
    ("FDDIMORDER4",    Float32,           27),
    ("FDF4OBS",        Float32,           28),
    ("FDF4SW",         Float32,           29),
    ("FDF4ORIG",       Float32,           30),
    ("FDF4FTFLAG",     Float32,           31),
    ("FDF4SIZE",       Float32,           32),
    ("FDF3APOD",       Float32,           50),
    ("FDF3QUADFLAG",   Float32,           51),
    ("FDF4APOD",       Float32,           53),
    ("FDF4QUADFLAG",   Float32,           54),
    ("FDF1QUADFLAG",   Float32,           55),
    ("FDF2QUADFLAG",   Float32,           56),
    ("FDPIPEFLAG",     Float32,           57),
    ("FDF3UNITS",      Float32,           58),
    ("FDF4UNITS",      Float32,           59),
    ("FDF3P0",         Float32,           60),
    ("FDF3P1",         Float32,           61),
    ("FDF4P0",         Float32,           62),
    ("FDF4P1",         Float32,           63),
    ("FDF2AQSIGN",     Float32,           64),
    ("FDPARTITION",    Float32,           65),
    ("FDF2CAR",        Float32,           66),
    ("FDF1CAR",        Float32,           67),
    ("FDF3CAR",        Float32,           68),
    ("FDF4CAR",        Float32,           69),
    ("FDUSER1",        Float32,           70),
    ("FDUSER2",        Float32,           71),
    ("FDUSER3",        Float32,           72),
    ("FDUSER4",        Float32,           73),
    ("FDUSER5",        Float32,           74),
    ("FDPIPECOUNT",    Float32,           75),
    ("FDFIRSTPLANE",   Float32,           77),
    ("FDLASTPLANE",    Float32,           78),
    ("FDF2CENTER",     Float32,           79),
    ("FDF1CENTER",     Float32,           80),
    ("FDF3CENTER",     Float32,           81),
    ("FDF4CENTER",     Float32,           82),
    ("FDF2APOD",       Float32,           95),
    ("FDF2FTSIZE",     Float32,           96),
    ("FDREALSIZE",     Float32,           97),
    ("FDF1FTSIZE",     Float32,           98),
    ("FDSIZE",         Float32,           99),
    ("FDF2SW",         Float32,           100),
    ("FDF2ORIG",       Float32,           101),
    ("FDQUADFLAG",     Float32,           106),
    ("FDF2ZF",         Float32,           108),
    ("FDF2P0",         Float32,           109),
    ("FDF2P1",         Float32,           110),
    ("FDF2LB",         Float32,           111),
    ("FDF2OBS",        Float32,           119),
    ("FDMCFLAG",       Float32,           135),
    ("FDF2UNITS",      Float32,           152),
    ("FDNOISE",        Float32,           153),
    ("FDTEMPERATURE",  Float32,           157),
    ("FDRANK",         Float32,           180),
    ("FDTAU",          Float32,           199),
    ("FDF3FTSIZE",     Float32,           200),
    ("FDF4FTSIZE",     Float32,           201),
    ("FDF1OBS",        Float32,           218),
    ("FDSPECNUM",      Float32,           219),
    ("FDF2FTFLAG",     Float32,           220),
    ("FDTRANSPOSED",   Float32,           221),
    ("FDF1FTFLAG",     Float32,           222),
    ("FDF1SW",         Float32,           229),
    ("FDF1UNITS",      Float32,           234),
    ("FDF1LB",         Float32,           243),
    ("FDF1P0",         Float32,           245),
    ("FDF1P1",         Float32,           246),
    ("FDMAX",          Float32,           247),
    ("FDMIN",          Float32,           248),
    ("FDF1ORIG",       Float32,           249),
    ("FDSCALEFLAG",    Float32,           250),
    ("FDDISPMAX",      Float32,           251),
    ("FDDISPMIN",      Float32,           252),
    ("FD2DPHASE",      Float32,           256),
    ("FDF2X1",         Float32,           257),
    ("FDF2XN",         Float32,           258),
    ("FDF1X1",         Float32,           259),
    ("FDF1XN",         Float32,           260),
    ("FDF3X1",         Float32,           261),
    ("FDF3XN",         Float32,           262),
    ("FDF4X1",         Float32,           263),
    ("FDF4XN",         Float32,           264),
    ("FDHOURS",        Float32,           283),
    ("FDMINS",         Float32,           284),
    ("FDSECS",         Float32,           285),
    ("FDSRCNAME",      FixedString{16},   286),
    ("FDUSERNAME",     FixedString{16},   290),
    ("FDMONTH",        Float32,           294),
    ("FDDAY",          Float32,           295),
    ("FDYEAR",         Float32,           296),
    ("FDTITLE",        FixedString{60},   297),
    ("FDCOMMENT",      FixedString{160},  312),
    ("FDLASTBLOCK",    Float32,           359),
    ("FDCONTBLOCK",    Float32,           360),
    ("FDBASEBLOCK",    Float32,           361),
    ("FDPEAKBLOCK",    Float32,           362),
    ("FDBMAPBLOCK",    Float32,           363),
    ("FDHISTBLOCK",    Float32,           364),
    ("FD1DBLOCK",      Float32,           365),
    ("FDF2TDSIZE",     Float32,           386),
    ("FDF1TDSIZE",     Float32,           387),
    ("FDF3TDSIZE",     Float32,           388),
    ("FDF4TDSIZE",     Float32,           389),
    ("FD2DVIRGIN",     Float32,           399),
    ("FDF3APODCODE",   Float32,           400),
    ("FDF3APODQ1",     Float32,           401),
    ("FDF3APODQ2",     Float32,           402),
    ("FDF3APODQ3",     Float32,           403),
    ("FDF3C1",         Float32,           404),
    ("FDF4APODCODE",   Float32,           405),
    ("FDF4APODQ1",     Float32,           406),
    ("FDF4APODQ2",     Float32,           407),
    ("FDF4APODQ3",     Float32,           408),
    ("FDF4C1",         Float32,           409),
    ("FDF2APODCODE",   Float32,           413),
    ("FDF1APODCODE",   Float32,           414),
    ("FDF2APODQ1",     Float32,           415),
    ("FDF2APODQ2",     Float32,           416),
    ("FDF2APODQ3",     Float32,           417),
    ("FDF2C1",         Float32,           418),
    ("FDF1APODQ1",     Float32,           420),
    ("FDF1APODQ2",     Float32,           421),
    ("FDF1APODQ3",     Float32,           422),
    ("FDF1C1",         Float32,           423),
    ("FDF1APOD",       Float32,           428),
    ("FDF1ZF",         Float32,           437),
    ("FDF3ZF",         Float32,           438),
    ("FDF4ZF",         Float32,           439),
    ("FDFILECOUNT",    Float32,           442),
    ("FDSLICECOUNT",   Float32,           443),
    ("FDOPERNAME",     FixedString{32},   464),
    ("FDF1AQSIGN",     Float32,           475),
    ("FDF3AQSIGN",     Float32,           476),
    ("FDF4AQSIGN",     Float32,           477),
    ("FDF2OFFPPM",     Float32,           480),
    ("FDF1OFFPPM",     Float32,           481),
    ("FDF3OFFPPM",     Float32,           482),
    ("FDF4OFFPPM",     Float32,           483),
]

function write_nmrpipe_data(io, signal, indirect_axis, direct_axis)
    @assert ndims(signal) == 2
    for i = size(signal, indirect_axis()):-1:1
        fid = signal[indirect_axis(i), direct_axis(:)]
        write(io, Float32.(real.(fid)))
        write(io, Float32.(imag.(fid)))
    end
end

function write_nmrpipe_data(io, signal, direct_axis)
    @assert ndims(signal) == 1
    write(io, Float32.(real.(signal)))
    write(io, Float32.(imag.(signal)))
end

#=
# TODO unified 1D & 2D stuff
function add_nmrpipe_header_axis_info(header, ax, obs_freq, car_ppm)
    # Single axis shared stuff
end

function add_nmrpipe_header_axes(header, signal, indirect_axis, direct_axis)
    # 1D dimension-specific header stuff
end

function add_nmrpipe_header_axes(header, signal, indirect_axis, direct_axis)
    # 2D dimension-specific header stuff
end
=#

"""
    save_nmrpipe(io, signal, axes; frequency, ref_freq_offset)

Convert twix into NMR Pipe format to be loaded into third party analytical
software, such as MNova. This currently supports up to 2 dimension.

At minimal, for each dimension we should have:
- Signal data
- Axes
- Observation frequency (frequency)
- Carrier frequency (`ref_freq_offset`)

`ref_freq_offset` is a tuple containing the relative frequency offset between
obs_freq and frequency which will be set to 0 on the ppm scale.
"""
function save_nmrpipe(io::IO, signal, axes; frequency, ref_freq_offset)
    ndim = ndims(signal)

    if ndim != 2
        error("Only 2 dimensional data is currently supported")
    end

    # f1 and f2 bandwidth
    indirect_axis, direct_axis = axes
    f1_bandwidth = ustrip(u"Hz", 1.0/step(AxisArrays.axes(signal, indirect_axis).val))
    f2_bandwidth = ustrip(u"Hz", 1.0/step(AxisArrays.axes(signal, direct_axis).val))

    # frequency
    obs_freq       = ustrip(u"MHz", frequency)
    car_f1, car_f2 = ustrip.(u"ppm", ref_freq_offset)
    is_complex_data = true # Our data is complex data
    encoding       = "magnitude"
    time_domain    = true # Our data is always time domain

    # Construct dictionary of NMR Pipe header
    data = Dict()

    # Assign default values
    data["FDF1CENTER"] = 1.
    data["FDF2CENTER"] = 1.
    data["FDF3CENTER"] = 1.
    data["FDF4CENTER"] = 1.

    data["FDF3SIZE"] = 1.
    data["FDF4SIZE"] = 1.

    data["FDF1QUADFLAG"] = 1.
    data["FDF2QUADFLAG"] = 1.
    data["FDF3QUADFLAG"] = 1.
    data["FDF4QUADFLAG"] = 1.

    data["FDSPECNUM"]   = 1.
    data["FDFILECOUNT"] = 1.
    data["FD2DVIRGIN"]  = 1.
    # dimention ordering

    data["FDDIMORDER1"] = 2.0
    data["FDDIMORDER2"] = 1.0
    data["FDDIMORDER3"] = 3.0
    data["FDDIMORDER4"] = 4.0

    data["FDF1LABEL"] = "Y"
    data["FDF2LABEL"] = "X"
    data["FDF3LABEL"] = "Z"
    data["FDF4LABEL"] = "A"

    # misc values
    # fdatap.h says this "Indicates IEEE floating point format"
    data["FDFLTFORMAT"] = read(IOBuffer(b"\xef\xeenO"), Float32)
    data["FDFLTORDER"]  = 2.345f0  # Byte order mark (as a Float32!)

    data["FDDIMCOUNT"]  = ndim

    data["FDF1SW"] = f1_bandwidth
    data["FDF2SW"] = ndim > 1 ? f2_bandwidth : nothing

    if ndim > 1
        for ax in 1:ndim
            data["FDF$(ax)OBS"] = obs_freq
            data["FDF$(ax)CAR"] = ax == 1 ? car_f1 : car_f2
            data["FDF$(ax)QUADFLAG"] = is_complex_data ? 0.0 : 1.0
            psize = is_complex_data && (ax != 1) ?
                    size(signal, indirect_axis)/2 : size(signal, direct_axis)/1
            osize = psize
            data["FDF$(ax)TDSIZE"] = psize
            data["FDF$(ax)FTFLAG"] = time_domain == true ? 0.0 : 1.0
            data["FDF$(ax)APOD"] = data["FDF$(ax)TDSIZE"]

            if encoding == "tppi"
                data["FD$(ax)DPHASE"] = 1.0
            elseif (encoding == "complex") || (encoding == "states") || (encoding == "states-tppi")
                data["FD$(ax)DPHASE"] = 2.0
            else
                data["FD$(ax)DPHASE"] = 0.0
            end

            data["FDF$(ax)CENTER"] = (ax == 1) || (data["FD$(ax)DPHASE"] != 1) ?
                                     Int(psize/2.0) + 1.0 : Int(psize/4.0) + 1
            data["FDF$(ax)ORIG"] = data["FDF$(ax)CAR"] * data["FDF$(ax)OBS"] - data["FDF$(ax)SW"] * (osize - data["FDF$(ax)CENTER"]) / osize

            if ax == 1
                data["FDSIZE"] = psize
                data["FDREALSIZE"] = psize
            end
            if ax == 2
                data["FDSPECNUM"] = size(signal, indirect_axis)
                data["FDF1AQSIGN"] = (encoding == "complex") || (encoding == "states") ? 0 : 16
            end
        end
    end

    # Write binary header
    nmrpipe_word_size = 4
    buf = IOBuffer(zeros(UInt8, 512*nmrpipe_word_size), write=true, truncate=false, read=true)
    for (key,type,offset) in nmrpipe_header_fields
        if haskey(data, key)
            seek(buf,  nmrpipe_word_size*offset)
            write(buf, type(data[key]))
        end
    end
    seek(buf, 0)
    write(io, buf)
    write_nmrpipe_data(io, signal, axes...)
    nothing
end

function save_nmrpipe(filename::AbstractString, args...; kws...)
    open(filename, "w") do io
        save_nmrpipe(io, args...; kws...)
    end
end

_open(f, io::IO, args...) = f(io)
_open(f, name::AbstractString, args...; kws...) = open(f, name, args...; kws...)

"""
    twix_to_nmrpipe(twix, nmrpipe_path)

`twix` is a path or open IO stream to the original twix data file.
`nmrpipe` is a path or open IO stream to the result file to be written.
"""
function twix_to_nmrpipe(twix, nmrpipe)
    _open(twix) do twix_io
        twix_data = load_twix(twix_io)
        expt = mr_load(twix_data)
        signal = simple_averaging(expt, downsample=2)

        frequency = standard_metadata(twix_data).frequency

        _open(nmrpipe, "w") do nmrpipe_io
            save_nmrpipe(nmrpipe_io, signal, (Axis{:time1}, Axis{:time2});
                         frequency=frequency,
                         ref_freq_offset=(_water_tms_offset, _water_tms_offset))
        end
    end
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
immutable Pade{T<:Number,S<:Number}
    p::Poly{T}
    q::Poly{S}
    var::Symbol
    function Pade(p::Poly{T},q::Poly{S})
        if p.var != q.var
            error("Polynomials must have same variable")
        end
        new(p, q, p.var)
    end
end
Pade{T<:Number,S<:Number}(p::Poly{T}, q::Poly{S}) = Pade{T,S}(p,q)

function Pade{T}(c::Poly{T},m::Int,n::Int)
    @assert m+n < length(c)
    rold,rnew = Poly([zeros(T,m+n+1);one(T)],c.var),Poly(c.a[1:m+n+1],c.var)
    uold,vold = Poly([one(T)],c.var),Poly([zero(T)],c.var)
    unew,vnew = vold,uold
    for i=1:n
        temp0,temp1,temp2 = rnew,unew,vnew
        q,rnew = divrem(rold,rnew)
        unew,vnew = uold-q*unew,vold-q*vnew
        rold,uold,vold = temp0,temp1,temp2

    end
    if vnew[0] == 0
        d = gcd(rnew,vnew)
        rnew = rnew ÷ d
        vnew = vnew ÷ d
    end
    Pade(rnew/vnew[0],vnew/vnew[0])
end
padeval{T}(PQ::Pade{T},x) = polyval(PQ.p,x)./polyval(PQ.q,x)


text:
include("header.jl")
include("dataset.jl")
include("network.jl")
include("callback.jl")

epochs = ProgressBar(iter:n_epoch);
loss_epoch = zeros(Float64, n_exp);
grad_norm = zeros(Float64, n_exp);
for epoch in epochs
    global p
    for i_exp in randperm(n_exp)
        if i_exp in l_val
            continue
        end
        grad = ForwardDiff.gradient(x -> loss_neuralode(x, i_exp), p)
        grad_norm[i_exp] = norm(grad, 2)
        if grad_norm[i_exp] > grad_max
            grad = grad ./ grad_norm[i_exp] .* grad_max
        end
        update!(opt, p, grad)
    end
    for i_exp = 1:n_exp
        loss_epoch[i_exp] = loss_neuralode(p, i_exp)
    end
    loss_train = mean(loss_epoch[l_train])
    loss_val = mean(loss_epoch[l_val])
    grad_mean = mean(grad_norm[l_train])
    set_description(
        epochs,
        string(
            @sprintf(
                "Loss train: %.2e val: %.2e grad: %.2e lr: %.1e",
                loss_train,
                loss_val,
                grad_mean,
                opt[1].eta
            )
        ),
    )
    cb(p, loss_train, loss_val, grad_mean)
end

conf["loss_train"] = minimum(l_loss_train)
conf["loss_val"] = minimum(l_loss_val)
YAML.write_file(config_path, conf)

for i_exp in randperm(n_exp)
    cbi(p, i_exp)
end

max_stars_repo_name: DENG-MIT/CRNN-Pyrolysis

max_issues_repo_name: DENG-MIT/CRNN-Pyrolysis

max_forks_repo_name: DENG-MIT/CRNN-Pyrolysis


text:
#=##############################################################################
# DESCRIPTION
Testing rotor-on-wing and wing-on-rotor modeling using a scaled-down Bertin's
wing (45-deg, swept-back, planar wing in Example 7.2, pp. 343 of Bertin's
Aerodynamics for Engineers) with two APC 10x7 propellers (one on each side).

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Dec 2019
  * License   : MIT
=###############################################################################

# ------------ MODULES ---------------------------------------------------------
# Load simulation engine
# import FLOWUnsteady
reload("FLOWUnsteady")
uns = FLOWUnsteady
vlm = uns.vlm

import GeometricTools
gt = GeometricTools

using PyPlot

noise = nothing
try
    # FLOWNoise https://github.com/byuflowlab/FLOWNoise (it requires PSU-WOPWOP)
    import FLOWNoise
    noise = FLOWNoise
catch e
    @warn("FLOWNoise not loaded: $e")
    println("FLOWNoise not loaded: $e")
end

# ------------ GLOBAL VARIABLES ------------------------------------------------
# Default path where to save data
extdrive_path = "/media/edoalvar/MyExtDrive/simulationdata7/"
# extdrive_path = "temps/"

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["singleprop", "isolatedwing", "blownwing", "postprocessing",
                                                                        "noise"]
    include("blownwing_"*header_name*".jl")
end

# ------------ DRIVERS ---------------------------------------------------------
function run_singleprop(; xfoil=true, prompt=true)
    singleprop(; xfoil=xfoil, save_path=extdrive_path*"fvs_singleprop04/",
                                                                  prompt=prompt)
end

function run_isolatedwing(; prompt=true)
    isolatedwing(; save_path=extdrive_path*"fvs_isolatedwing10/", prompt=prompt)
end

function run_blownwing(; xfoil=true, prompt=true)
    blownwing(; xfoil=xfoil, save_path=extdrive_path*"blownwing_test05/",
                                     VehicleType=uns.VLMVehicle, nrevs=10, prompt=prompt)
end

function run_noise(; )

    read_path = "/media/edoalvar/MyExtDrive/simulationdata5/fvs_blownwing11_alloutputs_noise/aero_rawoutputs/"
    save_path = "/media/edoalvar/MyExtDrive/simulationdata5/fvs_blownwing11_alloutputs_noise/test_noise06/"
    postprocessing_noise(  save_path;
                            # ---------- AERO SIMULATION INFO ----------
                            num0=2,                         # First time step to analyze
                            nrevs=29,                       # Number of revs to analyze
                            # ---------- OBSERVERS -------------------------
                            loading=true,                   # Include loading pressure
                            sph_R=1.2*1.693/2, sph_nR=24, sph_ntht=24, # Sphere definition
                            sph_nphi=24, sph_phimax=360,
                            sph_rotation=[90, 0, 0],
                            sph_C=[1.693/4, 0, 0],
                            microphoneX=nothing,
                            ww_nummin=72*15,
                            ww_numless=72*11-2,
                            # ---------- INPUT OPTIONS ---------------------
                            read_path=read_path,
                            # ---------- OUTPUT OPTIONS --------------------
                            prompt=true, debug_paraview=true,
                            debuglvl=1,                     # WW debug level
                        )
end


text:
# Monte Carlo estimation of the out-of-sample cost of a given first-stage decision

function estimateSolnQuality(z_soln::Array{Float64},covariate_obs::Array{Float64},degree::Float64,numMCScenarios::Int64,numMCReplicates::Int64,coeff_true::Array{Float64,2})

	objEstimates = zeros(Float64,numMCReplicates)

	for mc = 1:numMCReplicates
		srand(randomSeeds_MC[mc])
		demand_scen_MC = generateTrueCondScenarios(numMCScenarios,covariate_obs,degree,coeff_true)
		objEstimates[mc] = estimateCostOfSoln(z_soln,demand_scen_MC)
	end

	return objEstimates
end

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
include("test_interpolations.jl")

try
    import CUDA
    global CUDA_INSTALLED = true
catch
    global CUDA_INSTALLED = false
end
CUDA_FUNCTIONAL = CUDA.functional()
if CUDA_INSTALLED && CUDA_FUNCTIONAL
    @info """
    Found a working CUDA.jl package, running GPU tests
    """
    include("test_cuda.jl")
else
    msg = """
    Skipping CUDA tests. Got
    CUDA_INSTALLED = $CUDA_INSTALLED
    CUDA_FUNCTIONAL = $CUDA_FUNCTIONAL
    """
    @warn msg
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m3(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m3!(A_copy)
end

@inline function exp_ps_m3!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=5
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 A3 P1 C0 P0
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A3 with operation: mult
    mul!(memslots3,memslots1,memslots2)
    # Computing P1 = x*I+x*A+x*A2+x*A3
    coeff1=0.16666666666666666
    coeff2=0.041666666666666664
    coeff3=0.008333333333333333
    coeff4=0.001388888888888889
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing C0 with operation: mult
    mul!(memslots5,memslots4,memslots3)
    # Deallocating P1 in slot 4
    # Deallocating A3 in slot 3
    # Computing P0 = x*I+x*A+x*A2+x*C0
    coeff1=1.0
    coeff2=1.0
    coeff3=0.5
    coeff4=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating A2 in slot 2
    # Deallocating C0 in slot 5
    return memslots1 # Returning P0
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# TODO add quantile estimation as indices
function buildarray_annual(C::ClimateTools.ClimGrid, dataout, numYears)
    lonsymbol = Symbol(C.dimension_dict["lon"])
    latsymbol = Symbol(C.dimension_dict["lat"])
    FD = AxisArray(dataout, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]), Axis{:time}(Dates.year.(DateTime.(numYears))))
    return FD
end

function buildarray_climato(C::ClimateTools.ClimGrid, dataout)
    lonsymbol = Symbol(C.dimension_dict["lon"])
    latsymbol = Symbol(C.dimension_dict["lat"])
    FD = AxisArray(dataout, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]))
    return FD
end

function buildarray_resample(C::ClimateTools.ClimGrid, dataout, newtime)
    lonsymbol = Symbol(C.dimension_dict["lon"])
    latsymbol = Symbol(C.dimension_dict["lat"])
    FD = AxisArray(dataout, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]), Axis{:time}(newtime))
    return FD
end

# function buildarray_orig(C::ClimateTools.ClimGrid, dataout)
#     lonsymbol = Symbol(C.dimension_dict["lon"])
#     latsymbol = Symbol(C.dimension_dict["lat"])
#     FD = AxisArray(dataout, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]), Axis{:time}(get_timevec(C)))
#     return FD
# end

"""
    prcp1(C::ClimGrid)

Annual number with preciptation >= 1 mm. This function returns a ClimGrid. Input data should be in mm.
"""
function prcp1(C::ClimGrid)
  @argcheck C[9] == "pr"
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  dataout  = zeros(Int64, (size(C.data, 1), size(C.data, 2), length(numYears)))
  datain   = (C.data.data)

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.mapreducedim!(t -> t >= 1, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable="prcp1", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    frostdays(C::ClimGrid)

FD, Number of frost days: Annual count of days when TN (daily minimum temperature) < 0 Celsius.

Let TN[i,j] be daily minimum temperature on day i in year j. Count the number of days where:

  TN[i,j] < 0 Celsius.
"""
function frostdays(C::ClimGrid)
  @argcheck C[9] == "tasmin"
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  dataout  = zeros(Int64, (size(C.data, 1), size(C.data, 2), length(numYears)))
  datain   = C.data.data

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.mapreducedim!(t -> t < 0, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable="frostdays", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    summerdays(C::ClimGrid)

SD, Number of summer days: Annual count of days when TX (daily maximum temperature) > 25 degree Celsius.

Let TX[i,j] be daily maximum temperature on day i in year j. Count the number of days where:

  TX[i,j] > 25 Celsius.
"""
function summerdays(C::ClimGrid)
  @argcheck C[9] == "tasmax"
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  dataout  = zeros(Int64, (size(C.data, 1), size(C.data, 2), length(numYears)))
  datain   = C.data.data

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.mapreducedim!(t -> t > 25, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable="summerdays", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    icingdays(C::ClimGrid)

ID, Number of summer days: Annual count of days when TX (daily maximum temperature) < 0 degree Celsius.

Let TX[i,j] be daily maximum temperature on day i in year j. Count the number of days where:

  TX[i,j] < 0 Celsius.
"""
function icingdays(C::ClimGrid)
  @argcheck C[9] == "tasmax"
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  dataout  = zeros(Int64, (size(C.data, 1), size(C.data, 2), length(numYears)))
  datain   = C.data.data

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.mapreducedim!(t -> t < 0, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable="icingdays", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    tropicalnights(C::ClimGrid)

TropicalNights, Number of tropical nights: Annual count of days when TN (daily maximum temperature) > 20 degree Celsius.

Let TN[i,j] be daily minimum temperature on day i in year j. Count the number of days where:

  TN[i,j] > 20 Celsius.
"""
function tropicalnights(C::ClimGrid)
  @argcheck C[9] == "tasmin"
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  dataout  = zeros(Int64, (size(C.data, 1), size(C.data, 2), length(numYears)))
  datain   = C.data.data

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.mapreducedim!(t -> t > 20, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable="tropicalnights", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    customthresover(C::ClimGrid, thres)

customthresover, annual number of days over a specified threshold.

Let TS[i,j] be a daily time serie value on day i in year j. Count the number of days where:

  TS[i,j] > thres.

"""
function customthresover(C::ClimGrid, thres)
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  dataout  = fill(NaN, (size(C.data, 1), size(C.data, 2), length(numYears)))
  datain   = C.data.data

  dataout_rshp = reshape(dataout, (size(dataout, 1)*size(dataout, 2), size(dataout, 3)))
  datain_rshp = reshape(datain, (size(datain, 1)*size(datain, 2), size(datain, 3)))

  # Indice calculation
  Threads.@threads for k in 1:size(datain_rshp, 1)
    for i in 1:length(numYears)
      idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
      dataout_rshp[k, i] = sum(datain_rshp[k, idx] .> thres)
      # Base.mapreducedim!(t -> t > thres, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
    end
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable=string("Days over ", thres), typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    customthresunder(C::ClimGrid, thres)

customthresunder, annual number of days under a specified threshold.

Let TS[i,j] be a daily time serie value on day i in year j. Count the number of days where:

  TS[i,j] < thres.

"""
function customthresunder(C::ClimGrid, thres)
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  datain   = C.data.data
  dataout  = zeros(Int64, (size(C.data, 1), size(C.data, 2), length(numYears)))

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.mapreducedim!(t -> t < thres, +, view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits="days", latunits=C.latunits, lonunits=C.lonunits, variable=string("Days under ", thres), typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    annualmax(C::ClimGrid)

Annual maximum of array data.

Let data[i,j] be daily time serie on day i in year j. Extract the highest value for year j.
"""
function annualmax(C::ClimGrid)
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  datain   = C.data.data
  dataout  = zeros(typeof(datain[1]), (size(C.data, 1), size(C.data, 2), length(numYears)))

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Statistics.maximum!(view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable="annualmax", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    annualmin(C::ClimGrid)

Annual minimum of array data.

Let data[i,j] be daily time serie on day i in year j. Extract the lowest value for year j.
"""
function annualmin(C::ClimGrid)
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  datain   = C.data.data
  dataout  = zeros(typeof(datain[1]), (size(C.data, 1), size(C.data, 2), length(numYears)))

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Statistics.minimum!(view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable="annualmin", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    annualmean(C::ClimGrid)

Annual mean of array data.

Let data[i,j] be daily time serie on day i in year j. Calculate the mean value for year j.
"""
function annualmean(C::ClimGrid)
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  datain   = C.data.data
  dataout  = zeros(typeof(datain[1]), (size(C.data, 1), size(C.data, 2), length(numYears)))

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Statistics.mean!(view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable="annualmean", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    annualsum(C::ClimGrid)

Annual sum of array data.

Let data[i,j] be daily time serie on day i in year j. Sums daily values for year j.
"""
function annualsum(C::ClimGrid)
  years    = Dates.year.(C.data[Axis{:time}][:])
  numYears = unique(years)
  datain   = C.data.data
  dataout  = zeros(typeof(datain[1]), (size(C.data, 1), size(C.data, 2), length(numYears)))

  # Indice calculation
  Threads.@threads for i in 1:length(numYears)
    idx = searchsortedfirst(years, numYears[i]):searchsortedlast(years, numYears[i])
    Base.sum!(view(dataout, :, :, i:i), view(datain, :,:, idx))
  end

  # Apply mask
  dataout = applymask(dataout, C.msk)

  # Build output AxisArray
  FD = buildarray_annual(C, dataout, numYears)

  # Return climGrid type containing the indice
  return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="year", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable="annualsum", typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end


text:
using Pkg
#  If NetABM is not available then it should be installed via
#  Uncomment if not installed
#  ] add https://github.com/ollin18/NetABM.jl

function useit(list::Array{Symbol})
    installed = [key for key in keys(Pkg.installed())]
    strpackages = @. string(list)
    uninstalled = setdiff(strpackages,installed)

    map(Pkg.add,uninstalled)
    for package ∈ list
        @eval using $package
    end
end

thep = [:LightGraphs, :StatsBase, :Plots, :GraphPlot, :Statistics, :LaTeXStrings, :Bootstrap, :NetABM, :DelimitedFiles]
useit(thep)

"""
Creates a LightGraph's SimpleGraph object from an adjacency list array. Returns the SimpleGraph and the ordered node sequence,
this way node names can be strings.
...
# Arguments
- `red::Array`: Adjacenc list.
...
"""
function lectura_uw(red)
    Nodes = union(unique(red[:,1]),unique(red[:,2]))
    g = SimpleGraph()
    last_node = Int64(length(Nodes))
    add_vertices!(g,last_node)
    for n in 1:size(red)[1]
        add_edge!(g,red[n,1],red[n,2])
        add_edge!(g,red[n,2],red[n,1])
    end
    return g, Nodes
end

"""
Generates a LFR network
...
# Arguments
- `N::Int`: Number of nodes.
- `k::Float`: Average degree.
- `maxk::Float`: Maximum degree.
- `mu::Float`: Mixing parameter.
- `t1::Float`:: Minus exponent for the degree sequence.
- `t2::Float`:: Minus exponent for the community size distribution.
- `minc::Int`: Minimum for the community sizes.
- `maxc::Int`: Maximum for the community sizes.
...
"""
function get_network(;N = 100, k = 20, maxk = 40, mu = 0.1, t1 = 2, t2 = 1, minc = 10, maxc = 40)
        the_root = pwd()
        run(`$the_root/src/binary_networks/benchmark -N $N -k $k -maxk $maxk
                -mu $mu -t1 $t1 -t2 $t2 -minc $minc -maxc $maxc`)
end

function get_network(;N = 100, k = 20, maxk = 40, mu = 0.1, t1 = 2, t2 = 1, minc = 10, maxc = 40)
        the_root = pwd()
        try
            run(`$the_root/src/binary_networks/benchmark -N $N -k $k -maxk $maxk -mu $mu`)
        catch
            N
            #  println("OK")
        end
end

get_network(N = 100, mu=0.1);

the_root = pwd()
the_net = readdlm(the_root*"/src/binary_networks/network.dat")
true_com = readdlm(the_root*"/src/binary_networks/community.dat")
true_com = Int64.(true_com)

g, Nodes = lectura_uw(the_net)

global num_agents = 100
global n_boot = 1000
global cil = 0.95
global j
global adoptstep = (1,[0.1, 0.1])
global maxI = Array{Float64}(undef,0)
global sumI = []

mixing = 0.1
@time let
    for mixing in 0.0:0.1:1
        for copera in 0.0:0.1:1
            for adap in 0.0:0.1:1
                global meantotI = Array{Float64}(undef,0)
                global Imean = Array{Float64}(undef,0)
                global Ilow = Array{Float64}(undef,0)
                global Ihigh = Array{Float64}(undef,0)
                global Itotall = Array{Float64}(undef,0)
                global Iall = Array{Array}(undef,0)
                global Sall = Array{Array}(undef,0)
                global Rall = Array{Array}(undef,0)
                for a in 1:10
                    if mixing > 0
                        get_network(N = num_agents, mu=mixing);
                        the_root = pwd()
                        the_net = readdlm(the_root*"/network.dat")
                        true_com = readdlm(the_root*"/community.dat")
                        true_com = Int64.(true_com)
                        g, Nodes = lectura_uw(the_net)
                    else
                        g = erdos_renyi(num_agents, 0.3)
                    end
                    agents = [Agent(i) for i = 1:num_agents];
                    init_demographics!(agents;states=["S","I"],initial=[1.0,0.0]);
                    pat_0 = rand(1:nv(g))
                    agents[pat_0].state = "I"
                    set_coop_agents!(agents;p_cop=copera);
                    set_adapt_agents!(agents;p_cop=adap);
                    map(x -> assign_contacts!(g,x), agents);
                    num_S = []
                    num_I = []
                    num_R = []
                    global oldS = ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents
                    global oldI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                    global totI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                    for i in 1:100
                        push!(num_S, ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
                        push!(num_I, ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
                        push!(num_R, ([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)

                        next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=true);
                        update_state!(agents);
                        update_effect_given_distance_coop!(agents,g,adoptstep[1],1,adoptstep[2]);

                        newI = oldS - num_S[end]
                        global totI = totI + newI
                        oldS = num_S[end]

                    end
                    push!(Iall,num_I)
                    push!(Itotall,totI)
                    GC.gc()
                end
                for step in eachindex(Iall[1])
                    bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
                    bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
                    bci1 = confint(bs1, BasicConfInt(cil))
                    bci2 = confint(bs2, BasicConfInt(cil))
                    push!(Imean,bci1[1][1])
                    push!(Ilow,bci2[1][2])
                    push!(Ihigh,bci2[1][3])
                end
                bs1 = bootstrap(mean, Itotall,BasicSampling(n_boot))
                bci1 = confint(bs1, BasicConfInt(cil))
                push!(sumI,bci1[1][1])
                push!(maxI,maximum(Imean))
            end
            GC.gc()
        end
    end
end

theme(:default)
using DelimitedFiles

plt = heatmap(
        xlabel = "Cooperative Agents",
        ylabel = "Adaptive Agents",
        title = "Maximum Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(maxI, (11, 11, 11))
z1 = z[:,:,1]
z1 = z[:,:,2]
plt = heatmap!(xs, ys, z1, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
plt
png(plt,"figs/heatmap_scaled_maximum_coop.png")
savefig(plt,"figs/heatmap_scaled_maximum_coop.pdf")
writedlm("data/max_inf_coop_lfr_er.csv",maxI)

plt = heatmap(
        xlabel = "Cooperative Agents",
        ylabel = "Adaptive Agents",
        title = "Total Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(sumI, (11, 11))
plt = heatmap!(xs, ys, z, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
plt
png(plt,"figs/heatmap_scaled_total_coop.png")
savefig(plt,"figs/heatmap_scaled_total_coop.pdf")
writedlm("data/tot_inf_coop_lfr_er.csv",sumI)


text:
include("../models/StochasticProcess.jl")

struct DeterministicModel <: StochasticProcess
    domAlias
    domCurve
    forAliases
    forAssetSpots
    forCurves
    index
end

function DeterministicModel(domAlias,domCurve,forAliases,forAssetSpots,forCurves)
    if !isnothing(forAliases)
        index = Dict([ (forAliases[k],k) for k = 1:size(forAliases)[1] ])
    else
        index = nothing
    end
    return DeterministicModel(domAlias,domCurve,forAliases,forAssetSpots,forCurves,index)
end

stateSize(self::DeterministicModel) = 0

factors(self::DeterministicModel) = 0

initialValues(self::DeterministicModel) = zeros(0)

evolve(self::DeterministicModel, t0, X0, dt, dW, X1) = nothing

function shortRateOverPeriod(self::DeterministicModel, t0, dt, X0, X1)
    B_d = discount(self.domCurve,t0) / discount(self.domCurve,t0 + dt)
    return log(B_d) / dt
end

zeroBondVolatility(self::DeterministicModel, t, T) = zeros(0)

zeroBondVolatilityPrime(self::DeterministicModel, t, T) = zeros(0)

numeraire(self::DeterministicModel, t, X) = 1.0 / discount(self.domCurve,t)

function asset(self::DeterministicModel, t, X, alias)
    k = self.index[alias]  # this should throw an exception if alias is unknown
    return self.forAssetSpots[k] * discount(self.forCurves[k],t) / discount(self.domCurve,t)
end

function zeroBond(self::DeterministicModel, t, T, X, alias)
    if isnothing(alias) || alias==self.domAlias
        return discount(self.domCurve,T) / discount(self.domCurve,t)
    end
    k = self.index[alias]  # this should throw an exception if alias is unknown
    return discount(self.forCurves[k],T) / discount(self.forCurves[k],t)
end

DcfModel(curve) = DeterministicModel(nothing,curve,nothing,nothing,nothing)

DcfModel(curve, alias) = DeterministicModel(alias,curve,nothing,nothing,nothing)

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
###################
# Dynamics Models #
###################

export fderiv_earth_orbit
"""
Compute the state derivative.

Arguments:
- `epc::Epoch`: Current epoch
- `x::Array{<:Real, 1}`: Satellite state vector
- `mass::Real`: Satellite mass [kg]
- `area_drag`: Velocity-facing area affected by drag. [m^2]
- `coef_drag`: Coefficient of drag [dimensionless]
- `area_srp`: Velocity-facing area affected by drag. [m^2]
- `coef_srp`: Coefficient of reflectivity [dimensionless]
- `n_grav::Integer`: Gravity model degree (Default: 20)
- `m_grav::Integer`: Gravity model order (Default: 20)
- `drag::Bool`: Include cannonball atomospheric drag in force model (Default: `true`)
- `srp::Bool`: Include flat-plate solar radiation pressure in force model (Default: `true`)
- `moon::Bool`: Include thridbody lunar gravity in force model (Default: `true`)
- `sun::Bool`: Include thirdbody solar in force model (Default: `true`)
- `relativity::Bool`: Include relativistic effects in force model (Default: `true`)

Returns:
- `dx::Array{<:Float64, 1}`: Satellite state derivative, velocity and accelerations [m; m/s]
"""
function fderiv_earth_orbit(epc::Epoch, x::Array{<:Real} ;
             mass::Real=1.0, area_drag::Real=1.0, coef_drag::Real=2.3,
             area_srp::Real=1.0, coef_srp::Real=1.8,
             n_grav::Integer=20, m_grav::Integer=20,
             drag::Bool=true, srp::Bool=true, moon::Bool=true, sun::Bool=true,
             relativity::Bool=true)

    # Extract position and velocity
    r = x[1:3]
    v = x[4:6]

    # Compute ECI to ECEF Transformation -> IAU2010 Theory
    PN = bias_precession_nutation(epc)
    E  = earth_rotation(epc)
    W  = polar_motion(epc)
    R  = W * E * PN

    # Compute geolocation
    geod = sECEFtoGEOD(R*r)

    # Compute sun and moon position
    r_sun  = sun_position(epc)
    r_moon = moon_position(epc)

    # Compute acceleration
    a = zeros(Float64, 3)

    a += accel_gravity(x, R, n_grav, m_grav)

    # Drag
    if drag
        # Use PN*x[1:3] to compute the satellite position in the true-of-date inertial frame
        rho = density_nrlmsise00(epc, geod)
        a  += accel_drag(x, rho, mass, area_drag, coef_drag, Array{Real, 2}(PN))
    end

    # Solar Radiation Pressure
    if srp
        nu = eclipse_conical(x, r_sun)
        a += nu*accel_srp(x, r_sun, mass, area_srp, coef_srp)
    end

    # Thirdbody
    if sun
        a += accel_thirdbody_sun(x, r_sun)
    end

    if moon
        a += accel_thirdbody_moon(x, r_moon)
    end

    # Relativity
    if relativity
        a += accel_relativity(x)
    end

    return vcat(v, a)
end

##############
# Propagator #
##############

export EarthInertialState
"""
Satellite orbit state vector using an Earth inertial state representation and
dynamics model.

If an initial state transition matrix is provided it will be used for propagation
if there is no state transition matrix then only the state will be propagated.

Attributes:
- `rk4::RK4` Internal numerical integrator used for state propagation
- `dt::Real` Default propagation time step. All steps will be this size unless
the state vector is requested to propagate to a time smaller than this step size,
which it will do.
- `epc::Epoch` Epoch of state
- `x::Array{Float64, 1}` State vector. Earth-centered inertial Cartesian state.
- `phi::Union{Nothing, Array{Float64, 2}}` State transition matrix, or the matrix
of partial derivatives of the state at the current time with respect to the
start of propagation.

The following force model parametters can be set as keyword arguments

Parameters:
- `mass::Real`: Satellite mass [kg]
- `area_drag`: Velocity-facing area affected by drag. [m^2]
- `coef_drag`: Coefficient of drag [dimensionless]
- `area_srp`: Velocity-facing area affected by drag. [m^2]
- `coef_srp`: Coefficient of reflectivity [dimensionless]
- `n_grav::Integer`: Gravity model degree (Default: 20)
- `m_grav::Integer`: Gravity model order (Default: 20)
- `drag::Bool`: Include cannonball atomospheric drag in force model (Default: `true`)
- `srp::Bool`: Include flat-plate solar radiation pressure in force model (Default: `true`)
- `moon::Bool`: Include thridbody lunar gravity in force model (Default: `true`)
- `sun::Bool`: Include thirdbody solar in force model (Default: `true`)
- `relativity::Bool`: Include relativistic effects in force model (Default: `true`)

"""
mutable struct EarthInertialState
    rk4::RK4
    dt::Real
    epc::Epoch
    x::Array{Float64, 1}
    phi::Union{Nothing, Array{Float64, 2}}
end

function EarthInertialState(epc::Epoch, x::Array{<:Real, 1}, phi::Union{Nothing, Array{Float64, 2}}=nothing; dt::Real=30.0, kwargs...)
    rk4 = RK4(fderiv_earth_orbit; kwargs...)
    return EarthInertialState(rk4, dt, epc, x, phi)
end

export step!
"""
Step dynamics state the specified time.

Arguments:
- `state::EarthInertialState` State vector to propagate
- `dt::Real` Step size. If no value is provided will use default state stepsize
"""
function step!(state::EarthInertialState, dt::Real=0.0)
    # Set step size
    if dt == 0.0
        dt = state.h
    end

    if state.phi == nothing
        # If no state transition matrix is provided integrate
        # without propagating variational equations
        state.x = istep(state.rk4, state.epc, dt, state.x)
    else
        state.x, state.phi = istep(state.rk4, state.epc, dt, state.x, state.phi)
    end

    # Advance state propagator time
    state.epc += dt
end

export stepto!
"""
Step dynamics state to the specified time. Will take as many internal steps as
required to advance the propagator to the correct time. No internal step will
exceed the size specified in the state propagator.

Arguments:
- `state::EarthInertialState` State vector to propagate
- `time::Union{Real, Epoch}` Time to propagate internal state too. Can be either
a real number to advance the state by or the Epoch
"""
function stepto!(state::EarthInertialState, time::Union{Real, Epoch}=0.0)
    if typeof(time) <: Real
        time = state.epc + time
    end

    if time < state.epc
        @warn "Propagation time is before current state epoch. Will integrate dyanmics in reverse."
    end

    # Propagate state
    while state.epc < time
        # Set propagation state size
        dt = min(state.dt, time - state.epc)

        # Propagate state
        step!(state, dt)
    end
end

export sim!
"""
Simulate the sttate dynamics to the specified time. Takes same inputs as `stepto!`
but instead of just updating the state vector to the specified time, this
function also returns the timesteps, state values, and state transition
matrices for each time step.

Arguments:
- `state::EarthInertialState` State vector to propagate
- `time::Union{Real, Epoch}` Time to propagate internal state too. Can be either
a real number to advance the state by or the Epoch

Returns:
- `t::Array{Float64, 1}` Elapsed time as a scalar from the initial simulation epoch
- `epc::Array{Epoch, 1}` Epoch at each timestep
- `x::Array{Float64, 2}` State vectors at each time step. Time is along second axis
- `Phi::Array{Float64, 2}` Stacked array of state transition matrices
"""
function sim!(state::EarthInertialState, time::Union{Real, Epoch}=0.0)
    if typeof(time) <: Real
        time = state.epc + time
    end

    if time < state.epc
        @warn "Propagation time is before current state epoch. Will integrate dyanmics in reverse."
    end

    # Number of state variables
    n = length(state.x)

    # Compute number of time steps in propagation
    n_steps = ceil(Int, (time - state.epc)/state.dt) + 1

    # Initialize containers to store output
    t   = zeros(Float64, n_steps)
    epc = Array{Epoch}(undef, n_steps)
    x   = zeros(Float64, n, n_steps)
    A   = zeros(Float64, n*n_steps, n)

    # Save initial state
    idx = 1

    epc[idx]  = state.epc
    t[idx]    = epc[idx] - epc[1]
    x[:, idx] = state.x

    if state.phi != nothing
        A[(1+n*(idx-1)):(n*idx), :] = state.phi
    end

    # Propagate state
    while state.epc < time
        # Increase storage index
        idx += 1

        # Set propagation state size
        dt = min(state.dt, time - state.epc)

        # Propagate state
        step!(state, dt)

        # Save output
        epc[idx]  = state.epc
        t[idx]    = epc[idx] - epc[1]
        x[:, idx] = state.x
        if state.phi != nothing
            A[(1+n*(idx-1)):(n*idx), :] = state.phi
        end
    end

    # Output
    if state.phi != nothing
        # Return with state transition matrix if used
        return t, epc, x, A
    end

    # Default return without STM
    return t, epc, x
end

export reinit!
"""
Reinitialize State transition matrix to identity at the current time step.

Used to reinitialize the state transition matrix value to identity
"""
function reinit!(state::EarthInertialState)
    state.phi = diagm(0 => ones(Float64, length(state.x)))
end


text:
# TODO: Remove when done with this.

import Pkg; Pkg.activate("../")
using CytofDiffDensity; const cdd = CytofDiffDensity
using Distributions
using MCMC
using StatsPlots
import Random

function compute_post_density(chain::AbstractVector, ygrid)
  function compute_pdf(c::NamedTuple)
    comps = SkewT.(c.mu, c.sigma, c.nu, c.phi)
    return pdf.(MixtureModel(comps, c.eta), ygrid)
  end
  return compute_pdf.(chain)
end

# Normal Ordered Prior.
K = 3
m1 = MixtureModel([Normal(-2, .3), Normal(0, .2), Normal(1, .4)], [.3, .4, .3])
y = rand(m1, 1000)
histogram(y, label=nothing, bins=100)
onm_prior = cdd.make_ordered_prior(y, K)
model1= MixSkewT(y, K, mu=onm_prior, tdist=false, skew=false)
spl = make_sampler(model1)
chain, metrics = mcmc(spl, 5000, nburn=100, thin=1)
@assert isapprox(mean.(m1.components), mean(getindex.(chain, :mu)), atol=0.2)
@assert isapprox(std.(m1.components), mean(getindex.(chain, :sigma)), atol=0.1)
plot(hcat(getindex.(chain, :mu)...)')

# SkewT Ordered Prior.
K = 3
Random.seed!(1)
m2 = MixtureModel([SkewT(-5, 1, 10, -7),
                   SkewT(0, 2, 10, -12),
                   SkewT(3, 3, 10, -10)], [.3, .4, .3])
y = rand(m2, 10000)
histogram(y, bins=100, label=nothing)

Random.seed!(0)
# onm_prior = cdd.make_ordered_prior(y, K, s=(maximum(y)-minimum(y))/2K)
onm_prior = cdd.make_ordered_prior(y, K)
model2 = MixSkewT(y, K, mu=onm_prior, skew=true, tdist=true)
spl = make_sampler(model2)
nburn, nsamps, thin = (3000, 1000, 1)
callback = make_callback(model2, nburn=nburn, nsamps=nsamps, thin=thin)
chain, metrics = mcmc(spl, nsamps, nburn=nburn, thin=thin, callback=callback)
plot(metrics[:loglike])

mean(getindex.(chain, :mu))
mean(getindex.(chain, :eta))
plot(hcat(getindex.(chain, :mu)...)', label=nothing)
hline!(getfield.(m2.components, :loc), label=nothing, ls=:dot, color=:black)
plot(hcat(getindex.(chain, :sigma)...)')
hline!(getfield.(m2.components, :scale), label=nothing, ls=:dot, color=:black)
plot(hcat(getindex.(chain, :nu)...)')
hline!(getfield.(m2.components, :df), label=nothing, ls=:dot, color=:black)
plot(hcat(getindex.(chain, :phi)...)')
hline!(getfield.(m2.components, :skew), label=nothing, ls=:dot, color=:black)
plot(hcat(getindex.(chain, :eta)...)')
hline!(m2.prior.p, label=nothing, ls=:dot, color=:black)

begin
  ygrid = cdd.make_ygrid(y)
  dens = compute_post_density(chain, ygrid)
  plot(ygrid, pdf.(m2, ygrid), label=nothing, color=:black)
  dens_lower = vec(MCMC.quantiles(hcat(dens...), 0.025, dims=2))
  dens_upper = vec(MCMC.quantiles(hcat(dens...), 0.975, dims=2))
  plot!(ygrid, dens_lower, fillrange=dens_upper, label=nothing, alpha=0.3)
  histogram!(y, normalize=true, la=0, alpha=.3)
end


text:
function generate_traindataset1d(datatyp, n_meshes, polydeg, data_input_size)
    n_dims = 1
    #datatyp = "NNPP"

    X = zeros(data_input_size, 0)
    # if datatyp == "NNPP"
    #     X = zeros(2,0)
    # elseif datatyp == "NNPPh"
    #     X = zeros(3,0)
    # elseif datatyp == "NNRH"
    #     X = zeros(5,0)
    # end
    Y = zeros(2,0)

    n_traindata = zeros(26)
    n_troubledcells = 0

    # loop over meshs
    for i in 1:n_meshes
        println("Mesh $i")
        #include("../meshes/$(n_dims)d/mesh$(i).jl")
        mesh = get_mesh_1d(i)
        leaf_cell_ids = leaf_cells(mesh.tree)
        n_elements = length(leaf_cell_ids)

        #loop over polydeg
        for r in polydeg #2:7
            n_nodes = r + 1
            nodes, _ = gauss_lobatto_nodes_weights(n_nodes)
            node_coordinates = zeros(1,n_nodes, n_elements)
            c2e = zeros(Int, length(mesh.tree))

            for element_id in 1:n_elements
                cell_id = leaf_cell_ids[element_id]
                dx = length_at_cell(mesh.tree, cell_id)
                # Calculate node coordinates
                for i in 1:n_nodes
                    node_coordinates[1, i, element_id] = (mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
                end
                c2e[leaf_cell_ids[element_id]] = element_id
            end

            #loop over functions
            for func in 1:26
                #println("Funktion $func")
                u(x) = trainfunction1d(func,x)
                if datatyp == "NNPP"
                    Xi = zeros(2,1)
                    input_data = zeros(2, n_elements)
                elseif datatyp == "NNPPh"
                    Xi = zeros(3,1)
                    input_data = zeros(3, n_elements)
                elseif datatyp == "NNRH"
                    Xi = zeros(5,1)
                    input_data = zeros(3, n_elements)
                end

                for element_id in 1:n_elements
                    data = node_coordinates[1, :, element_id]
                    input_data[:, element_id] .= generate_inputdata1d(u, data, r, datatyp)
                end

                for element_id in 1:n_elements
                    cell_id = leaf_cell_ids[element_id]
                    n_traindata[func] += 1
                    Yi = zeros(2,1)

                    # Create Xi
                    Xi .= generate_inputXi1d(element_id, input_data, c2e, mesh, datatyp)

                    h = length_at_cell(mesh.tree, cell_id)
                    xi = mesh.tree.coordinates[1, cell_id]

                    # Create label Yi
                    if (func > 18 && func <=26) && (  0 >= xi-(1/2)*h && 0 <= xi+(1/2)*h )
                        Yi=[1,0]
                    else
                        Yi=[0,1]
                    end

                    # append to X and Y
                    X = cat(X, Xi, dims=2)
                    Y = cat(Y, Yi, dims=2)
                end

            end
            #troubled cells

            for t in 1:1100
                ul = rand(Uniform(-1,1))
                ur = rand(Uniform(-1,1))
                x0 = rand(Uniform(-0.75, 0.75))
                u(x) = troubledcellfunctionstep1d(x, ul, ur, x0)

                if datatyp == "NNPP"
                    Xi = zeros(2,1)
                    input_data = zeros(2, n_elements)
                elseif datatyp == "NNPPh"
                    Xi = zeros(3,1)
                    input_data = zeros(3, n_elements)
                elseif datatyp == "NNRH"
                    Xi = zeros(5,1)
                    input_data = zeros(3, n_elements)
                end

                for element_id in 1:n_elements
                    data = node_coordinates[1, :, element_id]
                    input_data[:, element_id] .= generate_inputdata1d(u, data, r, datatyp)
                end

                for element_id in 1:n_elements
                    cell_id = leaf_cell_ids[element_id]
                    # Get cell length
                    dx = length_at_cell(mesh.tree, cell_id)
                    # for last functions just troubled cells
                    if good_cell1d(node_coordinates[1, 1, element_id], dx, x0)
                        continue
                    end
                    n_troubledcells += 1
                    Yi = zeros(2,1)

                    # Create Xi
                    Xi .= generate_inputXi1d(element_id, input_data, c2e, mesh, datatyp)

                    # Create label Yi
                    Yi = [1,0]

                    # append to X and Y
                    X = cat(X, Xi, dims=2)
                    Y = cat(Y, Yi, dims=2)

                end
            end
            println(size(X))
        end
        println(size(X))
    end

    println(size(X))

    # Scale data
    println("Scale data")
    for i in 1:size(X)[2]
        X[:,i]=X[:,i]./max(maximum(abs.(X[:,i])),1)
    end

    println("Safe data")
    h5open("datasets/$(n_dims)d/traindata1d$(datatyp).h5", "w") do file
        write(file, "X", X)
        write(file, "Y", Y)
    end
    return size(X)
end
#generate_traindataset()


text:
using BenchmarkTools, Test

grid=range(rand(), stop=randn(), length=10)

d=NamedD(grid,Srcs(10), [:a,:b, :c]);
d1=deepcopy(d);
	v=vec(d)

	@btime copyto!(v,d)

	randn!(v)

	@btime copyto!(d,v)
	@btime copyto!(d1,v)

	@test isequal(d,d1)
	@btime reverse!(d1)
	@btime reverse!(d1)
	@test isequal(d,d1)

	@btime randn!(d)
	@btime randn!(d1)

	@btime dot(d,d1)
	@test dot(d,d1)==dot(d1,d)

	@btime fill!(d, 0.0)
	@test iszero(d)

	println("##################################")

d=VNamedD(grid,SSrcs(10), Srcs(10), [:a,:b, :c]);
d1=deepcopy(d);
	v=vec(d)

	@btime copyto!(v,d)

	randn!(v)

	@btime copyto!(d,v)
	@btime copyto!(d1,v)

	@test isequal(d,d1)
	@btime reverse!(d1)
	@btime reverse!(d1)
	@test isequal(d,d1)

	@btime randn!(d)
	@btime randn!(d1)

	@btime dot(d,d1)
	@test dot(d,d1)==dot(d1,d)

	@btime fill!(d, 0.0)
	@test iszero(d)


text:
# See https://en.wikipedia.org/wiki/Complex_number

expressions = [
    :(sqrt(-1+0im))
    :(exp(π * im) + 1)
    :(sin(im))
    :(acos(2+0im))
    :(asin(2+0im))
    :(atan(1+2im))
    :(log(-1+0im))
    :(log(-im))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
using Sloshing
using Test

@testset "Sloshing.jl" begin
    # Write your tests here.
end


text:
module ParallelDifferentialEvolution
export diffevo, nullcb, logcb

import StatsBase
import Logging
import Statistics: mean, std
import Dates
import LatinHypercubeSampling

function nullcb(gen, x, im, f, con, etime)
end

function logcb(gen, x, im, f, con, etime)
    @info "generation: $gen; convergence: $con; elapsed time: $etime" x[im], f[im], mean(f), std(f)
    flush(stdout)
    flush(stderr)
end

# fo is the objective function (must be @everywhere for distributed use)
# we minimise the fitness
# bounds translates each parameter from [0, 1] to [min, max]
# d is dimensionality of problem
# F should be 0.4 < F < 1.0; small F speeds convergence but can cause premature convergence
# CR should be 0.1 < CR < 1.0; higher CR speeds convergence
# np should be between 5d and 10d, and must be at least 4
# see https://doi.org/10.1023/A:1008202821328
# rtol and atol set convergence tolerance: std(f) < (atol + rtol*abs(mean(f)))
# cb is a callback to use after eacg generation: defaulst to null, logcb prints stats
# fmap is the map to use to evaluate fitness: either `map` or `pmap`
# seeds can be a list of individuals to insert into the initial population
function diffevo(fo, d; F=0.8, CR=0.6, np=d*10,
                 maxiter=1000, rtol=1e-3, atol=1e-14, cb=nullcb, fmap=map,
                 seeds=nothing)
    plan, _ = LatinHypercubeSampling.LHCoptim(np, d, 1000)
    plan = LatinHypercubeSampling.scaleLHC(plan, repeat([(0.0,1.0)], d))
    x = mapslices(x->[x], plan, dims=2)
    if !isnothing(seeds)
        for (i,seed) in enumerate(seeds)
            x[i] = seed
        end
    end
    f = fmap(fo, x)
    im = argmin(f)
    m = x[im]
    trials = similar(x)
    stime = Dates.now()
    gen = 1
    while gen <= maxiter
        for j in 1:np
            ii = [k for k in 1:np if k != j]
            a, b, c = x[StatsBase.sample(ii, 3, replace=false)]
            dither = F + rand()*(1 - F)
            mut = clamp.(a .+ dither .* (b .- c), 0.0, 1.0)
            cp = rand(d) .< CR
            if !any(cp)
                cp[rand(1:d)] = true
            end
            trials[j] = ifelse.(cp, mut, x[j])
        end
        tf = fmap(fo, trials)
        for j in 1:np
            if tf[j] < f[j]
                f[j] = tf[j]
                x[j] = trials[j]
                if f[j] < f[im]
                    im = j
                    m = trials[j]
                end
            end
        end
        global con = std(f) / (atol + rtol*abs(mean(f)))
        global etime = floor(Dates.now() - stime, Dates.Second)
        cb(gen, x, im, f, con, etime)
        if con < 1.0
           @info "converged on generation $gen"
           break
        end
        gen += 1
    end
    if con >= 1.0
        @info "maximum number of $gen iterations reached"
    end
    return (F=F, CR=CR, np=np, maxiter=maxiter, gen=gen, rtol=rtol, atol=atol,
            m=m, f=f, fm=f[im], etime=etime, x=x, im=im, con=con)
end

end