text:
import sbp
using MAT
using LinearAlgebra
include("../common/wave_equation.jl")
include("../common/bound.jl")
include("../common/discretization.jl")
disc = Discretization
pde = WaveEquation

datafile = ARGS[1]
nx = 100
nsamples = 100
order = 4
alpha = 10
dR = range(-1.0, stop=1.0, length=nsamples)

vars = pde.sizes(nx)
boundary_points = [1, nx]
boundary_conditions = [pde.non_reflective, () -> pde.leakage(alpha)]
sides = [pde.LEFT, pde.RIGHT]
boundaries = [disc.boundary(vars, [1]), disc.boundary(vars, [nx])]
nl = disc.normal(nx)

ops = disc.OperatorsType(nx, order)
h = ops.h
M0 = disc.discretization(pde.A(), ops)

sigma_d = LinearAlgebra.svd(Matrix(ops.Dx) * h)
rd = sigma_d.S[2]
sd = sigma_d.S[1]
rl = ops.Hxi[1,1] * h
sl = rl
println("estimated constants: $rd <= h\\|d\\|_2 <= $sd")
println("estimated constants: $rl <= h\\|l\\|_2 <= $sl")
i = 1

Mnorm = zeros(nsamples)
sigma_norm_left = zeros(nsamples)
sigma_norm_right = zeros(nsamples)
sigma_norm_upper = zeros(nsamples)
sigma_norm_lower = zeros(nsamples)
for dRi  in dR

        coef = [0, dRi]

        M = M0

        """
         Go through each boundary and construct projection matrices, and sat
         terms
        """
        for (dRi, bnd, bc, bp, side) in zip(coef, boundaries, boundary_conditions,
                                      boundary_points, sides)
                r = bc()
                dr = pde.reflection_matrix(dRi)
                Ai = pde.A(nl[bp])
                P = disc.projection_matrix(Ai, r, dr)
                sat_term, sat_term_data = pde.sat(Ai, side, ops, P)
                M = M + sat_term
                if side == pde.RIGHT
                 sigma_norm_right[i] = LinearAlgebra.opnorm(Ai * P)
                end
                if side == pde.LEFT
                 sigma_norm_left[i] = LinearAlgebra.opnorm(Ai * P)
                end
        end

        sigma_norm_upper[i] = sigma_norm_left[i] + sigma_norm_right[i]
        sigma_norm_lower[i] = sigma_norm_right[i] - sigma_norm_left[i]

        # Spectral radius
        eig_max = maximum(disc.check_stability(M))
        if (eig_max > 0)
                println("Warning: Unstable. eig_max = $eig_max")
        end
        Mnorm[i] = h * LinearAlgebra.opnorm(M)
        println("\\|M\\|_2 = ", Mnorm[i], " ", sigma_norm_lower[i], " <= \\|Sigma\\|_2 <= ", sigma_norm_upper[i])
        global i += 1
end
matwrite(datafile, Dict("dR" => collect(dR),
                                   "Mnorm" => Mnorm,
                                   "rd" => rd,
                                   "sd" => sd,
                                   "rl" => rl,
                                   "sl" => sl,
                                   "sigma_norm_upper" => sigma_norm_upper,
                                   "sigma_norm_lower" => sigma_norm_lower)
         )

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


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

module Visualize

export trajPlot
export pgfplot_hor, pgfplot_alt, pgfplot_heading, pgfplot_vrate

import Compat.ASCIIString

using ..DefineSave
using ..SaveHelpers
using ..AddSupplementary

include("../helpers/TikzUtils.jl")
include("vis_captions.jl")
using .TikzUtils
using .VisCaptions

using Iterators
using TikzPictures
import PGFPlots
import PGFPlots: Plots, Axis, GroupPlot

using RLESUtils, MathUtils, PGFPlotUtils

const RA_STYLE_MAP = [
    ((ra, h_d) -> ra && abs(h_d) < 5, "mark options={color=gray}, mark=*"),
    ((ra, h_d) -> ra && 5 <= h_d < 30, "mark options={color=orange}, mark=*"),
    ((ra, h_d) -> ra && 30 <= h_d, "mark options={color=red}, mark=*"),
    ((ra, h_d) -> ra && -30 < h_d <= -5, "mark options={color=cyan}, mark=*"),
    ((ra, h_d) -> ra && h_d <= -30, "mark options={color=violet}, mark=*")
    ]

const RESPONSE_STYLE_MAP = [
    (r -> r == "stay", "mark options={color=black}, mark=-"),
    (r -> r == "follow", "mark options={color=black}, mark=asterisk")
    ]

#xy
function pgfplot_hor(d::TrajLog)
    plotArray = vcat(pplot_line(d, "WorldModel", :y, :x),
        pplot_startpoint(d, "WorldModel", :y, :x, "top"), # label start point
        pplot_aircraft_num(d, "WorldModel", :y, :x, startdist = 3.4)) # label aircraft numbers
    xmin, xmax, ymin, ymax = manual_axis_equal(d, "WorldModel", :y, :x)
    ax = Axis(plotArray,
        ylabel = "N ($(get_unit(d, "WorldModel", 1, :x)))",
        xlabel = "E ($(get_unit(d, "WorldModel", 1, :y)))",
        title = "Horizontal Position",
        style = "xmin=$xmin,xmax=$xmax,ymin=$ymin,ymax=$ymax,enlarge x limits=true," *
            "enlarge y limits=true,axis equal,clip mode=individual")
    ax
end

#altitude vs. time
function pgfplot_alt(d::TrajLog)
    plotArray = vcat(pplot_z_label270s(d), # label270 short
        pplot_line(d, "WorldModel", :t, :z),
        pplot_startpoint(d, "WorldModel", :t, :z, "side", overrideangle = 0), # label start point
        pplot_aircraft_num(d, "WorldModel", :t, :z, startdist = 3.4)) # label aircraft numbers

    ax = Axis(plotArray,
            xlabel = "time ($(get_unit(d, "WorldModel", 1, :t)))",
            ylabel = "h ($(get_unit(d, "WorldModel", 1, :z)))",
            title = "Altitude vs. Time",
            style = "clip=false,clip mode=individual")
    ax
end

#heading vs time
function pgfplot_heading(d::TrajLog)
    plotArray = vcat(pplot_aircraft_num(d, "Dynamics", :t, :psi, fy=psidot_from_psi), # label aircraft numbers
        pplot_line(d, "Dynamics", :t, :psi, fy=psidot_from_psi))
    ax = Axis(plotArray,
        xlabel = "time ($(get_unit(d, "Dynamics", 1, :t)))",
        ylabel = "psidot ($(get_unit(d, "Dynamics", 1, :psi))/s)",
        title = "Heading Rate vs. Time")
    ax
end

#vertical rate vs time
function pgfplot_vrate(d::TrajLog)
    plotArray = vcat(pplot_aircraft_num(d, "WorldModel", :t, :vz), # label aircraft numbers
        pplot_line(d, "WorldModel", :t, :vz))
    ax = Axis(plotArray,
        xlabel = "time ($(get_unit(d, "WorldModel", 1, :t)))",
        ylabel = "vh ($(get_unit(d, "WorldModel", 1, :vz)))",
        title = "Vertical Rate vs. Time")
    ax
end

function pgfplotLog(d::TrajLog)
    tps = TikzPicture[]

    #xy and tz group
    g = GroupPlot(2, 2, groupStyle = "horizontal sep = 2.2cm, vertical sep = 2.2cm")

    #xy
    push!(g, pgfplot_hor(d))

    #altitude vs time
    push!(g, pgfplot_alt(d))

    #heading rate vs time
    push!(g, pgfplot_heading(d))

    #vertical rate vs time
    push!(g, pgfplot_vrate(d))

    tp = PGFPlots.plot(g)
    use_geometry_package!(tp, landscape = true)
    use_aircraftshapes_package!(tp)

    push!(tps, tp)

    tps
end

function manual_axis_equal(d::TrajLog, field::AbstractString, xname::Symbol,
    yname::Symbol; fx::Function = identity, fy::Function = identity)

    xmin = ymin = realmax(Float64)
    xmax = ymax = -realmax(Float64)

    for i = 1:get_num_aircraft(d)
        df = get_log(d, field, i)

        #plot trajectories
        xvals = convert(Array, df[xname])
        yvals = convert(Array, df[yname])

        xvals = fx(xvals)
        yvals = fy(yvals)

        tempmin, tempmax = extrema(xvals)
        xmin = min(xmin, tempmin)
        xmax = max(xmax, tempmax)
        tempmin, tempmax = extrema(yvals)
        ymin = min(ymin, tempmin)
        ymax = max(ymax, tempmax)
    end

    xrange = xmax - xmin
    yrange = ymax - ymin
    plot_range = max(xrange, yrange)

    delta = abs(yrange - xrange) / 2
    if yrange > xrange
        xmin -= delta
        xmax += delta
    else
        ymin -= delta
        ymax += delta
    end

    (xmin, xmax, ymin, ymax)
end

#xname = field name of x variable
#yname = field name of y variable
function pplot_aircraft_num(d::TrajLog, field::AbstractString, xname::Symbol,
    yname::Symbol; ind_start::Int64 = 1,
                            displaystart::Bool = true,
                            displayend::Bool = true,
                            startdist::Float64 = 2.0,
                            enddist::Float64 = 2.0,
                            scale::Float64 = 0.55,
                            fx::Function = identity,
                            fy::Function = identity)

    plotArray = Plots.Plot[]

    for i = 1:get_num_aircraft(d)
        df = get_log(d, field, i)

        #plot trajectories
        xvals = convert(Array, df[xname])
        yvals = convert(Array, df[yname])

        # apply user-supplied transformations
        xvals = fx(xvals)
        yvals = fy(yvals)

        if displaystart
            # determine where to put the label for the first: left or right
            x1 = xvals[ind_start]
            y1 = yvals[ind_start]
            x2 = xvals[ind_start + 1]

            #mark aircraft number
            dir = x2 < x1  ? "right" : "left"
            push!(plotArray, Plots.Node("$i", x1, y1,
                                  style = "$(dir)=$(startdist)mm,scale=$(scale),rotate=0"))
        end

        if displayend
            # determine where to put the label for the last: left or right
            xend = xvals[end]
            yend = yvals[end]
            xend_1 = xvals[end - 1]

            #mark aircraft number
            dir = xend_1 < xend  ? "right" : "left"
            push!(plotArray, Plots.Node("$i", xend, yend,
                                  style = "$(dir)=$(enddist)mm,scale=$(scale),rotate=0"))
        end
    end
    plotArray
end

#xname = field name of x variable
#yname = field name of y variable
#view = "top" or "side" view of aircraft
#angle = angle of aircraft in degrees.  Pointing right is 0.  nothing = auto-determine from first and second points.
#minwidth = minimum width of aircraft icon in cm

function pplot_startpoint(d::TrajLog, field::AbstractString, xname::Symbol,
    yname::Symbol, view::AbstractString;
                          overrideangle::Union{Void, Real} = nothing,
                          ind_start::Int64 = 1,
                          minwidth::Float64 = 0.65,
                          fx::Function = identity,
                          fy::Function = identity)
    plotArray = Plots.Plot[]
    for i = 1:get_num_aircraft(d)
        df = get_log(d, field, i)

        #plot trajectories
        xvals = convert(Array, df[xname])
        yvals = convert(Array, df[yname])

        # apply user-supplied transformations
        xvals = fx(xvals)
        yvals = fy(yvals)

        x1 = xvals[ind_start]
        y1 = yvals[ind_start]

        # determine angle of aircraft if not given
        if overrideangle == nothing
            x2 = xvals[ind_start + 1]
            y2 = yvals[ind_start + 1]

            angle = atan2(y2 - y1, x2 - x1) |> rad2deg
        else
            angle = overrideangle
        end

        #mark aircraft start point
        push!(plotArray, Plots.Node("", x1, y1,
            style = "aircraft $view,draw=white,thin,fill=black,minimum width=$(minwidth)cm,rotate=$angle"))
        #\node [aircraft top,fill=black,minimum width=1cm,rotate=30] at (0,0) {};
        #\node at (axis cs:-14746.707634583, 8514.015622487) [left=2.0mm,scale=0.55,rotate=0] {2};
    end

    plotArray
end

function pplot_z_label270s(d::TrajLog; start_time::Int64 = 0, end_time::Int64 = 51,
    label_scale::Float64 = 0.45)

    @assert start_time < end_time

    plotArray = Plots.Plot[]
    for i = 1:get_num_aircraft(d)
        wm = d["WorldModel_$i"]
        cas = d["CAS_$i"]
        ts = convert(Array, wm[:t])
        filter!(t -> start_time <= t <= end_time, ts) #filter based on start/end times

        #short labels
        prev_label = ""

        dirs = closest_is_above(d, i)
        xs = convert(Array, wm[:t])
        ys = convert(Array, wm[:z])
        for t in ts
            label = cas[t, :label270_short]
            if prev_label != label
                label_ = replace(label, "_", "\\_") #convert underscores to latex escape sequence
                label_ = "\\textbf{$(label_)}"
                dir = dirs[t] ? "left" : "right" #left=below, right=above
                push!(plotArray,Plots.Node(label_, xs[t], ys[t],
                    style="rotate=90,scale=$(label_scale),$dir=2mm,fill=white,rectangle,rounded corners=3pt"))
                prev_label = label

            end
        end
    end
    plotArray
end

function closest_is_above(d::TrajLog, own_id::Int64)
    wm = d["WorldModel_$(own_id)"]
    h1 = convert(Array, wm[:z])
    hs = copy(h1)
    for i = 1:get_num_aircraft(d)
        if i != own_id
            wm_i = d["WorldModel_$i"]
            h2 = convert(Array, wm_i[:z])
            hs = hcat(hs, h2-h1)
        end
    end
    v, i = findmin(abs(hs), 2)
    minvals = squeeze(hs[i], 2)

    minvals .>= 0.0
end

function pplot_line(d::TrajLog, field::AbstractString,
                         xname::Symbol,
                         yname::Symbol;
                         mark_ra::Bool = true,
                         fx::Function = identity,
                         fy::Function = identity)

    plotArray = Plots.Plot[]
    for i = 1:get_num_aircraft(d)
        df = get_log(d, field, i)

        #plot trajectories
        xvals = convert(Array, df[xname])
        yvals = convert(Array, df[yname])

        # apply user function transforms
        xvals = fx(xvals)
        yvals = fy(yvals)

        # apply time filters
        # TODO:...

        push!(plotArray,Plots.Linear(xvals, yvals,
            style="mark options={color=blue}", mark="*"))

        # RA markings
        if mark_ra
            #mark times of RA active
            t_style_array = get_ra_style(d, i)

            for (times, style) = t_style_array
                if !isempty(times)
                    push!(plotArray, Plots.Scatter(xvals[times], yvals[times], style = style))
                end
            end

            #mark times where pilot was following RA
            t_style_array = get_response_style(d, i)

            for (times, style) = t_style_array
                if !isempty(times)
                    push!(plotArray,Plots.Scatter(xvals[times], yvals[times], style = style))
                end
            end
        end
    end
    plotArray
end

function psidot_from_psi(psi::Vector{Float64})
    # Numerical differentiation
    psid = [angle_diff_rad(x1, x0) for (x0,x1) in partition(psi, 2, 1)]
    vcat(psid, psid[end]) #repeat last element so output is same length
end

function get_ra_style(d::TrajLog, aircraft_number::Int64)
    i = aircraft_number
    cas = d["CAS_$i"]
    ra_active = cas[:ra_active]
    target_rate = cas[:ownOutput_target_rate]
    t_style_array = Tuple{Vector{Int64}, ASCIIString}[]
    for (f, s) in RA_STYLE_MAP
        #find times where f is valid and tag it with style s
        times = find(x->f(x[1], x[2]), zip(ra_active, target_rate))
        push!(t_style_array,(times, s))
    end
    t_style_array #vector of tuples.  each tuple = (times::Vector, style::AbstractString)
end

function get_response_style(d::TrajLog, aircraft_number::Int64)
    i = aircraft_number
    pr = get_log(d, "Response", i)

    #stochastic linear case
    if haskey(pr, :response)
        t_style_array = (Vector{Int64}, string)[]

        for (f, s) in RESPONSE_STYLE_MAP
            times = find(f, pr[:response])
            push!(t_style_array,(times, s))
        end
        return t_style_array
    end

    #deterministic PR case
    if haskey(pr, :state)
        t_style_array = Tuple{Vector{Int64}, ASCIIString}[]

        for (f, s) in RESPONSE_STYLE_MAP
        times = find(f, pr[:state])
        push!(t_style_array,(times, s))
        end
        return t_style_array
    end
end

function trajPlot{T<:AbstractString}(savefiles::Vector{T}; format::Symbol=:TEXPDF)
    map(f->trajPlot(f, format=format), savefiles)
end

function trajPlot(savefile::AbstractString; format::Symbol=:TEXPDF)
    # add suppl info and reload.  This avoids adding suppl info to all files
    add_supplementary(savefile)
    d = trajLoad(savefile)

    outfileroot = getSaveFileRoot(savefile)
    trajPlot(outfileroot, d, format=format)

    savefile
end

function trajPlot(outfileroot::AbstractString, d::TrajLog; format::Symbol=:TEXPDF)
    #workaround, Tikzpictures/lualatex doesn't handle backslashes properly in path
    outfileroot = replace(outfileroot, "\\", "/")

    td = TikzDocument()
    tps = pgfplotLog(d)

    cap = string(vis_runtype_caps(d),
               vis_sim_caps(d),
               vis_runinfo_caps(d))

    add_to_document!(td, tps, cap)
    plot_tikz(outfileroot, td, format)
    td
end

end #module

ext: jl

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using NeXLMatrixCorrection, NeXLSpectrum
using PrettyTables

det = BasicEDS(4096, 0.0, 5.0, 132.0, 1, Dict(KShell => n"He", LShell => n"K", MShell => n"Cs", NShell => n"Pu"))

open(joinpath(homedir(), "Desktop", "suitable_references.tex"),"w") do io
    for el in elements[5:92]
        println(io,"\\subsection{Element: $(name(el))}")
        pretty_table(io,
            suitability(el, det, minC=0.01, latex=true),
            nosubheader=true, backend=:latex,
            label="The suitability of references for $(name(el))",
            wrap_table=false, table_type=:longtable )
        println(io)
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
# CUDA stream management

immutable CuStream
    handle::Ptr{Void}
    blocking::Bool
    priority::Int
end

function null_stream()
    CuStream(convert(Ptr{Void}, 0), true, 0)
end

function destroy(s::CuStream)
    @cucall(cuStreamDestroy, (Ptr{Void},), s.handle)
end

function synchronize(s::CuStream)
    @cucall(cuStreamSynchronize, (Ptr{Void},), s.handle)
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

import ACE

import ACE: ScalarBasis, OneParticleBasis, OnepBasisFcn,
              PIBasis, PIPotential, PIBasisFcn,
              VecOrTup, get_basis_spec,
              AbstractDegree, degree,
              gensparse,
              site_alloc_B, site_alloc_dB,
              site_evaluate!, site_evaluate_d!,
              set_Aindices!, add_into_A!, add_into_A_dA!,
              scaling


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Distributions.call, (Type{Distributions.MvNormal}, Array{Float64, 1}, PDMats.PDMat{Float64, Array{Float64, 2}},))
    precompile(Distributions.call, (Type{Distributions.MvNormal}, Array{Float64, 1}, PDMats.PDiagMat{Float64, Array{Float64, 1}},))
    precompile(Distributions.call, (Type{Distributions.MvNormal}, Array{Float64, 1}, Array{Float64, 2},))
end

max_stars_repo_name: waTeim/gpsFlow

max_issues_repo_name: waTeim/flying-monkey

max_forks_repo_name: waTeim/flying-monkey


text:
## Loading packages.
using MicroHH
using Tullio
using Statistics
using GLMakie

## Loading settings.
include("drycbl_settings.jl")

## Initialize the model.
n_domains = 1
m = Model("drycbl", n_domains, settings, Float32)
load_model!(m)
in_progress = prepare_model!(m)

## Make the figure.
f = m.fields[1]; g = m.grid[1]
x = @view g.x[g.is:g.ie]
y = @view g.y[g.js:g.je]
z = @view g.z[g.ks:g.ke]
s_bot = @view f.s_bot[g.is:g.ie, g.js:g.je]
s = @view f.s[g.is:g.ie, g.js:g.je, g.ks:g.ke]

fig = Figure(resolution=(800, 500))
node1 = Node(s_bot .- mean(s_bot))
h1 = heatmap(fig[1, 1], x, y, node1, colorrange=(-0.4, 0.4))
node2 = Node(s[:, 1, :] .- mean(s, dims=(1, 2))[:, 1, :])
h2 = heatmap(fig[1, 2], x, z, node2, colorrange=(-0.4, 0.4))
fig

## Run the model.
while in_progress
    global in_progress = step_model!(m)
    node1[] = s_bot .- mean(s_bot)
    node2[] = s[:, 1, :] .- mean(s, dims=(1, 2))[:, 1, :]
end


text:
using JuMP
const MOIT = MOI.Test

function bridged_mock(mock_optimize!::Function...;
                      model = JuMP._MOIModel{Float64}())
    mock = MOI.Utilities.MockOptimizer(model)
    bridged = MOI.Bridges.full_bridge_optimizer(mock, Float64)
    MOI.Utilities.set_mock_optimize!(mock, mock_optimize!...)
    return bridged
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.7.1
#   kernelspec:
#     display_name: Julia 1.2.0
#     language: julia
#     name: julia-1.2
# ---

# # School network model
# This model uses matrix expressions of school network to provide a samplewise likelihood function for primary school influenza data in Matsumoto city, Japan.

# ## Load and prepare data

# +
# School incidence data
using Dates, Distributions
const diZero=Dates.value(Date("2014-09-30")) # day integer for t=0
const diEnd=Dates.value(Date("2015-05-01")) # day integer for t=end
const diLen=diEnd-diZero
function d2di(d::Date) Dates.value(d)-diZero end
function di2d(di::Int) Date(Dates.UTInstant(Day(di+diZero))) end

winterbreak=d2di(Date("2014-12-27")):d2di(Date("2015-1-7")) # 27 Dec - 7 Jan was winter break in school in Matsumoto
duringbreak=(x->x.+1).(winterbreak);
logispar_est=[64.21977, 113.60699, 10.15022]
# -

# ## Construct school data

# +
include("src/SchoolOutbreak.jl")

# override function for recall bias adjustment
function SchoolOutbreak.updatesusinf!(student::Student,parameters::NT where NT<:NamedTuple)
    if student.isinfected
        student.susceptibility.=linregval(student.suscovar,parameters.suscoef)
        student.infectiousness.=linregval(student.infcovar,parameters.infcoef)
    else
        @views student.susceptibility.= exp(student.suscovar[1]*parameters.suscoef[1][]+
            sum(log.( (1 .-student.suscovar[2:end]) .+ (student.suscovar[2:end].*exp.(first.(parameters.suscoef[2:end]))) )))
        @views student.infectiousness.= exp(student.infcovar[1]*parameters.infcoef[1][]+
            sum(log.( (1 .-student.infcovar[2:end]) .+ (student.infcovar[2:end].*exp.(first.(parameters.infcoef[2:end]))) )))

    end
end
# -

using JLD2
@load "Matsumoto_studentsdata.jld2" students
@time schooldata=SchoolOutbreak.StudentsData.([filter(x->x.stratum[1]==schoolid,students) for schoolid in 1:29],Ref(parameters),Ref((schoolclosure=duringbreak,)));

# ## MLE (for sensitivity analysis and model selection only)

# +
function nll(lparms,schooldata=schooldata)
    parms=copy(lparms)
    parms[1:8].=exp.(lparms[1:8])
    @views SchoolOutbreak.updateArray.(schooldata[1].parameters.β,[0;parms[1:4]])
    @views SchoolOutbreak.updateArray.(schooldata[1].parameters.γ,[0;parms[[5:7;7]]])
    @views SchoolOutbreak.updateArray.(schooldata[1].parameters.γ,[0;parms[[8:10;10]]])
    @views SchoolOutbreak.updateArray.(schooldata[1].parameters.rcom,parms[11])
    @views SchoolOutbreak.updateArray(schooldata[1].parameters.suscoef,length(parms)>12&&parms[12:15])
    @views SchoolOutbreak.updateArray(schooldata[1].parameters.infcoef,parms[end-(length(parms)-12)÷2:end])
    -SchoolOutbreak.llfunc!(schooldata)
end
@time @show nll(fill(-1.,20),schooldata)

# Baseline
parmlen=11+length(sus_covlabels)+length(inf_covlabels)

# No covariate effects
#parmlen=12

# MLE
#@time opt=optimize(nll,fill(-20.0,parmlen),fill(0.0,parmlen),fill(-4.0,parmlen))
#exp.(opt.minimizer)
# -

# LAME
#using Calculus
#@time hess=Calculus.hessian(nll,opt.minimizer)
#@show eigen(hess).values
#logdetH=log(det(hess))
#@show LAME = 2opt.minimum-length(opt.minimizer)*log(2pi)+logdetH

# -

# ## MCMC

using Mamba
function schoolmodel(schooldata,
        usecovars=(sus=zeros(Bool,length(sus_covlabels)),inf=[zeros(Bool,length(inf_covlabels)-1);true]),
        ;likelihood! = SchoolOutbreak.llfunc!,
        regularize=(x)->0.0,
        invtemp=1.0,
        dimγ=3,dimδ=3,
        γfix=1.0,δfix=1.0)

    dotter=[100]
    inputs=Dict{Symbol,Any}(
        :schooldata=>schooldata,
        #:parameters=>parameters,
        :zerotrick=>0.0,
        :invtemp=>invtemp,
        :counter=>([0],dotter),
        :usecovars=>usecovars
    )
    parms=Dict{Symbol,Any}(
        :logβs=>.-ones(4).*5,
        :γs=>ones(dimγ)./2,
        :δs=>ones(dimδ)./2,
        :logrcom=>fill(-4.0),
        :suscoef=>zeros(sum(usecovars.sus)),
        :infcoef=>zeros(sum(usecovars.inf))
    )
    if(length(parms[:suscoef])==0) push!(inputs,:suscoef=>pop!(parms,:suscoef)) end
    if(length(parms[:infcoef])==0) push!(inputs,:infcoef=>pop!(parms,:infcoef)) end
    if(length(parms[:γs])==0) push!(inputs,:γs=>push!(pop!(parms,:γs),γfix)) end
    if(length(parms[:δs])==0) push!(inputs,:δs=>push!(pop!(parms,:δs),δfix)) end
    priors=Dict{Symbol,Any}()
    for parname in keys(parms)
        priors[parname]=Stochastic(length(size(parms[parname])),()->Uniform(-20,20))
    end

    inits=merge(parms,inputs)
    inits=[inits]

    model=Model(
        βs=Logical(1,(logβs)->exp.(logβs)),
        rcom=Logical((logrcom)->exp.(logrcom)),
        llfunc=Logical((βs,γs,δs,rcom,suscoef,infcoef,usecovars,schooldata)->begin
            parameters=schooldata[1].parameters
            SchoolOutbreak.updateArray.(parameters.γ,[γs[1];γs.+zeros(3);γs[end]]) # to enable a single common γs
            SchoolOutbreak.updateArray.(parameters.δ,[γs[1];δs.+zeros(3);δs[end]]) # to enable a single common γs
            SchoolOutbreak.updateArray.(parameters.β,[0;βs])
            SchoolOutbreak.updateArray.(parameters.rcom,rcom)
            SchoolOutbreak.updateArray(parameters.suscoef[usecovars.sus],suscoef)
            SchoolOutbreak.updateArray(parameters.suscoef[.!usecovars.sus],0.0)
            SchoolOutbreak.updateArray(parameters.infcoef[usecovars.inf],infcoef)
            SchoolOutbreak.updateArray(parameters.infcoef[.!usecovars.inf],0.0)
            ll=likelihood!(schooldata)
            ll-regularize(sum((suscoef).^2)+sum((infcoef).^2))
        end, false),
        llvalue=Logical((llfunc,invtemp)->llfunc*invtemp),
        zerotrick=Stochastic((llvalue)->Poisson(-llvalue),false),
        count=Logical((counter,rcom)->0);
        priors...
    )
    countup=Sampler([:count],
        (count,counter)->begin
            counter[1].+=1
            if counter[1][1]==counter[2][1]
                counter[1].=0
                print(".")
            end
            count+=1
        end
    )
    setsamplers!(model,[AMM(collect(keys(parms)),Matrix{Float64}(I,fill(length(reduce(vcat,values(parms))),2)...).*0.0005),countup])
    return(model=model,inputs=inputs,inits=inits)
end

# +
usecovars=(sus=ones(Bool,length(sus_covlabels)),inf=ones(Bool,length(inf_covlabels))) # include all covars
schoolmodel1 = schoolmodel(schooldata[Not([3,5,21])],usecovars) # exclude school 3, 5, 21
mcmclen=200000
burn=div(mcmclen,2)

chain1 = mcmc(schoolmodel1..., mcmclen, burnin=burn, thin=max(1,div(mcmclen,2000)), chains=1, verbose=true);
# -

Base.exp(chain::ModelChains)=begin newchain=deepcopy(chain);newchain.value.=exp.(chain.value);newchain end
showparam=[:logβs,:γs,:logrcom,:infcoef]
@show Mamba.draw(Mamba.plot(chain1[:,showparam,:]))
describe(exp(chain0[:,[:logβs,:logrcom,:infcoef],:]))
describe(chain1[:,[:γs,:llvalue],:])

write("Matsumoto_chain1.jls", chain1)


text:
using PDMats  # Positive Definite Matrices
using QuadGK  # 1-d integration
using Optim

# Read dataset
data = readdlm("../data/rvs_0001.txt");
times = data[:,1]
vels = data[:,2]
sigma_obs = data[:,3]
nobs = length(vels)

#
function compute_Fischer_matrix_offset{T<:Real}(data::Vector{T}, model::Array{T,2}, Sigma::PDMat{T})
  @assert length(data) == size(model,1) == size(Sigma,1)
  nparam = size(model,2)
  FIM = Array{Float64}(nparam,nparam)
  FIM[1,1] = sum(1.0./diag(Sigma))
  return FIM
end

function compute_Fischer_matrix_1sin{T<:Real}(data::Vector{T}, model::Array{T,2}, Sigma::PDMat{T})
  @assert length(data) == size(model,1) == size(Sigma,1)
  nparam = size(model,2)
  FIM = Array{Float64}(nparam,nparam)
  FIM[1,1] = sum(1.0./diag(Sigma))
  return FIM
end

const Cmax = 1000.0
function logprior_offset(param)
  -log(2*Cmax)
end

const Kmax = 999.0
const K0 = 1.0
function logprior_sinusoids(param)
  nparam = length(param)
  nsin = convert(Int64,floor((nparam-1)//2))
  logp = 0
  for i in 0:(nsin-1)
    K = sqrt(param[1+2*i]^2+param[2+2*i]^2)
    logp -= log1p(K/K0)+log(log1p(Kmax/K0))
  end
  C = param[nparam]
  if !(-Cmax<=C<=Cmax) logp -=Inf
  else  logp -= log(2*Cmax) end
  return logp
end

function calc_best_fit(data::Vector{T}, model::Array{T,2}, Sigma::PDMat{T}) where T<:Real
  inv_chol_Sigma = inv(Sigma.chol[:L])
  param_bf_linalg = (inv_chol_Sigma*model) \ (inv_chol_Sigma*data)
end

function calc_best_fit(data::Vector{T}, model::Array{T,2}, Sigma::PDiagMat{T}) where T<:Real
  inv_chol_Sigma = sqrt.(Sigma.inv_diag)
  param_bf_linalg = (inv_chol_Sigma.*model) \ (inv_chol_Sigma.*data)
end

function compute_Laplace_Approx(data::Vector{T}, model::Array{T,2}, Sigma::SigmaT; FIM::PDMat{T} = PDMat(Xt_invA_X(Sigma,X)), calclogprior::Function = x->0.0 ) where {T<:Real, SigmaT<:AbstractPDMat{T} }
  @assert length(data) == size(model,1) == size(Sigma,1)
  nobs = length(data)
  nparam = size(model,2)

  param_bf_linalg = calc_best_fit(data,model,Sigma)
  predict = X*param_bf_linalg

  #FIM = compute_Fischer_matrix_offset(data,model,Sigma)
  #FIM = compute_Fischer_matrix_1sin(data,model,Sigma)
  #chol_FIM = chol(FIM)
  #chol_FIM = chol(full(FIM))
  #sigma_param = inv(chol_FIM)
  sigma_param = inv(FIM).chol[:U]

  delta = data.-predict
  chisq = invquad(Sigma,delta) # sum(delta'*Sigma^-1*delta)

  loglikelihood_mode = 0.5*(-chisq-logdet(Sigma)-nobs*log(2pi))
  logprior_mode = calclogprior(param_bf_linalg)
  log_mode = logprior_mode + loglikelihood_mode
  LaplaceApprox = 0.5*nparam*log(2pi)-0.5*logdet(FIM)+log_mode
  #LaplaceApprox = 0.5*nparam*log(2pi)-logdet(chol_FIM)+log_mode
  return (param_bf_linalg,sigma_param,LaplaceApprox)
end

function generate_K_samples(bf::Vector, sigma, nsamples::Integer = 1)
  nparam = length(bf)
  theta = bf .+ sigma*rand(nparam,nsamples)
  #theta = param_bf'.+rand(nsamples,nparam)*sigma_param
  namps = convert(Int64,floor((nparam-1)//2))

  Ks = zeros(nsamples,namps)
  for i in 1:namps
    Ks[:,i] = sqrt.(theta[1+2*(i-1),:].^2+theta[2+2*(i-1),:].^2)
  end
  return Ks
end

nparam = 1
X = ones(nobs,nparam);
Sigma=PDiagMat(data[:,3].^2)
compute_Laplace_Approx(vels,X,Sigma)

# Simplified zero planet model to test linear algebra
# Brute force best-fit C and sigma for comparison
C_bf = sum(data[:,2]./data[:,3].^2)/sum(data[:,3].^(-2))
sigma_bf_0 = sqrt(1/sum(data[:,3].^(-2)))

# Now with covariances
function kernel_quaesiperiodic(dt::Float64)
  alpha_sq     = 3.0
  tau_evolve   = 50.0
  lambda_p     = 0.5
  tau_rotate   = 20.0
  alpha_sq*exp(-0.5*((dt/tau_evolve)^2 + (sin(pi*dt/tau_rotate)/lambda_p)^2) )
end

function make_Sigma(t::Vector{Float64}, sigma_obs::Vector{Float64}, sigma_j::Float64)
  PDMat( [ kernel_quaesiperiodic(t[i]-t[j]) + (i==j ? sigma_obs[i]^2 + sigma_j^2 : 0.0) for i in 1:length(t), j in 1:length(t)] )
end

function make_dSigmadsigmaj(t::Vector{Float64}, sigma_obs::Vector{Float64}, sigma_j::Float64)
  PDMat( [ 0.0 + (i==j ? 2*sigma_j : 0.0) for i in 1:length(t), j in 1:length(t)] )
end

const sigma_j_min = 0.0
const sigma_j_max = 99.0
const sigma_j_0 = 1.0
function logprior_jitter(sigma_j::Real)
  -log1p(sigma_j/sigma_j_0)-log(sigma_j_0)-log(log1p(sigma_j_max/sigma_j_0))
end

function compute_log_evidence(sigma_j::Real)
  Sigma = make_Sigma(times,sigma_obs,sigma_j)
  log_prior = logprior_jitter(sigma_j) + logprior_offset([])
  log_likelihood = compute_Laplace_Approx(vels,X,Sigma)[3]
  log_prior+log_likelihood
end

function compute_evidence(sigma_j::Real; log_normalization::Real=0.0)
  exp(compute_log_evidence(sigma_j)-log_normalization)
end

#lognorm = compute_log_evidence(sigma_j_0)  # needed if underflow
integrand(x) = compute_evidence(x) # ,log_normalization=0.0)
result = QuadGK.quadgk(integrand,sigma_j_min,sigma_j_max)./log(10)

function calc_opt_sigmaj(data::Vector{T}, model::Array{T,2}) where T<:Real
  function calc_laplace_approx_arg_jitter(sigmaj::Real)
  Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
    FIM = PDMat(Xt_invA_X(Sigma,X))
    (param_bf,sigma_param,evid) = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
	evid += logprior_jitter(sigmaj)
    -evid
  end
  result = optimize(calc_laplace_approx_arg_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-6)
end
res = calc_opt_sigmaj(vels,X)
dSigmadsigmaj = make_dSigmadsigmaj(times,sigma_obs,res.minimizer)
evid = (compute_log_evidence(res.minimizer)+0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2))))/log(10)

# Brute force 2-d integration
using Cubature

function log_target(x::Vector)
  offset = x[1]
  sigma_j = x[2]
  logprior = logprior_offset([offset]) + logprior_jitter(sigma_j)
  Sigma = make_Sigma(times,sigma_obs,sigma_j)
  delta = vels.-offset
  chisq = invquad(Sigma,delta)
  loglikelihood = 0.5*(-chisq -logdet(Sigma)-nobs*log(2pi))
  logprior + loglikelihood
end

function integrand(x::Vector; log_normalization::Real=0.0)
  exp(log_target(x)-log_normalization)
end
@time hcubature(integrand, [-Cmax/100,sigma_j_min],[Cmax/100,sigma_j_max],reltol=1e-2)
#@time pcubature(integrand, [-Cmax/100,sigma_j_min],[Cmax/100,sigma_j_max],reltol=1e-2)

using Cuba
function integrand_vegas(x::Vector, output::Vector; log_normalization::Real=0.0)
  xx = zeros(2)
  xx[1] = x[1]*2*Cmax/100-Cmax/100
  xx[2] = x[2]*(sigma_j_max-sigma_j_min)+sigma_j_min
  output[1] = exp(log_target(xx)-log_normalization)*(2*Cmax/100)*(sigma_j_max-sigma_j_min)
end
@time vegas(integrand_vegas,2,1,reltol=1e-3,maxevals=1e5)
@time suave(integrand_vegas,2,1,reltol=1e-3,maxevals=1e5)
@time divonne(integrand_vegas,2,1,reltol=1e-3,maxevals=1e5)
@time cuhre(integrand_vegas,2,1,reltol=1e-3,minevals=1e3,maxevals=1e5)

#=
# Brute Force:  Sequence of Nested 1-d integrations
function compute_log_target_fixed_offset(offset::Real)
  delta = vels.-offset
  function integrand(sigma_j::Real)
    logprior = logprior_offset([offset]) + logprior_jitter(sigma_j)
    Sigma = make_Sigma(times,sigma_obs,sigma_j)
    chisq = invquad(Sigma,delta)
    loglikelihood = 0.5*(-chisq -logdet(Sigma)-nobs*log(2pi))
    logtarget = logprior + loglikelihood
	exp(logtarget) # +log(1e210))
  end
  result = QuadGK.quadgk(integrand,sigma_j_min,sigma_j_0,sigma_j_max,reltol=1e-2)[1]
end
@time result = QuadGK.quadgk(x->begin p=compute_log_target_fixed_offset(x); println(x," ",p); p end,-Cmax/100,Cmax/100,reltol=1e-2)./log(10)

# Brute Force: 2-d integration
function log_target(x::Vector)
  offset = x[1]
  sigma_j = x[2]
  logprior = logprior_offset([offset]) + logprior_jitter(sigma_j)
  Sigma = make_Sigma(times,sigma_obs,sigma_j)
  delta = vels.-offset
  chisq = invquad(Sigma,delta)
  loglikelihood = 0.5*(-chisq -logdet(Sigma)-nobs*log(2pi))
  logprior + loglikelihood
end
=#

# Fit planet models
function make_design_matrix(times::Vector{T},periods::Vector{T}) where T<:Real
  nparam = 1+2*length(periods)
  freqs = 2pi./periods
  X = Array{T}(length(times),nparam)
  for i in 1:length(periods)
     X[:,1+(i-1)*2] = sin(freqs[i]*times)
	 X[:,2+(i-1)*2] = cos(freqs[i]*times)
  end
  X[:,nparam] = 1.0
  return X
end

# Fit one planet model
periods = [42.0]
nparam = 1+2*length(periods)
freqs = 2pi./periods
X = [sin(freqs[1]*times) cos(freqs[1]*times) ones(times)];
Sigma = PDMat(make_Sigma(times,sigma_obs,sigma_j_0))
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)

function calc_opt_sigmaj(data::Vector{T}, model::Array{T,2}) where T<:Real
  function calc_laplace_approx_arg_jitter(sigmaj::Real)
	Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
    FIM = PDMat(Xt_invA_X(Sigma,X))
    (param_bf,sigma_param,evid) = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
    evid += logprior_jitter(sigmaj)
	-evid
  end
  result = optimize(calc_laplace_approx_arg_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-4)
end
res = calc_opt_sigmaj(vels,X)
Sigma = PDMat(make_Sigma(times,sigma_obs,res.minimizer))
dSigmadsigmaj = make_dSigmadsigmaj(times,sigma_obs,res.minimizer)
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
evid += 0.5*(log(2pi)-log(0.5*trace((Sigma \ full(dSdj))^2)))

Ks = generate_K_samples(param_bf,sigma_param,1000);
m = mean(Ks,1);
sig = std(Ks,1,mean=m);
m-sig,m,m+sig

function profile_sigmaj(data::Vector{T}, model::Array{T,2}) where T<:Real
  function calc_laplace_approx_arg_jitter(sigmaj::Real) #param::Vector{T})
    #sigmaj = param[1]
	Sigma = PDMat(make_Sigma(times,sigma_obs,sigmaj))
    FIM = PDMat(Xt_invA_X(Sigma,X))
    (param_bf,sigma_param,evid) = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_sinusoids)
	evid += logprior_jitter(sigmaj)
    evid
  end
  #result = optimize(calc_laplace_approx_arg_jitter,[sigma_j_0])
  result = optimize(calc_laplace_approx_arg_jitter,sigma_j_min,sigma_j_max,rel_tol=1e-6)
end
res = calc_opt_sigmaj(vels,X)

# Fit two planet model
periods = [42.0,12.0]
X = make_design_matrix(times,periods)
Sigma = PDMat(make_Sigma(times,sigma_obs,0.68442348767604010984))
FIM = PDMat(Xt_invA_X(Sigma,X))
(param_bf,sigma_param,evid) = compute_Laplace_Approx(vels,X,Sigma,FIM=FIM,calclogprior=logprior_1sin)

Ks = generate_K_samples(param_bf,sigma_param,1000)
m = mean(Ks,1)
sig = std(Ks,1,mean=m)
m-sig,m+sig


text:
using MCMCChains
using Turing
using Test

@testset "get tests" begin
    @model model(x) = begin
        m = Vector{Real}(undef, 10)
        [m[i] ~ Normal(i, 0.1) for i in 1:10]
        s ~ InverseGamma(2,3)
        x ~ Normal(0, s)
    end

    n_samples = 1000
    model = model(2.0)
    sampler = MH(n_samples)
    chn = sample(model, sampler)

    get1 = get(chn, :m)
    get2 = get(chn, [:m, :s])

    @test length(get1.m) == 10
    @test length(get1.m[5]) == n_samples
    @test round(mean(get1.m[1])) ≈ 1.0
    @test round(mean(get1.m[10])) ≈ 10.0
    @test length(get2.s) == n_samples
    @test collect(keys(get2)) == [:m, :s]

    getsection = get(chn, section = :internals)
    @test length(keys(getsection)) == 2

    getall = get_params(chn)
    @test getall == get(chn, section=[:internals, :parameters])
    @test length(keys(getall)) == 4
end


text:
using ACE, Test, ACEbase, ACEbase.Testing, StaticArrays
using ACE.Random: rand_rot, rand_refl
using Random: shuffle

##
#=
@info("Rudimentary tests for sparse basis selectors and intersections of such")
r0cut = 2.0
rcut = 1.0
zcut = 2.0
env = ACE.EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=false, λ= .5)

maxorder = 3
dmaxdeg = 4
Bsel_p2 = ACE.SparseBasis(maxorder; p = 2, default_maxdeg = dmaxdeg)
Bsel_p1 = ACE.SparseBasis(maxorder; p = 1, default_maxdeg = dmaxdeg)

Bsel_intesect = ACE.intersect(Bsel_p1,Bsel_p2)
B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=dmaxdeg)
basis_intersect_inv = ACE.SymmetricBasis(ACE.Invariant(), B1p, Bsel_intesect)
=#

@info("Test bond basis selectors")

r0cut = 2.0
rcut = 1.0
zcut = 2.0
env = ACE.EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=false, λ= .5)

maxorder = 3
Bsel = ACE.SparseBasis(; maxorder=maxorder, p = 2, default_maxdeg = 4)

##

@info("Test invariance")

basis_inv = ACE.Utils.SymmetricBond_basis(ACE.Invariant(), env, Bsel; )
@show length(basis_inv)

rr0 = SVector{3}(rand(Float64,3))
cfg = [ ACE.State(rr = SVector{3}(rand(Float64,3)), rr0 = rr0,
                  be = rand([:bond,:env]))
        for _ = 1:10 ] |> ACEConfig
B1_inv = ACE.evaluate(basis_inv, cfg)

for ntest = 1:30
    Q = rand_refl() * rand_rot()
    Xs2 = ACE.shuffle([ ACE.State(rr = Q * X.rr, rr0 = Q * X.rr0, be = X.be)  for X in cfg.Xs ])
    B2_inv = ACE.evaluate(basis_inv, ACEConfig(Xs2))
    print_tf(@test isapprox(B1_inv, B2_inv, rtol=1e-10))
end
println()
@info("Test Euclidian covariance")

basis_cov = ACE.Utils.SymmetricBond_basis(ACE.EuclideanVector(), env, Bsel; )
@show length(basis_cov)

rr0 = SVector{3}(rand(Float64,3))
cfg = [ ACE.State(rr = SVector{3}(rand(Float64,3)), rr0 = rr0,
                  be = rand([:bond,:env]))
        for _ = 1:10 ] |> ACEConfig
B1_cov = ACE.evaluate(basis_cov, cfg)

for ntest = 1:30
    Q = rand_refl() * rand_rot()
    Xs2 = ACE.shuffle([ ACE.State(rr = Q * X.rr, rr0 = Q * X.rr0, be = X.be)  for X in cfg.Xs ])
    B2_cov = ACE.evaluate(basis_cov, ACEConfig(Xs2))
    print_tf(@test isapprox( map(x->Q*x, B1_cov), B2_cov, rtol=1e-10))
end
println()


text:
using FastTransforms, FastGaussQuadrature, Test

hermitepoints(n) = FastGaussQuadrature.unweightedgausshermite( n )[1]

@testset "Hermite" begin
    @test hermitepoints(1) == [0.0]
    @test hermitepoints(100_000)[end] ≈ 446.9720305443094

    @test weightedhermitetransform([1.0]) == [1.0]
    @test weightedhermitetransform(exp.(-hermitepoints(2).^2/2)) == [1.0,0.0]
    @test weightedhermitetransform(exp.(-hermitepoints(3).^2/2)) ≈ [1.0,0.0,0.0]
    @test weightedhermitetransform(exp.(-hermitepoints(1000).^2/2)) ≈ [1.0; zeros(999)]
    @test weightedhermitetransform(exp.(-hermitepoints(3000).^2/2)) ≈ [1.0; zeros(2999)]

    for n in (1, 5,100)
        x = randn(n)
        @test iweightedhermitetransform(weightedhermitetransform(x)) ≈ x
        @test weightedhermitetransform(iweightedhermitetransform(x)) ≈ x
    end

    x = hermitepoints(100)
    @test iweightedhermitetransform([0.0; 1.0; zeros(98)]) ≈ (exp.(-x.^2 ./ 2) .* 2x/sqrt(2))
    @test iweightedhermitetransform([0.0; 0; 1.0; zeros(97)]) ≈ (exp.(-x.^2 ./ 2) .* (4x.^2 .- 2)/(sqrt(2)*2^(2/2)))
    @test iweightedhermitetransform([0.0; 0; 0; 1.0; zeros(96)]) ≈ (exp.(-x.^2 ./ 2) .* (-12x + 8x.^3) / (sqrt(2*3)*2^(3/2)))
    @test iweightedhermitetransform([0.0; 0; 0; 0; 1.0; zeros(95)]) ≈ (exp.(-x.^2 ./ 2) .* (12 .- 48x.^2 .+ 16x.^4) / (sqrt(2*3*4)*2^(4/2)))
end


text:
const FI_NOT_IMPLEMENTED = [Weibull; BetaFourParameters; NegativeBinomial]

struct FakeDist{T<:Real} <: Distributions.ContinuousUnivariateDistribution
    foo::T
end

function instantiate_dist(D::Type{<:Distribution})
    if D == ScoreDrivenModels.TDistLocationScale
        params = [1.0 2.0 10]
        return ScoreDrivenModels.update_dist(D, params, 1)
    elseif D == ScoreDrivenModels.TDist
        params = [10][:, :]
        return ScoreDrivenModels.update_dist(D, params, 1)
    elseif D == ScoreDrivenModels.BetaFourParameters
        params = [0.0 10.0 2.0 3]
        return ScoreDrivenModels.update_dist(D, params, 1)
    else
        params = 0.5 * ones(1, ScoreDrivenModels.num_params(D))
        return ScoreDrivenModels.update_dist(D, params, 1)
    end
end

function test_score_mean(D::Type{<:Distribution}; n::Int = 10^7, seed::Int = 10,
                            atol::Float64 = 1e-3, rtol::Float64 = 1e-3)
    Random.seed!(seed)
    dist = instantiate_dist(D)
    pars = permutedims([ScoreDrivenModels.params_sdm(dist)...])
    n_params = ScoreDrivenModels.num_params(D)
    avg  = zeros(1, n_params)
    obs = rand(dist, n)
    score_aux = ones(1, n_params)
    for i = 1:n
        ScoreDrivenModels.score!(score_aux, obs[i], D, pars, 1)
        avg .+= score_aux
    end
    avg ./= n
    @test avg ≈ zeros(1, ScoreDrivenModels.num_params(D)) atol = atol rtol = rtol
end

function test_fisher_information(D::Type{<:Distribution}; n::Int = 10^6, seed::Int = 10,
                                 atol::Float64 = 1e-2, rtol::Float64 = 1e-2)
    Random.seed!(seed)
    dist = instantiate_dist(D)
    pars = permutedims([ScoreDrivenModels.params_sdm(dist)...])
    n_params = ScoreDrivenModels.num_params(D)
    var_terms  = zeros(n_params, n_params)
    obs = rand(dist, n)
    score_aux = ones(1, n_params)
    for i = 1:n
        ScoreDrivenModels.score!(score_aux, obs[i], D, pars, 1)
        var_terms .+= score_aux' * score_aux
    end
    var_terms ./= n
    aux_lin_alg = AuxiliaryLinAlg{Float64}(n_params)

    # Some distributions might not have the fisher information yet available
    if D in FI_NOT_IMPLEMENTED
        return
    else
        ScoreDrivenModels.fisher_information!(aux_lin_alg, D, pars, 1)
    end

    @test var_terms ≈ aux_lin_alg.fisher atol = atol rtol = rtol
end

function test_loglik(D::Type{<:Distribution}; atol::Float64 = 1e-3, rtol::Float64 = 1e-3,
                     seed::Int = 13, n::Int = 100)
    Random.seed!(seed)
    dist = instantiate_dist(D)
    y = rand(dist, n)
    pars = Matrix{Float64}(undef, n, ScoreDrivenModels.num_params(D))
    pars_dist = ScoreDrivenModels.params_sdm(dist)
    for t in axes(pars, 1), p in axes(pars, 2)
        pars[t, p] = pars_dist[p]
    end
    log_lik = ScoreDrivenModels.log_likelihood(D, y, pars, n)
    @test log_lik ≈ -loglikelihood(dist, y) atol = atol rtol = rtol
    return
end

function test_link_interfaces(D::Type{<:Distribution})
    n_params = ScoreDrivenModels.num_params(D)
    t = 1
    param_tilde = zeros(1, n_params)
    param = ones(1, n_params)
    ScoreDrivenModels.link!(param_tilde, D, param, t)
    ScoreDrivenModels.unlink!(param, D, param_tilde, t)
    @test param == ones(1, n_params)

    # No tests for jacobian_link! just pass through the function
    aux = AuxiliaryLinAlg{Float64}(n_params)
    ScoreDrivenModels.jacobian_link!(aux, D, param, t)
end

function test_dist_utils(D::Type{<:Distribution})
    # num_params
    dist = instantiate_dist(D)
    n_pars = ScoreDrivenModels.num_params(D)
    @test n_pars == length(ScoreDrivenModels.params_sdm(dist))

    # update_dist
    t = 1
    pars = permutedims([ScoreDrivenModels.params_sdm(dist)...])
    updated_dist = ScoreDrivenModels.update_dist(D, pars, t)
    @test typeof(updated_dist) <: D
end

function test_distribution_common_interface()
    score_til = ones(1, 1)
    param = ones(1, 1)
    param_tilde = ones(1, 1)
    aux = AuxiliaryLinAlg{Float64}(1)
    t = 1
    n = 1
    y_f64 = 1.0
    y_i64 = 1
    @test_throws ErrorException ScoreDrivenModels.score!(score_til, y_f64, FakeDist, param, t)
    @test_throws ErrorException ScoreDrivenModels.score!(score_til, y_i64, FakeDist, param, t)
    @test_throws ErrorException ScoreDrivenModels.fisher_information!(aux, FakeDist, param, t)
    @test_throws ErrorException ScoreDrivenModels.log_likelihood(FakeDist, [y_f64], param, n)
    @test_throws ErrorException ScoreDrivenModels.link!(param_tilde, FakeDist, param, t)
    @test_throws ErrorException ScoreDrivenModels.unlink!(param, FakeDist, param_tilde, t)
    @test_throws ErrorException ScoreDrivenModels.jacobian_link!(aux, FakeDist, param, t)
    @test_throws ErrorException ScoreDrivenModels.update_dist(FakeDist, param, t)
    @test_throws ErrorException ScoreDrivenModels.params_sdm(FakeDist(1))
    @test_throws ErrorException ScoreDrivenModels.num_params(FakeDist)
end

function test_dynamic(sigma::Float64, lags::Int; seed::Int = 12, atol::Float64 = 1e-1, rtol::Float64 = 1e-3, n::Int = 100)
    Random.seed!(seed)
    dist = Normal(0, sigma^2)
    obs_dynamic = kron(ones(n), collect(1:lags)) + rand(dist, n*lags)
    gas_lag_lag = Model(lags, lags, Normal, 0.0)
    initial_params = dynamic_initial_params(obs_dynamic, gas_lag_lag)
    for i in axes(initial_params, 1)
        @test initial_params[i, 1] ≈ i atol = atol rtol = rtol
        @test initial_params[i, 2] ≈ sigma^2 atol = atol rtol = rtol
    end
end

function simulate_GAS_1_1(D::Type{<:Distribution}, scaling::Float64, ω::Vector{T}, A::Matrix{T},
                            B::Matrix{T}, seed::Int) where T
    Random.seed!(seed)
    gas = ScoreDrivenModels.Model(1, 1, D, scaling)

    gas.ω = ω
    gas.A[1] = A
    gas.B[1] = B
    series, param = simulate_recursion(gas, 5000)

    return series
end

function simulate_GAS_1_12(D::Type{<:Distribution}, scaling::Float64, seed::Int)
    Random.seed!(seed)
    v = [0.1, 0.1]
    gas = Model([1, 12], [1, 12], D, scaling)

    gas.ω = v
    gas.A[1]  = convert(Matrix{Float64}, Diagonal(3*v))
    gas.A[12] = convert(Matrix{Float64}, Diagonal(-3*v))
    gas.B[1]  = convert(Matrix{Float64}, Diagonal(3*v))
    gas.B[12] = convert(Matrix{Float64}, Diagonal(-3*v))

    series, param = simulate_recursion(gas, 5000)

    return series
end

function test_coefficients_GAS_1_1(gas::Model{D, T}, ω::Vector{T}, A::Matrix{T}, B::Matrix{T};
                                    atol = 1e-1, rtol = 1e-1) where {D <: Distribution, T}
    @test gas.ω[1] ≈ ω[1] atol = atol rtol = rtol
    @test gas.ω[2] ≈ ω[2] atol = atol rtol = rtol
    @test gas.A[1][1, 1] ≈ A[1, 1] atol = atol rtol = rtol
    @test gas.A[1][2, 2] ≈ A[2, 2] atol = atol rtol = rtol
    @test gas.B[1][1, 1] ≈ B[1, 1] atol = atol rtol = rtol
    @test gas.B[1][2, 2] ≈ B[2, 2] atol = atol rtol = rtol
    return
end

function test_coefficients_GAS_1_12(gas::Model{D, T}; atol = 1e-1, rtol = 1e-1) where {D <: Distribution, T}
    @test gas.ω[1] ≈ 0.1 atol = atol rtol = rtol
    @test gas.ω[2] ≈ 0.1 atol = atol rtol = rtol
    @test gas.A[1][1, 1] ≈ 0.3 atol = atol rtol = rtol
    @test gas.A[1][2, 2] ≈ 0.3 atol = atol rtol = rtol
    @test gas.A[12][1, 1] ≈ -0.3 atol = atol rtol = rtol
    @test gas.A[12][2, 2] ≈ -0.3 atol = atol rtol = rtol
    @test gas.B[1][1, 1] ≈ 0.3 atol = atol rtol = rtol
    @test gas.B[1][2, 2] ≈ 0.3 atol = atol rtol = rtol
    @test gas.B[12][1, 1] ≈ -0.3 atol = atol rtol = rtol
    @test gas.B[12][2, 2] ≈ -0.3 atol = atol rtol = rtol
    return
end

function normality_quantile_and_pearson_residuals(D, n::Int, lags::Int; seed::Int = 11)
    Random.seed!(seed)

    gas = Model(1, 1, D, 0.0)
    gas.ω .= [0.0; 0.0]
    gas.A[1][[1; 4]] .= 0.2
    gas.B[1][[1; 4]] .= 0.2

    y, params = simulate_recursion(gas, n)
    quant_res = quantile_residuals(y, gas; initial_params = params[1:1, :])
    pearson = pearson_residuals(y, gas, initial_params = params[1:1, :])

    # quantile residuals
    jb = JarqueBeraTest(quant_res)
    @test pvalue(jb) >= 0.05
    Ljb = LjungBoxTest(quant_res, lags)
    @test pvalue(Ljb) >= 0.05

    # pearson
    Ljb = LjungBoxTest(pearson, lags)
    @test pvalue(Ljb) >= 0.05
    return
end

function test_GARCH_1_1(y, seed::Int, optimizer; atol = 1e-4, rtol = 1e-4)
    Random.seed!(seed)
    ini = [mean(y) var(y)]
    ub = [1.0; 1.0; 1.0; 1.0]
    lb = [-1.0; 0.0; 0.0; 0.0]
    gas = Model(1, 1, Normal, 1.0, time_varying_params = [2])
    f = fit!(gas, y; initial_params = ini, verbose = 2, opt_method = optimizer(gas, 10; ub = ub, lb = lb))
    show(stdout, fit_stats(f))

    @test gas.ω[1] - -0.00616637237701241 ≈ 0 atol = atol rtol = rtol
    @test gas.ω[2] - 0.010760592759725487 ≈ 0 atol = atol rtol = rtol
    @test gas.A[1][2, 2] - 0.15341133818189595 ≈ 0 atol = atol rtol = rtol
    @test gas.B[1][2, 2] - (0.15341133818189595 + 0.8058745318161223) ≈ 0 atol = atol rtol = rtol
    return
end


text:
"""
Generates various SearchLight files.
"""
module Generator

import Unicode, Logging
using SearchLight

include("FileTemplates.jl")
using .FileTemplates

"""
    newmodel(name::Union{String,Symbol}; path::Union{String,Nothing} = nothing, pluralize::Bool = true) :: Nothing

Generates a new SearchLight model file and persists it to the resources folder.
"""
function newmodel(name::Union{String,Symbol}; path::Union{String,Nothing} = nothing, pluralize::Bool = true) :: Nothing
  name = string(name) |> uppercasefirst
  model_name = SearchLight.Inflector.is_singular(name) ? SearchLight.Inflector.to_plural(name) : name

  model_path = setup_resource_path(model_name, path)
  mfn = model_file_name(model_name)
  write_resource_file(model_path, mfn, model_name, :model, pluralize = pluralize) &&
    @info "New model created at $(abspath(joinpath(model_path, mfn)))"

  nothing
end

"""
    newresource(resource_name::Union{String,Symbol}) :: Nothing

Generates all the files associated with a new resource and persists them to the resources folder.
"""
function newresource(resource_name::Union{String,Symbol}; path::Union{String,Nothing} = nothing, pluralize::Bool = true) :: Nothing
  resource_name = string(resource_name)

  sf = SearchLight.Inflector.tosingular(resource_name)

  model_name = (sf === nothing ? resource_name : sf) |> uppercasefirst
  newmodel(model_name, path = path, pluralize = pluralize)
  new_table_migration(resource_name, pluralize = pluralize)

  resource_name = uppercasefirst(resource_name)
  pluralize && SearchLight.Inflector.is_singular(resource_name) &&
    (resource_name = SearchLight.Inflector.to_plural(resource_name))

  resource_path = setup_resource_path(resource_name, path)
  for (resource_file, resource_type) in [(validator_file_name(resource_name), :validator)]
    write_resource_file(resource_path, resource_file, resource_name, resource_type, pluralize = pluralize) &&
      @info "New $resource_type created at $(abspath(joinpath(resource_path, resource_file)))"
  end

  isdir(SearchLight.TEST_PATH) || mkpath(SearchLight.TEST_PATH)
  test_file = resource_name * SearchLight.TEST_FILE_IDENTIFIER |> lowercase
  write_resource_file(SearchLight.TEST_PATH, test_file, resource_name, :test, pluralize = pluralize) &&
    @info "New unit test created at $(abspath(joinpath(SearchLight.TEST_PATH, test_file)))"

  nothing
end

function new_table_migration(migration_name::Union{String,Symbol}; pluralize::Bool = true) :: Nothing
  resource_name = uppercasefirst(string(migration_name))

  SearchLight.Inflector.is_singular(resource_name) && pluralize &&
    (resource_name = SearchLight.Inflector.to_plural(resource_name))

  migration_name = "create_table_" * lowercase(resource_name)
  SearchLight.Migration.new_table(migration_name, lowercase(resource_name))

  nothing
end

function newmigration(migration_name::Union{String,Symbol}) :: Nothing
  migration_name = replace(uppercasefirst(string(migration_name)) |> lowercase, " "=>"_")

  SearchLight.Migration.new(migration_name)

  nothing
end

"""
    setup_resource_path(resource_name::String) :: String

Computes and creates the directories structure needed to persist a new resource.
"""
function setup_resource_path(resource_name::String, path::Union{String,Nothing} = nothing) :: String
  resource_path = path === nothing ?
                  joinpath(SearchLight.RESOURCES_PATH, lowercase(resource_name)) :
                  path

  if ! isdir(resource_path)
    mkpath(resource_path)
    push!(LOAD_PATH, resource_path)
  end

  resource_path
end

"""
    write_resource_file(resource_path::String, file_name::String, resource_name::String) :: Bool

Generates all resouce files and persists them to disk.
"""
function write_resource_file(resource_path::String, file_name::String, resource_name::String, resource_type::Symbol; pluralize::Bool = true) :: Bool
  resource_name = (pluralize ? SearchLight.Inflector.tosingular(resource_name) : resource_name) |> SearchLight.Inflector.from_underscores

  try
    if resource_type == :model
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, FileTemplates.newmodel(resource_name, pluralize = pluralize))
      end
    end
  catch ex
    @error ex
  end

  try
    if resource_type == :validator
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, FileTemplates.newvalidator(resource_name, pluralize = pluralize))
      end
    end
  catch ex
    @error ex
  end

  try
    if resource_type == :test
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        uname = SearchLight.Inflector.from_underscores(resource_name)
        uname = pluralize ? SearchLight.Inflector.to_plural(uname) : uname

        write(f, FileTemplates.newtest(resource_name))
      end
    end
  catch ex
    @error ex
  end

  true
end

function newconfig(path::String = SearchLight.DB_PATH; filename = SearchLight.SEARCHLIGHT_DB_CONFIG_FILE_NAME) :: Nothing
  ispath(path) || mkpath(path)
  filepath = joinpath(path, filename) |> abspath
  open(filepath, "w") do io
    write(io, FileTemplates.newconfig())
  end

  @info "New config file create at $(filepath)"

  nothing
end

function resource_does_not_exist(resource_path::String, file_name::String) :: Bool
  if isfile(joinpath(resource_path, file_name))
    @warn "File already exists, $(joinpath(resource_path, file_name)) - skipping"

    return false
  end

  true
end

function model_file_name(resource_name::Union{String,Symbol})
  "$resource_name.jl"
end

function validator_file_name(resource_name::Union{String,Symbol})
  string(resource_name, SearchLight.SEARCHLIGHT_VALIDATOR_FILE_POSTFIX)
end

end


text:
function NaiveGauss(A::Array,B::Vector)

  #A = the matrix of coefficients 'A' must be squared
  #B = The Vector 'B' right hand side vector Vertical
    m,n =size(A)
    v = size(B,1)
    if m != n
        error("Matrix must be squared")
    elseif m != v
        error("Vector Length must be equal to Matrix")
    end

    # System of equations
    S=[A B]
    X=zeros(m)
    L=Matrix{Float64}(I,m,m)
    #Forward Elimination
    for j=1:m-1
       for i=j+1:m

            f=S[i,j]/S[j,j]
            S[i,:]=S[i,:].-f.*S[j,:]
            L[i,j]=f
        end
    end

    #Back Substitution
    X[m]=S[end,end]/S[end,end-1]

    for i=m-1:-1:1
        X[i]=(S[i,end] - (S[i,i+1:m]'*X[i+1:m]))/S[i,i]
    end

    return X
end


text:
# X and Z components in fT
# conducting 100 Ωm
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   20, 1, 100]
# operator
ntimesperdecade = 10
nfreqsperdecade = 5
nkᵣeval = 60
zTx = -120
zRx = -80
x_rx = -115.
y_rx = 0.
rx_roll = 0.
rx_pitch = 0.
rx_yaw = 10. #ACHTUNG!
tx_roll = 0.
tx_pitch = 0.
tx_yaw = 0.
addprimary = true
tempest = transD_GP.TEMPEST1DInversion.Bfield(
  				  ntimesperdecade = ntimesperdecade,
            nfreqsperdecade = nfreqsperdecade,
					  nkᵣeval = nkᵣeval,
            times  = times,
            ramp   = ramp,
					  zTx    = zTx,
					  zRx    = zRx,
					  x_rx   = x_rx,
					  y_rx   = y_rx,
					  rx_roll = rx_roll,
					  rx_pitch = rx_pitch,
					  rx_yaw = rx_yaw,
					  tx_roll = tx_roll,
					  tx_pitch = tx_pitch,
					  tx_yaw = tx_yaw,
					  addprimary = addprimary)
# X and Z components in fT from GA-AEM
ross_primary = [25.7847     -18.7313]
ross_secondary=[
   6.762     -8.0374
  4.7924     -6.6042
  4.0079     -5.9113
  3.4924     -5.3952
  3.1035     -4.9777
  2.7646     -4.608
  2.3721     -4.1681
  1.9079     -3.6087
  1.4085     -2.9386
0.92505      -2.1946
0.52168      -1.4595
0.24675      -0.84605
0.097237     -0.42402
0.032526     -0.18624
0.008896     -0.069836]

if addprimary
	ross_secondary .+= ross_primary
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
module PeriodicBCTests

using Gridap.Arrays
using Gridap.Geometry

using Test

## 2D mesh with periodic BC in x
model = CartesianDiscreteModel((0,1,0,1),(3,3);isperiodic=(true,false))

# Test cells to vertices
@test model.grid_topology.n_m_to_nface_to_mfaces[3,1] ==
  [[1, 2, 4, 5],
   [2, 3, 5, 6],
   [3, 1, 6, 4],
   [4, 5, 7, 8],
   [5, 6, 8, 9],
   [6, 4, 9, 7],
   [7, 8, 10, 11],
   [8, 9, 11, 12],
   [9, 7, 12, 10]]

# Test cells to faces
@test model.grid_topology.n_m_to_nface_to_mfaces[3,2] ==
  [[1, 2, 3, 4],
   [5, 6, 4, 7],
   [8, 9, 7, 3],
   [2, 10, 11, 12],
   [6, 13, 12, 14],
   [9, 15, 14, 11],
   [10, 16, 17, 18],
   [13, 19, 18, 20],
   [15, 21, 20, 17]]

# Test faces to vertices
@test model.grid_topology.n_m_to_nface_to_mfaces[2,1] ==
  [[1, 2], [4, 5], [1, 4], [2, 5], [2, 3], [5, 6], [3, 6], [3, 1], [6, 4],
   [7, 8], [4, 7], [5, 8], [8, 9], [6, 9], [9, 7],
   [10, 11], [7, 10], [8, 11], [11, 12], [9, 12], [12, 10]]

# Test vertices to faces
@test model.grid_topology.n_m_to_nface_to_mfaces[1,2] ==
  [[1, 3, 8],
   [1, 4, 5],
   [5, 7, 8],
   [2, 3, 9, 11],
   [2, 4, 6, 12],
   [6, 7, 9, 14],
   [10, 11, 15, 17],
   [10, 12, 13, 18],
   [13, 14, 15, 20],
   [16, 17, 21],
   [16, 18, 19],
   [19, 20, 21]]

# Test vertices to cells
@test model.grid_topology.n_m_to_nface_to_mfaces[1,3] ==
  [[1, 3],
   [1, 2],
   [2, 3],
   [1, 3, 4, 6],
   [1, 2, 4, 5],
   [2, 3, 5, 6],
   [4, 6, 7, 9],
   [4, 5, 7, 8],
   [5, 6, 8, 9],
   [7, 9],
   [7, 8],
   [8, 9]]

# Test faces to cells
@test model.grid_topology.n_m_to_nface_to_mfaces[2,3] ==
  [[1], [1, 4], [1, 3], [1, 2], [2], [2, 5], [2, 3], [3], [3, 6],
   [4, 7], [4, 6], [4, 5], [5, 8], [5, 6], [6, 9],
   [7], [7, 9], [7, 8], [8], [8, 9], [9]]

## 2D mesh with periodic BC in x & y
model = CartesianDiscreteModel((0,1,0,1),(3,3);isperiodic=(true,true))

 # Test cells to vertices
@test model.grid_topology.n_m_to_nface_to_mfaces[3,1] ==
  [[1, 2, 4, 5],
   [2, 3, 5, 6],
   [3, 1, 6, 4],
   [4, 5, 7, 8],
   [5, 6, 8, 9],
   [6, 4, 9, 7],
   [7, 8, 1, 2],
   [8, 9, 2, 3],
   [9, 7, 3, 1]]

# Test cells to faces
@test model.grid_topology.n_m_to_nface_to_mfaces[3,2] ==
  [[1, 2, 3, 4],
   [5, 6, 4, 7],
   [8, 9, 7, 3],
   [2, 10, 11, 12],
   [6, 13, 12, 14],
   [9, 15, 14, 11],
   [10, 1, 16, 17],
   [13, 5, 17, 18],
   [15, 8, 18, 16]]

# Test faces to vertices
@test model.grid_topology.n_m_to_nface_to_mfaces[2,1] ==
  [[1, 2], [4, 5], [1, 4], [2, 5], [2, 3], [5, 6], [3, 6], [3, 1], [6, 4],
   [7, 8], [4, 7], [5, 8], [8, 9], [6, 9], [9, 7],
   [7, 1], [8, 2], [9, 3]]

# Test vertices to faces
@test model.grid_topology.n_m_to_nface_to_mfaces[1,2] ==
  [[1, 3, 8, 16],
   [1, 4, 5, 17],
   [5, 7, 8, 18],
   [2, 3, 9, 11],
   [2, 4, 6, 12],
   [6, 7, 9, 14],
   [10, 11, 15, 16],
   [10, 12, 13, 17],
   [13, 14, 15, 18]]

# Test vertices to cells
@test model.grid_topology.n_m_to_nface_to_mfaces[1,3] ==
  [[1, 3, 7, 9],
   [1, 2, 7, 8],
   [2, 3, 8, 9],
   [1, 3, 4, 6],
   [1, 2, 4, 5],
   [2, 3, 5, 6],
   [4, 6, 7, 9],
   [4, 5, 7, 8],
   [5, 6, 8, 9]]

# Test faces to cells
@test model.grid_topology.n_m_to_nface_to_mfaces[2,3] ==
  [[1, 7], [1, 4], [1, 3], [1, 2], [2, 8], [2, 5], [2, 3], [3, 9], [3, 6],
   [4, 7], [4, 6], [4, 5], [5, 8], [5, 6], [6, 9],
   [7, 9], [7, 8], [8, 9]]

end #module


text:
## RHEOLOGY ###

@kwdef struct Rheology{T, D, V, E, P}
    damage::D = nothing
    viscosity::V = nothing
    elasticity::E = nothing
    plasticity::P = nothing
    function Rheology(damage, viscosity, elasticity, plasticity)
        fields = (damage, viscosity, elasticity, plasticity)
        cond1 = .!(<:).(typeof.(fields),AbstractBehavior{Function})
        cond2 = (<:).(typeof.(fields),AbstractBehavior)
        idxs = findall(cond1)
        if length(idxs) == 4
            if any(cond1 .& cond2)
                no_func_field = fields[findfirst(cond1 .& cond2)]
                T = typeof(getproperty(no_func_field,propertynames(no_func_field)[1]))
            else
                @error "Rheology instance can not be empty"
            end
        else
            T = Function
        end
        new{T,typeof(damage),typeof(viscosity),typeof(elasticity),typeof(plasticity)}(damage, viscosity, elasticity, plasticity)
    end
end
Rheology{T,D,V,E,P}(damage, viscosity, elasticity, plasticity) where {T,D,V,E,P} =
        Rheology(damage, viscosity, elasticity, plasticity)
Rheology{D,V,E,P}(damage, viscosity, elasticity, plasticity) where {D,V,E,P} =
        Rheology{Float64,D,V,E,P}(damage, viscosity, elasticity, plasticity)
Rheology{T}(damage, viscosity, elasticity, plasticity) where {T} =
        Rheology(damage, viscosity, elasticity, plasticity)
#Rheology(; damage, viscosity, elasticity, plasticity) = Rheology(damage, viscosity, elasticity, plasticity)

" Evaluate the functionaly defined rheology at coordinates x"
function (r::Rheology{Function})(x)
    behaviors = []
    for prop in propertynames(r)
        behavior = getproperty(r,prop)
        properties = []
        if typeof(behavior) <: AbstractBehavior{Function}
            for b_prop in propertynames(behavior)
                bp = getproperty(behavior,b_prop)
                if isa(bp, Function)
                    try
                    push!(properties,bp(x))
                    catch err
                        @error "An error occured during the evaluation of a field function. Make sure that the length of the coordinate argument of your defined functions matchs the dimensions of your grid"
                        throw(err)
                    end
                elseif b_prop == :Dᵉ
                    push!(properties,get_elastic_stiffness_tensor(properties[3],properties[5]))
                else
                    push!(properties,bp)
                end
            end
            push!(behaviors,updatetypeparameter(behavior,properties))
        else
            push!(behaviors,behavior)
        end
    end
    return Rheology(behaviors...)
end

(r::Rheology{T})(x) where {T<:Real} = @info "Rheology instance is already filled with Real values"

function Base.show(io::IO, ::MIME"text/plain", rheology::Rheology{T,D,V,E,P}) where {T,D,V,E,P}
    if T <: Real
        println(io, "⌗  $(typeof(rheology)) instance ")
    elseif T == Function
        println(io, "⌗  $(typeof(rheology)) instance with functional parameters ")
    end
    print(io, "    -> damage  : ")
    show(io,MIME"text/plain"(),rheology.damage); print(io,"\n")
    print(io, "    -> viscosity  : ")
    show(io,MIME"text/plain"(),rheology.viscosity); print(io,"\n")
    print(io, "    -> elasticity : ")
    show(io,MIME"text/plain"(),rheology.elasticity);
    print(io, "    -> plasticity : ")
    show(io,MIME"text/plain"(),rheology.plasticity);
end

rheology_summary(r::Rheology{T,Nothing,Nothing,E,Nothing}) where{T,E<:Elasticity} = "elastic"
rheology_summary(r::Rheology{T,Nothing,V,E,Nothing}) where{T,V<:Viscosity,E<:Elasticity} = "visco-elastic"
rheology_summary(r::Rheology{T,Nothing,V,Nothing,Nothing}) where{T,V<:Viscosity} = "viscous"
rheology_summary(r::Rheology{T,Nothing,Nothing,E,P}) where{T,E<:Elasticity,P<:Plasticity} = "elasto-plastic"
rheology_summary(r::Rheology{T,Nothing,V,E,P}) where{T,V<:Viscosity,E<:Elasticity,P<:Plasticity} = "visco-elasto-plastic"
rheology_summary(r::Rheology{T,D,Nothing,E,Nothing}) where{T,D<:Damage,E<:Elasticity} = "damaged-elastic"
rheology_summary(r::Rheology{T,D,Nothing,E,P}) where{T,D<:Damage,E<:Elasticity,P<:Plasticity} = "damaged-elasto-plastic"


text:
function contact_point_box(ps, xc, qc, kx, ky, kz)
    # corner kinematics
    cx = contact_point_origin(szeros(3), qc, kx)
    cy = contact_point_origin(szeros(3), qc, ky)
    cz = contact_point_origin(szeros(3), qc, kz)

    v1 = 2.0 * cx
    v2 = 2.0 * cy
    v3 = 2.0 * cz

    # origin is opposite corner
    origin = xc - cx - cy - cz

    @show tx = dot(ps - origin, v1) / dot(v1, v1)
    @show ty = dot(ps - origin, v2) / dot(v2, v2)
    @show tz = dot(ps - origin, v3) / dot(v3, v3)

    coc = origin

    if tx >= 1.0
        coc += v1
    elseif tx < 0.0
        nothing
    else
        coc += tx * v1
    end

    if ty >= 1.0
        coc += v2
    elseif ty < 0.0
        nothing
    else
        coc += ty * v2
    end

    if tz >= 1.0
        coc += v3
    elseif tz < 0.0
        nothing
    else
        coc += tz * v3
    end

    return coc
end

function ∂contact_point_box∂p(ps, xc, qc, kx, ky, kz)
    cx = contact_point_origin(szeros(3), qc, kx)
    cy = contact_point_origin(szeros(3), qc, ky)
    cz = contact_point_origin(szeros(3), qc, kz)

    v1 = 2.0 * cx
    v2 = 2.0 * cy
    v3 = 2.0 * cz

    origin = xc - cx - cy - cz

    tx = dot(ps - origin, v1) / dot(v1, v1)
    ty = dot(ps - origin, v2) / dot(v2, v2)
    tz = dot(ps - origin, v3) / dot(v3, v3)

    P = szeros(3, 3)

    if tx >= 1.0
        # coc += v1
        # X +=
    elseif tx < 0.0
        # nothing
        # X +=
    else
        # coc += tx * v1
        P += v1 * 1.0 / dot(v1, v1) * v1'
    end

    if ty >= 1.0
        # coc += v2
    elseif ty < 0.0
        # nothing
    else
        # coc += ty * v2
        P += v2 * 1.0 / dot(v2, v2) * v2'
    end

    if tz >= 1.0
        # coc += v3
    elseif tz < 0.0
        # nothing
    else
        # coc += tz * v3
        P += v3 * 1.0 / dot(v3, v3) * v3'
    end

    return P
    # FiniteDiff.finite_difference_jacobian(p -> contact_point_box(p, xc, qc, kx, ky, kz), ps)
end

function ∂contact_point_box∂x(ps, xc, qc, kx, ky, kz)
    cx = contact_point_origin(szeros(3), qc, kx)
    cy = contact_point_origin(szeros(3), qc, ky)
    cz = contact_point_origin(szeros(3), qc, kz)

    v1 = 2.0 * cx
    v2 = 2.0 * cy
    v3 = 2.0 * cz

    origin = xc - cx - cy - cz

    tx = dot(ps - origin, v1) / dot(v1, v1)
    ty = dot(ps - origin, v2) / dot(v2, v2)
    tz = dot(ps - origin, v3) / dot(v3, v3)

    X = 1.0 * I(3)

    if tx >= 1.0
        # coc += v1
    elseif tx < 0.0
        # nothing
    else
        # coc += tx * v1
        X += v1 * 1.0 / dot(v1, v1) * v1' * -1.0
    end

    if ty >= 1.0
        # coc += v2
    elseif ty < 0.0
        # nothing
    else
        # coc += ty * v2
        X += v2 * 1.0 / dot(v2, v2) * v2' * -1.0
    end

    if tz > 1.0
        # coc += v3
    elseif tz < 0.0
        # nothing
    else
        # coc += tz * v3
        X += v3 * 1.0 / dot(v3, v3) * v3' * -1.0
    end

    return X
    # FiniteDiff.finite_difference_jacobian(x -> contact_point_box(ps, x, qc, kx, ky, kz), xc)
end

function ∂contact_point_box∂q(ps, xc, qc, kx, ky, kz)
    cx = contact_point_origin(szeros(3), qc, kx)
    cy = contact_point_origin(szeros(3), qc, ky)
    cz = contact_point_origin(szeros(3), qc, kz)

    v1 = 2.0 * cx
    v2 = 2.0 * cy
    v3 = 2.0 * cz

    ∂v1∂q = 2.0 * ∂contact_point_origin∂q(szeros(3), qc, kx)
    ∂v2∂q = 2.0 * ∂contact_point_origin∂q(szeros(3), qc, ky)
    ∂v3∂q = 2.0 * ∂contact_point_origin∂q(szeros(3), qc, kz)

    origin = xc - cx - cy - cz

    tx = dot(ps - origin, v1) / dot(v1, v1)
    ty = dot(ps - origin, v2) / dot(v2, v2)
    tz = dot(ps - origin, v3) / dot(v3, v3)

    Q = -0.5 * (∂v1∂q + ∂v2∂q + ∂v3∂q)

    if tx >= 1.0
        # coc += v1
        Q += ∂v1∂q
    elseif tx < 0.0
        # nothing
    else
        # coc += tx * v1
        Q += tx * ∂v1∂q
        Q += v1 * (ps - xc)' ./ dot(v1, v1) * ∂v1∂q
        Q += v1 * -1.0 * dot(ps - xc, v1) ./ dot(v1, v1)^2 * 2.0 * v1' * ∂v1∂q
        Q += v1 * 1.0 / dot(v1, v1) * v1' * -0.5 * ∂v1∂q
    end

    if ty >= 1.0
        # coc += v2
        Q += ∂v2∂q
    elseif ty < 0.0
        # nothing
    else
        # coc += ty * v2
        Q += tx * ∂v2∂q
        Q += v2 * (ps - xc)' ./ dot(v2, v2) * ∂v2∂q
        Q += v2 * -1.0 * dot(ps - xc, v2) ./ dot(v2, v2)^2 * 2.0 * v2' * ∂v2∂q
        Q += v2 * 1.0 / dot(v2, v2) * v2' * -0.5 * ∂v2∂q
    end

    if tz >= 1.0
        # coc += v3
        Q += ∂v3∂q
    elseif tz < 0.0
        # nothing
    else
        # coc += tz * v3
        Q += tx * ∂v3∂q
        Q += v3 * (ps - xc)' ./ dot(v3, v3) * ∂v3∂q
        Q += v3 * -1.0 * dot(ps - xc, v3) ./ dot(v3, v3)^2 * 2.0 * v3' * ∂v3∂q
        Q += v3 * 1.0 / dot(v3, v3) * v3' * -0.5 * ∂v3∂q
    end

    return Q
    # FiniteDiff.finite_difference_jacobian(q -> contact_point_box(ps, xc, Quaternion(q..., false), kx, ky, kz), vector(qc))
end

# x_min = -1.0
# x_max = 1.0
# z_min = -1.0
# z_max = 1.0

# p = [-1.0; 1.0] #.- 1.0e-6

# function klamp(t, a, b)
#     if t <= a
#         return a
#     elseif b <= t
#         return b
#     else
#         da = abs(t - a)
#         db = abs(t - b)

#         if da < db
#             return a
#         else
#             return b
#         end
#     end
# end

# function nearest_point(p, x_min, x_max, z_min, z_max)
#     return [klamp(p[1], x_min, x_max), klamp(p[2], z_min, z_max)]
# end

# nearest_point(p, x_min, x_max, z_min, z_max)


text:
module FilterCircuit

using TinyModia

using DifferentialEquations
using ModiaPlot

setLogMerge(false)

include("../models/Electric.jl")

Filter = Model(
    R = Resistor | Map(R=0.5u"Ω"),
    C = Capacitor | Map(C=2.0u"F", v=Var(init=0.1u"V")),
    V = ConstantVoltage | Map(V=10.0u"V"),
    ground = Ground,
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

model = @instantiateModel(Filter, log=false, aliasReduction=true, logCode=false)
@time simulate!(model, Tsit5(), stopTime = 10)
#plot(model, [("R.v", "C.v")])

println("Simulate once more with different R.R")
@time simulate!(model, Tsit5(), stopTime = 10, merge = Map(R = Map(R = 5u"Ω")), requiredFinalStates = [6.3579935215716095])
plot(model, [("R.v", "C.v")])

# setLogMerge(true)
println("Filter without ground and parameter propagation")
Filter2 = Model(
    r = 2.0u"Ω",
    c = 1.0u"F",
    v = 10u"V",
    R = Resistor | Map(R=:r),
    C = Capacitor | Map(C=:c),
    V = ConstantVoltage | Map(V=:v),
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n)
    ]
)

# @showModel(Filter2)

model = @instantiateModel(Filter2, log=false, aliasReduction=true, logCode=false)
simulate!(model, Tsit5(), stopTime = 10, requiredFinalStates = [9.932620374719848])
plot(model, [("R.v", "C.v")])

println("Voltage divider by redeclaring capacitor to resistor")
Cpar = Map(C = 5.0u"F")

TwoFilters = Model( f1 = Filter | Map( R = Map(R = 10.0u"Ω"), C = Cpar), f2 = Filter)

VoltageDividerAndFilter = TwoFilters | Map(f1 = Map(C = Redeclare | Resistor | (R = 20.0u"Ω", v = Var(start = 0u"V"))))

model = @instantiateModel(VoltageDividerAndFilter, log=false, aliasReduction=true, logCode=false)
simulate!(model, Tsit5(), stopTime = 10, requiredFinalStates = [9.999550454584188])
plot(model, [("f1.R.v", "f1.C.v"), ("f2.R.v", "f2.C.v")])

println("Build array of filters")
Filters = Model(
    filters = [Filter | Map( R = Map(R = (10.0+5*i)*u"Ω")) for i in 1:10]
)

setLogMerge(false)

model = @instantiateModel(Filters, log=false, aliasReduction=true, logCode=false)

simulate!(model, Tsit5(), stopTime = 10, requiredFinalStates =
[2.9063400246452358, 2.2898722474751163, 1.8945655444974656, 1.6198309235728083, 1.4179087924692246, 1.2632806644107324, 1.1410907635368692, 1.0421095614435398, 0.9603029088053439, 0.8915602951695468])
plot(model, [("filters_1.R.v", "filters_1.C.v"), ("filters_2.R.v", "filters_2.C.v")])

end


text:
using ConstrainedDynamics

path = "examples/examples_files/doublependulum.urdf"
mech = Mechanism(path)

simulate!(mech, save = true)
visualize!(mech)


text:
# ============================================================
# Phase Space Methods
# ============================================================

export VSpace1D, VSpace2D, VSpace3D, UnstructVSpace
export MVSpace1D, MVSpace2D, MVSpace3D
export newton_cotes
export legendre_quadrature, octa_quadrature, quadrature_weights

include("creamer.jl")
include("velocity.jl")
include("quadrature.jl")


text:
function approximate(f, X)
    method = KernelInterpolation(f, X)
    train!(method)
    return x -> evaluate(x, method)
end

function approximate(method)
    if length(method.coeffs) == 0
        train!(method)
    end

    return x -> evaluate(x, method)
end


text:
""
function objective_min_fuel_cost(pm::_PM.AbstractPowerModel)
    model = _PM.check_cost_models(pm)
    if model == 1
        return objective_min_pwl_fuel_cost(pm)
    elseif model == 2
        return objective_min_polynomial_fuel_cost(pm)
    else
        error("Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model)")
    end

end

""
function objective_min_polynomial_fuel_cost(pm::_PM.AbstractPowerModel)
    order = _PM.calc_max_cost_index(pm.data)-1

    if order == 1
        return _objective_min_polynomial_fuel_cost_linear(pm)
    elseif order == 2
        return _objective_min_polynomial_fuel_cost_quadratic(pm)
    else
        error("cost model order of $(order) is not supported")
    end
end

function _objective_min_polynomial_fuel_cost_linear(pm::_PM.AbstractPowerModel)
    from_idx = Dict()
    for (n, nw_ref) in nws(pm)
        from_idx[n] = Dict(arc[1] => arc for arc in nw_ref[:arcs_from_dc])
    end

    return JuMP.@objective(pm.model, Min,
        sum(
            sum(   gen["cost"][1]*sum( var(pm, n, c, :pg, i) for c in conductor__PM.ids(pm, n))+
                   gen["cost"][2] for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in nws(pm))
    )
end
""
function _objective_min_polynomial_fuel_cost_quadratic(pm::_PM.AbstractPowerModel)
    from_idx = Dict()
    for (n, nw_ref) in nws(pm)
        from_idx[n] = Dict(arc[1] => arc for arc in nw_ref[:arcs_from_dc])
    end

    return JuMP.@objective(pm.model, Min,
        sum(
            sum(   gen["cost"][1]*sum( var(pm, n, c, :pg, i) for c in conductor__PM.ids(pm, n))^2 +
                   gen["cost"][2]*sum( var(pm, n, c, :pg, i) for c in conductor__PM.ids(pm, n))+
                   gen["cost"][3] for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in nws(pm))
    )
end
""
function objective_min_polynomial_fuel_cost(pm::_PM.AbstractConicModel)
    _PM.check_polynomial_cost_models(pm)
    pg_sqr = Dict()
    for (n, nw_ref) in _PM.nws(pm)
        for cnd in _PM.conductor__PM.ids(pm, n)
            pg_sqr = _PM.var(pm, n, cnd)[:pg_sqr] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, n, :gen)], base_name="$(n)_$(cnd)_pg_sqr",
                lower_bound = _PM.ref(pm, n, :gen, i, "pmin", cnd)^2,
                upper_bound = _PM.ref(pm, n, :gen, i, "pmax", cnd)^2
            )
            for (i, gen) in nw_ref[:gen]
                JuMP.@constraint(pm.model, [pg_sqr[i], var(pm, n, cnd, :pg, i)/sqrt(2), var(pm, n, cnd, :pg, i)/sqrt(2)] in JuMP.SecondOrderCone())
            end

        end
    end
    return JuMP.@objective(pm.model, Min,
        sum(
            sum( gen["cost"][1]*sum(_PM.var(pm, n, cnd, :pg_sqr, i) for cnd in _PM.conductor__PM.ids(pm, n)) +
                 gen["cost"][2]*sum(_PM.var(pm, n, cnd, :pg, i) for cnd in _PM.conductor__PM.ids(pm, n)) +
                 gen["cost"][3] for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in _PM.nws(pm))
    )
end
""
function objective_min_pwl_fuel_cost(pm::_PM.AbstractPowerModel)

    for (n, nw_ref) in _PM.nws(pm)
        pg_cost = _PM.var(pm, n)[:pg_cost] = JuMP.@variable(pm.model,
            [i in _PM.ids(pm, n, :gen)], base_name="$(n)_pg_cost"
        )

        # pwl cost
        gen_lines = _PM.get_lines(nw_ref[:gen])
        for (i, gen) in nw_ref[:gen]
            for line in gen_lines[i]
                JuMP.@constraint(pm.model, pg_cost[i] >= line["slope"]*sum(_PM.var(pm, n, cnd, :pg, i) for cnd in _PM.conductor__PM.ids(pm, n)) + line["intercept"])
            end
        end

    end

    return JuMP.@objective(pm.model, Min,
        sum(
            sum( _PM.var(pm, n,:pg_cost, i) for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in _PM.nws(pm))
        )
end

##################### TNEP Objective   ###################
function objective_min_cost(pm::_PM.AbstractPowerModel)
        gen_cost = Dict()
        for (n, nw_ref) in _PM.nws(pm)
            for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        return JuMP.@objective(pm.model, Min,
            sum(
                sum(conv["cost"]*_PM.var(pm, n, :conv_ne, i) for (i,conv) in nw_ref[:convdc_ne])
                +
                sum(branch["cost"]*_PM.var(pm, n, :branchdc_ne, i) for (i,branch) in nw_ref[:branchdc_ne])
                +
                sum( gen_cost[(n,i)] for (i,gen) in nw_ref[:gen] )
                for (n, nw_ref) in _PM.nws(pm)
                    )
        )
end

function objective_min_cost_acdc(pm::_PM.AbstractPowerModel)
        gen_cost = Dict()
        for (n, nw_ref) in _PM.nws(pm)
            for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        return JuMP.@objective(pm.model, Min,
            sum(
                sum(conv["cost"]*_PM.var(pm, n, :conv_ne, i) for (i,conv) in nw_ref[:convdc_ne])
                +
                sum(branch["construction_cost"]*_PM.var(pm, n, :branch_ne, i) for (i,branch) in nw_ref[:ne_branch])
                +
                sum(branch["cost"]*_PM.var(pm, n, :branchdc_ne, i) for (i,branch) in nw_ref[:branchdc_ne])
                +
                sum( gen_cost[(n,i)] for (i,gen) in nw_ref[:gen] )
                for (n, nw_ref) in _PM.nws(pm)
                    )
        )
end


text:
# using Revise
using Test, BifurcationKit, LinearAlgebra, Setfield, ForwardDiff
const BK = BifurcationKit

norminf = (x)->norm(x, Inf)
source_term(x; a = 0.5, b = 0.01) = 1 + (x + a*x^2)/(1 + b*x^2)
limits(x,i,N, b) = (i<1||i>N) ? b : x[i]

function F_chan(x, p)
	α, β = p
	N = length(x)
	f = similar(x)
	N = length(x)
	ind(x,i) = limits(x,i,N,β)
	for i=1:N
		f[i] = (ind(x,i-1) - 2 * x[i] + ind(x,i+1)) * (N-1)^2 + α * source_term(x[i], b = β)
	end
	return f
end

J_chan(u0, p) = ForwardDiff.jacobian(u -> F_chan(u,p), u0)
jet = BK.getJet(F_chan, J_chan)
par_chan = (α = 3.3, β = 0.01)

n = 101
	sol = [(i-1)*(n-i)/n^2+0.1 for i=1:n]
	opt_newton = NewtonPar(tol = 1e-9, verbose = false)
	out, hist, flag = newton(jet[1], jet[2], sol, par_chan,
							opt_newton, normN = norminf)

# _J = J_chan(sol, par_chan)
# heatmap(_J - _J', yflip=true)

# test with secant continuation
opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, pMax = 4.1, maxSteps = 150, newtonOptions = opt_newton, detectBifurcation = 3)
	br, = continuation(jet[1], jet[2],
		out, par_chan, (@lens _.α), opts_br0; )
####################################################################################################
# deflation newton solver, test of jacobian expression
deflationOp = DeflationOperator(2.0, dot, 1.0, [out])

# for testing
show(deflationOp)
length(deflationOp)
deflationOp(2out, out)
push!(deflationOp, rand(n))
deleteat!(deflationOp, 2)

chanDefPb = DeflatedProblem(jet[1], jet[2], deflationOp)

opt_def = setproperties(opt_newton; tol = 1e-10, maxIter = 1000)
outdef1, = newton(chanDefPb, out .* (1 .+0.01*rand(n)), par_chan, opt_def)

# we now compare the jacobians for the deflated problem either using finite differences or the explicit jacobian
rhs = rand(n)
J_def_fd = BK.finiteDifferences(u->chanDefPb(u, par_chan),1.5*out)
res_fd =  J_def_fd \ rhs

Jacdf = (u0, pb::DeflatedProblem, ls = opt_def.linsolve ) -> (return (u0, par_chan, pb, ls))
Jacdfsolver = DeflatedLinearSolver()

res_explicit = Jacdfsolver(Jacdf(1.5out, chanDefPb, opt_def.linsolver),rhs)[1]

# Test jacobian expression for deflated problem
@test norm(res_fd - res_explicit,Inf64) < 1e-4

opt_def = setproperties(opt_newton; tol = 1e-10, maxIter = 1000)
outdef1, = newton(
		jet[1], jet[2],
		out.*(1 .+ 0.1*rand(n)), par_chan,
		opt_def, deflationOp)

####################################################################################################
# Fold continuation, test of Jacobian expression
outfold = newtonFold(jet[1], jet[2], br, 2; startWithEigen = true)
@test  outfold[3] && outfold[4] == 2
outfold = newtonFold(jet[1], jet[2], br, 2; startWithEigen = true, issymmetric = true)
@test  outfold[3] && outfold[4] == 2
outfold = newton(jet[1], jet[2], br, 2; startWithEigen = true, issymmetric = true)
@test  outfold[3] && outfold[4] == 2

optcontfold = ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds= 0.01, pMax = 4.1, pMin = 0., newtonOptions = NewtonPar(verbose=false, tol = 1e-8), maxSteps = 50, detectBifurcation = 2)
outfoldco, = continuationFold(jet[1], jet[2], br, 2, (@lens _.β), optcontfold; startWithEigen = true, updateMinAugEveryStep = 1, plot = false)
outfoldco, = continuationFold(jet[1], jet[2], br, 2, (@lens _.β), optcontfold; startWithEigen = true, updateMinAugEveryStep = 1, issymmetric = true, plot = false)
outfoldco, = continuationFold(jet[1], jet[2], br, 2, (@lens _.β), optcontfold; startWithEigen = true, updateMinAugEveryStep = 1, issymmetric = true, d2F = jet[3], plot = false)

# manual handling
indfold = 1
foldpt = FoldPoint(br, indfold)
foldpb = FoldProblemMinimallyAugmented(
		jet[1], jet[2],
		nothing, nothing,
		(@lens _.α),
		br.specialpoint[indfold].x,
		br.specialpoint[indfold].x,
		opts_br0.newtonOptions.linsolver)
foldpb(foldpt, par_chan) |> norm

outfold, = newtonFold(jet[1], jet[2], foldpt, par_chan, (@lens _.α), br.specialpoint[indfold].x, br.specialpoint[indfold].x, NewtonPar(verbose=false))
	# println("--> Fold found at α = ", outfold.p, " from ", br.specialpoint[indfold].param)

# user defined Fold Problem
indfold = 1

# we define the following wrappers to be able to use ForwardDiff
Bd2Vec(x) = vcat(x.u, x.p)
Vec2Bd(x) = BorderedArray(x[1:end-1], x[end])
foldpbVec(x,p) = Bd2Vec(foldpb(Vec2Bd(x),p))

rhs = rand(n+1)
Jac_fold_fdMA(u0) = ForwardDiff.jacobian( u -> foldpbVec(u, par_chan), u0)
J_fold_fd = Jac_fold_fdMA(Bd2Vec(foldpt))
res_fd =  J_fold_fd \ rhs

Jac_fold_MA(u0, p, pb::FoldProblemMinimallyAugmented) = (return (x=u0, params=p, fldpb = pb))
jacFoldSolver = BK.FoldLinearSolverMinAug()
debugTmpForσ = zeros(n+1,n+1) # temporary array for debugging σ
res_explicit = jacFoldSolver(Jac_fold_MA(foldpt, par_chan, foldpb), Vec2Bd(rhs); debugArray = debugTmpForσ)

Jac_fold_MA(foldpt, 0.01, foldpb)[2]

# test whether the Jacobian Matrix for the Fold problem is correct
@test debugTmpForσ[1:end-1,1:end-1] - J_fold_fd[1:end-1,1:end-1] |> norminf == 0

# we test the expression for σp
σp_fd = J_fold_fd[end,end]
σp_fd_ana = debugTmpForσ[end,end]
@test σp_fd ≈ σp_fd_ana rtol = 1e-5

# we test the expression for σx
σx_fd = J_fold_fd[end,1:end-1]
σx_ana = debugTmpForσ[end,1:end-1]
@test σx_fd ≈ σx_ana rtol = 1e-2

σx_fd - σx_ana |> norminf
debugTmpForσ[1:end-1,1:end-1] - J_fold_fd[1:end-1,1:end-1] |> norminf

# check our solution of the bordered problem
res_exp = debugTmpForσ \ rhs
@test norm(res_exp - Bd2Vec(res_explicit[1]), Inf64) < 1e-10
#############################################
# same with foldpb.issymmetric = true
@set! foldpb.issymmetric = true
rhs = rand(n+1)
Jac_fold_fdMA(u0) = ForwardDiff.jacobian( u -> foldpbVec(u, par_chan), u0)
J_fold_fd = Jac_fold_fdMA(Bd2Vec(foldpt))
res_fd =  J_fold_fd \ rhs

Jac_fold_MA(u0, p, pb::FoldProblemMinimallyAugmented) = (return (x=u0, params=p, fldpb = pb))
jacFoldSolver = BK.FoldLinearSolverMinAug()
debugTmpForσ = zeros(n+1,n+1) # temporary array for debugging σ
res_explicit = jacFoldSolver(Jac_fold_MA(foldpt, par_chan, foldpb), Vec2Bd(rhs); debugArray = debugTmpForσ)

Jac_fold_MA(foldpt, 0.01, foldpb)[2]

# test whether the Jacobian Matrix for the Fold problem is correct
@test debugTmpForσ[1:end-1,1:end-1] - J_fold_fd[1:end-1,1:end-1] |> norminf == 0

# we test the expression for σp
σp_fd = J_fold_fd[end,end]
σp_fd_ana = debugTmpForσ[end,end]
@test σp_fd ≈ σp_fd_ana rtol = 1e-5

# we test the expression for σx
σx_fd = J_fold_fd[end,1:end-1]
σx_ana = debugTmpForσ[end,1:end-1]
@test σx_fd ≈ σx_ana rtol = 1e-2

σx_fd - σx_ana |> norminf
debugTmpForσ[1:end-1,1:end-1] - J_fold_fd[1:end-1,1:end-1] |> norminf

# check our solution of the bordered problem
res_exp = debugTmpForσ \ rhs
@test norm(res_exp - Bd2Vec(res_explicit[1]), Inf64) < 1e-10
####################################################################################################
# Use of different eigensolvers
opt_newton = NewtonPar(tol = 1e-8, verbose = false, eigsolver = EigKrylovKit())
opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, pMax = 4.1, maxSteps = 250, newtonOptions = opt_newton, detectFold = true, detectBifurcation = 1, nev = 15)

br, = continuation(jet[1], jet[2], out, par_chan, (@lens _.α), opts_br0, recordFromSolution = (x,p)->norm(x,Inf64), plot = false, verbosity = 0)

opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, pMax = 4.1, maxSteps = 250, newtonOptions = NewtonPar(tol =1e-8), detectFold = true, detectBifurcation = 1, nev = 15)

br, = continuation(jet[1], jet[2], out, par_chan, (@lens _.α),opts_br0,plot = false, verbosity = 0)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ff9de46a-7cb3-11eb-3d0a-41c763671f20
using Luxor, Colors, PlutoUI

# ╔═╡ 0972cbf2-7cb4-11eb-3309-a772deb54308
function color_circles(x)
	b = blend(Point(-200, 0), Point(200, 0), "orange", "blue")
	# hues = colormap("Oranges", 11)
	hues = range(colorant"purple", colorant"white", length=11)
	for (i, colr) in zip(110:-10:0, hues)
		# randomhue()
		# sethue(1 - 1 * i / 110, 0.4 - 1 * i / 110, 1 - 0.4 * i / 110)
		sethue(colr)
		circle(i*x, 0, i, :fill)
	end
end

# ╔═╡ 2b5cedb2-7cb9-11eb-33a2-4f2e02ee742a
@bind x Slider(-1.0:0.05:1.0)

# ╔═╡ 4123ccf6-7cb9-11eb-0a81-87b0cc299d39
@svg color_circles(x) 500 250 "color_circles_slide"

# ╔═╡ 0e403470-7cb7-11eb-24b2-4bee851633b4
function backdrop(scene, framenumber)
    background("white")
end

# ╔═╡ a642bac6-7cb8-11eb-0f12-450f251ea610
function frame(scene, framenumber)
	n = scene.easingfunction(framenumber, -1, 4, scene.framerange.stop)
	if n > 1
		n = 2 - n
	end
	color_circles(n)
end

# ╔═╡ bb178434-7cb8-11eb-2199-472c37a336eb
begin
	demo = Movie(500, 250, "color_circles")
	animate(demo, [
	    Scene(demo, backdrop, 0:250),
	    Scene(demo, frame, 0:250,
	        easingfunction=easingflat)
	    ],
	    creategif=true,
		pathname="color_circles_easingflat_onion.gif")
end

# ╔═╡ 625687c0-7cc0-11eb-187e-6919e50e5c34
range(colorant"red", colorant"green", length=11)

# ╔═╡ Cell order:
# ╠═ff9de46a-7cb3-11eb-3d0a-41c763671f20
# ╠═0972cbf2-7cb4-11eb-3309-a772deb54308
# ╠═2b5cedb2-7cb9-11eb-33a2-4f2e02ee742a
# ╠═4123ccf6-7cb9-11eb-0a81-87b0cc299d39
# ╟─0e403470-7cb7-11eb-24b2-4bee851633b4
# ╟─a642bac6-7cb8-11eb-0f12-450f251ea610
# ╠═bb178434-7cb8-11eb-2199-472c37a336eb
# ╠═625687c0-7cc0-11eb-187e-6919e50e5c34


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
            arg_type = String
        "--process_type"
            help = "The type of process to execute."
            arg_type = String
        "--log_dir"
            help = "The directory where to store the logs of the process."
            arg_type = String
        "--config_key"
            help = "Key corresponding to the keyword arguments used in this process."
            arg_type = String
        "--config_path"
            help = "Path to a julia file containing a dict with the keyword arguments."
            arg_type = String
            default = "data/configs/kwargs.jl"
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    include("../$(parsed_args["config_path"])")
    kwargs = process_kwargs[parsed_args["config_key"]]
    kwargs = Dict(Symbol(k) => v for (k, v) in kwargs)
    execute_process_mpi(db, parsed_args["process_type"], parsed_args["log_dir"];
                        kwargs...)
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using Decimals
using Test

@testset "Decimal constructor" begin

@test isa(d, Array{Decimal,1})

end


text:
"""

$(TYPEDEF)

Auxiliary structure used to update counters.

$(TYPEDFIELDS)

"""
mutable struct DminMol
  d::Float64
  jmol::Int
  iat::Int
  jat::Int
end