text:
import Base: mean, std, median, var, cov
export corrcoef

#Mean
mean(a::AFAbstractArray) = af_mean(a)
mean{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_mean(a, dim-1))
mean{T<:Real}(a::AFAbstractArray{T}, weights::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_mean(a,weights,dim-1))

#Median
median(a::AFAbstractArray) = af_median(a)
median{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_median(a, dim-1))

#Standard Deviation
std(a::AFAbstractArray) = af_stdev(a)
std{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_stdev(a,dim-1))

#Variance
var(a::AFAbstractArray) = af_var(a)
var{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_var(a,dim-1))
var{T<:Real}(a::AFAbstractArray{T}, weights::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_var(a,weights,dim-1))

#Covariance
cov{T}(a::AFAbstractArray{T}, b::AFAbstractArray{T}) = AFArray{T}(af_cov(a,b))

#Correlation coefficient
corrcoef{T}(a::AFAbstractArray{T}, b::AFAbstractArray{T}) = af_corrcoef(a,b)


text:
using HTTP, CSV, JSON3, StructTypes
using JSONStat

"""Downloads a all datasets in a domain"""
function downloaddataset(domain::Int, directory::String;
                         dataset::Union{Nothing,Int}=nothing)

    info = domaindetail(domain; lang="PT")
    isnothing(dataset) ? datasetnum = info.num_datasets : datasetnum = dataset
    pagesize = min(100, info.num_series) # API does not accept more than 100 series per page_size
    pagenum = divrem(info.num_series, 100)
    f = info.label

    isdir(directory) || mkdir(directory)  # Si no hay directorio lo creamos

    for dᵢ in 1:datasetnum # A domain can have several datasets
        if datasetnum == 1
            io = open(joinpath(directory, string(f, ".csv")), "w") # Open file
        else
            io = open(joinpath(directory, string(f, "_ds", dᵢ, ".csv")), "w") # Open file
        end
        url = domaindatasets(domain).link.item[dᵢ].href
        ds = HTTP.get(string(url, "&page_size=", pagesize, "&page=1")).body
        pᵢ = 1
        while true # A big dataset is served in several 'pages'
            js = JSON3.read(ds)
            dt = JSONStat.read(ds)
            (pᵢ == 1) ? CSV.write(io, dt; delim=";", decimal=',', transform=(col, val) -> something(val, missing)) : CSV.write(io, dt; delim=";", decimal=',', append=true, transform=(col, val) -> something(val, missing)) # First page opens, next appends
            haskey(js.extension, :next_page) ? begin url = js.extension.next_page; pᵢ += 1 end : break
            ds = HTTP.get(url).body
        end
        close(io)
    end
end

"""Downloads a specfic set of series provided as a character-separated-values string. The series MUST belong to the same dataset"""
function downloadseries(series_id::String, directory::String; file::String="file")

    series = seriesdetails(series_id; lang="EN")
    dataset = series[1].dataset_id
    domain = series[1].domain_ids[1]

    for i in series
       (i.dataset_id != dataset) && error("Series must belong to the same dataset")
    end

    dataset = HTTP.get(bpstat_API_URL * "/domains/"*string(domain)*"/datasets/"*dataset*
                    "/?lang=EN&series_ids=" * series_id).body

    rw = JSONStat.read(dataset)

    isdir(directory) || mkdir(directory) # Si no hay directorio lo creamos

    io = open(joinpath(directory, file * ".csv"), "w")
    CSV.write(io, rw; delim=";",decimal=',',transform=(col, val) -> something(val, missing))
    return close(io)
end

# Friendly domain id, domain name printout
function listdomains()
    d = domains()
    for dᵢ in d
        println(dᵢ["id"], " ", dᵢ["label"])
    end
end

# Interface for BPStat Data API

const bpstat_API_URL = "https://bpstat.bportugal.pt/data/v1"

#-------------------------------------------------------------------------------------------------------------------------
# List Domains Endpoint
#-------------------------------------------------------------------------------------------------------------------------

function domains(; lang::String="EN")
    return JSON3.read(HTTP.get(bpstat_API_URL * "/domains/?lang=" * lang).body)
end

#-------------------------------------------------------------------------------------------------------------------------
# Domain Detail Endpoint
#-------------------------------------------------------------------------------------------------------------------------

function domaindetail(id::Int; lang::String="EN")
    return JSON3.read(HTTP.get(bpstat_API_URL * "/domains/" * string(id) * "/?lang=" * lang).body)
end

#-------------------------------------------------------------------------------------------------------------------------
# List domain datasets Endpoint
#-------------------------------------------------------------------------------------------------------------------------

function domaindatasets(id::Int; lang::String="EN")
    return JSON3.read(HTTP.get(bpstat_API_URL * "/domains/" * string(id) * "/datasets/?lang=" * lang*"&page_size=25").body)
end

#-------------------------------------------------------------------------------------------------------------------------
# Get Dataset (Observations) Endpoint
#-------------------------------------------------------------------------------------------------------------------------

function dataset(domain::Int; lang::String="EN", pagesize::Int=10, dataset::Int=1)

    domaindetail(domain)[Symbol("has_series")] == false &&
        return error("Selected Dataset has no series")

    url = domaindatasets(domain; lang).link.item[dataset].href
    return JSON3.read(HTTP.get(url * "&page_size=" * string(pagesize)).body)
end

#-------------------------------------------------------------------------------------------------------------------------
# List series details Endpoint
#-------------------------------------------------------------------------------------------------------------------------

function seriesdetails(series::String; lang::String="EN")
    return JSON3.read(HTTP.get(bpstat_API_URL * "/series/?lang=" * lang*"&series_ids="*series).body)
end

ext: jl


text:
using ITensors
using ITensorNetworks
using ITensorsVisualization

using ITensorNetworks: inds_network, project_boundary, contract_approx, insert_projectors
using ITensorNetworks: Models

model = Models.Model"ising"()
βc = Models.critical_point(model)
β = 1.001 * βc
@show β / βc

N = (3, 3)
ndims = length(N)

#space = [QN(0) => 2]
space = 2
tn_inds = inds_network(N...; linkdims=space, addtags="S=1/2")
A = Models.local_boltzmann_weight("ising", Val(ndims); β)

tn = itensor.((A,), tn_inds)

# Project periodic boundary indices
# onto the 1 state
state = 1
tn = project_boundary(tn, state)

cutoff = 1e-15
_maxdim = 100
dir = "top_to_bottom"

# Approximately contract the tensor network.
# Outputs a Vector of boundary MPS.
boundary_mps_top = contract_approx(tn; alg="boundary_mps", dir=dir, cutoff=cutoff, maxdim=_maxdim)

# Insert approximate projectors into rows of the network
tn_projected = insert_projectors(tn, boundary_mps_top; dir=dir, center=(:, 1))

# Outputs a tuple of the original tensor network
# and the tensors making up the projectors
tn_split, Pl, Pr = tn_projected

Pl_flat = reduce(vcat, Pl)
Pr_flat = reduce(vcat, Pr)
tn_projected_flat = mapreduce(vec, vcat, (tn_split, Pl_flat, Pr_flat))

@show noncommoninds(tn_projected_flat...)
@visualize *(tn_projected_flat...) contract=false
#@visualize *(tn_split...) contract=false

@disable_warn_order begin
  @show contract(tn_projected_flat)[] / contract(vec(tn))[]
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã ``2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using ExprSearch
using GrammarCounts
using Base.Test

#function create_grammar()
#  @grammar grammar begin
#    start = ex
#    ex = sum | product | (ex) | value
#    sum = Expr(:call, :+, ex, ex)
#    product = Expr(:call, :*, ex, ex)
#    value = :x | :y | digit
#    digit = 0:9
#  end
#  return grammar
#end

function create_grammar()
  @grammar grammar begin
    start = bin

    bin = always | eventually | implies | count
    always = Expr(:call, :G, bin_vec) #global
    eventually = Expr(:call, :F, bin_vec) #future
    implies = Expr(:call, :Y, bin_vec, bin_vec)
    count = Expr(:call, :ctlt, bin_vec, timestep) | Expr(:call, :ctle, bin_vec, timestep) | Expr(:call, :ctgt, bin_vec, timestep) | Expr(:call, :ctge, bin_vec, timestep) | Expr(:call, :cteq, bin_vec, timestep)

    #produces a bin_vec
    bin_vec = bin_feat | and | or | not  | eq | lt | lte | abseq | abslt | abslte | diff_eq | diff_lt | diff_lte | sign | absdiff_eq | absdiff_lt | absdiff_lte
    and = Expr(:call, :&, bin_vec, bin_vec)
    or = Expr(:call, :|, bin_vec, bin_vec)
    not = Expr(:call, :!, bin_vec)

    #equal
    eq = vrate_eq | altdiff_eq | chi_angle_eq | psi_angle_eq | sr_eq | timer_eq | psid_eq | v_eq | alt_eq | abs_altdiff_eq
    vrate_eq = Expr(:call, :eq, vrate_feat, vrate_val) | Expr(:call, :eq, vrate_feat, vrate_feat)
    altdiff_eq = Expr(:call, :eq, altdiff_feat, altdiff_val)
    chi_angle_eq = Expr(:call, :eq, chi_angle_feat, angle_val) | Expr(:call, :eq, chi_angle_feat, chi_angle_feat)
    psi_angle_eq = Expr(:call, :eq, psi_angle_feat, psi_angle_feat)
    sr_eq = Expr(:call, :eq, sr_feat, sr_val)
    timer_eq = Expr(:call, :eq, timer_feat, timer_val) | Expr(:call, :eq, timer_feat, timer_feat)
    psid_eq = Expr(:call, :eq, psid_feat, psid_val) | Expr(:call, :eq, psid_feat, psid_feat)
    v_eq = Expr(:call, :eq, v_feat, v_val) | Expr(:call, :eq, v_feat, v_feat)
    alt_eq = Expr(:call, :eq, alt_feat, alt_val) | Expr(:call, :eq, alt_feat, alt_feat)
    abs_altdiff_eq = Expr(:call, :eq, abs_altdiff_feat, abs_altdiff_val)

    #absolute equals
    abseq = vrate_abseq | altdiff_abseq | chi_angle_abseq | psi_angle_abseq | sr_abseq | timer_abseq | psid_abseq | v_abseq | alt_abseq
    vrate_abseq = Expr(:call, :abeq, vrate_feat, vrate_val) | Expr(:call, :abeq, vrate_feat, vrate_feat)
    altdiff_abseq = Expr(:call, :abeq, altdiff_feat, altdiff_val)
    chi_angle_abseq = Expr(:call, :abeq, chi_angle_feat, angle_val) | Expr(:call, :abeq, chi_angle_feat, chi_angle_feat)
    psi_angle_abseq = Expr(:call, :abeq, psi_angle_feat, psi_angle_feat)
    sr_abseq = Expr(:call, :abeq, sr_feat, sr_val)
    timer_abseq = Expr(:call, :abeq, timer_feat, timer_val) | Expr(:call, :abeq, timer_feat, timer_feat)
    psid_abseq = Expr(:call, :abeq, psid_feat, psid_val) | Expr(:call, :abeq, psid_feat, psid_feat)
    v_abseq = Expr(:call, :abeq, v_feat, v_val) | Expr(:call, :abeq, v_feat, v_feat)
    alt_abseq = Expr(:call, :abeq, alt_feat, alt_val) | Expr(:call, :abeq, alt_feat, alt_feat)

    #less than
    lt = vrate_lt | altdiff_lt  | chi_angle_lt | psi_angle_lt | sr_lt | timer_lt | psid_lt | v_lt | alt_lt | abs_altdiff_lt
    vrate_lt = Expr(:call, :lt, vrate_feat, vrate_val) | Expr(:call, :lt, vrate_feat, vrate_feat)
    altdiff_lt = Expr(:call, :lt, altdiff_feat, altdiff_val)
    chi_angle_lt = Expr(:call, :lt, chi_angle_feat, angle_val) | Expr(:call, :lt, chi_angle_feat, chi_angle_feat)
    psi_angle_lt = Expr(:call, :lt, psi_angle_feat, psi_angle_feat)
    sr_lt = Expr(:call, :lt, sr_feat, sr_val)
    timer_lt = Expr(:call, :lt, timer_feat, timer_val) | Expr(:call, :lt, timer_feat, timer_feat)
    psid_lt = Expr(:call, :lt, psid_feat, psid_val) | Expr(:call, :lt, psid_feat, psid_feat)
    v_lt = Expr(:call, :lt, v_feat, v_val) | Expr(:call, :lt, v_feat, v_feat)
    alt_lt = Expr(:call, :lt, alt_feat, alt_val) | Expr(:call, :lt, alt_feat, alt_feat)
    abs_altdiff_lt = Expr(:call, :lt, abs_altdiff_feat, abs_altdiff_val)

    #absolute less than
    abslt = vrate_abslt | altdiff_abslt  | chi_angle_abslt | psi_angle_abslt | sr_abslt | timer_abslt | psid_abslt | v_abslt | alt_abslt
    vrate_abslt = Expr(:call, :ablt, vrate_feat, vrate_val) | Expr(:call, :ablt, vrate_feat, vrate_feat)
    altdiff_abslt = Expr(:call, :ablt, altdiff_feat, altdiff_val)
    chi_angle_abslt = Expr(:call, :ablt, chi_angle_feat, angle_val) | Expr(:call, :ablt, chi_angle_feat, chi_angle_feat)
    psi_angle_abslt = Expr(:call, :ablt, psi_angle_feat, psi_angle_feat)
    sr_abslt = Expr(:call, :ablt, sr_feat, sr_val)
    timer_abslt = Expr(:call, :ablt, timer_feat, timer_val) | Expr(:call, :ablt, timer_feat, timer_feat)
    psid_abslt = Expr(:call, :ablt, psid_feat, psid_val) | Expr(:call, :ablt, psid_feat, psid_feat)
    v_abslt = Expr(:call, :ablt, v_feat, v_val) | Expr(:call, :ablt, v_feat, v_feat)
    alt_abslt = Expr(:call, :ablt, alt_feat, alt_val) | Expr(:call, :ablt, alt_feat, alt_feat)

    #less then or equal
    lte = vrate_lte | altdiff_lte | chi_angle_lte | psi_angle_lte | sr_lte | timer_lte | psid_lte | v_lte | alt_lte | abs_altdiff_lte
    vrate_lte = Expr(:call, :lte, vrate_feat, vrate_val) | Expr(:call, :lte, vrate_feat, vrate_feat)
    altdiff_lte = Expr(:call, :lte, altdiff_feat, altdiff_val)
    chi_angle_lte = Expr(:call, :lte, chi_angle_feat, angle_val) | Expr(:call, :lte, chi_angle_feat, chi_angle_feat)
    psi_angle_lte = Expr(:call, :lte, psi_angle_feat, psi_angle_feat)
    sr_lte = Expr(:call, :lte, sr_feat, sr_val)
    timer_lte = Expr(:call, :lte, timer_feat, timer_val) | Expr(:call, :lte, timer_feat, timer_feat)
    psid_lte = Expr(:call, :lte, psid_feat, psid_val) | Expr(:call, :lte, psid_feat, psid_feat)
    v_lte = Expr(:call, :lte, v_feat, v_val) | Expr(:call, :lte, v_feat, v_feat)
    alt_lte = Expr(:call, :lte, alt_feat, alt_val) | Expr(:call, :lte, alt_feat, alt_feat)
    abs_altdiff_lte = Expr(:call, :lte, abs_altdiff_feat, abs_altdiff_val)

    #abs less then or equal
    abslte = vrate_abslte | altdiff_abslte | chi_angle_abslte | psi_angle_abslte | sr_abslte | timer_abslte | psid_abslte | v_abslte | alt_abslte
    vrate_abslte = Expr(:call, :able, vrate_feat, vrate_val) | Expr(:call, :able, vrate_feat, vrate_feat)
    altdiff_abslte = Expr(:call, :able, altdiff_feat, altdiff_val)
    chi_angle_abslte = Expr(:call, :able, chi_angle_feat, angle_val) | Expr(:call, :able, chi_angle_feat, chi_angle_feat)
    psi_angle_abslte = Expr(:call, :able, psi_angle_feat, psi_angle_feat)
    sr_abslte = Expr(:call, :able, sr_feat, sr_val)
    timer_abslte = Expr(:call, :able, timer_feat, timer_val) | Expr(:call, :able, timer_feat, timer_feat)
    psid_abslte = Expr(:call, :able, psid_feat, psid_val) | Expr(:call, :able, psid_feat, psid_feat)
    v_abslte = Expr(:call, :able, v_feat, v_val) | Expr(:call, :able, v_feat, v_feat)
    alt_abslte = Expr(:call, :able, alt_feat, alt_val) | Expr(:call, :able, alt_feat, alt_feat)

    #sign
    sign = vrate_sign | chi_angle_sign | psid_sign
    vrate_sign = Expr(:call, :sn, vrate_feat, vrate_feat)
    chi_angle_sign = Expr(:call, :sn, chi_angle_feat, chi_angle_feat)
    psid_sign = Expr(:call, :sn, psid_feat, psid_feat)

    #difference is equal
    diff_eq = vrate_diff_eq | chi_angle_diff_eq | psi_angle_diff_eq | timer_diff_eq | psid_diff_eq | v_diff_eq
    vrate_diff_eq = Expr(:call, :dfeq, vrate_feat, vrate_feat, vrate_val)
    chi_angle_diff_eq = Expr(:call, :dfeq, chi_angle_feat, chi_angle_feat, angle_val)
    psi_angle_diff_eq = Expr(:call, :dfeq, psi_angle_feat, psi_angle_feat, angle_val)
    timer_diff_eq = Expr(:call, :dfeq, timer_feat, timer_feat, timer_val)
    psid_diff_eq = Expr(:call, :dfeq, psid_feat, psid_feat, psid_val)
    v_diff_eq = Expr(:call, :dfeq, v_feat, v_feat, v_val)

    #abs difference is equal
    absdiff_eq = vrate_absdiff_eq | chi_angle_absdiff_eq | psi_angle_absdiff_eq | timer_absdiff_eq | psid_absdiff_eq | v_absdiff_eq
    vrate_absdiff_eq = Expr(:call, :adfeq, vrate_feat, vrate_feat, vrate_val)
    chi_angle_absdiff_eq = Expr(:call, :adfeq, chi_angle_feat, chi_angle_feat, angle_val)
    psi_angle_absdiff_eq = Expr(:call, :adfeq, psi_angle_feat, psi_angle_feat, angle_val)
    timer_absdiff_eq = Expr(:call, :adfeq, timer_feat, timer_feat, timer_val)
    psid_absdiff_eq = Expr(:call, :adfeq, psid_feat, psid_feat, psid_val)
    v_absdiff_eq = Expr(:call, :adfeq, v_feat, v_feat, v_val)

    #difference is less than some value.  Values may not be appropriate after taking difference
    diff_lt = vrate_diff_lt | chi_angle_diff_lt | psi_angle_diff_lt | timer_diff_lt | psid_diff_lt | v_diff_lt
    vrate_diff_lt = Expr(:call, :dflt, vrate_feat, vrate_feat, vrate_val)
    chi_angle_diff_lt = Expr(:call, :dflt, chi_angle_feat, chi_angle_feat, angle_val)
    psi_angle_diff_lt = Expr(:call, :dflt, psi_angle_feat, psi_angle_feat, angle_val)
    timer_diff_lt = Expr(:call, :dflt, timer_feat, timer_feat, timer_val)
    psid_diff_lt = Expr(:call, :dflt, psid_feat, psid_feat, psid_val)
    v_diff_lt = Expr(:call, :dflt, v_feat, v_feat, v_val)

    #abs difference is less than some value.  Values may not be appropriate after taking difference
    absdiff_lt = vrate_absdiff_lt | chi_angle_absdiff_lt | psi_angle_absdiff_lt | timer_absdiff_lt | psid_absdiff_lt | v_absdiff_lt
    vrate_absdiff_lt = Expr(:call, :adflt, vrate_feat, vrate_feat, vrate_val)
    chi_angle_absdiff_lt = Expr(:call, :adflt, chi_angle_feat, chi_angle_feat, angle_val)
    psi_angle_absdiff_lt = Expr(:call, :adflt, psi_angle_feat, psi_angle_feat, angle_val)
    timer_absdiff_lt = Expr(:call, :adflt, timer_feat, timer_feat, timer_val)
    psid_absdiff_lt = Expr(:call, :adflt, psid_feat, psid_feat, psid_val)
    v_absdiff_lt = Expr(:call, :adflt, v_feat, v_feat, v_val)

    #difference is less than or equal to some value.  Values may not be appropriate after taking difference
    diff_lte = vrate_diff_lte | chi_angle_diff_lte | psi_angle_diff_lte | timer_diff_lte | psid_diff_lte | v_diff_lte
    vrate_diff_lte = Expr(:call, :dfle, vrate_feat, vrate_feat, vrate_val)
    chi_angle_diff_lte = Expr(:call, :dfle, chi_angle_feat, chi_angle_feat, angle_val)
    psi_angle_diff_lte = Expr(:call, :dfle, psi_angle_feat, psi_angle_feat, angle_val)
    timer_diff_lte = Expr(:call, :dfle, timer_feat, timer_feat, timer_val)
    psid_diff_lte = Expr(:call, :dfle, psid_feat, psid_feat, psid_val)
    v_diff_lte = Expr(:call, :dfle, v_feat, v_feat, v_val)

    #absolute difference is less than or equal to some value.
    absdiff_lte = vrate_absdiff_lte | chi_angle_absdiff_lte | psi_angle_absdiff_lte | timer_absdiff_lte | psid_absdiff_lte | v_absdiff_lte
    vrate_absdiff_lte = Expr(:call, :adfle, vrate_feat, vrate_feat, vrate_val)
    chi_angle_absdiff_lte = Expr(:call, :adfle, chi_angle_feat, chi_angle_feat, angle_val)
    psi_angle_absdiff_lte = Expr(:call, :adfle, psi_angle_feat, psi_angle_feat, angle_val)
    timer_absdiff_lte = Expr(:call, :adfle, timer_feat, timer_feat, timer_val)
    psid_absdiff_lte = Expr(:call, :adfle, psid_feat, psid_feat, psid_val)
    v_absdiff_lte = Expr(:call, :adfle, v_feat, v_feat, v_val)

    #read features
    bin_feat = Expr(:ref, :D, :(:), bin_feat_id)
    vrate_feat = Expr(:ref, :D, :(:), vrate_feat_id)
    altdiff_feat = Expr(:ref, :D, :(:), altdiff_feat_id)
    abs_altdiff_feat = Expr(:ref, :D, :(:), abs_altdiff_feat_id)
    angle_feat = Expr(:ref, :D, :(:), angle_feat_id)
    psi_angle_feat = Expr(:ref, :D, :(:), psi_angle_feat_id)
    chi_angle_feat = Expr(:ref, :D, :(:), chi_angle_feat_id)
    sr_feat = Expr(:ref, :D, :(:), sr_feat_id)
    timer_feat = Expr(:ref, :D, :(:), timer_feat_id)
    psid_feat = Expr(:ref, :D, :(:), psid_feat_id)
    v_feat = Expr(:ref, :D, :(:), v_feat_id)
    alt_feat = Expr(:ref, :D, :(:), alt_feat_id)

    #indices of each type
    bin_feat_id = 1 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 23 | 24 | 25 | 26 | 27 | 28 | 30 | 31 | 32 | 38 | 44 | 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 60 | 61 | 62 | 63 | 64 | 65 | 67 | 68 | 69 | 75
    vrate_feat_id = 2 | 22 | 34 | 39 | 59 | 71
    altdiff_feat_id = 3 | 40
    abs_altdiff_feat_id = 76
    angle_feat_id = psi_angle_feat_id | chi_angle_feat_id
    psi_angle_feat_id = 4 | 41
    chi_angle_feat_id = 6 | 43
    sr_feat_id = 5 | 42 | 77
    timer_feat_id = 33 | 70
    psid_feat_id = 35 | 72
    v_feat_id = 36 | 73
    alt_feat_id = 37 | 74

    #values of each type
    vrate_val = -50 | -40 | -30 | -20 | -10 | -1 | 0 | 1 | 10 | 20 | 30 | 40 | 50
    altdiff_val = -2000 | -1500 | -1000 | -500 | -250 | -100 | -50 | -25 | -10 | -5 | -1 | 0 | 1 | 5 | 10 | 25 | 50 | 100 | 250 | 500 | 1000 | 1500 | 2000
    abs_altdiff_val = 0 | 1 | 5 | 10 | 25 | 50 | 100 | 250 | 500 | 1000 | 1500 | 2000
    angle_val = -180 | -135 | -90 | -45 | 0 | 45 | 90 | 135 | 180
    sr_val = 30000 | 25000 | 20000 | 15000 | 10000 | 7500 | 5000 | 2500 | 1000 | 500 | 250 | 100 | 50 | 25 | 10 | 1 | 0
    timer_val = 0:5
    psid_val = -10:10
    v_val = 50 | 75 | 100 | 150 | 200 | 250 | 300 | 350 | 400 | 450 | 500
    alt_val = 1000 | 2000 | 3000 | 5000 | 10000 | 15000 | 18000 | 20000 | 25000 | 29000
    timestep = 0:50
  end

  return grammar
end

function test(; depth::Int64=20)
  grammar = create_grammar()
  rcnts = rulecounts(grammar, depth)
  dcnts = decisioncounts(grammar, rcnts)

  rcnts, dcnts
end

test()

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
"""
    sample_neighbors(g, nodes, K=-1; dir=:in, replace=false, dropnodes=false)

Sample neighboring edges of the given nodes and return the induced subgraph.
For each node, a number of inbound (or outbound when `dir = :out``) edges will be randomly chosen.
If `dropnodes=false`, the graph returned will then contain all the nodes in the original graph,
but only the sampled edges.

The returned graph will contain an edge feature `EID` corresponding to the id of the edge
in the original graph. If `dropnodes=true`, it will also contain a node feature `NID` with
the node ids in the original graph.

# Arguments

- `g`. The graph.
- `nodes`. A list of node IDs to sample neighbors from.
- `K`. The maximum number of edges to be sampled for each node.
       If -1, all the neighboring edges will be selected.
- `dir`. Determines whether to sample inbound (`:in`) or outbound (``:out`) edges (Default `:in`).
- `replace`. If `true`, sample with replacement.
- `dropnodes`. If `true`, the resulting subgraph will contain only the nodes involved in the sampled edges.

# Examples

```julia
julia> g = rand_graph(20, 100)
GNNGraph:
    num_nodes = 20
    num_edges = 100

julia> sample_neighbors(g, 2:3)
GNNGraph:
    num_nodes = 20
    num_edges = 9
    edata:
        EID => (9,)

julia> sg = sample_neighbors(g, 2:3, dropnodes=true)
GNNGraph:
    num_nodes = 10
    num_edges = 9
    ndata:
        NID => (10,)
    edata:
        EID => (9,)

julia> sg.ndata.NID
10-element Vector{Int64}:
  2
  3
 17
 14
 18
 15
 16
 20
  7
 10

julia> sample_neighbors(g, 2:3, 5, replace=true)
GNNGraph:
    num_nodes = 20
    num_edges = 10
    edata:
        EID => (10,)
```
"""
function sample_neighbors(g::GNNGraph{<:COO_T}, nodes, K=-1;
        dir=:in, replace=false, dropnodes=false)
    @assert dir ∈ (:in, :out)
    _, eidlist = adjacency_list(g, nodes; dir, with_eid=true)
    for i in 1:length(eidlist)
        if replace
            k = K > 0 ? K : length(eidlist[i])
        else
            k = K > 0 ? min(length(eidlist[i]), K) : length(eidlist[i])
        end
        eidlist[i] = StatsBase.sample(eidlist[i], k; replace)
    end
    eids = reduce(vcat, eidlist)
    s, t = edge_index(g)
    w = get_edge_weight(g)
    s = s[eids]
    t = t[eids]
    w = isnothing(w) ? nothing : w[eids]

    edata = getobs(g.edata, eids)
    edata = (edata..., EID = eids)

    num_edges = length(eids)

    if !dropnodes
        graph = (s, t, w)

        gnew = GNNGraph(graph,
                    g.num_nodes, num_edges, g.num_graphs,
                    g.graph_indicator,
                    g.ndata, edata, g.gdata)
    else
        nodes_other = dir == :in ? setdiff(s, nodes) : setdiff(t, nodes)
        nodes_all = [nodes; nodes_other]
        nodemap = Dict(n => i for (i, n) in enumerate(nodes_all))
        s = [nodemap[s] for s in s]
        t = [nodemap[t] for t in t]
        graph = (s, t, w)
        graph_indicator = g.graph_indicator !== nothing ? g.graph_indicator[nodes_all] : nothing
        num_nodes = length(nodes_all)
        ndata = getobs(g.ndata, nodes_all)
        ndata = (ndata..., NID = nodes_all)

        gnew = GNNGraph(graph,
                    num_nodes, num_edges, g.num_graphs,
                    graph_indicator,
                    ndata, edata, g.gdata)
    end
    return gnew
end


text:
################################################################################
# Copyright (C) Marcus Greiff 2020
#
# @file example_FOF_observer_sim.jl
# @author Marcus Greiff
# @date July 2020
#
# This example simply demonstrates the integration of the observer dynamics in
# the FOF attitude feedback. This is a bit special, as it does not adhere to
# the regular attitude kinematics, but involves a set of innovation terms. In
# addition, the numerical integration needs to be done in such a way that the
# attitude remains on SU(2), which is the internal representation of the filter
# memory in the numerical integration. As such, the function which is
# exemplified here considers the simulation of the dynamics
#
#     d(Q(t))/dt = Q(t) * [ deltaQ(t) / 2 ]_{SU(2)}^{\land}
#     d(W(t))/dt = deltaW(t)
#
# one step ahead with a time-step dt [s]. This is done by evaluating
#
#     Q(t + dt) = Q(t) * exp([ deltaQ(t) / 2 ]_{SU(2)}^{\land})
#     W(t + dt) = W(t) + dt * deltaW(t)
#
# whih is nontrivial to implement in C, hence verified extrnally by comparison
# to a simulation of the same dynamics using the Tsit5 method in the
# DifferentialEquations.jl stack.
#
# This example computes these to separate solutions and subsequently plots them
# on-top of eachother, showing that we compute the exact same solution in C as
# we do in Julia, and verifying that there are no discernable bugs in the
#
#     attitude_FOF_SO3_continuous_simulate_observer()
#
# function located in cont_attitude_FOF_SO3_continuous.c.
#
################################################################################
using DifferentialEquations
using LinearAlgebra
using LaTeXStrings
using AerialVehicleControl

showPlot = 1;     # Set to 0 for no plot, or any number in (1,2)
savePlot = false; # Set to true if the plot is to be saved

# Recompile the C-code
recompile()

# Wrapper for the calling of the control update law from C
function update_estimator(
    Qm::matrix_double_t,
    Wm::matrix_double_t,
    deltaQm::matrix_double_t,
    deltaWm::matrix_double_t,
    dt::Cdouble)
    status = ccall((:attitude_FOF_SO3_continuous_simulate_observer, AerialVehicleControl.CONT_LIB_PATH),
        Cint,
        (Ref{matrix_double_t},
         Ref{matrix_double_t},
         Ref{matrix_double_t},
         Ref{matrix_double_t},
         Cdouble,),
        Qm,
        Wm,
        deltaQm,
        deltaWm,
        dt)
    return status;
end

function odefun!(dx,x,param,t)
    # Extract the states
    q  = x[1:4];
    w  = x[5:7];

    we      = zeros(4,1);
    dw      = zeros(3,1);
    we[2:4] = w;
    dq      = quaternion_left_product(q) * we / 2.0;
    dw[1]   = cos(t)+1;
    dw[2]   = sin(2*t);
    dw[3]   = sin(t)/2;

    # Externalize state derivatives
    dx[1:4]   = dq;
    dx[5:7]   = dw;
end

function initialize_example()

    # Initial states of the physical system
    q0 = rand(4,1).-0.5;
    q0 = q0./norm(q0);
    w0 = 1*randn(3,1);

    x0    = vcat(q0, w0)
    tspan = (0.0,20.0)
    return x0, tspan
end

# Compute solution using DifferentialEquations.jl
x0, tspan    = initialize_example()
prob         = ODEProblem(odefun!, x0, tspan, []);
sol          = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12)

# Compute solution using the implementation in C by the first order Euler method
# implemented in the C-file
time  = collect(0:0.01:20)
Qtraj = zeros(4, length(time));
Wtraj = zeros(3, length(time));
Qtraj[:,1] = x0[1:4,1];
Wtraj[:,1] = x0[5:7,1];

# Run a constant time-step simulation using the update_estimator funciton
# implemented in the FOF attitude controller
for i = 2:length(time);

    dw      = zeros(3,1);
    dw[1]   = cos(time[i-1])+1;
    dw[2]   = sin(2*time[i-1]);
    dw[3]   = sin(time[i-1])/2;

    Q,  _,  Qm = get_matrix_and_copy(Qtraj[:,i-1]);
    W,  _,  Wm = get_matrix_and_copy(Wtraj[:,i-1]);
    dQ, _, dQm = get_matrix_and_copy(Wtraj[:,i-1]);
    dW, _, dWm = get_matrix_and_copy(dw);
    dt = Cdouble(time[i] - time[i - 1]);

    update_estimator(Qm, Wm, dQm, dWm, dt)

    Qtraj[:,i] = Q;
    Wtraj[:,i] = W;

end

# Using the Tsit5() method in Julia as a benchmark,  the solutions in computed
# in DifferentialEquations.jl can be compared to the solution computed using
# the numerical CG-intergration cheme implemented in the FOF controller.
using Plots
pq = plot(sol, vars = (0,1:4), color=:black, linewidth=2, xaxis="Time (t)", yaxis="Attitude [.]", label=[L"q(t)\;\;(computed\;in\;Julia)" nothing nothing nothing])
plot!(time, Qtraj',  color=:blue, style=:dash,  linewidth=2, xaxis="Time (t)", yaxis="Attitude [.]", label=[L"q(t)\;\;(computed\;in\;C)" nothing nothing nothing])
pw = plot(sol, vars = (0,5:7), color=:black, linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\omega(t)\;\;(computed\;in\;Julia)" nothing nothing])
plot!(time, Wtraj', color=:blue, style=:dash, linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\omega(t)\;\;(computed\;in\;C)" nothing nothing])
plot(pq, pw, layout=(2,1), size=(1000,750))


text:
using QuantumInfo
using Test

import LinearAlgebra
import QuantumInfo.eye

@testset "partial trace" begin
    for i=1:100
        da = round(Int,round(rand()*10))+1
        db = round(Int,round(rand()*10))+1
        ra = partialtrace(projector(randn(da*3)+1im*randn(da*3)),[da,3],2)
        rb = partialtrace(projector(randn(db*3)+1im*randn(db*3)),[db,3],2)
        @test isapprox(partialtrace(kron(ra,rb),[da,db],1),rb,atol=1e-15)
        @test isapprox(partialtrace(kron(ra,rb),[da,db],2),ra,atol=1e-15)
    end
end

@testset "representation conversion" begin
    id_choi = 1/2 .* [1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1];
    deph_choi = 1/2 .* [1 0 0 0.99;0 0 0 0;0 0 0 0;0.99 0 0 1];
    deph_liou = [1 0 0 0;0 0.99 0 0;0 0 0.99 0;0 0 0 1];

    # Kraus operators for a dephasing channel
    ρ = 0.01;
    M_0 = sqrt(1-ρ)*eye(2);
    M_1 = sqrt(ρ)*Matrix{ComplexF64}([1.0 0.0; 0.0 0.0]);
    M_2 = sqrt(ρ)*Matrix{ComplexF64}([0.0 0.0;0.0 1.0]);

    @test liou(eye(2),eye(2)) == eye(4)
    @test liou(eye(2)) == liou(eye(2),eye(2))
    @test mat(vec(eye(2))) == eye(2)
    @test mat([1,2,3,4]) == [1 3;2 4]
    @test vec([1 2; 3 4]) == [1,3,2,4]

    @test choi_liou_involution(choi_liou_involution(eye(4)))==eye(4)

    @test isapprox(choi2liou(id_choi),eye(4))
    @test isapprox(liou2choi(eye(4)),id_choi)
    @test isapprox(choi2kraus(id_choi),Matrix{ComplexF64}[zeros(2,2), zeros(2,2), zeros(2,2), eye(2)],atol=1e-7)
    @test isapprox(kraus2liou(Matrix{ComplexF64}[eye(2)]),eye(4))
    @test isapprox(liou2kraus(eye(4)),Matrix{ComplexF64}[zeros(2,2), zeros(2,2), zeros(2,2), eye(2)],atol=1e-7)
    @test isapprox(kraus2choi(Matrix{ComplexF64}[eye(2)]),id_choi)

    @test isapprox(kraus2choi([M_0, M_1, M_2]),deph_choi)
    @test isapprox(kraus2liou([M_0, M_1, M_2]),deph_liou)

    @test isapprox(liou2pauliliou(eye(4)), eye(4))
    @test isapprox(pauliliou2liou(eye(4)), eye(4))
end

@testset "depolarization" begin
    @test isapprox(depol(2),projector(vec(eye(2))),atol=1e-15)
    @test isapprox(depol(2),depol(2,1.))
    @test isapprox(depol(2,0.),eye(4),atol=1e-15)
end

@testset "density matrix properties" begin
  for i=1:100
    rrho = partialtrace(projector(randn(3^3)+1im*randn(3^3)),[3,9],2)
    @test isapprox(LinearAlgebra.tr(rrho),1.)
    @test QuantumInfo.ishermitian(rrho)
    @test ispossemidef(rrho)

    ru = LinearAlgebra.svd(randn(3,3)+1im*randn(3,3)).U
    @test isapprox(ru*ru',eye(3),atol=1e-13)
    @test isapprox(ru'*ru,eye(3),atol=1e-13)

    rv = LinearAlgebra.svd(randn(3,3)+1im*randn(3,3)).U
    re = (liou(ru)+liou(rv))/2
    @test ispossemidef(liou2choi(re),tol=1e-13)
    @test isapprox(partialtrace(liou2choi(re),[3,3],2),eye(3)/3,atol=1e-15)
  end
end


text:
module HnswAnn

using Printf

include("hnswrs.jl")

export
    Neighbour,
    Neighbourhood,
    HnswApi,
    createHnswApi,
    setRustlibPath,
    insert,
    search,
    fileDump,
    HnswDescription,
    getDescription,
    loadHnsw

"""

# HnswApi

This is the structure encapsulting the rust api.

## FIELDS
 -------

* rust : a reference pointer coming from rust library. An opaque structure not to be manipulated.

* type : The type of Vector

* maxNbConn : The number of connection to use in Hnsw structure. Between 8 and 64.

* efConstruction : The with of neighbours search used when constructing links between nodes.
    A rule of thumb is between maxNbConn and 64.

* distname : the name o the distance to use. It can be "DistL1", "DistL2", "DistCosine", DistDot",
        "DistHamming", "DistJaccard" or "DistPtr"

* A pointer to a function computing distances created by the macro @cfunction.
    The signature of the function must be Cfloat, (Ptr{Cfloat}, Ptr{Cfloat}, Culonglong)).

    See the testdistptr function in test.jl and the documentation of @cfunction.

"""
mutable struct HnswApi
    rust :: Ref{Hnswrs}
    type :: DataType
    maxNbConn::Int64
    efConstruction::Int64
    distname::String
    distfunctPtr :: Union{ Some{Ptr{Cvoid}} , Nothing }
end

"""
  #  `function createHnswApi(type::DataType, maxNbConn::Int64, efConstruction::Int64, distname::String)`

  This function creates a standard HnswApi with pre-defined distances.
  It benefits from parallel insert/search method

"""
function createHnswApi(type::DataType, maxNbConn::Int64, efConstruction::Int64, distname::String)
    rust_type_name = checkForImplementedType(type)
    rust = hnswInit(type, maxNbConn, efConstruction, distname)
    HnswApi(rust, type, maxNbConn, efConstruction, distname, nothing)
end

"""
# ` function createHnswApi(type::DataType, maxNbConn::Int64, efConstruction::Int64, distance::Ptr{Cvoid}) `

Create a HnswApi with a custom distance function pointer. In this case the name of the distance
is set to "DistPtr"

"""
function createHnswApi(type::DataType, maxNbConn::Int64, efConstruction::Int64, distance::Ptr{Cvoid})
    rust_type_name = checkForImplementedType(type)
    rust = hnswInit(type, maxNbConn, efConstruction, distance)
    # we store distance to avoid garbage collection of pointer as it is referenced in rust library.
    HnswApi(rust, type, maxNbConn, efConstruction, "DistPtr", Some(distance))
end

"""
 #  insert

 `function insert(hnsw::HnswApi, data::Vector{T}, id::UInt64) where {T <: Number}`

## Args

- data: the data vector to insert in struct hnsw
- id: the (unique) id of data. Can be the rank of insertion or any hash value enabling
        identification of point for possible dump/restore of the whole hnsw structure
"""
function insert(hnsw::HnswApi, data::Vector{T}, id::UInt64) where {T <: Number}
    one_insert(hnsw.rust, data, id)
end

"""
# parallel insertion

`function insert(hnsw::HnswApi, datas::Vector{Tuple{Vector{T}, UInt}}) where {T <: Number}`

This function does a block parallel insertion of datas in structure.

## Args
----
- datas: a vector of insertion request as a Tuple associating the point to insert and its id.
"""
function insert(hnsw::HnswApi, datas::Vector{Tuple{Vector{T}, UInt}}) where {T <: Number}
    parallel_insert(hnsw.rust, datas)
end

"""
# search

` function search(hnsw::Hnsw, vector::Vector{T}, knbn::Int64, ef_search ::Int64) where {T<:Number}`

 ## Args

- The vector for which we search neighbours
- knbn : the number of neughbour we search
- ef_search : the search parameter.

"""
function search(hnsw::HnswApi, vector::Vector{T}, knbn::Int64, ef_search ::Int64) where {T<:Number}
    one_search(hnsw.rust, vector, knbn, ef_search)
end

"""

# parallel search

`function search(hnsw::Hnsw , datas::Vector{Vector{T}}, knbn::Int64, ef_search:: Int64) where {T<:Number}`

As for insertion this function is interesting if we batch sufficiently many requests. (some hundreds)

##  Args

- datas : The vector of point we search the neighbours of.
- knbn the number of neighbours wanted.
- ef_search : parameter search.
"""
function search(hnsw::HnswApi , datas::Vector{Vector{T}}, knbn::Int64, ef_search:: Int64) where {T<:Number}
    parallel_search(hnsw.rust, datas, knbn, ef_search)
end

"""
# `function fileDump(hnsw::HnswApi,  filename::String)`

Returns 1 if OK , -1 if failure.
`function fileDump(hnsw::HnswApi,  filename::String)`

This dumps the whole data in the structure in 2 binary files:
- one file has name obtained by appending ".hnsw.graph" to filename and contains the graph structure of Hnsw
- the other file has name obtained by appending ".hnsw.data" to filename and contains for each data vector
    inserted in the hnsw structure the couple id , data vector

* NOTE:
    It is not possible to dump/reload a HnswApi with a custom distance generated by
    the macro @cfunction in the current version.
"""
function fileDump(hnsw::HnswApi,  filename::String)
    res = filedump(hnsw.rust, hnsw.type, filename)
    res
end

end # module


text:
using ForwardDiff
using LinearAlgebra
include("backtracking_line_search.jl")
include("newton_method.jl")
include("gradient_algorithm_l2.jl")
include("coordinate_algorithm_l1.jl")
include("gradient_algorithm_l2.jl")
include("coordinate_algorithm_l1.jl")
dimension = 10 #defult
epsilon = 0.001 #default
objective(x) =  x' * x   #default
x = rand(dimension) #default, point selection is also an important step
alpha = 0.1 #default
beta = 0.5 #default
step = 100 #default

x = newton_algo(epsilon,objective,x)
x = coordinate_algo(dimension,epsilon,objective,x,alpha,beta,step)

print(x)
print(objective(x))

#brief readme: rough underestimator. for learning only

max_stars_repo_name: zhuhs-q/basic-optimal-method

max_issues_repo_name: Jacob-jinxy/basic-optimal-method

max_forks_repo_name: Jacob-jinxy/basic-optimal-method


text:
module JlrsTests
mutable struct MutableStruct
    x
    y::UInt64
end

function inlinetuple()::Tuple{UInt32, UInt16, Int64}
    (1, 2, 3)
end

function symbol()::Symbol
    :x
end

function base()::Module
    Base
end

function datatype()::DataType
    Bool
end

function callrust(ptr::Ptr)::Bool
    ccall(ptr, Bool, ())
end

function callrustwitharr(ptr::Ptr, arr::Array{Float64, 1})::Bool
    ccall(ptr, Bool, (Array,), arr)
end

function vecofmodules()::Vector{Module}
    [Base; Core; Main]
end

function anothervecofmodules()::Vector{Module}
    [Base; Core; Main]
end
end

max_stars_repo_name: nikhilmitrax/jlrs

max_issues_repo_name: nikhilmitrax/jlrs

max_forks_repo_name: nikhilmitrax/jlrs


text:
using BioSymbols
using BioSequences
using ArgParse
using Glob

using Rifraf

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--multi-reference"
        help = "each sequence is followed by its reference"
        action = :store_true

        "--log-p"
        help = "log error probability"
        arg_type = Float64
        default = -1.0

        "--bandwidth"
        help = "alignment bandwidth. If < 0, choose dynamically."
        arg_type = Int
        default = -1

        "--errors"
        help = "comma-seperated reference error ratios - mm, codon ins, codon del"
        arg_type = String
        default = "10,0.00001,0.00001,1,1"

        "--verbose", "-v"
        help = "print progress"
        arg_type = Int
        default = 0

        # required arguments
        "input"
        help = "input fasta file, with sequence/reference alternating pairs"
        required = true

        "output"
        help = "output fasta file of corrected sequences"
        required = true
    end
    return parse_args(s)
end

function main()
    args = parse_commandline()

    score_args = map(x -> parse(Float64, x), split(args["errors"], ","))
    scores = Scores(ErrorModel(score_args...))

    input = args["input"]
    input_sequences = Rifraf.read_fasta_records(input)

    multi = args["multi-reference"]
    sequences = multi ? input_sequences[1:2:end-1] : input_sequences[2:end]
    references = multi ? input_sequences[2:2:end] : repeated(input_sequences[1])

    stream = open(FASTAWriter, args["output"])
    for (seq, ref) in zip(sequences, references)
        result = Rifraf.correct_shifts(DNASeq(seq.seq), DNASeq(ref.seq);
                                       log_p=args["log-p"],
                                       bandwidth=args["bandwidth"],
                                       scores=scores)
        record = Seq.FASTASeqRecord(seq.name, result)
        write(stream, record)
    end
    close(stream)
end

main()


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

# This is a toolbox for integrate general grid to gaussian grid using ESMF_RegridWeiGen
# The transformation refers to
# [Section 12.9.3] of http://www.earthsystemmodeling.org/esmf_releases/public/ESMF_8_0_0/ESMF_refdoc/node3.html

mutable struct ESMFWeight

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
module MTIWrapper

"""
    function install_web_api(root_dir)
Installs MTI Java jar to specified location
"""
function install_web_api(root_dir; api_version="SKR_Web_API_V2_3")

    module_dir = dirname(@__FILE__)
    curr_dir = Base.pwd()

    @info "MTIWrapper: Current Directory $module_dir"

    api_dir = "$root_dir/$api_version"

    if !isdir(root_dir)
        mkdir(root_dir)
    else
        if isdir(api_dir)
           rm(api_dir, recursive = true, force = true)
           rm("$api_dir.jar")
        end
    end

    println("Dowloading fresh copy of sources")
    # get sources and expand
    cd(root_dir)
    download("https://ii.nlm.nih.gov/Web_API/$api_version.jar", "$api_version.jar")
    run(`java sun.tools.jar.Main xf $api_version.jar`)

    # We have modified our class to pass credentials as an input - we'll put it in examples folder
    # so that MTI default compilation scripts work
    cp("$module_dir/GenericBatchCustom.java", "$api_dir/examples/GenericBatchCustom.java", force=true)
    cp("$module_dir/generic_batch.sh", "$root_dir/generic_batch.sh", force=true)

    # compile
    # cd("SKR_Web_API_V2_3")
    run(`chmod +x $api_dir/compile.sh $api_dir/run.sh $api_dir/build.sh`)
    run(`$api_dir/compile.sh $api_dir/examples/GenericBatchCustom.java`)

    cd(curr_dir)
    println(pwd())

end

"""
    function mti_batch_query(mti_java_dir, email, username, password, in_file, out_file)
Send a batch query to MTI. Use function `abstracts_to_request_file` to generate `in_file`
"""
function mti_batch_query(mti_java_dir, email, username, password, in_file, out_file)

    run(`$mti_java_dir/generic_batch.sh $mti_java_dir $email $username $password $in_file $out_file`)
end

"""
    abs_to_request_file()

Write all abstracts in a dataframe, to a file to be used for MTI batch query.
The format is:

UI - pmid
AB - abstract_text
"""
function abstracts_to_request_file(df, out_file; uid_column::Symbol = :pmid, abstract_column::Symbol = :abstract)

    #call MTI
    open(out_file, "w") do file

        for row in eachrow(df)
            uid = row[uid_column]
            abstract_text = row[abstract_column]

            if ismissing(abstract_text)
                println( "Skipping empty abstract for uid: ", uid)
                continue
            end

            # convert to ascii - all unicode caracters to " "
            abstract_ascii = replace(abstract_text, r"[^\u0000-\u007F]"=>"")
            write(file, "UI  - $uid \n")
            write(file, "AB  - $abstract_ascii \n \n")
        end
    end

end

"""
    abs_to_request_file()

Write all abstracts in a dataframe, to a file to be used for MTI batch query.
The format is:

UI - pmid
AB - abstract_text
"""
function title_abstracts_to_request_file(df, out_file; uid_column::Symbol = :pmid,
                                         title_column::Symbol = :title, abstract_column::Symbol = :abstract)

    #call MTI
    open(out_file, "w") do file

        for row in eachrow(df)
            uid = row[uid_column]
            abstract_text = row[abstract_column]
            title_text = row[title_column]

            write(file, "UI  - $uid \n")

            if !ismissing(title_text)
                # convert to ascii - all unicode caracters to " "
                title_ascii = replace(title_text, r"[^\x00-\x7F]+"=>"")
                write(file, "TI  - $title_ascii \n")
            end

           if !ismissing(abstract_text)
                # convert to ascii - all unicode caracters to " "
                abstract_ascii = replace(abstract_text, r"[^\x00-\x7F]+"=>"")
                abstract_ascii = replace(abstract_ascii, "Abstract: "=>"")
                write(file, "AB  - $abstract_ascii \n \n")
            end

        end
    end

end

end # module


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.6.0
#   kernelspec:
#     display_name: Julia 1.5.1
#     language: julia
#     name: julia-1.5
# ---

using Plots

x = range(-pi,pi,length=100)
y = sin.(x)
plot(aspect_ratio=:equal,xlim=[-pi,pi])
plot!(sin,color=:red)
plot!(x, y, fill=y .+ 1,color=:green)
plot!(x, y, fill=y .- 1,color=:blue)

x = range(-pi,pi,length=100)
y = sin.(x)
l = -ones(length(y))
u = ones(length(y))
plot(aspect_ratio=:equal,xlim=[-pi,pi])
plot!(sin, ribbons=[l, u], color=:blue,label=:none)
plot!(sin,color=:red,label=:none)

x=range(-pi,pi,length=200)
plot(x,sin.(x),fill=(x-> sin(x)>0 ? sin(x) : 0),alpha=0.5)


text:
struct LTCCell{W,SENS,SYN,SENSP,SENSRE,SYNP,SYNRE,MI,MO,V,S,SOLVER,SENSEALG}

  wiring::W

  mapin::MI
  mapout::MO

  sens_f::SENS
  sens_p::SENSP
  sens_pl::Int
  sens_re::SENSRE
  syn_f::SYN
  syn_p::SYNP
  syn_pl::Int
  syn_re::SYNRE

  cm::V
  Gleak::V
  Eleak::V

  state0::S

  solver::SOLVER
  sensealg::SENSEALG

  function LTCCell(wiring, mapin, mapout, sens_f, sens_p, sens_pl, sens_re, syn_f, syn_p, syn_pl, syn_re, cm, Gleak, Eleak, state0, solver, sensealg)
    new{typeof(wiring),typeof(sens_f),typeof(syn_f),typeof(sens_p),typeof(sens_re),typeof(syn_p),typeof(syn_re),typeof(mapin),typeof(mapout),typeof(cm),typeof(state0),typeof(solver),typeof(sensealg)}(
      wiring,mapin,mapout,sens_f,sens_p,sens_pl,sens_re,syn_f,syn_p,syn_pl,syn_re,cm,Gleak,Eleak,state0,solver,sensealg)
  end
end

function LTCCell(wiring, solver, sensealg)
  n_in = wiring.n_in
  out = wiring.out
  n_total = wiring.n_total
  cm    = rand_uniform(1.6, 2.5, n_total)
  Gleak = rand_uniform(0.001, 1, n_total)
  Eleak = rand_uniform(-0.2, 0.2, n_total)
  state0 = fill(Float32(0.01), n_total, 1)

  mapin = Mapper(n_in)
  mapout = Mapper(n_total)

  sens_f = gNN(n_in, n_total)
  syn_f = gNN(n_total, n_total)
  sens_p, sens_re = Flux.destructure(sens_f)
  syn_p, syn_re = Flux.destructure(syn_f)

  LTCCell(wiring,mapin,mapout,sens_f,sens_p,length(sens_p),sens_re,syn_f,syn_p,length(syn_p),syn_re, cm,Gleak,Eleak,state0,solver,sensealg)
end

Flux.@functor LTCCell (mapin, mapout, sens_p, syn_p, cm, Gleak, Eleak, state0,)
Flux.trainable(m::LTCCell)=(Flux.trainable(m.mapin), Flux.trainable(m.mapout), m.sens_p, m.syn_p, m.cm, m.Gleak, m.Eleak, m.state0)

NCP(args...) = Flux.Recur(LTCCell(args...))
Flux.Recur(m::LTCCell) = Flux.Recur(m, m.state0)

#(m::LTCCell)(h::AbstractVector, x::AbstractVecOrMat) == m(repeat(h, 1, size(x,2)-size(h,2)+1), x)

function (m::LTCCell)(h::AbstractVecOrMat, x::AbstractVecOrMat)

  x = m.mapin(x)
  I = x#::Matrix{Float32}

  # TODO: sens_f() here?

  u0 = repeat(h, 1, size(x,2)-size(h,2)+1)

  # condition(u,t,integrator) = any(x -> (x < -100) || (x > 100), u)
  # affect!(integrator) = terminate!(integrator)
  # cb = DiscreteCallback(condition,affect!)

  # sens_f = m.sens_re(m.sens_p)
  # syn_f  = m.syn_re(m.syn_p)
  # p = [m.cm; m.Gleak; m.Eleak; m.sens_p; m.syn_p]::Vector{Float32}
  prob = ODEProblem{true}((dx,x,p,t)->dltcdt!(dx,x,p,t, m, I),u0,Float32.((0, 1)),[m.cm; m.Gleak; m.Eleak; m.sens_p; m.syn_p])
  #prob = ODEProblem((x,p,t)->oop(x,p,t, m, m.sens_f, m.syn_f, I),u0,Float32.((0, 1)),[m.cm; m.Gleak; m.Eleak; m.sens_p; m.syn_p])
  #prob = ODEProblem((x,p,t)->oop(x,p,t, m),u0,Float32.((0, 1)),[m.cm; m.Gleak; m.Eleak; m.sens_p; m.syn_p; x[:]])
  sol = solve(prob, m.solver; dt=0.1, adaptive=false, sensealg=m.sensealg, save_everystep=false, save_start=false)#, isoutofdomain=(u,p,t)->any(x -> (x < -1e40) || (x > 1e40), u))
  #solve(prob, m.solver, alias_u0=true, dtmax=1.0, sensealg=m.sensealg, save_everystep=false, save_start=false)#, isoutofdomain=(u,p,t)->any(x -> (x < -10) || (x > 10), u))
  #sol = solve(prob, m.solver; sensealg=m.sensealg, save_everystep=false, save_start=false)

  #hh = copy(u0)#[:,:]
  hh = sol[:,:,end]#::typeof(h)
  #h =
  out = m.mapout(hh)
  return hh, @view out[end-m.wiring.n_motor+1:end, :]
end

#function dltcdt!(dx,x,p,t, m, I::Matrix{Float32}, sens_f, syn_f)
function dltcdt!(dx,x,p,t, m, I)

  cm_pl = size(x,1)::Int
  Gleak_pl = size(x,1)::Int
  Eleak_pl = size(x,1)::Int
  sens_pl = m.sens_pl::Int
  syn_pl = m.syn_pl::Int

  cm    = @view p[1 : cm_pl]#::AbstractVecOrMat{eltype(x)}
  Gleak = @view p[cm_pl + 1 : cm_pl + Gleak_pl]#::AbstractVecOrMat{eltype(x)}
  Eleak = @view p[cm_pl + Gleak_pl + 1 : cm_pl + Gleak_pl + Eleak_pl]#::AbstractVecOrMat{eltype(x)}

  sens_p = @view p[cm_pl + Gleak_pl + Eleak_pl + 1 : cm_pl + Gleak_pl + Eleak_pl + sens_pl]#::AbstractVecOrMat{eltype(x)}
  syn_p = @view p[cm_pl + Gleak_pl + Eleak_pl + sens_pl + 1 : cm_pl + Gleak_pl + Eleak_pl + sens_pl + syn_pl]#::AbstractVecOrMat{eltype(x)}

  # I = reshape((@view p[cm_pl + Gleak_pl + Eleak_pl + sens_pl + syn_pl + 1 : end]), m.wiring.n_in, :)

  sens_f = m.sens_re(sens_p)
  syn_f  = m.syn_re(syn_p)

  #@unpack cm, Gleak, Eleak, sens_p, syn_p, I = p

  I_sens = sens_f(x, I, m.wiring.sens_mask, m.wiring.sens_pol)#::AbstractVecOrMat{eltype(x)}
  I_syn = syn_f(x, x, m.wiring.syn_mask, m.wiring.syn_pol)#::AbstractVecOrMat{eltype(x)}

  @. dx = cm * (-(Gleak * (x - Eleak)) + I_syn + I_sens)
  nothing
end

# function oop(x,p,t, m, I::Matrix{Float32})
function oop(x,p,t, m, sens_f, syn_f, I)

  cm_pl = size(x,1)::Int
  Gleak_pl = size(x,1)::Int
  Eleak_pl = size(x,1)::Int
  sens_pl = m.sens_pl::Int
  syn_pl = m.syn_pl::Int

  cm    = @view p[1 : cm_pl]#::AbstractVecOrMat{eltype(x)}
  Gleak = @view p[cm_pl + 1 : cm_pl + Gleak_pl]#::AbstractVecOrMat{eltype(x)}
  Eleak = @view p[cm_pl + Gleak_pl + 1 : cm_pl + Gleak_pl + Eleak_pl]#::AbstractVecOrMat{eltype(x)}

  sens_p = @view p[cm_pl + Gleak_pl + Eleak_pl + 1 : cm_pl + Gleak_pl + Eleak_pl + sens_pl]#::AbstractVecOrMat{eltype(x)}
  syn_p = @view p[cm_pl + Gleak_pl + Eleak_pl + sens_pl + 1 : cm_pl + Gleak_pl + Eleak_pl + sens_pl + syn_pl]#::AbstractVecOrMat{eltype(x)}

  # I = reshape(p[cm_pl + Gleak_pl + Eleak_pl + sens_pl + syn_pl + 1 : end], m.wiring.n_in, :)

  sens_f = m.sens_re(sens_p)
  syn_f  = m.syn_re(syn_p)
  #@unpack cm, Gleak, Eleak, sens_p, syn_p, I = p

  I_sens = sens_f(x, I, m.wiring.sens_mask, m.wiring.sens_pol)#::AbstractVecOrMat{eltype(x)}
  I_syn = syn_f(x, x, m.wiring.syn_mask, m.wiring.syn_pol)#::AbstractVecOrMat{eltype(x)}

  @. cm * (-(Gleak * (x - Eleak)) + I_syn + I_sens)
end

struct Mapper{V<:AbstractArray}
  W::V
  b::V
  Mapper(W::V,b::V) where {V<:AbstractArray} = new{V}(W,b)
end
Mapper(in::Integer) = Mapper(ones(Float32,in), zeros(Float32,in))

(m::Mapper)(x::AbstractVecOrMat) = m.W .* x .+ m.b

Flux.@functor Mapper

struct gNN{A<:AbstractArray}
  G::A
  μ::A
  σ::A
  E::A
  gNN(G::A, μ::A, σ::A, E::A) where {A<:AbstractArray} = new{A}(G, μ, σ, E)
end
function gNN(in::Integer, out::Integer)
  G = rand_uniform(0.001, 1, in, out)
  μ = rand_uniform(0.3, 0.8, in, out)
  σ = rand_uniform(3, 8, in, out)
  E = rand_uniform(-0.3, 0.3, in, out)
  gNN(G,μ,σ,E)
end

function (m::gNN)(h::AbstractVecOrMat, x::AbstractVecOrMat, bitmask::Matrix{Float32}, polmask::Matrix{Float32})
  G, μ, sig, E = m.G, m.μ, m.σ, m.E

  #fun4e(bitmask, G, E, x, polmask, h, μ, σ) = @tullio out[j] := bitmask[i,j] * G[i,j] * sigmoid((x[i] - μ[i,j]) * σ[i,j]) * (polmask[i,j] * E[i,j] - h[j])

  #fun1(bitmask, G, E, x, polmask, h, μ, σ) = reshape(reduce(vcat, [reshape(sum(bitmask .* G .* Flux.sigmoid.(((x[:,b]) .- μ) .* σ) .* (polmask .* E .- reshape((h[:,b]),1,:)),dims=1), :) for b in 1:size(h,2)]), size(h,1),:)
  function fun1(bitmask, G, E, x, polmask, h, μ, σ)
    out = reduce(hcat, [reshape(sum(bitmask .* G .* sigmoid.((x[:,b] .- μ) .* σ) .* ( polmask .* E .- reshape(h[:,b], 1,:)),dims=1),:) for b in 1:size(x,2)])
    return out
  end

  function fun2(bitmask, G, E, x, polmask, h, μ, σ)
    @views out = hcat([reshape(sum(bitmask .* G .* sigmoid.((x[:,b] .- μ) .* σ) .* (polmask .* E .- reshape(h[:,b], 1,:)),dims=1),:) for b in 1:size(x,2)]...)
    return out
  end

  # First function call produced NaNs. Exiting.
  function fun3(bitmask, G, E, x, polmask, h, μ, σ)
  	@views out = hcat([reshape(sum(bitmask .* G .* (1 ./ (1 .+ exp.(-((x[:,b] .- μ) .* σ)))) .* (polmask .* E .- reshape(h[:,b], 1,:)),dims=1),:) for b in 1:size(x,2)]...)
  	return out
  end

  # First function call produced NaNs. Exiting.
  # function fun4(bitmask, G, E, x, polmask, h, μ, σ)
  #   @tullio out[j,b] := bitmask[i,j] * G[i,j] * (1 / (1 + exp(-((x[i,b] - μ[i,j]) * σ[i,j])))) * (polmask[i,j] * E[i,j] - h[j,b]) #threads=false avx=false verbose=true
  #   return out
  # end
  # no gradient
  function fun5(bitmask, G, E, x, polmask, h, μ, sig)
    @tullio et[i,j,b] := (x[i,b] - μ[i,j]) * sig[i,j] verbose=true
    @tullio t[i,j,b] := exp(-abs(et[i,j,b])) verbose=true
    @tullio s[i,j,b] := ifelse(et[i,j,b] ≥ 0, inv(1 + t[i,j,b]), t[i,j,b] / (1 + t[i,j,b])) verbose=true

    @tullio out[j,b] := bitmask[i,j] * G[i,j] * s[i,j,b] * (polmask[i,j] * E[i,j] - h[j,b]) verbose=true #threads=false avx=false
    return out
  end
  function fun6(bitmask, G, E, x, polmask, h, μ, sig)
    @tullio out[j,b] := bitmask[i,j] * G[i,j] * sigmoid(-((x[i,b] - μ[i,j]) * sig[i,j])) * (polmask[i,j] * E[i,j] - h[j,b]) grad=Dual #threads=false avx=false verbose=true
    return out
  end

  return fun1(bitmask, G, E, x, polmask, h, μ, sig)

  # @views reshape(reduce(vcat, [reshape(sum(bitmask .* G .* Flux.sigmoid.(((x[:,b]) .- μ) .* sig) .* (polmask .* E .- reshape((h[:,b]),1,:)),dims=1), :) for b in 1:size(h,2)]), size(h,1),:)
  # Flux.stack([reshape(sum(bitmask .* G .* Flux.sigmoid.(((x[:,b]) .- μ) .* sig) .* (polmask .* E .- reshape((h[:,b]),1,:)),dims=1), :) for b in 1:size(h,2)], 2)
  #reshape(reduce(vcat, [reshape(sum(bitmask .* G .* Flux.sigmoid.(((@view x[:,b]) .- μ) .* σ) .* (polmask .* E .- reshape((@view h[:,b]),1,:)),dims=1), :) for b in 1:size(h,2)]), size(h,1),:)

  #Flux.stack([reshape(sum(bitmask .* G .* Flux.sigmoid.(((@view x[:,b]) .- μ) .* σ) .* (polmask .* E .- reshape((@view h[:,b]),1,:)),dims=1), :) for b in 1:size(h,2)], 2)
end

function (m::gNN)(p::AbstractVector, h::AbstractVecOrMat, x::AbstractVecOrMat, bitmask::Matrix{Float32}, polmask::Matrix{Float32})
  Gl = length(m.G)
  μl = length(m.μ)
  σl = length(m.σ)
  El = length(m.E)
  @views G = p[1 : Gl]
  @views μ = p[Gl + 1 : Gl + μl]
  @views σ = p[Gl + μl + 1 : Gl + μl + σl]
  @views E = p[Gl + μl + σl + 1 : Gl + μl + σl + El]

  G = reshape(G, size(m.G))
  μ = reshape(μ, size(m.μ))
  σ = reshape(σ, size(m.σ))
  E = reshape(E, size(m.E))

  #fun4e(bitmask, G, E, x, polmask, h, μ, σ) = @tullio out[j] := bitmask[i,j] * G[i,j] * sigmoid((x[i] - μ[i,j]) * σ[i,j]) * (polmask[i,j] * E[i,j] - h[j])

  #fun1(bitmask, G, E, x, polmask, h, μ, σ) = reshape(reduce(vcat, [reshape(sum(bitmask .* G .* Flux.sigmoid.(((x[:,b]) .- μ) .* σ) .* (polmask .* E .- reshape((h[:,b]),1,:)),dims=1), :) for b in 1:size(h,2)]), size(h,1),:)
  function fun1(bitmask, G, E, x, polmask, h, μ, σ)
    @views out = reduce(hcat, [reshape(sum(bitmask .* G .* sigmoid.((x[:,b] .- μ) .* σ) .* ( polmask .* E .- reshape(h[:,b], 1,:)),dims=1),:) for b in 1:size(x,2)])
    return out
  end

  return fun1(bitmask, G, E, x, polmask, h, μ, σ)

end

Flux.@functor gNN


text:
module AlertPushover

export pushover_alert!

using Alert
using HTTP

"""
    pushover_alert!(;token,user)

Use the [Pushover](https://pushover.net) notificaiton service to send notifications when
calling `alert`. You must provide an API and user token, which are documented in the
[description of pushover's api](https://pushover.net/api).

Call this funciton with no arguments to return to the default alert notification backend.
"""
function pushover_alert!(;token=nothing,user=nothing)
    if isnothing(token) && isnothing(user)
        Alert.set_alert_backend!()
    end

    Alert.set_alert_backend!() do message
        @info "Sending message '"*message*"' to https://pushover.net."
        HTTP.post("https://api.pushover.net/1/messages.json",[],
            readtimeout = 15,
            retry_non_indempotent = true,
            retry = true,
            retries = 8,
            HTTP.URIs.escapeuri(Dict(
                :token => string(token),
                :user => string(user),
                :message => string(message)
            ))
        )

        nothing
    end
    @info "Using remote pushover service as `alert` backend."

    nothing
end

end


text:
module DCEGM

using JSON
using Interpolations
using Distributions: Normal
using FastGaussQuadrature
using Roots
using Plots
using Colors
using Statistics
using Random
using Interact
using DataStructures: OrderedDict
using InvertedIndices: Not
using Blink
using BenchmarkTools

import Base: +, -, *, /, ==, isless, in
import Base: iterate, typemin, isapprox, convert
import Base.size, Base.getindex, Base.setindex!, Base.eltype
import Base.prepend!,
       Base.append!,
       Base.insert!,
       Base.delete!,
       Base.sort!,
       Base.length,
       Base.show,
       Base.IndexStyle

# Types
export MLine, Point, Envelope, yisless

# methods
export interp, splitat, getx, gety, gets, splitLine, getr

# includes for MultiLine
include("line.jl")
include("envelope.jl")

# # includes for dcegm
include("param.jl")
include("utils.jl")
include("minimal_EGM.jl")
include("dc_algo.jl")
include("simulate.jl")

include("plotting.jl")
include("bench.jl")
include("bk.jl")
include("interact.jl")

end # module


text:
function iterated_local_search(tuning_run::Run,
                               channel::RemoteChannel;
                               acceptance_criterion::Function = metropolis,
                               temperature::Function          = log_temperature,
                               subsidiary_iterations::Int     = 30,
                               threshold::AbstractFloat       = 2.)
    iteration  = 1

    function iterate(tuning_run::Run)
        for i = 1:subsidiary_iterations
            p      = temperature(i)
            result = probabilistic_improvement(tuning_run,
                                               threshold = p)
        end

        p          = temperature(i)
        iteration += 1
        return probabilistic_improvement(tuning_run,
                                         threshold = p)
    end

    technique(tuning_run,
              channel,
              iterate,
              name = "Iterated Local Search")
end


text:
# Shared methods between surface meshers
# * pushback
# * isoriented
# * volume (because isoriented could be made more robust)
# * subdivision

"Calculation of precise volume. Needed for convergence check"
function volume(points,faces)

    normal0 = [0,0,1]

    s = 0

    for tri in 1:size(faces,2)
        face = faces[:,tri]
        y1 = points[:,face[1]]
        y2 = points[:,face[2]]
        y3 = points[:,face[3]]

        normaly = cross(y2-y1,y3-y1)
        normaly /= norm(normaly)

        area = norm(cross(y2-y1,y3-y1))/2
        areaproj = dot(normaly,normal0)*area
        volume = dot(y1 + y2 + y3,normal0)/3*areaproj

        s += volume
    end

    return s

    # Calculate face normal
    # Calculate projected area
    # Calculate ordinary volume
    # Calculate volume between projected and real area
    # (+) if normal is outwards
end

using ForwardDiff               # Can be made optional

function pushback(sdist::Function,x::Array{Float64,1})
    for i in 1:100
        grad = ForwardDiff.gradient(sdist,x)
        x = x - sdist(x)*grad/norm(grad)^2
        # I allow the root to be only on gradient
    end
    return x
end

### Loading mesh by filename in the /preproc directgory
### At this point it should be unimportant if it is either generated with cgal or matlab

### Creation routines

function isoriented(points::Array{T,2},triangles::Array{S,2}) where {T <: AbstractFloat, S <: Integer}

    tri = triangles[:,1]
    p1 = points[:,tri[1]]
    p2 = points[:,tri[2]]
    p3 = points[:,tri[3]]

    a1 = p2 - p1
    a2 = p3 - p1

    normal = cross(a1,a2)
    normal /= sqrt(dot(normal,normal))

    direction = p1
    direction /= sqrt(dot(p1,p1))

    angle = acos(dot(normal,direction))

    if angle<pi/2
        return true
    else
        return false
    end
end

function subdivision(faces)
    edges = []

    for ti in 1:size(faces,2)
        v1,v2,v3 = faces[:,ti]

        v2<v1 || push!(edges,(v1,v2))
        v3<v2 || push!(edges,(v2,v3))
        v1<v3 || push!(edges,(v3,v1))
    end

    rfaces = Array(Int,3,size(faces,2)*4)
    N = maximum(faces)

    for ti in 1:size(faces,2)

        v1,v2,v3 = faces[:,ti]
        e3 = findfirst(edges,v2>v1 ? (v1,v2) : (v2,v1)) + N
        e1 = findfirst(edges,v3>v2 ? (v2,v3) : (v3,v2)) + N
        e2 = findfirst(edges,v1>v3 ? (v3,v1) : (v1,v3)) + N

        rfaces[:,4*ti - 3] = [v1,e3,e2]
        rfaces[:,4*ti - 2] = [v2,e1,e3]
        rfaces[:,4*ti - 1] = [v3,e2,e1]
        rfaces[:,4*ti] = [e1,e2,e3]
    end

    return rfaces
end

function subdivision(points,faces; n=1,method=:linear)
    if method==:linear
        rpoints,rfaces = LinearSubdivision(points,faces)
    elseif method==:paraboloid
        rpoints,rfaces = ParaboloidSubdivision(points,faces)
    end

    return rpoints,rfaces
end

function LinearSubdivision(points,faces)

    rfaces = subdivision(faces)
    rpoints = Array(Float64,3,maximum(rfaces))

    for ti in 1:size(faces,2)

        v1,v2,v3 = faces[:,ti]
        e1,e2,e3 = rfaces[:,4*ti]

        rpoints[:,v1] = points[:,v1]
        rpoints[:,v2] = points[:,v2]
        rpoints[:,v3] = points[:,v3]
        rpoints[:,e1] = (points[:,v2] + points[:,v3])/2
        rpoints[:,e2] = (points[:,v1] + points[:,v3])/2
        rpoints[:,e3] = (points[:,v1] + points[:,v2])/2
    end

    return rpoints,rfaces
end

function ParaboloidSubdivision(points,faces)

    normals = Array(Float64,size(points)...)
    NormalVectors!(normals,points,faces,i->FaceVRing(i,faces))

    rfaces = subdivision(faces)
    rpoints = Array(Float64,3,maximum(rfaces))

    vproperties = []
    for xkey in 1:size(points,2)

        x = points[:,xkey]
        nx = normals[:,xkey]

        vects = Float64[]
        for i in VertexVRing(xkey,faces)
            vi = points[:,i] - x
            vects = [vects..., vi...]
        end
        vects = reshape(vects,3,div(length(vects),3))

        R = eye(3)
        C,D,E = SurfaceGeometry.ZinchenkoDifertential!(vects,nx,R)

        push!(vproperties,(C,D,E,R))
    end

    z(x,y,C,D,E) = C*x^2 + D*x*y + E*y^2
    pushon(x,C,D,E,R) = ((xp,yp,zp) = R*x; inv(R)*[xp,yp,z(xp,yp,C,D,E)])

    for ti in 1:size(faces,2)
        v1,v2,v3 = faces[:,ti]
        e1,e2,e3 = rfaces[:,4*ti]

        rpoints[:,v1] = points[:,v1]
        rpoints[:,v2] = points[:,v2]
        rpoints[:,v3] = points[:,v3]

        ### It gets overwritten once
        xe1 = (points[:,v2] + points[:,v3])/2
        rpoints[:,e1] = xe1 + (pushon(xe1-points[:,v2],vproperties[v2]...) + pushon(xe1-points[:,v3],vproperties[v3]...))/2

        xe2 = (points[:,v3] + points[:,v1])/2
        rpoints[:,e2] = xe2 + (pushon(xe2-points[:,v1],vproperties[v1]...) + pushon(xe2-points[:,v3],vproperties[v3]...))/2

        xe3 = (points[:,v1] + points[:,v2])/2
        rpoints[:,e3] = xe3 + (pushon(xe3-points[:,v1],vproperties[v1]...) + pushon(xe3-points[:,v2],vproperties[v2]...))/2
    end

    return rpoints,rfaces
end

"Subdivision of surface and pushback"
function subdivision(points,faces,sdist::Function)
    rpoints,rfaces = LinearSubdivision(points,faces)

    for i in 1:size(rpoints,2)
        rpoints[:,i] = pushback(sdist,rpoints[:,i])
    end

    return rpoints,rfaces
end


text:
# ---
# title: 1491. Average Salary Excluding the Minimum and Maximum Salary
# id: problem1491
# author: zhwang
# date: 2022-03-18
# difficulty: Easy
# categories: Array, Sort
# link: <https://leetcode.com/problems/average-salary-excluding-the-minimum-and-maximum-salary/description/>
# hidden: true
# ---
#
# Given an array of **unique** integers `salary` where `salary[i]` is the salary
# of the employee `i`.
#
# Return the average salary of employees excluding the minimum and maximum
# salary.
#
#
#
# **Example 1:**
#
#
#
#     Input: salary = [4000,3000,1000,2000]
#     Output: 2500.00000
#     Explanation: Minimum salary and maximum salary are 1000 and 4000 respectively.
#     Average salary excluding minimum and maximum salary is (2000+3000)/2= 2500
#
#
# **Example 2:**
#
#
#
#     Input: salary = [1000,2000,3000]
#     Output: 2000.00000
#     Explanation: Minimum salary and maximum salary are 1000 and 3000 respectively.
#     Average salary excluding minimum and maximum salary is (2000)/1= 2000
#
#
# **Example 3:**
#
#
#
#     Input: salary = [6000,5000,4000,3000,2000,1000]
#     Output: 3500.00000
#
#
# **Example 4:**
#
#
#
#     Input: salary = [8000,9000,2000,3000,6000,1000]
#     Output: 4750.00000
#
#
#
#
# **Constraints:**
#
#   * `3 <= salary.length <= 100`
#   * `10^3 <= salary[i] <= 10^6`
#   * `salary[i]` is unique.
#   * Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

function average_salary_excluding_mami(salary::Vector{Int})
    return (sum(salary) - minimum(salary) - maximum(salary)) / (length(salary) - 2)
end
## @lc code=end


text:
"""
    direct time marching refers essentially to explicit Euler integration.
"""
function direct_time_marching!(Δt,U,Z,L)
    dU = (I - U*U')*L*Z*pinv(Z'*Z)
    dZ = L'*U
    Z .+= Δt * dZ
    U .+= Δt * dU
end

function direct_time_marching!(Δt,U,Z,L_U,L_Z)
    dU = (I - U*U')*L_U*(L_Z'*Z*pinv(Z'*Z))
    dZ = L_Z*(L_U'*U)
    U .+= Δt * dU
    Z .+= Δt * dZ
end

"""
    explicit time step with subsequent projection into low rank manifold
"""
step_and_project!(Δt, U, Z, L) = project_to_low_rank!(U, Z, U*Z' + Δt*L)
step_and_project!(Δt, U, Z, L_U, L_Z) = project_to_low_rank!(U, Z, cat(U, L_U, dims=2), cat(Z, Δt*L_Z, dims=2))

"""
    integration along geodesics of the low rank matrix manifold
    --> High order integration schemes
    --> Direct time marching approximates this only to first order
"""
function exponential_map_integration!(Δt, U, Z, dU, dZ)
# TBD
# requires additional implementation of numerical schemes to solve the
# geodesic equations (ODE):
# d²Z/dt² - Z (dU/dt)ᵀ dU/dt = 0
# d²U/dt² + U (dU/dt)ᵀ dU/dt + 2 dU/dt (dZ/dt)ᵀ Z inv(ZᵀZ) = 0
# with IC: dU/dt(0) = dU, U(0) = U
#          dZ/dt(0) = dZ, Z(0) = Z
# The ODE system is solved for [0, Δt] and the update is given by U(Δt), Z(Δt)
# positive: We can write tailored code to integrate these equations
end


text:
using Graphs

function create_graph(start_node, end_node)
    @assert Base.length(start_node)==Base.length(end_node)

    no_node = max(maximum(start_node), maximum(end_node))
    no_arc = Base.length(start_node)

    graph = simple_inclist(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin
        parent = state.parents[current]
        link_idx = link_dic[parent,current]
        if link_idx != 0
            x[link_idx] = 1
        end
        current = parent
    end
    return x
end

# preparing a graph
graph = create_graph(start_node, end_node)
link_dic = sparse(start_node, end_node, 1:numLinks)

function BPR(flowVec, fcoeffs)
    bpr = similar(flowVec)
    for a = 1:length(bpr)
        bpr[a] = free_flow_time[a] * sum([fcoeffs[i] * (flowVec[a]/capacity[a])^(i-1) for i = 1:length(fcoeffs)])
    end
    return bpr
end

function BPRSocial(flowVec, fcoeffs)
    bpr = similar(flowVec)
    # refer to [Page 50; Patriksson 1994, 2015]
    for a = 1:length(bpr)
        bpr[a] = free_flow_time[a] * sum([fcoeffs[i] *
        (flowVec[a]/capacity[a])^(i-1) for i = 1:length(fcoeffs)])
        + free_flow_time[a] * sum([fcoeffs[i] * (i-1) *
        (flowVec[a]/capacity[a])^(i-1) for i = 2:length(fcoeffs)])
    end
    return bpr
end

function all_or_nothing(travel_time, demands)
    state = []
    path = []
    x = zeros(size(start_node))

    for r=1:numZones
        # for each origin node r, find shortest paths to all destination nodes
        state = dijkstra_shortest_paths(graph, travel_time, r)

        for s=1:numZones
            # for each destination node s, find the shortest-path vector
            # load travel demand
            x = x + demands[(r,s)] * get_vector(state, r, s, link_dic)
        end
    end

    return x
end

function tapMSA(demands, fcoeffs, numIter=500, tol=1e-6)
    # Finding a starting feasible solution
    travel_time = BPR(zeros(numLinks), fcoeffs)
    xl = all_or_nothing(travel_time, demands)

    l = 1

    while l < numIter
        l += 1

        xl_old = xl

        # Finding yl
        travel_time = BPR(xl, fcoeffs)

        yl = all_or_nothing(travel_time, demands)

        # assert(yl != xl)

        xl = xl + (yl - xl)/l

        xl_new = xl

        relative_gap = norm(xl_new - xl_old, 1) / norm(xl_new, 1)

        if relative_gap < tol
            break
        end

    end

    tapFlows = Dict()

    for i = 1:length(ta_data.start_node)
        key = (ta_data.start_node[i], ta_data.end_node[i])
        tapFlows[key] = xl[i]
    end

    tapFlowVect = xl

    return tapFlows, tapFlowVect
end

function tapMSASocial(demands, fcoeffs, numIter=1000, tol=1e-6)
    # Finding a starting feasible solution
    travel_time = BPRSocial(zeros(numLinks), fcoeffs)
    xl = all_or_nothing(travel_time, demands)

    l = 1

    while l < numIter
        l += 1

        xl_old = xl

        # Finding yl
        travel_time = BPRSocial(xl, fcoeffs)

        yl = all_or_nothing(travel_time, demands)

        # assert(yl != xl)

        xl = xl + (yl - xl)/l

        xl_new = xl

        relative_gap = norm(xl_new - xl_old, 1) / norm(xl_new, 1)

        if relative_gap < tol
            break
        end

    end

    tapFlows = Dict()

    for i = 1:length(ta_data.start_node)
        key = (ta_data.start_node[i], ta_data.end_node[i])
        tapFlows[key] = xl[i]
    end

    tapFlowVect = xl

    return tapFlows, tapFlowVect
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.992552298143351 <= q <= 4.750344613921084)

                     add_NL_constraint(m, :(swish(0.47930083070023954 + -0.2811056889538768*swish(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])) + -0.8918376604733185*swish(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])) + -0.8083246377220759*swish(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))) + swish(-0.740550943643274 + -0.007926005977438866*swish(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])) + 0.328386397895855*swish(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])) + 0.6266373732109729*swish(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))) + swish(-0.3322291763851517 + -0.3879865017411772*swish(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])) + 0.13951112377587105*swish(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])) + -0.16669374723112762*swish(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#= -------------------------------------------------------------------------
# @file solution.jl
#
# @date 10/28/16 23:51:31
# @author Martin Noblia
# @email martin.noblia@openmailbox.org
#
# @brief
#
# @detail
#
  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#
module Reverse

"""
Function for returning the lines of text file
"""
function read_file(str::AbstractString)
    file = open(str)
    lines = eachline(file)
    close(file)
    return lines
end

function lines(io::IO)
    lines_array = Array{String, 1}()
    for line in eachline(io)
        push!(lines_array, chomp(line))
    end
    return lines_array
end

lines(path::AbstractString) = open(lines, path)

"""
Function for generate the cases
"""
function create_cases(file_name::String)

   data = lines(file_name)
   number_of_cases = parse(Int, data[1])

   S = [reverse(split(s)) for s in data[2:end]] # S Strings in reverse
   L = Array{String, 1}()
   for s in S
      l = length(s)
      if l > 1
         # concatenate with products
         aux = prod(s .* " ")
         aux = strip(aux) # strip the last whitespace introduce
         push!(L, aux)
      end
      if l == 1
         push!(L, s[1])
      end
   end
   return L, S
end

"""
Function for print the cases results

"""
function print_results(L::Array{String, 1})

   for (num, s) in enumerate(L)
      println("Case #$num: $s")
   end

end

end

max_stars_repo_name: elsuizo/CodeJam

max_issues_repo_name: elsuizo/CodeJam

max_forks_repo_name: elsuizo/CodeJam


text:
export
    q2args_sqth_fno,
    q2ρ_sqth_fno,
    q2args_sqth_res

function q2args_sqth_fno()
    modes = (12, )
    ch = 64
    σ = leakyrelu

    return Chain(
        Dense(2, ch),
        FourierOperator(ch=>ch, modes, σ),
        FourierOperator(ch=>ch, modes, σ),
        FourierOperator(ch=>ch, modes, σ),
        FourierOperator(ch=>ch, modes),
        Dense(ch, 4ch, σ),
        Dense(4ch, 1),

        flatten,
        Dense(4096, 3),
    )
end

function q2ρ_sqth_fno()
    modes = (12, )
    ch = 64
    σ = leakyrelu

    return Chain(
        Dense(2, ch),
        FourierOperator(ch=>ch, modes, σ),
        FourierOperator(ch=>ch, modes, σ),
        FourierOperator(ch=>ch, modes, σ),
        FourierOperator(ch=>ch, modes),
        Dense(ch, 4ch, σ),
        Dense(4ch, 1),

        flatten,
        Dense(4096, 100*100),
        x -> cholesky2ρ(x)
    )
end

function q2args_sqth_res()
    σ = relu

    return Chain(
        # stage 0
        Conv((31, ), 1=>8, pad=15),
        BatchNorm(8, σ),

        # res
        res_block((8, 4, 4, 16), (1, 15, 7), (0, 7, 3), 1, 0, 0, σ),
        res_block((16, 8, 8, 32), (1, 15, 7), (0, 7, 3), 1, 0, 0, σ),
        res_block((32, 16, 16, 64), (1, 15, 7), (0, 7, 3), 1, 0, 0, σ),
        res_block((64, 16, 16, 32), (1, 15, 7), (0, 7, 3), 1, 0, 4, σ),
        res_block((32, 8, 8, 16), (1, 15, 7), (0, 7, 3), 1, 0, 4, σ),
        res_block((16, 4, 4, 8), (1, 15, 7), (0, 7, 3), 1, 0, 4, σ),

        # stage 1
        flatten,
        Dense(8*64, 64, σ),
        Dense(64, 16, σ),
        Dense(16, 3, relu)
    )
end


text:
# import stuff
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using Statistics
@everywhere using GRASS
@everywhere using SharedArrays
@everywhere using EchelleCCFs
using CSV
using DataFrames
using LaTeXStrings
using LsqFit

# define rms loop function
include(GRASS.moddir * "figures/fig_functions.jl")

# some global stuff
const N = 132
const Nt = 200
const Nloop = 200

# get command line args and output directories
run, plot = parse_args(ARGS)
grassdir, plotdir, datadir = check_plot_dirs()

function main()
    # set up stuff for lines
    lines = [5434.5]
    depths = [0.8]
    resolution = 700000.0
    top = NaN
    contiguous_only=false

    # get angles
    n_inc = 20
    incls = (π/180.0) .* range(0.0, 90.0, length=n_inc)
    ycomp = sin.(incls)
    zcomp = cos.(incls)

    # make pole vectors
    poles = Array{Tuple{Float64, Float64, Float64}, 1}(undef, n_inc)
    for i in eachindex(incls)
        pole = (0.0, ycomp[i], zcomp[i])
        poles[i] = pole
    end

    # allocate shared sharred arrays
    avg_avg_inc = SharedArray{Float64}(n_inc)
    std_avg_inc = SharedArray{Float64}(n_inc)
    avg_rms_inc = SharedArray{Float64}(n_inc)
    std_rms_inc = SharedArray{Float64}(n_inc)

    # loop over inclinations
    @sync @distributed for i in eachindex(poles)
        println("running pole " * string(i) * " of " * string(n_inc))

         # create spec and disk params instances
        spec = SpecParams(lines=lines, depths=depths, resolution=resolution,
                          extrapolate=true, contiguous_only=contiguous_only)
        disk = DiskParams(N=N, Nt=Nt, pole=poles[i])

        # synthesize spectra, get velocities and stats
        avg_avg1, std_avg1, avg_rms1, std_rms1 = spec_loop(spec, disk, Nloop, top=top)
        avg_avg_inc[i] = avg_avg1
        std_avg_inc[i] = std_avg1
        avg_rms_inc[i] = avg_rms1
        std_rms_inc[i] = std_rms1
    end

    # make data frame
    df = DataFrame()
    df[!,:inc] = incls
    df[!,:avg_avg_inc] = avg_avg_inc
    df[!,:std_avg_inc] = std_avg_inc
    df[!,:avg_rms_inc] = avg_rms_inc
    df[!,:std_rms_inc] = std_rms_inc

    # write results to CSV
    fname = datadir * "inclination_" *  string(N) * ".csv"
    CSV.write(fname, df)
    return nothing
end

# run the simulation
if run
    main()
end

# plotting code block
if plot
    # plotting imports
    import PyPlot; plt = PyPlot; mpl = plt.matplotlib; plt.ioff()
    using PyCall; animation = pyimport("matplotlib.animation")
    mpl.style.use(GRASS.moddir * "figures/fig.mplstyle")

    # read in the data
    fname = datadir * "inclination_" *  string(N) * ".csv"
    df = CSV.read(fname, DataFrame)

    # assign to variable names
    ang = df.inc .* (180.0/π)
    avg_avg_inc = df.avg_avg_inc
    std_avg_inc = df.std_avg_inc
    avg_rms_inc = df.avg_rms_inc
    std_rms_inc = df.std_rms_inc

    # get the errors
    err_avg_inc = std_avg_inc ./ sqrt(Nloop)
    err_rms_inc = std_rms_inc ./ sqrt(Nloop)

    # plot the results
    fig = plt.figure()
    ax1 = fig.add_subplot()
    ax1.errorbar(ang, avg_rms_inc, yerr=err_rms_inc, capsize=3.0, color="black", fmt=".")
    ax1.fill_between(ang, avg_rms_inc .- std_rms_inc, avg_rms_inc .+ std_rms_inc, color="tab:blue", alpha=0.3)

    # annotate the axes
    arrowprops = Dict("facecolor"=>"black", "shrink"=>0.05, "width"=>2.0,"headwidth"=>8.0)
    ax1.annotate(L"\textnormal{Pole-on}", xy=(69.8, 0.39), xytext=(0.0,0.388), arrowprops=arrowprops)
    ax1.annotate(L"\textnormal{Equator-on}", xy=(70.0, 0.388))

    # set labels, etc.
    ax1.set_xlabel(L"{\rm Inclination\ (deg)}")
    ax1.set_ylabel(L"{\rm RMS}_{\rm RV}\ {\rm (m s}^{-1})")
    ax1.set_xticks(range(0, 90, length=10))
    ax1.set_ylim(0.375, 0.675)

    # save the figure
    fig.savefig(plotdir * "fig6.pdf")
    plt.clf(); plt.close()
    println(">>> Figure written to: " * plotdir * "fig6.pdf")
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
struct ZeroField <: AbstractField{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, 3} end

@inline Base.getindex(::ZeroField, i, j, k) = 0


text:
function KolmogorovSmirnovTest()
    return KolmogorovSmirnovTest(())
end

function KolmogorovSmirnovTest(arg0::jlong)
    return KolmogorovSmirnovTest((jlong,), arg0)
end

function approximate_p(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jint, arg2::jint)
    return jcall(obj, "approximateP", jdouble, (jdouble, jint, jint), arg0, arg1, arg2)
end

function bootstrap(obj::KolmogorovSmirnovTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint)
    return jcall(obj, "bootstrap", jdouble, (Vector{jdouble}, Vector{jdouble}, jint), arg0, arg1, arg2)
end

function bootstrap(obj::KolmogorovSmirnovTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jint, arg3::jboolean)
    return jcall(obj, "bootstrap", jdouble, (Vector{jdouble}, Vector{jdouble}, jint, jboolean), arg0, arg1, arg2, arg3)
end

function cdf(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jint)
    return jcall(obj, "cdf", jdouble, (jdouble, jint), arg0, arg1)
end

function cdf(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jint, arg2::jboolean)
    return jcall(obj, "cdf", jdouble, (jdouble, jint, jboolean), arg0, arg1, arg2)
end

function cdf_exact(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jint)
    return jcall(obj, "cdfExact", jdouble, (jdouble, jint), arg0, arg1)
end

function exact_p(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jint, arg2::jint, arg3::jboolean)
    return jcall(obj, "exactP", jdouble, (jdouble, jint, jint, jboolean), arg0, arg1, arg2, arg3)
end

function kolmogorov_smirnov_statistic(obj::KolmogorovSmirnovTest, arg0::RealDistribution, arg1::Vector{jdouble})
    return jcall(obj, "kolmogorovSmirnovStatistic", jdouble, (RealDistribution, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_statistic(obj::KolmogorovSmirnovTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "kolmogorovSmirnovStatistic", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_test(obj::KolmogorovSmirnovTest, arg0::RealDistribution, arg1::Vector{jdouble})
    return jcall(obj, "kolmogorovSmirnovTest", jdouble, (RealDistribution, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_test(obj::KolmogorovSmirnovTest, arg0::RealDistribution, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(obj, "kolmogorovSmirnovTest", jdouble, (RealDistribution, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function kolmogorov_smirnov_test(obj::KolmogorovSmirnovTest, arg0::RealDistribution, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "kolmogorovSmirnovTest", jboolean, (RealDistribution, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function kolmogorov_smirnov_test(obj::KolmogorovSmirnovTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "kolmogorovSmirnovTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function kolmogorov_smirnov_test(obj::KolmogorovSmirnovTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(obj, "kolmogorovSmirnovTest", jdouble, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function ks_sum(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jdouble, arg2::jint)
    return jcall(obj, "ksSum", jdouble, (jdouble, jdouble, jint), arg0, arg1, arg2)
end

function pelz_good(obj::KolmogorovSmirnovTest, arg0::jdouble, arg1::jint)
    return jcall(obj, "pelzGood", jdouble, (jdouble, jint), arg0, arg1)
end


text:
module CUSOLVER

using CUDArt
using CUDArt.rt.cudaStream_t

import Base.one
import Base.zero

include("libcusolver_types.jl")

immutable CUSOLVERError <: Exception
    msg::AbstractString
    status::UInt32

    function CUSOLVERError(status)
        new(status,statusmessage(status))
    end
end

function statusmessage( status )
    if status == CUSOLVER_STATUS_SUCCESS
        return "cusolver success"
    elseif status == CUSOLVER_STATUS_NOT_INITIALIZED
        return "cusolver not initialized"
    elseif status == CUSOLVER_STATUS_ALLOC_FAILED
        return "cusolver allocation failed"
    elseif status == CUSOLVER_STATUS_INVALID_VALUE
        return "cusolver invalid value"
    elseif status == CUSOLVER_STATUS_ARCH_MISMATCH
        return "cusolver architecture mismatch"
    elseif status == CUSOLVER_STATUS_EXECUTION_FAILED
        return "cusolver execution failed"
    elseif status == CUSOLVER_STATUS_INTERNAL_ERROR
        return "cusolver internal error"
    elseif status == CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED
        return "cusolver matrix type not supported"
    end
end

function statuscheck( status )
    if status == CUSOLVER_STATUS_SUCCESS
        return nothing
    end
    warn("CUSOLVER error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    println()
    throw(CUSOLVERError( status ))
end

cuda_versions = ["80", "75", "70", "65"]
lib_list = vcat(["libcusolver", "cusolver"], "cusolver64_" .* cuda_versions, "cusolver32_" .* cuda_versions)
const libcusolver = Libdl.find_library(lib_list, ["/usr/local/cuda", "/usr/lib", ""])
if isempty(libcusolver)
    error("CUSOLVER library not found!")
end

include("libcusolver.jl")

#setup handler for cusolver

cusolverSphandle = cusolverSpHandle_t[0]
cusolverSpCreate( cusolverSphandle )
cusolverDnhandle = cusolverDnHandle_t[0]
cusolverDnCreate( cusolverDnhandle )

function cusolverDestroy()
    cusolverSpDestroy(cusolverSphandle[1])
    cusolverDnDestroy(cusolverDnhandle[1])
end
#clean up handle at exit
atexit( ()->cusolverDestroy() )

include("sparse.jl")
include("dense.jl")

end


text:
@testset "legendre correct values edge,origin" begin
    # legendre == jacobi(a,b=0,0); same function body
    for n in 1:10
        # design of jacobi polynomials is such that edge = 1
        @test legendre(n, [1.]) ≈ [1.]
        correct = [-1.]
        if iseven(n)
            correct = [1.]
        end
        @test legendre(n, [-1.]) ≈ correct
        if !iseven(n)
            @test legendre(n, [0.]) ≈ [0.]
        end
    end
end

@testset "legendre sum correct values" begin
    # legendre == jacobi(a,b=0,0), same function body
    terms = [1,3,5,7,9]
    weights = ones(size(terms)...)
    @test legendre_sum(terms, weights, [0.]) ≈ [0.]
    terms = [1,2,3,4,5]
    weights = ones(size(terms)...)
    # -1 accumulates 3 times and +1 two times on left edge => expect -1.
    # 1 accumulates 5 times on right edge => expect 5.
    @test legendre_sum(terms, weights, [-1., 1.]) ≈ [-1., 5.]
end

@testset "legendre series correct values contiguous" begin
    terms = [1,2,3,4,5]
    ledge = [-1.,1.,-1.,1.,-1.]
    center = [0.,-0.5,0.,0.375,0.]
    redge = [1.,1.,1.,1.,1.]
    series = legendre_series(terms, [-1., 0., 1.])
    @test series[1,:] ≈ ledge
    @test series[2,:] ≈ center
    @test series[3,:] ≈ redge
end

@testset "legendre series correct values sparse" begin
    terms = [1,3,5]
    ledge = [-1.,-1.,-1.]
    center = [0.,0.,0.]
    redge = [1.,1.,1.]
    series = legendre_series(terms, [-1., 0., 1.])
    @test series[1,:] ≈ ledge
    @test series[2,:] ≈ center
    @test series[3,:] ≈ redge
end


text:
using CImGui
using CImGui.LibCImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(1280, 720, "Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsDark()
# CImGui.StyleColorsClassic()
# CImGui.StyleColorsLight()

# load Fonts
# - If no fonts are loaded, dear imgui will use the default font. You can also load multiple fonts and use `CImGui.PushFont/PopFont` to select them.
# - `CImGui.AddFontFromFileTTF` will return the `Ptr{ImFont}` so you can store it if you need to select the font among multiple.
# - If the file cannot be loaded, the function will return C_NULL. Please handle those errors in your application (e.g. use an assertion, or display an error and quit).
# - The fonts will be rasterized at a given size (w/ oversampling) and stored into a texture when calling `CImGui.Build()`/`GetTexDataAsXXXX()``, which `ImGui_ImplXXXX_NewFrame` below will call.
# - Read 'fonts/README.txt' for more instructions and details.
#fonts_dir = joinpath(@__DIR__, "..", "fonts")
#fonts = CImGui.GetIO().Fonts
# default_font = CImGui.AddFontDefault(fonts)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Cousine-Regular.ttf"), 15)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "DroidSans.ttf"), 16)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Karla-Regular.ttf"), 10)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "ProggyTiny.ttf"), 10)
#CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)
# @assert default_font != C_NULL

#fonts_dir = joinpath(pathof(CImGui), "..","..","fonts")
fonts_dir = raw"c:\Windows\Fonts";
fonts = CImGui.GetIO().Fonts

builder=ImFontGlyphRangesBuilder_ImFontGlyphRangesBuilder()
ImFontGlyphRangesBuilder_AddText(builder, "♐", C_NULL)
ImFontGlyphRangesBuilder_AddChar(builder,'☽')
ImFontGlyphRangesBuilder_AddRanges(builder, ImFontAtlas_GetGlyphRangesDefault(fonts))

struct GlyphRange
	first::ImWchar
	last::ImWchar
	zero::ImWchar
	function GlyphRange(first::ImWchar,last::ImWchar)
		new(first,last,0x0)
	end
end
greekRange=GlyphRange(0x0300, 0x03ff)
greekRange_ptr=Core.Ptr{ImWchar}(Base.pointer_from_objref(Ref(greekRange)))
ImFontGlyphRangesBuilder_AddRanges(builder, greekRange_ptr)

ranges=ImVector_ImWchar_create()
ImVector_ImWchar_Init(ranges)
ImFontGlyphRangesBuilder_BuildRanges(builder, ranges)
r = unsafe_wrap(Vector{ImVector_ImWchar}, ranges, 1)
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "JuliaMono-Regular.ttf"), 16, C_NULL, r[1].Data )

CImGui.Build(fonts)
default_font = CImGui.AddFontDefault(fonts)
#@assert default_font != C_NULL

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    show_demo_window = true
    show_another_window = false
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        # show the big demo window
        show_demo_window && @c CImGui.ShowDemoWindow(&show_demo_window)

        # show a simple window that we create ourselves.
        # we use a Begin/End pair to created a named window.
        @cstatic f=Cfloat(0.0) counter=Cint(0) begin
            CImGui.Begin("Hello, world!")  # create a window called "Hello, world!" and append into it.
            CImGui.Text("This is some useful text.")  # display some text
            @c CImGui.Checkbox("Demo Window", &show_demo_window)  # edit bools storing our window open/close state
            @c CImGui.Checkbox("Another Window", &show_another_window)

            @c CImGui.SliderFloat("float", &f, 0, 1)  # edit 1 float using a slider from 0 to 1
            CImGui.ColorEdit3("clear color", clear_color)  # edit 3 floats representing a color
            CImGui.Button("Button") && (counter += 1)

            CImGui.SameLine()
            CImGui.Text("counter = $counter")
            CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))

            CImGui.End()
        end

        # show another simple window.
        if show_another_window
            @c CImGui.Begin("Another Window", &show_another_window)  # pass a pointer to our bool variable (the window will have a closing button that will clear the bool when clicked)
            CImGui.Text("Hello from another window!")
            CImGui.Button("Close Me") && (show_another_window = false;)
            CImGui.End()
        end

		CImGui.Begin("Special characters")
		CImGui.Text("π ϕ θ ♐ ☽")
		CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end
catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    ImVector_ImWchar_destroy(ranges)
    ImFontGlyphRangesBuilder_destroy(builder)
	CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end