text:
################################################
#
# Functions to tear systems of equations
#
# Author: Martin Otter, DLR-SR (first version: Jan. 14, 2017)
#
# Details are described in the paper:
#   Otter, Elmqvist (2017): Transformation of Differential Algebraic Array Equations to
#                           Index One Form. Modelica'2017 Conference.
#
# The following utility algorithm is used below to incrementally added edges to a
# DAG (Directed Acyclic Graph). This algorithm leads to an O(n*m) worst time complexity of the
# tearing (instead of O(m*m)) where n is the number of equations and m is the number of
# variable incidences. Note, the traversals of the DAG are not performed with recursive function
# calls but with while loops and an explicit stack, in order to avoid function stack overflow
# for large algebraic loops.
#
#   Bender, Fineman, Gilbert, Tarjan:
#      A New Approach to Incremental Cycle Detection and Related Problems.
#      ACM Transactions on Algorithms, Volume 12, Issue 2, Feb. 2016
#      http://dl.acm.org/citation.cfm?id=2756553
#
#   Text excerpt from this paper (the advantage of Algorithm N is that it
#   is simple to implement and needs no sophisticated data structures)
#
#       3. A ONE-WAY-SEARCH ALGORITHM FOR DENSE GRAPHS
#           [..]
#       To develop the algorithm, we begin with a simple algorithm and then modify it to
#       improve its running time. We call the simple algorithm Algorithm N (for “naive”). The
#       algorithm initializes k(v) = 1 for each vertex v. The initial vertex levels are a weak
#       topological numbering since the initial graph contains no arcs. To add an arc (v,w),
#       if k(v) < k(w), merely add the arc. If, on the other hand, k(v) ≥ k(w), add the arc and
#       then do a selective forward search that begins by traversing (v,w) and continues until
#       the search traverses an arc into v (there is a cycle) or there are no more candidate
#       arcs to traverse. To traverse an arc (x, y), if y = v, stop (there is a cycle); otherwise, if
#       k(x) ≥ k(y), increase k(y) to k(x)+1 and add all arcs (y, z) to the set of candidate arcs to
#       traverse.
#       It is easy to show that (1) after each arc addition that does not form a cycle
#       the vertex levels are a weak topological numbering, (2) Algorithm N is correct, and
#       (3) 1 ≤ k(v) ≤ size(v) ≤ n for all v. Since an arc (x, y) notEqual (v,w) is only traversed as a
#       result of k(x) increasing, each arc is traversed at most n times, resulting in a total time
#       bound of O(nm) for all arc additions.
#
################################################

const Undefined = typemin(Int)

"""
    td = TraverseDAG(G,nv)

Generate an object td to traverse a set of equations that are
represented as a DAG (Directed Acyclic Graph).
G is the bipartite graph of all relevant equations
and nv is the largest variable number used in G (or larger).
"""
mutable struct TraverseDAG
   minlevel::Int
   curlevel::Int
   level::Vector{Int}
   lastlevel::Vector{Int}
   levelStack::Vector{Int}
   visitedStack::Vector{Int}
   vActive::Vector{Bool}
   visited::Vector{Bool}
   check::Vector{Bool}
   stack::Vector{Int}
   eSolved::Vector{Int}
   vSolved::Vector{Int}
   G                      # Vector{ Vector{Int} }
   assign::Vector{Int}
   es::Vector{Int}
   vs::Vector{Int}

   function TraverseDAG(G,nv::Int)
      visited      = fill(false,length(G))
      check        = fill(false,length(G))
      vActive      = fill(false,nv)
      level        = fill(Undefined,nv)
      lastlevel    = fill(Undefined,nv)
      levelStack   = fill(0,0)
      stack        = fill(0,0)
      visitedStack = fill(0,0)
      eSolved      = fill(0,0)
      vSolved      = fill(0,0)
      assign       = fill(0,nv)

      new(0, Undefined, level, lastlevel, levelStack, visitedStack, vActive,
          visited, check, stack, eSolved, vSolved, G, assign)
   end
end

"""
    initAlgebraicSystem(td::TraverseDAG,es,vs)

Define the set of equations and the set variables for which the equations shall be solved for
(equations es shall be solved for variables vs) and re-initialize td.

eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1]
"""
function initAlgebraicSystem(td::TraverseDAG,es::Vector{Int},vs::Vector{Int},
                             eSolvedFixed::Vector{Int},vSolvedFixed::Vector{Int},vTearFixed::Vector{Int})
   # check arguments
   for i in eachindex(es)
      if es[i] <= 0
	     error("\n\n... Internal error in Tearing.jl: es[",i,"] = ", es[i],".\n")
      end
   end
   for i in eachindex(vs)
      if vs[i] <= 0
	     error("\n\n... Internal error in Tearing.jl: vs[",i,"] = ", vs[i],".\n")
	  end
   end

   # check that all elements of eSolvedFixed are in es, vSolvedFixed in vs,
   # vTearFixed in vs and that vSolvedFixed and vTearFixed have no variables in common
   @assert( length(eSolvedFixed) == length(vSolvedFixed) )
   ediff = setdiff(eSolvedFixed, es)
   @assert(length(ediff) == 0)
   vdiff = setdiff(vSolvedFixed, vs)
   @assert(length(vdiff) == 0)
   vdiff2 = setdiff(vTearFixed,vs)
   @assert(length(vdiff2) == 0)
   vdiff3 = intersect(vSolvedFixed,vTearFixed)
   @assert(length(vdiff3) == 0)

   # Re-initialize td
   td.minlevel = 0
   td.curlevel = Undefined
   for i in eachindex(td.visited)
      td.visited[i] = false
      td.check[i]   = false
   end
   for i in eachindex(td.vActive)
      td.vActive[i]   = false
      td.assign[i]    = 0
      td.level[i]     = Undefined
      td.lastlevel[i] = Undefined
   end
   for i in eachindex(vs)
      td.vActive[ vs[i] ] = true
   end
   empty!(td.levelStack)
   empty!(td.stack)
   empty!(td.visitedStack)
   empty!(td.eSolved)
   empty!(td.vSolved)

   # Define initial DAG
   vs2 = Int[]
   for i in eachindex(vSolvedFixed)
      vFixed = vSolvedFixed[i]
      td.assign[vFixed] = eSolvedFixed[i]
      td.level[ vFixed] = i
      push!(vs2, vFixed)
   end
   for i in eachindex(vTearFixed)
      td.vActive[ vTearFixed[i] ] = false   # vTearFixed shall not be assigned
   end

   # Store es, vs in td
   td.es = es
   td.vs = vs

   return vs2
end

in_vs(td,v) = td.vActive[v]

function setlevel(td::TraverseDAG, v::Int, parentLevel::Int)
   td.lastlevel[v] = td.level[v]
   td.level[v]     = parentLevel + 1
   push!(td.visitedStack,v)
end

"""
    success = visit!(td::TraverseDAG, v)

Traverse potential DAG starting from new variable node v.
If no cycle is detected return true, otherwise return false.
"""
function visit!(td::TraverseDAG, vcheck::Int)
   empty!(td.stack)
   empty!(td.levelStack)
   empty!(td.visitedStack)
   td.curlevel = td.level[vcheck]
   push!(td.levelStack,td.curlevel)
   push!(td.stack,vcheck)
   first = true

   while length(td.stack) > 0
      parentLevel = pop!(td.levelStack)
      veq = pop!(td.stack)
      eq  = td.assign[veq]
      if first
         first = false
      else
         if td.level[veq] == td.curlevel
            # cycle detected
            return false
         elseif td.level[veq] == Undefined || td.level[veq] <= parentLevel
            setlevel(td,veq,parentLevel)
         end
      end

      if eq > 0
         # Push all child nodes on stack
         parentLevel = td.level[veq]
         for v in td.G[eq]
            if in_vs(td,v) && v!=veq   # v is an element of td.vs and is not the variable to solve for
               eq2 = td.assign[v]
               if eq2 == 0 || td.level[v] <= parentLevel
                  push!(td.levelStack,parentLevel)
                  push!(td.stack,v)
               end
            end
         end
      end
   end

   return true
end

"""
    visit2!(td::TraverseDAG,v)

Traverse DAG starting from variable v and store visited equations and variables in stacks
eSolved, vSolved. If a cycle is deteced, raise an error (signals a programming error).
"""
function visit2!(td::TraverseDAG, vVisit::Int)
   push!(td.stack,vVisit)
   while length(td.stack) > 0
      veq = td.stack[end]
      eq  = td.assign[veq]
      if !td.visited[eq]
         td.visited[eq] = true
         td.check[eq]   = true
         for v in td.G[eq]
            if in_vs(td,v) && v!=veq  # v is an element of td.vs and is not the variable to solve for
               eq2 = td.assign[v]
               if eq2 != 0
                  if !td.visited[eq2]   # visit eq2 if not yet visited
                     push!(td.stack,v)
                  elseif td.check[eq2]  # cycle detected
                     error("... error in Tearing.jl code: \n",
                           "    cycle detected (should not occur): eq = ", eq, ", veq = ", veq, ", eq2 = ", eq2, ", v = ", v)
                  end
               end
            end
         end
      else
         td.check[eq] = false
         push!(td.eSolved,eq)
         push!(td.vSolved,veq)
         pop!(td.stack)
      end
   end
   nothing
end

"""
    (eSolved, vSolved) = sortDAG!(td::TraverseDAG, vs)

Sort the equations that are assigned by variables vs using object td of type TraverseDAG
and return the sorted equations eSolved and assigned variables vSolved.
"""
function sortDAG!(td::TraverseDAG, vs::Vector{Int})
   # initialize data structure
   empty!(td.stack)
   empty!(td.eSolved)
   empty!(td.vSolved)
   for i in eachindex(td.visited)
      td.visited[i] = false
      td.check[i]   = false
   end

   # visit all assigned variables and equations
   for veq in vs
      if !td.visited[ td.assign[veq] ]
         visit2!(td,veq)
      end
   end

   return (td.eSolved, td.vSolved)
end

"""
    (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs; eSolvedFixed=Int[], vSolvedFixed=Int[], vTearFixed=Int[])

Equations es shall be solved with respect to variables vs. The function returns
the teared equation so that if vTear is given, vSolved can be computed from eSolved
in a forward sequence (so solving eSolved[1] for vSolved[1], eSolved[2] for vSolved[2],
and so on). vTear must be selected, so that the equations eResidues are fulfilled.
Equations es are the union of eSolved and eResidue.
Variables vs are the union of vSolved and vTear.

Input argument td is an object of type TraverseDAG. Gsolvable defines the variables
that can be explicitly solved in every equation without influencing the solution space
(= rank preserving operation).

eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1]
"""
function tearEquations!(td::TraverseDAG, Gsolvable, es::Vector{Int}, vs::Vector{Int};
                        eSolvedFixed::Vector{Int}=Int[], vSolvedFixed::Vector{Int}=Int[], vTearFixed::Vector{Int}=Int[])
   vs2 = initAlgebraicSystem(td,es,vs,eSolvedFixed,vSolvedFixed,vTearFixed)
   #eResidue = fill(0,0)
   residue  = true
   esReduced = setdiff(es, eSolvedFixed)
   #println("    es = ", es, ", eSolvedFixed = ", eSolvedFixed, ", esReduced = ", esReduced)
   #println("    vs = ", vs, ", vSolvedFixed = ", vSolvedFixed)
   for eq in esReduced  # iterate only over equations that are not in eSolvedFixed
      residue = true
      for vj in Gsolvable[eq]
         if td.assign[vj]==0 && in_vs(td,vj)
            # vj is an element of vs that is not yet assigned
            # Add equation to graph
            td.assign[vj] = eq

            # Check for cycles
            if td.level[vj] == Undefined
               # (eq,vj) cannot introduce a cycle
               # Introduce a new level (the smallest level that exists yet)
               td.minlevel += -1
               td.level[vj] = td.minlevel

               # Inspect all childs and use level+1, if child has no level yet
               for v in td.G[eq]
                  if in_vs(td,v) && v!=vj &&
                     (td.level[v]==Undefined || td.level[v] <= td.level[vj]) # v is an element of td.vs and is not the variable to solve for and no level yet defined
                     setlevel(td,v,td.level[vj])
                  end
               end
               push!(vs2,vj)
               residue = false
               break # continue with next equation

            else # Traverse DAG starting from eq
               if visit!(td,vj)
                  # accept vj
                  push!(vs2,vj)
                  residue = false
                  break   # continue with next equation
               else
                  # cycle; remove vj from DAG and undo its changes
                  for vv in td.visitedStack
                     td.level[vv] = td.lastlevel[vv]
                  end
                  td.assign[vj] = 0
                  # continue with next variable in equation eq
               end
            end
         end
      end
      #if residue
      #   push!(eResidue, eq)
      #end
   end

   # Determine solved equations and variables
   (eSolved, vSolved) = sortDAG!(td, vs2)
   vTear = setdiff(vs,vSolved)
   eResidue = setdiff(es,eSolved)
   return (eSolved, vSolved, eResidue, vTear)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6be32f5e-05ce-4be3-b7af-64821076fda1
#hide
import Pkg

# ╔═╡ 0a803a48-8ad1-4d3c-9272-c0d1270d8e92
#hide
Pkg.activate(Base.current_project())

# ╔═╡ 00758500-175a-11ec-173c-670141a85b02
#export
begin
include(normpath(joinpath(@__FILE__,"..","..", "src/Export.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/Documenter.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/Makedocs.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/CodeRunner.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/Make.jl")))
end

# ╔═╡ 1ab8fabc-1b46-4cc3-9b58-03ac6b0fc9f9
#export
begin
notebook2script = Export.notebook2script
showdoc = Documenter.showdoc
build = Make.build
export notebook2script, showdoc, build
end

# ╔═╡ 9724fdac-8ffd-483c-9d7e-a4cdc8d0dd72
#hide
Export.notebook2script(joinpath("..", "nbs"), joinpath("..", "src"))

# ╔═╡ Cell order:
# ╠═6be32f5e-05ce-4be3-b7af-64821076fda1
# ╠═0a803a48-8ad1-4d3c-9272-c0d1270d8e92
# ╠═00758500-175a-11ec-173c-670141a85b02
# ╠═1ab8fabc-1b46-4cc3-9b58-03ac6b0fc9f9
# ╠═9724fdac-8ffd-483c-9d7e-a4cdc8d0dd72

ext: jl


text:
# LOAD IN PACKAGES AND DATA

using DataFrames, CSV, Plots, StatsPlots, StatsBase, Missings

# setup filepaths
path_source = string(@__DIR__,"\\..\\source");
path_dev = string(@__DIR__,"\\..\\dev");

df = DataFrame(CSV.read(joinpath(path_dev, "19-21VAERSCOMB_clean.csv", DataFrame));

size(df)

vax_ct = countmap(df[!, :VAX_NAME])

#SERIOUS_EVENT =1 is serious and SERIOUS_EVENT=0 is non-serious adverse event

#find most frequent serious adverse event - pull from SYMPTOM1, SYMPTOM2, SYMPTOM3, SYMPTOM4, SYMPTOM5
#find most frequent non-serious adverse event - pull from SYMPTOM1, SYMPTOM2, SYMPTOM3, SYMPTOM4, SYMPTOM5

df_nonserious = df[df.SERIOUS_EVENT .== 0, :]

df_serious = df[df.SERIOUS_EVENT .== 1, :]

# Frequency table of serious events
freq_serious = sort(freqtable(df, :SERIOUS_EVENT), rev=true);
freq_serious

# Frequency table of vaccines
freq_vax = sort(freqtable(df, :VAX_NAME), rev=true);
freq_vax

#function for aims analysis, to calc PRR from contingency tables created below

function get_freqtable(df, vax_name_str, symptom_str)
    test = select(df, ["VAX_NAME", "SYMPTOMS"])
    # Create dummy variable for vax_name
    test.VAX_IND = (test.VAX_NAME .== vax_name_str)
    # Create dummy variable for symptom
    test.SYMPTOM_IND = zeros(size(test, 1))
    for rownumber in 1:size(test, 1)
        if symptom_str in test[rownumber,:SYMPTOMS]
            test[rownumber,:SYMPTOM_IND] = 1
        end
    end
    # Create frequency table
    tbl = freqtable(test, :VAX_IND, :SYMPTOM_IND)
    return tbl
end;
function get_PRR(tbl)
    # Calculate PRR
    a = tbl[2,2]
    b = tbl[2,1]
    c = tbl[1,2]
    d = tbl[1,1]
    PRR = (a/(a+b))/(c/(c+d))
    return PRR
end;

# PRR for #1 freq VAX = "COVID19 (COVID19 (PFIZER-BIONTECH))" and serious report
df.COVID19_PFIZER = (df.VAX_NAME .== "COVID19 (COVID19 (PFIZER-BIONTECH))");
tbl = freqtable(df, :COVID19_PFIZER, :SERIOUS_EVENT)

#PRR for Pfizer COVID-19
get_PRR(tbl)

#PRR for #2 freq VAX = "Zoster" and serious report
df.ZOSTER = (df.VAX_NAME .== "ZOSTER (SHINGRIX)");
tbl = freqtable(df, :ZOSTER, :SERIOUS_EVENT)

#PRR for Zoster
get_PRR(tbl)

#PRR for #3 freq VAX = "Moderna" and serious report
df.COVID19_MODERNA = (df.VAX_NAME .== "COVID19 (COVID19 (MODERNA))");
tbl = freqtable(df, :COVID19_MODERNA, :SERIOUS_EVENT)

#PRR for Moderna COVID-19
get_PRR(tbl)

#PRR for #4 freq VAX = "Pneumo" and serious report
df.PNEUMO = (df.VAX_NAME .== "PNEUMO (PNEUMOVAX)");
tbl = freqtable(df, :PNEUMO, :SERIOUS_EVENT)

#PRR for Pneumo
get_PRR(tbl)

#PRR for #5 freq VAX = "JANSSEN COVID-19" and serious report
df.COVID19_JANSSEN = (df.VAX_NAME .=="COVID19 (COVID19 (JANSSEN))");
tbl = freqtable(df, :COVID19_JANSSEN, :SERIOUS_EVENT)

#PRR for Janssen COVID-19
get_PRR(tbl)


text:
facts("Testing only challenger reward - challenge accepted") do
    srand(1234)
    challenger = Agent(100)
    other = Agent(100)
    result = [false, [other], [challenger]]
    TcrCore.onlyChallengerRewardRedistribution(result, challenger, 10)
    @fact challenger.balance --> 70 #Received double deposit
end

facts("Testing punishment redistribution") do
    srand(1234)
    challenger = Agent(100)
    pro = [Agent(100), Agent(100)]
    cons = [Agent(100), Agent(100)]
    result = [false, pro, cons]
    TcrCore.punishmentRedistribution(result, challenger, 10)
    @fact challenger.balance --> 70 #Received double deposit
    @fact pro[1].balance --> 45
    @fact pro[2].balance --> 45
    @fact cons[1].balance --> 50
    @fact cons[2].balance --> 50
end

facts("Testing reward redistribution") do
    srand(1234)
    challenger = Agent(100)
    pro = [Agent(100), Agent(100)]
    cons = [Agent(100), Agent(100)]
    result = [false, pro, cons]
    TcrCore.rewardRedistribution(result, challenger, 10)
    @fact challenger.balance --> 60 #Received double deposit
    @fact pro[1].balance --> 50
    @fact pro[2].balance --> 50
    @fact cons[1].balance --> 55
    @fact cons[2].balance --> 55
end

facts("Testing punishment & reward redistribution") do
    srand(1234)
    challenger = Agent(100)
    pro = [Agent(100), Agent(100)]
    cons = [Agent(100), Agent(100)]
    result = [false, pro, cons]
    TcrCore.punishmentAndRewardRedistribution(result, challenger, 10)
    @fact challenger.balance --> 70 #Received double deposit
    @fact pro[1].balance --> 45
    @fact pro[2].balance --> 45
    @fact cons[1].balance --> 55
    @fact cons[2].balance --> 55
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
import PredictMD
import Test

Test.@test length(PredictMD.registry_url_list()) > 0


text:
#=
This file is auto-generated. Do not edit.
=#

"""Data Structure Operational Cost Data in two parts fixed and variable cost."""
mutable struct TwoPartCost <: OperationalCost
    variable::VariableCost
    fixed::Float64
    internal::InfrastructureSystemsInternal
end

function TwoPartCost(variable, fixed, )
    TwoPartCost(variable, fixed, InfrastructureSystemsInternal())
end

function TwoPartCost(; variable, fixed, )
    TwoPartCost(variable, fixed, )
end

# Constructor for demo purposes; non-functional.

function TwoPartCost(::Nothing)
    TwoPartCost(;
        variable=VariableCost((0.0, 0.0)),
        fixed=0.0,
    )
end

"""Get TwoPartCost variable."""
get_variable(value::TwoPartCost) = value.variable
"""Get TwoPartCost fixed."""
get_fixed(value::TwoPartCost) = value.fixed
"""Get TwoPartCost internal."""
get_internal(value::TwoPartCost) = value.internal


text:
using Distributions
using LinearAlgebra
using KalmanFilter
using Test

@testset "Dimensions check" begin
    M = 1:3
    P = 1:3
    for p in P
        for m in M
            y = rand(10, p)
            H = rand(p, m)
            A = rand(m, m)
            Q = Matrix(Diagonal(fill(1.0, m)))
            R = Matrix(Diagonal(fill(1.0, p)))
            x = rand(m)
            P = Matrix(Diagonal(fill(1.0, m)))
            model = Model(y, H, A, Q, R, x, P)

            # Test p dimension
            @test p == size(model.H, 1) == size(model.R, 1) == size(model.R, 2)
            # Test m dimension
            @test m == size(model.H, 2) == size(model.Q, 1) == size(model.Q, 2) == length(model.x) == size(model.P, 1) == size(model.P, 2) == size(model.A,1) == size(model.A, 2)
        end
    end
end

@testset "Univairate Local Level" begin
    T = 100
    Q = 0.1
    R = 4.0
    x = 0.0
    P = 1000.0
    y = rand(Normal(0.0, sqrt(R)), T)

    model = LocalLevel(y, Q, R, x, P)
    kf = kalmanfilter(model)
    ks = kalmansmoother(model, kf.priors, kf.posteriors)

    # Test posterior = smoother at time T
    @test isapprox(kf.posteriors.x[end, 1], ks.x[end, 1])
    @test isapprox(kf.posteriors.P[end, 1, 1], ks.P[end, 1, 1])

    # Test H and A
    @test isapprox(model.H[1, 1], 1.0)
    @test isapprox(model.A[1, 1], 1.0)
end


text:
using Dierckx

type BicubicSpline
  spline::Dierckx.Spline2D
end

BicubicSpline{T <: Real}(x::Vector{T}, y::Vector{T}, z::Matrix{T}) = BicubicSpline(Dierckx.Spline2D(x, y, z))

type NaturalCubicSpline{T <: Number} <: Interpolation
  x_vert::Vector{T}
  y_vert::Vector{T}
  b::Vector{T}
  c::Vector{T}
  d::Vector{T}
end

# adapted from here: http://sepwww.stanford.edu/sep/sergey/128A/answers6.pdf
function NaturalCubicSpline{T <: Real}(x_vert::Vector{T}, y_vert::Vector{T})
  n = length(x_vert)
  h = zeros(n - 1)
  b = zeros(n - 1)
  a = zeros(n - 2)
  g = zeros(n - 2)
  c = zeros(n)
  d = zeros(n - 1)

  for i = 1:n - 1
    h[i] = x_vert[i + 1] - x_vert[i]
    b[i] = (y_vert[i + 1] - y_vert[i]) / h[i]
  end

  for i = 1:n - 2
    a[i] = 2.0 * (h[i] + h[i + 1])
    # u[i] = 6 * (b[i + 1] - b[i])
    g[i] = b[i + 1] - b[i]
  end

  Alu = lufact(Tridiagonal(h[2:end-1], a[1:end], h[2:end-1]))
  c[2:end - 1] = Alu \ g

  for i=1:n - 1
    d[i] = (c[i+1] - c[i])/h[i]
    b[i] -= (2.0 * c[i] + c[i + 1]) * h[i]
    c[i] *= 3.0
  end

  return NaturalCubicSpline(x_vert, y_vert, b, c, d)
end

function call{T <: Real}(spl::NaturalCubicSpline, my_x::T)
  # get coefficients
  # b, c, d = gen_splines(x_vert, y_vert)

  x_idx = searchsortedlast(spl.x_vert, my_x)

  if x_idx == 0
    x_idx = 1
  elseif x_idx == length(spl.x_vert)
    x_idx = length(spl.x_vert) - 1
  end

  diff = my_x - spl.x_vert[x_idx]

  return spl.y_vert[x_idx] + diff*(spl.b[x_idx] + diff * (spl.c[x_idx] + (diff * spl.d[x_idx])))
end

## Cubic Interpolation ##
type CubicInterpolation{D <: DerivativeApprox, B1 <: BoundaryCondition, B2 <: BoundaryCondition} <: Interpolation
  derivativeApprox::D
  leftBoundaryCondition::B1
  rightBoundaryCondition::B2
  leftValue::Float64
  rightValue::Float64
  monotonic::Bool
  x_vals::Vector{Float64}
  y_vals::Vector{Float64}
  a::Vector{Float64}
  b::Vector{Float64}
  c::Vector{Float64}
  tmp::Vector{Float64}
  dx::Vector{Float64}
  S::Vector{Float64}
  n::Int
  L::TridiagonalOperator
end

function CubicInterpolation(dApprox::DerivativeApprox, leftBoundary::BoundaryCondition, rightBoundary::BoundaryCondition, x_vals::Vector{Float64}, y_vals::Vector{Float64},
                            leftValue::Float64 = 0.0, rightValue::Float64 = 0.0, monotonic::Bool = true)
  n = length(x_vals)
  ci = CubicInterpolation(dApprox, leftBoundary, rightBoundary, leftValue, rightValue, monotonic, x_vals, y_vals, zeros(n - 1), zeros(n - 1), zeros(n - 1), zeros(n), zeros(n - 1), zeros(n - 1),
                          n, TridiagonalOperator(n))

  update!(ci)
  return ci
end

call(interp::CubicInterpolation, x::Float64) = value(interp, x)

# type aliases #
typealias SplineCubicInterpolation{D, B1, B2} CubicInterpolation{Spline, B1, B2} # First derivative approximation

function cubicInterpolationPolynomialDerivative(x_vals::Vector{Float64}, y_vals::Vector{Float64}, x::Float64)
  a, b, c, d = x_vals
  u, v, w, z = y_vals

  return (-((((a-c)*(b-c)*(c-x)*z-(a-d)*(b-d)*(d-x)*w)*(a-x+b-x)+((a-c)*(b-c)*z-(a-d)*(b-d)*w)*(a-x)*(b-x))*(a-b)+
          ((a-c)*(a-d)*v-(b-c)*(b-d)*u)*(c-d)*(c-x)*(d-x)+((a-c)*(a-d)*(a-x)*v-(b-c)*(b-d)*(b-x)*u)
          *(c-x+d-x)*(c-d)))/((a-b)*(a-c)*(a-d)*(b-c)*(b-d)*(c-d))
end

function update_left_boundary_condition!(::Lagrange, interp::SplineCubicInterpolation)
  set_first_row!(interp.L, 1.0, 0.0)
  interp.tmp[1] = cubicInterpolationPolynomialDerivative(interp.x_vals[1:4], interp.y_vals[1:4], interp.x_vals[1])

  return interp
end

function update_right_boundary_condition!(::Lagrange, interp::SplineCubicInterpolation)
  set_last_row!(interp.L, 0.0, 1.0)
  interp.tmp[interp.n] = cubicInterpolationPolynomialDerivative(interp.x_vals[interp.n - 3:interp.n], interp.y_vals[interp.n - 3:interp.n], interp.x_vals[interp.n])

  return interp
end

function update!(interp::SplineCubicInterpolation)
  n = interp.n

  for i = 1:n - 1
    interp.dx[i] = interp.x_vals[i + 1] - interp.x_vals[i]
    interp.S[i] = (interp.y_vals[i + 1] - interp.y_vals[i]) / interp.dx[i]
  end

  for i = 2:n - 1
    set_mid_row!(interp.L, i, interp.dx[i], 2.0 * (interp.dx[i] + interp.dx[i - 1]), interp.dx[i - 1])
    interp.tmp[i] = 3.0 * (interp.dx[i] * interp.S[i - 1] + interp.dx[i - 1] * interp.S[i])
  end

  # update left boundary
  update_left_boundary_condition!(interp.leftBoundaryCondition, interp)

  # update right boundary
  update_right_boundary_condition!(interp.rightBoundaryCondition, interp)

  # solve
  solve_for!(interp.L, interp.tmp, interp.tmp)
  # interp.tmp = solve_for(interp.L, interp.tmp)

  if interp.monotonic
    corretion = 0.0
    pm = pu = pd = M = 0.0
    for i = 1:n
      if i == 1
        if interp.tmp[i] * interp.S[i] > 0.0
          correction = interp.tmp[i] / abs(interp.tmp[i]) * min(abs(interp.tmp[i]), abs(3.0 * interp.S[1]))
        else
          correction = 0.0
        end

        if correction != interp.tmp[i]
          interp.tmp[i] = correction
        end
      elseif i == n
        if interp.tmp[i] * interp.S[n - 1] > 0.0
          correction = interp.tmp[i] / abs(interp.tmp[i]) * min(abs(interp.tmp[i]), abs(3.0 * interp.S[n - 1]))
        else
          correction = 0.0
        end

        if correction != interp.tmp[i]
          interp.tmp[i] = correction
        end
      else
        pm = (interp.S[i - 1] * interp.dx[i] + interp.S[i] * interp.dx[i - 1]) / (interp.dx[i - 1] + interp.dx[i])
        M = 3.0 * min(min(abs(interp.S[i - 1]), abs(interp.S[i])), abs(pm))

        if i > 2
          if (interp.S[i - 1] - interp.S[i - 2]) * (interp.S[i] - interp.S[i - 1]) > 0.0
            pd = (interp.S[i - 1] * (2.0 * interp.dx[i - 1] + interp.dx[i - 2]) - interp.S[i - 2] * interp.dx[i - 1]) / (interp.dx[i - 2] + interp.dx[i - 1])

            if pm * pd > 0.0 && pm * (interp.S[i - 1] - interp.S[i - 2]) > 0.0
              M = max(M, 15 * min(abs(pm), abs(pd)))
            end
          end
        end
        if i < n - 1
          if (interp.S[i] - interp.S[i - 1]) * (interp.S[i + 1] - interp.S[i]) > 0.0
            pu = (interp.S[i] * (2.0 * interp.dx[i] + interp.dx[i + 1]) - interp.S[i + 1] * interp.dx[i]) / (interp.dx[i] + interp.dx[i + 1])
            if pm * pu > 0.0 && -pm * (interp.S[i] - interp.S[i - 1]) > 0.0
              M = max(M, 1.5 * min(abs(pm), abs(pu)))
            end
          end
        end

        if interp.tmp[i] * pm > 0.0
          correction = interp.tmp[i] / abs(interp.tmp[i]) * min(abs(interp.tmp[i]), M)
        else
          correction = 0.0
        end

        if correction != interp.tmp[i]
          interp.tmp[i] = correction
        end
      end
    end
  end

  # coefficients
  for i = 1:interp.n - 1
    interp.a[i] = interp.tmp[i]
    interp.b[i] = (3.0 * interp.S[i] - interp.tmp[i + 1] - 2.0 * interp.tmp[i]) / interp.dx[i]
    interp.c[i] = (interp.tmp[i + 1] + interp.tmp[i] - 2.0 * interp.S[i]) / (interp.dx[i] * interp.dx[i])
  end

  return interp
end

function value(interp::CubicInterpolation, x::Float64)
  j = locate(interp, x)
  dx = x - interp.x_vals[j]
  return interp.y_vals[j] + dx * (interp.a[j] + dx * (interp.b[j] + dx * interp.c[j]))
end


text:
abstract type FiniteElement{dim,shape,FuncOrder,GeomOrder} end
getshape(::FiniteElement{dim,shape}) where {dim,shape} = shape

####################
# Template for Finite Elements
####################

struct GenericFiniteElement{dim,shape,qorder,gorder} <: FiniteElement{dim,shape,qorder,gorder}
    func_basis::Interpolation{dim,shape,qorder}
    geom_basis::Interpolation{dim,shape,gorder}
    topology::Dict{Int,Int}
    geom_topology::Dict{Int,Int}
    M::Matrix{Float64}
end

function GenericFiniteElement(func_basis::Interpolation{dim,shape,qorder},
    geom_basis::Interpolation{dim,shape,gorder} = get_default_geom_interpolator(shape, Val{dim})) where {dim,shape<:AbstractRefShape, qorder,gorder}
    topology = gettopology(func_basis)
    geom_topology = gettopology(geom_basis)
    n_geom_basefuncs = getnbasefunctions(geom_basis)
    if gorder == 1
        M = one(Matrix{Float64}(undef,n_geom_basefuncs,n_geom_basefuncs))
    else
        M = _nodal_geom_data(geom_basis)
    end
    GenericFiniteElement{dim, shape, qorder, gorder}(func_basis, geom_basis, topology, geom_topology, M)
end

function _nodal_geom_data(geom_interpol::Lagrange{dim,shape,order}) where {dim,shape,order}
    # Matrix to get spacial coordinates
    nodal_points, topology = get_nodal_points(shape, Val{dim}, order)
    T = eltype(nodal_points[1])
    qrs = QuadratureRule{dim,shape,T}(fill(T(NaN), length(nodal_points)), nodal_points) # weights will not be used
    n_qpoints = length(getweights(qrs))
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints)
    for (qp, ξ) in enumerate(qrs.points)
        for i in 1:n_geom_basefuncs
            M[i, qp] = value(geom_interpol, i, ξ)
        end
    end
    M
end

@inline getnbasefunctions(fe::GenericFiniteElement) = getnbasefunctions(fe.func_basis)
@inline getngeombasefunctions(fe::GenericFiniteElement) = getnbasefunctions(fe.geom_basis)
@inline gettopology(fe::GenericFiniteElement) = fe.topology
@inline getgeomtopology(fe::GenericFiniteElement) = fe.geom_topology

"""
value(ip::ContinuousLagrange{dim,shape<:AbstractRefShape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Continuous Lagrange Finite Element basis `j` at point ξ on shape `shape`
"""
function value(ip::GenericFiniteElement{dim}, k::Int, ξ::Vec{dim,T}) where {dim, T}
    value(ip.func_basis, k, ξ)
end

"""
geom_value(ip::ContinuousLagrange{dim,shape<:AbstractRefShape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Continuous Lagrange Finite Element basis `j` at point ξ on shape `shape`
"""
function geom_value(ip::GenericFiniteElement{dim}, k::Int, ξ::Vec{dim,T}) where {dim, T}
    value(ip.geom_basis, k, ξ)
end

"""
gradient_value(ip::FiniteElement{dim,shape<:AbstractRefShape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Finite Element basis `j` derivative at point ξ
on the reference shape `shape`
"""
function gradient_value(ip::FiniteElement{dim}, k::Int, ξ::Vec{dim,T}) where {dim, T}
    gradient(ξ -> value(ip, k, ξ), ξ)
end

function spatial_nodal_coordinate(fe::FiniteElement{dim}, n_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_base_funcs = getngeombasefunctions(fe)
    @assert length(x) == n_base_funcs
    vec = zero(Vec{dim,T})
    @inbounds for i in 1:n_base_funcs
        vec += fe.M[i, n_point] * x[i]
    end
    return vec
end

function spatial_nodal_coordinate(fe::FiniteElement{dim,shape,1,1}, n_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,shape,T}
    @assert length(x) ==  getngeombasefunctions(fe)
    return x[n_point]
end


text:
using DashBootstrapComponents

inputs = html_div([
    dbc_input(placeholder = "Valid input...", valid = true, className = "mb-3"),
    dbc_input(placeholder = "Invalid input...", invalid = true),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
include("euler/euler.jl")
using .Numbers: get_digits

function compute(n::Float64)::Int
    bouncy = 0
    i = 1
    while bouncy / i < n
        i += 1
        digits = get_digits(i)
        sorted_digits = sort(digits)
        if !(digits == sorted_digits || digits == reverse(sorted_digits)) bouncy += 1 end
    end
    return i
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
export gripenberg

struct BranchState{MT<:AbstractMatrix, ET}
    A::MT
    seq::Vector{ET}
    mode::Int
    p::Float64 # min_j ||prod_{i=1}^j A_i||^{1/j}
end

"""
    gripenberg(s::AbstractDiscreteSwitchedSystem; δ=1e-2,
                    max_eval = 10000, max_ρ_eval = max_eval,
                    max_norm_eval = max_eval, max_length = 50,
                    matrix_norm = A -> opnorm(A, 2), verbose = 1)

Gripenberg algorithm [G96] for computing an upper bound `ub` and a lower bound
`lb` to the Joint Spectral Radius such that `ub - lb ≤ 1e-2`.

[G96] Gripenberg, G. Computing the joint spectral radius.
*Linear Algebra and its Applications*, *Elsevier*, **1996**, *234*, 43-60
"""
function gripenberg(s::AbstractDiscreteSwitchedSystem; δ=1e-2,
                    max_eval = 10000, max_ρ_eval = max_eval,
                    max_norm_eval = max_eval, max_length = 50,
                    matrix_norm = A -> opnorm(A, 2), verbose = 1)
    MT = typeof(dynamicfort(s, first(out_transitions(s, 1))))
    ET = transitiontype(s)
    branches = [BranchState{MT, ET}(Matrix(LinearAlgebra.I, HybridSystems.statedim(s, mode), HybridSystems.statedim(s, mode)), ET[], mode, Inf) for mode in modes(s)]
    smp = nothing
    ub = Inf
    n_ρ_eval = 0
    n_norm_eval = 0
    cur_length = 0
    while !isempty(branches) && n_ρ_eval < max_ρ_eval &&
        n_norm_eval < max_norm_eval && cur_length < max_length &&
        (smp === nothing || ub > smp.growthrate + δ)
        @assert cur_length == length(first(branches).seq)
        cur_length += 1
        new_branches = eltype(branches)[]
        for branch in branches
            for t in out_transitions(s, branch.mode)
                A = dynamicfort(s, t) * branch.A
                seq = [branch.seq; t]
                if source(s, first(seq)) == target(s, last(seq))
                    new_smp = periodicswitching(s, seq, A)
                    n_ρ_eval += 1
                    notifyperiodic!(s, new_smp)
                    if smp === nothing || isbetter(new_smp, smp)
                        smp = new_smp
                    end
                end
                p = min(matrix_norm(A)^(1/length(seq)), branch.p)
                n_norm_eval += 1
                if smp === nothing || p > smp.growthrate + δ
                    push!(new_branches, BranchState(A, seq, target(s, t), p))
                end
            end
        end
        branches = new_branches
        β = 0.0
        if !isempty(branches)
            β = max(β, maximum(branch -> branch.p, branches))
        end
        if smp !== nothing
            β = max(smp.growthrate + δ, β)
        end
        ub = min(ub, β)
    end
    if verbose ≥ 1
        function _keyword_feedback(cur, max)
            if cur < max
                print(cur, " < ", max)
            else
                printstyled(cur, " ≥ ", max, bold=true, color=:red)
            end
        end
        print("ρ evaluations   : ")
        _keyword_feedback(n_ρ_eval, max_ρ_eval)
        println(" = max_ρ_eval.")
        print("norm evaluations: ")
        _keyword_feedback(n_norm_eval, max_norm_eval)
        println(" = max_norm_eval.")
        print("switch length   : ")
        _keyword_feedback(cur_length, max_length)
        println(" = max_length.")
    end
    return smp, ub
end


text:
# #
# Reproduce the figures in
# [V. Duchêne and C. Klein](http://dx.doi.org/10.3934/dcdsb.2021300)
# on the Serre-Green-Naghdi and its fully dispersive counterpart.
# #
export PlotSolitaryWaveWGN1,PlotSolitaryWaveWGN2,PlotSolitaryWaveWGN3,PlotJacobianWGN,IntegrateSolitaryWaveWGN,StabilitySolitaryWaveWGN,IntegrateWGN
using WaterWaves1D,FFTW,Plots,LinearAlgebra,ProgressMeter;
#using JLD #(uncomment if using function save)

#---- Figures 1 and 2
"""
	PlotSolitaryWaveWGN1(;kwargs)

First method to compute the WGN solitary wave.
Use GMRES-iterative method for the inversion of the Jacobian.

Arguments are all optional:
- `c` the velocity,
- `N` the number of collocation points,
- `L` the half-length of the mesh,
- `verbose` provides details on residuals if `true` (default is `false`)
- `name`: a string used to save the figure as a name.pdf.

Return `(η,u,v,mesh)`, where `mesh.x` are the collocation points.
"""
function PlotSolitaryWaveWGN1(;c=2,N=2^10,L=10*π,verbose=false,name=nothing)
	param = ( μ  = 1,
			ϵ  = 1,
        	N  = N,
            L  = L,
			c  = c
				)

	(η,u,v,mesh) = SolitaryWaveWhithamGreenNaghdi(
				param;
				method=2, α = 0,
				tol =  1e-14, max_iter=10,
				ktol =0*1e-11, gtol = 1e-16,
				iterative = true, q=1,
				verbose = verbose)

	(ηGN,uGN) = SolitaryWaveSerreGreenNaghdi(param)

	plt = plot(layout=(1,2))
	plot!(plt[1,1], mesh.x, [u uGN];
	  title=string("c=",c),
	  xlabel = "x",
	  ylabel = "u",
	  label=["WGN" "SGN"])
	plot!(plt[1,2], fftshift(mesh.k),
	  [log10.(abs.(fftshift(fft(u)))) log10.(abs.(fftshift(fft(uGN))))];
	  title="frequency",
	  label=["WGN" "SGN"])
	display(plt)
	if name != nothing savefig(string(name,".pdf")); end
	return (η,u,v,mesh)
end

#---- Figure 3
"""
	PlotSolitaryWaveWGN2(;kwargs)

Second method to compute the WGN solitary wave. Same as `PlotSolitaryWaveWGN1`, but
use non-iterative method for the inversion of the Jacobian.
To be used with higher values of the velocity (default is `c=20`).
"""
function PlotSolitaryWaveWGN2(;c=20,L=10*π,N=2^10,verbose=false,name=nothing)
	param = ( μ  = 1,
		ϵ  = 1,
      	N  = N,
        L  = L,
		c = c )

	(η,u,v,mesh) = SolitaryWaveWhithamGreenNaghdi(
				param;
				method=2, α = 1, #ici α = 1 évite des oscillations important si c = 3 ou c = 20
				tol =  1e-14, max_iter=15,
				iterative = false, q=1,
				verbose = verbose)

	(ηGN,uGN) = SolitaryWaveSerreGreenNaghdi(param)

	plt = plot(layout=(1,2))
	plot!(plt[1,1], mesh.x, [u uGN];
	  title=string("c=",c),
	  xlabel = "u",
	  ylabel = "u/c",
	  label=["WGN" "SGN"])
	plot!(plt[1,2], fftshift(mesh.k),
	  [log10.(abs.(fftshift(fft(u)))) log10.(abs.(fftshift(fft(uGN))))];
	  title="frequency",
	  label=["WGN" "SGN"])
    display(plt)
	if name != nothing savefig(string(name,".pdf")); end
	return (η,u,v,mesh)
end
#---- Figure 4
"""
	PlotSolitaryWaveWGN3(;kwargs)

Third method to compute the WGN solitary wave. Same as `PlotSolitaryWaveWGN1`, but
- use non-iterative method for the inversion of the Jacobian; and
- use a rescaled equation.
To be used with highest values of the velocity (default is `c=100`).
"""
function PlotSolitaryWaveWGN3(;c=100,L=10*π,N=2^10,verbose=false,name=nothing)
	param = ( μ  = 1,
		ϵ  = 1,
      	N  = N,
        L  = L,
		c = c )

	(η,u,v,mesh) = SolitaryWaveWhithamGreenNaghdi(
				param;
				method=3, α = 1,
				tol =  1e-10, max_iter=10,
				iterative = false, q=1,
				verbose = verbose)
	(ηGN,uGN) = SolitaryWaveSerreGreenNaghdi(param)

	plt = plot(layout=(1,2))
	plot!(plt[1,1], mesh.x, [u/c uGN/c];
	  title=string("c=",c),
	  xlabel = "x",
	  ylabel = "u/c",
	  label=["WGN" "SGN"])
	plot!(plt[1,2], fftshift(mesh.k),
	  [log10.(abs.(fftshift(fft(u/c)))) log10.(abs.(fftshift(fft(uGN/c)))) ];
	  title="frequency",
	  label=["WGN" "SGN"])
	display(plt)
	if name != nothing savefig(string(name,".pdf")); end
	return (η,u,v,mesh)
end

#------ Figure 6
"""
	PlotJacobianWGN(;kwargs)

Compute the Jacobian matrix to be inverted in SGN or WGN equation.

Arguments are all optional:
- `c` the velocity,
- `N` the number of collocation points,
- `L` the half-length of the mesh,
- `SGN`: use SGN if `true`, and WGN if `false` (default is `false`),
- `name`: a string used to save the figure as a name.pdf.

"""
function PlotJacobianWGN(;c=20,L=10*π,N=2^10,SGN=false,verbose=false,name=nothing)
	ϵ,μ,α=1,1,0

	(η,u,v,mesh) = SolitaryWaveWhithamGreenNaghdi(
				(c=c,L=L,N=N,μ=μ,ϵ=ϵ); SGN = SGN,
				method=2, α = 1, #ici α = 1 évite des oscillations important si c = 3 ou c = 20
				tol =  1e-12, max_iter=15,
				iterative = false, q=1,
				verbose = verbose)

	k,x=mesh.k,mesh.x
	if SGN == true
		F₀ = sqrt(μ)*1im * k
	else
		F₁ = tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k))
		F₁[1] = 1
		F₀ = 1im * sqrt.(3*(1 ./F₁ .- 1)).*sign.(k)
	end
	x₀ = mesh.x[1]
	FFT = exp.(-1im*k*(x.-x₀)');
	IFFT = exp.(1im*k*(x.-x₀)')/length(x);
	M₀ = IFFT * Diagonal( F₀ )* FFT
	M(v) = Diagonal( v )

	dxu = real.(ifft(1im*k.*fft(u)))
	dxu ./= norm(dxu,2)
	hu = c ./(c .- ϵ*u)
	Fu = hu.* real.(ifft(F₀.*fft(u)))
	F2u = real.(ifft(F₀.*fft(hu.^2 .* Fu ) ))
	Du = c ./hu .+ 2*ϵ/3 * F2u ./ hu .+ ϵ^2/c * hu .* Fu.^2 .- (hu.^2)/c
	Jac = (-1/3 *M(1 ./ hu.^2 )* M₀ * M(hu.^3)* (c*M₀ .+ 3*ϵ * M( Fu ))
					.+ ϵ * M( hu .* Fu ) * M₀
					.+ M( Du ) .+ α*M( hu.^2 )*dxu*dxu' *M(1 ./ hu.^2 ) )
	Jacstar = -1/3 *M(1 ./ hu.^2 )* M₀ * M(hu.^3)* c*M₀

	plt = plot(layout=(1,2))
	surface!(plt[1,1],fftshift(k),fftshift(k)[N:-1:1],
		log10.(abs.(FFT*Jac*IFFT)),
		title = "Jacobian")
	surface!(plt[1,2],fftshift(k),fftshift(k)[N:-1:1],
		log10.(abs.(FFT*Jacstar*IFFT)),
		title = "non-diagonal part")
	display(plt)

	if name != nothing savefig(string(name,".pdf")); end
	return (Jac,Jacstar,FFT,IFFT)
end

#---- Figures 7 and 8
"""
	IntegrateSolitaryWaveWGN(;kwargs)

Integrate in time the SGN or WGN solitary wave.

Arguments are all optional:
- `c` the velocity,
- `N` the number of collocation points,
- `L` the half-length of the mesh,
- `T` the final time of integration,
- `dt` the timestep,
- `SGN`: use SGN if `true`, and WGN if `false` (default is `false`),
- `name`: a string used to save raw data and the figure.

Return `problem` of type `Problem`, containing all the information.
"""
function IntegrateSolitaryWaveWGN(;SGN=false,c=2,N=2^10,L=10*π,T=1,dt=1/2000,name=nothing)
	if name != nothing ns=floor(Int,max(1,T/dt/100)) else ns=1 end

	param = ( μ  = 1, ϵ  = 1, c = c,
				N  = N,
	            L  = L,
	            T  = T,
	            dt = dt, ns=ns)
	if SGN == true
		(ηinit,uinit,vinit,mesh) = SolitaryWaveSerreGreenNaghdi(param)
		(ηfin,ufin,vfin) = SolitaryWaveSerreGreenNaghdi(param; x₀ = c*T)
	else
		(ηinit,uinit,vinit,mesh) = SolitaryWaveWhithamGreenNaghdi(
				param; method = 3, max_iter=10,α=1,verbose=true, tol=1e-12)
	  	(ηfin,ufin,vfin) = SolitaryWaveWhithamGreenNaghdi(
				param; x₀ = c*T, method = 3, max_iter=10,α=1,verbose=true, tol=1e-12)
	end
	init     = Init(mesh,ηinit,vinit)
	model = WhithamGreenNaghdi(param;SGN=SGN, ktol=0, iterate=true, precond = true)
	problem = Problem(model, init, param)
	solve!( problem )

	(ηcomp,vcomp,ucomp)   =  model.mapfrofull(last(problem.data.U))

	E(η,u,v) = sum(η.^2 .+ (1 .+ param.ϵ*η).*u.*v)
	dE(η1,u1,v1,η2,u2,v2) = sum(η1.^2-η2.^2) + sum((1 .+ param.ϵ*η1).*u1.*v1 - (1 .+ param.ϵ*η2).*u2.*v2)

	print(string("absolute energy difference: ",dE(ηinit,uinit,vinit,ηcomp,ucomp,vcomp),"\n"))
	print(string("relative energy difference: ",dE(ηinit,uinit,vinit,ηcomp,ucomp,vcomp)/E(ηinit,uinit,vinit),"\n"))

	plt = plot(layout=(1,2))
	plot!(plt[1,1],fftshift(mesh.k),fftshift(log10.(abs.(fft(uinit))));
			title = "frequency",
			label = "initial",
			xlabel="x",
			ylabel="u")
	plot!(plt[1,1],fftshift(mesh.k),fftshift(log10.(abs.(fft(ucomp))));
			label = "final")
	plot!(plt[1,2],mesh.x,ufin-ucomp;
			title = string("error at time t=",problem.times.tfin),
			label="difference",
			xlabel="x",
			ylabel="δu")
	display(plt)

	if name != nothing
		savefig(string(name,".pdf"));
		create_animation(problem;name=string(name,"-anim.pdf"))
		plot_solution(problem)
		savefig(string(name,"-final.pdf"));
		#save(problem,name);
	end
	return problem
end

#------ Figures 9 to 15

"""
	StabilitySolitaryWaveWGN(;kwargs)

Integrate in time a pertubed SGN or WGN solitary wave.

Arguments are all optional:
- `p` (1,2, or a real) is the type of perturbation:
	1. If `p=1`, then `u` is multiplied by `λ = 0.99`
	2. If `p=2`, then `u` is multiplied by `λ = 1.01`
	3. Otherwise, one adds `p exp(-x^2)` to `u`
- `c` the velocity,
- `N` the number of collocation points,
- `L` the half-length of the mesh,
- `T` the final time of integration,
- `dt` the timestep,
- `SGN`: use SGN if `true`, and WGN if `false` (default is `false`),
- `iterate`: use GMRES if `true`, and LU decomposition otherwise (default is `true`),
- `precond` gives some choice in the preconditioner for GMRES,
- `dealias` use 2/3 dealiasing rule if `1`, (default is `0`, i.e. no dealiasing),
- `name`: a string used to save raw data and the figures.

Return `problem` of type `Problem`, containing all the information.
"""
function StabilitySolitaryWaveWGN(;p=2,c=2,N=2^10,L=10*π,T=10,dt=10/10^4,SGN=false,precond=true,iterate=true,dealias=0,name=nothing)
	if name != nothing ns=floor(Int,max(1,T/dt/100)) else ns=1 end
	if p == 1
		λ = 0.99
	elseif p == 2
		λ = 1.01
	else
		λ = p
	end

	param = ( μ = 1, ϵ = 1, c = c, λ = λ,
				N  = N,
	            L  = L,
	            T  = T,
	            dt = dt, ns=ns,
				dealias = dealias,
				SGN=SGN,
				ktol=0*1e-6,
				gtol=1e-12,
				iterate=iterate)

	if SGN == true
		(η,u,v,mesh) = SolitaryWaveSerreGreenNaghdi(param)
	else
		(η,u,v,mesh) = SolitaryWaveWhithamGreenNaghdi(
				param; method = 3, tol=1e-14, max_iter=10,α=1,verbose=true)
	end
	k = mesh.k

	if precond > 0
		precond = Diagonal( 1 .+ μ/3*(precond^2*k).^2 )
	elseif precond < 0
		precond = Diagonal( (1 .+ μ/3000*k.^8)  )
	end
	if p == 1 || p == 2
		u= λ*u
	else
		u .+= λ*exp.(-mesh.x.^2)
	end
	if SGN == true
		F₀ = sqrt(param.μ)*1im * k
	else
		F₁ = tanh.(sqrt(param.μ)*abs.(k))./(sqrt(param.μ)*abs.(k))
		F₁[1] = 1
		F₀ = 1im * sqrt.(3*(1 ./F₁ .- 1)).*sign.(k)
	end
	DxF(v) = real.(ifft(F₀ .* fft(v)))
	h = 1 .+ param.ϵ*η
	v = u - 1/3 ./h .* (DxF(h.^3 .*DxF(u)))

	init = Init(mesh,η,v)
	model = WhithamGreenNaghdi(param;SGN=param.SGN, dealias = param.dealias, ktol=param.ktol, gtol=param.gtol, iterate=param.iterate, precond = precond)
	problem = Problem(model, init, param)
	solve!( problem )

	#if name != nothing save(problem,name); end

	(ηfin,vfin,ufin)   =  model.mapfrofull(last(problem.data.U))

	E(η,u,v) = sum(η.^2 .+ (1 .+ param.ϵ*η).*u.*v)
	dE(η1,u1,v1,η2,u2,v2) = sum(η1.^2-η2.^2) + sum((1 .+ param.ϵ*η1).*u1.*v1 - (1 .+ param.ϵ*η2).*u2.*v2)

	print(string("normalized error: ",dE(η,u,v,ηfin,ufin,vfin)/E(η,u,v),"\n"))

	plt = plot(layout=(1,2))
	plot!(plt[1,1],fftshift(mesh.k),fftshift(log10.(abs.(fft(η))));
			title = "frequency",
			label = "initial")
	plot!(plt[1,1],fftshift(mesh.k),fftshift(log10.(abs.(fft(ηfin))));
			label = "final")
	plot!(plt[1,2],mesh.x,[η ηfin];
			title = string("at time t=",problem.times.tfin),
			label=["zeta initial" "zeta final"])
	plot!(plt[1,2],mesh.x,[u ufin];
			label=["u initial" "u final"])
	display(plt)
	if name != nothing savefig(string(name,"-final.pdf")); end

	ts = problem.times.ts
	x = mesh.x
	k = mesh.k
	X=interpolate(mesh,real.(ifft(problem.data.U[1][:,1])))[1].x
	zs=zeros(length(X),length(ts));
	for i in 1:length(ts)
		zs[:,i].=interpolate(mesh,real.(ifft(problem.data.U[i][:,1])))[2]
	end

	plt = plot()
	scatter!(plt,ts,maximum(abs.(zs),dims=1)',
			title="maximum of surface deformation",
			label="",
			xlabel="time t")
	display(plt)

	if name != nothing
		savefig(string(name,"-znorm.pdf"));

		us=zeros(length(X),length(ts));
		@showprogress 1 "Computing v..." for i in 1:length(ts)
			us[:,i].=interpolate(mesh,model.mapfro(problem.data.U[i])[2])[2]
		end
		plt = plot()
		scatter!(plt,ts,maximum(abs.(us),dims=1)',
				title="maximum of velocity",
				label="",
				xlabel="time t")
		display(plt)
		savefig(string(name,"-vnorm.pdf"));

		plt=plot()
		my_cg = cgrad([:blue,:green])
		#surface!(plt,X,ts,zs',view_angle=(20,30), color = my_cg)
		display(plt)
		savefig(string(name,"-evol.pdf"));

		create_animation(problem;name=string(name,"-anim"))
	end
	return problem
end

#---- Figures 16 to end
"""
	IntegrateWGN(scenario;kwargs)

Integrate in time SGN or WGN with an initial data depending on a given `scenario`
- If `scenario = 1`, produces a dispersive shock wave for a unidirectional wave constructed from the Saint-Venant model.
- If `scenario = 2`, produces a dispersive shock wave for a unidirectional wave constructed from the Camassa-Holm model.
- If `scenario = 3`, the initial data challenges the non-cavitation assumption.

Other arguments are optional:
- `δ` is either
	- the square root of the shallowness parameter (the nonlinear parameter `ϵ=1`) if `scenario=1` or `2`;
	- the minimal depth if `scenario=3`
- `N` the number of collocation points,
- `L` the half-length of the mesh,
- `T` the final time of integration,
- `dt` the timestep,
- `SGN`: use SGN if `true`, and WGN if `false` (default is `false`),
- `iterate`: use GMRES if `true`, and LU decomposition otherwise (default is `true`),
- `precond` gives some choice in the preconditioner for GMRES,
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `name`: a string used to save raw data and the figures.

Return `(problem,plt)` where `problem` contains all the information and `plt` a plot of the final time solution.
"""
function IntegrateWGN(scenario;δ=0.1,N=2^11,L=3*π,x₀=-3,T= 5,dt = 5/10^4,SGN=false,dealias=0,iterate=true,precond=true,name=nothing)

	if name != nothing ns=floor(Int,max(1,T/dt/100)) else ns=1 end
	if scenario == 1 || scenario == 2
		μ  = δ^2
	else
		μ = 1
	end
	param = ( μ  = μ, ϵ  = 1,
				N  = N, L  = L,
	            T  = T, dt = dt,
				ns=ns )

	mesh=Mesh(param)
	k=mesh.k
	if precond > 0
		precond = Diagonal( 1 .+ μ/3*(precond^2*k).^2 )
	elseif precond < 0
		precond = Diagonal( (1 .+ μ/3000*k.^8)  )
	end

	if scenario == 1
		η= exp.(-(mesh.x .-x₀).^2)
		u= 2*sqrt.(1 .+ param.ϵ*η) .-2
	elseif scenario == 2
		function krasny!(v)
			v[abs.(v).<1e-14].=0
			return v
		end
		Dx(v) = real.(ifft( 1im*mesh.k.* krasny!(fft(v))))
		Dx2(v) = Dx(Dx(v))
		ϵ = param.ϵ;μ=param.μ;
		w = - (mesh.x.-x₀).* exp.(-(mesh.x.-x₀).^2)
		u = w .+ μ/12 *Dx2(w) .+ μ*ϵ/6* w.*Dx2(w)
		η = u .+ ϵ/4*u.^2 .- μ/6* Dx2( u+3*ϵ/4* u.^2) .- μ*ϵ/6 * u .* Dx2(u) .- 5*μ*ϵ/48 * Dx(u).^2
	elseif scenario == 3
		η = -(1-δ)*exp.(-mesh.x.^2)
		u = -mesh.x.*exp.(-mesh.x.^2)
	else
		error("the first argument must be 1, 2 or 3")
	end
	if SGN == true
		F₀ = sqrt(param.μ)*1im*mesh.k
	else
		F₁ 	= tanh.(sqrt(param.μ)*abs.(mesh.k))./(sqrt(param.μ)*abs.(mesh.k))
		F₁[1]= 1                 # Differentiation
		F₀   = 1im * sqrt.(3*(1 ./F₁ .- 1)).*sign.(mesh.k)
	end
	DxF(v) = real.(ifft(F₀ .* fft(v)))
	h=1 .+ param.ϵ*η
	v= u - 1/3 ./h .* (DxF(h.^3 .*DxF(u)))

	init     = Init(mesh,η,v)
	model = WhithamGreenNaghdi(param;SGN=SGN, ktol=0, gtol=1e-12, iterate=iterate, precond = precond, dealias = dealias)
	problem = Problem(model, init, param)
	solve!( problem )

	(ηfin,vfin,ufin)   =  model.mapfrofull(last(problem.data.U))
	E(η,u,v) = sum(η.^2 .+ (1 .+ param.ϵ*η).*u.*v)
	dE(η1,u1,v1,η2,u2,v2) = sum(η1.^2-η2.^2) + sum((1 .+ param.ϵ*η1).*u1.*v1 - (1 .+ param.ϵ*η2).*u2.*v2)

	print(string("normalized error: ",dE(η,u,v,ηfin,ufin,vfin)/E(η,u,v),"\n"))

	fftηfin=last(problem.data.U)[:,1]

	plt = plot(layout=(1,2))
	plot!(plt[1,1],fftshift(mesh.k),fftshift(log10.(abs.(fftηfin)));
			title = "frequency",label="")
	plot!(plt[1,2],mesh.x,real.(ifft(fftηfin));
			title = string("surface deformation at time t=",problem.times.tfin),
			label="")
	display(plt)

	if name != nothing
		savefig(string(name,".pdf"));
		#save(problem,name);
		ts = problem.times.ts
		x = mesh.x
		k = mesh.k
		zs=zeros(param.N,length(ts));
		@showprogress 1 for i in 1:length(ts)
			zs[:,i].=real.(ifft(problem.data.U[i][:,1]))
		end
		plt3=plot()
		my_cg = cgrad([:blue,:green])
		surface!(plt3,x,ts,zs',view_angle=(20,30), color = my_cg)
		display(plt3)
		savefig(string(name,"-evol.pdf"));
		create_animation(problem;ylims=false,name=string(name,"-anim"))

	end
	display(plt)
	return problem,plt
end
nothing


text:
@testset "Norm of polynomial" begin
    Mod.@polyvar x y
    @test norm(x) == norm(x, 2)  == 1.0
    @test norm(2x) == norm(2x, 2) == 2.0
    @test norm(2x+1) == norm(2x+1, 2) ≈ sqrt(5)
    @test norm(2x*y+x^2+3) ≈ norm(2x*y+x^2+3,2) ≈ sqrt(14)
end


text:
module TestTasks

using SimpleTasks.Types

import SimpleTasks.Tasks.DaemonTask
import SimpleTasks.Tasks.BasicTask
import SimpleTasks.Services.Datasource

export TEST_ID, TEST_TASK_NAME, TEST_BASE_DIRECTORY, TEST_INPUTS
export make_valid_basic_info

export MockTaskNoExecute, MockTaskExecute
export make_valid_task_execute, make_valid_task_no_execute

const TEST_ID = 1
const TEST_TASK_NAME = "test_name"
const TEST_BASE_DIRECTORY = "base_directory"
const TEST_INPUTS = ["input_1", "input_2"]
function make_valid_basic_info()
    return BasicTask.Info(TEST_ID, TEST_TASK_NAME, TEST_BASE_DIRECTORY,
        TEST_INPUTS)
end

type MockTaskNoExecute <: DaemonTaskDetails
    basic_info::BasicTask.Info
    payload_info::AbstractString
end

type MockTaskExecute <: DaemonTaskDetails
    basic_info::BasicTask.Info
    payload_info::AbstractString
end
# register a registered test task
function DaemonTask.execute(task::MockTaskExecute,
        datasource::DatasourceService)
    println("executing")
end

function make_valid_task_no_execute()
    return MockTaskNoExecute(make_valid_basic_info(),
        "Mock Task No Execute")
end

function make_valid_task_execute()
    return MockTaskExecute(make_valid_basic_info(),
        "Mock Task Execute")
end

end # module MockTasks

max_forks_repo_name: seung-lab/SimpleTasks


text:
using Random
using Parameters
using Debugger

Random.seed!(1)

GOAL = "goal"
COLLISION = "collision"
ONGOING = "ongoing"

abstract type Mdp end

@with_kw mutable struct ActMdp <: Mdp
	s_max::Float64 = 200.0
	sd_max::Float64 = 30.0
	gamma::Float64 = 1.0

	dt::Float64 = 0.250 # ms
	tMax::Float64 = 25 # secs

	goal::Vector{Float64} = [200.0, 30.0] # column vector
	obstacles::Array{Tuple{Float64, Float64}} = [] # list of (s,t) tuples
	dcol::Float64 = 10.0 # meters

	start_state::Vector{Float64} = [0.0, 20.0, 0.0] # s, sd (column vector)
	state::Vector{Float64} = [0.0, 20.0] # s, sd (column vector)
	actions::Vector{Float64} = [-4, -2, -1, 0, 1, 2]
	# Transtion model
	Ts::Array{Float64,2} = [1.0 dt 0.0; 0  1 0; 0 0 1]
	Ta::Array{Float64,1} = [dt^2/2; dt; 0]
	Tt::Array{Float64,1} = [0.0; 0; dt]

end

"""
mdp init
"""
function init!(m::ActMdp;
				s_max::Float64 = 200.0,
				sd_max::Float64 = 30.0,
				gamma::Float64 = 1.0,
				dt::Float64 = 0.25,
				tMax::Float64 = 30.0,
				goal::Vector{Float64} = [200.0, 30.0],
				obstacles=[],
				dcol::Float64 = 10.0,
				start_state::Vector{Float64} = [0.0, 20.0, 0.0],
				actions::Vector{Float64} = [-4.0, -2, -1, 0, 1, 2])

	m.s_max = s_max
	m.sd_max = sd_max
	m.gamma = gamma

	m.dt = dt
	m.tMax = tMax

	m.goal[1], m.goal[2] = goal[1], goal[2]
	m.obstacles = deepcopy(obstacles)
	m.dcol = dcol

	m.start_state = start_state
	m.state = start_state
	m.actions = actions
end

"""
mdp startState
"""
function startState(m::ActMdp)
	#return Tuple(m.start_state)
	return m.start_state
end

"""
mdp isEnd
"""
function isEnd(m::ActMdp, s::Array{Float64,1})::Int64
	ret = 0
	if s[1] >= m.goal[1]
		ret = 2 # true, "goal"
	elseif s[2] < 0 || s[3] > m.tMax
		ret = 1 # true, "collision" # handle driving backward as collision
	elseif s[3] > m.obstacles[end][1] # small optim
		ret = 0 # false, "ongoing"
	else
		for obstacle in m.obstacles
			#s_cross, t_cross = obstacle
			t_cross, s_cross = obstacle
			if s[3] == t_cross
				#println("==========> DBG: t_cross=$t_cross s_cross=$s_cross s_ego=$(s[1])")
				if abs(s_cross - s[1]) <= m.dcol
					ret = 1 # true, "collision"
					return ret
				end
			end
		end
	end
	return ret #false, "ongoing"
end

"""
mdp actions
"""
function actions(m::ActMdp, s::Array{Float64,1})::Vector{Float64}
	return m.actions
end

"""
mdp discount
"""
function discount(m::ActMdp)
	return m.gamma
end

"""
mdp pi0
"""
#function pi0(m::ActMdp, s)
#	return rand(actions(m, s)::Array)
#end

# use a generic function
pi0 = (m::ActMdp,s)-> rand(actions(m, s)::Array)::Float64

"""
mdp succProbReward
	returns sp, proba=T(s,a,sp), R(s,a,sp)
"""
function succProbReward(m::ActMdp, s::Array{Float64,1}, a::Float64)
	# UNUSED
	results = []
	return results
end

"""
mdp sampleSuccReward
"""
#function sampleSuccReward(m::ActMdp, s::Array{Float64,1}, a::Float64)::Tuple{Vector{Float64},Float64,Tuple{Bool, String}}
function sampleSuccReward(m::ActMdp, s::Array{Float64,1}, a::Float64)::Tuple{Vector{Float64},Float64,Int64}
	sp = m.Ts * s + m.Ta*a + m.Tt
	done = isEnd(m, sp)
	if sp[2] < 0 || 1==done #"collision" == done[2]
		r = -1.0
	elseif a <= -4
		r = -0.002
	else
		r = -0.001
	end
	#return Tuple(sp), r
	return sp, r, done
end

function succAndCost(m::ActMdp, s::Array{Float64, 1})::Array{Tuple{Float64, Vector{Float64}, Float64}}
	res = Tuple{Float64, Vector{Float64}, Float64}[] # (action, nextState, cost)
	for a in actions(m, s)
		sp, r, _ = sampleSuccReward(m, s, a)
		push!(res, (a, sp, -r))
	end
	return res
end


text:
using Documenter, Literate
using IndividualDisplacements

#download data dependencies if needed
IndividualDisplacements.get_ecco_velocity_if_needed();
IndividualDisplacements.get_occa_velocity_if_needed();

# generate tutorials and how-to guides using Literate
src = joinpath(@__DIR__, "src/")
lit = joinpath(@__DIR__, "../examples/")
notebooks = joinpath(src, "notebooks")

execute = true # Set to true for executing notebooks and documenter!
nb = true      # Set to true to generate the notebooks

lst1 = ["solid_body_rotation","random_flow_field","detailed_look","particle_cloud","global_ocean_circulation","three_dimensional_ocean"]
lst2 = ["solid_body_rotation","random_flow_field","detailed_look","particle_cloud","global_ocean_circulation","three_dimensional_ocean"]
#lst2 = ["none"]
tst1(x) = !isempty(lst1) && Bool(sum(isequal.(x, lst1)))
tst2(x) = !isempty(lst2) && Bool(sum(isequal.(x, lst2)))

for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue
	tst1(splitext(file)[1]) || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    Literate.markdown(ipath, opath, documenter = execute)
    nb && Literate.notebook(ipath, notebooks, execute = execute*tst2(splitext(file)[1]))
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) = [joinpath(folder, f) for f in readdir(joinpath(src, folder)) if ismd(f)]

p=pages("basics"); np=length(p)
i=findall((occursin).("solid_body_rotation",p));
j=findall((occursin).("random_flow_field",p));
p_tu=[p[j];p[i]]
i=findall((occursin).("detailed_look",p));
j=findall((occursin).("particle_cloud",p));
p_mi=[p[i];p[j]]

makedocs(
    sitename = "IndividualDisplacements",
    format = Documenter.HTML(),
    pages = [
		"Introduction" => "index.md",
        "User Guide" => "workflow.md",
		"Tool Box" => "API.md",
        "Example Guide" => "examples.md",
        "Tutorial Examples" => p_tu,
		"Real Ocean Cases" => pages("worldwide"),
        "MITgcm Examples" => p_mi],
    modules = [IndividualDisplacements]
)

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


text:
# ---
# title: 954. Array of Doubled Pairs
# id: problem954
# author: Jimmy Shen
# date: 2021-08-11
# difficulty: Medium
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/array-of-doubled-pairs/description/>
# hidden: true
# ---
#
# Given an array of integers `A` with even length, return `true` if and only if
# it is possible to reorder it such that `A[2 * i + 1] = 2 * A[2 * i]` for every
# `0 <= i < len(A) / 2`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [3,1,3,6]
#     Output: false
#
#
# **Example 2:**
#
#
#
#     Input: A = [2,1,2,6]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: A = [4,-2,2,-4]
#     Output: true
#     Explanation: We can take two groups, [-2,-4] and [2,4] to form [-2,-4,2,4] or [2,4,-2,-4].
#
#
# **Example 4:**
#
#
#
#     Input: A = [1,2,4,16,8,4]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `0 <= A.length <= 3 * 104`
#   * `A.length` is even.
#   * `-105 <= A[i] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
function can_reorder_doubled(arr::Vector{Int})::Bool
    cnt = counter(arr)
    for a in sort!(arr; by=abs)
        cnt[a] == 0 && continue
        cnt[2a] == 0 && return false
        cnt[2a] -= 1
        cnt[a] -= 1
    end
    return true
end
## @lc code=end


text:
# Sample script for stream tracing from a given starting point.
#
# Hongyang Zhou, hyzhou@umich.edu

# FieldTracer is not a registered pkg yet
# using Pkg
# Pkg.add(url="https://github.com/henry2004y/FieldTracer.jl.git", rev="master")
using PyPlot, FieldTracer, Vlasiator

Re = Vlasiator.Re

file = "bulk.0000501.vlsv"
nameρ = "rho"
nameV = "rho_v"

meta = load(file)

pcolormesh(meta, nameρ)

v = readvariable(meta, nameV)
vx = reshape(v[1,:], meta.ncells[1], meta.ncells[2])
vy = reshape(v[2,:], meta.ncells[1], meta.ncells[2])
# tracing starting point
xstart, ystart = 12Re, 0Re
# regular Cartesian mesh
x = range(meta.coordmin[1], meta.coordmax[1], length=meta.ncells[1])
y = range(meta.coordmin[2], meta.coordmax[2], length=meta.ncells[2])

# RK4 scheme by default
x1, y1 = trace2d(vx, vy, xstart, ystart, x, y;
   ds=0.5, maxstep=3000, gridType="ndgrid")
x1 ./= Re
y1 ./= Re

plot(x1, y1)

streamplot(meta, nameV, comp="xy", color="w", density=1.0)


text:
# This file was generated, do not modify it. # hide
y, X = unpack(caravan, ==(:Purchase), col->true)

mstd = machine(Standardizer(), X)
fit!(mstd)
Xs = MLJ.transform(mstd, X)

var(Xs[:,1]) |> r3

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# return the enclosure of a polynomial over a box
#=
function enclose(p::AbstractPolynomialLike, X::IntervalBox)

end
=#

max_stars_repo_name: JuliaReach/BernsteinExpansions

max_issues_repo_name: JuliaReach/BernsteinExpansions

max_forks_repo_name: JuliaReach/BernsteinExpansions


text:
using MANA
using JuliaDB, JWAS, DataFrames

function main()
	dir="~/MANA/example/"
        snpfilename=dir .* "plink.example"
        phenofilename=dir .* "phen.example2.txt"
        pedfilename=dir .* "ped.example2.csv"
        repeated_measures=false
        outputfilename=dir .* "out.example2.txt"
        Rvar=0.7 #total, not a proportion
        Gvar=0.3 #total, not a proportion
	Pvar=0
        y_tbl = loadtable(phenofilename; spacedelim = true)
        pedigree = get_pedigree(pedfilename)
        g_tbl, markerInfo = make_genotype_table(snpfilename)
        g_tbl = generate_genotypeclass_columns(g_tbl, markerInfo)
        tbl = join(y_tbl, g_tbl; lkey=:anml_key, rkey=:anml_key)
        tbl = MANA.center_phenotype(tbl)
        covariates = []
        if repeated_measures
                tbl = IndexedTables.transform(tbl, :anml_keyx => :anml_key)
                df = DataFrame(tbl)
                run_pedigree_model_multirecord(df, markerInfo, pedigree, Rvar, Gvar, Pvar; chainlength = 50000, outputfilename = outputfilename, covariates = covariates)
        else
                df = DataFrame(tbl)
                run_pedigree_model(df, markerInfo, pedigree, Rvar, Gvar; chainlength = 50000, outputfilename = outputfilename, covariates = covariates)
        end
end

main()

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
# https://www.youtube.com/watch?v=m2DvmHusyk4

import Pluto

session = Pluto.ServerSession()

path = tempname()
original = joinpath(pathof(Pluto) |> dirname |> dirname, "sample", "Tower of Hanoi.jl")
# so that we don't overwrite the file:
Pluto.readwrite(original, path)

@info "Loading notebook"
nb = Pluto.load_notebook(Pluto.tamepath(path));
session.notebooks[nb.notebook_id] = nb;

@info "Running notebook"
Pluto.update_save_run!(session, nb, nb.cells; run_async=false, prerender_text=true)

# nice! we ran the notebook, so we already precompiled a lot

@info "Starting HTTP server"
# next, we'll run the HTTP server which needs a bit of nasty code
t = @async Pluto.run(40404; session=session)

sleep(5)
download("http://localhost:40404/")

# this is async because it blocks for some reason
@async Base.throwto(t, InterruptException())
sleep(2) # i am pulling these numbers out of thin air

@info "Warmup done"


text:
# Management of a set of portfolios with overlapping stocks and constraints on different risk measures
# Choose how much to invest in each of a limited number of stocks, to maximize payoff

# source: https://github.com/JuliaOpt/Pajarito.jl/blob/master/examples/portfoliorisk/portfoliorisk.jl

mutable struct Portfolio
    id::Int
    size::Int
    stocks::Vector{String}
    returns::Dict{String,Float64}
    Sigmahalf::Dict{Tuple{String,String},Float64}
    risk::Symbol
    gamma::Float64
    Delta::Dict{Tuple{String,String},Float64}
end

function Base.show(io::IO, p::Portfolio)
    print(io, p.id)
end

#=========================================================
JuMP model functions
=========================================================#

# If using entropy constraints, requires a special branch of JuMP to allow modeling exponential cone (https://github.com/chriscoey/JuMP.jl/tree/coneconstr)
using JuMP

function portfoliorisk(solver, Portfolios, Stocks, Smax)
    m = Model(solver=solver)

    # Define investment fraction variables for stocks in each portfolio and limit sum to 1
    @variable(m, x[p in Portfolios, s in p.stocks] >= 0)
    @constraint(m, sum(x) <= 1)

    # Maximize total returns over all stocks
    @objective(m, Max, sum(p.returns[s]*x[p,s] for p in Portfolios, s in p.stocks))

    # Define stock choice indicator variables and limit number of chosen stocks to Smax and force zero investment if stock is not chosen
    @variable(m, y[s in Stocks], Bin)
    @constraint(m, sum(y) <= Smax)
    @constraint(m, [p in Portfolios, s in p.stocks], x[p,s] <= y[s])

    # Add risk constraints on Sigmahalf_p*x_p for each portfolio p
    @variable(m, Shx[p in Portfolios, s in p.stocks])
    @constraint(m, [p in Portfolios, s1 in p.stocks], Shx[p,s1] == sum(p.Sigmahalf[(s1,s2)]*x[p,s2] for s2 in p.stocks))

    for p in Portfolios
        if p.risk == :norm2
            @constraint(m, norm(Shx[p,s] for s in p.stocks) <= p.gamma)
        elseif p.risk == :robustnorm2
            lambda = @variable(m, lowerbound=0)
            @SDconstraint(m, [p.gamma [Shx[p,s] for s in p.stocks]' [x[p,s] for s in p.stocks]'; [Shx[p,s] for s in p.stocks] (diagm(fill(p.gamma, p.size)) - [lambda.*p.Delta[(s1,s2)] for s1 in p.stocks, s2 in p.stocks]) zeros(p.size, p.size); [x[p,s] for s in p.stocks] zeros(p.size, p.size) diagm(fill(lambda, p.size))] >= 0)
        # elseif p.risk == :entropy
        #     ent1 = @variable(m, [s in p.stocks])
        #     ent2 = @variable(m, [s in p.stocks])
        #     @constraint(m, sum(ent1[s] + ent2[s] for s in p.stocks) <= p.gamma^2)
        #     for s in p.stocks
        #         @Conicconstraint(m, [-ent1[s], 1 + Shx[p,s], 1] >= 0, :ExpPrimal)
        #         @Conicconstraint(m, [-ent2[s], 1 - Shx[p,s], 1] >= 0, :ExpPrimal)
        #     end
        else
            error("Invalid risk type $(p.risk)")
        end
    end

    return (m, x, y)
end

#=========================================================
Data generation functions
=========================================================#

# Generate model data from basic model options, reading portfolio data from datafiles
function generatedata(risks, counts, maxstocks, gammas, datadir, datafiles)
    N = sum(counts)
    @printf "\n\nGenerating data for %d portfolios\n" N
    @printf "\n%6s %6s %12s" "ID" "Size" "Risk type"

    Portfolios = Vector{Portfolio}(N)

    k = 0
    for b in 1:length(risks), bp in 1:counts[b] # For each risk type, each portfolio of the risk type
        k += 1

        pid = k
        pgamma = gammas[b]
        prisk = risks[b]

        (pstocks, preturns, pSigmahalf) = loadportfolio(datadir, datafiles[k], maxstocks[b]) # Read raw data for portfolio
        psize = length(pstocks)

        pDelta = Dict{Tuple{String,String},Float64}()
        if prisk == :robustnorm2
            # Generate random matrix and scale and clean zeros
            Deltahalf = randn(psize, psize)
            scalefactor = 1/10*norm([v for v in values(pSigmahalf)])/norm(vec(Deltahalf))
            @assert 1e-3 < scalefactor < 1e2
            for i in 1:psize, j in 1:psize
                val = scalefactor*Deltahalf[i,j]
                if val < 1e-3
                    Deltahalf[i,j] = 0.
                else
                    Deltahalf[i,j] = val
                end
            end

            # Fill Delta = Deltahalf*Deltahalf'
            for i in 1:psize, j in i:psize
                (s1, s2) = (pstocks[i], pstocks[j])
                pDelta[(s1,s2)] = pDelta[(s2,s1)] = vecdot(Deltahalf[i,:], Deltahalf[:,j])
            end
        end

        @printf "\n%6d %6d %12s" pid psize string(prisk)

        Portfolios[k] = Portfolio(pid, psize, pstocks, preturns, pSigmahalf, prisk, pgamma, pDelta)
    end

    stockset = Set{String}()
    for p in Portfolios
        union!(stockset, p.stocks)
    end
    Stocks = collect(stockset)

    return (Portfolios, Stocks)
end

# Load data from a .por file, returning returns, sqrt of covariance matrix, ticker names; take at most maxstocks stocks
function loadportfolio(datadir, datafile, maxstocks::Int)
    @assert endswith(datafile, ".por")
    file = open(joinpath(datadir, datafile), "r")

    n = parse(Int, chomp(readline(file)))

    if n > maxstocks
        taken = maxstocks
        takestocks = randperm(n)[1:taken]
    else
        taken = n
        takestocks = randperm(n)
    end

    rawreturns = split(chomp(readline(file)))
    @assert length(rawreturns) == n

    rawSigmahalf = zeros(n, n)
    for i in 1:n
        data = split(chomp(readline(file)))
        @assert length(data) == n
        for j in 1:n
            rawSigmahalf[i,j] = parse(Float64, data[j])
        end
    end

    @assert chomp(readline(file)) == ""

    line = chomp(readline(file))
    @assert startswith(line, "['") && endswith(line, "']")
    rawnames = split(line[3:end-2], "', '")
    @assert length(rawnames) == n

    pstocks = [String(rawnames[s]) for s in takestocks]

    preturns = Dict{String,Float64}(String(rawnames[s]) => parse(Float64, rawreturns[s]) for s in takestocks)

    pSigmahalf = Dict{Tuple{String,String},Float64}((String(rawnames[s1]), String(rawnames[s2])) => rawSigmahalf[s1,s2] for s1 in takestocks, s2 in takestocks)

    return (pstocks, preturns, pSigmahalf)
end

#=========================================================
Specify MICP solver
=========================================================#

using Pajarito

mip_solver_drives = false
log_level = 3
rel_gap = 1e-4

using CPLEX
mip_solver = CplexSolver(
    CPX_PARAM_SCRIND=(mip_solver_drives ? 1 : 0),
    # CPX_PARAM_SCRIND=1,
    CPX_PARAM_EPINT=1e-9,
    CPX_PARAM_EPRHS=1e-9,
    CPX_PARAM_EPGAP=(mip_solver_drives ? rel_gap : 0.)
)

# using SCS
# cont_solver = SCSSolver(eps=1e-3, max_iters=100000, verbose=0, warm_start=true)

# using ECOS
# cont_solver = ECOSSolver(verbose=false)

using Mosek
cont_solver = MosekSolver(LOG=0)

solver = PajaritoSolver(
    mip_solver_drives=mip_solver_drives,
    log_level=log_level,
    rel_gap=rel_gap,
    mip_solver=mip_solver,
    cont_solver=cont_solver,
    # solve_subp=false,
    # solve_relax=false,
    # init_sdp_soc=false,
    # sdp_soc=false,
    sdp_eig=false,
    # prim_cuts_only=true,
    # prim_cuts_always=true,
    # prim_cuts_assist=true,
    # dump_subproblems = true,
    # dump_basename = joinpath(pwd(), "portcbf/port")
)

#=========================================================
Specify model options and generate data
=========================================================#

srand(102)

risks = [:norm2, :entropy, :robustnorm2]
counts = [0, 0, 1]
maxstocks = [100, 20, 20]
gammas = [0.05, 0.05, 0.05]

datadir = joinpath(pwd(), "data")
datafiles = readdir(datadir)

(Portfolios, Stocks) = generatedata(risks, counts, maxstocks, gammas, datadir, datafiles)

Smax = round(Int, length(Stocks)/3)

@printf "\n\nChoose %d of %d unique stocks (sum of %d portfolio sizes is %d)\n\n" Smax length(Stocks) length(Portfolios) sum(p.size for p in Portfolios)

#=========================================================
Solve and print solution
=========================================================#

(m, x, y) = portfoliorisk(solver, Portfolios, Stocks, Smax)

# Save as CBF file
# name = "8,100_8,20_8,8"
# path = "/home/coey/portcbf/socexppsd/"
# ConicBenchmarkUtilities.jump_to_cbf(m, name, joinpath(path, "$name.cbf"))

status = solve(m)

@printf "\nStatus = %s\n" status

@printf "\nReturns (obj) = %8.4f\n" getobjectivevalue(m)
@printf "\nTotal number chosen = %d\n" sum(round(Int, getvalue(y[s])) for s in Stocks)
@printf "\nTotal fraction invested = %8.4f\n" getvalue(sum(x))

for p in Portfolios
    @printf "\nPortfolio %d investments\n" p.id

    for s in p.stocks
        if getvalue(y[s]) > 0.1
            @printf "%6s %8.4f\n" s getvalue(x[p,s])
        end
    end
end
println()


text:
using Plots, LaTeXStrings, Measures; pyplot()

f(x,y) = x^2 + y^2
f0(x) = f(x,0)
f2(x) = f(x,2)

xVals, yVals = -5:0.1:5 , -5:0.1:5
plot(xVals, [f0.(xVals), f2.(xVals)],
	c=[:blue :red], xlims=(-5,5), legend=:top,
	ylims=(-5,25), ylabel=L"f(x,\cdot)", label=[L"f(x,0)" L"f(x,2)"])
p1 = annotate!(0, -0.2, text("(0,0) The minimum\n of f(x,0)", :left, :top, 10))

z = [ f(x,y) for y in yVals, x in xVals ]
p2 = surface(xVals, yVals, z, c=cgrad([:blue, :red]),legend=:none,
	ylabel="y", zlabel=L"f(x,y)")

M = z[1:10,1:10]
p3 = heatmap(M, c=cgrad([:blue, :red]), yflip=true, ylabel="y",
	xticks=([1:10;], xVals), yticks=([1:10;], yVals))

plot(p1, p2, p3, layout=(1,3), size=(1200,400), xlabel="x", margin=5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using LinearAlgebra

const libblas = Base.libblas_name
const liblapack = Base.liblapack_name
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, chkstride1, Givens

function dggsvp3!(A::AbstractMatrix{Float64}, B::AbstractMatrix{Float64})
    chkstride1(A, B)
    jobu = 'U'
    jobv = 'V'
    jobq = 'Q'
    m, n = size(A)
    p = size(B,1)
    # if BlasInt(size(B, 2)) != n
    #     throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs $n"))
    # end
    tola = max(m, n) * opnorm(A, 1) * eps(Float64)
    tolb = max(p, n) * opnorm(B, 1) * eps(Float64)
    k = Ref{BlasInt}()
    l = Ref{BlasInt}()
    ldu = max(1, m)
    U = similar(A, Float64, ldu, m)
    ldv = max(1, p)
    V = similar(A, Float64, ldv, p)
    ldq = max(1, n)
    Q = similar(A, Float64, ldq, n)
    iwork = Vector{BlasInt}(undef, n)
    tau = Vector{Float64}(undef, n)
    work = Vector{Float64}(undef, 1)
    lwork = BlasInt(-1)
    info  = Ref{BlasInt}()

    for i in 1:2  # first call returns lwork as work[1]
        ccall((@blasfunc(dggsvp3_), liblapack), Cvoid,
        (Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
        Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
        Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
        Ref{Float64}, Ref{Float64}, Ref{BlasInt}, Ref{BlasInt},
        Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
        Ptr{Float64}, Ref{BlasInt},
        Ptr{BlasInt}, Ptr{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt}
        ),
        jobu, jobv, jobq,
        m, p, n, A,
        max(1, stride(A, 2)), B, max(1, stride(B, 2)),
        tola, tolb, k, l,
        U, ldu, V, ldv,
        Q, ldq,
        iwork, tau, work, lwork, info)
        # chklapackerror(info[])
        if i == 1
            lwork = BlasInt(real(work[1]))
            resize!(work, lwork)
        end
    end
    return U, V, Q, k[], l[], A, B
end


text:
import Compat: @static, is_windows

"""
    extract(archive, directory, [files]) -> Void

Extracts files from a compressed tar `archive` to the specified `directory`. If `files` is
specified only the files given will be extracted.
"""
function extract(archive, directory, files=String[])
    @static if is_windows()
        run(pipeline(`7z x $archive -y -so`, `7z x -si -y -ttar -o$directory $files`))
    else
        run(`tar xvf $archive --directory=$directory $files`)
    end
end

"""
    @optional(expr)

Creates multiple method signatures to allow optional arguments before required arguments.
For example:

    f(a=1, b=2, c) = ...

becomes:

    f(a, b, c) = ...
    f(a, c) = f(a, 2, c)
    f(c) = f(1, 2, c)
"""
macro optional(ex)
    esc(Expr(:block, optional(ex)...))
end

function optional(ex::Expr)
    funcs = Expr[]
    if ex.head == :function || ex.head == :(=)
        sig, body = ex.args
    else
        throw(ArgumentError("Expected expression to be a function"))
    end

    # Determine the optional arguments that occur prior to the last required argument
    dynamic = Int[]
    last_required = 0
    for (i, arg) in enumerate(sig.args)
        if isa(arg, Expr) && arg.head == :kw
            push!(dynamic, i)
        end

        if isa(arg, Symbol) || arg.head == :(::)
            last_required = i
        end
    end
    dynamic = filter(i -> i < last_required, dynamic)

    # Return early if we don't have to generate additional methods
    isempty(dynamic) && return Expr[ex]

    # Generate a new signature which converts the dynamic arguments to required arguments
    default = Array{Any}(length(dynamic))
    new_sig = Array{Any}(length(sig.args))
    for (i, arg) in enumerate(sig.args)
        if i in dynamic
            name, value = arg.args
            new_sig[i] = name
            default[findin(dynamic,i)] = value
        else
            new_sig[i] = arg
        end
    end

    # Generate the primary function
    push!(funcs, Expr(ex.head, Expr(:call, new_sig...), body))

    # Generate the
    func_call = copy(sig)
    for (i, argument) in enumerate(func_call.args)
        isa(argument, Expr) || continue

        if argument.head == :parameters
            for keyword in argument.args
                keyword.args[2] = keyword.args[1]
            end
        elseif argument.head == :(::)
            func_call.args[i] = argument.args[1]
        elseif argument.head == :kw
            if isa(argument.args[1], Symbol)
                func_call.args[i] = argument.args[1]
            else
                func_call.args[i] = argument.args[1].args[1]
            end
        end
    end

    # Generate additional methods which call the primary function
    while !isempty(dynamic)
        i = pop!(dynamic)

        deleteat!(new_sig, i)
        func_call.args[i] = pop!(default)

        push!(funcs, Expr(:(=), Expr(:call, new_sig...), Expr(:block, copy(func_call))))
    end

    return funcs
end


text:
using Revise
# using Plots
using Pkg
# Pkg.activate(".")
using BenchmarkTools
using RadiativeTransfer
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM

# Load parameters from file
# parameters = vSmartMOM.parameters_from_yaml("RadiativeTransfer/src/vSmartMOM/ModelParameters/DefaultParameters.yaml")

# params_dict = YAML.load_file("RadiativeTransfer/src/vSmartMOM/ModelParameters/DefaultParameters.yaml")

# vSmartMOM.validate_yaml_parameters(params_dict)

# file = "RadiativeTransfer/src/vSmartMOM/ModelParameters/SampleProfile.yaml"
# profile_nc4 = vSmartMOM.read_atmos_profile(parameters.file, parameters.lat, parameters.lon, parameters.timeIndex);
# profile_yaml = vSmartMOM.read_atmos_profile(file);

# Sets all the "specific" parameters
parameters = vSmartMOM.default_parameters();

# Generates all the derived attributes from above parameters
model = model_from_parameters(parameters);
#model.params.architecture = RadiativeTransfer.Architectures.GPU();

R = vSmartMOM.rt_run(model)


text:
# TODO: Connect with GitLink config
const _LOCK_FORCE_TIME = 360.0
const _LOOP_FREC_FAIL_PENALTY = 5.0
const _LOOP_FREC_IDLE_PENALTY = 0.5

"""
    sync_link(gl::GitLink; verbose = true, force = false, tries = 1)

Try to pull/merge state/push the `GitLink`.
It is a lazy method, if no action is require no action will be made (use `force` to avoid it).
This method will sleep till (or timeout `tout`) the GitLink lock is free (which must by must of the time, but...).
The method will attempt `tries` tries till success.
Returns `true` if the action was successful.
"""
function sync_link(gl::GitLink;
        verbose = true, force = true,
        before_push::Function = () -> nothing,
        tout = _LOCK_FORCE_TIME,
        tries = 1
    )

    for t in 1:max(tries, 0)

        # Globals
        lf = lock_file(gl)
        lid = ""
        pull_token = _get_pull_token(gl)
        push_token = _get_push_token(gl)
        stage_token = _get_stage_token(gl)

        try

            ## ---------------------------------------------------
            # ACQUIRE LOCK
            unlock(lf, lid) # In case it is mine
            lid, ttag = acquire_lock(lf; tout, force = true)

            ## ---------------------------------------------------
            verbose && @info("Lock acquired", lid, ttag)

            ## ---------------------------------------------------
            # RESOLVE ACTION FlAG
            pull_flag = force || is_pull_required(gl)
            push_flag = force || is_push_required(gl)
            doloop = pull_flag || push_flag
            if !doloop # Handle idle
                add_loop_frec!(gl, _LOOP_FREC_IDLE_PENALTY)
                continue
            end

            ## ---------------------------------------------------
            verbose && @info("Doing", pull_flag, push_flag)

            ## ---------------------------------------------------
            # HARD PULL (Loop)
            pull_ok = hard_pull(gl; verbose, clearwd = true)
            if !pull_ok # Handle fail
                add_loop_frec!(gl, _LOOP_FREC_FAIL_PENALTY)
                continue
            end
            pull_token = _set_pull_token(gl) # Aknowlage successful pull

            ## ---------------------------------------------------
            verbose && @info("Pull info", pull_token, chash = _HEAD_hash(gl))

            ## ---------------------------------------------------
            if push_flag

                ## ---------------------------------------------------
                # MERGE STAGE
                _merge_stage(gl)

                ## ---------------------------------------------------
                # CALLBACK
                before_push()

                ## ---------------------------------------------------
                # SOFT PUSH
                push_ok = soft_push(gl; verbose)
                if !push_ok # Handle fail
                    add_loop_frec!(gl, _LOOP_FREC_FAIL_PENALTY)
                    continue
                end
                push_token = _set_push_token(gl) # Aknowlage successful pull

                ## ---------------------------------------------------
                verbose && @info("Push info", push_token, chash = _HEAD_hash(gl))

                ## ---------------------------------------------------
                stage_token = _set_stage_token(gl, push_token) # Aknowlage successful upload

            end

            ## ---------------------------------------------------
            # HANDLE SUCCEFUL LOOP
            loop_frec!(gl, _MIN_LOOP_FREC) # Reset loop frec
            is_stage_sync = _is_stage_up_to_day(gl)

            ## ---------------------------------------------------
            verbose && @info("Success!", stage_token, is_stage_sync)
            verbose && println()

            return true

        finally
            # Free lock
            unlock(lf, lid)
        end

    end #tries

    return false

end


text:
# GMM estimation for a sample from Chi^2(theta)
# compare to two method of moments estimators (see chi2mm.m)
using Distributions, Plots
function main()
n = 30
theta = 3.0
reps = 1000
results = zeros(reps,2)
for i = 1:1000
    y = rand(Chisq(theta), n)
    moments = theta -> [theta.-y  theta.-0.5.*(y .- mean(y)).^2.0]
    thetahat, junk, junk, ms, junk = gmm(moments, [3.0], eye(2))
    results[i,1] = sqrt(n)*(thetahat[1]-theta)
    W = inv(cov(ms))
    thetahat, junk, ms, junk = gmm(moments, thetahat, W)
    results[i,2] = sqrt(n)*(thetahat[1]-theta)
end
p1 = npdensity(results[:,1])
plot!(p1, title="Inefficient")
p2 = npdensity(results[:,2])
plot!(p2, title="Efficient")
plot(p1,p2,layout=(2,1),xlims=(-10,20))
savefig("Efficient.svg")
gui()
println("Monte Carlo covariance: (1,1) is inefficient GMM, (2,2) is efficient")

prettyprint(cov(results))
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# Estimates 5 separate models for the Nerlove Cobb-Douglas model,
# and does a Chow test for pooled coefficients
using Econometrics, Plots, DelimitedFiles, LinearAlgebra
function main()
data = readdlm("nerlove.data")
data = data[:,2:6]
data = log.(data)
n = size(data,1)
y = data[:,1]
x = data[:,2:end]
x = [ones(n,1) x]
#b, junk, junk, junk = ols(y, x, names=names)
k = size(x,2)

# create the block diagonal X matrix corresponding to separate coefficients
big_x = zeros(n,5*k)
for i=1:k
	startrow = (i-1)*29+1
	endrow = i*29
	startcol =(i-1)*k + 1
	endcol = i*k
	big_x[startrow:endrow,startcol:endcol] += x[startrow:endrow,:]
end
x = big_x

names = ["constant", "output","labor", "fuel", "capital"]
names = [names; names; names; names; names] # copy 5 times

println("Nerlove model: 5 separate regressions for different output levels")
b, junk, junk, junk = ols(y, x, names=names)
output = 1:5
output = output .*5 .+2 .- 5
output = b[output,:]
rts = 1 ./ output

# gset term X11
group = 1:5
group = group
plot(group, rts, label = "", yaxis="RTS", xaxis = "Output Group", show=true)
#savefig("rts.png")

# Chow test
R = Matrix{Float64}(I, 5, 5)
Z = zeros(5,5)
R = [
	R -R Z Z Z;
	R Z -R Z Z;
	R Z Z -R Z;
	R Z Z Z -R]
r = zeros(20,1)

println("Chow test: note that the restricted model")
println("gives the same results as the original model")
ols(y, x, R=R, r=r, names=names)
TestStatistics(y, x, R, r)
return
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics