text:
export get_JuMP_cartesian_model, get_JuMP_solution

"""
    get_JuMP_cartesian_model(problem_poly::Problem, mysolver)

Return JuMP cartesian model associated to `mysolver` defined by `problem_poly`, a polynomial problem in complex numbers

# Example
```jldoctest
V1 = Variable("VOLT_1",Complex)
V2 = Variable("VOLT_2",Complex)
p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
p_ctr1 = abs2(V1)
p_ctr2 = 3im * V1*conj(V2) + abs2(V1)

problem_poly=Problem()
add_variable!(problem_poly,V1)
add_variable!(problem_poly,V2)
add_constraint!(problem_poly, "ctr1", 0.95^2 << p_ctr1 << 1.05^2)
add_constraint!(problem_poly, "ctr2", p_ctr2==0)
set_objective!(problem_poly, p_obj)
print(problem_poly)
mysolver = ClpSolver()
m = get_JuMP_cartesian_model(problem_poly, mysolver)
print(m)
```
"""
function get_JuMP_cartesian_model(problem_poly::Problem, mysolver)
    pb_poly_real = problem_poly
    m = Model(solver = mysolver)
    variables_jump = SortedDict{String, JuMP.Variable}()

    ## Define JuMP variables
    startReal = 1.1
    startBool = 0
    for (varname, vartype) in pb_poly_real.variables
        if vartype==Real
            variables_jump["$varname"] = @variable(m, basename="$varname", start=startReal)
        elseif vartype==Bool
            variables_jump["$varname"] = @variable(m, category=:Bin, basename="$varname", start=startBool)
        else
            error(LOGGER, "$varname must be of type Real or Bool, here type is $vartype")
        end
    end

    ## Define JuMP constraints
    ctr_jump = SortedDict{String,JuMP.ConstraintRef}()
    ctr_exp = SortedDict{String,Tuple{Any,Float64,Float64}}()
    for (ctrname, modeler_ctr) in pb_poly_real.constraints
        polynome = modeler_ctr.p

        modeler_ctr.lb == -Inf-Inf*im || imag(modeler_ctr.lb) == 0 || error(LOGGER, "get_JuMP_cartesian_model(): Found constraint with complex bound.\nOnly real models are valid at this point.")
        modeler_ctr.ub == +Inf+Inf*im || imag(modeler_ctr.ub) == 0 || error(LOGGER, "get_JuMP_cartesian_model(): Found constraint with complex bound.\nOnly real models are valid at this point.")
        lb = real(modeler_ctr.lb)
        ub = real(modeler_ctr.ub)

        precond = modeler_ctr.precond
        for value in values(polynome)
            if imag(value)!=0
                error(LOGGER, "Polynom coefficients have to be real numbers")
            end
        end
        my_timer = @elapsed s_ctr, ispolylinear = poly_to_NLexpression(m, variables_jump,polynome)
        ctr_exp[ctrname] = (s_ctr, lb, ub)
        # @printf("%-35s%10.6f s\n", "poly_to_NLexpression for $ctr", my_timer)
        if ispolylinear
            @constraint(m, lb <= s_ctr <= ub)
        else
            if precond == :sqrt
                if lb == -Inf && ub >0
                    ##TODO: sqrt constraints for Smax
                    ctr_jump[ctrname] = @NLconstraint(m, -Inf <= s_ctr <= ub)
                else
                    error(LOGGER, "sqrt non applicable to $lb <= $s_ctr <= $ub")
                end
            else
                ctr_jump[ctrname] = @NLconstraint(m, lb <= s_ctr <= ub)
            end
        end
    end

    ## Define JuMP objective
    polynome_obj = pb_poly_real.objective
    my_timer = @elapsed s_obj, ispolylinear = poly_to_NLexpression(m, variables_jump,polynome_obj)
    # @printf("%-35s%10.6f s\n", "poly_to_NLexpression for objective", my_timer)
    if ispolylinear
        @objective(m, Min, s_obj)
    else
        @NLobjective(m,Min,s_obj)
    end
    return m, variables_jump, ctr_jump, ctr_exp
end

function poly_to_NLexpression(m::JuMP.Model, variables_jump::SortedDict{String, JuMP.Variable},polynome::Polynomial)
    s = 0
    ispolylinear = true
    d = compute_degree(polynome) #max degree on one variable, not total degree
    for expo in keys(polynome.poly)
        explsum, conjsum = get_sumdegs(expo)
        d.explvar = max(d.explvar, explsum)
        d.conjvar = max(d.conjvar,conjsum)
    end
    if d.explvar + d.conjvar > 1
        ispolylinear = false
    end
    for (monome,coeff) in polynome.poly
        prod = 1
        for (varname, degree) in monome.expo
            if !ispolylinear
                if degree.explvar > 1
                prod = @NLexpression(m, prod * variables_jump["$varname"]^degree.explvar)
                elseif degree.explvar == 1
                prod = @NLexpression(m, prod * variables_jump["$varname"])
                end
            else
                if degree.explvar != 1
                    error(LOGGER, "Exponent is supposed to be degree 1 since polynome is linear. ")
                end
                prod = @expression(m, variables_jump["$varname"])
            end
        end
        if ispolylinear
            s = @expression(m, s + coeff * prod)
        else
            s = @NLexpression(m, s + coeff * prod)
        end
    end
    return s, ispolylinear
end

"""
    sol = get_JuMP_solution(m, variables_jump, rosenbrock)

    Return a `Point` storing the candidate point in `m`.
"""
function get_JuMP_solution(m, variables_jump, pb::Problem)
    sol = Point()
    for (varname, vartype) in pb.variables
        sol[MathProgComplex.Variable(varname, vartype)] = convert(vartype, getvalue(variables_jump[varname]))
    end
    return sol
end

###TEST
#
# ###########################################
# ###polynomial problem with complex variables
# V1 = Variable("VOLT_1",Complex)
# p_obj = 0.5*(V1+conj(V1))-0.5*im*(V1-conj(V1))
# p_ctr1 = abs2(V1)
# problem_poly=Problem()
# add_variable!(problem_poly,V1)
# add_constraint!(problem_poly, "ctr1", 1 << p_ctr1 << 1 )
# set_objective!(problem_poly, p_obj)
# print(problem_poly)
# ###########################################
# pb_poly_real = pb_cplx2real(problem_poly)
# println(pb_poly_real)
#
# mysolver = KnitroSolver(KTR_PARAM_MIP_INTVAR_STRATEGY=0, KTR_PARAM_OUTLEV=4, #=KTR_PARAM_HESSOPT=2=#)
# m = get_JuMP_cartesian_model(problem_poly, mysolver)
# print(m)
# solve(m)
#
#
# #################################################
# ## the same problem with JuMP
# mj = Model(solver=mysolver)
# @variable(mj, v1_re)
# @variable(mj, v1_im)
# @objective(mj, Min, v1_re+v1_im)
# @NLconstraint(mj, 1<= 1*v1_re^2 + 1*v1_im^2 <=1)
# print(mj)
# solve(mj)
# #################################################
#
# ###########################################
# ###polynomial problem with real variables
# V1_Re = Variable("VOLT_1_Re",Real)
# V1_Im = Variable("VOLT_1_Im",Real)
# p_obj = V1_Re + V1_Im
# p_ctr1 = V1_Re^2 + V1_Im^2
# problem_poly=Problem()
# add_variable!(problem_poly,V1_Re)
# add_variable!(problem_poly,V1_Im)
# add_constraint!(problem_poly, "ctr1", 1 << p_ctr1 << 1 )
# set_objective!(problem_poly, p_obj)
# print(problem_poly)
# mysolver = KnitroSolver(KTR_PARAM_MIP_INTVAR_STRATEGY=0, KTR_PARAM_OUTLEV=4)
# m = get_JuMP_cartesian_model(problem_poly, mysolver)
# print(m)
# solve(m)
#
# ###########################################
#

"""
    get_JuMP_polar_model(pb::Problem, mysolver)

Return JuMP polar model associated to `mysolver` defined by `pb`, a polynomial problem in complex numbers

# Example
```jldoctest
V1 = Variable("VOLT_1",Complex)
V2 = Variable("VOLT_2",Complex)
p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
p_ctr1 = abs2(V1)
p_ctr2 = Polynomial((3+4im) * V1*conj(V2))
pb=Problem()
add_variable!(pb,V1)
add_variable!(pb,V2)
add_constraint!(pb, "ctr1", 0 << p_ctr1 << 100)
# add_constraint!(pb, "ctr2", p_ctr2==3+4im)
set_objective!(pb, p_obj)
print(pb)
mysolver = IpoptSolver()
m = get_JuMP_polar_model(pb, mysolver)
print(m)
```
"""
function get_JuMP_polar_model(pb::Problem, mysolver)
    m = Model(solver = mysolver)
    jump_vars = SortedDict{String, JuMP.Variable}()
    for (varname, vartype) in pb.variables
        mod = "ρ_$varname"
        jump_vars["ρ_$varname"] = @variable(m, basename="ρ_$varname")
        jump_vars["θ_$varname"] = @variable(m, basename="θ_$varname")
    end
    ctr_jump = SortedDict{String,JuMP.ConstraintRef}()
    for (ctr, modeler_ctr) in pb.constraints
        p = modeler_ctr.p
        lb = modeler_ctr.lb
        ub = modeler_ctr.ub
        ps = SortedDict{Exponent, Tuple{Real, Real}}(exp=>(real(coeff), imag(coeff)) for (exp, coeff) in p)
        ## Real part of constraint
        @NLconstraint(m, real(lb) <= sum( coeff[1] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * cos(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        - coeff[2] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * sin(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        for (mon, coeff) in ps) <= real(ub))

        ## Imag part of constraint
        @NLconstraint(m, imag(lb) <= sum( coeff[2] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * cos(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        + coeff[1] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * sin(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        for (mon, coeff) in ps) <= imag(ub))

    polynome_obj = pb.objective
    ps = SortedDict{Exponent, Tuple{Real, Real}}(exp=>(real(coeff), imag(coeff)) for (exp, coeff) in polynome_obj)
    lb_real, ub_real = real(lb), real(ub)
    @NLobjective(m, :Min, lb_real <= sum( coeff[1] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * cos(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        - coeff[2] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * sin(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        for (mon, coeff) in ps) <= ub_real)
    end
    return m
end

####TEST
# V1 = Variable("VOLT_1",Complex)
# V2 = Variable("VOLT_2",Complex)
# p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
# p_ctr1 = abs2(V1)
# p_ctr2 = Polynomial((3+4im) * V1*conj(V2))
# pb=Problem()
# add_variable!(pb,V1)
# add_variable!(pb,V2)
# add_constraint!(pb, "ctr1", 0 << p_ctr1 << 100)
# # add_constraint!(pb, "ctr2", p_ctr2==3+4im)
# set_objective!(pb, p_obj)
# print(pb)
# mysolver = IpoptSolver()
# m = get_JuMP_polar_model(pb, mysolver)
# print(m)


text:
# Default implementations of basic operators

# cpdf and logcpdf have default operators in terms of the other. Sfuncs should implement one or the other.

@impl begin
    struct SFuncCpdf end
    function cpdf(sf::SFunc{I,O}, i::I, o::O)::AbstractFloat where {I,O}
        exp(logcpdf(sf, i, o))
    end
end

@impl begin
    struct SFuncLogcpdf end
    function logcpdf(sf::SFunc{I,O}, i::I, o::O)::AbstractFloat where {I,O}
        log(cpdf(sf, i, o))
    end
end

ext: jl


text:
module MaxHelpingHandFloatySpaceBlockWithAttachedSidewaysJumpthruSupport

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/FloatySpaceBlockWithAttachedSidewaysJumpthruSupport" FloatySpaceBlockWithAttachedSidewaysJumpthruSupport(x::Integer, y::Integer,
    width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, tiletype::String="3", disableSpawnOffset::Bool=false)

const placements = Ahorn.PlacementDict(
    "Floaty Space Block (supporting Sideways Jump Throughs)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        FloatySpaceBlockWithAttachedSidewaysJumpthruSupport,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    )
)

Ahorn.editingOptions(entity::FloatySpaceBlockWithAttachedSidewaysJumpthruSupport) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::FloatySpaceBlockWithAttachedSidewaysJumpthruSupport) = 8, 8
Ahorn.resizable(entity::FloatySpaceBlockWithAttachedSidewaysJumpthruSupport) = true, true

Ahorn.selection(entity::FloatySpaceBlockWithAttachedSidewaysJumpthruSupport) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FloatySpaceBlockWithAttachedSidewaysJumpthruSupport, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# Concurrent Call Macro
# =====================

"""
    @cc expr

Take an expression `expr` and override function calls in the concurrent style.
"""
macro cc(ex)
    if ex isa Expr
        return cc_transform(ex)
    else
        return ex
    end
end

function cc_transform(ex::Expr)
    scheduler = gensym(:scheduler)
    quote
        $(scheduler) = $(Scheduler)()
        try
            $(Base.map)($(Base.fetch), $(rewrite_funccalls(ex, scheduler)))
        catch
            $(Base.rethrow)()
        finally
            $(teardown)($(scheduler))
            # println(sprint(show, $(scheduler).stats))
            # println(STDOUT, $(scheduler).trace)
        end
    end
end

function rewrite_funccalls(ex::ANY, scheduler::Symbol)
    if !(ex isa Expr)
        return ex isa Symbol ? esc(ex) : ex
    end
    rewrite(ex) = rewrite_funccalls(ex, scheduler)
    head = ex.head
    args = ex.args
    if head == :call  # function call
        if length(args) > 0 && args[1] == :pmap
            return Expr(:call, ConcurrentCalls.pmap, scheduler, map(rewrite, args[2:end])...)
        else
            return Expr(:call, ConcurrentCalls.call, scheduler, map(rewrite, args)...)
        end
    elseif head == :generator
        @assert length(args) == 2
        @assert args[2].head == :(=)
        iter = args[2].args[2]
        closure = Expr(:->, esc(args[2].args[1]), esc(args[1]))
        return Expr(:call, ConcurrentCalls.ChunkingGenerator, scheduler, rewrite(iter), closure)
    elseif head == :comprehension
        @assert length(args) == 1
        return Expr(:call, collect, rewrite(args[1]))
    elseif head == :->  # closure
        return esc(ex)
    else
        return Expr(head, map(rewrite, args)...)
    end
end


text:
module TestModelStructure

# tests the framework of components and connections

using Test
using Mimi

import Mimi:
    connect_param!, unconnected_params, set_dimension!,  build,
    get_connections, internal_param_conns, dim_count,  dim_names,
    modeldef, modelinstance, compdef, getproperty, setproperty!, dimension, compdefs

@defcomp A begin
    varA = Variable{Int}(index=[time])
    parA = Parameter{Int}()

    function run_timestep(p, v, d, t)
        v.varA[t] = p.parA
    end
end

@defcomp B begin
    varB = Variable{Int}()

    function run_timestep(p, v, d, t)
        if t.t < 10
            v.varB = 1
        else
            v.varB = 10
        end
    end
end

@defcomp C begin
    varC = Variable{Int}()
    parC = Parameter{Int}()

    function run_timestep(p, v, d, t)
        v.varC = p.parC
    end
end

m = Model()

# TBD: This is not necessarily an error with composites.
# make sure you can't add a component before setting time dimension
# @test_throws ErrorException add_comp!(m, A)

set_dimension!(m, :time, 2015:5:2100)

add_comp!(m, A)
add_comp!(m, B, before=:A)
add_comp!(m, C, after=:B)
# Component order is B -> C -> A.

connect_param!(m, :A, :parA, :C, :varC)

unconns = unconnected_params(m)
@test length(unconns) == 1
c = compdef(m, :C)
uconn = unconns[1]
@test uconn.comp_name == :C
@test uconn.datum_name == :parC

connect_param!(m, :C => :parC, :B => :varB)

@test_throws ErrorException add_comp!(m, C, after=:A, before=:B)

@test length(m.md) == 3

@test length(internal_param_conns(m)) == 2
c = compdef(m, :C)
@test get_connections(m, :A, :incoming)[1].src_comp_path == c.comp_path

@test length(get_connections(m, :B, :incoming)) == 0
@test get_connections(m, :B, :outgoing)[1].dst_comp_path == c.comp_path

@test length(get_connections(m, :A, :all)) == 1

@test length(unconnected_params(m)) == 0

run(m)

#############################################
#  Tests for model def and instance         #
#############################################

@test modeldef(m) == m.md
@test modelinstance(m) == m.mi

#############################################
#  Tests for connecting scalar parameters   #
#############################################

@test all([m[:A, :varA][t] == 1 for t in 1:9])

@test all([m[:A, :varA][t] == 10 for t in 10:dim_count(m.md, :time)])

##########################
#   tests for indexing   #
##########################

@test dim_count(m.md, :time) == 18

@test m[:A, :parA] == 10
@test_throws ErrorException m[:A, :xx]

time = dimension(m, :time)
a = collect(keys(time))
@test all([a[i] == 2010 + 5*i for i in 1:18])

@test dim_names(m, :A, :varA)[1] == :time
@test length(dim_names(m, :A, :parA)) == 0

################################
#  tests for delete! function  #
################################

@test_throws ErrorException delete!(m, :D)
@test length(internal_param_conns(m.md)) == 2
delete!(m, :A)
@test length(internal_param_conns(m.md)) == 1
@test !(:A in compdefs(m))
@test length(compdefs(m)) == 2

#######################################
#   Test check for unset parameters   #
#######################################

@defcomp D begin
  varD = Variable(index=[time])
  parD = Parameter()
end

add_comp!(m, D)
@test_throws ErrorException build(m)

##########################################
#   Test init function                   #
##########################################

@defcomp E begin
    varE = Variable{Int}()
    parE1 = Parameter{Int}()
    parE2 = Parameter{Int}()

    function init(p, v, d)
        v.varE = p.parE1
    end

    function run_timestep(p, v, d, t)
        if !is_first(t)
            v.varE = p.parE2
        end
    end
end

m = Model()

# run for several timesteps
set_dimension!(m, :time, 2015:5:2100)

add_comp!(m, E)
set_param!(m, :E, :parE1, 1)
set_param!(m, :E, :parE2, 10)

run(m)
@test m[:E, :varE] == 10

# run for just one timestep, so init sets the value here
set_dimension!(m, :time, [2015])

run(m)
@test m[:E, :varE] == 1

end # module


text:
module DoubleDummyImperfectInfo
  export Game, Board
  include("game.jl")
  module Training
    using AlphaZero
    include("params.jl")
  end
  include("solver.jl")
end

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
# ---
# title: 145. Binary Tree Postorder Traversal
# id: problem145
# author: zhwang
# date: 2022-02-19
# difficulty: Medium
# categories: Stack, Tree
# link: <https://leetcode.com/problems/binary-tree-postorder-traversal/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _the postorder traversal of its
# nodes ' values_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/pre1.jpg)
#
#
#
#     Input: root = [1,null,2,3]
#     Output: [3,2,1]
#
#
# **Example 2:**
#
#
#
#     Input: root = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#
#
# **Example 4:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/pre3.jpg)
#
#
#
#     Input: root = [1,2]
#     Output: [2,1]
#
#
# **Example 5:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/pre2.jpg)
#
#
#
#     Input: root = [1,null,2]
#     Output: [2,1]
#
#
#
#
# **Constraints:**
#
#   * The number of the nodes in the tree is in the range `[0, 100]`.
#   * `-100 <= Node.val <= 100`
#
#
#
# **Follow up:**
#
# Recursive solution is trivial, could you do it iteratively?
#
#
#
#
## @lc code=start
using LeetCode
postorder_traversal(::Nothing) = Int[]
function postorder_traversal(root::TreeNode)::Vector{Int}
    res, stack = Int[], [(root, -1)] ## -1 for left subtree, 1 for right subtree
    while !isempty(stack)
        cur = last(stack)
        if last(cur) == -1 ## search left subtree
            stack[end] = (cur[1], 1)
            !isnothing(cur[1].left) && push!(stack, (cur[1].left, -1))
        elseif last(cur) == 1 ## search right subtree
            stack[end] = (cur[1], 0)
            !isnothing(cur[1].right) && push!(stack, (cur[1].right, -1))
        else
            cur = pop!(stack)
            push!(res, cur[1].val) ## postorder traversal: put codes here
        end
    end
    return res
end

##### code template for inorder/preorder/postorder traversal #####
## function traversal(root::TreeNode)::Vector{Int}
##     res, stack = Int[], [(root, -1)] ## -1 for left subtree, 1 for right subtree
##     while !isempty(stack)
##         cur = last(stack)
##         if last(cur) == -1 ## search left subtree
##             stack[end] = (cur[1], 1)
##             ## preorder traversal: put codes here
##             !isnothing(cur[1].left) && push!(stack, (cur[1].left, -1))
##         elseif last(cur) == 1 ## search right subtree
##             stack[end] = (cur[1], 0)
##             ## inorder traversal: put codes here
##             !isnothing(cur[1].right) && push!(stack, (cur[1].right, -1))
##         else
##             cur = pop!(stack)
##             ## postorder traversal: put codes here
##         end
##     end
##     return res
## end

## @lc code=end


text:
cell1 = cell()
cost1 = cost_default()
cell1_converted,cost1_converted = convert_all(cell1,cost1,mult)

cost1.general_costs.no_units_mfg = converter(80.0, mult.units_mfg)
cost_kwh,thing = cost_calc(cell1,cost1)
println("Cost/kWH: ",cost_kwh)
breakdown_gen = cost_calc(cell1,cost1,breakdown=true)

(total_unit_cost_all_materials  ,
total_direct_labor_per_unit           ,
variable_overhead_per_unit            ,
GSA_per_unit                          ,
R_and_D_per_unit                      ,
depreciation_per_unit                 ,
profits_after_taxes_per_unit          ,
warrant_cost_per_unit                 ,
cell_energy_storage) = breakdown_gen

breakdown_vec = breakdown_gen[1:end-1]
total_cost = sum(breakdown_vec)
breakdown_percentages = breakdown_vec./total_cost
breakdown_big_categories = zeros(4)
breakdown_big_categories[1] = total_unit_cost_all_materials #Materials
breakdown_big_categories[2] = total_direct_labor_per_unit #Direct Labor
breakdown_big_categories[3] = variable_overhead_per_unit+GSA_per_unit+R_and_D_per_unit+profits_after_taxes_per_unit+warrant_cost_per_unit
breakdown_big_categories[4] = depreciation_per_unit
breakdown_big_categories = breakdown_big_categories./total_cost

mat_breakdown= cost_calc(cell1,cost1,materials_breakdown = true)
cathode_cost,anode_cost,electrolyte_cost,others_cost = mat_breakdown
mat_normalized_breakdown = mat_breakdown./sum(mat_breakdown)

using PyPlot

figure(1)
clf()
pie(breakdown_big_categories,autopct="%1.1f%%")
legend(["Materials","Labor","Overhead","Depreciation"])

figure(2)
clf()
pie(mat_normalized_breakdown,autopct="%1.1f%%")
legend(["Cathode","Anode","Electrolyte","Other"])

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
module CartPendulum

global const g = 9.82

using DifferentialEquations
using ControlSystems
mutable struct Pendulum{T<:Real}
    M::T
    m::T
    l::T
    b::T
    ulims::Tuple{T,T}
    nx::T
    nθ::T
    sgn::T
end

function Pendulum()
    Pendulum{Float64}(1., 0.2, 0.6, 0.1, (-10., 10.), 0.1, 0.01,1)
end

function I(pen::Pendulum{T}) where T<:Real
    return pen.m*pen.l^2/3
end

function linearize(pen::Pendulum{<:Real})
    (M, m, l, b) = (pen.M, pen.m, pen.l, pen.b)
    denom = I(pen) * (M + m) + M * m * l^2
    a22 = -(I(pen) + m * l^2) * b / denom
    a23 = m^2 * g * l^2 / denom;
    a42 = -m * l * b / denom;
    a43 = m * g * l * (M + m) / denom;
    b2 = (I(pen) + m * l^2) / denom;
    b4 = m * l / denom;

    A = [0 1 0 0
    0 a22 a23 0
    0 0 0 1
    0 a42 a43 0]

    B = [0;b2;0;b4];
    # For now, assume full state feedback
    # C = diagm(ones(4));
    C = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
    return ss(A, B, C, 0)
end

# Deterministic part
function f!(dx::AbstractVector{T},
    x::AbstractVector{T},
    pen::Pendulum{T}, t::T) where T<:Real
    (M, m, l, b, ulims) = (pen.M, pen.m, pen.l, pen.b, pen.ulims)

    lump1 = m^2 * l^2 / (I(pen) + m * l^2)
    c = cos(x[3])
    s = sin(x[3])

    dx[1] = x[2];
    dx[2] = (-b * x[2] + g * lump1 * c * s + m * l * s * x[4]^2 + pen.sgn*sat(x[5], ulims)) / (M + m - lump1 * c^2)
    dx[3] = x[4];
    dx[4] = (-m * l * c * dx[2] - m * g * l * s) / (I(pen) + m * l^2);
    dx[5:end] .= 0; # Controller state is set by callback and held constant otherwise
end

# Stochastic noise
function g!(dx::AbstractVector{T},
    x::AbstractVector{T},
    pen::Pendulum{T}, t::T) where T<:Real

    (M, m, l, b) = (pen.M, pen.m, pen.l, pen.b)
    (nx, nθ) = (pen.nx, pen.nθ)
    lump1 = m^2 * l^2 / (I(pen) + m * l^2)
    c = cos(x[3])
    s = sin(x[3])

    dx[1] = 0;
    dx[2] = pen.nx/ (M + m - lump1 * c^2)
    dx[3] = 0;
    dx[4] = nθ/(I(pen) + m * l^2)
    dx[5:end] .= 0
end

function sat(val::T, lims::Tuple{T,T}) where T<:Real
    return min(max(val, lims[1]), lims[2])
end
end


text:
using DifferentialEquations
using LinearAlgebra
using Plots

struct Wave
    f
    g
    c
end

# Exercise

function solve_wave(T, L, wave::Wave; n_t=100, n_x=100)
    ts = range(0, T; length=n_t)
    xs = range(0, L; length=n_x)
    Δt = ts[2] - ts[1]
    Δx = xs[2] - xs[1]
    y = zeros(n_t, n_x)

    y[:,1] .= wave.f(0)
    y[:,end] .= wave.f(L)

    y[1,2:end-1] = wave.f.(xs[2:end-1])
    y[2,2:end-1] = y[1,2:end-1] + Δt*wave.g.(xs[2:end-1])

    for t in 2:n_t-1, x in 2:n_x-1
        ∂y_xx = (y[t, x+1] - 2*y[t, x] + y[t, x-1])/Δx^2
        y[t+1, x] = c^2 * Δt^2 * ∂y_xx  + 2*y[t, x] - y[t-1, x]
    end

    return y
end

# Visualization

function plot_wave(y, file_name; fps = 60, kwargs...)
    anim = @animate for (i, y_row) in enumerate(eachrow(y))
        plot(
            y_row;
            title = "t = $(i-1)Δt",
            xlabel = "x",
            ylabel = "y(t, x)",
            legend = false,
            linewidth = 2,
            kwargs...
        )
    end
    gif(anim, file_name; fps, show_msg = false)

    return nothing
end

# Exercise

f(x,L) = 2*exp(-(x-L/2)^2) + x/L
g(x) = 0

L = 1.5*pi
T = 240
c = 0.02

wave = Wave(x -> f(x,L), g, c)

y1 = solve_wave(T, L, wave; n_t=241, n_x=101)
plot_wave(y1, "wave1.gif"; ylims=(-2,3), label="")

y2 = solve_wave(T, L, wave; n_t=241, n_x=7)
plot_wave(y2, "wave2.gif"; ylims=(-2,3), label="")

# DiferentialEquations

f(u,p,t) = 0.98*u

u0 = 1.0
tspan = (0.0, 1.0)

prob = ODEProblem(f, u0, tspan)

sol = solve(prob)

plot(sol; label="")

sol(0.8)

# Exercise

sol2 = solve(prob, dense=false)
sol3 = solve(prob, RK4())

ts = range(tspan...; length = 100)

plot(ts, t->exp(0.98*t), label="True solution", legend=:topleft)
plot!(ts, t->sol(t), label="Default")
plot!(ts, t->sol2(t), label="Linear")
plot!(ts, t->sol3(t), label="Runge-Kutta")

# Lorenz system

function lorenz(u, p, t)
    σ, ρ, β = p
    x_t = σ*(u[2]-u[1])
    y_t = u[1]*(ρ-u[3]) - u[2]
    z_t = u[1]*u[2] - β*u[3]
    return [x_t; y_t; z_t]
end

u0 = [1.0; 0.0; 0.0]
p = [10; 28; 8/3]

tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan, p)

sol = solve(prob)

plot(sol)

plt1 = plot(sol, vars=(1,2,3), label="")

plot(sol, vars=(1,2,3), denseplot=false; label="")

traj = hcat(sol.u...)
plot(traj[1,:], traj[2,:], traj[3,:]; label="")

# Exercise

p0 = (nextfloat(p[1]), p[2:end]...)

prob0 = ODEProblem(lorenz, u0, tspan, p0)
sol0 = solve(prob0)

plt0 = plot(sol0, vars=(1,2,3), label="")

plot(plt1, plt0; layout=(1,2))

# Solution comparison

hcat(sol(tspan[2]), sol0(tspan[2]))

# Bisection setup

function bisection(f, a, b; tol=1e-6)
    fa = f(a)
    fb = f(b)
    fa == 0 && return a
    fb == 0 && return b
    fa*fb > 0 && error("Wrong initial values for bisection")
    while b-a > tol
        c = (a+b)/2
        fc = f(c)
        fc == 0 && return c
        if fa*fc > 0
            a = c
            fa = fc
        else
            b = c
            fb = fc
        end
    end
    return (a+b)/2
end

# PMSM structure

struct PMSM{T<:Real}
    ρ::T
    ω::T
    A::Matrix{T}
    invA::Matrix{T}

    function PMSM(ρ, ω)
        A = -ρ*[1 0; 0 1] -ω*[0 -1; 1 0]
        return new{eltype(A)}(ρ, ω, A, inv(A))
    end
end

function expA(p::PMSM, t)
    ρ, ω = p.ρ, p.ω
    return exp(-ρ*t)*[cos(ω*t) sin(ω*t); -sin(ω*t) cos(ω*t)]
end

ρ = 0.1
ω = 2
x0 = [0;-0.5]
q = [1;0]

ps = PMSM(ρ, ω)

# Exercise

t = 5
exp0 = exp.(t*ps.A)
exp1 = exp(t*ps.A)
exp2 = expA(ps, t)

norm(exp1 - exp0) >= 1e-10 || error("Matrices are wrong")
norm(exp1 - exp2) <= 1e-10 || error("Matrices are wrong")

# Exercise

function trajectory_fin_diff(p::PMSM, x0, ts, q)
    xs = zeros(length(x0), length(ts))
    xs[:, 1] = x0

    for i in 1:length(ts)-1
        xs[:, i+1] = xs[:, i] + (ts[i+1]-ts[i])*(p.A * xs[:, i] + q)
    end
    return xs
end

function trajectory_exact(p::PMSM, x0, ts, q)
    xs = zeros(length(x0), length(ts))

    for (i, t) in enumerate(ts)
        xs[:, i] = expA(p, t)*(x0 + p.invA * (I - expA(p, -t))*q)
    end
    return xs
end

ts = 0:0.01:10

xs1 = trajectory_fin_diff(ps, x0, ts, q)
xs2 = trajectory_exact(ps, x0, ts, q)

plot(xs1[1,:], xs1[2,:], label="Finite differences")
plot!(xs2[1,:], xs2[2,:], label="True value")

# Plotting trajectory

ts = 0:0.0001:10

xs1 = trajectory_fin_diff(ps, x0, ts, q)
xs2 = trajectory_exact(ps, x0, ts, q)

plot(xs1[1,:], xs1[2,:], label="Finite differences")
plot!(xs2[1,:], xs2[2,:], label="True value")

# Exercise

U_max = 0.1
x_t = [0.25;-0.5]

f(t) = norm(expA(ps, -t)*x_t - x0 - ps.invA*(I-expA(ps, -t))*q) - U_max/ps.ρ*(exp(ps.ρ*t)-1)

τ = bisection(f, minimum(ts), maximum(ts))

# Plotting trajectory

function trajectory_control(p::PMSM, x0, ts, q, U_max, p0)
    xs = zeros(length(x0), length(ts))

    for (i, t) in enumerate(ts)
        eAt  = expA(p, t)
        emAt = expA(p, -t)
        xs[:, i] = eAt*(x0 + p.invA * (I - emAt)*q + U_max/ρ*(exp(ρ*t) - 1)*p0)
    end
    return xs
end

p0 = ps.ρ/(U_max*(exp(ps.ρ*τ)-1))*(expA(ps, -τ)*x_t - x0 - ps.invA*(I-expA(ps, -τ))*q)
p0 /= norm(p0)

ts = range(0, τ; length=100)

traj = trajectory_control(ps, x0, ts, q, U_max, p0)

plot(traj[1,:], traj[2,:], label="Optimal trajectory")
scatter!([x0[1]], [x0[2]], label="Starting point")
scatter!([x_t[1]], [x_t[2]], label="Target point")

# BONUS

ts = 0:0.01:10

plt = plot()
for α = 0:π/4:2*π
    trj = trajectory_control(ps, x0, ts, q, U_max, [sin(α); cos(α)])
    plot!(plt, trj[1,:], trj[2,:], label="")
end
display(plt)

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
function get_pitch(M::AbstractEquilibrium, c::HamiltonianCoordinate, r::T, z::T) where {T<:Number}
    psi = M(r,z)
    g = poloidal_current(M,psi)
    babs = norm(Bfield(M,r,z))
    phi = electric_potential(M,psi)
    KE = c.energy - 1e-3*phi
    f = -babs/(sqrt(2e3*e0*KE*c.m)*g*B0Ip_sign(M))
    pitch = f*(c.p_phi - c.q*e0*psi)
    #pitchabs = sqrt(max(1.0-(c.mu*babs/(1e3*e0*KE)), 0.0))
    #if !isapprox(abs(pitch), pitchabs, atol=1.e-1)
    #    warn("abs(pitch) != abspitch: ",pitchabs," ",pitch)
    #end
    return clamp(pitch,-1.0,1.0)
end

function get_pitch(M::AbstractEquilibrium, gcp::GCParticle, p_para::T, mu::T, r::T, z::T) where {T<:Number}
    Babs = norm(Bfield(M,r,z))
    m = gcp.m
    p = sqrt(p_para^2 + 2*m*Babs*mu)
    pitch = B0Ip_sign(M)*p_para/p
    return pitch
end

function get_pitch(M::AbstractEquilibrium, c::T, path::OrbitPath) where {T<:AbstractOrbitCoordinate}
    return get_pitch.(M, c, path.r, path.z)
end

function get_pitch(M::AbstractEquilibrium, c::T, r::S, z::S) where {T<:AbstractOrbitCoordinate, S<:Number}
    hc = HamiltonianCoordinate(M, c)
    return get_pitch(M, hc, r, z)
end

function get_pitch(M::AbstractEquilibrium, o::Orbit)
    return get_pitch.(M, o.coordinate, o.path)
end

function get_kinetic_energy(M::AbstractEquilibrium, c::HamiltonianCoordinate, r::T, z::T) where {T<:Number}
    psi = M(r,z)
    return c.energy - 1e-3*electric_potential(M,psi)
end

function get_kinetic_energy(M::AbstractEquilibrium, gcp::GCParticle, p_para::T, mu::T, r::T, z::T) where {T<:Number}
    Babs = norm(Bfield(M,r,z))
    m = gcp.m
    mc2 = m*c0^2
    p = sqrt(p_para^2 + 2*m*Babs*mu)
    KE = 1e-3*(hypot(p*c0, mc2) - mc2)/e0 #keV
end

function get_kinetic_energy(M::AbstractEquilibrium, c::T, path::OrbitPath) where {T<:AbstractOrbitCoordinate}
    return get_kinetic_energy.(M, c, path.r, path.z)
end

function get_kinetic_energy(M::AbstractEquilibrium, c::T, r::S, z::S) where {T<:AbstractOrbitCoordinate, S<:Number}
    hc = HamiltonianCoordinate(M, c)
    return get_kinetic_energy(M, hc, r, z)
end

function get_kinetic_energy(M::AbstractEquilibrium, o::Orbit)
    path = o.path
    c = o.coordinate
    return get_kinetic_energy(M, o.coordinate, o.path)
end

function classify(r, z, pitch, axis; n=length(r))

    op = Boundary(r,z)

    sign_p = sign.(pitch)
    if in_boundary(op, axis)
        if all(sign_p .== sign_p[1])
            if sign_p[1] > 0.0
                class = :co_passing
            else
                class = :ctr_passing
            end
        else
            if sign_p[argmax(r)] >= 0.0
                class = :potato
            else
                class = :ctr_passing #really ctr_potato
            end
        end
    else
        if all(sign_p .== sign_p[1])
            class = :stagnation
        else
            class = :trapped
        end
    end

    return class
end

function classify(path::OrbitPath, axis)
    return classify(path.r,path.z,path.pitch, axis, n=length(path))
end

"""
    _eigmax(A) -> λ_max

Returns approximate largest eigval for a 3x3 matrix using power iteration
"""
function _eigmax(A::SMatrix{3,3}; N=10_000_000)
    b_k = SVector{3}(1.0,0.0,0.0)
    for i=1:N
        # calculate the matrix-by-vector product Ab
        b_k = A*b_k
        b_k = b_k/norm(b_k)
    end
    λ_k = (b_k'*A*b_k)/(b_k'*b_k)
    return λ_k
end

"""

gcde_check(M, o; verbose=false)

For each point along the orbit path, compute the gyro-radius r_g, the magnetic field
amplitude B and the square root of the maximum eigenvalue of the matrix ̂M. Use these quantities
to determine whether the usage of guiding-center drift equations (GCDE) could be expected to produce
a significant discrepancy compared to full-orbit equations (FOE). The matrix ̂M is the matrix
found in equation (2) of D. Pfefferlé et al (2015) (https://doi.org/10.1088/0741-3335/57/5/054017)
and the condition to be evaluated is the condition found in equation (3) of the same paper.
In section 3 of the paper, the abritrary threshold of 0.073 is used. The default is the same here.
If the condition is satisfied for the threshold, no significant discrepancy between GCDE
and FOE is to be expected. The equations to compute ̂M are from Appendix B of the same paper (B.2).

Return true if the criterion is fulfilled, and it is ok to use GCDE.
Return false if the criterion is not fulfilled, and FOE should be used instead.
"""
function gcde_check(M::AbstractEquilibrium, gcp::GCParticle, path::OrbitPath; threshold=0.073, verbose=false)

    g_rz = (x) -> poloidal_current(M, M(x[1],x[2])) # Poloidal current function (F=R*Bt) as a function of R,Z
    dpsi_rz = (x) -> psi_gradient(M, x[1], x[2]) # [dψ/dR,dψ/dZ] as a function of R,Z
    m = gcp.m # Mass of particle. kg
    KE = gcp.energy # Kinetic energy. keV
    mc2 = m*c0*c0 # Rest energy. Joule
    KE_j = e0*KE*1e3 # Kinetic energy. Joule
    p_rel2 = ((KE_j + mc2)^2 - mc2^2)/(c0*c0) # Square of relativistic momentum. kg^2 m^2 s^-2
    q = gcp.q*e0 # Particle charge. Coulomb

    maxcrit = 0.0 # To keep track of the maximum of criterion violation
    #Mhat = zeros(3,3) # Alternative function execution approach
    for i=1:length(path)
        R = path.r[i] # Major radius position of particle. m
        φ = path.phi[i] # Toroidal angle position of particle. -
        Z = path.z[i] # Vertical position of particle. m
        pitch = path.pitch[i] # Pitch of particle. -

        p_perp2 = p_rel2*(1-pitch^2) # Square of relativistic perpendicular momentum. kg^2 m^2 s^-2
        B = Equilibrium.Bfield(M, R, Z) # Magnetic field vector at particle position.
        Babs = norm(B) # Magnetic field magnitude. Tesla

        r_g = sqrt(p_perp2) / (abs(q)*Babs) # Gyroradius at particle position. m

        F = g_rz(SVector{2,Float64}(R,Z)) # Poloidal current function in meter*Tesla
        grad_psi = SVector{2,Float64}(psi_gradient(M,R,Z)) # [dψ/dR,dψ/dZ]. SVector for efficiency
        cc = cocos(M) # Cocos factor
        cocos_factor = cc.sigma_RpZ*cc.sigma_Bp/((2pi)^cc.exp_Bp)

        grad_F = SVector{2,Float64}(ForwardDiff.gradient(g_rz,SVector{2}(R,Z))) # [∂F/∂r,∂F/∂z] SVector for efficiency
        J_grad_psi = SMatrix{2,2}(ForwardDiff.jacobian(dpsi_rz,SVector{2}(R,Z))) # [d^2ψ/dr^2 d^2ψ/dzdr;d^2ψ/drdz d^2ψ/dz^2] SMatrix for efficiency

        D = SMatrix{3,3}(cocos_factor*inv(R)*(J_grad_psi[2,1]-grad_psi[2]*inv(R)),grad_F[1]-F*inv(R),cocos_factor*inv(R)*(-J_grad_psi[1,1]+inv(R)*grad_psi[1]),-F*inv(R),cocos_factor*grad_psi[2],0.0,cocos_factor*inv(R)*J_grad_psi[2,2],grad_F[2],-cocos_factor*inv(R)*J_grad_psi[1,2]) # Matrix D as in equation (2) in D. Pfefferlé et al (2015). Note D_{ij}=B_{i;j} in appendix B.

        G = SMatrix{3,3}(1.0,0.0,0.0,0.0,inv(R*R),0.0,0.0,0.0,1.0)

        P = G .- inv(Babs^2)*SMatrix{3,3}(grad_psi[2]*grad_psi[2]*inv(R^2),
                                          cocos_factor*grad_psi[2]*F*inv(R^3),
                                          -grad_psi[1]*grad_psi[2]*inv(R^2),
                                          cocos_factor*grad_psi[2]*F*inv(R^3),
                                          (F^2)*inv(R^4),
                                          -cocos_factor*F*grad_psi[1]*inv(R^3),
                                          -grad_psi[1]*grad_psi[2]*inv(R^2),
                                          -cocos_factor*F*grad_psi[1]*inv(R^3),
                                          grad_psi[1]*grad_psi[1]*inv(R^2))

        Λ = SMatrix{3,3}(cos(φ), sin(φ), 0.0, -R*sin(φ), R*cos(φ), 0.0, 0.0, 0.0,1.0)

        DPL = D*P*Λ

        #Mhat .= transpose(DPL)*G*DPL
        #λmax = maximum(eigvals(Mhat)) # Could use λmax = tr(Mhat) instead to speed up. But less accurate.
        λmax = _eigmax(transpose(DPL)*G*DPL,N=100)

        criterion = sqrt(λmax)*r_g / Babs

        if criterion > maxcrit
            maxcrit = criterion
        end

        if criterion > threshold # criterion violated
            if verbose
                println("--- Criterion violated! ---")
                println("Criterion > Threshold: $(criterion) > $(threshold)")
                println("- Orbit path index at bad position: $(i)")
                println("- Gyroradius (r_g) at bad position: $(r_g) m")
                println("- |B| at bad position: $(Babs) T")
                println("- √(λmax) at bad position: $(sqrt(λmax)) T/m")
                println("√(λmax) * r_g / |B| > Threshold")
            end
            return false
        end
    end
    if verbose
        println("- Criterion was not violated at any point along the orbit path.")
        println("maximum(criterion) < threshold: $(maxcrit) < $(threshold)")
        println("- gcde ok to use!")
    end
    return true
end

function gcde_check(M::AbstractEquilibrium, o::Orbit; kwargs...)
    gcde_check(M, GCParticle(o.coordinate), o.path; kwargs...)
end

function hits_wall_path(path::OrbitPath, wall::Wall)

    not_in_vessel = [~in_vessel(wall,p) for p in zip(path.r,path.z)]
    ind = findfirst(not_in_vessel)
    if ind != nothing
        return false, path
    else
        ind = ind-1
        return true, OrbitPath(path.vacuum, path.drift,
                               path.energy[1:ind], path.pitch[1:ind],
                               path.r[1:ind], path.z[1:ind], path.phi[1:ind],
                               path.dt[1:ind])
    end
end

function hits_wall(path::OrbitPath, wall::Wall)
    not_in_vessel = [~in_vessel(wall,p) for p in zip(path.r,path.z)]
    return any(not_in_vessel)
end

function lorentz_factor(K_keV,m)
    KE_j = e0*1e3*K_keV
    mc2 = m*c0^2
    p_rel2 = ((KE_j + mc2)^2 - mc2^2)/(c0*c0)
    return sqrt(1 + p_rel2/((m*c0)^2))
end

function lorentz_factor(p::GCParticle)
    return lorentz_factor(p.energy, p.m)
end

function lorentz_factor(p::Particle)
    v = hypot(p.vr,p.vt,p.vz)
    beta = v/c0
    return inv(sqrt(1 - beta^2))
end

function cyclotron_frequency(M,p::T) where T <: AbstractParticle
    gamma = lorentz_factor(p)
    Babs = norm(Bfield(M, p.r, p.z))
    return abs(p.q*e0)*Babs/(gamma*p.m)
end

function cyclotron_period(M,p::T) where T <: AbstractParticle
    return 2*pi/cyclotron_frequency(M,p)
end

function gyroradius(M,p::GCParticle)
    γ = lorentz_factor(p)
    Babs = norm(Bfield(M, p.r, p.z)) # Magnetic field magnitude. Tesla
    m = p.m # Mass of particle
    KE = p.energy # Kinetic energy. keV
    mc2 = m*c0*c0 # Rest energy. Joule
    KE_j = e0*KE*1e3 # Kinetic energy. Joule
    p_rel2 = ((KE_j + mc2)^2 - mc2^2)/(c0*c0) # Relativistic momentum
    p_perp2 = p_rel2*(1-p.pitch^2) # Square of relativistic perpendicular momentum

    return sqrt(p_perp2) / (abs(p.q)*B)
end

function gyroradius(M,p::Particle)
    B = Bfield(M, p.r, p.z)
    vpara, vperp = parallel_perpendicular([p.vr p.vt p.vz], B)
    γ = lorentz_factor(p)
    return γ*(p.m)*vperp / (abs(p.q)*norm(B))
end

function normalize(M::AbstractEquilibrium, hc::HamiltonianCoordinate)
    E = hc.energy
    B0 = norm(Bfield(M, magnetic_axis(M)...))
    mu = (abs(B0)/(E*e0*1e3))*hc.mu
    flux = abs(-(psi_limits(M)...))
    Pphi = (B0Ip_sign(M)/(e0*flux))*hc.p_phi

    return E, Pphi, mu
end

function normalize(M::AbstractEquilibrium, KE, p, r, z; amu=H2_amu, q=1)
    normalized(M,HamiltonianCoordinate(M,KE,p,r,z,amu*mass_u,q))
end

"""
    v_para, v_perp = parallel_perpendicular(v, B)

Returns parallel velocity component (`v_para`) w.r.t to B and the perpendicular vector (`v_perp`)
"""
function parallel_perpendicular(v, B)
    Bhat = B/norm(B)
    v_para = dot(v,Bhat)
    v_perp = cross(-Bhat,cross(Bhat,v))
    return v_para, v_perp
end

"""
    a, c = perpendicular_vectors(b)

Calculates normalized vectors that are perpendicular to `b` such that `a x c = b_norm`
"""
function perpendicular_vectors(B::T) where T
    Bhat = B/norm(B)

    if abs(Bhat[3]) == 1
        a = convert(T, [1, 0, 0])
        c = convert(T, [0, 1, 0])
    else
        if Bhat[3] == 0
            a = convert(T, [0, 0, 1])
            c = convert(T, [Bhat[2], -Bhat[1], 0])
        else
            a = convert(T, [Bhat[2], -Bhat[1], 0]/norm(Bhat[1:2]))
            c = convert(T, -[a[2] , -a[1] , (a[1]*Bhat[2] - a[2]*Bhat[1])/Bhat[3]])
            c = c/norm(c)
            if Bhat[3] < 0
                c = -c
            end
        end
    end
    return a, c
end

"""
    velocity(gcp::GCParticle, gamma)

Returns the velocity of the particle in cartesian coordinates. Assumes toroidal angle = 0
"""
function velocity(M::AbstractEquilibrium, gcp::GCParticle, gamma)
    mc2 = gcp.m*c0^2

    p0 = sqrt(((1e3*e0*gcp.energy + mc2)^2 - mc2^2)/(c0*c0)) # The initial particle momentum
    if abs(gcp.pitch) == 1.0
        pitch0 = sign(gcp.pitch)*prevfloat(abs(gcp.pitch))
    else
        pitch0 = gcp.pitch
    end

    p_para0 = p0*pitch0*B0Ip_sign(M) # The initial parallel momentum
    p_perp0 = p0*sqrt(1.0 - pitch0^2) # The initial perpendicular momentum

    B = Bfield(M,gcp.r,0.0,gcp.z)
    Babs = norm(B)
    b = B/Babs
    a, c = perpendicular_vectors(B)

    pvec = p_perp0*cos(gamma)*a .+ p_perp0*sin(gamma)*c .+ p_para0*b

    v = pvec/(gcp.m*lorentz_factor(gcp))
    return v
end

"""
    gyro_step(M::AbstractEquilibrium, gcp, gamma) -> r_gyro

Returns gyro-step such that the r_p = r_gc - r_gyro in cartesian coordinates. Assumes toroidal angle=0
"""
function gyro_step(M::AbstractEquilibrium, gcp::GCParticle, gamma)

    v = velocity(M, gcp, gamma)

    B = Bfield(M,gcp.r,0.0,gcp.z)
    Babs = norm(B)
    b = B/Babs

    # First order
    Ω_c = cyclotron_frequency(M,gcp)
    r_gyro = cross(v,b)/Ω_c

    # Second order
    # Belova, E. V., N. N. Gorelenkov, and C. Z. Cheng. "Self-consistent equilibrium model of low aspect-
    # ratio toroidal plasma with energetic beam ions." Physics of Plasmas (1994-present) 10.8 (2003):
    vpara = dot(v,b)
    cB = curlB(M,gcp.r,0.0,gcp.z)/Babs
    term1 = vpara*dot(b,cB)/Ω_c
    gB = gradB(M,gcp.r,0.0,gcp.z)
    term2 = -dot(r_gyro,gB)/(2*Babs)

    correction = (1 - term1 - term2)
    if (correction <= 0.0) || (correction >= 2.0)
        @warn "Gyro correction results in negative distances or too large of a shift: $(correction)"
    end

    r_gyro = r_gyro*correction

    return r_gyro
end

"""
    larmor_radius(M::AbstractEquilibrium, gcp::GCParticle)

Returns the average larmor radius
"""
function larmor_radius(M::AbstractEquilibrium, gcp::GCParticle; N=10)
    return mean(norm(gyro_step(M, gcp, g)) for g in range(0,2pi,length=N))
end


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const config = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                               optimal_status=MOI.LOCALLY_SOLVED)
const config_no_duals = MOIT.TestConfig(atol=1e-4, rtol=1e-4, duals=false,
                                        optimal_status=MOI.LOCALLY_SOLVED)

@testset "MOI utils" begin
    optimizer = MadNLP.Optimizer()
    @testset "SolverName" begin
        @test MOI.get(optimizer, MOI.SolverName()) == "MadNLP"
    end
    @testset "supports_default_copy_to" begin
        @test MOIU.supports_default_copy_to(optimizer, false)
        @test !MOIU.supports_default_copy_to(optimizer, true)
    end
    @testset "MOI.Silent" begin
        @test MOI.supports(optimizer, MOI.Silent())
        @test MOI.get(optimizer, MOI.Silent()) == false
        MOI.set(optimizer, MOI.Silent(), true)
        @test MOI.get(optimizer, MOI.Silent()) == true
    end
    @testset "MOI.TimeLimitSec" begin
        @test MOI.supports(optimizer, MOI.TimeLimitSec())
        my_time_limit = 10.
        MOI.set(optimizer, MOI.TimeLimitSec(), my_time_limit)
        @test MOI.get(optimizer, MOI.TimeLimitSec()) == my_time_limit
    end
    @testset "MOI.MaxIter" begin
        MOI.set(optimizer,MOI.RawParameter("max_iter"),1)
        @test MOI.get(optimizer,MOI.RawParameter("max_iter")) == 1
    end
end

@testset "Testing getters" begin
    MOIT.copytest(MOI.instantiate(()->MadNLP.Optimizer(print_level=MadNLP.ERROR),
                                  with_bridge_type=Float64), MOIU.Model{Float64}())
end

@testset "Bounds set twice" begin
    optimizer = MadNLP.Optimizer(print_level=MadNLP.ERROR)
    MOIT.set_lower_bound_twice(optimizer, Float64)
    MOIT.set_upper_bound_twice(optimizer, Float64)
end

@testset "MOI Linear tests" begin
    optimizer = MadNLP.Optimizer(print_level=MadNLP.ERROR)
    exclude = ["linear1", # modify constraints not allowed
               "linear5", # modify constraints not allowed
               "linear6", # constraint set for l/q not allowed
               "linear7",  # VectorAffineFunction not supported.
               "linear8a", # Behavior in infeasible case doesn't match test.
               "linear8b", # Behavior in unbounded case doesn't match test.
               "linear8c", # Behavior in unbounded case doesn't match test.
               "linear10", # Interval not supported yet
               "linear10b", # Interval not supported yet
               "linear11", # Variable cannot be deleted
               "linear12", # Behavior in infeasible case doesn't match test.
               "linear14", # Variable cannot be deleted
               "linear15", # VectorAffineFunction not supported.
               ]
    MOIT.contlineartest(optimizer, config_no_duals,exclude)
end

@testset "MOI NLP tests" begin
    optimizer = MadNLP.Optimizer(print_level=MadNLP.ERROR)
    exclude = [
        "feasibility_sense_with_objective_and_no_hessian", # we need Hessians
        "feasibility_sense_with_no_objective_and_no_hessian", # we need Hessians
        "hs071_no_hessian", # we need Hessians
        "hs071_hessian_vector_product_test", # Hessian-vector product is needed
    ]
    MOIT.nlptest(optimizer,config,exclude)
end

@testset "Unit" begin
    bridged = MOIB.full_bridge_optimizer(MadNLP.Optimizer(print_level=MadNLP.ERROR),Float64)
    exclude = ["delete_variable", # Deleting not supported.
               "delete_variables", # Deleting not supported.
               "getvariable", # Variable names not supported.
               "solve_zero_one_with_bounds_1", # Variable names not supported.
               "solve_zero_one_with_bounds_2", # Variable names not supported.
               "solve_zero_one_with_bounds_3", # Variable names not supported.
               "getconstraint", # Constraint names not suported.
               "variablenames", # Variable names not supported.
               "solve_with_upperbound", # loadfromstring!
               "solve_with_lowerbound", # loadfromstring!
               "solve_integer_edge_cases", # loadfromstring!
               "solve_affine_lessthan", # loadfromstring!
               "solve_affine_greaterthan", # loadfromstring!
               "solve_affine_equalto", # loadfromstring!
               "solve_affine_interval", # loadfromstring!
               "get_objective_function", # Function getters not supported.
               "solve_constant_obj",  # loadfromstring!
               "solve_blank_obj", # loadfromstring!
               "solve_singlevariable_obj", # loadfromstring!
               "solve_objbound_edge_cases", # ObjectiveBound not supported.
               "solve_affine_deletion_edge_cases", # Deleting not supported.
               "solve_unbounded_model", # `NORM_LIMIT`
               "number_threads", # NumberOfThreads not supported
               "delete_nonnegative_variables", # get ConstraintFunction n/a.
               "update_dimension_nonnegative_variables", # get ConstraintFunction n/a.
               "delete_soc_variables", # VectorOfVar. in SOC not supported
               "solve_result_index", # DualObjectiveValue not supported
               "time_limit_sec", #time limit given as Flaot64?
               "solve_farkas_interval_lower",
               "solve_farkas_lessthan",
               "solve_farkas_interval_upper",
               "solve_farkas_greaterthan",
               "solve_farkas_variable_lessthan_max",
               "solve_farkas_variable_lessthan",
               "solve_farkas_equalto_lower",
               "solve_farkas_equalto_upper",
               "solve_qp_edge_cases"
               ]
    MOIT.unittest(bridged, config, exclude)
end

@testset "MOI QP/QCQP tests" begin
    optimizer = MadNLP.Optimizer(print_level=MadNLP.ERROR)
    qp_optimizer = MOIU.CachingOptimizer(MOIU.Model{Float64}(), optimizer)
    MOIT.qptest(qp_optimizer, config)
    exclude = ["qcp1", # VectorAffineFunction not supported.
              ]
    MOIT.qcptest(qp_optimizer, config_no_duals, exclude)
end


text:
# %% Standard map island shrinking
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot
grid = 0.0:0.5:2π
sm = Systems.standardmap()
ks = [0.1, 1.0, 4.0]

fig, axs = subplots(1,3; sharey = true)

for i in 1:3
    set_parameter!(sm, 1, ks[i])
    ax = axs[i]
    for θ in grid, p in grid
        u0 = SVector(θ, p)
        tr = trajectory(sm, 1000, u0)
        c = DynamicalSystems.lyapunovspectrum(sm, 4000, 1; u0 = u0)[1] > 0.01 ? "C0" : "C2"
        ax.plot(columns(tr)...; ls = "None", marker = "o", ms = 0.2, c = c)
    end
    ax.set_xlim(0, 2π)
    ax.set_ylim(0, 2π)
    ax.text(4.5, 5.5, "\$k=$(ks[i])\$"; bbox = bbox)
    ax.set_xticks(0:2:6)
end

axs[1].set_ylabel("\$p\$", labelpad = -5)
axs[2].set_xlabel("\$\\theta\$", labelpad = -15)
fig.tight_layout(pad=0.3)
# fig.subplots_adjust(left = 0.05, bottom = 0.14, right = 0.95, top = 0.97, wspace = 0.06)
wsave(plotsdir("8","kam_sm"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI   = point(6,5,3)
const pII  = point(5,2,4)
const pIII = point(7,1,0)
const pIV  = point(3,4,0)

const P = simplex(pI,pII,pIII,pIV)
const Q =  simplex(pI,pII,pIII,pIV)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 18
ct_ref = CommonVolume6D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
   ((x-pI)'*(y-pI))*return(exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(u,v)
   if u[3]< 0.0 || u[3]>1-u[1]-u[2] || u[2]< 0.0 || u[2]>1-u[1] || u[1]<0.0 || u[1]>1.0
      println("u: ",u)
   end
   if v[3]< 0.0 || v[3]>1-v[1]-v[2] || v[2]< 0.0 || v[2]>1-v[1] || v[1]<0.0 || v[1]>1.0
      println("v: ",v)
   end
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, ct_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm =[]
n1 = []
n2 = []
n3 = []
for i in 2:1:15
   Accuracy = i
   ct = CommonVolume6D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, ct)

   num_pts = 18*length(ct.qps)^6
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

for i in 2:1:6
   Accuracy = i
   ct_s = CommonVolume6D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._shunnham2D(Accuracy),
                          SauterSchwab3D._shunnham4D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, ct_s)

   num_pts = 2*length(ct_s.qps[1])^2*length(ct_s.qps[3])+
            16*length(ct_s.qps[3])*length(ct_s.qps[2])
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   ct_gm = CommonVolume6D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1),
                          SauterSchwab3D._grundmannMoeller4D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, ct_gm)

   num_pts = 2*length(ct_gm.qps[1])^2*length(ct_gm.qps[3])+
            16*length(ct_gm.qps[3])*length(ct_gm.qps[2])
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times" )
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex-Product SH",markershape=:rect)
plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e2,1e7),ylims=(1e-8,1))
plot!(xlabel="#Quad. Pts./Func. Evals", ylabel="Rel. Error.", title="Identical Tetrahedron 6D",legend=:bottomleft)

using BenchmarkTools

ref = 1.749281958454714 - 4.409866361065269im

ct = CommonVolume6D(sing,SauterSchwab3D._legendre(5,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, ct)

num_pts = 18*length(ct.qps)^6
println("#Pts: ",num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println("Rell Err; ",err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, ct) end

ct_s = CommonVolume6D_S(sing,(SauterSchwab3D._legendre(5,0.0,1.0),
                         SauterSchwab3D._shunnham2D(5),
                         SauterSchwab3D._shunnham4D(5)))

int_sp = sauterschwab_parameterized(INTEGRAND, ct_s)

num_pts = 2*length(ct_s.qps[1])^2*length(ct_s.qps[3])+
          16*length(ct_s.qps[3])*length(ct_s.qps[2])
println("#Pts: ",num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println("Rel Err: ",err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, ct_s) end

#=
ct_gm = CommonVolume6D_S(sing,(SauterSchwab3D._legendre(7,0.0,1.0),
                         SauterSchwab3D._grundmannMoeller2D(2*7-1),
                         SauterSchwab3D._grundmannMoeller4D(2*7-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, ct_gm)

num_pts = 2*length(ct_gm.qps[1])^2*length(ct_gm.qps[3])+
          16*length(ct_gm.qps[3])*length(ct_gm.qps[2])
println("#Pts: ",num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println("Rel Err: ",err_gm)

=#


text:
export LocalOperator, X, Y, Z
using MacroTools

"""
    LocalOperator{T, N, O, MT}

Local operator on `N` qubits/sites with `O` eigen states.
"""
struct LocalOperator{T <: Number, N, O, MT <: AbstractMatrix{T}} <: AbstractMatrix{T}
    val::MT
end

LocalOperator(N::Int, O::Int, val::MT) where {T, MT <: AbstractMatrix{T}} = LocalOperator{T, N, O, MT}(val)

const X = LocalOperator(1, 2, Pauli.X{ComplexF64}())
const Y = LocalOperator(1, 2, Pauli.Y{ComplexF64}())
const Z = LocalOperator(1, 2, Pauli.Z{ComplexF64}())

Base.summary(io::IO, op::LocalOperator{T, N, O}) where {T, N, O} = print(io, "LocalOperator{$T, $N, $O}")

# inherit other array interface
MacroTools.@forward LocalOperator.val Base.getindex,
                                      Base.size, Base.length
                                      Base.stride, Base.strides

Base.eltype(::LocalOperator{T}) where T = T
Base.:(==)(lhs::LocalOperator{T, N, O}, rhs::LocalOperator{T, N, O}) where {T, N, O} = lhs.val == rhs.val

Base.getindex(A::LocalOperator{<:Number, N}, lhs::NTuple{N, <:Integer}, rhs::NTuple{N, <:Integer}) where N =
    getindex(A.val, to_indices(A, (lhs, rhs))...)

@generated function Base.to_index(A::LocalOperator{<:Number, N, O}, i::NTuple{N, <:Integer}) where {N, O}
    ex = :(i[1] - 1)
    for k = 2:N
        factor = O^(k-1)
        ex = :($ex + $factor * (i[$k] - 1))
    end
    :($ex + 1)
end

Base.:(*)(lhs::LocalOperator{T, N, O, MT1}, rhs::LocalOperator{T, N, O, MT2}) where {T, N, O, MT1, MT2} =
    LocalOperator(N, O, lhs.val * rhs.val)

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
endpoint(::GitHubAPI, ::typeof(is_member), org::AStr, user::AStr) =
    Endpoint(:GET, "/orgs/$org/members/$user"; allow_404=true)
@not_implemented(::GitHubAPI, ::typeof(is_member), ::String, ::Int64)
postprocessor(::GitHubAPI, ::typeof(is_member)) = DoSomething(ismemberorcollaborator)
into(::GitHubAPI, ::typeof(is_member)) = Bool

@not_implemented(::GitHubAPI, ::typeof(groups), )


text:
using AeroMDAO

##
fuse = Fuselage(0.2:0.1:1., 1.:10)

##
projected_area(fuse)

##
length(fuse)

##
cosine_spacing(fuse, 40)

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
const crouch_height = -0.08
const normal_height = -0.16

function step_script(s::mjSim, robot)
    elapsed_time = round(Int, s.d.time * 1000)  # elapsed time in milliseconds (non-paused simulation)

    # check every 100 milliseconds for another action to take
    if !s.paused && elapsed_time % 100 == 0 && elapsed_time > 0
        #println(elapsed_time, ": ", elapsed_time, "\tframecount: ", round(Int, s.framecount))

        # After he's done falling and getting up, we return to a normal height and pitch
        if elapsed_time == 2000 && robot.command.height > -0.1
            robot.command.height = normal_height
            robot.command.pitch = 0.0
            # We begin trotting here for a few seconds
            toggle_trot(robot)
            println("Standing up and beginning march with velocity", robot.command.horizontal_velocity)
        end
    end
end


text:
RD.@autodiff struct Cartpole{T} <: RobotDynamics.ContinuousDynamics
    mc::T
    mp::T
    l::T
    g::T
end

Cartpole() = Cartpole{Float64}(1.0, 0.2, 0.5, 9.81)

function RobotDynamics.dynamics(model::Cartpole, x, u)
    mc = model.mc  # mass of the cart in kg (10)
    mp = model.mp   # mass of the pole (point mass at the end) in kg
    l = model.l   # length of the pole in m
    g = model.g  # gravity m/s^2

    q = x[ @SVector [1,2] ]
    qd = x[ @SVector [3,4] ]

    s = sin(q[2])
    c = cos(q[2])

    H = @SMatrix [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = @SMatrix [0 -mp*qd[2]*l*s; 0 0]
    G = @SVector [0, mp*g*l*s]
    B = @SVector [1, 0]

    qdd = -H\(C*qd + G - B*u[1])
    return [qd; qdd]
end

function RobotDynamics.dynamics!(model::Cartpole, xdot, x, u)
    mc = model.mc  # mass of the cart in kg (10)
    mp = model.mp   # mass of the pole (point mass at the end) in kg
    l = model.l   # length of the pole in m
    g = model.g  # gravity m/s^2

    q = x[ @SVector [1,2] ]
    qd = x[ @SVector [3,4] ]

    s = sin(q[2])
    c = cos(q[2])

    H = @SMatrix [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = @SMatrix [0 -mp*qd[2]*l*s; 0 0]
    G = @SVector [0, mp*g*l*s]
    B = @SVector [1, 0]

    qdd = -H\(C*qd + G - B*u[1])
    xdot[1] = qd[1]
    xdot[2] = qd[2]
    xdot[3] = qdd[1]
    xdot[4] = qdd[2]
    return nothing
end

function RD.jacobian!(model::Cartpole, J, xdot, x, u, t)
    mc = model.mc  # mass of the cart in kg (10)
    mp = model.mp   # mass of the pole (point mass at the end) in kg
    l = model.l   # length of the pole in m
    g = model.g  # gravity m/s^2

    q = x[ @SVector [1,2] ]
    qd = x[ @SVector [3,4] ]

    s = sin(q[2])
    c = cos(q[2])

    H = @SMatrix [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = @SMatrix [0 -mp*qd[2]*l*s; 0 0]
    G = @SVector [0, mp*g*l*s]
    B = @SVector [1, 0]
    qdd = -H\(C*qd + G - B*u[1])

    ∂H∂z = SA[
        0 -mp*l*s*qdd[2] 0 0 0
        0 -mp*l*s*qdd[1] 0 0 0
    ]
    ∂C∂z = SA[
        0 -mp*l*c*qd[2]^2 0 -2*mp*l*qd[2]*s 0
        0 0 0 0 0
    ]
    ∂G∂z = SA[
       0 0 0 0 0
       0 mp*g*l*c 0 0 0
    ]
    ∂B∂z = SA[
        0 0 0 0 1
        0 0 0 0 0
    ]
    J[1:2,:] .= 0
    J[1,3] = 1.0
    J[2,4] = 1.0
    J[3:4,:] .= H\(-∂H∂z - ∂C∂z - ∂G∂z + ∂B∂z)
    return nothing
end

RobotDynamics.state_dim(::Cartpole) = 4
RobotDynamics.control_dim(::Cartpole) = 1


text:
bread(D, V) = D' * V^-1 * D
meat(D, V, resid) = D' * V^-1 * resid * resid' * V^-1 * D
sandwich(bread::Vector{Matrix{T}}, meat::Vector{Matrix{T}}) where {T<:Real} = sum(bread)^-1 * sum(meat) * sum(bread)^-1

function vcov(m::GeneralizedEstimatingEquationModel)
    V = [updateV(m.R, μ) for μ in m.GEE.μᵢ]
    D = [updateD(m.d, X, μ) for (X, μ) in zip(m.GEE.Xᵢ, m.GEE.μᵢ)]
    rᵢ = [y - μ for (y, μ) in zip(m.GEE.yᵢ, m.GEE.μᵢ)]
    return sandwich(bread.(D, V), meat.(D, V, rᵢ))
end

function stderror(m::GeneralizedEstimatingEquationModel)
    [sqrt(d) for d in diag(vcov(m))]
end


text:
"""
Submodule of `SirenSeq`.
Standard musical scales are defined here.
"""
module Scales

## a scale is a function that takes (ocv::Int, val::Int) as input and outputs a midi pitch value {0..127}
## `ocv` is an octave offset while `val` is the note interval on the scale
## scales can be any length between 1 and 11

export

## standard 7 note scales

aLyd,aMaj,aMix,aDor,aMin,aPhr,aLoc,
bLyd,bMaj,bMix,bDor,bMin,bPhr,bLoc,
cLyd,cMaj,cMix,cDor,cMin,cPhr,cLoc,
dLyd,dMaj,dMix,dDor,dMin,dPhr,dLoc,
eLyd,eMaj,eMix,eDor,eMin,ePhr,eLoc,
fLyd,fMaj,fMix,fDor,fMin,fPhr,fLoc,
gLyd,gMaj,gMix,gDor,gMin,gPhr,gLoc,

## a scale which ignores `ocv` and uses `val` as the midi pitch value
## useful for percussion tracks

noScale,

## returns the pitch value for a given scale and given (octave,interval) pair

getPitch

abstract Scale

type ChromaticMode <: Scale
	pmap::Dict{Int,Int}
	basev::Int
	name::AbstractString
end

type PercussionScale <: Scale
	name::AbstractString
end

function Base.isequal(x::ChromaticMode, y::ChromaticMode)
	if !isequal(x.pmap,y.pmap) ; return false ; end
	if !isequal(x.basev,y.basev) ; return false ; end
	if !isequal(x.name,y.name) ; return false ; end
	true
end

function Base.isequal(x::PercussionScale, y::PercussionScale)
	true
end

function Base.isequal(x::PercussionScale, y::ChromaticMode)
	false
end

function toPitch(ocv::Int, val::Int, scaleDict::Dict{Int,Int}, ofs::Int)
	ssz = length(scaleDict) ;  @assert  0 < ssz < 12
	while val < 1
		val += ssz
		ocv -= 1
	end
	while val > ssz
		val -= ssz
		ocv += 1
	end
	clamp(ofs+scaleDict[val]+60+12*(ocv-4),0,127)
end

## returns the pitch value for a given scale and given (octave,interval) pair
function getPitch(sca::ChromaticMode, ocv::Int, val::Int)
	toPitch(ocv,val,sca.pmap,sca.basev)
end

function getPitch(sca::PercussionScale, ocv::Int, val::Int)
	clamp(val,0,127)
end

const cLydDc = Dict{Int,Int}( 1=>0, 2=>2, 3=>4, 4=>6, 5=>7, 6=>9, 7=>11 )
const cMajDc = Dict{Int,Int}( 1=>0, 2=>2, 3=>4, 4=>5, 5=>7, 6=>9, 7=>11 )
const cMixDc = Dict{Int,Int}( 1=>0, 2=>2, 3=>4, 4=>5, 5=>7, 6=>9, 7=>10 )
const cDorDc = Dict{Int,Int}( 1=>0, 2=>2, 3=>4, 4=>5, 5=>7, 6=>8, 7=>10 )
const cMinDc = Dict{Int,Int}( 1=>0, 2=>2, 3=>3, 4=>5, 5=>7, 6=>8, 7=>10 )
const cPhrDc = Dict{Int,Int}( 1=>0, 2=>1, 3=>3, 4=>5, 5=>7, 6=>8, 7=>10 )
const cLocDc = Dict{Int,Int}( 1=>0, 2=>1, 3=>3, 4=>5, 5=>6, 6=>8, 7=>10 )

noScale = PercussionScale("perc")

aLyd = ChromaticMode(cLydDc,-3,"cLyd")
aMaj = ChromaticMode(cMajDc,-3,"cMaj")
aMix = ChromaticMode(cMixDc,-3,"cMix")
aDor = ChromaticMode(cDorDc,-3,"aDor")
aMin = ChromaticMode(cMinDc,-3,"aMin")
aPhr = ChromaticMode(cPhrDc,-3,"aPhr")
aLoc = ChromaticMode(cLocDc,-3,"aLoc")

bLyd = ChromaticMode(cLydDc,-1,"bLyd")
bMaj = ChromaticMode(cMajDc,-1,"bMaj")
bMix = ChromaticMode(cMixDc,-1,"bMix")
bDor = ChromaticMode(cDorDc,-1,"bDor")
bMin = ChromaticMode(cMinDc,-1,"bMin")
bPhr = ChromaticMode(cPhrDc,-1,"bPhr")
bLoc = ChromaticMode(cLocDc,-1,"bLoc")

cLyd = ChromaticMode(cLydDc,0,"cLyd")
cMaj = ChromaticMode(cMajDc,0,"cMaj")
cMix = ChromaticMode(cMixDc,0,"cMix")
cDor = ChromaticMode(cDorDc,0,"cDor")
cMin = ChromaticMode(cMinDc,0,"cMin")
cPhr = ChromaticMode(cPhrDc,0,"cPhr")
cLoc = ChromaticMode(cLocDc,0,"cLoc")

dLyd = ChromaticMode(cLydDc,2,"dLyd")
dMaj = ChromaticMode(cMajDc,2,"dMaj")
dMix = ChromaticMode(cMixDc,2,"dMix")
dDor = ChromaticMode(cDorDc,2,"dDor")
dMin = ChromaticMode(cMinDc,2,"dMin")
dPhr = ChromaticMode(cPhrDc,2,"dPhr")
dLoc = ChromaticMode(cLocDc,2,"dLoc")

eLyd = ChromaticMode(cLydDc,4,"eLyd")
eMaj = ChromaticMode(cMajDc,4,"eMaj")
eMix = ChromaticMode(cMixDc,4,"eMix")
eDor = ChromaticMode(cDorDc,4,"eDor")
eMin = ChromaticMode(cMinDc,4,"eMin")
ePhr = ChromaticMode(cPhrDc,4,"ePhr")
eLoc = ChromaticMode(cLocDc,4,"eLoc")

fLyd = ChromaticMode(cLydDc,6,"fLyd")
fMaj = ChromaticMode(cMajDc,6,"fMaj")
fMix = ChromaticMode(cMixDc,6,"fMix")
fDor = ChromaticMode(cDorDc,6,"fDor")
fMin = ChromaticMode(cMinDc,6,"fMin")
fPhr = ChromaticMode(cPhrDc,6,"fPhr")
fLoc = ChromaticMode(cLocDc,6,"fLoc")

gLyd = ChromaticMode(cLydDc,7,"gLyd")
gMaj = ChromaticMode(cMajDc,7,"gMaj")
gMix = ChromaticMode(cMixDc,7,"gMix")
gDor = ChromaticMode(cDorDc,7,"gDor")
gMin = ChromaticMode(cMinDc,7,"gMin")
gPhr = ChromaticMode(cPhrDc,7,"gPhr")
gLoc = ChromaticMode(cLocDc,7,"gLoc")

end


text:
module TensorAlgebra

using LinearAlgebra

export Tensor, Covector, VectorSpace, DualSpace, ProductSpace, TensorSpace, TensorProduct
export spaces, tensors, label, dimensions, field, degree, dual, domain, ×, ⊗

abstract type AbstractSpace{K,N} end

struct VectorSpace{K,L,D} <: AbstractSpace{K,1} end

struct DualSpace{K,L,D} <: AbstractSpace{K,1} end

abstract type AbstractProductSpace{K,N,S} <: AbstractSpace{K,N} end

struct ProductSpace{K,N,S} <: AbstractProductSpace{K,N,S} end

struct TensorSpace{K,N,S} <: AbstractProductSpace{K,N,S} end

abstract type AbstractTensor{K,N} <: AbstractArray{K,N} end

struct Tensor{K,N,D} <: AbstractTensor{K,N}
    array::Array{K,N}
end

const Covector{K,L,D} = Tensor{K,1,VectorSpace{K,L,D}}

struct TensorProduct{K,N,S} <: AbstractTensor{K,N}
    scalar::K
    tensors::S
end

Base.collect(t::Tensor) = t.array

field(::AbstractSpace{K}) where {K} = K

field(::AbstractTensor{K}) where {K} = K

field(::Vector{K}) where {K} = K

degree(::AbstractSpace{K,N}) where {K,N} = N

degree(::AbstractTensor{K,N}) where {K,N} = N

spaces(as::AbstractSpace{K,1}) where {K} = (as,)

spaces(::AbstractProductSpace{K,N,S}) where {K,N,S} = S

label(::Union{VectorSpace{K,L},DualSpace{K,L}}) where {K,L} = L

dimensions(::Union{VectorSpace{K,L,D},DualSpace{K,L,D}}) where {K,L,D} = D

dual(::VectorSpace{K,L,D}) where {K,L,D} = DualSpace{K,L,D}()

dual(::DualSpace{K,L,D}) where {K,L,D} = VectorSpace{K,L,D}()

function Covector(vs::VectorSpace{K,L,D}, a) where {K,L,D}
    size(a)[1] === D || throw(DimensionMismatch("Tried to create a covector of dimension $(size(a)[1]) in the vector space $(label(vs)) of dimension $(D)."))
    Tensor{K,1,TensorSpace{K,1,(VectorSpace{K,L,D}(),)}}(a)
end

function Vector(vs::VectorSpace{K,L,D}, a) where {K,L,D}
    size(a)[1] === D || throw(DimensionMismatch("Tried to create a vector of dimension $(size(a)[1]) in the vector space $(label(vs)) of dimension $(D)."))
    Tensor{K,1,TensorSpace{K,1,(DualSpace{K,L,D}(),)}}(a)
end

VectorSpace(L::Symbol,::Type{K},dim::Int) where {K,D} = VectorSpace{K,L,dim}()

ProductSpace(args::AbstractSpace{K,1}...) where {K} = ProductSpace{K,length(args),(args...,)}()

TensorSpace(args::Vararg{<:AbstractSpace{K,1}}) where {K} = TensorSpace{K,length(args),(args...,)}()

TensorSpace(::ProductSpace{K,N,S}) where {K,N,S} = TensorSpace{K,N,S}()

TensorSpace(ts::TensorSpace{K,N,S}) where {K,N,S} = ts

(::Type{PS})(args...) where {PS <: ProductSpace} = ProductSpace(args...)

(::Type{TS})(args...) where {TS <: TensorSpace} = TensorSpace(args...)

dual(ps::PS) where {PS <: AbstractProductSpace} = PS(dual.(spaces(ps))...)

Tensor(dom::D,a) where {K,N,D <: AbstractSpace{K,N}} = Tensor{K,N,typeof(TensorSpace(dom))}(a)

Tensor(::D,a) where {K,N,D <: TensorSpace{K,N}} = Tensor{K,N,D}(a)

TensorProduct(args::Tensor{K}...) where {K} = TensorProduct{K,sum(degree.(args)),typeof((args...,))}(one(K), (args...,))

TensorProduct(t::Tensor) = t

domain(::Tensor{K,N,D}) where {K,N,D} = D()

domain(tp::TensorProduct{K,N,S}) where {K,N,S} = TensorSpace{K,N,((spaces.(domain.(tensors(tp)))...)...,)}()

dimensions(at::AbstractTensor) = ((dimensions.(spaces(domain(at)))...)...,)

scalar(::Tensor{K}) where {K} = one(K)

scalar(tp::TensorProduct) = tp.scalar

tensors(t::Tensor) = t

tensors(tp::TensorProduct) = tp.tensors

function (f::AbstractTensor{K,2})(::typeof(-), x::Tensor{K,1}) where {K}
    dom = spaces(domain(f))
    codom = spaces(domain(x))
    dual(dom[2]) === codom[1] || throw(DomainError(x,"Domain mismatch: Expected $(dual(dom[2])). Got $(codom[1])."))
    Tensor(dom[1], collect(f) * collect(x))
end

function (f::AbstractTensor{K,2})(x::Tensor{K,1}, ::typeof(-)) where {K}
    dom = spaces(domain(f))
    codom = spaces(domain(x))
    dual(dom[1]) === codom[1] || throw(DomainError(x,"Domain mismatch: Expected $(dual(dom[1])). Got $(codom[1])."))
    Tensor(dom[2], (collect(x)' * collect(f)).parent)
end

function (f::AbstractTensor{K,N})(x::AbstractTensor{K,N}) where {K,N}
    dual(domain(f)) === domain(x) || throw(DomainError(x,"Domain mismatch: Expected $(dual(domain(f))). Got $(domain(x))."))
    dot(f, x)
end

(f::Tensor{K,N})(::Vararg{typeof(-),N}) where {K,N} = f

function (f::Tensor{K,1})(x::Tensor{K,1}) where {K}
    dual(domain(f)) === domain(x) || throw(DomainError(x,"Domain mismatch: Expected $(dual(domain(f))). Got $(domain(x))."))
    dot(f, x)
end

function (f::Tensor{K,N})(xs::Vararg{Tensor{K,1},N}) where {K,N}
    x = TensorProduct(xs...)
    dual(domain(f)) === domain(x) || throw(DomainError(x,"Domain mismatch: Expected $(dual(domain(f))). Got $(domain(x))."))
    dot(f, x)
end

function (f::AbstractTensor{K,N})(xs::Vararg{Union{typeof(-),AbstractTensor{K,1}},N}) where {K,N}
    isarg = xs .!== -
    indims = findall(isarg)
    x = TensorProduct(xs[indims]...)
    any(.!isarg) || return dot(f, x)
    outdims = findall(.!isarg)
    a = reshape(mapslices(f, dims=indims) do slice
        dot(x, slice)
    end,size(f)[outdims]...)
    d = TensorSpace(spaces(domain(f))[outdims]...)
    Tensor(d, a)
end

Base.size(t::Tensor) = size(t.array)

Base.size(t::TensorProduct) = ((size.(t.tensors)...)...,)

Base.getindex(t::Tensor,ix::Vararg{Int}) = getindex(t.array, ix...)

function Base.getindex(tp::TensorProduct, ix::Vararg{Int})
    value = one(field(tp))
    offset = 0
    for (i, tensor) in enumerate(tp.tensors)
        it = degree(tensor) === 1 ? ix[offset + 1] : ix[offset + 1:offset + degree(tensor)]
        value *= getindex(tensor, it...)
        offset += degree(tensor)
    end
    value * scalar(tp)
end

Base.:^(v::AbstractSpace,::typeof(*)) = dual(v)

Base.first(ts::AbstractProductSpace) = first(spaces(ts))

Base.last(ts::AbstractProductSpace) = last(spaces(ts))

Base.in(t::AbstractTensor{K}, vs::AbstractSpace{K}) where {K} = dual(domain(t)) === TensorSpace(vs)

function Base.:+(at1::AbstractTensor{K,N}, at2::AbstractTensor{K,N}) where {K,N}
    domain(at1) === domain(at2) || throw(DomainError(at2,"Domain mismatch: Expected $(domain(at1)). Got $(domain(at2))."))
    Tensor{K,N,typeof(domain(at1))}(collect(at1) + collect(at2))
end

function Base.:-(at1::AbstractTensor{K,N}, at2::AbstractTensor{K,N}) where {K,N}
    domain(at1) === domain(at2) || throw(DomainError(at2,"Domain mismatch: Expected $(domain(at1)). Got $(domain(at2))."))
    Tensor{K,N,typeof(domain(at1))}(collect(at1) - collect(at2))
end

function Base.:-(t::Tensor{K,N,D}) where {K,N,D}
    Tensor{K,N,D}(-t.array)
end

function Base.:-(tp::TensorProduct{K,N,S}) where {K,N,S}
    TensorProduct{K,N,S}(-tp.scalar, tp.tensors)
end

Base.:*(x::Number,t::Tensor{K,N,D}) where {K,N,D} = Tensor{K,N,D}(x * t.array)

Base.:*(t::Tensor{K,N,D},x::Number) where {K,N,D} = Tensor{K,N,D}(x * t.array)

Base.:*(t1::Tensor{K,1,D},t2::Tensor{K,1,D}) where {K,D} = Tensor{K,1,D}(t1.array .* t2.array)

Base.:*(x::Number,tp::TensorProduct{K,N,S}) where {K,N,S} = TensorProduct{K,N,S}(x * scalar(tp), tensors(tp))

Base.:*(tp::TensorProduct{K,N,S},x::Number) where {K,N,S} = x * tp

×(vs1::AbstractSpace{K,1},vs2::AbstractSpace{K,1}) where {K} = ProductSpace(vs1, vs2)

×(vs::AbstractSpace{K,1},ts::ProductSpace{K}) where {K} = ProductSpace(vs, spaces(ts)...)

×(ts::ProductSpace{K},vs::AbstractSpace{K,1}) where {K} = ProductSpace(spaces(ts)..., vs)

×(ts1::ProductSpace{K},ts2::ProductSpace{K}) where {K} = ProductSpace((spaces(ts1)..., spaces(ts2)...))

⊗(vs1::AbstractSpace{K,1},vs2::AbstractSpace{K,1}) where {K} = TensorSpace(vs1, vs2)

⊗(vs::AbstractSpace{K,1},ts::TensorSpace{K}) where {K} = TensorSpace(vs, spaces(ts)...)

⊗(ts::TensorSpace{K},vs::AbstractSpace{K,1}) where {K} = TensorSpace(spaces(ts)..., vs)

⊗(ts1::TensorSpace{K},ts2::TensorSpace{K}) where {K} = TensorSpace((spaces(ts1)..., spaces(ts2)...))

⊗(t1::Tensor{K,R1},t2::Tensor{K,R2}) where {K,R1,R2} = TensorProduct{K,R1 + R2,typeof((t1, t2))}(one(K), (t1, t2))

⊗(tp::TensorProduct{K,R1},t::Tensor{K,R2}) where {K,R1,R2} = TensorProduct{K,R1 + R2,typeof((tensors(tp)..., t))}(scalar(tp), (tensors(tp)..., t))

⊗(t::Tensor{K,R1},tp::TensorProduct{K,R2}) where {K,R1,R2} = TensorProduct{K,R1 + R2,typeof((t, tensors(tp)...))}(scalar(tp), (t, tensors(tp)...))

⊗(t1::TensorProduct{K,R1},t2::TensorProduct{K,R2}) where {K,R1,R2} = TensorProduct{K,R1 + R2,typeof((tensors(t1)..., tensors(t2)...))}(scalar(t1) * scalar(t2), (tensors(t1)..., tensors(t2)...))

⊗(t1::AbstractTensor,t2::AbstractTensor,ts::Vararg{<:AbstractTensor}) = TensorProduct(scalar(t1) * scalar(t2), t1 ⊗ t2, ts...)

Base.show(io::IO, ::VectorSpace{K,L}) where {K,L} = print(io, L)

# Base.show(io::IO, ::Type{VectorSpace{K,L}}) where {K,L} = print(io, L)

Base.show(io::IO, ::DualSpace{K,L}) where {K,L} = print(io, L, "^*")

# Base.show(io::IO, ::Type{DualSpace{K,L}}) where {K,L} = print(io, L, "^*")

# Base.show(io::IO, ::ProductSpace{K,N,S}) where {K,N,S} = print(io, join(S, " × "))

# Base.show(io::IO, ::Type{ProductSpace{K,N,S}}) where {K,N,S} = print(io, join(S, " × "))

# Base.show(io::IO, ::TensorSpace{K,N,S}) where {K,N,S} = print(io, join(S, " ⊗ "))

# Base.show(io::IO, ::Type{TensorSpace{K,N,S}}) where {K,N,S} = print(io, join(S, " ⊗ "))

end # module


text:
module OpticSim

import Unitful
using LinearAlgebra: eigen, svd, I, qr, dot, cross, norm, det, normalize, inv
using StaticArrays
using FiniteDifferences
using BenchmarkTools
using DataFrames: DataFrame
using Images
using Base: @.
using ForwardDiff
using StringEncodings

unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))

include("GlassCat/GlassCat.jl")
import OpticSim.GlassCat: plot_indices, index, polyfit_indices, absairindex, absorption, info, glassid, glassname, glassforid, isair, findglass, modelglass, glassfromMIL, GlassID

include("Matrix.jl")
include("Constants.jl")
include("Utilities.jl")
include("Geometry/Geometry.jl")
include("Optical/Optical.jl")
include("Visualization.jl")
include("TestData.jl")
include("Examples.jl")
include("Optimization/Optimizable.jl")
include("Diagnostics.jl")

include("Benchmarks.jl")
include("Data/Spectra.jl")

include("Sysimage.jl")

################################################

# This can be used to track NaN, particularly in ForwardDiff gradients, causing problems
# e.g. Diagnostics.testoptimization(lens = Examples.doubleconvex(NaNCheck{Float64}), samples = 1)

# struct NaNCheck{T<:Real} <: Real
#     val::T
#     function NaNCheck{T}(a::S) where {T<:Real, S<:Real}
#         @assert !(T <: NaNCheck)
#         new{T}(T(a))
#     end
# end
# export NaNCheck
# Base.isnan(a::NaNCheck{T}) where{T} = isnan(a.val)
# Base.isinf(a::NaNCheck{T}) where{T} = isinf(a.val)
# Base.typemin(::Type{NaNCheck{T}}) where{T} = NaNCheck{T}(typemin(T))
# Base.typemax(::Type{NaNCheck{T}}) where{T} = NaNCheck{T}(typemax(T))
# Base.eps(::Type{NaNCheck{T}}) where {T} = NaNCheck{T}(eps(T))
# Base.decompose(a::NaNCheck{T}) where {T} = Base.decompose(a.val)
# Base.round(a::NaNCheck{T}, m::RoundingMode) where {T} = NaNCheck{T}(round(a.val, m))

# struct NaNException <: Exception end

# # (::Type{Float64})(a::NaNCheck{S}) where {S<:Real} = NaNCheck{Float64}(Float64(a.val))
# (::Type{T})(a::NaNCheck{S}) where {T<:Integer,S<:Real} = T(a.val)
# (::Type{NaNCheck{T}})(a::NaNCheck{S}) where {T<:Real,S<:Real} = NaNCheck{T}(T(a.val))
# Base.promote_rule(::Type{NaNCheck{T}}, ::Type{T}) where {T<:Number} = NaNCheck{T}
# Base.promote_rule(::Type{T}, ::Type{NaNCheck{T}}) where {T<:Number} = NaNCheck{T}
# Base.promote_rule(::Type{S}, ::Type{NaNCheck{T}}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}
# Base.promote_rule(::Type{NaNCheck{T}}, ::Type{S}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}
# Base.promote_rule(::Type{NaNCheck{S}}, ::Type{NaNCheck{T}}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}

# for op = (:sin, :cos, :tan, :log, :exp, :sqrt, :abs, :-, :atan, :acos, :asin, :log1p, :floor, :ceil, :float)
#     eval(quote
#         function Base.$op(a::NaNCheck{T}) where{T}
#             temp = NaNCheck{T}(Base.$op(a.val))
#             if isnan(temp)
#                 throw(NaNException())
#             end
#             return temp
#         end
#     end)
# end

# for op = (:+, :-, :/, :*, :^, :atan)
#     eval(quote
#         function Base.$op(a::NaNCheck{T}, b::NaNCheck{T}) where{T}
#             temp = NaNCheck{T}(Base.$op(a.val, b.val))
#             if isnan(temp)
#                 throw(NaNException())
#             end
#             return temp
#         end
#     end)
# end

# for op =  (:<, :>, :<=, :>=, :(==), :isequal)
#     eval(quote
#         function Base.$op(a::NaNCheck{T}, b::NaNCheck{T}) where{T}
#             temp = Base.$op(a.val, b.val)
#             return temp
#         end
#     end)
# end

################################################

end # module


text:
include("bpMCSetup.jl")

function makeRecycleLTD(ys::Vector{Float64}, N::Int64)
  rpio::RPIO{Int64Particle} = RPIO{Int64Particle}(N, length(ys), 1)
  function ltd(in::SVector{d, Float64})
    θ::BPθ = toBPθ(in)
    lp::Float64 = lglogprior(θ)
    if lp == -Inf return -Inf end
    model::Model = makeBPModel(θ, ys)
    return lp + rpEstimateDiscrete(model, rpio)
  end
end

logtargetRecycle = makeRecycleLTD(ys, N)
PRecycle = makeAMKernel(logtargetRecycle, d)

Random.seed!(12345)
setRNGs(54321)

recycleElapsed = @elapsed (chainRecycle = simulateChainProgress(PRecycle, t0, chainLength))
rar = PRecycle(:acceptanceRate)
recycleCovEstimate = PRecycle(:covEstimate)

@save "data/bpMCRecycle.jld" chainRecycle rar recycleCovEstimate recycleElapsed


text:
import Obj
import Obj.Objective

"""
    bt([; c, α_0, η, maxiter])

    Returns a LineSearch strategy for Backtracking (Armijo) Line Search
    with given parametrization:

          c: Armijo condition constant (default 1e-1)
        α_0: Initial step size value (default 1.0)
          η: Decrease rate (default: 0.1)
    maxiter: Max number of LS iterations (default: 10)
"""
function bt(; c::Float64=1e-1, α_0::Float64=1.0, η::Float64=0.1, maxiter::Int=10)

    return LineSearch(
        # name
        "bt",

        # strategy
        function (w::Array{Float64,1}, s::Array{Float64,1}, obj::Union{Objective,Void},
                  b::Array{Int,1}, fw::Float64, gw::Array{Float64,1})
            return _bt(w, s, obj, fw, gw, c, α_0, η, maxiter)   # b is ignored
        end
    )

end

@fastmath function _bt(w::Array{Float64,1}, s::Array{Float64,1}, obj::Objective,
                       fw::Float64, gw::Array{Float64,1},
                       c::Float64, α_0::Float64, η::Float64, maxiter::Int)

    @inbounds gws = (gw'*s)[1]
    for i = 1:maxiter

        if Obj.f(obj, w + α_0*s) <= fw + c*α_0*gws # Armijo satisfied?
            return α_0, i-1
        else
            α_0 = η*α_0
        end

    end
    return α_0, maxiter

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
# Find All Occurrences of a Pattern in a String

function find_pattern_occurrences(text, pattern)
    k = length(pattern)
    matching_indices = Int64[]
    for ix in 1:(length(text)-k+1)
        substring = text[ix:ix+k-1]
        if substring == pattern
            push!(matching_indices, ix - 1)
        end
    end
    matching_indices
end

function main()
    pattern = readline()
    genome = readline()
    result = find_pattern_occurrences(genome, pattern)
    println(join(result, " "))
end

main()


text:
#=
Snippets for testing struct types.
=#

struct Foo
    field1::Int
end

isconcretetype(Foo) # true
# Foo is concrete because it is a struct.

fieldtype(Foo, 1)  # Int64
# While declared an Int, the actual type is Int64
# TODO because why?

Int64 <: Int  # true
# Int64 is a subtype of Int

Int64 <: Integer  # true
# Int64 is also a subtype of Integer
# Integer is Abstract supertype for all integers.

Int <: Integer  # true
Integer <: Int  # false
# Int is a subtype of Integer

isconcretetype(Int) # true
isconcretetype(Integer) # false

Int <: Signed   # true
Integer <: Signed   # true

isconcretetype(Int) # true
# Int

x::Int = 1  # ERROR: syntax: type declarations on global variables are not yet supported

x = 1
typeof(x)   # Int64

Int64 :< Int

typeof(Int) # DataType
# the type of a type is DataType

#=
Tests that all fields of given type expr are concrete.
=#
macro isconcretestruct(expr)
    @assert expr.head == :struct
    S = expr.args[2]
    return quote
        $(esc(expr))

        for n in fieldnames($S)
            if !isconcretetype(fieldtype($S, n))
                error("field $n is not concrete")
            end
        end
    end
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
n_values  = [200, 300, 400]
k_values  = [200, 300, 400]
m_values  = [200, 300, 400]

testset_name_suffix = "(main)"

for T ∈ (Float64,Float32,Int64,Int32)
  @time test_complex(T, m_values, k_values, n_values, testset_name_suffix)
  @time test_real(T, m_values, k_values, n_values, testset_name_suffix)
end


text:
# This file was generated, do not modify it. # hide
for m in models(matching(X, y))
    if m.prediction_type == :probabilistic
        println(rpad(m.name, 30), "($(m.package_name))")
    end
end

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
function sqlite3_errcode(dbptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_errcode, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
    @unix_only return ccall(
        (:sqlite3_errcode, sqlite3_lib),
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
end

function sqlite3_extended_errcode(dbptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_extended_errcode, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
    @unix_only return ccall(
        (:sqlite3_extended_errcode, sqlite3_lib),
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
end

function sqlite3_errmsg(dbptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_errmsg, sqlite3_lib),
        stdcall,
        Ptr{Uint8},
        (Ptr{Void}, ),
        dbptr
    )
    @unix_only return ccall(
        (:sqlite3_errmsg, sqlite3_lib),
        Ptr{Uint8},
        (Ptr{Void}, ),
        dbptr
    )
end

function sqlite3_errstr(errcode::Cint)
    @windows_only return ccall(
        (:sqlite3_errstr, sqlite3_lib),
        stdcall,
        Ptr{Uint8},
        (Cint, ),
        errcode
    )
    @unix_only return ccall(
        (:sqlite3_errstr, sqlite3_lib),
        Ptr{Uint8},
        (Cint, ),
        errcode
    )
end


text:
using Base.Test

const LOG_FILE_OPTS = tempname()
println("Path to LOG_FILE_OPTS: $LOG_FILE_OPTS")

configure(; modes = ["debug", "info", "warn", "error", "crazy"])
add_truck(Lumberjack.LumberjackTruck(LOG_FILE_OPTS, nothing, Dict{Any,Any}(:is_colorized => true, :uppercase => true)), "colorlumberjacklogfile")

log("debug", "some-msg")
log("info", "some-msg")
log("warn", "some-msg")
log("error", "some-msg")
log("crazy", "some-msg")

remove_truck("optslumberjacklogfile")

# custom colors
add_truck(Lumberjack.LumberjackTruck(LOG_FILE_OPTS, nothing, Dict{Any,Any}(:colors => Dict{Any,Any}("debug" => :black, "info" => :red, "crazy" => :green), :uppercase => true)), "optslumberjacklogfile")

log("debug", "some-msg")
log("info", "some-msg")
log("crazy", "some-msg")

remove_truck("optslumberjacklogfile")

log_lines = readlines(open(LOG_FILE_OPTS, "r"))

# default colors: {"debug" => :cyan, "info" => :blue, "warn" => :yellow, "error" => :red}
# test with default colors
X = VERSION >= v"0.6.0" ? "" : "\n"
Y = VERSION >= v"0.6.0" ? :default : :normal
@test log_lines[1] == "$(Base.text_colors[:cyan])DEBUG: some-msg$(X)"
@test log_lines[2] == "$(Base.text_colors[Y]Base.text_colors[:blue])INFO: some-msg$(X)"
@test log_lines[3] == "$(Base.text_colors[Y]Base.text_colors[:yellow])WARN: some-msg$(X)"
@test log_lines[4] == "$(Base.text_colors[Y]Base.text_colors[:red])ERROR: some-msg$(X)"
@test log_lines[5] == "$(Base.text_colors[Y])CRAZY: some-msg$(X)"

# test with custom colors
@test log_lines[6] == "$(Base.text_colors[:black])DEBUG: some-msg$(X)"
@test log_lines[7] == "$(Base.text_colors[Y]Base.text_colors[:red])INFO: some-msg$(X)"
@test log_lines[8] == "$(Base.text_colors[Y]Base.text_colors[:green])CRAZY: some-msg$(X)"

# clean up
@test success(`rm $LOG_FILE_OPTS`)


text:
## getFundamentals
##----------------

## calculate fundamental helping values
function getFundamentals(mod::MuSigmaModel)
    nAss = size(mod.mu, 1)

    ## get return moments
    mus = mod.mu
    covMatr = mod.sigma

    ## get helper matrices
    onesN = ones(nAss,1)
    a = ones(1, nAss)/covMatr*ones(nAss, 1)
    b = ones(1, nAss)/covMatr*mus
    c = mus'/covMatr*mus
    d = a.*c-b.^2

    return a[1], b[1], c[1], d[1]
end

## get global minimum variance portfolio with short selling
function getGMVSS(mod::MuSigmaModel)
    nAss = size(mod.mu, 1)

    ## get return moments
    covMatr = mod.sigma

    a, b, c, d = getFundamentals(mod)

    ## compute GMV weights
    xGMV = covMatr\ones(nAss, 1)/a

    return xGMV[:]
end

## efficient frontier
##-------------------

function getEffFrontier(mod::MuSigmaModel;
                        nPoints = 50)

    ## get GMV portfolio
    xGMV = AssetMgmt.getGMVSS(mod)

    ## get associated portfolio mu
    muGMV = AssetMgmt.getPMean(xGMV, mod.mu)

    ## find maximum asset mu
    maxMu = maximum(mod.mu)

    ## find maximum asset sigma
    ## sigmas = AssetMgmt.getVolas(mod)
    ## maxSigma = maximum(sigmas)

    ## define mu grid
    muGrid = linspace(muGMV, maxMu, nPoints)

    ## get fundamentals
    a, b, c, d = AssetMgmt.getFundamentals(mod)

    ## get sigmas for given mu grid
    effSigmas = sqrt((a*muGrid.^2 - 2*b*muGrid + c)/d)

    return muGrid, effSigmas
end

function getEffPortfolios(mod::MuSigmaModel;
                          nPoints = 50)

    ## get GMV portfolio
    xGMV = AssetMgmt.getGMVSS(mod)

    ## get associated portfolio mu
    muGMV = AssetMgmt.getPMean(xGMV, mod.mu)

    ## find maximum asset mu
    maxMu = maximum(mod.mu)

    ## define mu grid
    muGrid = linspace(muGMV, maxMu, nPoints)

    ## get number of assets
    nAss = size(mod.mu, 1)

    ## preallocation
    nPfs = length(muGrid)
    effWgts = ones(nPfs, nAss)

    for ii=1:nPfs
        muTarget = muGrid[ii]

        ## get weights
        wgts = getEffPfGivenMu(mod, muTarget)
        effWgts[ii, :] = wgts
    end

    return effWgts
end

function getEffPfGivenMu(mod, muTarget)
    covMatr = mod.sigma

    ## get number of assets
    nAss = size(mod.mu, 1)

    ## create helping matrices
    R = [ones(nAss) mod.mu]

    ## get constraints
    muConstr = [1; muTarget]

    ## get weights
    wgts = inv(covMatr)*R*inv(R'*inv(covMatr)*R)*muConstr

    return wgts'
end


text:
# fbp2_back_fan.jl
# obsolete: fbp2_back now covers both fan and par case

export fbp2_back_fan

using LazyGrids: ndgrid

"""
    img = fbp2_back_fan(sg, ig, sino; ia_skip)

2D backprojection for fan-beam FBP.

in
- `sg::SinoGeom`
- `ig::ImageGeom`
- `sino::AbstractArray{<:Number}`      sinogram(s) (line integrals)

options
- `ia_skip::Int`                        downsample in angle to save time for quick tests (default: 1)

out
- `img::AbstractArray{<:Number}`       reconstructed image(s)

"""
function fbp2_back_fan(sg::SinoGeom, ig::ImageGeom, sino::AbstractMatrix{<:Number}; ia_skip::Int=1)

    sg isa SinoFan || throw("need fan type")

    if sg.dfs == 0
        is_arc=true
    elseif isinf(dsf)
        is_arc=false
    else
        throw("bad dsf")
    end

    return fbp2_back_fan(sino, sg.orbit, sg.orbit_start,
    sg.dsd, sg.dso, sg.dfs, sg.ds, sg.offset,
    sg.source_offset,
    ig.nx, ig.ny, ig.dx, ig.dy, ig.offset_x, ig.offset_y,
    is_arc, ig.mask, ia_skip)

end

function fbp2_back_fan(sino::AbstractMatrix{<:Number}, orbit::Union{Symbol,Real}, orbit_start::Real,
    dsd::RealU, dso::Real, dfs::RealU, ds::RealU, offset::Real, source_offset::Real,
    nx::Int, ny::Int, dx::RealU, dy::RealU, offset_x::Real, offset_y::Real,
     is_arc::Bool, mask::AbstractMatrix{Bool}, ia_skip::Int,
)
    rmax=[]

    na,nb=size(sino)

    # trick: extra zero column saves linear interpolation indexing within loop!
    sino = [sino;zeros(eltype(sino),size(sino,2),2)']

    # precompute as much as possible
    wx = (nx+1)/2 - offset_x
    wy = (ny+1)/2 - offset_y
    xc, yc = ndgrid(dx * ((1:nx) .- wx), dy * ((1:ny) .- wy))
    rr = @.(sqrt(abs2(xc) + abs2(yc))) # [nx,ny]

    smax = ((nb-1)/2 - abs(offset)) * ds

    #todo: rmax possibly extracable from SinoGeom
    if is_arc
        gamma_max = smax / dsd
    else # flat
        gamma_max = atan(smax / dsd)
    end
    rmax = dso * sin(gamma_max)

    mask = mask .& (rr .< rmax)
    xc = xc[vec(mask)] # [np] pixels within mask
    yc = yc[vec(mask)]
    #clear wx wy rr smax

    betas = @.(deg2rad(orbit_start + orbit * (0:na-1) / na)) # [na]
    wb = (nb+1)/2 + offset

    img = 0

    for ia=1:ia_skip:na
    #ticker(mfilename, ia, na)

        beta = betas[ia]
        d_loop = @.(dso + xc * sin(beta) - yc * cos(beta)) # dso - y_beta
        r_loop = @.(xc * cos(beta) + yc * sin(beta) - source_offset) # x_beta-roff

        if is_arc
            sprime_ds = (dsd/ds) * atan.(r_loop, d_loop) # s' / ds
            w2 = @.(dsd^2 / (abs2(d_loop) + abs2(r_loop))) # [np] image weighting
        else # flat
            mag = dsd ./ d_loop
            sprime_ds = mag .* r_loop ./ ds
            w2 = abs2.(mag) # [np] image-domain weighting
        end

        bb = sprime_ds .+ wb # [np] bin "index"

        # nearest neighbor interpolation:
        #=
    %    ib = round(bb);
    %    if any(ib < 1 | ib > nb), error 'bug', end
    %    % trick: make out-of-sinogram indices point to those extra zeros
    %%    ib(ib < 1 | ib > nb) = nb+1;
    %    img = img + sino(ib, ia) ./ L2;
        =#

        # linear interpolation:
        il = floor.(Int, bb) # left bin
        ir = 1 .+ il # right bin

        # deal with truncated sinograms
        ig = (il .>= 1) .& (ir .<= nb)
        il[.!ig] .= nb+1
        ir[.!ig] .= nb+1
    #    if any(il < 1 | il >= nb), error 'bug', end

    # Temporary:
        il=min.(il,na)
        ir=min.(ir,na)

        wr = bb .- il # left weight
        wl = 1 .- wr # right weight

        img = @.(img + (wl * sino[il, ia] + wr * sino[ir, ia]) * w2)
    end

    return π / (na/ia_skip) * embed(img,mask)
end

function fbp2_back_fan(sg::SinoGeom, ig::ImageGeom, sino::AbstractArray{<:Number}; kwargs...)
    return mapslices(sino -> fbp2_back_fan(sg, ig, sino; kwargs...), sino, [1,2])
end


text:
using Agents, LinearAlgebra

"""
Wrapper function which updates agent velocities at each simulation time step. The function modifies the model in place so this function doesn't return anything.

The wrapper function determines which pairs are proximal to each other based on model parameters then computes the updated velocity for each agent accordingly.
"""
function FMP(model::AgentBasedModel)

    # get list of interacting_pairs within some radius
    agent_iter = interacting_pairs(model, model.r, :all)

    # construct interaction_array which is (num_agents x num_agents)
    #   array where interaction_array[i,j] = 1 implies that
    #   agent_i and agent_j are within the specified interaction radius
    interaction_array = falses( nagents(model), nagents(model))
    agents = agent_iter.agents
    for pair in agent_iter.pairs

        i, j = pair
        if agents[i].type == :A && agents[j].type == :A
            interaction_array[i, j] = true
            interaction_array[j, i] = true
        end

    end

    # determine object ids

    # loop through agents and update velocities
   for i in keys(agents)
        Ni = findall(x->x==1, interaction_array[i, :])
        # move_this_agent_to_new_position(i) in FMP paper
        UpdateVelocity(model, i, Ni, agents)
    end
end

"""
Sub-wrapper function that updates the agent velocity based on the velocity subcomponents described in the FMP paper. The three forces included are:

- Repulsive force: analogous to a "magnetic" repulsion based on proximity of agent and other agents in the state space.
- Navigational force: an attractive force drawing an agent to its goal position
- Obstacle force: similar to repulsive force but generated by proximity to objects in the state space

After computing the resultant vectors from each component, an overall resultant vector is computed. This is then capped based on the global max velocity constraint.
"""
function UpdateVelocity(model::AgentBasedModel, i, Ni, agents)

    # compute forces and resultant velocities
    fiR = RepulsiveForce(model, agents, i, Ni)
    fiGamma = NavigationalFeedback(model, agents, i)
    fiObject = ObstactleFeedback(model, agents, i)
    ui = fiR .+ fiGamma .+ fiObject
    vi = agents[i].vel .+ ui .* model.dt
    vi = CapVelocity(model.vmax, vi)

    # update agent velocities
    agents[i].vel = vi

end

"""
Function to calculate the resultant velocity vector from the repulsive component of the FMP
algorithm.
"""
function RepulsiveForce(model::AgentBasedModel, agents, i, Ni)
    # compute repulsive force for each agent
    # note the "." before most math operations, required for component wise tuple math
    f = ntuple(i->0, length(agents[i].vel))
    for j in Ni
        dist = norm(agents[j].pos .- agents[i].pos)
        if dist < model.r
            force = -model.rho * (dist - model.r)^2
            distnorm = (agents[j].pos .- agents[i].pos) ./dist
            f = f .+ (force .* distnorm)
        end
    end

    # targets/objects do not experience repulsive feedback
    if agents[i].type == :O || agents[i].type == :T
        return  ntuple(i->0, length(agents[i].vel))
    else
        return f
    end
end

"""
Function to calculate the resultant velocity vector from the navigational component of the FMP
algorithm.
"""
function NavigationalFeedback(model::AgentBasedModel, agents, i)
    # compute navigational force for each agent
    # note the "." before most math operations, required for component wise tuple math
    f = (-model.c1 .* (agents[i].pos .- agents[i].tau)) .+ (- model.c2 .* agents[i].vel)
    if agents[i].type == :T
        return  ntuple(i->0, length(agents[i].vel))  # targets to not experience navigational feedback
    else
        return f
    end
end

"""
Function to calculate the resultant velocity vector from the obstacle avoidance component
of the FMP algorithm.
"""
function ObstactleFeedback(model::AgentBasedModel, agents, i)
    # determine obstacle avoidance feedback term
    # note the "." before most math operations, required for component wise tuple math
    f = ntuple(i->0, length(agents[i].vel))

    for id in model.obstacle_list
        # the original paper defines z as p_j-r_j-p_i in equation 17/18
        #   in the paper r_j is treated a vector, however it makes more sense to
        #   treat as a scalar quantity so we take the norm, then subtract the radius
        #   (j is obstacle (id) and i is agent (i))
        dist = norm(agents[id].pos  .- agents[i].pos) - agents[id].radius
        if dist < agents[i].radius
            force = -model.rho_obstacle * (dist - agents[id].radius)^2
            distnorm = (agents[id].pos .- agents[i].pos) ./ norm(agents[id].pos .- agents[i].pos)
            f = f .+ (force .* distnorm)
        end
    end
    if agents[i].type == :O || agents[i].type == :T
        return ntuple(i->0, length(agents[i].vel))
    else
        return f
    end
end

"""
Function to bound computed velocities based on globally set vmax parameter.
"""
function CapVelocity(vmax, vel)
    # bound velocity by vmax
    # note the "." before most math operations, required for component wise tuple math
    if norm(vel) > vmax
        vi = (vel ./ norm(vel)) .* vmax
        return vi
    else
        return vel
    end
end

max_stars_repo_name: rmcsqrd/antfarm

max_issues_repo_name: rmcsqrd/antfarm

max_forks_repo_name: rmcsqrd/antfarm