text:
using BioMedQuery.PubMed
using BioServices.EUtils
using SQLite
using MySQL
using XMLDict
using LightXML

"""
    pubmed_search_and_save!(email, search_term::String, article_max,
    conn, verbose=false)

Search pubmed and save the results into a database connection. The database is expected to exist and
have the appriate pubmed related tables. You can create such tables using
`PubMed.create_tables(conn)`

## Arguments

* `email` : valid email address (otherwise pubmed may block you)
* `search_term` : search string to submit to PubMed e.g `(asthma[MeSH Terms]) AND ("2001/01/29"[Date - Publication] : "2010"[Date - Publication])` see http://www.ncbi.nlm.nih.gov/pubmed/advanced for help constructing the string
* `article_max` : maximum number of articles to return
* `conn` : database connection
* `verbose` : if true, the NCBI xml response files are saved to current directory
"""
function pubmed_search_and_save!(email, search_term::String, article_max,
    conn, verbose=false)

    retstart = 0
    retmax = 10000  #e-utils only allows 10,000 at a time
    article_max = article_max

    if article_max < retmax
        retmax = article_max
    end

    db = nothing
    article_total = 0

    for rs=retstart:retmax:(article_max- 1)

        rm = rs + retmax
        if rm > article_max
            retmax = article_max - rs
        end

        println("Getting ", retmax, " articles, starting at index ", rs)

        println("------ESearch--------")

        esearch_response = esearch(db = "pubmed", term = search_term,
        retstart = rs, retmax = retmax, tool = "BioJulia",
        email = email)

        # if verbose
        #     xdoc = parse_string(esearch_response)
        #     save_file(xdoc, "./esearch.xml")
        # end

        #convert xml to dictionary
        esearch_dict = parse_xml(String(esearch_response.body))

        if !haskey(esearch_dict, "IdList")
            println("Error with esearch_dict:")
            println(esearch_dict)
            error("Response esearch_dict does not contain IdList")
        end

        println("------EFetch--------")

        #get the list of ids and perfom a fetch
        ids = [parse(Int64, id_node) for id_node in esearch_dict["IdList"]["Id"]]

        efetch_response = efetch(db = "pubmed", tool = "BioJulia", retmode = "xml", rettype = "null", id = ids)

        # if verbose
        #     xdoc = parse_string(efetch_response)
        #     save_file(xdoc, "./efetch.xml")
        # end

        #convert xml to dictionary
        efetch_doc = root(parse_string(String(efetch_response.body)))

        #save the results of an entrez fetch
        println("------Save to database--------")

        save_efetch!(conn, efetch_doc, verbose, true)

        free(efetch_doc)

        article_total+=length(ids)

        if (length(ids) < retmax)
            break
        end

    end

    println("Finished searching, total number of articles: ", article_total)
    return db
end

function pubmed_pmid_search(email, search_term::String, article_max, verbose=false)

     retstart = 0
     retmax = 10000  #e-utils only allows 10,000 at a time
     article_max = article_max

     if article_max < retmax
         retmax = article_max
     end

     all_pmids = Array{Int64,1}()
     article_total = 0

     for rs=retstart:retmax:(article_max- 1)

        rm = rs + retmax
        if rm > article_max
            retmax = article_max - rs
        end

        println("Getting ", retmax, " articles, starting at index ", rs)

        println("------ESearch--------")

        esearch_response = esearch(db = "pubmed", term = search_term,
        retstart = rs, retmax = retmax, tool = "BioJulia",
        email = email)

        if verbose
            xdoc = parsexml(esearch_response)
            save_file(xdoc, "./esearch.xml")
        end

        #convert xml to dictionary
        esearch_dict = parse_xml(String(esearch_response.body))

        if !haskey(esearch_dict, "IdList")
            println("Error with esearch_dict:")
            println(esearch_dict)
            error("Response esearch_dict does not contain IdList")
        end

        #get the list of ids and perfom a fetch
        ids = [parse(Int64, id_node) for id_node in esearch_dict["IdList"]["Id"]]

        append!(all_pmids, ids)

        article_total+=length(ids)

        if (length(ids) < retmax)
            break
        end

     end

     println("Finished searching, total number of articles: ", article_total)

     return all_pmids
 end

function pubmed_pmid_search_and_save!(email, search_term::String, article_max,
    conn, verbose=false)

    pmids = pubmed_pmid_search(email, search_term, article_max, verbose)
    db = save_pmids!(conn, pmids, verbose)

    narticles = length(pmids)
    @info "Finished saving $narticles articles"

    return nothing
end

"""
    pubmed_search_and_parse(email, search_term::String, article_max, verbose=false)

Search pubmed and parse the results into a dictionary of DataFrames.  The dataframes have the same names
and fields as the pubmed database schema. (e.g. df_dict["basic"] returns a dataframe with the basic article info)

## Arguments

* `email` : valid email address (otherwise pubmed may block you)
* `search_term` : search string to submit to PubMed e.g `(asthma[MeSH Terms]) AND ("2001/01/29"[Date - Publication] : "2010"[Date - Publication])` see http://www.ncbi.nlm.nih.gov/pubmed/advanced for help constructing the string
* `article_max` : maximum number of articles to return
* `verbose` : if true, the NCBI xml response files are saved to current directory
"""
function pubmed_search_and_parse(email, search_term::String, article_max, verbose=false)

    retstart = 0
    retmax = 10000  #e-utils only allows 10,000 at a time
    article_max = article_max

    dfs = Dict{String, DataFrame}()

    if article_max < retmax
        retmax = article_max
    end

    db = nothing
    article_total = 0

    for rs=retstart:retmax:(article_max- 1)

        rm = rs + retmax
        if rm > article_max
            retmax = article_max - rs
        end

        println("Getting ", retmax, " articles, starting at index ", rs)

        println("------ESearch--------")

        esearch_response = esearch(db = "pubmed", term = search_term,
        retstart = rs, retmax = retmax, tool = "BioJulia",
        email = email)

        # if verbose
        #     xdoc = parse_string(esearch_response)
        #     save_file(xdoc, "./esearch.xml")
        # end

        #convert xml to dictionary
        esearch_dict = parse_xml(String(esearch_response.body))

        if !haskey(esearch_dict, "IdList")
            println("Error with esearch_dict:")
            println(esearch_dict)
            error("Response esearch_dict does not contain IdList")
        end

        println("------EFetch--------")

        #get the list of ids and perfom a fetch
        ids = [parse(Int64, id_node) for id_node in esearch_dict["IdList"]["Id"]]

        efetch_response = efetch(db = "pubmed", tool = "BioJulia", retmode = "xml", rettype = "null", id = ids)

        # if verbose
        #     xdoc = parse_string(efetch_response)
        #     save_file(xdoc, "./efetch.xml")
        # end

        #convert xml to dictionary
        efetch_doc = root(parse_string(String(efetch_response.body)))

        #save the results of an entrez fetch
        println("------Save to dataframes--------")
        this_dfs = PubMed.parse_articles(efetch_doc)
        for (table, df) in this_dfs
            if haskey(dfs, table)
                dfs[table] = vcat(dfs[table], df)
            else
                dfs[table] = df
            end
        end

        free(efetch_doc)

        article_total+=length(ids)

        if (length(ids) < retmax)
            break
        end
    end

    return dfs
end

max_stars_repo_name: bcbi/NLM


text:
module DrMZ

# Load required packages
using FFTW: fft, ifft
using Flux
using Flux.Data: DataLoader
using Flux: mse
using ProgressMeter: @showprogress
using Distributions: MvNormal
using LinearAlgebra
using Random: randperm, AbstractRNG
using Random
using DifferentialEquations
using Interpolations
using BSON: @save
using BSON: @load
using Printf
using ToeplitzMatrices
using Statistics: mean, std
using FastGaussQuadrature
using ForwardDiff
using LegendrePolynomials

# Export the functions for General.jl
export error_se, error_rel, mse_error, norm_rel_error, norm_rel_error_continuous, norm_infinity_error, ic_error, average_ic_error, average_error,
       periodic_fill_domain, periodic_fill_solution, solution_interpolation, reduced_initial_condition, solution_spatial_sampling, solution_temporal_sampling,
       fft_norm, ifft_norm,
       fourier_diff,
       trapz,
       gauss_legendre,
       clenshaw_curtis, cheby_grid, cheby_diff_matrix, cheby_diff,
       trapezoid,
       orthonormal_check,
       legendre_norm, legendre_norm_collect

# Export the functions for OperatorNN.jl
export train_model, loss_all, predict,
       build_dense_model, build_branch_model, build_trunk_model,
       exp_kernel_periodic, generate_sinusoidal_functions_2_parameter, generate_periodic_train_test_initial_conditions, generate_periodic_train_test, generate_periodic_train_test_initial_condition_load, generate_periodic_functions, solution_extraction,
       save_model, load_model, load_branch, load_trunk, save_data, load_data, load_data_initial_conditions, load_data_train_test, load_data_initial_conditions, save_data_initial_conditions

# Export the functions for PDESolve.jl
export advection_pde!, advection_diffusion_pde!, inviscid_burgers_pde!, viscous_burgers_pde!,
       quadratic_nonlinear,
       generate_fourier_solution, generate_periodic_train_test_muscl,
       central_difference,
       minmod, ub, fl, ulpl, urpl, ulnl, urnl,
       muscl_minmod_RHS!, muscl_minmod_viscous_RHS!, generate_muscl_minmod_solution, generate_muscl_reduced,
       get_1D_energy_fft,
       get_1D_energy_custom, get_1D_energy_custom_coefficients,
       mode_extractor, get_1D_energy_upwind,
       spectral_approximation_fourier

# Export the functions for DBasis.jl
export trunk_build, trunk_ortho_build, build_basis,
       basis_eval, expansion_coefficients, expansion_approximation,
       basis_derivative,
       save_basis, load_basis

# Load functions
include("General.jl")
include("OperatorNN.jl")
include("PDESolve.jl")
include("DBasis.jl")

end

ext: jl


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
Formulation described in:

    Garver, L. L. (1962). Power generation scheduling by integer
    programming-development of theory. Transactions of the American Institute
    of Electrical Engineers. Part III: Power Apparatus and Systems, 81(3), 730-734.
    DOI: https://doi.org/10.1109/AIEEPAS.1962.4501405

"""
module Gar1962

import ..PiecewiseLinearCostsFormulation
import ..ProductionVarsFormulation
import ..StatusVarsFormulation

struct ProdVars <: ProductionVarsFormulation end
struct PwlCosts <: PiecewiseLinearCostsFormulation end
struct StatusVars <: StatusVarsFormulation end

end


text:
export newton, linop

struct DualFun{F,T}
    f::F
    J::T
end
DualFun(f::Fun) = DualFun(f,Operator(I,space(f)))

domain(df::DualFun) = domain(df.f)

differentiate(d::DualFun) = DualFun(d.f',Derivative(rangespace(d.J))*d.J)
integrate(d::DualFun) = DualFun(integrate(d.f),Integral(rangespace(d.J))*d.J)
function cumsum(d::DualFun)
    Q = Integral(rangespace(d.J))*d.J
    rs = rangespace(Q)
    DualFun(cumsum(d.f),Q-ones(rs)*Evaluation(rs,leftendpoint)*Q)
end

# total definite integral of a function
function sum(df::DualFun)
    Q = Integral(rangespace(df.J))*df.J
    return DualFun(sum(df.f), Evaluation(rangespace(Q), rightendpoint)*Q)
end

adjoint(d::DualFun) = differentiate(d)

^(d::DualFun, k::Number) = DualFun(d.f^k,k*d.f^(k-1)*d.J)

# from DualNumbers
for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
    @eval function $(funsym)(z::DualFun)
        x = z.f
        xp = z.J
        DualFun($(funsym)(x),$exp*xp)
    end
end

for OP in (:+,:-)
    @eval begin
        $OP(a::DualFun,b::Union{Number,Fun}) = DualFun($OP(a.f,b),a.J)
        $OP(a::Union{Number,Fun},b::DualFun) = DualFun($OP(a,b.f),$OP(b.J))
        $OP(a::DualFun,b::DualFun) = DualFun($OP(a.f,b.f),$OP(a.J,b.J))
    end
end
-(a::DualFun)=DualFun(-a.f,-a.J)

*(a::Union{Number,Fun},b::DualFun)=DualFun(a*b.f,a*b.J)
*(a::DualFun,b::Union{Number,Fun})=DualFun(b*a.f,b*a.J)
*(a::DualFun,b::DualFun)=DualFun(a.f*b.f,a.f*b.J+b.f*a.J)

/(a::Union{Number,Fun},b::DualFun)=DualFun(a/b.f,-a/b.f^2*b.J)
/(a::DualFun,b::Union{Number,Fun})=DualFun(a.f/b,a.J/b)
/(a::DualFun,b::DualFun)=DualFun(a.f/b.f,a.J/b.f-a.f/b.f^2*b.J)

(d::DualFun)(x) = DualFun(d.f(x),Evaluation(rangespace(d.J),x)*d.J)
first(d::DualFun) = DualFun(first(d.f),Evaluation(rangespace(d.J),leftendpoint)*d.J)
last(d::DualFun) = DualFun(last(d.f),Evaluation(rangespace(d.J),rightendpoint)*d.J)

jacobian(d::DualFun) = d.J
jacobian(a::Number) = zero(a)
jacobian(f::Fun) = Operator(I,space(f))

promote_rule(::Type{DF},::Type{T}) where {DF<:DualFun,T<:Number}=DualFun
convert(::Type{DualFun},b::Number) = DualFun(b,0)

function linop(f::Function,ds::Space)
    if (isgeneric(f) && applicable(f,0)) || (!isgeneric(f)&&arglength(f)==1)
        df=f(DualFun(zeros(ds)))
    elseif (isgeneric(f) && applicable(f,0,0)) || (!isgeneric(f)&&arglength(f)==2)
        df=f(Fun(ds),DualFun(zeros(ds)))
    else
        error("Not implemented")
    end

    if isa(df,Array)
        map(u->u.J,df)
    else
        df.J
    end
end

linop(f::Function,d) = linop(f,Space(d))
linop(f::Function) = linop(f,Chebyshev())  #TODO: UnsetSpace

# full operator should be
# N=u->[B*u-bcs;...]
function newton(N,u0::Fun;maxiterations=15,tolerance=1E-15)
    u=u0
    err=Inf
    for k=1:maxiterations
        DF=N(DualFun(u))
        J=map(jacobian,DF)
        F=map(d->d.f,DF)
        unew=u-J\F
        err=norm(unew-u)
        if err≤10tolerance
            return unew
        else
            u=chop(unew,tolerance)
        end
    end
    @warn "Maximum number of iterations $maxiterations reached, with approximate accuracy of $err."
    return u
end

# Newton iteration for a system of nonlinear equations
# INPUTS:
#   N - a nonlinear function N(u1, u2, u3, ...) = [BCs;...]
#   u0 - an array of initial Funs
#       NOTE: these initial Funs should be "chopped" as much as possible for optimal speed
#
# KSS and ZJS
function newton(N, u0::Array{T,1}; maxiterations=15, tolerance=1E-15) where T <: Fun
    u = copy(u0)
    err = Inf

    numf = length(u0)   # number of functions

    if numf == 0
        error("u0 must contain at least 1 function")
    end

    Js = Array{Any,2}(undef,1,numf) # jacobians

    for k = 1:maxiterations
        # ------ calculate Jacobian for each function - O(numf^2) ----

        # use the first call to populate the value of the functions
        # this saves one call to N()
        F = N([ (j == 1 ? DualFun(u[j]) : u[j]) for j = 1:numf ]...)
        @inbounds Js[1] = jacobian.(F)
        F = map(d -> typeof(d) <: DualFun ? d.f : d, F)

        for i = 2:numf
            @inbounds Js[i] = jacobian.(N([ (j == i ? DualFun(u[j]) : u[j]) for j = 1:numf ]...))
        end

        F = N(u...)
        J = Base.typed_hcat(Operator, Js...)

        # ------- update step ---------
        unew = u .- J\F
        err = maximum(norm.(unew-u)) # use the max norm error as the overall error

        if err ≤ 10*tolerance
            return unew
        else
            u = map(q -> chop(q, tolerance), unew)
        end
    end

    @warn "Maximum number of iterations $maxiterations reached, with approximate accuracy of $err."

    return u
end

newton(N,d::Domain;opts...) =
    newton(N,zeros(d);opts...)

newton(N,d;opts...) =
    newton(N,Domain(d);opts...)


text:
# This file contains operations that modify the presolve problem data

"
--- Row and Column Operations ---
There are four functions for each - add, enque, deque and remove.
There are separate functions for enque and deque intentionally.
There are times we want to remove a row from the active list but not remove it from the problem.
Julia doesnt have NULL. Thus, self-referencing indicates end of the list in either direction.
Initially, the aij pointers are assigned the value nothing. They are updated later by the add_aij functions.
"

# Creates row i with b[i] value b_val and adds it to Presolve_Problem p.
function add_row!(verbose::Bool, p::Presolve_Problem, i::Int, lb::Float64, ub::Float64)
    v = verbose
    v && println("Adding row : $(i)")

    row = Presolve_Row()
    row.i = i
    row.lb = lb
    row.ub = ub
    if(lb == ub)
        row.b_val = lb
    else
        row.b_val = nothing
    end
    row.aij = nothing
    row.is_active = false
    if(i!=1)
        row.prev = p.dictrow[i-1]
        row.prev.next = row
    else
        row.prev = row
        p.rowptr = row
    end
    row.next = row
    enque_row!(v,p,row)
    p.dictrow[i] = row
end

# places the specified row in the active list.
function enque_row!(verbose::Bool, p::Presolve_Problem, row::Presolve_Row)
    v = verbose
    v && println("Queueing row : $(row.i)")

    if(row.is_active == false)
        row.is_active = true
        row.active_prev = row.prev
        row.active_next = row
        if(p.rowque.i == -1)
            p.rowque = row
        end
        if(row.active_prev != row)
            row.active_prev.active_next = row
        end
    end
end

# removes the specified row from the active list.
function deque_row!(verbose::Bool, p::Presolve_Problem, row::Presolve_Row)
    v = verbose
    v && println("Dequeueing row : $(row.i)")

    if(row.is_active == true)
        row.is_active = false
        if(row.active_prev == row)
            p.rowque = row.active_next
        else
            row.active_prev.active_next = row.active_next
        end
        if(row.active_next != row)
            row.active_next.active_prev = row.active_prev
        end
    end
end

# Creates col j with c[i] value c_val and adds it to Presolve_Problem p.
function add_col!(verbose::Bool, p::Presolve_Problem, j::Int, c_val::Float64, l::Float64, u::Float64)
    v = verbose
    v && println("Adding col : $(j)")

    col = Presolve_Col()
    col.j = j
    col.c_val = c_val
    col.l = l
    col.u = u
    col.aij = nothing
    col.is_independent = false
    if(j!=1)
        col.prev = p.dictcol[j-1]
        col.prev.next = col
    else
        col.prev = col
        p.colptr = col
    end
    col.next = col
    enque_col!(v,p,col)
    p.dictcol[j] = col
end

# places the specified col in the independent list.
function enque_col!(verbose::Bool, p::Presolve_Problem, col::Presolve_Col)
    v = verbose
    v && println("Queueing col : $(col.j)")

    if(col.is_independent == false)
        col.ind_prev = col.prev
        col.is_independent = true
        col.ind_next = col
        if(p.colque.j == -1)
            p.colque = col
        end
        if(col.ind_prev != col)
            col.ind_prev.ind_next = col
        end
    end
end

# removes the specified col from the independent list.
function deque_col!(verbose::Bool, p::Presolve_Problem, col::Presolve_Col)
    v = verbose
    v && println("Dequeueing col : $(col.j)")

    if(col.is_independent == true)
        col.is_independent = false
        if(col.ind_prev == col)
            p.colque = col.ind_next
        else
            col.ind_prev.ind_next = col.ind_next
        end
        if(col.ind_next != col)
            col.ind_next.ind_prev = col.ind_prev
        end
    end
end

function remove_row!(verbose::Bool, p::Presolve_Problem, row::Presolve_Row)
    v = verbose
    deque_row!(v,p,row)

    while(row.aij != nothing)
        v && println("INSIDE------ and aij is $(row.aij.row.i),$(row.aij.col.j)")
        tmp = row.aij
        key = rc(tmp.row.i,tmp.col.j,p.originaln)
        #enque_col!(p,tmp.col)
        row.aij = tmp.row_next

        if(tmp.col_prev == tmp)
            if(tmp.col_next == tmp)
                tmp.col.aij = nothing
            else
                tmp.col.aij = tmp.col_next
                tmp.col_next.col_prev = tmp.col_next
            end
        else
            if(tmp.col_next == tmp)
                tmp.col_prev.col_next = tmp.col_prev
            else
                tmp.col_prev.col_next = tmp.col_next
            end
        end

        delete!(p.dictaij,key)

        if(row.aij == tmp)
            row.aij = nothing
        end
        tmp = nothing
    end

    if(row.prev == row)
        if(row.next != row)
            p.rowptr = row.next
            row.next.prev = row.next
        else
            p.rowptr = Presolve_Row()
        end
    else
        if(row.next != row)
            row.prev.next = row.next
            row.next.prev = row.prev
        else
            row.prev.next = row.prev
        end
    end

    delete!(p.dictrow,row.i)
end

function remove_col!(verbose::Bool, p::Presolve_Problem, col::Presolve_Col)
    v = verbose
    v && println("REMOVE COLUMN CALL")
    deque_col!(v,p,col)

    while(col.aij != nothing)
        tmp = col.aij
        key = rc(tmp.row.i,tmp.col.j,p.originaln)
        #enque_row!(p,tmp.row)
        col.aij = tmp.col_next

        if(tmp.row_prev == tmp)
            if(tmp.row_next == tmp)
                tmp.row.aij = nothing
            else
        #        println("here 1 ")
                tmp.row.aij = tmp.row_next
                tmp.row_next.row_prev = tmp.row_next
            end
        else
            if(tmp.row_next == tmp)
                tmp.row_prev.row_next = tmp.row_prev
            else
                tmp.row_prev.row_next = tmp.row_next
            end
        end

        delete!(p.dictaij,key)

        if(col.aij == tmp)
            col.aij = nothing
        end
        tmp = nothing
    end

    if(col.prev == col)
        if(col.next != col)
            p.colptr = col.next
            col.next.prev = col.next
        else
            p.colptr = Presolve_Col()
        end
    else
        if(col.next != col)
            col.prev.next = col.next
            col.next.prev = col.prev
        else
            col.prev.next = col.prev
        end
    end
    delete!(p.dictcol,col.j)
end

"
--- Matrix Element Operations ---
Matrix elements have only two functions - add_aij_normal and add_aij_transpose
They are never separately removed outside of removing rows or columns.
The input is a sparse matrix stored in the CSC format.
It can efficiently be accessed only in a column major order.
We need our matrix elements to act as two doubly linked lists.
One along the row and one along the column.
add_aij_normal creates the links along the column.
add_aij_transpose creates the links along the row.

We traverse the CSC matrix twice. First in the regular column major order and call add_aij_normal.
The second time we traverse the transpose(A) in column major order and call add_aij_transpose.
Note that the matrix element is already created by the time we call the second function.
"
function add_aij_normal!(verbose::Bool, p::Presolve_Problem, row_id::Int, col_id::Int, row_prev_id::Int, val::Float64)
    v = verbose
    v && println("Adding mat element (normal): $(row_id),$(col_id)")

    aij = Presolve_Matrix()
    aij.row = p.dictrow[row_id]
    aij.col = p.dictcol[col_id]
    p.rowcounter[row_id] += 1
    p.colcounter[col_id] += 1
    aij.val = val
    aij.row_prev = aij
    aij.row_next = aij
    aij.col_next = aij
    aij.col_prev = aij

    if(row_prev_id != -1)
        prev_key = rc(row_prev_id,col_id,p.originaln)
        p.dictaij[prev_key].col_next = aij
        aij.col_prev = p.dictaij[prev_key]
    end
    #can also be done by checking if row_prev == -1
    if(p.dictcol[col_id].aij == nothing)
        p.dictcol[col_id].aij = aij
    end

    key = rc(row_id,col_id,p.originaln)
    p.dictaij[key] = aij
end

function add_aij_transpose!(verbose::Bool, p::Presolve_Problem, row_id::Int, col_id::Int, col_prev_id::Int, val::Float64)
    v = verbose
    v && println("Adding mat element (transpose): $(row_id),$(col_id)")

    aij = p.dictaij[rc(row_id,col_id,p.originaln)]
    if(col_prev_id != -1)
        prev_key = rc(row_id,col_prev_id,p.originaln)
        p.dictaij[prev_key].row_next = aij
        aij.row_prev = p.dictaij[prev_key]
    end

    if(p.dictrow[row_id].aij == nothing)
        p.dictrow[row_id].aij = aij
    end
end

"
--- Presolve Setup and Cleanup ---
make_presolve       : Sets up the linked list connections.
print_info          : Prints all the linked list information. Useful for debugging.
make_new            : Converts the final presolve problem data into the format of the original problem"
function make_presolve!(verbose::Bool, p::Presolve_Problem, A::SparseMatrixCSC{Float64,Int64}, collb::Array{Float64,1},colub::Array{Float64,1}, c::Array{Float64,1}, rowlb::Array{Float64,1},rowub::Array{Float64,1})
    v = verbose
    m,n = size(A)

    # TODO.. make changes everywhere with collb,colub and rowlb,rowub

    # checks to ensure input problem is valid.
    p.originalm != m && error("Wrong size of b wrt A")
    p.originalm != length(rowlb) && error("Wrong size of rowlb wrt A")
    p.originalm != length(rowub) && error("Wrong size of rowub wrt A")
    p.originaln != n && error("Wrong size of c wrt A")
    p.originaln != length(collb) && error("Wrong size of collb wrt A")
    p.originaln != length(colub) && error("Wrong size of colub wrt A")

    v && println("Row SETUP ----- ")
    for i in 1:p.originalm
        add_row!(v,p,i,rowlb[i],rowub[i])
    end

    v && println("COL SETUP -----")
    for j in 1:p.originaln
        add_col!(v,p,j,c[j],collb[j],colub[j])
    end

    # Iterating through the non-zeros of sparse matrix A to construct the dictionary
    Arows = rowvals(A)
    v && println("MAT ELEMENT SETUP -----")
    Avals = nonzeros(A)
    for j = 1:p.originaln
        tmp = -1
        for i in nzrange(A,j)
            r = Arows[i]
            rcval = rc(r,j,p.originaln)
            #dictA[rcval] = Avals[i]
            #p.rowcounter[r] += 1
            #p.colcounter[j] += 1
            add_aij_normal!(v,p,r,j,tmp,Avals[i])
            tmp = r
        end
    end

    B = A'
    Arows = rowvals(B)
    Avals = nonzeros(B)
    for i = 1:p.originalm
        tmp = -1
        for c in nzrange(B,i)
            j = Arows[c]
            rcval = rc(i,j,p.originaln)
            add_aij_transpose!(v,p,i,j,tmp,Avals[c])
            tmp = j
        end
    end
end

# Constructs the reduced problem
function make_new(verbose::Bool, p::Presolve_Problem)
    v = verbose

    currentn = 0
    newc = Array{Float64,1}()
    newcollb = Array{Float64,1}()
    newcolub = Array{Float64,1}()

    col = p.colptr
    if(col.j != -1)
        v && println("Constructing newc,newlb,newub")
        while(col != nothing)
            v && @show col.j
            push!(newc,col.c_val)
            push!(newlb,col.l)
            push!(newub,col.u)

            currentn = currentn + 1
            p.finalcols[col.j] = currentn

            if(col.next == col)
                col = nothing
            else
                col = col.next
            end
        end
    end
    v && @show currentn

    currentm = 0
    newrowlb = Array{Float64,1}()
    newrowub = Array{Float64,1}()
    row = p.rowptr
    if(row.i != -1)
        v && println("Constructing newb")
        while(row != nothing)
            push!(newrowlb,row.lb)
            push!(newrowub,row.ub)
            currentm = currentm + 1
            p.finalrows[row.i] = currentm
            if(row.next == row)
                row = nothing
            else
                row = row.next
            end
        end
    end
    v && @show currentm

    v && println(p.finalcols)
    v && println(p.finalrows)

    I = Array{Int64,1}()
    J = Array{Int64,1}()
    Val = Array{Float64,1}()
    col = p.colptr
    if(col.j != -1)
        v && println("Constructing the new A matrix")
        while(col != nothing)
            tmp = col.aij
            while(tmp != nothing)
                v && @show tmp.row.i , tmp.col.j
                v && @show tmp.val
                push!(J,p.finalcols[tmp.col.j])
                push!(I,p.finalrows[tmp.row.i])
                push!(Val,tmp.val)
                if(tmp.col_next == tmp)
                    tmp = nothing
                else
                    tmp = tmp.col_next
                end
            end
            if(col.next == col)
                col = nothing
            else
                col = col.next
            end
        end
    end
    newA = sparse(I,J,Val,currentm,currentn)
    return newA,newcollb,newcolub,newc,newrowlb,newrowub
end

# For Debugging.
function print_info(p::Presolve_Problem)
    println("Row Information--------------------------------------")
    for key in keys(p.dictrow)
        println("-----------")
        row = p.dictrow[key]
        @show row.i
        @show row.lb
        @show row.ub
        @show row.b_val
        if(row.aij != nothing)
            @show row.aij.row.i
            @show row.aij.col.j
            @show row.aij.val
        end
        @show row.prev.i
        @show row.next.i
        @show row.is_active
        @show row.active_prev.i
        @show row.active_next.i
    end

    println("Col Information-------------------------------------------")
    for key in keys(p.dictcol)
        println("-----------")
        col = p.dictcol[key]
        @show col.j
        @show col.c_val
        @show col.l
        @show col.u
        if(col.aij != nothing)
            @show col.aij.row.i
            @show col.aij.col.j
            @show col.aij.val
        end
        @show col.prev.j
        @show col.next.j
        @show col.is_independent
        @show col.ind_prev.j
        @show col.ind_next.j
    end

    println("MAT ELEMENT Information---------------------------")
    for key in keys(p.dictaij)
        println("-----------")
        aij = p.dictaij[key]
        @show aij.row.i
        @show aij.col.j
        @show aij.val

        @show aij.row_prev.row.i, aij.row_prev.col.j, aij.row_prev.val
        @show aij.row_next.row.i, aij.row_next.col.j, aij.row_next.val
        @show aij.col_prev.row.i, aij.col_prev.col.j, aij.col_prev.val
        @show aij.col_next.row.i, aij.col_next.col.j, aij.col_next.val
    end
end


text:
using Test

include("../../src/utils/chordalextension.jl")

@testset "utils/chordalextension" begin
    #TODO
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
include("globals.jl")

module update
using globals

export update_electric_field!
export update_magnetic_field!

#
# Vacuum
#
function update_magnetic_field!(ez, hy)

    # Hack
    size = length(hy)
    if length(hy) == length(ez)
        size -= 1
    end
    #
    #
    for i = 1:size
        hy[i] = hy[i] + (ez[i + 1] - ez[i]) / globals.imp0
    end
end
function update_electric_field!(ez, hy)
    for i = 2:length(ez)
        ez[i] = ez[i] + (hy[i] - hy[i-1]) * globals.imp0
    end
end

#
# Material
#

function update_magnetic_field!(ez, hy, mu)
    for i = 1:length(hy)
        hy[i] = hy[i] + (ez[i + 1] - ez[i]) / globals.imp0 / mu[i]
    end
end
function update_electric_field!(ez, hy, eps)
    for i = 2:length(ez)-1
        ez[i] = ez[i] + (hy[i] - hy[i-1]) * globals.imp0 / eps[i]
    end
end

#
# Lossy
#

function update_magnetic_field!(ez, hy, mu, chyh, chye)
    for i = 1:length(hy)
        hy[i] = chyh[i] * hy[i] + chye[i] * (ez[i + 1] - ez[i]) / globals.imp0 / mu[i]
    end
end

function update_electric_field!(ez, hy, eps, cezh, ceze)
    for i = 2:length(ez)-1
        ez[i] = ceze[i] * ez[i] + cezh[i] * (hy[i] - hy[i-1]) * globals.imp0 / eps[i]
    end
end

end


text:
"""
    \\lead{text}

A simple lead class.
"""
lx_lead(com, _) = "@@lead $(proc(com)) @@"

"""
    \\figure{path=..., alt=..., width=..., ...}

Insert a figure with specifications.
"""
function lx_figure(com, _)
    content = strip(proc(com))
    _, kwargs = lxargs(content, "figure")
    return _figure(; kwargs...)
end

function _figure(; path="", alt="", width="", style="", caption="", imgclass="", captionclass="")
    style = ifelse(isempty(style), "", "style=\"$style\"")
    if !isempty(caption)
        isempty(alt) && (alt = caption)
        caption = """

            $caption

            """
    end
    return html("""

          $caption

        """)
end

"""
    \\table{..., opt1=...}

Insert a table with specifications.
"""
function lx_table(com, _)
    content = strip(proc(com))
    args, kwargs = lxargs(content, "table")
    length(args) == 1 ||
        error("Expected a single un-named argument for `\\table`.")
    return _table(args[1]; kwargs...)
end

function _table(md; style="", caption="", class="")
    caption = ifelse(isempty(caption), "", "
$caption
")
    parser = CommonMark.Parser()
    CommonMark.enable!(parser, CommonMark.TableRule())
    h = CommonMark.html(parser(md))
    h = replace(h, "
"=>"")
    return html("""
$caption $h""")
end

function lx_alert(com, _)
    content = strip(proc(com))
    return html("""

        """) * content * html("
")
end


text:
query = """
{
  neomatrix{
    nombre
    linkedin
  }
}
"""

r = Queryclient("https://neomatrix.herokuapp.com/graphql",query,headers=Dict("header"=>"value"))

@test r.Info.status == 200
@test r.Data == "{\"data\":{\"neomatrix\":{\"nombre\":\"Acevedo Maldonado Josue\",\"linkedin\":\"https://www.linkedin.com/in/acevedo-maldonado-josue/\"}}}"

#=client = GraphQLClient("https://api.graph.cool/simple/v1/movies",auth="Bearer my-jwt-token",headers=Dict("header"=>"value"))

query2 = """
query getMovie(\$title: String!) {
  Movie(title:\$title) {
    releaseDate
    actors {
      name
    }
  }
}
"""
r = client.Query(query2,vars=Dict("title" => "Inception"))

@test r.Info.status == 200
@test r.Data == "{\"data\":{\"Movie\":{\"releaseDate\":\"2010-08-28T20:00:00.000Z\",\"actors\":[{\"name\":\"Leonardo DiCaprio\"},{\"name\":\"Ellen Page\"},{\"name\":\"Tom Hardy\"},{\"name\":\"Joseph Gordon-Levitt\"},{\"name\":\"Marion Cotillard\"}]}}}"

r = client.Query(query2,check=true)
@test r == "ok"

query2 = """
{
  Movie(title: "Inception"){
    actors{
      name
    }
  }
}
"""
r = client.Query(query2)

@test r.Info.status == 200
@test r.Data == "{\"data\":{\"Movie\":{\"actors\":[{\"name\":\"Leonardo DiCaprio\"},{\"name\":\"Ellen Page\"},{\"name\":\"Tom Hardy\"},{\"name\":\"Joseph Gordon-Levitt\"},{\"name\":\"Marion Cotillard\"}]}}}"

r = client.Query(query2,getlink=true)

@test r == "https://api.graph.cool/simple/v1/movies?query=%7B%0A%20%20Movie%28title%3A%20%22Inception%22%29%7B%0A%20%20%20%20actors%7B%0A%20%20%20%20%20%20name%0A%20%20%20%20%7D%0A%20%20%7D%0A%7D%0A"

@test client.headers(Dict("header"=>"value")) == Dict("header"=>"value")

@test client.serverAuth("Bearer my-jwt-token") == "Bearer my-jwt-token"

=#

query="https://neomatrix.herokuapp.com/graphql?query=%7B%0A%20%20neomatrix%7B%0A%20%20%20%20nombre%0A%20%20%20%20linkedin%0A%20%20%7D%0A%7D"
r = Queryclient(query)
@test r.Info.status == 200
@test r.Data == "{\"data\":{\"neomatrix\":{\"nombre\":\"Acevedo Maldonado Josue\",\"linkedin\":\"https://www.linkedin.com/in/acevedo-maldonado-josue/\"}}}"

query = """
{
  neomatrix{
    nombre
    linkedin
  }
}
"""

r = Queryclient("https://neomatrix.herokuapp.com/graphql",query,getlink=true)
@test r == "https://neomatrix.herokuapp.com/graphql?query=%7B%0A%20%20neomatrix%7B%0A%20%20%20%20nombre%0A%20%20%20%20linkedin%0A%20%20%7D%0A%7D%0A"

query = """
       query consulta{
          neomatrix{
            nombre
            linkedin
          }
        }

       query hola{
          neomatrix{
            nombre
          }
       }
       """
r = Queryclient("https://neomatrix.herokuapp.com/graphql",query,operationName="hola")
@test r.Data == "{\"data\":{\"neomatrix\":{\"nombre\":\"Acevedo Maldonado Josue\"}}}"

#=
query = """
query consulta{
  Movie(title: "Inception"){
    actors{
      name
    }
  }
}
query hola{
  Movie(title: "Inception"){
    actors{
      name
    }
  }
}
"""
r = client.Query(query,operationName="consulta")

@test r.Data == "{\"data\":{\"Movie\":{\"actors\":[{\"name\":\"Leonardo DiCaprio\"},{\"name\":\"Ellen Page\"},{\"name\":\"Tom Hardy\"},{\"name\":\"Joseph Gordon-Levitt\"},{\"name\":\"Marion Cotillard\"}]}}}"
=#
query = """
query consulta{
  neomatrix{
    nombre
  }
}
query hola{
  neomatrix{
    nombre
    linkedin
  }
}
"""
r = Queryclient("https://neomatrix.herokuapp.com/graphql",query,getlink=true,operationName="consulta")

@test r == "https://neomatrix.herokuapp.com/graphql?query=query%20consulta%7B%0A%20%20neomatrix%7B%0A%20%20%20%20nombre%0A%20%20%7D%0A%7D%0Aquery%20hola%7B%0A%20%20neomatrix%7B%0A%20%20%20%20nombre%0A%20%20%20%20linkedin%0A%20%20%7D%0A%7D%0A&operationName=consulta"


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

"""
    FullTopology(connectivities)

A data structure that stores *all* `connectivities` of a mesh.

### Notes

This data structure is sometimes referred to as the "soup of geometries".
It does *not* support topological relations and is therefore incompatible
with algorithms that rely on neighborhood search. It is still useful for
mesh visualization and IO operations.
"""
struct FullTopology{C<:Connectivity} <: Topology
  # input fields
  connec::Vector{C}

  # state fields
  ranks::Vector{Int}
  elms::Vector{Int}

  function FullTopology{C}(connec) where {C}
    ranks = [paramdim(c) for c in connec]
    elms  = findall(isequal(maximum(ranks)), ranks)
    new(connec, ranks, elms)
  end
end

FullTopology(connec) = FullTopology{eltype(connec)}(connec)

==(t1::FullTopology, t2::FullTopology) = t1.connec == t2.connec

# ---------------------
# HIGH-LEVEL INTERFACE
# ---------------------

nvertices(t::FullTopology) = maximum(i for c in t.connec for i in indices(c))

function faces(t::FullTopology, rank)
  cs, rs = t.connec, t.ranks
  (cs[i] for i in 1:length(cs) if paramdim(cs[i]) == rank)
end

element(t::FullTopology, ind) = t.connec[t.elms[ind]]

nelements(t::FullTopology) = length(t.elms)

facets(t::FullTopology) = faces(t, maximum(t.ranks) - 1)

nfacets(t::FullTopology) = count(==(maximum(t.ranks) - 1), t.ranks)

# ------------
# CONVERSIONS
# ------------

function Base.convert(::Type{FullTopology}, t::Topology)
  # TODO: add all faces, not just the elements
  FullTopology(collect(elements(t)))
end


text:
function ABCDQR_boots_s(y,u,A,B,C,D,Q,R,x10,nb)
	#
	# computes nb STATIONARY bootstrap replications for state space model
	#
	# x_{t+1} = A*x_{t} + B*u_{t} + w_{t}
	# y_{t}   = C*x_{t} + D*u_{t} + v_{t}
	#
	# javier.cara@upm.es, 2016-04

	ny,nt = size(y)
	nx = size(A,1)

	# step 0: compute the innovations
	# -----------------------------------------------------
	xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = ABCDQR_kfilter_s(y,u,A,B,C,D,Q,R,x10)

	# step 1: standardized innovations
	# -----------------------------------------------------
	# St^(1/2)
	Std,V = eig(St)
	St_05 = V*diagm(sqrt(Std))*V' # St is symmetric = inv(V)=V'
	# St^(-1/2)*et
	et_st = St_05\et

	# bootstrap replications
	ybo = zeros(ny,nt,nb)
	xtt1_bo = zeros(nx,nt)
	for i in 1:nb
		# step 2: sampling et_st WITH replacement
   	# ------------------------------------------------------
   	orden = rand(1:nt,nt)
   	et_st_bo = et_st[:,orden]

   	# step 3: bootstrap data set
   	# ------------------------------------------------------
		for t in 1:nt-1
			xtt1_bo[:,t+1] = A*xtt1_bo[:,t] + B*u[:,t] + A*Kt*St_05*et_st_bo[:,t]
		end
		ybo[:,:,i] = C*xtt1_bo + D*u + St_05*et_st_bo

	end # nb

	return ybo

end


text:
# This package deals with complex spherical coordinates and complex spherical harmonics. Many of the identities used are not well known, so we pay special attention to tests for these complex functions.

using Test, LinearAlgebra

@testset "Tests for complex transmission direction" begin

    k_eff = rand(-1:0.11:1.0) + rand(-1:0.11:1.0) * im
    incident_wavevector = rand(-1:0.11:1.0,3) + rand(-1:0.11:1.0,3) .* im
    surface_normal = rand(-1:0.11:1.0,3)
    surface_normal = surface_normal / norm(surface_normal)

    # make sure this surface_normal is the outward normal
    if real(dot(surface_normal, incident_wavevector)) > 0
        surface_normal = - surface_normal
    end

    normal_eff = transmission_direction(k_eff, incident_wavevector, surface_normal)

    @test sum(normal_eff .^2) ≈ 1.0 + 0.0im
    @test imag(k_eff * dot(-surface_normal, normal_eff)) > 0

    k_eff = rand(-1:0.1:1.0) + rand(-1:0.1:1.0) * im
    incident_wavevector = rand(-1:0.1:1.0,2) + rand(-1:0.1:1.0,2) .* im
    surface_normal = rand(-1:0.1:1.0,2)
    surface_normal = surface_normal / norm(surface_normal)

    if real(dot(surface_normal, incident_wavevector)) > 0
        surface_normal = - surface_normal
    end

    normal_eff = transmission_direction(k_eff, incident_wavevector, surface_normal)

    @test sum(normal_eff .^2) ≈ 1.0 + 0.0im
    @test imag(k_eff * dot(-surface_normal, normal_eff)) > 0

    k_eff = rand(-1:0.1:1.0) + rand() * im
    k = rand()
    θin = rand(-pi/2.0:0.2:pi/2.0)
    θeff = transmission_angle_wiener(k, k_eff, θin)

    normal2 = [cos(θeff+pi/2.0),sin(θeff+pi/2.0)]

    surface_normal = [0.0,-1.0];
    incident_wavevector = k*[cos(θin+pi/2.0),sin(θin+pi/2.0)]
    normal_eff = transmission_direction(k_eff, incident_wavevector, surface_normal)

    @test norm(normal_eff - normal2) < 1e-10

end

@testset "Tests for complex transmission angle" begin

    surface_normal = [0.0,0.0,-1.0]
    xs = [rand(-1.01:0.1:1.0,3) + rand(-1.01:0.1:1.0,3)*im for i = 1:100];

    # make sure these are incident upon the surface
    xs = [real(dot(surface_normal, x)) > 0 ? -x : x for x in xs];

    @test all(cartesian_to_radial_coordinates(x)[2] == transmission_angle(x, surface_normal) for x in xs)

    surface_normal = [-1.0,0.0]
    xs = [rand(-1.01:0.1:1.0,2) + rand(-1.01:0.1:1.0,2)*im for i = 1:100];
    xs = [real(dot(surface_normal, x)) > 0 ? -x : x for x in xs];

    @test all(cartesian_to_radial_coordinates(x)[2] == transmission_angle(x, surface_normal) for x in xs)

    N = 100
    is = 1:N
    ns = [rand(-1.01:0.1:1.0,2) for i = is];
    ns = ns ./ norm.(ns);

    v_ins = [rand(-1.01:0.1:1.0,2) for i = is]
    v_ins = [real(dot(ns[i], v_ins[i])) > 0 ? -v_ins[i] : v_ins[i] for i in is]

    θ_ins = transmission_angle.(v_ins, ns)

    k_effs = rand(-1.01:0.1:1.0,N) + rand(-1.01:0.1:1.0,N) .* im
    v_effs = k_effs .* transmission_direction.(k_effs, v_ins, ns)
    θ_effs = transmission_angle.(v_effs, ns)

    # Check snells law, i.e. that the components orthogonal to the surface are the same for the two wavevectors
    @test maximum(norm(v_effs[i] - v_ins[i] - ns[i] .* dot(ns[i], v_effs[i] - v_ins[i])) for i in is) < 1e-14
    @test maximum(abs.((k_effs  .* sin.(θ_effs)) .^2 - (norm.(v_ins)  .* sin.(θ_ins)) .^2)) < 1e-10

    # Check that the effective wave attenuates as it moves into the material
    @test all(imag.(dot.(- ns,v_effs)) .> 0)

    # NOTE sum(v_effs[i] .^2) ≈ k_effs[i] ^2, however sqrt(sum(v_effs[i] .^2)) is not necessarily equal to k_effs[i], as the direction of v_effs[i] is chosen so that the wave attenuates when transmitted into the material.

    # When real part of dot(-n,v) is positive then the wave propagates into the material, in which case if real(k_effs) > 0 then we can recover k_effs from the v_effs.
    ns = [ -rand(0.2:0.1:1.0,2) for i = is]
    ns = ns ./ norm.(ns)

    v_ins = [rand(0.2:0.1:1.0,2) for i = is]
    θ_ins = transmission_angle.(v_ins, ns)

    @test all(-pi/2 .< real.(θ_ins) .< pi/2)

    # choose positive real part
    k_effs = rand(0.01:0.1:1.0,N) + rand(-1.1:0.1:1.0,N) .* im;
    v_effs = k_effs .* transmission_direction.(k_effs, v_ins, ns);
    θ_effs = transmission_angle.(v_effs, ns);

    @test maximum(abs.(k_effs  .* sin.(θ_effs) - norm.(v_ins)  .* sin.(θ_ins)) ) < 1e-10

    #angle from x-axis
    θ_ns = [transmission_angle(-n, [-1.0,0.0]) for n in ns];

    @test maximum(norm(v_ins[i] - norm(v_ins[i]) .* [cos(θ_ins[i] + θ_ns[i]), sin(θ_ins[i] + θ_ns[i])]) for i in is) < 1e-10
    @test maximum(norm(v_effs[i] - k_effs[i] .* [cos(θ_effs[i] + θ_ns[i]), sin(θ_effs[i] + θ_ns[i])]) for i in is) < 1e-10

# test for 3-dimensions

    ns = [rand(-1.01:0.1:1.0,3) for i = is];
    ns = ns[is] ./ norm.(ns[is]);

    v_ins = [rand(-1.01:0.1:1.0,3) for i = is];
    v_ins = [real(dot(ns[i], v_ins[i])) > 0 ? -v_ins[i] : v_ins[i] for i in is];

    θ_ins = transmission_angle.(v_ins, ns);

    v_effs = k_effs .* transmission_direction.(k_effs, v_ins, ns);
    θ_effs = transmission_angle.(v_effs, ns);

    # Check snells law, i.e. that the components orthogonal to the surface are the same for the two wavevectors
    @test maximum(norm(v_effs[i] - v_ins[i] - ns[i] .* dot(ns[i], v_effs[i] - v_ins[i])) for i in is) < 1e-14
    @test maximum(abs.((k_effs  .* sin.(θ_effs)) .^2 - (norm.(v_ins)  .* sin.(θ_ins)) .^2)) < 1e-10

    # Check that the effective wave attenuates as it moves into the material
    @test all(imag.(dot.(- ns,v_effs)) .>= 0.0)
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
using Pkg
const HOME = ENV["HOME"]
using Printf
Pkg.activate("$HOME/HelioseismicKernels")

using MPI
MPI.Init()
using HelioseismicKernels
const comm = MPI.COMM_WORLD

const xobs1 = Point2D(pi/2, pi/3)
const xobs2 = Point2D(pi/2, pi/2)
const ℓ_range = 1:100

try
    HelioseismicKernels.kernel_uₛₜ(comm, TravelTimes(), xobs1, xobs2, los_radial(), save = false, ℓ_range = 1:1)
    HelioseismicKernels.kernel_uₛₜ(comm, TravelTimes(), xobs1, xobs2, los_earth(), save = false, ℓ_range = 1:1)

    s_max = 25

    tstart = time()
    HelioseismicKernels.kernel_uₛₜ(comm,
        TravelTimes(), xobs1, xobs2, los_earth(),
        s_max = s_max, t_max = s_max, ℓ_range = ℓ_range)

    tend = time()
    Δt = round(Int, tend - tstart)
    if MPI.Comm_rank(comm) == 0
        @printf "los\t%d\t%d\n" s_max Δt
    end

    tstart = time()

    HelioseismicKernels.kernel_uₛₜ(comm,
        TravelTimes(), xobs1, xobs2, los_radial(),
        s_max = s_max, t_max = s_max, ℓ_range = ℓ_range)

    tend = time()
    Δt = round(Int, tend - tstart)
    if MPI.Comm_rank(comm) == 0
        @printf "rad\t%d\t%d\n" s_max Δt
    end

finally
    MPI.Finalize()
end


text:
using Test, DataLoaders

include("dataset.jl")

@testset "DataLoder" begin
  batchsize = 16
  dataset = getdataset((10,10), 100)
  loader  = DataLoader(dataset, batchsize, ntasks = 8, shuffle=true, droplast=false)

  @test length(loader) == round(Int, 100 / 16, RoundDown)

  (x, y), index = iterate(loader)

  @test size(x) == (10, 10, batchsize)
  @test size(y) == (batchsize,)
  @test index   == batchsize + 1
end


text:
################################################################################
#
#  Lattice database
#
################################################################################

export number_of_lattices, lattice_name, lattice,
       lattice_automorphism_group_order, lattice_database

struct LatticeDB
  path::String
  max_rank::Int
  db::Vector{Vector{NamedTuple{(:name, :rank, :deg, :amb, :basis_mat, :min, :aut, :kissing),
                               Tuple{String, Int, Int, Vector{Rational{BigInt}}, Vector{Rational{BigInt}}, BigInt, BigInt, BigInt}}}}

  function LatticeDB(path::String)
    db = Meta.eval(Meta.parse(Base.read(path, String)))
    max_rank = length(db)
    return new(path, max_rank, db)
  end
end

# TODO: Write a parser for the data

function show(io::IO, L::LatticeDB)
  print(io, "Nebe-Sloan database of lattices (rank limit = ", L.max_rank, ")")
end

const default_lattice_db = Ref(joinpath(pkgdir, "data/lattices"))

################################################################################
#
#  For creating a lattice database
#
################################################################################

function lattice_database()
  if !isfile(joinpath(pkgdir, "data/lattices"))
    download_data(data = "quadratic_lattices")
  end
  return LatticeDB(default_lattice_db[])
end

function lattice_database(path::String)
  return LatticeDB(path)
end

################################################################################
#
#  Conversion from linear indicies
#
################################################################################

function from_linear_index(L::LatticeDB, i::Int)
  k = 1
  while i > length(L.db[k])
    i = i - length(L.db[k])
    k += 1
  end
  return (k, i)
end

################################################################################
#
#  Out of bounds error functions
#
################################################################################

@inline function _check_rank_range(L, r)
  r < 0 || r > L.max_rank &&
        throw(error("Rank ($(r)) must be between 1 and $(L.max_rank)"))
end

@inline function _check_range(L, r, i)
  r < 0 || r > L.max_rank &&
          throw(error("Rank ($(r)) must be between 1 and $(L.max_rank)"))
  j = number_of_lattices(L, r)
  i < 0 || i > j && throw(error("Index ($(i)) must be between 1 and $(j)"))
end

@inline function _check_range(L, i)
  j = number_of_lattices(L)
  i < 0 || i > j && throw(error("Index ($(i)) must be between 1 and $(j)"))
end

################################################################################
#
#  Access functions
#
################################################################################

function number_of_lattices(L::LatticeDB, r::Int)
  _check_rank_range(L, r)
  return length(L.db[r])
end

function number_of_lattices(L::LatticeDB)
  return sum(length.(L.db))
end

function lattice_name(L::LatticeDB, r::Int, i::Int)
  _check_range(L, r, i)
  return L.db[r][i].name
end

function lattice_name(L::LatticeDB, i::Int)
  _check_range(L, i)
  return lattice_name(L, from_linear_index(L, i)...)
end

function lattice_automorphism_group_order(L::LatticeDB, r::Int, i::Int)
  _check_range(L, r, i)
  return L.db[r][i].aut
end

function lattice_automorphism_group_order(L::LatticeDB, i::Int)
  _check_range(L, i)
  return lattice_automorphism_group_order(L, from_linear_index(L, i)...)
end

function lattice(L::LatticeDB, r::Int, i::Int)
  _check_range(L, r, i)
  d = L.db[r][i].deg
  A = matrix(FlintQQ, d, d, L.db[r][i].amb)
  B = matrix(FlintQQ, r, d, L.db[r][i].basis_mat)
  return Zlattice(B, gram = A)
end

function lattice(L::LatticeDB, i::Int)
  _check_range(L, i)
  return lattice(L, from_linear_index(L, i)...)
end

################################################################################
#
#  Quadratic lattices DB
#
################################################################################

const default_quad_lattice_db = Ref(joinpath(pkgdir, "data/quadratic_lattices"))

struct QuadLatDB
  path::String
  #db::Vector{Tuple{Vector{BigInt}, Vector{Vector{Rational{BigInt}}}, Vector{Vector{Rational{BigInt}}}, Int}}
  metadata
  head::Int
  length::Int

  function QuadLatDB(path::String)
    metadata = Dict()
    f = open(path)
    head = 0
    while true
      line = readline(f)
      if line[1] == '#'
        head += 1
        if !(':' in line)
          continue
        end
        i = 2
        while line[i] != ':'
          i += 1
        end
        key = strip(line[2:i-1])
        val = strip(line[i+1:end])
        metadata[key] = val
      else
        break
      end
    end
    seekstart(f)
    length = countlines(f)
    close(f)
    return new(path, metadata, head, length - head)
  end
end

function quadratic_lattice_database()
  if !isfile(joinpath(pkgdir, "data/quadratic_lattices"))
    download_data(data = "quadratic_lattices")
  end
  return QuadLatDB(default_quad_lattice_db[])
end

function quadratic_lattice_database(path::String)
  return QuadLatDB(path)
end

Base.length(D::QuadLatDB) = D.length

class_number(D::QuadLatDB, i::Int) = _lattice_data(D, i)[4]

function Base.show(io::IO, D::QuadLatDB)
  s = get(D.metadata, "Description", "Quadratic lattices database")
  print(io, s, "\n")
  if haskey(D.metadata, "Author")
    print(io, "Author: $(D.metadata["Author"])\n")
  end
  if haskey(D.metadata, "Source")
    print(io, "Source: $(D.metadata["Source"])\n")
  end
  if haskey(D.metadata, "Version")
    print(io, "Version: $(VersionNumber(D.metadata["Version"]))\n")
  end

  print(io, "Number of lattices: ", D.length)
end

function versioninfo(D::QuadLatDB)
  return VersionNumber(D.metadata["Version"])
end

function _lattice_data(D::QuadLatDB, i::Int)
  it = Iterators.drop(eachline(D.path), D.head + i - 1)
  line = IOBuffer(first(it))
  return _parse_quad(line, versioninfo(D))
end

function lattice(D::QuadLatDB, i::Int)
  return _get_lattice(_lattice_data(D, i))
end

function _get_lattice(data)
  f = Globals.Qx(data[1])
  K, a = NumberField(f, "a", cached = false)
  diag = map(K, data[2])
  gens = map(K, data[3])
  D = diagonal_matrix(diag)
  n = nrows(D)
  @assert iszero(mod(length(gens), n))
  gens_split = collect(Iterators.partition(gens, n))
  return quadratic_lattice(K, generators = gens_split, gram_ambient_space = D)
end

################################################################################
#
#  Hermitian lattices DB
#
################################################################################

const default_herm_lattice_db = Ref(joinpath(pkgdir, "data/hermitian_lattices"))

struct HermLatDB
  path::String
  #db::Vector{Tuple{Vector{BigInt}, Vector{Vector{Rational{BigInt}}}, Vector{Vector{Rational{BigInt}}}, Int}}
  metadata
  head::Int
  length::Int

  function HermLatDB(path::String)
    metadata = Dict()
    f = open(path)
    head = 0
    while true
      line = readline(f)
      if line[1] == '#'
        head += 1
        if !(':' in line)
          continue
        end
        i = 2
        while line[i] != ':'
          i += 1
        end
        key = strip(line[2:i-1])
        val = strip(line[i+1:end])
        metadata[key] = val
      else
        break
      end
    end
    seekstart(f)
    length = countlines(f)
    close(f)
    return new(path, metadata, head, length - head)
  end
end

function hermitian_lattice_database()
  if !isfile(joinpath(pkgdir, "data/hermitian_lattices"))
    download_data(data = "hermitian_lattices")
  end
  return HermLatDB(default_herm_lattice_db[])
end

function hermitian_lattice_database(path::String)
  return HermLatDB(path)
end

Base.length(D::HermLatDB) = D.length

class_number(D::HermLatDB, i::Int) = _lattice_data(D, i)[5]

function Base.show(io::IO, D::HermLatDB)
  s = get(D.metadata, "Description", "Hermratic lattices database")
  print(io, s, "\n")
  if haskey(D.metadata, "Author")
    print(io, "Author: $(D.metadata["Author"])\n")
  end
  if haskey(D.metadata, "Source")
    print(io, "Source: $(D.metadata["Source"])\n")
  end
  if haskey(D.metadata, "Version")
    print(io, "Version: $(VersionNumber(D.metadata["Version"]))\n")
  end

  print(io, "Number of lattices: ", D.length)
end

function versioninfo(D::HermLatDB)
  return VersionNumber(D.metadata["Version"])
end

function _lattice_data(D::HermLatDB, i::Int)
  it = Iterators.drop(eachline(D.path), D.head + i - 1)
  line = IOBuffer(first(it))
  return _parse_herm(line, versioninfo(D))
end

function lattice(D::HermLatDB, i::Int)
  return _get_hermitian_lattice(_lattice_data(D, i))
end

function _get_hermitian_lattice(data)
  f = Globals.Qx(data[1])
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t", cached = false)
  E, b = NumberField(Kt(map(K, data[2])), "b", cached = false)
  diag = map(E, map(K, data[3]))
  k = degree(K)
  #gens = [ E(map(K, collect(Iterators.partition(v, k))))  map(E, data[4])
  D = diagonal_matrix(diag)
  n = nrows(D)
  @assert iszero(mod(length(data[4]), n))
  gens_split = collect(Iterators.partition(data[4], n))
  gens = []
  for v in gens_split
    push!(gens, map(E, [map(K, collect(Vector.(Iterators.partition(w, k)))) for w in v]))
  end
  return hermitian_lattice(E, generators = gens, gram_ambient_space = D)
end

################################################################################
#
#  Parse
#
################################################################################

base_type(::Type{fmpq}) = fmpz

base_type(::Type{Rational{T}}) where {T} = T

function parse_int(io, delim = UInt8[UInt8(',')])
  n = IOBuffer()
  b = Base.read(io, UInt8)
  while !(b in delim)
    Base.write(n, b)
    if eof(io)
      break
    end
    b = Base.read(io, UInt8)
  end

  return b, parse(Int, String(take!(n)))
end

function parse_rational(::Type{T}, io, delim = [UInt8(',')]) where {T}
  n = IOBuffer()
  d = IOBuffer()
  read_den = false
  b = Base.read(io, UInt8)
  while !(b in delim)
    if b == UInt8('/')
      read_den = true
      skip(io, 1)
      @goto read
    end

    if !read_den
      Base.write(n, b)
    else
      Base.write(d, b)
    end

    @label read
    if eof(io)
      break
    end
    b = Base.read(io, UInt8)
  end

  nu = String(take!(n))

  if length(nu) < 18
    num = base_type(T)(parse(Int, nu))
  else
    num = parse(base_type(T), nu)
  end

  if !read_den
    return b, T(num)
  else
    de = String(take!(d))
    if length(de) < 18
      den = base_type(T)(parse(Int, de))
    else
      den = parse(base_type(T), de)
    end
    return b, num//den
  end
end

function parse_array(::Type{T}, io, ldelim = '[', rdelim = ']') where {T}
  res = Vector{T}()
  b = Base.read(io, UInt8)
  @assert b == UInt8(ldelim)
  mark(io)
  b = Base.read(io, UInt8)
  if b == UInt8(rdelim)
    return b, res
  else
    reset(io)
    unmark(io)
    while true
      b, z = parse_rational(T, io, map(UInt8, [',', rdelim]))
      push!(res, z)
      if b != UInt8(',')
        @assert b == UInt8(rdelim)
        break
      end
    end
    return b, res
  end
end

function parse_array(::Type{Vector{T}}, io, ldelim = '[', rdelim = ']') where {T}
  res = Vector{Vector{T}}()
  b = Base.read(io, UInt8)
  @assert b == UInt8(ldelim)
  mark(io)
  b = Base.read(io, UInt8)
  if b == rdelim
    return res
  else
    reset(io)
    unmark(io)
    while true
      b, z = parse_array(T, io)
      push!(res, z)
      b = Base.read(io, UInt8)
      if b == UInt8(rdelim)
        break
      end
      if b != UInt8(',')
        @assert b == UInt8(rdelim)
        break
      end
    end
    return res
  end
end

#db::Vector{Tuple{Vector{BigInt}, Vector{Vector{Rational{BigInt}}}, Vector{Vector{Rational{BigInt}}}, Int}}
function _parse_quad(io, version)
  @assert version == v"0.0.1"
  b, def_poly = parse_array(fmpq, io)
  @assert b == UInt8(']')
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  diagonal = parse_array(Vector{fmpq}, io)
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  gens = parse_array(Vector{fmpq}, io)
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  b, cl = parse_int(io)
  return def_poly, diagonal, gens, cl
end

function _parse_herm(io, version)
  @assert version == v"0.0.1"
  b, def_poly = parse_array(fmpq, io)
  @assert b == UInt8(']')
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  ext_poly = parse_array(Vector{fmpq}, io)
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  diagonal = parse_array(Vector{fmpq}, io)
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  gens = parse_array(Vector{fmpq}, io)
  b = Base.read(io, UInt8)
  @assert b == UInt8(',')
  b, cl = parse_int(io)
  return def_poly, ext_poly, diagonal, gens, cl
end


text:
module StringBuilders

export StringBuilder

mutable struct StringBuilder
    buffer::IOBuffer
    as_string::Union{String,Nothing}

    function StringBuilder()
        return new(IOBuffer(), nothing)
    end
end

function Base.String(sb::StringBuilder)

    if sb.as_string === nothing
        sb.as_string = String(take!(sb.buffer))
    end

    return sb.as_string
end

function Base.append!(sb::StringBuilder, s::AbstractString)
    if sb.as_string !== nothing
        print(sb.buffer, sb.as_string)
        sb.as_string = nothing
    end
    print(sb.buffer, s)
end

function Base.append!(sb::StringBuilder, ss::AbstractString...)
    for s in ss
        append!(sb, s)
    end
end

end # module


text:
using Interpolations
using StaticArrays
using LinearAlgebra
using Statistics
using QuadGK
using Plots
# using Makie

function Rot2(θ)
    sθ,cθ = sincos(θ)
    @SMatrix [cθ -sθ; sθ cθ]
end

abstract type AbstractPath end

struct CarPath{T} <: AbstractPath
    s::Vector{T}
    ϕ::Vector{T}
    κ::Vector{T}
    X::Vector{T}
    Y::Vector{T}
end

function curvature(path::CarPath, s)
    k = searchsortedfirst(path.s, s)
    return path.κ[k]
end

struct CirclePath <: AbstractPath
    r::Float64
end
curvature(path::CirclePath, s) = 1/path.r

""" ArcPath """
struct ArcPath <: AbstractPath
    ang0::Float64
    radius::Float64
    ang::Float64
    xc::Float64
    yc::Float64
    function ArcPath(ang0, radius, ang)
        xc = -radius*sin(ang0)
        yc =  radius*cos(ang0)
        new(ang0, radius, ang, xc, yc)
    end
end

function position_change(arc::ArcPath)
    localToGlobal(arc, total_length(arc), 0)
end
total_length(arc::ArcPath) = arc.radius * arc.ang
curvature(arc::ArcPath, s) = 1/arc.radius
final_angle(arc::ArcPath) = arc.ang0 + arc.ang
function localToGlobal(arc::ArcPath, s, e)
    r = arc.radius
    θ = s ./ r
    radius = r .- e
    # radius = r
    ϕ = θ .- π/2 .+ arc.ang0
    x = radius .* cos.(ϕ) .+ arc.xc
    y = radius .* sin.(ϕ) .+ arc.yc
    return x,y
end

""" Straight Path """
struct StraightPath <: AbstractPath
    len::Float64
    ang::Float64
end

function position_change(path::StraightPath)
    x = path.len * cos(path.ang)
    y = path.len * sin(path.ang)
    return x,y
end
total_length(path::StraightPath) = path.len
final_angle(path::StraightPath) = path.ang
curvature(path::StraightPath, s) = 0.0
function localToGlobal(path::StraightPath, s, e)
    Rot2(path.ang)*@SVector [s, e]
end

"Build a straight path that continues from an arc"
function StraightPath(arc::ArcPath, len)
    StraightPath(len, final_angle(arc))
end

"Build an arc that continues from a line"
function ArcPath(line::StraightPath, radius, ang)
    ArcPath(final_angle(line), radius, ang)
end

""" Dubins Path """
struct DubinsPath{P} <: AbstractPath
    paths::P
    lengths::Vector{Float64}
    ds::Vector{Float64}
    dx::Vector{Float64}
    dy::Vector{Float64}
    curvature::Vector{Float64}
end

function DubinsPath(paths)
    num_paths = length(paths)
    lengths = cumsum([total_length(p) for p in paths])
    ds = circshift(lengths, 1)
    ds[1] = 0
    dx,dy = zeros(num_paths), zeros(num_paths)
    for i = 2:num_paths
        dx[i],dy[i] = position_change(paths[i-1])
    end
    dx = cumsum(dx)
    dy = cumsum(dy)
    κ = [curvature(path,0.0) for path in paths]
    DubinsPath(paths, lengths, ds, dx, dy, κ)
end

total_length(path::DubinsPath) = path.lengths[end]
final_angle(path::DubinsPath) = final_angle(path.paths[end])
function curvature(path::DubinsPath, s)
    i = match_segment(path, s)::Int
    return path.curvature[i]
end

function match_segment(path::DubinsPath, s::AbstractVector)
    [searchsortedfirst(path.ds, s_) for s_ in s]
end

function match_segment(path::DubinsPath, s::Real)
    j = searchsortedlast(path.ds, s)
    min(j, length(path.ds))
end

function localToGlobal(path::DubinsPath, s, e)
    x,y = zero(s), zero(s)
    for i in eachindex(s)
        j = match_segment(path, s[i])
        s_ = s[i] - path.ds[j]
        x[i],y[i] = localToGlobal(path.paths[j], s_, e[i])
        x[i] += path.dx[j]
        y[i] += path.dy[j]
    end
    return x,y
end

function Plots.plot(path::AbstractPath; kwargs...)
    s = range(0, total_length(path), length=101)
    e = zero(s)
    x,y = localToGlobal(path, s, e)
    plot(x,y; kwargs...)
end

function Plots.plot!(path::AbstractPath; kwargs...)
    s = range(0, total_length(path), length=101)
    e = zero(s)
    x,y = localToGlobal(path, s, e)
    plot!(x,y; kwargs...)
end

"""
Get a trajectory of s, ϕ, κ from a set of E,N pairs defining the path, where
    s is the distance along the path
    ϕ is the heading wrt the X-axis (E)
    κ is the curvature of the road, in 1/m
"""
function pathToLocal(X,Y)
    # Parameterize the curve from 0 to 1
    N = length(X)
    t = range(0,1,length=N)

    # Create cubic interpolations
    itpX = CubicSplineInterpolation(t,X)
    itpY = CubicSplineInterpolation(t,Y)

    # Compute the gradient vectors
    g = map(t) do i
        dx = Interpolations.gradient(itpX, i)[1]
        dy = Interpolations.gradient(itpY, i)[1]
        @SVector [dx,dy]
    end

    # Compute the hessian vectors
    h = map(t) do i
        ddx = Interpolations.hessian(itpX, i)[1]
        ddy = Interpolations.hessian(itpY, i)[1]
        @SVector [ddx,ddy]
    end

    # Calculate the heading
    ϕ = map(x->atan(x[2],x[1]), g)

    # Calculate the curvature
    κ = map(1:N) do k
        r′ = g[k]
        r″ = h[k]
        k = norm(r′ × r″) / norm(r′)^3
    end

    # Calulcate the path length: ∫ₐᵇ ‖r′(t)‖ dt
    x,w = gauss(N,0,1)
    @show any(isnan,x)
    r = map(x) do i
        dx = Interpolations.gradient(itpX, i)[1]
        dy = Interpolations.gradient(itpY, i)[1]
        sqrt(dx^2 + dy^2)
    end
    s = cumsum(r .* w)

    # Make the first path length zero
    s[1] = 0.0

    return s, ϕ, κ
end

"""
Attempt to recover the E,N coordanates from a set of local coordinates.
Doesn't work super well due to integration drift
"""
function localToGlobal(s,ϕ,κ)
    N = length(s)
    @assert length(ϕ) == N
    @assert length(κ) == N

    X = zeros(N)
    Y = zeros(N)
    for k = 2:N
        r = 1/mean(κ[k-1:k])
        ds = s[k] - s[k-1]
        θ = ds/r
        psi = mean(ϕ[k-1:k])
        psi = ϕ[k-1]
        xc = X[k-1] - r*sin(psi)
        yc = Y[k-1] + r*cos(psi)
        X[k] = xc + r*sin(psi + θ)
        Y[k] = yc - r*cos(psi + θ)
    end
    return X,Y
end

"""
Get the E,N coordinate given the distance along the path (s) and the longitudinal error (e)
    given a nominal path described by X,Y,S,Φ
Inputs s,e can be vectors or scalars
"""
function localToGlobal(path::CarPath, s, e)
    plot(path.X,path.Y, aspect_ratio=:equal)
    X,Y,S,Φ = path.X, path.Y, path.s, path.ϕ
    itpX = interpolate((S,),X, Gridded(Linear()))
    itpY = interpolate((S,),Y, Gridded(Linear()))
    itpΦ = interpolate((S,),Φ, Gridded(Linear()))
    nomX = itpX.(s)
    nomY = itpY.(s)
    nomΦ = itpΦ.(s)
    x = @. nomX - e*sin(nomΦ)
    y = @. nomY + e*cos(nomΦ)
    return x,y
end

function localToGlobal(path::CirclePath, s, e)
    r = path.r
    t = range(0,2pi,length=101)
    X,Y = r*cos.(t), r*sin.(t)
    plot(X,Y, aspect_ratio=:equal)
    θ = s ./ r
    radius = r .+ e
    x = radius .* cos.(θ)
    y = radius .* sin.(θ)
    return x,y
end


text:
# Independent Component Analysis

#### FastICA type

abstract type AbstractICAAlg end

struct FastICA <: AbstractICAAlg end

struct ICA{T<:Real}
    mean::Vector{T}   # mean vector, of length m (or empty to indicate zero mean)
    W::Matrix{T}      # component coefficient matrix, of size (m, k)
end

indim(M::ICA) = size(M.W, 1)
outdim(M::ICA) = size(M.W, 2)
mean(M::ICA) = fullmean(indim(M), M.mean)

transform(M::ICA, x::AbstractVecOrMat) = transpose(M.W) * centralize(x, M.mean)

#### core algorithm

# the abstract type for all g functions:
#
# Let f be an instance of such type, then
#
#   evaluate(f, x) --> (v, d)
#
# It returns a function value v, and derivative d
#
abstract type ICAGDeriv end

struct Tanh{T} <: ICAGDeriv
    a::T
end

evaluate(f::Tanh{T}, x::T) where T<:Real = (a = f.a; t = tanh(a * x); (t, a * (1 - t * t)))

function update_UE!(f::Tanh{T}, U::AbstractMatrix{T}, E1::AbstractVector{T}) where T
    n,k = size(U)
    _s = zero(T)
    a = f.a
    @inbounds for j = 1:k
        @inbounds @fastmath for i = 1:n
            t = tanh(a * U[i,j])
            U[i,j] = t
            _s += a * (1 - t^2)
        end
        E1[j] = _s / n
    end
end

struct Gaus <: ICAGDeriv end

evaluate(f::Gaus, x::T) where T<:Real = (x2 = x * x; e = exp(-x2/2); (x * e, (1 - x2) * e))

function update_UE!(f::Gaus, U::AbstractMatrix{T}, E1::AbstractVector{T}) where T
    n,k = size(U)
    _s = zero(T)
    @inbounds for j = 1:k
        for i = 1:n
            u = U[i,j]
            u2 = u^2
            e = exp(-u2/2)
            U[i,j] = u * e
            _s += (1 - u2) * e
        end
        E1[j] = _s / n
    end
end

# Fast ICA
#
# Reference:
#
#   Aapo Hyvarinen and Erkki Oja
#   Independent Component Analysis: Algorithms and Applications.
#   Neural Network 13(4-5), 2000.
#
function ica!(::FastICA, W::DenseMatrix{T},      # initialized component matrix, size (m, k)
                  X::DenseMatrix{T},      # (whitened) observation sample matrix, size(m, n)
                  fun::ICAGDeriv,         # approximate neg-entropy functor
                  maxiter::Int,           # maximum number of iterations
                  tol::Real) where T<:Real# convergence tolerance

    # argument checking
    m = size(W, 1)
    k = size(W, 2)
    size(X, 1) == m || throw(DimensionMismatch("Sizes of W and X mismatch."))
    n = size(X, 2)
    k <= min(m, n) || throw(DimensionMismatch("k must not exceed min(m, n)."))

    @debug "FastICA Algorithm" m=m n=n k=k

    # pre-allocated storage
    Wp = Matrix{T}(undef, m, k)    # to store previous version of W
    U  = Matrix{T}(undef, n, k)    # to store w'x & g(w'x)
    Y  = Matrix{T}(undef, m, k)    # to store E{x g(w'x)} for components
    E1 = Vector{T}(undef, k)       # store E{g'(w'x)} for components

    # normalize each column
    for j = 1:k
        w = view(W,:,j)
        rmul!(w, one(T) / sqrt(sum(abs2, w)))
    end

    # main loop
    chg = T(NaN)
    t = 0
    converged = false
    @inbounds while !converged && t < maxiter
        t += 1
        copyto!(Wp, W)

        # apply W of previous step
        mul!(U, transpose(X), W) # u <- w'x

        # compute g(w'x) --> U and E{g'(w'x)} --> E1
        update_UE!(fun, U, E1)

        # compute E{x g(w'x)} --> Y
        rmul!(mul!(Y, X, U), one(T) / n)

        # update w: E{x g(w'x)} - E{g'(w'x)} w := y - e1 * w
        for j = 1:k
            w = view(W,:,j)
            y = view(Y,:,j)
            e1 = E1[j]
            @. w = y - e1 * w
        end

        # symmetric decorrelation: W <- W * (W'W)^{-1/2}
        copyto!(W, W * _invsqrtm!(W'W))

        # compare with Wp to evaluate a conversion change
        chg = maximum(abs.(abs.(diag(W*Wp')) .- one(T)))
        converged = (chg < tol)

        @debug "Iteration $t" change=chg tolerance=tol
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))
    return W
end

#### interface function

function fit(::Type{ICA}, X::AbstractMatrix{T},             # sample matrix, size (m, n)
                          k::Int;                           # number of independent components
                          alg::AbstractICAAlg,             # choice of algorithm
                          fun::ICAGDeriv=Tanh(one(T)),      # approx neg-entropy functor
                          do_whiten::Bool=true,             # whether to perform pre-whitening
                          maxiter::Integer=100,             # maximum number of iterations
                          tol::Real=1.0e-6,                 # convergence tolerance
                          mean=nothing,                     # pre-computed mean
                          winit::Matrix{T}=zeros(T,0,0)) where T<:Real  # init guess of W, size (m, k)

    # check input arguments
    m, n = size(X)
    n > 1 || error("There must be more than one samples, i.e. n > 1.")
    k <= min(m, n) || error("k must not exceed min(m, n).")
    maxiter > 1 || error("maxiter must be greater than 1.")
    tol > 0 || error("tol must be positive.")

    # preprocess data
    mv = preprocess_mean(X, mean)
    Z::Matrix{T} = centralize(X, mv)

    W0= zeros(T, 0, 0)  # whitening matrix
    if do_whiten
        C = rmul!(Z * transpose(Z), 1.0 / (n - 1))
        Efac = eigen(C)
        ord = sortperm(Efac.values; rev=true)
        (v, P) = extract_kv(Efac, ord, k)
        W0 = rmul!(P, Diagonal(1 ./ sqrt.(v)))
        Z = W0'Z
    end

    # initialize
    W = (isempty(winit) ? randn(T, size(Z,1), k) : copy(winit))

    # invoke core algorithm
    ica!(alg, W, Z, fun, maxiter, tol)

    # construct model
    if do_whiten
        W = W0 * W
    end
    return ICA(mv, W)
end

# function fastica!(W::DenseMatrix{T},      # initialized component matrix, size (m, k)
#                   X::DenseMatrix{T},      # (whitened) observation sample matrix, size(m, n)
#                   fun::ICAGDeriv{T},      # approximate neg-entropy functor
#                   maxiter::Int,           # maximum number of iterations
#                   tol::Real) where T<:Real# convergence tolerance
#
#     # argument checking
#     m = size(W, 1)
#     k = size(W, 2)
#     size(X, 1) == m || throw(DimensionMismatch("Sizes of W and X mismatch."))
#     n = size(X, 2)
#     k <= min(m, n) || throw(DimensionMismatch("k must not exceed min(m, n)."))
#
#     @debug "FastICA Algorithm" m=m n=n k=k
#
#     # pre-allocated storage
#     Wp = Matrix{T}(undef, m, k)    # to store previous version of W
#     U  = Matrix{T}(undef, n, k)    # to store w'x & g(w'x)
#     Y  = Matrix{T}(undef, m, k)    # to store E{x g(w'x)} for components
#     E1 = Vector{T}(undef, k)       # store E{g'(w'x)} for components
#
#     # normalize each column
#     for j = 1:k
#         w = view(W,:,j)
#         rmul!(w, one(T) / sqrt(sum(abs2, w)))
#     end
#
#     # main loop
#     chg = NaN
#     t = 0
#     converged = false
#     @inbounds while !converged && t < maxiter
#         t += 1
#         copyto!(Wp, W)
#
#         # apply W of previous step
#         mul!(U, transpose(X), W) # u <- w'x
#
#         # compute g(w'x) --> U and E{g'(w'x)} --> E1
#         _s = zero(T)
#         for j = 1:k
#             for i = 1:n
#                 u, v = evaluate(fun, U[i,j])
#                 U[i,j] = u
#                 _s += v
#             end
#             E1[j] = _s / n
#         end
#
#         # compute E{x g(w'x)} --> Y
#         rmul!(mul!(Y, X, U), one(T) / n)
#
#         # update w: E{x g(w'x)} - E{g'(w'x)} w := y - e1 * w
#         for j = 1:k
#             w = view(W,:,j)
#             y = view(Y,:,j)
#             e1 = E1[j]
#             @. w = y - e1 * w
#
#         end
#
#         # symmetric decorrelation: W <- W * (W'W)^{-1/2}
#         copyto!(W, W * _invsqrtm!(W'W))
#
#         # compare with Wp to evaluate a conversion change
#         chg = maximum(abs.(abs.(diag(W*Wp')) .- one(T)))
#         converged = (chg < tol)
#
#         @debug "Iteration $t" change=chg tolerance=tol
#     end
#     converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))
#     return W
# end


text:
import PipeLayout: nonzero, isapproxlt, isapproxgt, isapproxin, isSOS1,
       isSOS2, select_subset

@testset "numerical comparisons" begin
    @test nonzero(0.0) == false
    @test nonzero(1e-15) == false
    @test nonzero(-1e-15) == false
    @test nonzero(1e-3) == true
    @test nonzero(-1e-3) == true

    @test isapproxlt(1.5, 2.0)
    @test isapproxlt(2.0 + 1e-10, 2.0)
    @test !isapproxlt(2.1, 2.0)

    @test isapproxgt(2.5, 2.0)
    @test isapproxgt(2.0 - 1e-10, 2.0)
    @test !isapproxgt(1.9, 2.0)

    @test !isapproxin(0.5, 1.0, 2.0)
    @test isapproxin(1.0 - 1e-9, 1.0, 2.0)
    @test isapproxin(1.5, 1.0, 2.0)
    @test isapproxin(2.0 + 1e-9, 1.0, 2.0)
    @test !isapproxin(2.5, 1.0, 2.0)

    @test !isSOS1([-1., 0.0])
    @test isSOS1( [0.0, 0.0])
    @test isSOS1( [0.5, 0.0])
    @test isSOS1( [0.0, 0.5])
    @test !isSOS1([0.3, 0.5])

    @test !isSOS2([-1., 0.0, 0.0])
    @test isSOS2( [0.0, 0.0, 0.0])
    @test isSOS2( [0.5, 0.0, 0.0])
    @test isSOS2( [0.0, 0.5, 0.3])
    @test !isSOS2([0.3, 0.0, 0.5])
    @test !isSOS2([1.0, 0.5, 0.3])
end

@testset "random subset sampling" begin
    base = collect(1:6)
    sample0 = select_subset(base, 0)
    sample3 = select_subset(base, 3)
    sample6 = select_subset(base, 6)

    @test length(sample0) == 0
    @test length(sample3) == 3
    @test unique(sample3) == sample3
    @test length(sample6) == 6
    @test unique(sample6) == sample6
    @test_throws AssertionError select_subset(base, 9)
end


text:
id = 0;
counter = 0;
maxdiffreduvscriti = 0.0;
maxfluid = "";
uneq = Set();
critphasefail = Set();
for fluid in coolpropfluids
  id+=1;
  pure = get_fluid_param_string(fluid, "pure");
  tcrit = PropsSI("TCRIT", fluid);
  pcrit = PropsSI("PCRIT", fluid);
  try
    (PhaseSI("P", pcrit, "T", tcrit, fluid) != "critical_point") && push!(critphasefail, fluid);
  catch err
    push!(critphasefail, fluid);
  end
  @test PhaseSI("P", pcrit+50000, "T", tcrit+3, fluid)=="supercritical";
  @test PhaseSI("P", pcrit+50000, "T", tcrit-3, fluid)=="supercritical_liquid";
  @test PhaseSI("P", pcrit-50000, "T", tcrit+3, fluid)=="supercritical_gas";
  psat = PropsSI("P","T",tcrit-3,"Q",1, fluid);
  t3 = PropsSI("Ttriple", fluid);
  p3 = PropsSI("ptriple", fluid);
  @test PhaseSI("P", p3, "T", (t3+tcrit)/2, fluid)=="gas";
  @test PhaseSI("P", (pcrit+p3)/2, "T", (t3+tcrit)/2, fluid)=="liquid";
  (pure == "true") && (tcrit != PropsSI("T_REDUCING", fluid)) && push!(uneq, fluid);
  diffreduvscriti = abs(PropsSI("TCRIT", fluid) - PropsSI("T_REDUCING", fluid));
  if (diffreduvscriti > maxdiffreduvscriti)
    maxfluid = fluid;
    maxdiffreduvscriti = diffreduvscriti;
  end
end
println("max diff between reducing vs critical point T: $maxdiffreduvscriti for $maxfluid");
@test uneq == Set(fails_tcrit_eq_treducing);
println("different reducing vs critical point: $uneq");
@test issubset(critphasefail, Set(fails_critical_point));
println("fails to get phase for critical point: $critphasefail");


text:
export ingredients, merge!, similarname, alcohols, sides

using Lazy: @forward
@forward DrinkBook.recipes (Base.length, Base.iterate, Base.getindex)
@forward Systembolag.alcohols (Base.length, Base.iterate, Base.getindex)
@forward Supermarket.sides (Base.length, Base.iterate, Base.getindex)

ingredients(r::Recipe) = [r.alcohols; r.sides]
alcohols(v)            = [i for i in v if i isa Alcohol]
sides(v)               = [i for i in v if i isa Side]
Base.length(r::Recipe) = length(r.alcohols) + length(r.sides)

Base.:(==)(x::Systembolag, y::Systembolag) = x.alcohols == y.alcohols
Base.:(==)(x::Supermarket, y::Supermarket) = x.sides == y.sides
Base.:(==)(x::BarCabinet, y::BarCabinet) = x.sides == y.sides && x.alcohols == y.alcohols
Base.:(==)(x::DrinkBook, y::DrinkBook) = x.recipes == y.recipes

Base.:(==)(x::Ingredient, y::Ingredient) = x.name == y.name
function Base.:(==)(x::Recipe, y::Recipe)
    x.name == y.name &&
    x.ingredients == y.ingredients &&
    x.description == y.description
end

name(dc::DrinkComponent) = name(dc.ingredient)
name(x::ElementType) = x.name
name(s::String) = s

function Base.iterate(b::BarCabinet, state=1)
    state == length(b) && return nothing
    state <= length(b.alcohols) && return b.alcohols[state], state+1
    return b.sides[state - length(b.alcohols)], state+1
end
Base.length(b::BarCabinet) = length(b.alcohols) + length(b.sides)

function Base.push!(db::DrinkBook, x::Recipe)
    fuzzyin(name(x), db) && return false
    push!(db.recepies, x)
    return true
end

function Base.push!(db::Systembolag, x::Alcohol)
    fuzzyin(x, db) && return false
    push!(db.alcohols, x)
    return true
end

function Base.push!(db::Supermarket, x::Side)
    fuzzyin(x, db) && return false
    push!(db.sides, x)
    return true
end

function Base.push!(db::BarCabinet, x::Ingredient)
    fuzzyin(x, db) && return false
    x isa Alcohol && push!(db.alcohols, x)
    x isa Side && push!(db.sides, x)
    return true
end

function similarname(c1::T) where T <: ElementType
    function (c2::T)
        s = compare(Levenshtein(), lowercase(name(c1)), lowercase(name(c2)))
        s > SIMILARITY_THRESHOLD
    end
end
similarname(c1::T,c2::T) where T <: ElementType = similarname(c1)(c2)

function fuzzyin(item, collection)
    sim = similarname(item)
    any(sim(c) for c in collection)
end

"""
`push!`es all elements in `x2` into `x1`
"""
function Base.merge!(x1::T,x2::T) where T <: ContainerType
    foreach(e->push!(x1,e), x2)
end

"""
avaiable::DrinkBook, missing::DB = available_recipes(drinkbook::DrinkBook, cabinet::BarCabinet; ignore_missing_sides=false)

Return two `DrinkBook`s with the drink recipes that can be mixed with the ingredients available in `cabinet`, and those that can't be mixed.
"""
function available_recipes(drinkbook::DrinkBook, cabinet::BarCabinet; ignore_missing_sides=false)
    available = Recipe[]
    missing = Recipe[]
    for recipe in drinkbook
        mixable = all(!component.mandatory ||
            ignore_missing_sides && component isa Side ||
            component.ingredient ∈ cabinet
            for component in recipe)
        mixable ? push!(available, recipe) : push!(missing, recipe)
    end
    DrinkBook(available), DrinkBook(missing)
end

function which_ingredient_to_buy(drinkbook::DrinkBook, cabinet::BarCabinet, depth=1; ignore_missing_sides=false)
    error("Not yet implemented")
    # TODO: perform tree seach of depth `depth` to find ingredients not in cabinet that would make the largest number of new recipes mixable
    # Pseudo-code:
    _, candidates = available_recipes(drinkbook, cabinet; ignore_missing_sides=ignore_missing_sides)
    candidates = filter(recipe->countmissing(recipe, cabinet) <= depth, candidates)
    for ingredient in missing_ingredients(cabinet, systembolag)
        for recipe in drinkbook
        end
    end
end

countmissing(recipe, cabinet) = count(i->i ∉ cabinet, recipe.ingredients)


text:
using TableIOInterface
using Test

@testset "TableIOInterface.jl" begin

    fnames = ["test.csv", "test.jdf", "test.parquet", "test.arrow", "test.xlsx", "test.json",
                "test.zip", "test.sqlite"]

    @testset "File type ID" begin

        expected_file_types = [TableIOInterface.CSVFormat, TableIOInterface.JDFFormat, TableIOInterface.ParquetFormat,
        TableIOInterface.ArrowFormat, TableIOInterface.ExcelFormat, TableIOInterface.JSONFormat,
        TableIOInterface.ZippedFormat, TableIOInterface.SQLiteFormat]

        for (i, fname) ∈ enumerate(fnames)
            file_type = get_file_type(fname)
            @test file_type == expected_file_types[i]()
        end
    end

    @testset "Sample code generation" begin
        expected_sample_code = [
            """df_test_csv = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.csv")); copycols=false)
            end""",

            """df_test_jdf = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.jdf")); copycols=false)
            end""",

            """df_test_parquet = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.parquet")); copycols=false)
            end""",

            """df_test_arrow = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.arrow")); copycols=false)
            end""",

            """df_test_xlsx = let
                import TableIO
                import DataFrames
                df_test_xlsx = DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.xlsx"),
                    # "Sheet1", # uncomment this to specify the sheet name to be imported
                    ); copycols=false)
            end""",

            """df_test_json = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.json")); copycols=true)
            end""",

            """df_test_zip = let
                import TableIO
                import DataFrames
                df_test_zip = DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.zip"),
                    # "filename_in_archive", # uncomment this to specify a specific file name in the archive
                    ); copycols=false)
            end""",

            """df_test_sqlite = let
                import TableIO
                import DataFrames
                df_test_sqlite = DataFrames.DataFrame(TableIO.read_table(joinpath(dirname(@__FILE__), "my_notebook.assets", "test.sqlite"),
                    "tablename", # define the table name here
                    ); copycols=false)
            end""",

        ]

        for (i, fname) ∈ enumerate(fnames)
            sample_code = get_example_code("my_notebook.assets", fname)
            @test sample_code == expected_sample_code[i]
            #println(sample_code)
        end

    end

end


text:
module Jlo

function test()
    print("test")
end

function gaussian(x_px, mean_px, stdev_px)
    exp(-((x_px-mean_px)/(2*stdev_px)).^2)
end

#image has to be square
function propTF(E0, L_m, lambda_m, z_m)
    fx=Jlo.fs(linspace(0,L_m,size(E0,1)))
    Fx, Fy =  Jlo.meshgrid(fx,fx)
    H = exp((-im)*pi*lambda_m*z_m*(Fx.^2+Fy.^2))
    E_prop = Jlo.ft(H.*Jlo.ift(E0))
    return E_prop
end

function propIR
end

function propHF
end

function crop
end

function pad
end

function fs(t)
    ff = 1/(t[end]-t[1])
    Nf = length(t)
    df =   ff./Nf
    f = (0:df:(ff-df)) - (ff-mod(Nf,2)*df)/2

    return f
end

function ft(x)
    y = fftshift(ifft(ifftshift(x)))
    return y
end

"centered fourier transform"
function ift(x)
    y = fftshift(fft(ifftshift(x)))
    return y
end

function meshgrid(vx, vy)
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

end

max_stars_repo_name: awojdyla/jlo

max_issues_repo_name: awojdyla/jlo

max_forks_repo_name: awojdyla/jlo


text:
## circulant_embedding.jl : Gaussian random field generator using fft; only for uniformly spaced GRFs

"""
CirculantEmbedding <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) that uses FFT to compute samples of the Gaussian random field. Circulant embedding can only be applied if the points are specified on a structured grid.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,CirculantEmbedding(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a circulant embedding

julia> contourf(grf)
[...]

julia> plot_eigenvalues(grf)
[...]

```
With appropriate ordering, the covariance matrix of a Gaussian random field is a (nested block) Toeplitz matrix. This matrix can be embedded into a larger (nested block) circulant matrix, whose eigenvalues can be rapidly computed using FFT. A difficulty here is that although the covariance matrix is positive semi-definite, its circulant extension in general is not. As a remedy, one can add so-called *ghost points* outside the domain of interest using the optional flag `padding`.
```jldoctest
julia> m = Matern(1,1)
Matérn (λ=1.0, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(1,m)
1d Matérn covariance function (λ=1.0, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0,stop = 1,length = 256)
0.0:0.00392156862745098:1.0

julia> grf = GaussianRandomField(c,CirculantEmbedding(),pts)
WARNING: negative eigenvalue -0.001465931113950698 detected, Gaussian random field will be approximated (ignoring all negative eigenvalues)
WARNING: increase padding if possible
Gaussian random field with 1d Matérn covariance function (λ=1.0, ν=1.0, σ=1.0, p=2.0) on a 256-point structured grid, using a circulant embedding

julia> grf = GaussianRandomField(c,CirculantEmbedding(),pts,padding=5)
Gaussian random field with 1d Matérn covariance function (λ=1.0, ν=1.0, σ=1.0, p=2.0) on a 256-point structured grid, using a circulant embedding

```
For faster sampling (but slower initialization), use the optional argument `measure` (default=true).

See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`KarhunenLoeve`](@ref)
"""
struct CirculantEmbedding <: GaussianRandomFieldGenerator end

const CirculantGRF = GaussianRandomField{C,CirculantEmbedding} where {C}

function _GaussianRandomField(mean,cov::CovarianceFunction{d},method::CirculantEmbedding,pts...;padding=1,measure=true) where {d}

    # add ghost points by padding
    padded_pts = pad.(pts,padding)
    n = length.(padded_pts)
    n2 = n[2:d]
    padded_pts = normalize.(padded_pts)
    padded_pts = (padded_pts[1],mirror.(padded_pts[2:d])...) # do not mirror dimension 1

    # compute eigenvalues of circulant matrix
    c = zeros(n[1],(2 .*n2.-1)...)
    for (i,idx) in enumerate(Base.product(padded_pts...))
        c[i] = apply(cov.cov,collect(idx))
    end
    c̃ = zeros(n[1],(2 .* n2)...)
    c̃[:,broadcast(:,2,2 .*n2)...] = c
    c̃ = fftshift(c̃,2:d)
    Λ = irfft(c̃,2*size(c̃,1)-1)
    Λ⁺ = zeros(size(Λ))
    Λ⁻ = zeros(size(Λ))
    for i in 1:length(Λ)
        if Λ[i] < 0
            Λ⁻[i] = -Λ[i]
        else
            Λ⁺[i] = Λ[i]
        end
    end
    Λ_max = maximum(Λ⁻)
    Λ_max > 0. && ( warn("negative eigenvalue $(-Λ_max) detected, Gaussian random field will be approximated (ignoring all negative eigenvalues)"); warn("increase padding if possible") )

    # optimize
    Σ = sqrt.(Λ⁺)
    P = measure ? plan_fft(Σ) : plan_fft(Σ,flags=FFTW.MEASURE)

    GaussianRandomField{typeof(cov),CirculantEmbedding,typeof(pts)}(mean,cov,pts,(Σ,P))
end

normalize(pts::R where {R<:AbstractRange}) = pts.-pts[1]
mirror(pts::R where {R<:AbstractRange}) = -pts[end]:(pts[2]-pts[1]):pts[end]

function pad(x,n)
    x0 = x[1]
    xn = x[end]
    dx = x[2]-x[1]
    x0:dx:n*xn
end

# returns the required dimension of the random points
randdim(grf::CirculantGRF) = length(grf.data[1])

# sample function
function _sample(grf::CirculantGRF, xi)
    v = grf.data[1]
    y = v.*reshape(xi,size(v))
    P = grf.data[2]
    w = P*y # fft
    w = real(w) + imag(w)
    n = length.(grf.pts)
    z = w[broadcast(:,1,n)...] # select appropriate elements

    grf.mean + std(grf.cov)*z
end

show(io::IO,::CirculantEmbedding) = print(io,"circulant embedding")


text:
function init_observable_mappings!(m::PoolModel)
    observables = OrderedDict{Symbol,Observable}()

    identity_transform_805 = function(levels)
        return levels[!,:p805]
    end
    identity_transform_904 = function(levels)
        return levels[!,:p904]
    end

    if subspec(m) in ["ss1", "ss3", "ss5", "ss7", "ss9"]
        ############################################################################
        ## 1. Model 904 predictive density
        ############################################################################

        observables[:Model904] = Observable(:Model904, [:p904__wrongorigmatlab],
                                            identity_transform_904,
                                            identity_transform_904,
                                            "Model 904 predictive density",
                                            "Model 904 conditional predictive density scores")

        ############################################################################
        ## 2. Model 805 predictive density
        ############################################################################
        observables[:Model805] = Observable(:Model805, [:p805__wrongorigmatlab],
                                            identity_transform_805,
                                            identity_transform_805,
                                            "Model 805 predictive density",
                                            "Model 805 conditional predictive density scores")
    elseif subspec(m) in ["ss2", "ss4", "ss6", "ss8", "ss10"]
        ############################################################################
        ## 1. Model 904 predictive density
        ############################################################################
        observables[:Model904] = Observable(:Model904, [:p904__right805orig904],
                                            identity_transform_904,
                                            identity_transform_904,
                                            "Model 904 predictive density",
                                            "Model 904 conditional predictive density scores")

        ############################################################################
        ## 2. Model 805 predictive density
        ############################################################################
        observables[:Model805] = Observable(:Model805, [:p805__right805orig904],
                                            identity_transform_805,
                                            identity_transform_805,
                                            "Model 805 predictive density",
                                            "Model 805 conditional predictive density scores")

    end

    m.observable_mappings = observables
end


text:
@testset "CoreModel generic interface" begin
    model = load_model(CoreModel, model_paths["e_coli_core.mat"])

    @test reaction_stoichiometry(model, "EX_ac_e") == Dict("ac_e" => -1)
    @test reaction_stoichiometry(model, 44) == Dict("ac_e" => -1)
end

@testset "Conversion from and to StandardModel" begin
    cm = load_model(CoreModel, model_paths["e_coli_core.mat"])

    sm = convert(StandardModel, cm)
    cm2 = convert(CoreModel, sm)

    @test Set(reactions(cm)) == Set(reactions(sm))
    @test Set(reactions(cm)) == Set(reactions(cm2))

    @test reaction_gene_association(sm, reactions(sm)[1]) ==
          reaction_gene_association(cm, reactions(sm)[1])
end