text:
# ## [Tracers](@id tracer-model)
#
#md # !!! note
#md #
#md #     Usage: Enable tracers using a keyword argument in the AtmosModel specification\
#md #     `tracers = NoTracer()`\
#md #     `tracers = NTracers{N, FT}(δ_χ)` where N is the number of tracers required.\
#md #     FT is the float-type and $\delta_{\chi}$ is an SVector of diffusivity scaling coefficients
#
# In `tracers.jl`, we define the equation sets governing
# tracer dynamics. Specifically, we address the the equations
# of tracer motion in conservation form,
#
#
using DocStringExtensions
#
# ### [Equations](@id tracer-eqns)
# ```math
# \frac{\partial \rho\chi}{\partial t} +  \nabla \cdot ( \rho\chi u) = \nabla \cdot (-\rho\delta_{D\chi}\mathrm{D_{T}}\nabla\chi) + \rho \mathrm{S}
# ```
# where  $$\chi$$ represents the tracer species, $$\mathrm{S}$$ represents the tracer source terms and $$\delta_{D\chi} \mathrm{D_{T}}$$ represents the scaled turbulent eddy diffusivity for each tracer.
# Currently a default scaling of `1` is supported.
# The equation as written above corresponds to a single scalar tracer, but can be extended to include
# multiple independent tracer species.
#
# We first define an abstract tracer type, and define the
# default function signatures. Two options are currently
# supported. [`NoTracers`](@ref no-tracers),
# and [`NTracers`](@ref multiple-tracers).
#
# ### [Abstract Tracer Type](@id abstract-tracer-type)
#
# Default methods for a generic tracer type are defined here.
#

abstract type TracerModel end

export TracerModel, NoTracers, NTracers

vars_state(::TracerModel, ::AbstractStateType, FT) = @vars()

function atmos_init_aux!(
    ::TracerModel,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    nothing
end
function atmos_nodal_update_auxiliary_state!(
    ::TracerModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
function compute_gradient_flux!(
    ::TracerModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
function compute_gradient_argument!(
    ::TracerModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
function wavespeed_tracers!(
    ::TracerModel,
    wavespeed::Vars,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

# ### [NoTracers](@id no-tracers)
# The default tracer type in both the LES and GCM configurations is the
# no tracer model. (This means no state variables for tracers are being
# carried around). For the purposes of this model, moist variables are
# considered separately in `moisture.jl`.
#
"""
    NoTracers <: TracerModel
No tracers. Default model.
"""
struct NoTracers <: TracerModel end

# ### [NTracers](@id multiple-tracers)
# Allows users to specify an integer corresponding to the number of
# tracers required. This can be extended to provide a vector
# corresponding to a scaling factor for each individual tracer
# component. Note that tracer naming is not currently supported,
# i.e. the user must track each tracer variable based on its
# numerical index. Sources can be added to each tracer based on the
# same numerical index. Initial profiles must be specified using the
# `init_state_prognostic!` hook at the experiment level.

"""
    NTracers{N, FT} <: TracerModel
Currently the simplest way to get n-tracers in an AtmosModel run
using the existing machinery. Model input: SVector of
diffusivity scaling coefficients. Length of SVector allows number
of tracers to be inferred. Tracers are currently identified by indices.

# Fields
#
$(DocStringExtensions.FIELDS)
"""
struct NTracers{N, FT} <: TracerModel
    "N-component `SVector` with scaling ratios for tracer diffusivities"
    δ_χ::SVector{N, FT}
end

vars_state(tr::NTracers, ::Prognostic, FT) = @vars(ρχ::typeof(tr.δ_χ))
vars_state(tr::NTracers, ::Gradient, FT) = @vars(χ::typeof(tr.δ_χ))
vars_state(tr::NTracers, ::GradientFlux, FT) =
    @vars(∇χ::SMatrix{3, length(tr.δ_χ), FT, 3 * length(tr.δ_χ)})
vars_state(tr::NTracers, ::Auxiliary, FT) = @vars(δ_χ::typeof(tr.δ_χ))

function atmos_init_aux!(
    tr::NTracers,
    am::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.tracers.δ_χ = tr.δ_χ
end
function compute_gradient_argument!(
    tr::NTracers,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.tracers.χ = state.tracers.ρχ * ρinv
end
function compute_gradient_flux!(
    tr::NTracers,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.tracers.∇χ = ∇transform.tracers.χ
end

function wavespeed_tracers!(
    tr::NTracers{N, FT},
    wavespeed::Vars,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
) where {N, FT}

    ρinv = 1 / state.ρ
    u = ρinv * state.ρu
    uN = abs(dot(nM, u))

    wavespeed.tracers.ρχ = fill(uN, SVector{N, FT})

    return nothing
end

#####
##### Tendency specification
#####

eq_tends(pv::Tracers{N}, m::NTracers{N}, tt::Flux{SecondOrder}) where {N} =
    (Diffusion(),)


text:
# +
module PetriNets

abstract type AbstractNode end

const Option{T} = Union{T, Missing}

mutable struct Node <: AbstractNode
    name::Any
    value::Int32
end

struct Transition <: AbstractNode
    name::Any
    inputs::Vector{Tuple{Node,Bool}}
    outputs::Vector{Node}
end

struct PetriNet
    nodes::Vector{Node}
    transitions::Vector{Transition}
end

struct PetriError <: Exception
    msg::AbstractString
end

name(n::AbstractNode) = n.name

"""    ready(t::Tuple{Node,Bool})

check that the node is ready, the bool is used to indicate inhibitor arc
"""
function ready(t::Tuple{Node,Bool})
    n = t[1]
    b = t[2]
    return b ? value(n) > 0 : value(n) == 0
end

function ready(n::Transition)
    return length(filter(x->ready(x) == false, n.inputs)) == 0
end

function value(n::Node)
    return n.value
end

clear(n::Node) = begin
    n.value = 0
    return
end

function run!(n::Transition)
    #@show name(n)
    if ready(n)
        if length(findall(i->!xor(i[1] in n.outputs, i[2]), n.inputs)) == length(n.inputs)
            throw(PetriError("Infinite loop detected at transition $(name(n))"))
        end
        #@info "Decrease inputs"
        map(n.inputs) do i
            # if you know that i is ready, and it is empty, then it must be an inhibitor.
            i[1].value -= i[1].value > 0 ? 1 : 0
        end
        #@info "Increase outputs"
        map(n.outputs) do i
            i.value += 1
        end
    end
end

function run!(p::PetriNet)
    run = true
    while run
        run = false
        for n in p.transitions
            while ready(n)
                run = true
                run!(n)
            end
        end
    end
end

end

# +
n1 = PetriNets.Node(:n1, 10)
n2 = PetriNets.Node(:n2, 1)
n3 = PetriNets.Node(:n3, 0)
n4 = PetriNets.Node(:n4, 0)
n5 = PetriNets.Node(:n5, 0)
n6 = PetriNets.Node(:n6, 0)
n7 = PetriNets.Node(:n7, 1)
n8 = PetriNets.Node(:n8, 0)

t1 = PetriNets.Transition(:t1, [(n1,true)], [n3])
t2 = PetriNets.Transition(:t2, [(n2,true)], [n3])
t3 = PetriNets.Transition(:t3, [(n3,true)], [n4,n5])
t4 = PetriNets.Transition(:t4, [(n4,true)], [n6])
t5 = PetriNets.Transition(:t5, [(n5,true)], [n6])
t6 = PetriNets.Transition(:t6, [(n6,true)], [n8])
t7 = PetriNets.Transition(:t7, [(n7,true)], [n8])

pn1 = PetriNets.PetriNet([n1,n2,n3,n4,n5,n6,n7,n8], [t1, t2, t3, t4, t5, t6, t7])
PetriNets.run!(pn1)
@show n8
@assert pn1.nodes[end].value == 23 "pn1 did not compute correctly"

n9 = PetriNets.Node(:n9, 1)
n10 = PetriNets.Node(:n10, 0)
n11 = PetriNets.Node(:n11, 0)

t8 = PetriNets.Transition(:t8, [(n9, true), (n10, false)], [n11])

pn2 = PetriNets.PetriNet([n9, n10, n11], [t8])
PetriNets.run!(pn2)
@show n11

n12 = PetriNets.Node(:n12, 1)
n13 = PetriNets.Node(:n13, 0)
n14 = PetriNets.Node(:n14, 0)
t9 = PetriNets.Transition(:t9, [(n12, true), (n13, false)], [n12,n14])
pn3 = PetriNets.PetriNet([n12, n13], [t9])
try
    PetriNets.run!(pn3)
catch PetriError
    @warn "pn3 has an infinite loop, this warning is benign"
end

# -

import TikzPictures
using Catlab.Graphics
using Catlab.Doctrines
using Catlab.WiringDiagrams
import Catlab.WiringDiagrams.WiringDiagramCore: boxes, WiringDiagram

A, B, C, D, One = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C, :D, :1)

function boxes(p::PetriNets.PetriNet, mapping)
    j = 1
    map(p.nodes) do n
        @info "Constructing Homsets to $(name(n))"
        map(n.inputs) do i
            s,t = name(i), name(n)
            x,y = mapping[s], mapping[t]
            j += 1
            WiringDiagram(Hom(Symbol("f_$j"), x,y))
        end
    end
end

composel(x::Vector{Vector{WiringDiagram}}) = foldl.(compose, x)
composel(x::Vector{WiringDiagram}) = foldl(compose, x)

WiringDiagram(p::PetriNets.PetriNet, mappi lng) = composel(boxes(p, mapping)) |> composel

mapping = Dict(:A=>A,
              :B=>B,
              :C=>C,
              :D=>D,
              :C1=>One)
bx = boxes(p, mapping)

@assert length(bx) == PetriNets.nv(p)
d = WiringDiagram(p, mapping)

struct PetriEdge
    value
    input
    output
    meta
end

function PetriNet(d::WiringDiagram, mapping)
    edges = map(boxes(d)) do box
        PetriEdge(box.value, box.input_ports[1], box.output_ports[1],())
    end
end
# can we represent a petrinet as a presentation?
# This is like a database schema where you have Ob representing tables
# and Homs representing foreign keys. The equations are the business rules
# preserved by the DB.
# How do we represent a specific petri net?
# As a Petri-Instance, ie. a functor from Petri to what?
@present Petri(FreeCategory) begin
  # Primitive concepts.
  Node::Ob
  Symb::Ob
  Input::Ob
  Position::Ob

  name::Hom(Node, Symb)
  output::Hom(Node, Node)
  slot::Hom(Input, Position)
  source::Hom(Input, Node) # the place a connection came from
  target::Hom(Input, Node) # the place a connection goes to

  # Defined concepts.
  # second_level_manager := compose(manager, manager)
  # third_level_manager := compose(manager, manager, manager)

  # Abbreviations (no syntactic term for LHS).
  # boss = manager

  # Managers work in the same department as their employees.
  # compose(boss, works_in) == works_in
  # The secretary of a department works in that department.
  # compose(secretary, works_in) == id(Department)
end

struct Instance{C, D}
    index::C
    image::D
end

inst = Instance(Petri, p)

import Base: getindex, setindex
function getindex(inst::Instance{Presentation{Symbol},PetriNets.PetriNet}}, obj::Catlab.Doctrines.FreeCategory.Ob{:generator})
    instance.index
end

ext: jl


text:
using PyPlot
using Convex
#using CenterOfMass
include("CenterOfMass.jl")
include("setup_parameters.jl")

plt[:close]("all")

mass = min_mass + rand() * (max_mass - min_mass)
com_p = r_attach * rand_in_circle()

prior = initialize_prior(circle_ps, resolution, interior_q, masses)

critical_forces_by_point = map(attachment_ws) do wrench
  get_critical_values(circle_ws, wrench, prior, interior_q)
end

normals_matrices = map(critical_forces_by_point) do forces
  normal_matrix(forces[:], sigma^2)
end

pygui(false)
for ii = 1:16
  # selection of second point
  n_attach = length(critical_forces_by_point)
  csqmis = map(1:n_attach) do ii
    println("computing csqmi")
    compute_mutual_information(get_data(prior)[:], normals_matrices[ii])
  end

  #fig = figure(figsize=(6,6), dpi=600)

  # application at second point
  applied_p = attachment_ps[indmax(csqmis)]
  critical_forces = critical_forces_by_point[indmax(csqmis)]

  #@show critical_forces
  #fig, ax = plt[:subplots](1)
  #plot_field(reshape(critical_forces[1,:,:], (size(critical_forces, 2),
  #size(critical_forces, 3))))

  applied_f = critical_force_from_points(circle_ps, com_p, applied_p, mass)

  f_hat = applied_f + sigma * randn()
  update_prior!(prior, critical_forces, f_hat, sigma)

  cloud = to_cloud(prior, interior_q)

  fig, ax = plt[:subplots](1)

  plot_attachment_csqmis(circle_ps, attachment_ps, csqmis)
  plot_solution(applied_p)

  scaled_ps = (15*cloud[4,:]'/maximum(cloud[4,:]))
  scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', "z", scaled_ps.^2, "purple", alpha = 0.5)

  # plot point
  scatter3D([com_p[1]], [com_p[2]], [mass], "z", 200, "red", marker="*", alpha=1)
  #scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', cloud[4,:]')

  circle_polygon = hcat(circle_ps..., circle_ps[1])
  plot3D(circle_polygon[1,:][:], circle_polygon[2,:][:], circle_polygon[3,:][:],
    color = "k", alpha = 1, linestyle = "solid")

  #title("probability distribution of center-of-mass")
  xlabel("X")
  ylabel("Y")
  zlabel("M")
  fig[:axes][1][:get_yaxis]()[:set_visible](false)
  fig[:axes][1][:get_xaxis]()[:set_visible](false)
  xlim(-1.0,1.0)
  ylim(-1.0,1.0)
  #axis("off")
  savefig("fig/belief_3d$(ii).png", pad_inches=0.01, bbox_inches="tight")
end
pygui(true)

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
############################################################################
### Types
############################################################################

struct SemDiffNLopt{A, A2, B, B2, C} <: SemDiff
    algorithm::A
    local_algorithm::A2
    options::B
    local_options::B2
    equality_constraints::C
    inequality_constraints::C
end

Base.@kwdef mutable struct NLoptConstraint
    f
    tol = 0.0
end

############################################################################
### Constructor
############################################################################

function SemDiffNLopt(;
        algorithm = :LD_LBFGS,
        local_algorithm = nothing,
        options = Dict{Symbol, Any}(),
        local_options = Dict{Symbol, Any}(),
        equality_constraints = Vector{NLoptConstraint}(),
        inequality_constraints = Vector{NLoptConstraint}(),
        kwargs...)
    applicable(iterate, equality_constraints) || (equality_constraints = [equality_constraints])
    applicable(iterate, inequality_constraints) || (inequality_constraints = [inequality_constraints])
    return SemDiffNLopt(algorithm, local_algorithm, options, local_options, equality_constraints, inequality_constraints)
end

############################################################################
### Recommended methods
############################################################################

update_observed(diff::SemDiffNLopt, observed::SemObs; kwargs...) = diff

############################################################################
### additional methods
############################################################################

algorithm(diff::SemDiffNLopt) = diff.algorithm
local_algorithm(diff::SemDiffNLopt) = diff.local_algorithm
options(diff::SemDiffNLopt) = diff.options
local_options(diff::SemDiffNLopt) = diff.local_options
equality_constraints(diff::SemDiffNLopt) = diff.equality_constraints
inequality_constraints(diff::SemDiffNLopt) = diff.inequality_constraints

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SemDiffNLopt)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end


text:
# Add fallbacks for predict_* which require mean, mode, median.

const BadMeanTypes   = Union{AbstractArray{<:Finite},Table(Finite)}
const BadMedianTypes = Union{AbstractArray{<:Finite},Table(Finite)}

const err_wrong_target_scitype(actual_scitype) = ArgumentError(
    "Attempting to compute mode of predictions made "*
    "by a model with `$actual_scitype` targets. ")

# mode:
predict_mode(m, fitresult, Xnew) =
    mode.(predict(m, fitresult, Xnew))

# mean:
predict_mean(m, fitresult, Xnew) =
    predict_mean(m, fitresult, Xnew, target_scitype(m))
predict_mean(m, fitresult, Xnew, ::Any) =
    mean.(predict(m, fitresult, Xnew))
predict_mean(m, fitresult, Xnew, ::Type{<:BadMeanTypes}) =
    throw(err_wrong_target_scitype(Finite))

# median:
predict_median(m, fitresult, Xnew) =
    predict_median(m, fitresult, Xnew, target_scitype(m))
predict_median(m, fitresult, Xnew, ::Any) =
    median.(predict(m, fitresult, Xnew))
predict_median(m, fitresult, Xnew, ::Type{<:BadMedianTypes}) =
    throw(err_wrong_target_scitype(Finite))

# not in MLJModelInterface as methodswith requires InteractiveUtils
MLJModelInterface.implemented_methods(::FI, M::Type{<:MLJType}) =
    getfield.(methodswith(M), :name) |> unique

# serialization fallbacks:
MLJModelInterface.save(filename, model, fitresult; kwargs...) = fitresult
MLJModelInterface.restore(filename, model, serializable_fitresult) =
                          serializable_fitresult

# to suppress inclusion of abstract types in the model registry.
for T in (:Supervised, :Unsupervised,
          :Interval, :Static, :Deterministic, :Probabilistic)
    ex = quote
        MLJModelInterface.is_wrapper(::Type{$T}) = true
    end
    eval(ex)
end


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
##------------------------  ------------------------##
using Ipopt
using LinearAlgebra, ForwardDiff, Zygote
##------------------------  ------------------------##
# Implement hs071 from IPOPT/IPOPT.jl documentation using
# analytical and automatic differentiation.

# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

# solution/starting point = x
# objective function = f(x)
# constraints = g(x)

# number of variables
n = 4

# bounds
# 1 <= x1, x2, x3, x4 <= 5
# lower bounds
x_L = [1.0, 1.0, 1.0, 1.0]
# upper bounds
x_U = [5.0, 5.0, 5.0, 5.0]

# number of constraints, dim of g(x)
m = 2

# constraints lower bounds
# x1 * x2 * x3 * x4 >= 25
# x1^2 + x2^2 + x3^2 + x4^2 = 40
g_L = [25.0, 40.0]
# constraints upper bounds
g_U = [2.0e19, 40.0]# unbounded first term

# initial point
x_0 = [1.0, 5.0, 5.0, 1.0]

##---------------------- obj/constrain functions ----------------------##
# objective function
function eval_f(x)
    return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]
end

# constraint function
function eval_g(x, g)
    g[1] = x[1] * x[2] * x[3] * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
    return g
end

function eval_g_ad(x, g)
    g = eltype(x).(g)
    g[1] = x[1] * x[2] * x[3] * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
    return g
end

##------------------------ grad obj function ------------------------##

# grad of obj function ∇f(x)
grad_f = zeros(n)
# analytical grad of obj function
function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
    grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
    grad_f[2] = x[1] * x[4]
    grad_f[3] = x[1] * x[4] + 1
    return grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

# now with automatic differentiation
∇f = x -> ForwardDiff.gradient(eval_f, x)

# check both are the same
eval_grad_f(x_0, grad_f)
@show grad_f
@show ∇f(x_0)

function eval_grad_f_ad(x::Vector{Float64}, grad_f::Vector{Float64})
    grad_f[:] = ForwardDiff.gradient(eval_f, x)[:]
end
@show eval_grad_f_ad(x_0, grad_f)

##--------------------- jac of constraints ∇g(x)^T --------------------##
println("jacobian of constraints")
g = zeros(m)
∇g = x -> ForwardDiff.jacobian(x->eval_g_ad(x, g), x)
jac_temp = ∇g(x_0)
sparsity = jac_temp .!= 0.0
I = ((i, j) for i = 1:m, j = 1:n if sparsity[i,j] != 0.0)
rows = Int32.(getindex.(I, 1)[:])
cols = Int32.(getindex.(I, 2)[:])
values = zeros(length(rows))

# analytical version
function eval_jac_g(
    x::Vector{Float64},
    rows::Vector{Int32},
    cols::Vector{Int32},
    values::Union{Nothing,Vector{Float64}},
)
    if values === nothing
    # return the structure of the jacobian,
    # Sparse Triplet Matrix Format
    # i.e. the rows and cols of the non-zero entries
    # index style=FORTRAN STYLE (1-indexed)
        # Constraint (row) 1
        rows[1] = 1
        cols[1] = 1
        rows[2] = 1
        cols[2] = 2
        rows[3] = 1
        cols[3] = 3
        rows[4] = 1
        cols[4] = 4
        # Constraint (row) 2
        rows[5] = 2
        cols[5] = 1
        rows[6] = 2
        cols[6] = 2
        rows[7] = 2
        cols[7] = 3
        rows[8] = 2
        cols[8] = 4
    else
    # return the values of the jacobian of the constraints
        # Constraint (row) 1
        values[1] = x[2] * x[3] * x[4]  # 1,1
        values[2] = x[1] * x[3] * x[4]  # 1,2
        values[3] = x[1] * x[2] * x[4]  # 1,3
        values[4] = x[1] * x[2] * x[3]  # 1,4
        # Constraint (row) 2
        values[5] = 2 * x[1]  # 2,1
        values[6] = 2 * x[2]  # 2,2
        values[7] = 2 * x[3]  # 2,3
        values[8] = 2 * x[4]  # 2,4
    end
    @show values
    return nothing
end
println("analytical")
eval_jac_g(x_0, rows, cols, nothing)
@show rows; @show cols
eval_jac_g(x_0, rows, cols, values)
@show values
display(reshape(values, m, n))

println("\nAD")

function jac_structure!(
    m, n,
    rows::AbstractVector{<:Integer},
    cols::AbstractVector{<:Integer},
  )
    I = ((i, j) for i = 1:m, j = 1:n)
    rows .= getindex.(I, 1)[:]
    cols .= getindex.(I, 2)[:]
    return rows, cols
end

function jac_structure_sparse!(
    x, m, n,
    rows::AbstractVector{<:Integer},
    cols::AbstractVector{<:Integer},
  )
    jac_temp =  ForwardDiff.jacobian(x->eval_g_ad(x, g), x)
    sparsity = jac_temp .!= 0.0
    I = ((i, j) for i = 1:m, j = 1:n if sparsity[i,j] != 0.0)
    rows .= getindex.(I, 1)[:]
    cols .= getindex.(I, 2)[:]
    return rows, cols
end

∇g = x -> ForwardDiff.jacobian(x->eval_g_ad(x, g), x)

function eval_jac_g_ad(x, rows, cols, values)

    if values === nothing
        jac_structure_sparse!(x_0, m, n, rows, cols)
        # jac_structure!(m, n, rows, cols)
    else

        values[:] = ForwardDiff.jacobian(x->eval_g_ad(x, g), x)'[:]
    end
    return nothing
end

#jac_structure!(m, n, rows, cols)
@show jac_structure_sparse!(x_0, m, n, rows, cols)
@show rows; @show cols

eval_jac_g_ad(x_0, rows, cols, nothing)
eval_jac_g_ad(x_0, rows, cols, values)

##------------------------ Hessian of Langrangian  ------------------------##
# L = f(x) + g(x)^T * λ
# ∇L = σ_f  ∇^2 f(x_k) + Σ^m_i λ_i * ∇^2 g_i(x_k)
println("Hessian of Langrangian")
values = ones(10)

function eval_h(
    x::Vector{Float64},
    rows::Vector{Int32},
    cols::Vector{Int32},
    obj_factor::Float64,
    lambda::Vector{Float64},
    values::Union{Nothing,Vector{Float64}},
)
    if values === nothing
        # Symmetric matrix, fill the lower left triangle only
        idx = 1
        for row in 1:4
            for col in 1:row
                rows[idx] = row
                cols[idx] = col
                idx += 1
            end
        end
    else
        # Again, only lower left triangle
        # Objective
        values[1] = obj_factor * (2 * x[4])  # 1,1
        values[2] = obj_factor * (x[4])  # 2,1
        values[3] = 0                      # 2,2
        values[4] = obj_factor * (x[4])  # 3,1
        values[5] = 0                      # 3,2
        values[6] = 0                      # 3,3
        values[7] = obj_factor * (2 * x[1] + x[2] + x[3])  # 4,1
        values[8] = obj_factor * (x[1])  # 4,2
        values[9] = obj_factor * (x[1])  # 4,3
        values[10] = 0                     # 4,4

        # First constraint
        values[2] += lambda[1] * (x[3] * x[4])  # 2,1
        values[4] += lambda[1] * (x[2] * x[4])  # 3,1
        values[5] += lambda[1] * (x[1] * x[4])  # 3,2
        values[7] += lambda[1] * (x[2] * x[3])  # 4,1
        values[8] += lambda[1] * (x[1] * x[3])  # 4,2
        values[9] += lambda[1] * (x[1] * x[2])  # 4,3

        # Second constraint
        values[1] += lambda[2] * 2  # 1,1
        values[3] += lambda[2] * 2  # 2,2
        values[6] += lambda[2] * 2  # 3,3
        values[10] += lambda[2] * 2  # 4,4
    end
    return
end

println("analytical")
eval_h(x_0, Int32.(ones(10)), Int32.(ones(10)), 1., ones(2), nothing)
# Hessian of constraints
eval_h(x_0, Int32.(ones(10)), Int32.(ones(10)), 0., ones(2), values)
@show values

function L(x, λ, σ)
    λ = typeof(x)(λ)
    return σ * eval_f(x) + eval_g_ad(x,g)'*λ
end

λ = ones(2)
σ = 0.0
L(x_0, λ, σ)
∇H = x -> ForwardDiff.hessian(x->L(x, λ, σ), x)
sparsity = LowerTriangular(∇H(x_0)) .!= 0.0
I = ((i, j) for i = 1:size(sparsity,1), j = 1:size(sparsity,2) if sparsity[i,j] != 0.0)
rows_hess = Int32.(getindex.(I, 1)[:])
cols_hess = Int32.(getindex.(I, 2)[:])
values_hess = zeros(length(rows_hess))

function hess_structure_sparse!(x, rows, cols)
    ∇H = x -> ForwardDiff.hessian(x->L(x, λ, 1), x)
    sparsity = LowerTriangular(∇H(x_0)) .!= 0.0
    I = ((i, j) for i = 1:size(sparsity,1), j = 1:size(sparsity,2) if sparsity[i,j] != 0.0)
    rows .= Int32.(getindex.(I, 1)[:])
    cols .= Int32.(getindex.(I, 2)[:])
    return rows, cols
end

function eval_h_ad(
    x::Vector{Float64},
    rows,
    cols,
    σ::Float64,
    λ::Vector{Float64},
    values::Union{Nothing,Vector{Float64}},
)
    if values === nothing
        hess_structure_sparse!(x, rows, cols)
    else
        hess = ForwardDiff.hessian(x->L(x, λ, σ), x)
        values[:] = hess[tril(trues(size(hess)))]
    end
end
println("automatic differentiation")
eval_h_ad(x_0, rows_hess, cols_hess, σ, λ, values_hess)
@show rows_hess; @show cols_hess
eval_h_ad(x_0, rows_hess, cols_hess, σ, λ, values_hess)
@show values_hess

##------------------------ Solve problem analytical ------------------------##

n = 4
x_L = [1.0, 1.0, 1.0, 1.0]
x_U = [5.0, 5.0, 5.0, 5.0]

m = 2
g_L = [25.0, 40.0]
g_U = [2.0e19, 40.0]

prob = Ipopt.CreateIpoptProblem(
    n,
    x_L,
    x_U,
    m,
    g_L,
    g_U,
    8,
    10,
    eval_f,
    eval_g,
    eval_grad_f,
    eval_jac_g,
    eval_h,
)

prob.x = [1.0, 5.0, 5.0, 1.0]
solvestat = Ipopt.IpoptSolve(prob)
@show prob.x
##------------------------ solve with ad ------------------------##

n = 4
x_L = [1.0, 1.0, 1.0, 1.0]
x_U = [5.0, 5.0, 5.0, 5.0]

m = 2
g_L = [25.0, 40.0]
g_U = [2.0e19, 40.0]

g = zeros(m)
∇g = x -> ForwardDiff.jacobian(x->eval_g_ad(x, g), x)
jac_temp = ∇g(x_0)
sparsity = jac_temp .!= 0.0
I = ((i, j) for i = 1:m, j = 1:n if sparsity[i,j] != 0.0)
rows = Int32.(getindex.(I, 1)[:])
cols = Int32.(getindex.(I, 2)[:])
values = zeros(length(rows))
rows = nothing
cols = nothing
values = nothing

prob = Ipopt.CreateIpoptProblem(
    n,
    x_L,
    x_U,
    m,
    g_L,
    g_U,
    8,
    10,
    eval_f,
    eval_g,
    eval_grad_f_ad,
    eval_jac_g_ad,
    eval_h_ad,
    # function h() return false end
)
# Ipopt.AddIpoptStrOption(prob, "hessian_approximation", "limited-memory")
# Ipopt.AddIpoptStrOption(prob, "check_derivatives_for_naninf", "yes")
# Ipopt.AddIpoptIntOption(prob, "file_print_level", 10)
# Ipopt.AddIpoptIntOption(prob, "output_file", "ipopt.out")
# Ipopt.OpenIpoptOutputFile(prob, "blah.txt", 5)

prob.x = [1.0, 5.0, 5.0, 1.0]
solvestat = Ipopt.IpoptSolve(prob)
@show prob.x


text:
# ```@setup whitenoise
# # all examples are run on an example global.toml and downsized maps.
# ARGS = ["example.toml"]
# ```
#
# # [White Noise Levels (whitenoise.jl)](@id whitenoise)
# This file estimates the amplitude of the noise power spectrum from the
# pixel variance map, under the assumption that the noise is white.
# We want to compute
# ```math
# N_{\ell}^{\mathrm{white}} = \frac{1}{\sum_p m_p^2} \sum_p \sigma_{II}^2 \Omega_p m_p^2
# ```
# where ``m_p`` is mask at pixel ``p``, ``\Omega_p`` is the pixel area, and ``\sigma_{II}^2``
# is the pixel variance.
#
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia whitenoise.jl example.toml</code></pre>
# ```
#
# We first need to do the usual setup steps. We read the command-line arguments and load
# the packages we need.

configfile = first(ARGS)

## setup data
using Plots
using PowerSpectra
using Healpix
using TOML
using CSV
using DataFrames

config = TOML.parsefile(configfile)
nside = config["general"]["nside"]
run_name = config["general"]["name"]
spectrapath = joinpath(config["scratch"], "rawspectra")
lmax = min(2508,nside2lmax(nside))
npix = nside2npix(nside)
Ωp = 4π / npix

# Next, we loop over each half-mission frequency map, and estimate the white noise power
# spectrum.
df = DataFrame(freq = String[], split = String[], noiseT = Float64[], noiseP = Float64[])

for freq ∈ ("100", "143", "217"),  split ∈ ("1", "2")
    mapid = "P$(freq)hm$(split)"
    maskfileT = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid)_maskT.fits")
    maskfileP = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid)_maskP.fits")
    mapfile = joinpath(config["scratch"], "maps", config["map"][mapid])
    maskT = readMapFromFITS(maskfileT, 1, Float64)
    maskP = readMapFromFITS(maskfileP, 1, Float64)
    covII = nest2ring(readMapFromFITS(mapfile, 5, Float64)) * 1e12
    covQQ = nest2ring(readMapFromFITS(mapfile, 8, Float64)) * 1e12
    covUU = nest2ring(readMapFromFITS(mapfile, 10, Float64)) * 1e12
    N_white_T = sum(maskT.pixels.^2 .* (covII.pixels) .* Ωp) ./ (sum(maskT.pixels.^2))
    N_white_P = sum(maskP.pixels.^2 .* (covUU.pixels + covQQ.pixels) .* Ωp) ./ (2 * sum(maskP.pixels.^2))
    push!(df, (freq, split, N_white_T, N_white_P))
end

# Finally, we save this all to a CSV.
csvfile = joinpath(config["scratch"], "whitenoise.dat")
CSV.write(csvfile, df)

# Let's inspect the values that we saved.
df

max_forks_repo_name: sgiardie/PSpipe


text:
export simulate_unit!

# ----------------------
# Holds simulations
# ----------------------

type_to_sym(x) = x |> typeof |> Symbol |> string

const SimulationList = Tuple{DrillReward,AbstractUnitProblem,Vector}

"Create `DataFrame` of info given vector of `simulationPrimitives`"
function simulationPrimitives_information(simulations::Vector{<:SimulationList})
    df = DataFrame()
    df[!, :i] = 1:length(simulations)

    df[!, :tech]        = CategoricalArray( map(s -> s |> x -> getindex(x,1) |> revenue |> tech    |> type_to_sym, simulations) )
    df[!, :tax ]        = CategoricalArray( map(s -> s |> x -> getindex(x,1) |> revenue |> tax     |> type_to_sym, simulations) )
    df[!, :learn]       = CategoricalArray( map(s -> s |> x -> getindex(x,1) |> revenue |> learn   |> type_to_sym, simulations) )
    df[!, :royalty]     = CategoricalArray( map(s -> s |> x -> getindex(x,1) |> revenue |> royalty |> type_to_sym, simulations) )
    df[!, :state_space] = CategoricalArray( map(s -> s |> x -> getindex(x,2)                       |> type_to_sym, simulations) )
    df[!, :theta]       = CategoricalArray( map(s -> s |> x -> getindex(x,3) |> x -> round.(x; digits=3) |> string     , simulations) )
    return df
end

"construct `DataFrame`s to hold simulations"
function dataFramesForSimulations(simulations::Vector{<:SimulationList}, data::DataDrill, Tstop)

    nSim = length(simulations)
    nT = length(zchars(data))
    timestamps = _timestamp(zchars(data))
    nI = num_i(data)

    Dmx = _Dmax(getindex(first(simulations), 2))

    df_d_iter = product(1:nT,        1:nI, 1:nSim)
    df_D_iter = product(Tstop:Tstop, 1:nI, 1:nSim)

    ndf_d = length(df_d_iter)
    ndf_D = length(df_D_iter)

    df_d = DataFrame(
        sim     = vec(collect(Int32(ns)    for (t, i, ns,) in df_d_iter)),
        i       = vec(collect(Int32(i)     for (t, i, ns,) in df_d_iter)),
        mondate = vec(collect(timestamps[t] for (t, i, ns,) in df_d_iter)),
    )
    vars = [
        :d0, :d1, :d0psi, :d1psi, :d0eur, :d1eur, :d0eursq, :d1eursq, :d0eurcub, :d1eurcub,
        :epsdeq1, :epsdgt1, :Prdeq1, :Prdgt1,
        :Eeps, :profit, :surplus, :revenue, :drillcost, :extension,
    ]
    for v in vars
        df_d[!, v] = zeros(Float64, ndf_d)
    end

    df_D = DataFrame(
        sim     = vec(collect(Int32(ns)     for (t, i, ns,) in df_D_iter)),
        i       = vec(collect(Int32(i)      for (t, i, ns,) in df_D_iter)),
        mondate = vec(collect(timestamps[t]  for (t, i, ns,) in df_D_iter)),
    )
    for D in 0:Dmx
        df_D[!, Symbol("D$(D)")] = zeros(Float64, ndf_D)
    end

    return (df_d, df_D,)
end

@deprecate dataFrames_from_simulationPrimitives(args...) dataFramesForSimulations(args...)

# ----------------------
# Holds simulations
# ----------------------

"update `df_d` and `df_D` with simulation number `sim`, held in `drillsim`"
function update_sim_dataframes_from_simdata!(df_d::DataFrame, df_D::DataFrame, drillsim::SharedSimulations, sim::Integer)

    is_this_simulation_d = df_d[!,:sim] .== sim
    is_this_simulation_D = df_D[!,:sim] .== sim

    sum(is_this_simulation_d) == 0 && throw(error("no observations selected to update"))

    # update period-t simulations
    df_d[is_this_simulation_d, :d0       ] = vec(sdata(drillsim.d0))
    df_d[is_this_simulation_d, :d1       ] = vec(sdata(drillsim.d1))
    df_d[is_this_simulation_d, :d0psi    ] = vec(sdata(drillsim.d0psi))
    df_d[is_this_simulation_d, :d1psi    ] = vec(sdata(drillsim.d1psi))
    df_d[is_this_simulation_d, :d0eur    ] = vec(sdata(drillsim.d0eur))
    df_d[is_this_simulation_d, :d1eur    ] = vec(sdata(drillsim.d1eur))
    df_d[is_this_simulation_d, :d0eursq  ] = vec(sdata(drillsim.d0eursq))
    df_d[is_this_simulation_d, :d1eursq  ] = vec(sdata(drillsim.d1eursq))
    df_d[is_this_simulation_d, :d0eurcub ] = vec(sdata(drillsim.d0eurcub))
    df_d[is_this_simulation_d, :d1eurcub ] = vec(sdata(drillsim.d1eurcub))
    df_d[is_this_simulation_d, :epsdeq1  ] = vec(sdata(drillsim.epsdeq1))
    df_d[is_this_simulation_d, :epsdgt1  ] = vec(sdata(drillsim.epsdgt1))
    df_d[is_this_simulation_d, :Prdeq1   ] = vec(sdata(drillsim.Prdeq1))
    df_d[is_this_simulation_d, :Prdgt1   ] = vec(sdata(drillsim.Prdgt1))
    df_d[is_this_simulation_d, :Eeps     ] = vec(sdata(drillsim.Eeps))
    df_d[is_this_simulation_d, :profit   ] = vec(sdata(drillsim.profit))
    df_d[is_this_simulation_d, :surplus  ] = vec(sdata(drillsim.surplus))
    df_d[is_this_simulation_d, :revenue  ] = vec(sdata(drillsim.revenue))
    df_d[is_this_simulation_d, :drillcost] = vec(sdata(drillsim.drillcost))
    df_d[is_this_simulation_d, :extension] = vec(sdata(drillsim.extension))

    # update final period simulations
    df_D[is_this_simulation_D, :D0] = 1 .- vec(sum(sdata(drillsim.D_at_T);dims=1))
    for D in 1:size(drillsim.D_at_T, 1)
        @views df_D[is_this_simulation_D, Symbol("D$(D)")] = sdata(drillsim.D_at_T)[D,:]
    end
end

function doSimulations(datafull::DataSetofSets, simlist::Vector{<:SimulationList}, Tstop, M; do_royalty=true)

    N = num_i(datafull)
    wrkrs = CachingPool(workers())

    # generate smaller DataDrill to transfer to workers
    datadrill = drill(datafull)
    ddm_novf = DDM_NoVF(_model(datadrill))
    datadrill_dso = DataDrillStartOnly(datadrill)
    datadrill_bare= DataDrill(datadrill_dso, ddm_novf)

    # DataFull set for transfer
    data_bare = DataSetofSets(datafull, datadrill_bare)

    # generate common variables
    sharesim = SharedSimulations(datadrill_bare)
    simMat = SimulationDraws(data_bare, M)

    # DataFrames for results
    df_t, df_Tstop = dataFramesForSimulations(simlist, datadrill_bare, Tstop)

    # send to workers
    @eval @everywhere set_g_BaseDataSetofSets($data_bare)
    @eval @everywhere set_g_SimulationDrawsMatrix($simMat)
    @eval @everywhere set_g_SharedSimulations($sharesim)

    for (k, (rwrd, wp, theta)) in enumerate(simlist)

        check_theta(datafull, theta)
        fill!(sharesim, 0)
        rev = revenue(rwrd)
        thet_d = theta_drill(datafull, theta)
        thet_r = vw_revenue(rwrd, thet_d)
        println_time_flush("Simulation $k of $(length(simlist))")
        println("\tTech = $(tech(rev))")
        println("\tLearn = $(learn(rev))")
        println("\tRoyalty = $(royalty(rev))")
        println("\tStatespace = $(wp)")
        println("\tTheta_rev = $(thet_r)")

        @eval @everywhere set_g_SimulationPrimitives($rwrd, $wp, $Tstop, $theta)
        simprim = get_g_SimulationPrimitives()

        # map(i -> simulate_unit!(simprim, i, do_royalty), OneTo(N))
        pmap(i -> simulate_unit!(i, do_royalty), wrkrs, OneTo(N))
        update_sim_dataframes_from_simdata!(df_t, df_Tstop, sharesim, k)
    end

    return df_t, df_Tstop
end

function theta_cost(d::DataSetofSets, theta)
    length(theta) == _nparm(d) || throw(DimensionMismatch())
    rwrd = reward(_model(drill(d)))
    theta_d = theta_drill(d, theta)
    thet_c = vw_cost(rwrd, theta_d)
    return thet_c
end

function theta_revenue(d::DataSetofSets, theta)
    length(theta) == _nparm(d) || throw(DimensionMismatch())
    rwrd = reward(_model(drill(d)))
    theta_d = theta_drill(d, theta)
    thet_r = vw_revenue(rwrd, theta_d)
    return thet_r
end

function average_cost(dfull::DataSetofSets, theta, d::Int, x=length(statespace(_model(drill(dfull)))))
    theta_c = theta_cost(dfull, theta)
    ddrill = drill(dfull)
    zs = zchars(ddrill)
    m = _model(ddrill)
    C = cost(reward(m))

    sim = SimulationDraw(0,0,0,zeros(Int,0))
    grad = similar(theta_c)

    obs(z) = ObservationDrill(m, (Inf,Inf), z, 0, x)
    drillcost(z) = flow!(grad, C, d, obs(z), theta_c, sim, false)

    avgcost = collect(drillcost(z) / d for z in _timevars(zs))

    return avgcost
end

function average_cost_df(d::DataSetofSets, theta)

    zs = zchars(drill(d))
    cost1 = average_cost(d, theta, 1, 1)
    cost2 = average_cost(d, theta, 2)

    df = DataFrame(
        date = _timestamp(zs),
        cost1 = cost1,
        cost2 = cost2
    )
    return df
end

function Theta_NoTech(d::DataSetofSets, theta::Vector, TECH_YEAR_ZERO::Integer)
    length(theta) == _nparm(d) || throw(DimensionMismatch())

    theta_notech = copy(theta)

    ddrill = drill(d)
    m = _model(ddrill)
    R = revenue(reward(m))

    thet_r = theta_revenue(d, theta_notech)
    basey = baseyear(tech(R))
    alphat = thet_r[idx_t(R)]

    if tech(R) isa TimeTrend
        thet_r[idx_0(R)] -= alphat*(basey - TECH_YEAR_ZERO)
    end
    for i in idx_t(R)
        thet_r[i] = 0
    end

    return theta_notech
end


text:
using LispSyntax
using Test

#-------------------------------------------------------------------------------
# Setup
macro incr(x)
  esc(quote
    $x = $x + 1
    $x
  end)
end

#-------------------------------------------------------------------------------
@testset "Reader" begin
    @test LispSyntax.read("1.1f") == 1.1f0
    @test LispSyntax.read("1.2f") == 1.2f0
    @test LispSyntax.read("2f")   == 2f0

    @test LispSyntax.read("3.0d") == 3.0

    @test LispSyntax.read("4")    == 4

    @test LispSyntax.read("\\u2312") == '\u2312'

    @test LispSyntax.read("\\040") == ' '

    @test LispSyntax.read("\\c") == 'c'

    @test LispSyntax.read("\"test\"") == "test"

    @test LispSyntax.read("true") == true
    @test LispSyntax.read("false") == false

    @test LispSyntax.read("test") == :test

    @test LispSyntax.read("()") == sx()
    @test LispSyntax.read("(1.1f)") == sx(1.1f0)
    @test LispSyntax.read("(1.1f 2.2f)") == sx(1.1f0, 2.2f0)
    @test LispSyntax.read("(+ 1.1f 2)") == sx(:+, 1.1f0, 2)
    @test LispSyntax.read("(this (+ 1.1f 2))") == sx(:this, sx(:+, 1.1f0, 2))
    @test LispSyntax.read("(this (+ 1.1f 2) )") == sx(:this, sx(:+, 1.1f0, 2))

    @test LispSyntax.read("#{1 2 3 4}") == Set([1, 2, 3, 4])

    @test LispSyntax.read("""#{
                            1 2
                            3 4
                            }""") == Set([1, 2, 3, 4])

    @test LispSyntax.read("{a 2 b 3}") == Dict(:a => 2, :b => 3)

    @test LispSyntax.read("""{
                            a 2
                            b 3
                            }""") == Dict(:a => 2, :b => 3)

    @test LispSyntax.read("[1 2 3 4]")  == sx(1, 2, 3, 4)

    @test LispSyntax.read("""[
                            1 2
                            3 4
                            ]""")  == sx(1, 2, 3, 4)

    @test LispSyntax.read("[]")         == sx()
    @test LispSyntax.read("[1]")        == sx(1)

    @test LispSyntax.read("'test")      == sx(:quote, :test)

    @test LispSyntax.read("`test")      == sx(:quasi, :test)

    @test LispSyntax.read("~test")      == sx(:splice, :test)
    @test LispSyntax.read("~@(1 2 3)")  == sx(:splice_seq, sx(1, 2, 3))

    @test LispSyntax.read("`~test")     == sx(:quasi, sx(:splice, :test))

    @test desx(sx(:splice_seq, sx(1, 2, 3))) == Any[ :splice_seq, [1, 2, 3] ]
    @test desx(sx(:splice_seq, sx(1, 2, sx(3)))) == Any[ :splice_seq, Any[ 1, 2, [3] ] ]

    @test LispSyntax.read("""(defn multiline
                               [x]
                               (+ x 1))""") == sx(:defn, :multiline, sx(:x), sx(:+, :x, 1))

    @test LispSyntax.read("""
    (defn f1 [n]
       (if (< n 2)
           1
           (+ (f1 (- n 1))
              (f1 (- n 2)))))
    """) == sx(:defn, :f1, sx(:n),
               sx(:if, sx(:<, :n, 2),
                  1,
                  sx(:+, sx(:f1, sx(:-, :n, 1)), sx(:f1, sx(:-, :n, 2)))))

    assign_reader_dispatch(:sx, x -> sx(x.vector...))
    assign_reader_dispatch(:hash, x -> Dict(x.vector[i] => x.vector[i+1] for i = 1:2:length(x.vector)))
    @test LispSyntax.read("#sx[a b c]") == sx(:a, :b, :c)
    @test LispSyntax.read("#sx [ 1 2 3 ]") == sx(1, 2, 3)
end

#-------------------------------------------------------------------------------
@testset "Code generation" begin
    @test codegen(desx(LispSyntax.read("(if true a)"))) == :(true && a)
    @test codegen(desx(LispSyntax.read("(if true a b)"))) == :(true ? a : b)

    @test codegen(desx(LispSyntax.read("(call)"))) == :(call())
    @test codegen(desx(LispSyntax.read("(call a)"))) == :(call(a))
    @test codegen(desx(LispSyntax.read("(call a b)"))) == :(call(a,b))
    @test codegen(desx(LispSyntax.read("(call a b c)"))) == :(call(a,b,c))

    @test codegen(desx(LispSyntax.read("(lambda (x) (call x))"))) == Base.remove_linenums!(:(function (x) call(x) end))
    @test codegen(desx(LispSyntax.read("(def x 3)"))) == :(global x = 3)
    @test codegen(desx(LispSyntax.read("(def x (+ 3 1))"))) == :(global x = 3 + 1)

    construct_sexpr = LispSyntax.construct_sexpr
    @test codegen(desx(LispSyntax.read("test"))) == :test
    @test codegen(desx(LispSyntax.read("'test"))) == QuoteNode(:test)
    @test codegen(desx(LispSyntax.read("'(1 2)"))) == :($construct_sexpr(1, 2))
    @test codegen(desx(LispSyntax.read("'(1 x)"))) == :($construct_sexpr(1, :x))
    @test codegen(desx(LispSyntax.read("'(1 (1 2))"))) == :($construct_sexpr(1, $construct_sexpr(1, 2)))
    @test codegen(desx(LispSyntax.read("'(1 (test x))"))) == :($construct_sexpr(1, $construct_sexpr(:test, :x)))
    @test codegen(desx(LispSyntax.read("(call 1 '2)"))) == :(call(1,2))
end

#-------------------------------------------------------------------------------
@testset "Scope and variables" begin
    x = 10
    @test lisp"x" == 10

    let
        # In clojure, def affects global bindings
        lisp"(def w (+ 3 1))"
    end
    let
        @test w == 4
    end
end

#-------------------------------------------------------------------------------
@testset "Quoting and splicing" begin
    x = 10
    @test lisp"`~x" == 10
    @test lisp"'test" == :test
    @test lisp"'(1 2)" == Any[1, 2]
    @test lisp"'(1 x)" == Any[1, :x]
    @test lisp"'(1 (1 2))" == Any[1, Any[1, 2]]
    @test lisp"'(1 (test x))" == Any[1, Any[:test, :x]]
    @test lisp"`(test ~x)" == Any[ :test, 10 ]
    @test lisp"`(~x ~x)" == Any[ 10, 10 ]
    global y = Any[ 1, 2 ]
    @test lisp"`(~x ~@y)" == Any[ 10, 1, 2 ]
    @test lisp"`(~x ~y)" == Any[ 10, Any[1, 2] ]

    @test lisp"`(10 ~(+ 10 x))" == Any[10, 20]

    @test lisp"(quote (+ 1 2))" == Any[:+, 1, 2]
end

#-------------------------------------------------------------------------------
@testset "Functions" begin
    lisp"(defn xxx [a b] (+ a b))"
    @test lisp"(xxx 1 2)" == 3

    global z = 10
    lisp"(defn yyy [a] (+ a z))"
    @test lisp"(yyy 1)" == 11
    @test lisp"(yyy z)" == 20

    # recursion
    lisp"(defn fib [a] (if (< a 2) a (+ (fib (- a 1)) (fib (- a 2)))))"
    @test lisp"(fib 2)" == 1
    @test lisp"(fib 4)" == 3
    @test lisp"(fib 30)" == 832040
    @test lisp"(fib 40)" == 102334155

    # Note this version is slow due to the non-const global binding fib2
    lisp"(def fib2 (lambda [a] (if (< a 2) a (+ (fib2 (- a 1)) (fib2 (- a 2))))))"
    @test lisp"(fib2 2)" == 1
    @test lisp"(fib2 4)" == 3
    @test lisp"(fib2 30)" == 832040

    lisp"(defn dostuff [a] (@incr a) (@incr a) (@incr a))"
    @test lisp"(dostuff 3)" == 6
    @test lisp"(dostuff 6)" == 9

    lisp"(def dostuff2 (lambda [a] (@incr a) (@incr a) (@incr a)))"
    @test lisp"(dostuff2 3)" == 6
    @test lisp"(dostuff2 6)" == 9

    lisp"(def dostuff3 (fn [a] (@incr a) (@incr a) (@incr a)))"
    @test lisp"(dostuff3 3)" == 6
    @test lisp"(dostuff3 6)" == 9
    @test lisp"((lambda [x] (+ x 1)) 5)" == 6
    @test lisp"#{1 2 z}" == Set([1, 2, 10])
    @test lisp"{1 2 2 z}" == Dict(1 => 2, 2 => 10)
    @test lisp"#sx[+ 1 2]" == 3
    @test lisp"#hash['+ 1 '- z]" == Dict(:+ => 1, :- => 10)
end

#-------------------------------------------------------------------------------
# Macros. Note that LispSyntax currently uses the Julia macro system, which
# includes the hygenic-by-default hygiene rules. This is very different from
# Clojure style explicit hygiene.
lisp"(defmacro fapply [f a] (esc `(~f ~a)))"
fcount = 0
lisp"(defmacro fapply_trace [f a] (esc `(do (global fcount) (@incr fcount) (~f ~a))))"
@testset "Macros" begin
    lisp"(defn fact [a] (if (< a 1) 1 (* a (fact (- a 1)))))"
    @test @fapply(fib2, 2) == 1
    @test @fapply(fact, 3 + 1) == 24
    @test lisp"(@fapply fib2 2)" == 1
    @test lisp"(@fapply fact (+ 3 1))" == 24

    @test @fapply_trace(fib2, 2) == 1
    @test fcount == 1
    @test @fapply_trace(fact, 3 + 1) == 24
    @test fcount == 2
end

#-------------------------------------------------------------------------------
@testset "Loops" begin
    number = 0
    output = 0

    lisp"(while (< number 2) (@incr number) (@incr output))"
    @test number == 2
    @test output == 2
    r = output
    lisp"(for [i (: 1 10)] (@incr r))"
    @test r == 12

    r = 0
    lisp"(for [i (: 1 10) j (: 1 10)] (@incr r))"
    @test r == 100
end

#-------------------------------------------------------------------------------
@testset "Let and do" begin
    number = 2
    r = 100
    output = 2
    @test lisp"(let [x 10] x)" == 10
    @test lisp"(let [x 10 y 20] (+ x y))" == 30
    @test lisp"(let [x 10 y 20 z 20] (+ x y z))" == 50
    @test lisp"(let [x 10 y 20 z 20] (+ x y z number))" == 52
    @test lisp"(let [x 10 y 20 z 20 number 10] (+ x y z number))" == 60
    @test lisp"(let [x 10 y 20 z 20] (- (+ x y z number) output))" == 50

    lisp"(do (@incr r) (@incr number))"
    @test number == 3
    @test r == 101
end

#-------------------------------------------------------------------------------
lisp"(import ParserCombinator)"
@testset "Module import" begin
    @test lisp"(@E_str \"S\")" == E"S"
end

#-------------------------------------------------------------------------------
@testset "Include from file" begin
    # Return value is value of last expression
    @test include_lisp(@__MODULE__, "lisp.clj") == 100
    # Test objects defined in lisp.clj
    @test func_in_clj_file(1, 2) == "x = 1; y = 2"
    @test func_in_clj_file(10, 20) == "x = 10; y = 20"
    @test some_global === 1.23f0
    @test !isdefined(@__MODULE__, :not_a_global)
end

#-------------------------------------------------------------------------------
# Bug reports
@testset "Bug reports" begin
    @test lisp"""(def game_map (Dict
              (=> 'living_room
                  '((you are in the living room
                     of a wizards house - there is a wizard
                     snoring loudly on the couch -)
                    (west door garden)
                    (upstairs stairway attic)))))""" ==
                    Dict(:living_room =>
                         Any[ Any[ :you, :are, :in, :the, :living, :room, :of, :a, :wizards, :house, :-,
                                  :there, :is, :a, :wizard, :snoring, :loudly, :on, :the, :couch, :- ],
                             Any[ :west, :door, :garden ],
                             Any[ :upstairs, :stairway, :attic ] ])
end


text:
function anderson_mix(vin, vout, beta, ymat, smat, iter::Int64, mixdim)
    # function that computes the new potential usign the anderson mix
    # TODO: I think that the mixing can be efficiently be explained by a small
    # jupyter notebook

    n = length(vin);
    #
    # function evaluation overwrites vout
    #
    # residue
    res = vin - vout;
    #
    iterused = min(iter-1,mixdim);
    ipos = iter - 1 - round(Integer, floor((iter-2)/mixdim))*mixdim;
    #
    if (iter > 1)
        # compute the changes in function evaluations and the step (changes
        # in potentials)
        ymat[:,ipos] = res - ymat[:,ipos];
        smat[:,ipos] = vin - smat[:,ipos];
    end
    #
    vopt  = vin;
    ropt  = res;
    #
    if (iter > 1)
        # Solve the least square problem
        gammas = ymat[:,1:iterused] \ res;
        vopt   = vin - smat[:,1:iterused] * gammas;
        ropt   = res - ymat[:,1:iterused] * gammas;
    end

    inext = iter - round(Integer, floor((iter - 1) / mixdim)) * mixdim;
    ymat[:,inext] = res;
    smat[:,inext] = vin;

    vnew = vopt - beta*ropt;

    return (vnew, ymat, smat)
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
"""
    mutable struct SalientPoleQuadratic <: Machine
        base_machine::SalientPoleMachine
        saturation_coeffs::Tuple{Float64, Float64}

3-states salient-pole synchronous machine with exponential saturation:
IEEE Std 1110 §5.3.2 (Model 2.1). GENSAL in PSSE and PSLF.

# Arguments:
- `base_machine::SalientPoleMachine`: Salient Pole Machine model.
- `saturation_coeffs::Tuple{Float64, Float64}``: Saturation coefficients for quadratic model.
"""
mutable struct SalientPoleQuadratic <: Machine
    base_machine::SalientPoleMachine
    saturation_coeffs::Tuple{Float64, Float64}
end
IS.@forward((SalientPoleQuadratic, :base_machine), SalientPoleMachine)

function SalientPoleQuadratic(
    R::Float64,
    Td0_p::Float64,
    Td0_pp::Float64,
    Tq0_pp::Float64,
    Xd::Float64,
    Xq::Float64,
    Xd_p::Float64,
    Xd_pp::Float64,
    Xl::Float64,
    Se::Tuple{Float64, Float64},
)
    saturation_coeffs = get_quadratic_saturation(Se)
    return SalientPoleQuadratic(
        SalientPoleMachine(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se),
        saturation_coeffs,
    )
end

function SalientPoleQuadratic(; R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se)
    return SalientPoleQuadratic(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se)
end

function SalientPoleQuadratic(::Nothing)
    return SalientPoleQuadratic(;
        R = 0.0,
        Td0_p = 0.0,
        Td0_pp = 0.0,
        Tq0_pp = 0.0,
        Xd = 0.0,
        Xq = 0.0,
        Xd_p = 0.0,
        Xd_pp = 0.0,
        Xl = 0.0,
        Se = (0.0, 0.0),
    )
end

get_base_machine(value::SalientPoleQuadratic) = value.base_machine
get_saturation_coeffs(value::SalientPoleQuadratic) = value.saturation_coeffs
set_base_machine!(value::SalientPoleQuadratic, val::SalientPoleMachine) =
    value.base_machine = val
set_saturation_coeffs!(value::SalientPoleQuadratic, val::Tuple{Float64, Float64}) =
    value.saturation_coeffs = val

function IS.deserialize_struct(::Type{SalientPoleQuadratic}, data::Dict)
    vals = IS.deserialize_to_dict(SalientPoleQuadratic, data)
    return SalientPoleQuadratic(vals[:base_machine], vals[:saturation_coeffs])
end


text:
struct COSMOSAC02Param <: EoSParam
    Pi::SingleParam{Vector{Float64}}
    V::SingleParam{Float64}
    A::SingleParam{Float64}
end

abstract type COSMOSAC02Model <: ActivityModel end

struct COSMOSAC02{c<:EoSModel} <: COSMOSAC02Model
    components::Array{String,1}
    icomponents::UnitRange{Int}
    params::COSMOSAC02Param
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel COSMOSAC02
export COSMOSAC02

function COSMOSAC02(components; puremodel=PR,
    userlocations=String[],
     verbose=false)
    params = getparams(components, ["Activity/COSMOSAC/COSMOSAC02_like.csv"]; userlocations=userlocations, verbose=verbose)
    Pi  = COSMO_parse_Pi(params["Pi"])
    A  = params["A"]
    V  = params["V"]
    icomponents = 1:length(components)

    init_puremodel = [puremodel([components[i]]) for i in icomponents]
    packagedparams = COSMOSAC02Param(Pi,V,A)
    references = String[]
    model = COSMOSAC02(components,icomponents,packagedparams,init_puremodel,1e-12,references)
    return model
end

function activity_coefficient(model::COSMOSAC02Model,V,T,z)
    return exp.(@f(lnγ_comb) .+ @f(lnγ_res))
end

function lnγ_comb(model::COSMOSAC02Model,V,T,z)
    r0 = 66.69
    q0 = 79.53

    V = model.params.V.values
    A = model.params.A.values

    x = z ./ sum(z)

    r = V ./ r0
    q = A ./ q0

    Φ = r/sum(x[i]*r[i] for i ∈ @comps)
    θ = q/sum(x[i]*q[i] for i ∈ @comps)
    lnγ_comb = @. log(Φ)+(1-Φ)-5*q*(log(Φ/θ)+(1-Φ/θ))
    return lnγ_comb
end

function lnγ_res(model::COSMOSAC02Model,V,T,z)
    x = z ./ sum(z)

    aeff = 7.5
    A = model.params.A.values
    n = A ./ aeff
    Pi = model.params.Pi.values
    PS = sum(x[i]*Pi[i][:] for i ∈ @comps) ./ sum(x[i]*A[i] for i ∈ @comps)

    lnΓS = @f(lnΓ,PS)
    lnΓi = [@f(lnΓ,Pi[i]./A[i]) for i ∈ @comps]
    lnγ_res_ =  [n[i]*sum(Pi[i][v]/A[i]*(lnΓS[v]-lnΓi[i][v]) for v ∈ 1:51) for i ∈ @comps]

    return lnγ_res_
end

function lnΓ(model::COSMOSAC02Model,V,T,z,P)
    Γ0 = ones(length(P))
    σ  = -0.025:0.001:0.025
    Γold = exp.(-log.(sum(P[i]*Γ0[i]*exp.(-ΔW.(σ,σ[i])./T) for i ∈ 1:51)))
    Γnew = Γold
    tol = 1
    i = 1
    damp_factor = 0.5
    while tol>sqrt(model.absolutetolerance)
        Γnew = exp.(-log.(sum(P[i]*Γold[i]*exp.(-ΔW.(σ,σ[i])./T) for i ∈ 1:51)))
        Γnew .*= (1-damp_factor)
        Γnew .+= damp_factor .* Γold
        tol = cosmo_tol(Γnew,Γold)
        Γold .= Γnew
        i+=1
    end
    lnΓ = log.(Γold)
    return lnΓ
end

function ΔW(σm,σn)
    σdon,σacc = minmax(σm,σn)
    σhb  = 0.0084
    chb  = 85580
    α    = 16466.72
    R    = 0.001987
    return (α/2*(σm+σn)^2+chb*max(0,σacc-σhb)*min(0,σdon+σhb))/R
end


text:
# This file was generated, do not modify it. # hide
y = smarket.Direction
X = select(smarket, Not(:Direction));


text:
function settype!(b::Binding, type::Binding)
    push!(type.refs, b)
    b.type = type
end

function settype!(b::Binding, type)
    b.type = type
end

function infer_type(binding::Binding, scope, state)
    if binding isa Binding
        binding.type !== nothing && return
        if binding.val isa EXPR && CSTParser.defines_module(binding.val)
            settype!(binding, CoreTypes.Module)
        elseif binding.val isa EXPR && CSTParser.defines_function(binding.val)
            settype!(binding, CoreTypes.Function)
        elseif binding.val isa EXPR && CSTParser.defines_datatype(binding.val)
            settype!(binding, CoreTypes.DataType)
        elseif binding.val isa EXPR
            if isassignment(binding.val)
                if CSTParser.is_func_call(binding.val.args[1])
                    settype!(binding, CoreTypes.Function)
                else
                    infer_type_assignment_rhs(binding, state, scope)
                end
            elseif binding.val.head isa EXPR && valof(binding.val.head) == "::"
                infer_type_decl(binding, state, scope)
            elseif iswhere(parentof(binding.val))
                settype!(binding, CoreTypes.DataType)
            end
        end
    end
end

function infer_type_assignment_rhs(binding, state, scope)
    rhs = binding.val.args[2]
    if is_loop_iter_assignment(binding.val)
        settype!(binding, infer_eltype(rhs))
    elseif headof(rhs) === :ref && length(rhs.args) > 1
        ref = refof_maybe_getfield(rhs.args[1])
        if ref isa Binding && ref.val isa EXPR
            settype!(binding, infer_eltype(ref.val))
        end
    else
        if CSTParser.is_func_call(rhs)
            callname = CSTParser.get_name(rhs)
            if isidentifier(callname)
                resolve_ref(callname, scope, state)
                if hasref(callname)
                    rb = get_root_method(refof(callname), state.server)
                    if (rb isa Binding && (CoreTypes.isdatatype(rb.type) || rb.val isa SymbolServer.DataTypeStore)) || rb isa SymbolServer.DataTypeStore
                        settype!(binding, rb)
                    end
                end
            end
        elseif headof(rhs) === :INTEGER
            settype!(binding, CoreTypes.Int)
        elseif headof(rhs) === :HEXINT
            if length(rhs.val) < 5
                settype!(binding, CoreTypes.UInt8)
            elseif length(rhs.val) < 7
                settype!(binding, CoreTypes.UInt16)
            elseif length(rhs.val) < 11
                settype!(binding, CoreTypes.UInt32)
            else
                settype!(binding, CoreTypes.UInt64)
            end
        elseif headof(rhs) === :FLOAT
            settype!(binding, CoreTypes.Float64)
        elseif CSTParser.isstringliteral(rhs)
            settype!(binding, CoreTypes.String)
        elseif headof(rhs) === :TRUE || headof(rhs) === :FALSE
            settype!(binding, CoreTypes.Bool)
        elseif isidentifier(rhs) || is_getfield_w_quotenode(rhs)
            refof_rhs = isidentifier(rhs) ? refof(rhs) : refof_maybe_getfield(rhs)
            if refof_rhs isa Binding
                if refof_rhs.val isa SymbolServer.GenericStore && refof_rhs.val.typ isa SymbolServer.FakeTypeName
                    settype!(binding, maybe_lookup(refof_rhs.val.typ.name, state))
                elseif refof_rhs.val isa SymbolServer.FunctionStore
                    settype!(binding, CoreTypes.Function)
                elseif refof_rhs.val isa SymbolServer.DataTypeStore
                    settype!(binding, CoreTypes.DataType)
                else
                    settype!(binding, refof_rhs.type)
                end
            elseif refof_rhs isa SymbolServer.GenericStore && refof_rhs.typ isa SymbolServer.FakeTypeName
                settype!(binding, maybe_lookup(refof_rhs.typ.name, state))
            elseif refof_rhs isa SymbolServer.FunctionStore
                settype!(binding, CoreTypes.Function)
            elseif refof_rhs isa SymbolServer.DataTypeStore
                settype!(binding, CoreTypes.DataType)
            end
        end
    end
end

function infer_type_decl(binding, state, scope)
    t = binding.val.args[2]
    if isidentifier(t)
        resolve_ref(t, scope, state)
    end
    if iscurly(t)
        t = t.args[1]
        resolve_ref(t, scope, state)
    end
    if CSTParser.is_getfield_w_quotenode(t)
        resolve_getfield(t, scope, state)
        t = t.args[2].args[1]
    end
    if refof(t) isa Binding
        rb = get_root_method(refof(t), state.server)
        if rb isa Binding && CoreTypes.isdatatype(rb.type)
            settype!(binding, rb)
        else
            settype!(binding, refof(t))
        end
    elseif refof(t) isa SymbolServer.DataTypeStore
        settype!(binding, refof(t))
    end
end

# Work out what type a bound variable has by functions that are called on it.
function infer_type_by_use(b::Binding, env::ExternalEnv)
    b.type !== nothing && return # b already has a type
    possibletypes = []
    visitedmethods = []
    ifbranch = nothing
    for ref in b.refs
        new_possibles = []
        ref isa EXPR || continue # skip non-EXPR (i.e. used for handling of globals)
        # Some simple handling for :if blocks
        if ifbranch === nothing
            ifbranch = find_if_parents(ref)
        else
            newbranch = find_if_parents(ref)
            if !in_same_if_branch(ifbranch, newbranch)
                return
            end
            ifbranch = newbranch
        end
        check_ref_against_calls(ref, visitedmethods, new_possibles, env)
        if !isempty(new_possibles)
            if isempty(possibletypes)
                possibletypes = new_possibles
            else
                possibletypes = intersect(possibletypes, new_possibles)
            end
            if isempty(possibletypes)
                return
            end
        end
    end
    # Only do something if we're left with a singleton set at the end.
    if length(possibletypes) == 1
        type = first(possibletypes)
        if type isa Binding
            settype!(b, type)
        elseif type isa SymbolServer.DataTypeStore
            settype!(b, type)
        elseif type isa SymbolServer.VarRef
            settype!(b, SymbolServer._lookup(type, getsymbols(env))) # could be nothing
        elseif type isa SymbolServer.FakeTypeName && isempty(type.parameters)
            settype!(b, SymbolServer._lookup(type.name, getsymbols(env))) # could be nothing
        end
    end
end

function check_ref_against_calls(x, visitedmethods, new_possibles, env::ExternalEnv)
    if is_arg_of_resolved_call(x) && !call_is_func_sig(x.parent)
        sig = parentof(x)
        # x is argument of function call (func) and we know what that function is
        if CSTParser.isidentifier(sig.args[1])
            func = refof(sig.args[1])
        else
            func = refof(sig.args[1].args[2].args[1])
        end
        argi = get_arg_position_in_call(sig, x) # what slot does ref sit in?
        tls = retrieve_toplevel_scope(x)
        if func isa Binding
            for method in func.refs
                method = get_method(method)
                method === nothing && continue
                if method isa EXPR
                    if defines_function(method)
                        get_arg_type_at_position(method, argi, new_possibles)
                    # elseif CSTParser.defines_struct(method)
                        # Can we ignore this? Default constructor gives us no type info?
                    end
                else # elseif what?
                    iterate_over_ss_methods(method, tls, env, m -> (get_arg_type_at_position(m, argi, new_possibles);false))
                end
            end
        else
            iterate_over_ss_methods(func, tls, env, m -> (get_arg_type_at_position(m, argi, new_possibles);false))
        end
    end
end

function call_is_func_sig(call::EXPR)
    # assume initially called on a :call
    if call.parent isa EXPR
        if call.parent.head === :function || CSTParser.is_eq(call.parent.head)
            true
        elseif isdeclaration(call.parent) || iswhere(call.parent)
            call_is_func_sig(call.parent)
        else
            false
        end
    else
        false
    end
end

function is_arg_of_resolved_call(x::EXPR)
    parentof(x) isa EXPR && headof(parentof(x)) === :call && # check we're in a call signature
    (caller = parentof(x).args[1]) !== x && # and that x is not the caller
    ((CSTParser.isidentifier(caller) && hasref(caller)) || (is_getfield(caller) && headof(caller.args[2]) === :quotenode && hasref(caller.args[2].args[1])))
end

function get_arg_position_in_call(sig::EXPR, arg)
    for i in 1:length(sig.args)
        sig.args[i] == arg && return i
    end
end

function get_arg_type_at_position(method, argi, types)
    if method isa EXPR
        sig = CSTParser.get_sig(method)
        if sig !== nothing &&
            sig.args !== nothing && argi <= length(sig.args) &&
            hasbinding(sig.args[argi]) &&
            (argb = bindingof(sig.args[argi]); argb isa Binding && argb.type !== nothing) &&
            !(argb.type in types)
            push!(types, argb.type)
            return
        end
    elseif method isa SymbolServer.DataTypeStore || method isa SymbolServer.FunctionStore
        for m in method.methods
            get_arg_type_at_position(m, argi, types)
        end
    end
    return
end

function get_arg_type_at_position(m::SymbolServer.MethodStore, argi, types)
    if length(m.sig) >= argi && m.sig[argi][2] != SymbolServer.VarRef(SymbolServer.VarRef(nothing, :Core), :Any) && !(m.sig[argi][2] in types)
        push!(types, m.sig[argi][2])
    end
end

# Assumes x.head.val == "="
is_loop_iter_assignment(x::EXPR) = x.parent isa EXPR && ((x.parent.head == :for || x.parent.head == :generator) || (x.parent.head == :block && x.parent.parent isa EXPR && (x.parent.parent.head == :for || x.parent.parent.head == :generator)))

function infer_eltype(x::EXPR)
    if isidentifier(x) && hasref(x) # assume is IDENT
        r = refof(x)
        if r isa Binding && r.val isa EXPR
            if isassignment(r.val) && r.val.args[2] != x
                return infer_eltype(r.val.args[2])
            end
        end
    elseif headof(x) === :ref && hasref(x.args[1])
        r = refof(x.args[1])
        if r isa SymbolServer.DataTypeStore ||
            r isa Binding && CoreTypes.isdatatype(r.type)
            r
        end
    elseif headof(x) === :STRING
        return CoreTypes.Char
    elseif headof(x) === :call && length(x.args) > 2 && CSTParser.is_colon(x.args[1])
        if headof(x.args[2]) === :INTEGER && headof(x.args[3]) === :INTEGER
            return CoreTypes.Int
        elseif headof(x.args[2]) === :FLOAT && headof(x.args[3]) === :FLOAT
            return CoreTypes.Float64
        elseif headof(x.args[2]) === :CHAR && headof(x.args[3]) === :CHAR
            return CoreTypes.Char
        end
    elseif hasbinding(x) && isdeclaration(x) && length(x.args) == 2
        return maybe_get_vec_eltype(x.args[2])
    end
end

function maybe_get_vec_eltype(t)
    if iscurly(t)
        lhs_ref = refof_maybe_getfield(t.args[1])
        if lhs_ref isa SymbolServer.DataTypeStore && CoreTypes.isarray(lhs_ref) && length(t.args) > 1
            refof(t.args[2])
        end
    end
end


text:
# Autogenerated wrapper script for Notcurses_jll for x86_64-apple-darwin
export libnotcurses, libnotcurses_core, libnotcurses_ffi

using Ncurses_jll
using libunistring_jll
using libdeflate_jll
JLLWrappers.@generate_wrapper_header("Notcurses")
JLLWrappers.@declare_library_product(libnotcurses, "@rpath/libnotcurses.3.dylib")
JLLWrappers.@declare_library_product(libnotcurses_core, "@rpath/libnotcurses-core.3.dylib")
JLLWrappers.@declare_library_product(libnotcurses_ffi, "@rpath/libnotcurses-ffi.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Ncurses_jll, libunistring_jll, libdeflate_jll)
    JLLWrappers.@init_library_product(
        libnotcurses,
        "lib/libnotcurses.3.0.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnotcurses_core,
        "lib/libnotcurses-core.3.0.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnotcurses_ffi,
        "lib/libnotcurses-ffi.3.0.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ConstrainedDynamics
using ConstrainedControl
using LinearAlgebra
using Rotations

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0

pll = [0.0;0.0;length1/2] # joint connection point
pul = [0.0;0.0;length1/4]
pp = [0.0;0.0;length1/4*sqrt(2)]

# Links
origin = Origin{Float64}()
lowerlegl = Box(0.1, 0.1, length1, length1)
lowerlegr = deepcopy(lowerlegl)
upperlegl = Box(0.1, 0.1, length1/2, length1/2)
upperlegr = deepcopy(upperlegl)
platform = Box(0.1, 0.1, length1/2*sqrt(2), length1/2*sqrt(2))

# Constraints
floorlr = EqualityConstraint(Revolute(origin, lowerlegl, joint_axis; p2=-pll),Revolute(origin, lowerlegr, joint_axis; p2=-pll),FixedOrientation(origin,platform;qoffset = UnitQuaternion(RotX(pi/2))))
kneel = EqualityConstraint(Revolute(lowerlegl, upperlegl, joint_axis; p1=pll, p2=-pul))
kneer = EqualityConstraint(Revolute(lowerlegr, upperlegr, joint_axis; p1=pll, p2=-pul))
platl = EqualityConstraint(Revolute(platform, upperlegl, joint_axis; p2=pul, p1=pp))
platr = EqualityConstraint(Revolute(platform, upperlegr, joint_axis; p2=pul, p1=-pp))

links = [lowerlegl;lowerlegr;upperlegl;upperlegr;platform]
constraints = [platl;platr;floorlr;kneel;kneer]

mech = Mechanism(origin, links, constraints, g=-9.81, Δt = 0.01)
setPosition!(origin,lowerlegl,p2 = -pll,Δq = UnitQuaternion(RotX(pi/4)))
setPosition!(origin,lowerlegr,p2 = -pll,Δq = UnitQuaternion(RotX(-pi/4)))
setPosition!(lowerlegl,upperlegl,p1 = pll, p2 = -pul,Δq = UnitQuaternion(RotX(-pi/2)))
setPosition!(lowerlegr,upperlegr,p1 = pll, p2 = -pul,Δq = UnitQuaternion(RotX(pi/2)))
setPosition!(upperlegl,platform,p1 = pul, p2 = pp,Δq = UnitQuaternion(RotX(3pi/4)))

xd = [links[i].state.xc for i=1:5]
qd=[links[i].state.qc for i=1:5]

Q = [diagm(ones(12))*0.0 for i=1:5]
Q[5][2,2]=10.0
Q[5][3,3]=10.0
Q[5][5,5]=1.0
Q[5][6,6]=1.0
R = [ones(1,1)*0.1 for i=1:2]

lqr = LQR(mech, getid.(links), getid.(constraints[1:2]), Q, R, Inf, xd=xd, qd=qd, Fτd = [[[6.7879484]];[[-6.7879484]]])
@test true


text:
#gridtuneの評価関数の定義

function Evalfun(Ytraining,Xtraining,CV_n)
  Ytraining=Ytraining
  Xtraining=Xtraining

  return function evalfun(scores)
    if(typeof(Ytraining) == Array{Float64,1}||typeof(Ytraining) == Array{Float64,2})
      return CV_n/sum(scores)
    else
      return mean_and_std(scores)
    end
  end
end

max_stars_repo_name: wakakusa/RFCVGT

max_issues_repo_name: wakakusa/RFCVGT

max_forks_repo_name: wakakusa/RFCVGT