text:
# You descend into the ocean trench and encounter some snailfish. They say they
# saw the sleigh keys! They'll even tell you which direction the keys went if
# you help one of the smaller snailfish with his math homework.
#
# Snailfish numbers aren't like regular numbers. Instead, every snailfish
# number is a pair - an ordered list of two elements. Each element of the pair
# can be either a regular number or another pair.
#
# Pairs are written as [x,y], where x and y are the elements within the pair.
# Here are some example snailfish numbers, one snailfish number per line:
#
# [1,2]
# [[1,2],3]
# [9,[8,7]]
# [[1,9],[8,5]]
# [[[[1,2],[3,4]],[[5,6],[7,8]]],9]
# [[[9,[3,8]],[[0,9],6]],[[[3,7],[4,9]],3]]
# [[[[1,3],[5,3]],[[1,3],[8,7]]],[[[4,9],[6,9]],[[8,2],[7,3]]]]
#
# This snailfish homework is about addition. To add two snailfish numbers, form
# a pair from the left and right parameters of the addition operator. For
# example, [1,2] + [[3,4],5] becomes [[1,2],[[3,4],5]].
#
# There's only one problem: snailfish numbers must always be reduced, and the
# process of adding two snailfish numbers can result in snailfish numbers that
# need to be reduced.
#
# To reduce a snailfish number, you must repeatedly do the first action in this
# list that applies to the snailfish number:
#
# - If any pair is nested inside four pairs, the leftmost such pair explodes.
# - If any regular number is 10 or greater, the leftmost such regular number
#   splits.
#
# Once no action in the above list applies, the snailfish number is reduced.
#
# During reduction, at most one action applies, after which the process returns
# to the top of the list of actions. For example, if split produces a pair that
# meets the explode criteria, that pair explodes before other splits occur.
#
# To explode a pair, the pair's left value is added to the first regular number
# to the left of the exploding pair (if any), and the pair's right value is
# added to the first regular number to the right of the exploding pair (if
# any). Exploding pairs will always consist of two regular numbers. Then, the
# entire exploding pair is replaced with the regular number 0.
#
# Here are some examples of a single explode action:
#
# - [[[[[9,8],1],2],3],4] becomes [[[[0,9],2],3],4] (the 9 has no regular
#   number to its left, so it is not added to any regular number).
# - [7,[6,[5,[4,[3,2]]]]] becomes [7,[6,[5,[7,0]]]] (the 2 has no regular
#   number to its right, and so it is not added to any regular number).
# - [[6,[5,[4,[3,2]]]],1] becomes [[6,[5,[7,0]]],3].
# - [[3,[2,[1,[7,3]]]],[6,[5,[4,[3,2]]]]] becomes
#   [[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]] (the pair [3,2] is unaffected because the
#   pair [7,3] is further to the left; [3,2] would explode on the next action).
# - [[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]] becomes [[3,[2,[8,0]]],[9,[5,[7,0]]]].
#
# To split a regular number, replace it with a pair; the left element of the
# pair should be the regular number divided by two and rounded down, while the
# right element of the pair should be the regular number divided by two and
# rounded up. For example, 10 becomes [5,5], 11 becomes [5,6], 12 becomes
# [6,6], and so on.
#
# Here is the process of finding the reduced result of
# [[[[4,3],4],4],[7,[[8,4],9]]] + [1,1]:
#
# after addition: [[[[[4,3],4],4],[7,[[8,4],9]]],[1,1]]
# after explode:  [[[[0,7],4],[7,[[8,4],9]]],[1,1]]
# after explode:  [[[[0,7],4],[15,[0,13]]],[1,1]]
# after split:    [[[[0,7],4],[[7,8],[0,13]]],[1,1]]
# after split:    [[[[0,7],4],[[7,8],[0,[6,7]]]],[1,1]]
# after explode:  [[[[0,7],4],[[7,8],[6,0]]],[8,1]]
#
# Once no reduce actions apply, the snailfish number that remains is the actual
# result of the addition operation: [[[[0,7],4],[[7,8],[6,0]]],[8,1]].
#
# The homework assignment involves adding up a list of snailfish numbers (your
# puzzle input). The snailfish numbers are each listed on a separate line. Add
# the first snailfish number and the second, then add that result and the
# third, then add that result and the fourth, and so on until all numbers in
# the list have been used once.
#
# For example, the final sum of this list is [[[[1,1],[2,2]],[3,3]],[4,4]]:
#
# [1,1]
# [2,2]
# [3,3]
# [4,4]
#
# The final sum of this list is [[[[3,0],[5,3]],[4,4]],[5,5]]:
#
# [1,1]
# [2,2]
# [3,3]
# [4,4]
# [5,5]
#
# The final sum of this list is [[[[5,0],[7,4]],[5,5]],[6,6]]:
#
# [1,1]
# [2,2]
# [3,3]
# [4,4]
# [5,5]
# [6,6]
#
# Here's a slightly larger example:
#
# [[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]]
# [7,[[[3,7],[4,3]],[[6,3],[8,8]]]]
# [[2,[[0,8],[3,4]]],[[[6,7],1],[7,[1,6]]]]
# [[[[2,4],7],[6,[0,5]]],[[[6,8],[2,8]],[[2,1],[4,5]]]]
# [7,[5,[[3,8],[1,4]]]]
# [[2,[2,2]],[8,[8,1]]]
# [2,9]
# [1,[[[9,3],9],[[9,0],[0,7]]]]
# [[[5,[7,4]],7],1]
# [[[[4,2],2],6],[8,7]]
#
# The final sum [[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]] is found
# after adding up the above snailfish numbers:
#
#   [[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]]
# + [7,[[[3,7],[4,3]],[[6,3],[8,8]]]]
# = [[[[4,0],[5,4]],[[7,7],[6,0]]],[[8,[7,7]],[[7,9],[5,0]]]]
#
#   [[[[4,0],[5,4]],[[7,7],[6,0]]],[[8,[7,7]],[[7,9],[5,0]]]]
# + [[2,[[0,8],[3,4]]],[[[6,7],1],[7,[1,6]]]]
# = [[[[6,7],[6,7]],[[7,7],[0,7]]],[[[8,7],[7,7]],[[8,8],[8,0]]]]
#
#   [[[[6,7],[6,7]],[[7,7],[0,7]]],[[[8,7],[7,7]],[[8,8],[8,0]]]]
# + [[[[2,4],7],[6,[0,5]]],[[[6,8],[2,8]],[[2,1],[4,5]]]]
# = [[[[7,0],[7,7]],[[7,7],[7,8]]],[[[7,7],[8,8]],[[7,7],[8,7]]]]
#
#   [[[[7,0],[7,7]],[[7,7],[7,8]]],[[[7,7],[8,8]],[[7,7],[8,7]]]]
# + [7,[5,[[3,8],[1,4]]]]
# = [[[[7,7],[7,8]],[[9,5],[8,7]]],[[[6,8],[0,8]],[[9,9],[9,0]]]]
#
#   [[[[7,7],[7,8]],[[9,5],[8,7]]],[[[6,8],[0,8]],[[9,9],[9,0]]]]
# + [[2,[2,2]],[8,[8,1]]]
# = [[[[6,6],[6,6]],[[6,0],[6,7]]],[[[7,7],[8,9]],[8,[8,1]]]]
#
#   [[[[6,6],[6,6]],[[6,0],[6,7]]],[[[7,7],[8,9]],[8,[8,1]]]]
# + [2,9]
# = [[[[6,6],[7,7]],[[0,7],[7,7]]],[[[5,5],[5,6]],9]]
#
#   [[[[6,6],[7,7]],[[0,7],[7,7]]],[[[5,5],[5,6]],9]]
# + [1,[[[9,3],9],[[9,0],[0,7]]]]
# = [[[[7,8],[6,7]],[[6,8],[0,8]]],[[[7,7],[5,0]],[[5,5],[5,6]]]]
#
#   [[[[7,8],[6,7]],[[6,8],[0,8]]],[[[7,7],[5,0]],[[5,5],[5,6]]]]
# + [[[5,[7,4]],7],1]
# = [[[[7,7],[7,7]],[[8,7],[8,7]]],[[[7,0],[7,7]],9]]
#
#   [[[[7,7],[7,7]],[[8,7],[8,7]]],[[[7,0],[7,7]],9]]
# + [[[[4,2],2],6],[8,7]]
# = [[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]]
#
# To check whether it's the right answer, the snailfish teacher only checks the
# magnitude of the final sum. The magnitude of a pair is 3 times the magnitude
# of its left element plus 2 times the magnitude of its right element. The
# magnitude of a regular number is just that number.
#
# For example, the magnitude of [9,1] is 3*9 + 2*1 = 29; the magnitude of [1,9]
# is 3*1 + 2*9 = 21. Magnitude calculations are recursive: the magnitude of
# [[9,1],[1,9]] is 3*29 + 2*21 = 129.
#
# Here are a few more magnitude examples:
#
# [[1,2],[[3,4],5]] becomes 143.
# [[[[0,7],4],[[7,8],[6,0]]],[8,1]] becomes 1384.
# [[[[1,1],[2,2]],[3,3]],[4,4]] becomes 445.
# [[[[3,0],[5,3]],[4,4]],[5,5]] becomes 791.
# [[[[5,0],[7,4]],[5,5]],[6,6]] becomes 1137.
# [[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]] becomes 3488.
#
# So, given this example homework assignment:
#
# [[[0,[5,8]],[[1,7],[9,6]]],[[4,[1,2]],[[1,4],2]]]
# [[[5,[2,8]],4],[5,[[9,9],0]]]
# [6,[[[6,2],[5,6]],[[7,6],[4,7]]]]
# [[[6,[0,7]],[0,9]],[4,[9,[9,0]]]]
# [[[7,[6,4]],[3,[1,3]]],[[[5,5],1],9]]
# [[6,[[7,3],[3,2]]],[[[3,8],[5,7]],4]]
# [[[[5,4],[7,7]],8],[[8,3],8]]
# [[9,3],[[9,9],[6,[4,9]]]]
# [[2,[[7,7],7]],[[5,8],[[9,3],[0,2]]]]
# [[[[5,2],5],[8,[3,7]]],[[5,[7,5]],[4,4]]]
#
# The final sum is:
#
# [[[[6,6],[7,6]],[[7,7],[7,0]]],[[[7,7],[7,7]],[[7,8],[9,9]]]]
#
# The magnitude of this final sum is 4140.
#
# Add up all of the snailfish numbers from the homework assignment in the order
# they appear. What is the magnitude of the final sum?

# read input
lines = map(eachline("input.txt")) do l
  r = []
  stack = [r]

  # first character will always be an opening bracket - skip it
  for c in l[2:end]
    # skip commas - numbers are always 1 digit in input, so we don't need to do
    # anything special with commas
    c == ',' && continue
    if c == '['
      nextpair = []
      push!(stack[end], nextpair)
      push!(stack, nextpair)
    elseif c == ']'
      pop!(stack)
    else
      # must be a number
      push!(stack[end], c - '0')
    end
  end
  return r
end

# utility function to handle explosions
function do_explosion(l)
  last_number = nothing
  add_to_next = nothing
  has_explosion = false
  exploder = function (p, depth)
    has_explosion && isnothing(add_to_next) && return p
    if depth == 5 && !has_explosion
      if !isnothing(last_number)
        last_number[1][last_number[2]] += p[1]
      end
      add_to_next = p[2]
      has_explosion = true
      return 0
    end

    if isa(p[1], Int)
      if isnothing(add_to_next)
        last_number = (p, 1)
      else
        p[1] += add_to_next
        add_to_next = nothing
        return p
      end
    else
      p[1] = exploder(p[1], depth + 1)
    end
    if isa(p[2], Int)
      if isnothing(add_to_next)
        last_number = (p, 2)
      else
        p[2] += add_to_next
        add_to_next = nothing
        return p
      end
    else
      p[2] = exploder(p[2], depth + 1)
    end
    return p
  end
  exploder(l, 1)
  return has_explosion
end

# utility function for splits
function do_split(l)
  had_split = false
  splitter = function (p)
    had_split && return p
    if isa(p, Int)
      if p >= 10
        (d,r) = divrem(p, 2)
        had_split = true
        return Any[d, d+r]
      end
      return p
    end
    p[1] = splitter(p[1])
    p[2] = splitter(p[2])
    return p
  end
  splitter(l)
  return had_split
end

# sum
total_sum = foldl(lines) do a,b
  r = Any[a, b]
  while true
    do_explosion(r) || do_split(r) || break
  end
  return r
end

# find magnitude
function magnitude(p)
  isa(p, Int) && return p
  3 * magnitude(p[1]) + 2 * magnitude(p[2])
end

println("Result: ", magnitude(total_sum))


text:
@testset "Fitting data" begin
	x = collect(1.:1:10*365)./365;
	y = 1869.9 .- 782.0 .* exp.(-0.085 .* x);
	par,er = fitexp(x,y);
	@test round(sum(par*10)) ≈ round(18699-7820-0.85)
	@test sum(er) < 1e-2

	@test evalexp([3.],[10.,0.5,0.05]) ≈ [10. + 0.5*exp(0.05*3.)]

	x = collect(1.:1:365);
	y = 10.0 .+ 0.1 .* x + rand(length(x))./20;
	par,er = fitpoly(x,y,deg=1);
	@test par[1] ≈ 0.1 atol=0.001
	@test par[2] ≈ 10. atol=0.1
	@test sum(er) < 1e-2

	@test evalpoly([10.],[0.01, 0.1, 1.0]) ≈ [1+0.1*10+0.01*10*10]

	@test evalpoly([DateTime(1,1,1,0,1)],[0.01, 0.1, 1.0]) ≈ [1+0.1*86460000+0.01*86460000.0^2]

	## Test model creation
	m = ResampleAndFit.getpolymodel(0)
	@test m([10.,20.],2.0) == [2.0, 2.0]
	m = ResampleAndFit.getpolymodel(3)
	@test m(1,[4.,3.,2.,1]) == sum([4.,3.,2.,1])
	m = ResampleAndFit.getpolymodel(4)
	@test m(1,[5.,4.,3.,2.,1]) == sum([5.,4.,3.,2.,1])
	@test isnan(ResampleAndFit.getpolymodel(10))
end

ext: jl


text:
"""
    NodeType

Each node in the directed graph can be classified into the following types
- VARIABLE: Denotes a decision variable.
- PARAMETER: An adjustable parameter value (not a decision variable).
- CONSTANT: A constant value
- EXPRESSION: Any other expression that isn't a subexpression
- SUBEXPRESSION: Any expression referencing a different graph representation.
"""
@enum(NodeClass, CONSTANT, PARAMETER, VARIABLE, EXPRESSION, SUBEXPRESSION)
#=
- SELECT: These nodes do not store values rather they reference the output
          stored in the buffer assocatied with an EXPRESSION. These nodes
          only occur after a multiple output function. In the mimo JuMP
          extension they correspond to intermediate variables introduced
          by the user.
=#

struct Variable end
struct Parameter end
struct Constant end
struct Expression end
struct Subexpression end
struct User end

function user end
function usern end

const ATOM_TYPE_INSTANCES = instances(AtomType)

abstract type AbstractNode end

"""
Describes connectivity and expression represented by node.
"""
struct Node <: AbstractNode
    node_class::NodeClass
    ex_type::AtomType
    first_index::Int
    second_index::Int
    arity::Int
    children::Vector{Int}
end

for (t, s, a) in ((Variable, VARIABLE, VAR_ATOM),
                  (Parameter, PARAMETER, PARAM_ATOM),
                  (Constant, CONSTANT, CONST_ATOM),
                  (Subexpression, SUBEXPRESSION, SUBEXPR),)
    @eval Node(::$t, i::Int) = Node($s, $a, i, 0, 0, Int[])
end

for v in (PLUS, MINUS, MULT, POW, DIV, MAX, MIN)
    @eval Node(::Val{true}, ::Val{$v}, c::Vector{Int}) = Node(EXPRESSION, $v, 0, 0, length(c), c)
end
for d in ALL_ATOM_TYPES
    @eval Node(::Val{false}, ::Val{$d}, c::Vector{Int}) = Node(EXPRESSION, $d, 0, 0, 1, c)
end
Node(::Val{true}, ::Val{USER}, i::Int, c::Vector{Int}) = Node(EXPRESSION, USER, i, 0, 1, c)
Node(::Val{true}, ::Val{USERN}, i::Int, c::Vector{Int}) = Node(EXPRESSION, USERN, i, 0, length(c), c)

Node(::Val{:first_index}, n::Node, i::Int) = Node(node_class(n), ex_type(n), i, second_index(n), arity(n), children(n))

node_class(n::Node)   = n.node_class
ex_type(n::Node)      = n.ex_type
first_index(n::Node)  = n.first_index
second_index(n::Node) = n.node_second_index
arity(n::Node)        = n.arity
children(n::Node)     = n.children
child(n::Node, i)     = @inbounds getindex(n.children, i)

node_is_class(::Variable, n::Node) = node_class(n) == VARIABLE
node_is_class(::Parameter, n::Node) = node_class(n) == PARAMETER
node_is_class(::Constant, n::Node) = node_class(n) == CONSTANT

mv_eago_not_jump = setdiff(JuMP._Derivatives.operators,
                           union(Symbol[k for k in keys(REV_BIVARIATE_ATOM_DICT)],
                                 Symbol[k for k in keys(REV_NARITY_ATOM_DICT)]))
eago_mv_switch = quote end
for s in mv_eago_not_jump
    global eago_mv_switch = quote
        $eago_mv_switch
        (d == $s) && (return Node(Val(true), Val($s), v[c]))
    end
end
@eval function _create_call_node(i, v, c::UnitRange{Int}, op::OperatorRegistry)
    if i == 1
        return Node(Val(true), Val(PLUS), v[c])
    elseif i == 2
        return Node(Val(true), Val(MINUS), v[c])
    elseif i == 3
        return Node(Val(true), Val(MULT), v[c])
    elseif i == 4
        return Node(Val(true), Val(POW), v[c])
    elseif i == 5
        return Node(Val(true), Val(DIV), v[c])
    elseif i == 6
        error("If-else currently unsupported...")
    elseif i == 7
        return Node(Val(true), Val(MAX), v[c])
    elseif i == 8
        return Node(Val(true), Val(MIN), v[c])
    elseif i >= JuMP._Derivatives.USER_OPERATOR_ID_START
        i_mv = i - JuMP._Derivatives.USER_OPERATOR_ID_START + 1
        d = op.multivariate_id[i_mv]
        $eago_mv_switch
        return Node(Val(true), Val(USERN), i_mv, v[c])
    end
end

function binary_switch_typ(ids, exprs)
    if length(exprs) <= 3
        out = Expr(:if, Expr(:call, :(==), :i, ids[1]),
                   :(Node(Val(false), Val($(exprs[1])), v[c])))
        if length(exprs) > 1
            push!(out.args, binary_switch_typ(ids[2:end], exprs[2:end]))
        end
        return out
    else
        mid = length(exprs) >>> 1
        return Expr(:if, Expr(:call, :(<=), :i, ids[mid]),
                         binary_switch_typ(ids[1:mid], exprs[1:mid]),
                         binary_switch_typ(ids[mid+1:end], exprs[mid+1:end]))
    end
end

indx_JuMP = Int[]
indx_EAGO = AtomType[]
for k in univariate_operators
    if haskey(REV_UNIVARIATE_ATOM_DICT, k)
        k_EAGO = REV_UNIVARIATE_ATOM_DICT[k]
        push!(indx_JuMP, univariate_operator_to_id[k])
        push!(indx_EAGO, k_EAGO)
    end
end

uni_eago_not_jump = setdiff(univariate_operators, Symbol[k for k in keys(REV_UNIVARIATE_ATOM_DICT)])
uni_eago_not_jump = push!(uni_eago_not_jump, :-)
eago_uni_switch = quote end
for s in uni_eago_not_jump
    global eago_uni_switch = quote
        $eago_uni_switch
        (d == $s) && (return Node(Val(false), Val($s), v[c]))
    end
end
atom_switch = binary_switch_typ(indx_JuMP, indx_EAGO)
@eval function _create_call_node_uni(i::Int, v, c::UnitRange{Int}, op::OperatorRegistry)

    if i >= JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START
        j = i - JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START + 1
        dop = op.univariate_operator_id[j]
        d = op.univariate_operator_to_id[dop]
        $eago_uni_switch
        return Node(Val(true), Val(USER), j, v[c])
    end
    $atom_switch
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function set_wallpoint(nodes, bdcon, cellxmax, cellymax, cellzmax)
	# cal number of ponits
	nop = 0
	if Int(bdcon[1][1]) == 2 || Int(bdcon[1][1]) == 3
		nop = nop + (cellymax-1)*(cellzmax-1)
	end
	if Int(bdcon[2][1]) == 2 || Int(bdcon[2][1]) == 3
		nop = nop + (cellymax-1)*(cellzmax-1)
	end
	if Int(bdcon[3][1]) == 2 || Int(bdcon[3][1]) == 3
		nop = nop + (cellzmax-1)*(cellxmax-1)
	end
	if Int(bdcon[4][1]) == 2 || Int(bdcon[4][1]) == 3
		nop = nop + (cellzmax-1)*(cellxmax-1)
	end
	if Int(bdcon[5][1]) == 2 || Int(bdcon[5][1]) == 3
		nop = nop + (cellxmax-1)*(cellymax-1)
	end
	if Int(bdcon[6][1]) == 2 || Int(bdcon[6][1]) == 3
		nop = nop + (cellxmax-1)*(cellymax-1)
	end

	wallpoint = zeros(nop,3)
	ite = 1

	# x-
	if Int(bdcon[1][1]) == 2 || Int(bdcon[1][1]) == 3
		for j in 2:cellymax
			for k in 2:cellzmax
				wallpoint[ite,1] = nodes[2,j,k,1]
				wallpoint[ite,2] = nodes[2,j,k,2]
				wallpoint[ite,3] = nodes[2,j,k,3]
				ite = ite + 1
			end
		end
	end
	# x+
	if Int(bdcon[2][1]) == 2 || Int(bdcon[2][1]) == 3
		for j in 2:cellymax
			for k in 2:cellzmax
				wallpoint[ite,1] = nodes[cellzmax,j,k,1]
				wallpoint[ite,2] = nodes[cellzmax,j,k,2]
				wallpoint[ite,3] = nodes[cellzmax,j,k,3]
				ite = ite + 1
			end
		end
	end
	# y-
	if Int(bdcon[3][1]) == 2 || Int(bdcon[3][1]) == 3
		for k in 2:cellzmax
			for i in 2:cellxmax
				wallpoint[ite,1] = nodes[i,2,k,1]
				wallpoint[ite,2] = nodes[i,2,k,2]
				wallpoint[ite,3] = nodes[i,2,k,3]
				ite = ite + 1
			end
		end
	end
	# y+
	if Int(bdcon[4][1]) == 2 || Int(bdcon[4][1]) == 3
		for k in 2:cellzmax
			for i in 2:cellxmax
				wallpoint[ite,1] = nodes[i,cellymax,k,1]
				wallpoint[ite,2] = nodes[i,cellymax,k,2]
				wallpoint[ite,3] = nodes[i,cellymax,k,3]
				ite = ite + 1
			end
		end
	end
	# z-
	if Int(bdcon[5][1]) == 2 || Int(bdcon[5][1]) == 3
		for i in 2:cellxmax
			for j in 2:cellymax
				wallpoint[ite,1] = nodes[i,j,2,1]
				wallpoint[ite,2] = nodes[i,j,2,2]
				wallpoint[ite,3] = nodes[i,j,2,3]
				ite = ite + 1
			end
		end
	end
	# z+
	if Int(bdcon[6][1]) == 2 || Int(bdcon[6][1]) == 3
		for i in 2:cellxmax
			for j in 2:cellymax
				wallpoint[ite,1] = nodes[i,j,cellzmax,1]
				wallpoint[ite,2] = nodes[i,j,cellzmax,2]
				wallpoint[ite,3] = nodes[i,j,cellzmax,3]
				ite = ite + 1
			end
		end
	end
	return wallpoint, nop
end

function set_wally(cellcenter, wallpoint, nop, cellxmax, cellymax, cellzmax)
	wally = zeros(cellxmax, cellymax, cellzmax)

	distance_list = zeros(nop,2)
	distance = 0

	pickup_point = 5

	# 最大値を5つ抽出
	for i in 2:cellxmax-1
		for j in 2:cellymax-1
			for k in 2:cellzmax-1
				for n in 1:nop
						distance = ((cellcenter[i,j,k,1]-wallpoint[n,1])^2 +
									(cellcenter[i,j,k,2]-wallpoint[n,2])^2 +
									(cellcenter[i,j,k,3]-wallpoint[n,3])^2)^0.5
						distance_list[n,1] = n
						distance_list[n,2] = distance
				end

				qdistance_list = copy(distance_list)
				qdistance_list = quicksort(qdistance_list,1,nop,2,2)

				#println(distance_list[1,2])

				# 5wall_point
				wp = zeros(Int, pickup_point)
				for m in 1:pickup_point
					wp[m] = Int(qdistance_list[m,1])
				end

				# max(y)
				# 近傍5点を用いて一番遠いyが正解
				# 例1 1点はy+側、2点はy-側のとき、y-とy+で作れる三角形の高さが一番短いが間違い
				# 逆に一番遠いものが正解
				ytemp = zeros(10)
				ite1 = [1,1,1,1,1,1,2,2,2,3]
				ite2 = [2,2,2,3,3,4,3,3,4,4]
				ite3 = [3,4,5,4,5,5,4,5,5,5]

				# println(qdistance_list[5])
				for m in 1:10
					wp1 = wp[ite1[m]]
					wp2 = wp[ite2[m]]
					wp3 = wp[ite3[m]]

					dis1 = distance_list[wp1,2]
					dis2 = distance_list[wp2,2]
					dis3 = distance_list[wp3,2]

					wallpoint_dis1 = ((wallpoint[wp1,1]-wallpoint[wp2,1])^2 + (wallpoint[wp1,2]-wallpoint[wp2,2])^2 + (wallpoint[wp1,3]-wallpoint[wp2,3])^2)^0.5
					wallpoint_dis2 = ((wallpoint[wp2,1]-wallpoint[wp3,1])^2 + (wallpoint[wp2,2]-wallpoint[wp3,2])^2 + (wallpoint[wp2,3]-wallpoint[wp3,3])^2)^0.5
					wallpoint_dis3 = ((wallpoint[wp3,1]-wallpoint[wp1,1])^2 + (wallpoint[wp3,2]-wallpoint[wp1,2])^2 + (wallpoint[wp3,3]-wallpoint[wp1,3])^2)^0.5

					V = cal_tetrahedron_volume(dis1, dis2, dis3, wallpoint_dis1, wallpoint_dis2, wallpoint_dis3)

					S = cal_area(wallpoint_dis1, wallpoint_dis2, wallpoint_dis3)

					ytemp[m] = V/S
				end
				wally[i,j,k] = maximum(ytemp)
			end
		end
	end

	#=
	・壁面に下した垂線の計算について
	面積S = (s(s-a)(s-b)(s-c))^0.5 = hc/2  , h=垂線,c=底面(壁面の2点間距離)
	これで三辺の長さがわかれば垂線が算出できる
	=#
	return wally
end

function cal_tetrahedron_volume(a1,a2,a3,a4,a5,a6)
	temp1 = a1^2 * a5^2 * (a2^2 + a3^2 + a4^2 + a6^2 - a1^2 - a5^2)
	temp2 = a2^2 * a6^2 * (a1^2 + a3^2 + a4^2 + a5^2 - a2^2 - a6^2)
	temp3 = a3^2 * a4^2 * (a1^2 + a2^2 + a5^2 + a6^2 - a3^2 - a4^2)
	temp4 = -(a1*a2*a4)^2 - (a2*a3*a5)^2 - (a1*a3*a6)^2 - (a4*a5*a6)^2
	V2 = 1/144 * (temp1 + temp2 + temp3 + temp4)

	# ---------------------------------
	# 5つの近傍点で計算するため，構造格子だと，xが同じ値の場合がありエラーになるから
	# ---------------------------------
	volume = 0.0
	try volume = (V2)^0.5
	catch

	end

	return volume
end

function cal_area(a1,a2,a3)
	s = 0.5*(a1 + a2 + a3)

	# ---------------------------------
	# 5つの近傍点で計算するため，構造格子だと，xが同じ値の場合がありエラーになるから
	# ---------------------------------
	surface = s * (s-a1) * (s-a2) * (s-a3)

	if surface <= 0.0
		surface = 1.0e50
	end

	surface = surface^0.5

	return surface
end

function cal_yplus(Qbase, wally, mu, cellxmax, cellymax, cellzmax)
	yplus = zeros(cellxmax, cellymax, cellzmax)

	for i in 2:cellxmax-1
		for j in 2:cellymax-1
			for k in 2:cellzmax-1
				nu = mu[i,j,k] / Qbase[i,j,k,1]
				yplus[i,j,k], uplus = cal_wall_val_spalding(Qbase[i,j,k,2], wally[i,j,k], nu)

				if isequal(yplus[i,j,k], NaN) == true
					println("      ")
					println(Qbase[i,j,k,2])
					println(wally[i,j,k])
					println(nu)
					throw(UndefVarError(:x))
				end
			end
		end
	end

	return yplus
end

function cal_wall_val_spalding(u,y,nu)
	k = 0.4
	B = 5.5

	u_tau = 1.0
	for i in 1:20
		#=
		println("   zzz   ")
		println(u)
		println(y)
		println(nu)
		println(spalding(u_tau,u,y,nu,k,B))
		println(spalding_dash(u_tau,u,y,nu,k,B))
		println(u_tau)
		=#

		old = u_tau
		u_tau = u_tau-spalding(u_tau,u,y,nu,k,B)/spalding_dash(u_tau,u,y,nu,k,B) + 1.0e-50

		if abs(u_tau-old)/u_tau < 10^(-4)
			break
		end
	end

	yplus = u_tau*y/nu
	uplus = u/u_tau
	return yplus,uplus
end

function spalding(u_tau,u,y,nu,k,B)
	t = k*u/u_tau
	F = u/u_tau-u_tau*y/nu+exp(-k*B)*(exp(t)-1-t-t^2/2-t^3/6)
	return F
end

function spalding_dash(u_tau,u,y,nu,k,B)
	t = k*u/u_tau
	Fdash = -u/u_tau^2-y/nu+exp(-k*B)*(t/u_tau)*(-exp(t)+1+t+t^2/2)
	return Fdash
end

# test
#=
uu = 0.0
yy = 0.00014
nunu = 1.5e-5
println(cal_wall_val_spalding(uu,yy,nunu))
=#

function SGSstress()
	# 境界におけるtauが必要
end

function cal_Strain_rate_tensor()

end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
"""
    @hprofile

A macro which
- resets the default `TimerOutputs.get_defaulttimer` to zero
- execute the code block
- print the profiling details

This is useful as a coarse-grained profiling strategy to get a rough idea of
where time is spent. Note that this relies on `TimerOutputs` annotations
manually inserted in the code.
"""
macro hprofile(block)
    return quote
        TimerOutputs.enable_debug_timings(IFGF)
        reset_timer!()
        $(esc(block))
        print_timer()
    end
end

# fast invsqrt code taken from here
# https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-julia-8.html
@inline function invsqrt(x::Float64)
    y = @fastmath Float64(1 / sqrt(Float32(x)))
    # This is a Newton-Raphson iteration.
    return 1.5y - 0.5x * y * (y * y)
end

"""
    cone_domain_size_func(Δs₀::NTuple{N,T},k)

Returns an anonymous function `(node) -> Δs` that computes an appropriate size
`Δs` for the interpolation domain of `node` given an intial size `Δs₀` and a
wavenumber `k`. The function is constructed so as to scale `Δs₀` by the inverse
of the acoustic size of `node`.
"""
function cone_domain_size_func(k, ds)
    if k == 0
        func = (node) -> ds
    else
        # oscillatory case (e.g. Helmholtz, Maxwell)
        # k: wavenumber
        func = (node) -> begin
            bbox = IFGF.container(node)
            w    = maximum(high_corner(bbox) - low_corner(bbox))
            δ    = max(k * w / 2,1)
            ds ./ δ
        end
    end
    return func
end

function cone_domain_size_func(k, ds::Number)
    if k == 0
        func = (node) -> begin
            N = ambient_dimension(node)
            ntuple(i->ds,N)
        end
    else
        # oscillatory case (e.g. Helmholtz, Maxwell)
        # k: wavenumber
        func = (node) -> begin
            N  = ambient_dimension(node)
            ds_tup = ntuple(i->ds,N)
            bbox = IFGF.container(node)
            w    = maximum(high_corner(bbox) - low_corner(bbox))
            δ    = max(k * w / 2,1)
            ds_tup ./ δ
        end
    end
    return func
end

"""
    modified_admissible_condition(target,source,[η])

A target and source are admissible under the *modified admissiblility condition*
(MAC) if the target box lies farther than `r*η` away, where `r` is the radius of
the source box and `η >= 1` is an adjustable parameter. By default, `η = N /
√N`, where `N` is the ambient dimension.
"""
function modified_admissibility_condition(target,source,η)
    # compute distance between source center and target box
    xc = source |> container |> center
    h  = source |> container |> radius
    bbox = container(target)
    dc   = distance(xc, bbox)
    # if target box is outside a sphere of radius h*η, consider it admissible.
    return dc > η*h
end

function modified_admissibility_condition(target,source)
    N = ambient_dimension(target)
    η = N / sqrt(N)
    modified_admissibility_condition(target,source,η)
end

"""
    _density_type_from_kernel_type(T)

Helper function to compute the expected density `V` type associated with a kernel type
`T`. For `T<:Number`, simply return `T`. For `T<:SMatrix`, return an `SVector`
with elements of the type `eltype(T)` and of length `size(T,2)` so that
multiplying an element of type `T` with an element of type `V` makes sense.
"""
function _density_type_from_kernel_type(T)
    if T <: Number
        return T
    elseif T <: SMatrix
        m,n = size(T)
        return SVector{n,eltype(T)}
    else
        error("kernel type $T not recognized")
    end
end

"""
    cheb_error_estimate(coefs::AbstractArray{T,N},dim)

Given an `N` dimensional array of coefficients , estimate the relative error
commited by the Chebyshev interpolant along dimension `dim`.
"""
function cheb_error_estimate(coefs::AbstractArray{T,N},dim) where {T,N}
    sz = size(coefs)
    I  = ntuple(N) do d
        if d == dim
            sz[d]:sz[d]
        else
            1:sz[d]
        end
    end
    norm(view(coefs,I...),2) / norm(coefs,2)
end

"""
    wavenumber(K::Function)

For oscillatory kernels, return the characteristic wavenumber (i.e `2π` divided
by he wavelength). For non-oscillatory kernels, return `0`.
"""
function wavenumber end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
mutable struct Job
    id::JobId
    nworkers::Int32
    sample_rate::Int32
    locations::Dict{ValueId,Location}
    pending_requests::Vector{Request}
    # This is a `WeakKeyDict` so that futures can be GC-ed as long as all
    # references elsewhere are gone.
    futures_on_client::WeakKeyDict{ValueId,Future}
    cluster_name::String
    # To know the current status of the job, check if it is in the global
    # `jobs`. If it is, it's running. If it's not it either failed or
    # completed (`destroy_job` being called doesn't necessarily mean that the
    # job failed).

    # This struct just stores local state for the job.
    function Job(cluster_name::String, job_id::JobId, nworkers::Integer, sample_rate::Integer)::Job
        new(job_id, nworkers, sample_rate, Dict(), [], Dict(), cluster_name)
    end
end

# mutable struct Job
#     job_id::JobId
#     failed::Bool

#     # function Job(; kwargs...)
#     #     new_job_id = create_job(; kwargs...)
#     #     #new_job_id = create_job(;cluster_name="banyancluster", nworkers=2)
#     #     new_job = new(new_job_id)
#     #     finalizer(new_job) do j
#     #         destroy_job(j.job_id)
#     #     end

#     #     new_job
#     # end
# end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using ActuaryUtilities

@testset "ActuaryUtilities.jl integration tests" begin
    cfs = [5, 5, 105]
    times    = [1, 2, 3]

    discount_rates  = [0.03,Yields.Periodic(0.03,1), Yields.Constant(0.03)]

    for d in discount_rates
        @test present_value(d, cfs, times)           ≈ 105.65722270978935
        @test duration(Macaulay(), d, cfs, times)    ≈ 2.86350467067113
        @test duration(d, cfs, times)                ≈ 2.7801016220108057
        @test convexity(d, cfs, times)               ≈ 10.625805482685939
    end
end


text:
"""
Plot example file, and create a matlab file for good measure.  If you are not interested in the
matlab file comment out those sections.
"""

# Pkg.add("Plots")  # uncomment for first run
# Pkg.add("Revise")  # uncomment for first run
# using Revise
using MAT
using HDF5
using Plots
using Gantner

# load file data
name = joinpath(@__DIR__, "data.dat")
(t, time, fs, chanTextLegend) = gantnerread(name)
# t = range(0, length=size(Time, 1), step=1/fs)

plotly()       # use for general purpose plotting and looking for peaks
plot(t, time, label = reshape(chanTextLegend, 1, :))
xlabel!("Time [s]")
ylabel!("Amp")
title!(name)
ylims!(0, 1.2)
xlims!(0,0.3)

# base, ext = split(filename, ".")
# fnamemat = "$base.mat"
# command from .julia/packages/MAT/.../README.md
matwrite(joinpath(@__DIR__, "data.mat"), Dict(
	"t" => collect(t),
	"Time" => time);
	compress = true)

# delete data.mat when done
sleep(10)
rm(joinpath(@__DIR__, "data.mat"))


text:
# # How to plot only specific variables from the results

# See [How to set up plots](3.0_set_up_plots.md) to get started

# ## How to plot a single dataframe

# ```julia
# plots = plot_dataframe(dataframe, time_range)
# ```
# ![this one](plots-54/Generation.png)

# ## How to plot a single dataframe overlayed on an existing plot

# ```julia
# plot = plot_variable(dataframe_1, time_range)
# plot_2 = plot_variable(plot, dataframe_2, time_range)
# ```
# ![this one](plots-55/overlay.png)

# ## How to plot a single variable

# ```julia
# plots = plot_variable(results, "P__ThermalStandard")
# ```
# ![this one](plots-50/P__ThermalStandard.png)

# ## How to plot a single variable overlayed on an existing plot

# ```julia
# plot = plot_variable(results, "P__ThermalStandard")
# plot_2 = plot_variable(plot, results, "P__RenewableDispatch")
# ```
# ![this one](plots-53/overlay.png)

# ## To plot only a couple of variables from the collected results:
# ### Define the variables to be plotted

# ```julia
# variables = [Symbol("P__ThermalStandard")]
# plots = stack_plot(results, variables)
# ```

# ![this one](plots-51/P__ThermalStandard_Stack.png)
# ![this one](plots-51/Plot_with_Fewer_Variables.png)

# This will plot all of the generators for the listed variables.

# ## To collect only a subset of variables and generators:

# ```julia
# selected_variables = Dict(
#   Symbol("P__ThermalStandard") => [:Brighton, :Solitude],
# )
# results_subset = sort_data(results; Variables = selected_variables)
# plots = stack_plot(results_subset)
# ```

# ![this one](plots-52/P__ThermalStandard_Stack.png)
# ![this one](plots-52/Selected_Variables_Plot.png)

# `sort_data()` creates a new results object only containing the subset of variables and generators listed in the dictionary.
# If the Variables key word argument is not called, the default for `sort_data()` is to alphebatize the generators per variable.


text:
const drivers = AG.listdrivers()
const drivermapping = Dict(
    ".shp" => "ESRI Shapefile",
    ".gpkg" => "GPKG",
    ".geojson" => "GeoJSON",
)
const fieldmapping = Dict(v => k for (k, v) in AG._FIELDTYPE)

function read(fn::AbstractString, layer::Union{Integer,AbstractString}=0)
    ds = AG.read(fn)
    layer = AG.getlayer(ds, layer)
    table = AG.Table(layer)
    df = DataFrame(table)
    "" in names(df) && rename!(df, Dict(Symbol("") => :geom, ))  # needed for now
    df
end

function write(fn::AbstractString, table, layer_name::AbstractString="data", geom_column::Symbol=:geom, crs::GFT.GeoFormat=GFT.EPSG(4326), driver::Union{Nothing,AbstractString}=nothing)
    rows = Tables.rows(table)
    sch = Tables.schema(rows)

    # Geom type can't be inferred from here
    geom_type = AG.getgeomtype(rows[1][geom_column])

    # Find driver
    _, extension = splitext(fn)
    if extension in keys(drivermapping)
        driver = AG.getdriver(drivermapping[extension])
    elseif driver !== nothing
        driver = AG.getdriver(driver)
    else
        error("Couldn't determine driver for $extension. Please provide one of $(keys(drivermapping))")
    end

    # Figure out attributes
    fields = Vector{Tuple{Symbol,DataType}}()
    for (name, type) in zip(sch.names, sch.types)
        if type != AG.IGeometry && name != geom_column
            push!(fields, (Symbol(name), type))
        end
    end

    AG.create(
        fn,
        driver=driver
    ) do ds
        AG.createlayer(
            name=layer_name,
            geom=geom_type,
            spatialref=AG.importCRS(crs)
        ) do layer
            for (name, type) in fields
                AG.addfielddefn!(layer, String(name), fieldmapping[type])
            end
            for row in rows
                AG.createfeature(layer) do feature
                    AG.setgeom!(feature, getproperty(row, geom_column))
                    for (name, type) in fields
                        AG.setfield!(feature, AG.findfieldindex(feature, name), getproperty(row, name))
                    end
                end
            end
            AG.copy(layer, dataset=ds, name=layer_name)
        end
    end
    fn
end


text:
@testset "1095.find-in-mountain-array.jl" begin
    @test find_in_mountain_array(MountainArray([1, 2, 3, 4, 5, 3, 1]), 3) == 3
    @test find_in_mountain_array(MountainArray([0, 1, 2, 4, 2, 1]), 3) == -1
    @test find_in_mountain_array(MountainArray([1, 2, 3, 5, 4, 3, 1]), 4) == 5
end


text:
require("covar_matrix_self_adaptation_es.jl")

# Calls to ensure everything has been compiled on all processors.
@everywhere cmsa_es(2, rosenbrock; max_fevals = 2, trace = false)
@everywhere cmsa_es(2, rosenbrock; max_fevals = 2, covarMatrixSampler = CholeskyCovarSampler, trace = false)

# Parameter study
num_runs_per_dim = 25
dims = 2.^(1:6)
# Max evals per dim
max_evals_per_dim = [1000, 10000, 100000]
#max_evals_per_dim = [1000]
samplers = [EigenCovarSampler, CholeskyCovarSampler]
utils = [linear_utilities, log_utilities]
#utils = [linear_utilities]
#mus = [0.05, 0.10, 0.20, 0.30]
mus = [1, 0.05, 0.10, 0.20]
#mus = [1]
lambdas = [(n) -> int(4 + 3*floor(log(n))), (n) -> 4*n, (n) -> 4*n*n]
problems = [rosenbrock, sphere, cigar, cigtab]
#problems = [rosenbrock]

@everywhere machine1 = "macbook_pro_retina_cth"
@everywhere machine2 = "macbook_pro_ssd_bth"
@everywhere machine = machine1

# Now create parameter combinations
Params = Any[]
for(run in 1:num_runs_per_dim)
  for(prob in problems)
    for(d in dims)
      for(maxfevals in max_evals_per_dim)
        for(uf in utils)
          for(s in samplers)
            for(l in lambdas)

              if typeof(l) == Function
                l = l(d)
              end

              for(mu in mus)
                if typeof(mu) != Int64
                  mu = int(mu * l)
                end
                if mu < l
                  if mu < 1
                    mu = 1
                  end
                  push!(Params, (run, prob, d, maxfevals*d, uf, s, mu, l))
                end
              end

            end
          end
        end
      end
    end
  end
end

@everywhere function one_run(args)
  run, prob, d, maxfevals, uf, s, mu, l = args

  ts = strftime("%y%m%d,%H%M%S", time())

  tic()
  x, f, fevals = cmsa_es(d, prob;
    max_fevals = maxfevals, covarMatrixSampler = s,
    mu = mu, lambda = l,
    utilitiesFunc = uf, trace = false)
  exectime = toq()

  return "$(ts),$(run),$(prob),$(d),$(machine),$(s),$(uf),$(mu),$(l),$(fevals),$(exectime),$(f)"
end

@everywhere function one_run_and_save(args)
  res = one_run(args)
  println(STDERR, res)
  res
end

results = pmap(one_run_and_save, Params)

fh = open(strftime("parameter_studies/%y%m%d_%H%M%S_$(machine).csv", time()), "w")
println(fh, "Date,Time,RunID,Problem,N,Machine,Sampler,Utilities,Mu,Lambda,NumFuncEvals,ExecTime,Fitness")
for(i in 1:length(results))
  println(fh, results[i])
end
close(fh)


text:
@doc raw"""
	quantiledep(x,y,q)

Quantile Dependence at the `q`-th quantile of vector-like data `x` and `y`.

Let ``q_x`` and ``q_y`` be the `q`-th quantiles of `x` and `y`.

Returns ``E[xy | x \leq q_x, y \leq q_y] / q`` if `q<=0.5` and
``(1 -2q + E[xy | x \leq q_x, y \leq q_y]) / (1-q)`` otherwise.
"""
function quantiledep(x,y,q)
    qx = quantile(x,q)
    qy = quantile(y,q)
    tmp = mean((x.<=qx).*(y.<=qy))
    if q <= 0.5
        return tmp/q
    else
        return (1-2*q + tmp)/(1-q)
    end
end

"""
	exceedancecor(x,y,c;kind=:both,side=:lower)

Computes exceedance correlations between `x` and `y`, using exceedance threshold `c`.
Keyword `kind` refers to which type of conditioning: `:both` conditions on `x` and `y` simultaneously,
while `:onx` conditions only on `x`. Keyword `side` refers to look below (`:lower`) or above (`:upper`)
the threshold `c`. Inputs `x` and `y` are standardized before conditioning,
so `c` should be in standard deviation units.
"""
function exceedancecor(x,y,c;kind=:both,side=:lower)
    xs = standardize(x)
	ys = standardize(y)
    if side==:lower
        if kind==:both
            loc = (xs .<= c) .& (ys .<= c)
        elseif kind==:onx
            loc = (xs .<= c)
        end
    elseif side==:upper
        if kind==:both
            loc = (xs .>= c) .& (ys .>= c)
        elseif kind==:onx
            loc = (xs .>= c)
        end
    end
	return cor(xs[loc],ys[loc])
end

"""
	downsidebeta(x,y;kind=:both)

Computes the downside beta of `y` on `x` as defined in Ang & Chen.

`kind=:onx` allows to condition only on `x`.
"""
function downsidebeta(x,y;kind=:both)
	mux = mean(x)
	muy = mean(y)
	if kind==:both
		loc = (x .<= mux) .& (y .<= muy)
	elseif kind==:onx
		loc = (x .<= mux)
	end
	xtmp = x[loc]
	ytmp = y[loc]
	varx = var(xtmp)
	covxy = cov(xtmp,ytmp)
	return covxy/varx
end

"""
	upsidebeta(x,y;kind=:both)

Computes the upside beta of `y` on `x` as defined in Ang & Chen.

`kind=:onx` allows to condition only on `x`.
"""
function upsidebeta(x,y;kind=:both)
	mux = mean(x)
	muy = mean(y)
	if kind==:both
		loc = (x .>= mux) .& (y .>= muy)
	elseif kind==:onx
		loc = (x .>= mux)
	end
	xtmp = x[loc]
	ytmp = y[loc]
	varx = var(xtmp)
	covxy = cov(xtmp,ytmp)
	return covxy/varx
end

# x here is "the market return"

@doc raw"""
	coskewness(x,y)

Computes the coskewness between vectors `x` and `y`. The vector `x` receives the
powers in the expression.

``\frac{E[(x - E[x])^2 (y-E[y])]}{\sigma_x^2 \sigma_y}``

"""
function coskewness(x,y)
    dmx = x .- mean(x)
    dmy = y .- mean(y)
    top = mean( dmy .* (dmx.^2))
    bot = std(dmy) * var(dmx)
    return top/bot
end

@doc raw"""
	cokurtosis(x,y;kind::Symbol=:asymmetric)

Computes the cokurtosis between vectors `x` and `y`. The default is asymmetric cokurtosis.

`kind=:asymmetric`:
``\frac{E[(x - E[x])^3 (y-E[y])]}{\sigma_x^3 \sigma_y}``

`kind=:symmetric`:
``\frac{E[(x - E[x])^2 (y-E[y])^2]}{\sigma_x^2 \sigma_y^2}``
"""
function cokurtosis(x,y;kind::Symbol=:asymmetric)
    dmx = x .- mean(x)
    dmy = y .- mean(y)
    if kind == :symmetric
        top = mean( (dmy.^2) .* (dmx.^2))
        bot = var(dmy) * var(dmx)
    else
        top = mean( dmy .* (dmx.^3))
        bot = std(dmy)  * std(dmx)^3
    end
    return top/bot
end

#### TODO
# - better API

#### Correlation-based (Ang & Chen 2001 + statistical test of Hong, Tu, & Zhou 2007 + one-sided as in Schreindorfer 2019)
function downsidecor(x,y)
    qs = 0.15:0.025:1
    out = zeros(length(qs))
    for (i,qi) in enumerate(qs)
        qxi = quantile(x,qi)
        loc = x .<= qxi
        tmpxi = x[loc]
        tmpyi = y[loc]
        out[i] = cor(tmpxi,tmpyi)
        #println([length(tmpxi) length(tmpyi) length(tmpxj) length(tmpyj)])
    end
    return out
end

function upsidecor(x,y)
    qs = 0:0.025:0.85
    out = zeros(length(qs))
    for (i,qi) in enumerate(qs)
        qxi = quantile(x,qi)
        loc = x .>= qxi
        tmpxi = x[loc]
        tmpyi = y[loc]
        out[i] = cor(tmpxi,tmpyi)
        #println([length(tmpxi) length(tmpyi) length(tmpxj) length(tmpyj)])
    end
    return out
end

##### Asymmetry test from Guofu Zhou et al.
#standardize(x) = (x .- mean(x)) ./ std(x)
bartlett(z) = abs(z) < 1 ? (1-abs(z)) : zero(eltype(z))
# this only does 1 value of c
function buildEC(x,y,c)
	T = length(x)
	xs = standardize(x)
	ys = standardize(y)
	iL = (xs .< -c) #.& (ys .< -c)
	iU = (xs .> c) #.& (ys .> c)
	lower_xs = xs[iL]
	upper_xs = xs[iU]
	lower_ys = ys[iL]
	upper_ys = ys[iU]
	Tl = length(lower_xs)
	Tu = length(upper_xs)
	lower_cor = cor(lower_xs,lower_ys)
	upper_cor = cor(upper_xs,upper_ys)

	# standardize again
	xs[iL] = standardize(lower_xs)
	ys[iL] = standardize(lower_ys)
	xs[iU] = standardize(upper_xs)
	ys[iU] = standardize(upper_ys)

	EC = (((T-Tu)/T) .* (xs.*ys .- upper_cor) .* iU) .- (((T-Tl)/T) .* (xs.*ys .- lower_cor) .* iL)
	return (EC=EC,corL=lower_cor,corU=upper_cor)
end

# this wraps it together
# conditioning is standardized x .< -c[i] vs standardized x .> c[i]
# p represents some sort of bandwith? -- lag order?
function asymmetrictest(x,y,c,p)
	T = length(x)
	Nc = length(c)

	if Nc == 1
		outs=buildEC(x,y,c)
		ECmat=outs.EC
		uppercors = outs.corU
		lowercors = outs.corL
		S = 0.0
		for L=-p:p
			w = bartlett(L/p)
			g_l = 0.0
			for t=(abs(L)+1):T
				g_l += ECmat[t]*ECmat[t-abs(L)]
			end
			S += w*(1/T)*g_l
		end
	else
		ECmat = zeros(eltype(x),T,Nc)
		uppercors = zeros(eltype(x),Nc)
		lowercors = zeros(eltype(x),Nc)
		for i =1:Nc
			outs=buildEC(x,y,c[i])
			ECmat[:,i]=outs.EC
			uppercors[i] = outs.corU
			lowercors[i] = outs.corL
		end
		S = zeros(eltype(x),Nc,Nc)
		for L=-p:p
			w = bartlett(L/p)
			g_l = zeros(eltype(x),Nc,Nc)
			for ic=1:Nc,jc=1:Nc
				for t=(abs(L)+1):T
					g_l[ic,jc] += ECmat[t,ic]*ECmat[t-abs(L),jc]
				end
			end
			S += w*(1/T).*g_l
		end
	end

	J = T * (uppercors - lowercors)' * inv(S) * (uppercors - lowercors)
	pval = pdf(Chisq(Nc),J)
	return (pval=pval,J=J,S=S,corL=lowercors,corU=uppercors,T=T)
end

#### Entropy-based (Jiang, Wu, & Zhou 2018)
quadrantprobabilities(f,c;kwargs...) = (LQP=hcubature(f, [-10;-10], [-c;-c];kwargs...), UQP=hcubature(f, [c;c], [10;10];kwargs...))
function entropy_S(f,c;kwargs...)
    newf(x) = (sqrt(abs(f(x))) - sqrt(abs(f(-x))))^2
    int,se = hcubature(x->newf(x),[c;c],[20;20];kwargs...)
    return (int/2,se)
end
function DOWN_ASY_JWZ(f,c=0;kwargs...)
    lqp,uqp = quadrantprobabilities(f,c;kwargs...)
    S = entropy_S(f,c;kwargs...)
    val = sign(first(lqp)-first(uqp))*first(S)
    return val
end
# FF = kde((standardize(dfnew[!,:GDPg]),standardize(dfnew[!,:NDSg])))
# PF = InterpKDE(FF)
# quadrantprobabilities.(x->pdf(PF,x[1],x[2]),0:0.5:1.5)


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

using Markdown
using InteractiveUtils

# ╔═╡ 88af7dc6-70fd-11eb-3af2-f7a194340290
begin
	using PlutoUI, Logging
	# set global loglevel
	Logging.global_logger(ConsoleLogger(stderr, Logging.LogLevel(-5000)))
	Logging.disable_logging(Logging.LogLevel(-1001))
end

# ╔═╡ bb6ed424-70fc-11eb-2a2a-ff0556bb7cda
md"""
# Rational agents

## Our first agent
Consider the following example from the book: a simple vacuum-cleaner agent that cleans a square if it is dirty and moves to the other square if not. We consider the assume the following:
* Performance measure: we awards one point for each clean square at each time step, over a “lifetime” of 1000 time steps.
* Environment: the “geography” of the environment is known a priori but the dirt distribution and the initial location of the agent are not. Clean squares stay clean and sucking cleans the current square. The Left and Right actions move the agent left and right except when this would take the agent outside the environment, in which case the agent remains where it is.
* Actuators: the only available actions are Left, Right, and Suck.
* Sensors: the agent correctly perceives its location and whether that location contains dirt.

$(PlutoUI.LocalResource("./img/vacuum.png", :width => 500, :align => "middle"))

Questions:

1. How many world states are there?
2. Identify the goal states amongst those from q1.
3. Implement a performance-measuring environment simulator for the vacuum-cleaner world. Thinks about a modular implementation so that the sensors, actuators, and environment characteristics (size, shape, dirt placement, etc.) can be changed easily.
4. What is the value of the performance measure for each world state? Confirm your simulation returns the expected value.
"""

# ╔═╡ 9fd46218-70fe-11eb-1ddd-132f13e0611c
begin
	# build the world
	struct world
		locations::Dict{Tuple{Int64,Int64}, Bool}
	end

	function world(locations::Array{Tuple{Int64,Int64}})
		return world(Dict((location => rand(Bool) for location in locations)))
	end

	function world(locations::Array{Tuple{Int64,Int64}}, dirt::Array{Bool,1})
		if length(locations) == length(dirt)
			return world(Dict(zip(locations, dirt)))
		end
	end

	mutable struct reflex_vacuum
		w::world
		percept::Tuple{Tuple{Int64, Int64}, Bool}
	end

	# include the actions
	function update!(v::reflex_vacuum)
		pos = v.percept[1]
		status = v.percept[2]

		@debug "current pos: $(pos), status: $(status)"
		if status == true # if dirty => clean
			@debug "status is dirty, cleaning..."
			v.w.locations[pos] = false
			v.percept = (pos, false)
		else # if clean => move
			pos = rules[pos]
			v.percept = (pos, v.w.locations[pos]) # update sensor
		end
	end

	function runprog(v::reflex_vacuum, n::Int=1000)
		score = 0
		for _ in 1:n
			score += sum(.!(values(v.w.locations)))
			update!(v)
		end
		return score
	end

	# __main__

	# set the rules
	const rules = Dict((0,0)=>(1,0), (1,0)=>(0,0))

	@info "Starting..."
	for layout in Base.Iterators.product([true, false],[true, false])
		for startpos in [(0,0),(1,0)]
			# initiate the world
			w = world([(0,0),(1,0)],[i for i in layout])
			# initiate reflex_vacuum
			v = reflex_vacuum(w, (startpos, w.locations[startpos]))
			# scoring
			score = runprog(v)
			# output msg
			@info "layout: $(layout), startpos: $(startpos), score: $(score)"
		end
	end
	@info "Finished"
end

# ╔═╡ ab564aec-70ff-11eb-0d36-cd95fb8def23
md"""
## Modified agent
Consider a modified version of the vacuum environment, in which the agent is penalized one point for each movement.

1. Can a simple reflex agent be perfectly rational for this environment? Explain.
2. What about a reflex agent with state?
3. How do your answers to a and b change if the agent’s percepts give it the clean/dirty status of every square in the environment?
"""

# ╔═╡ f2915910-70ff-11eb-2f60-e9963db3aaeb
md"""
## Another modified agent
Consider a modified version of the simple vacuum environment in which the geography of the environment — its extent, boundaries, and obstacles — is unknown, as is the initial dirt configuration. (The agent can go Up and Down as well as Left and Right.)

1. Can a simple reflex agent be perfectly rational for this environment? Explain.
2. Can a simple reflex agent with a randomized agent function outperform a simple reflex agent? Design such an agent (principle) and measure its performance on several environments (if you have time).
3. Can you design an environment in which your randomized agent will perform poorly? Show your results.
"""

# ╔═╡ Cell order:
# ╠═88af7dc6-70fd-11eb-3af2-f7a194340290
# ╠═bb6ed424-70fc-11eb-2a2a-ff0556bb7cda
# ╠═9fd46218-70fe-11eb-1ddd-132f13e0611c
# ╟─ab564aec-70ff-11eb-0d36-cd95fb8def23
# ╟─f2915910-70ff-11eb-2f60-e9963db3aaeb

max_stars_repo_path: PS/PS01 - rational agents.jl

max_issues_repo_path: PS/PS01 - rational agents.jl

max_forks_repo_path: PS/PS01 - rational agents.jl


text:
include("bitstrings.jl")
include("graphs.jl")
include("sats.jl")


text:
abstract XMLError

immutable XMLParseError <: XMLError
    msg::String
end

immutable XMLNoRootError <: XMLError
end

immutable XMLAttributeNotFound <: XMLError
end

immutable XMLWriteError <: XMLError
    msg::String
end

immutable XMLTreeError <: XMLError
    msg::String
end

const dom_exception_causes = [
    "Index size error",        #  1
    "DOM string size error",   #  2
    "Hierarchy request error", #  3
    "Wrong document",          #  4
    "Invalid character",       #  5
    "No data allowed",         #  6
    "No modification allowed", #  7
    "Not found",               #  8
    "Not supported",           #  9
    "Inused attribute"         # 10
]

immutable DOMException <: XMLError
    code::Int
    cause::String

    DOMException(code::Int) = new(code, dom_exception_causes[code])
end


text:
using Test

using Revise
push!(LOAD_PATH, ".")
using Haptics

Haptics.get_version_info()
@test Haptics.compatible()

Haptics.get_error_msg.(-100:100)
Haptics.get_error_msg.(-1:28)

Haptics.last_error_msg()

Haptics.tick()

Haptics.count_devices()

# test existing
for dev in -1:(Haptics.count_devices()-2)
    Haptics.get_device_name(dev)
end

# test missing
Haptics.get_device_name(10)

# test valid (mock device)
Haptics.initialize(-1)
# test it twice: error already init
Haptics.initialize(-1)

Haptics.deinitialize()
# Try it twice
Haptics.deinitialize()

# Start without init
Haptics.start()

# Start after init ok
Haptics.initialize(-1)
Haptics.start()
Haptics.deinitialize()
# Deinit twice fail again
Haptics.deinitialize()

Haptics.initialize(-1)
Haptics.start()
# start twice fails
Haptics.start()
Haptics.deinitialize()

# Initialize when running
Haptics.initialize(-1)
Haptics.start()
Haptics.initialize(-1)
Haptics.stop()
Haptics.initialize(-1)
Haptics.start()
Haptics.deinitialize() # found a crash

# No loop frequency if thread is not initialized
@test Haptics.get_loop_frequency() < 0

# No loop frequency if thread is not running
Haptics.initialize(-1)
@test Haptics.get_loop_frequency() < 0
Haptics.deinitialize()

# Loop frequency if thread is running
Haptics.initialize(-1)
Haptics.start()
# let the thread start
sleep(1.0)
@test Haptics.get_loop_frequency() > 500
Haptics.deinitialize()

# Loop number -1 if not inited
@test Haptics.get_loops() == -1

# Loop number 0 if not running
Haptics.initialize(-1)
# @test Haptics.get_loops() == 0

# NOT WORKING!
# # Loop number strictly growing
# n = Haptics.get_loops()
# @show n
# Haptics.initialize(-1)
# Haptics.start()
# for i in 1:5
#     sleep(1.0)
#     n2 = Haptics.get_loops()
#     @test n2 > n
#     n = n2
# end
# Haptics.stop()
# Haptics.deinitialize()

Haptics.get_log_frame()
Haptics.deinitialize()

Haptics.initialize(-1)
Haptics.start()
sleep(1)
Haptics.get_log_frame()
Haptics.stop()
Haptics.deinitialize()

Haptics.is_initialized()
Haptics.is_running()

# Mock device, can't be pressed
@test all(Haptics.is_tool_button_pressed.(-100:100)) == false

Haptics.initialize(-1)
Haptics.start()
Haptics.is_initialized()
Haptics.is_running()

# Not logging
Haptics.get_log_frame()

Haptics.get_loops()

Haptics.stop()
Haptics.deinitialize()

## those tests require a real device!
@test Haptics.count_devices() == 2

@test Haptics.initialize(0) == 0
@test Haptics.start() == 0
@test Haptics.stop() == 0
@test Haptics.deinitialize() == 0

#=
init
- count
- get device <- fails
=#

for i in 1:10
    @show i
    @test Haptics.count_devices() == 2
    @test Haptics.initialize(0) == 0
    @test Haptics.start() == 0
    @test Haptics.stop() == 0
    @test Haptics.deinitialize() == 0
    @info "Before: $(Haptics.count_devices())"
    @test Haptics.initialize(-1) == 0 ## this is required. Why?
    sleep(1)
    @test Haptics.deinitialize() == 0
    @info "After: $(Haptics.count_devices())"
end

for i in 1:100
    @test Haptics.initialize(-1) == 0
    @test Haptics.deinitialize() == 0
    @test Haptics.count_devices() == 2
end

for i in 1:100
    @show i
    @test Haptics.initialize(0) == 0
    @test Haptics.deinitialize() == 0
    @test Haptics.count_devices() == 2
end

max_stars_repo_name: HapticPlugin/HPGE-wrappers

max_issues_repo_name: HapticPlugin/HPGE-wrappers

max_forks_repo_name: HapticPlugin/HPGE-wrappers


text:
module general
using Unitful

"""
    P_Antoine(CA, T, i)

```math
\\log_{10}{P(mmHg)}=A-\\frac{B}{C+T(°C)}
```

**Los campos de entrada son:**

- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `T :: Float` Es la temperatura del sistema en °C
- `i :: Int` Componente `i` al que se quiere conocer su presión

**Salida:**

Regresa la presión `P` en mmHg

- `P :: m-element Array{Float64,1}` presión
"""
function P_Antoine(CA::Array, T::Unitful.Temperature, i::Int)::Unitful.Pressure
    A,B,C = CA[i, :]
    T = T |> u"K"
    P = u"Torr"*10^(A-B/((T/1u"K" -273.15)+C))
    return P
end

"""
    T_Antoine(CA, P, i)

```math
T(°C)=-C-\\frac{B}{log_{10}{P(mmHg)}-A}
```

**Los campos de entrada son:**

- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `P :: Float` Es la presión del sistema en mmHg
- `i :: Int` Componente `i` al que se quiere conocer su temperatura

**Salida:**

Regresa la temperatura `T` en °C

- `T :: m-element Array{Float64,1}` temperatura
"""
function T_Antoine(CA::Array, P::Unitful.Pressure, i::Int)::Unitful.Temperature
    A,B,C = CA[i, :]
    P = P |> u"Torr"
    (-B/(log10(P/1u"Torr")-A) -C +273.15)*u"K"
end

function P_Antoine(CA::Array, T, i::Int)::Float64
    A,B,C = CA[i, :]
    P = 10^(A - B / ((T - 273.15) + C))
    return P
end

function T_Antoine(CA::Array, P, i::Int)::Float64
    A,B,C = CA[i, :]
    (-B / (log10(P) - A) - C + 273.15)
end

function nelemen(n::Int, ñ::Int)::Int
    m = Array{Int}(undef,ñ,n)
    for i=1:n, j=1:ñ
        if i == 1 || j == 1
            m[j, i] = 1
        else
            m[j, i] = m[j, i - 1] + m[j - 1, i]
        end
    end
    return m[ñ, n]
end
function nc(xx::StepRangeLen...)
    kk = collect(Iterators.product(xx...))
    sk = sum.(kk)
    m = length(xx)
    ñ = length(xx[1])
    n = length(sk)
    nn = nelemen(m + 1,ñ)
    x = Array{Float64}(undef,m + 1,nn)
    k = 1
    for i=1:n
        if sk[i] <= 1
            for j=1:m
                x[j, k] = kk[i][j]
            end
            x[m + 1, k] = 1 - sk[i]
            k +=1
        end
    end
    return x
end
function nc(xx::Array...)
    m = length(xx)
    ñ = length(xx[1])
    x = Array{Float64}(undef,m + 1,ñ)
    for i=1:ñ
        sk = 0
        for j=1:m
            x[j, i] = xx[j][i]
            sk += xx[j][i]
        end
        if sk <= 1
            x[m + 1, i] = 1 - sk
        else
            println("La suma de las primeras variables es mayor a 1")
        end
    end
    return x
end
function nc_ord(xx::StepRangeLen...)
    kk = collect(Iterators.product(xx...))
    sk = sum.(kk)
    m = length(xx)
    ñ = length(xx[1])
    n = length(sk)
    nn = nelemen(m + 1,ñ)
    x = Array{Float64}(undef,m + 1,nn)
    k = 1
    iter = 1
    for i=1:ñ
        if iter == 1
            for h=1:ñ
                if sk[i, h] <= 1
                    for j=1:m
                        x[j, k] = kk[i, h][j]
                    end
                    x[m + 1, k] = 1 - sk[i, h]
                    k +=1
                end
            end
            iter = 0
        else
            for h=reverse(1:ñ)
                if sk[i, h] <= 1
                    for j=1:m
                        x[j, k] = kk[i, h][j]
                    end
                    x[m + 1, k] = 1 - sk[i, h]
                    k +=1
                end
            end
            iter = 1
        end
    end
    return x
end

end  # module

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
include("src/init.jl")

"""
This file contains the script to execute learning containing
the Value and Critic function in one Network!
"""

θ = init_θ_combined(params.structure_net_pl)
res = zeros(params.n_iter)
precision = zeros(params.n_iter)
all_traj  = init_all_traj(params)

function actor_critic_combined(params::Params, θ::Any, res::Array{Float64,1}, precision::Array{Float64,1}, all_traj::Traj_Data)
    #Init Part:
    #------------------------
    #Initialize the Network
    ĝ = copy(θ)
    opt_para = optimizers(θ,   Adam; lr = params.lr_pl)

    #Init Progress bar
    progress = Progress(params.n_iter,1)

    # Preallocate trajectory
    rms_state = init_rms(params.s_init)

    a_t_old = all_traj.traj[1].a_t
    weights_old = zeros(length(net_to_vec(θ)),2)
    #FIM = Symmetric(eye(length(net_to_vec(θ))))*0 #preallocate FisherMatrix

    for i_iter = 1:params.n_iter
        ĝ *= 0

        #Make a Video
        if floor(i_iter/400) == ceil(i_iter/400)
            figure(1)
            #t_horizon = (all_traj.traj[1].len)*(params.dt)
            t_plot = 0:params.dt:params.t_horizon;
            x1 = all_traj.traj[1].s_t
            u1 = all_traj.traj[1].a_t'
            filename = string("Results/Videos/Actor_Critic_",Dates.today(),"_iter_",i_iter,".mp4")
            visualise(:SegwayRobot, x1', t_plot, filename; options = Dict(:U => u1, :size => (1920, 1080)))
        end

        #-----------------------------------------------------------------------
        #Solve it in Parallel! (How to do asynchron updates..?)
        ĝ, res[i_iter], precision[i_iter] = parallel_episodes_combined(params,θ,all_traj,rms_state)

        #Reorganice Data in minibatches
        mini_batch = build_minibatches(all_traj,params)

        #Do the Policy / and Value updates
        θ = parallel_PPO_combined(params, θ, mini_batch, opt_para)

        #Visualization_stuff
        #-----------------------------------------------------------------------
        ProgressMeter.update!(progress,i_iter) #update progress bar
            a_t_old ,weights_old = plot_stuff(params, all_traj, a_t_old,
            weights_old, i_iter,θ, θ, res, precision/params.n_traj, ĝ, ĝ)
        #-----------------------------------------------------------------------
    end
end

actor_critic_combined(params, θ, res, precision, all_traj)
#save_data(res, precision, θa, θv, all_traj)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using MP3

if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

# encoding and decoding can introduce some delay
DELAY_THRESHOLD = 1160

# write your own tests here
@testset "Loading MP3" begin
    reference = load(joinpath(dirname(@__FILE__), "Sour_Tennessee_Red_Sting.mp3"))
    @test typeof(reference) == SampledSignals.SampleBuf{MP3.PCM16Sample,2,MP3.Hertz}
    @test size(reference, 2) == 2
    @test abs(size(reference, 1) - 245376) <= DELAY_THRESHOLD
    @test reference.samplerate == 44100Hz
    @test reference.samplerate == 44.1kHz
end

@testset "Saving MP3" begin
    reference = load(joinpath(dirname(@__FILE__), "Sour_Tennessee_Red_Sting.mp3"))

    outpath = "$(tempname()).mp3"
    save(outpath, reference)
    audio = load(outpath)
    @test size(audio) == size(reference)
    @test audio.samplerate == reference.samplerate

    for samplerate in [8000, 11025, 12000, 16000, 22050, 24000, 32000, 44100, 48000]
        save(outpath, reference; samplerate = samplerate)
        audio = load(outpath)
        @test audio.samplerate == samplerate * Hz
    end

    for bitrate in [96, 128, 160, 192, 224, 256, 320]
        save(outpath, reference; bitrate = bitrate)
        audio = load(outpath)
        @test size(audio, 2) == size(reference, 2)
        @test abs(size(audio, 1) - size(reference, 1)) <= DELAY_THRESHOLD
    end

    save(outpath, reference[:, 1])
    audio = load(outpath)
    @test size(audio, 2) == 1

    save(outpath, reference[:, 2])
    audio = load(outpath)
    @test size(audio, 2) == 1

    save(outpath, reference[:, 1]; nchannels = 2)
    audio = load(outpath)
    @test size(audio, 2) == 2

    f32 = map(Float32, reference)
    save(outpath, f32)
    audio = load(outpath)
    @test audio.samplerate == 44100Hz
    @test size(audio, 2) == size(reference, 2)
    @test abs(size(audio, 1) - size(reference, 1)) <= DELAY_THRESHOLD

    f64 = map(Float32, reference)
    save(outpath, f64)
    @test audio.samplerate == 44100Hz
    @test size(audio, 2) == size(reference, 2)
    @test abs(size(audio, 1) - size(reference, 1)) <= DELAY_THRESHOLD

end


text:
module TestExamples
include("preamble.jl")

using Bifurcations: resolved_points
using Bifurcations.Codim1: resolve_point
using Bifurcations.Codim2: BackReferencingAS, NormalizingAS
using Bifurcations.Continuations: find_errors, print_errors
using Bifurcations: examples
using Bifurcations.Examples.Calcium
using Bifurcations.Examples.Calcium: CalciumParam
using Bifurcations.Examples.Bazykin85: Bazykin85Param
using Bifurcations.Examples.BazykinKhibnik81: BazykinKhibnik81Param
using Bifurcations.Examples.DuffingVanDerPol: DuffingVanDerPolParam
using Bifurcations.Examples.Bautin: BautinParam
using Bifurcations.Examples.Cusp: CuspParam

struct Example{P}
    prob::P
end

make_codim2(::Any, ::Any, ::Any) = []
make_codim2(::Real, point, solver1) = []
make_codim2(::Tuple, point, solver1) = []  # PredatorPrey
make_codim2(::DuffingVanDerPolParam, point, solver1) = []
make_codim2(::CuspParam, point, solver1) = []  # TODO: define something

make_codim2(::CalciumParam, point, solver1) = [
    BifurcationProblem(
        point,
        solver1,
        (@lens _.gca),
        (0.0, 8.0),
    )
]

function make_codim2(::Bazykin85Param, point, solver1)
    if (point.point_type == Codim1.PointTypes.saddle_node ||
        isapprox(point.u[2], 0; atol=1e-3))
        return []
    end
    return [BifurcationProblem(
        point,
        solver1,
        (@lens _.δ),
        (0.0, 10.0),
    )]
end

function make_codim2(::BazykinKhibnik81Param, point, solver1)
    if point.point_type !== Codim1.PointTypes.hopf
        return []
    end
    resolved = resolve_point(point, solver1)
    if any(@. resolved.u[1:2] < 1e-3)
        return []
    end
    return [BifurcationProblem(
        resolved,
        solver1,
        (@lens _.n),
        (-1e-3, 2.0),
    )]
end

function make_codim2(::BautinParam, point, solver1)
    args = (
        point,
        solver1,
        (@lens _.β₂),
        (-2.0, 2.0),
    )
    return [
        BifurcationProblem(
            args...;
            augmented_system = BackReferencingAS(),
        ),
        BifurcationProblem(
            args...;
            augmented_system = NormalizingAS(),
        ),
    ]
end

EXAMPLES = [
    examples()...,
    :CalciumSTDArray => Example(Calcium.make_prob(
        f = generalized_f(Calcium),
        u0 = Array(Calcium.u0),
    )),
]

@testset "example $name" for (name, ex) in EXAMPLES
    prob1 = ex.prob
    solver1 = init(prob1)
    solve!(solver1)
    errors = find_errors(solver1)
    print_errors(errors)
    @test isempty(errors)

    @testset "$i-th $(point.point_type); $j-th problem" for
            (i, point) in enumerate(resolved_points(solver1)),
            (j, prob2) in enumerate(
                make_codim2(prob1.p.de_prob.p, point, solver1))
        solver2 = init(prob2)
        solve!(solver2)
        errors = find_errors(solver2)
        print_errors(errors)
        @test isempty(errors)
    end
end

end  # module


text:
"""
`SqExponentialKernel([ρ=1.0])`

The squared exponential kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-ρ²‖x-y‖²)
```
See also [`ExponentialKernel`](@ref) for a
related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
"""
struct SqExponentialKernel{Tr} <: Kernel{Tr}
    transform::Tr
end

@inline kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
@inline iskroncompatible(::SqExponentialKernel) = true

metric(::SqExponentialKernel) = SqEuclidean()

## Aliases ##
const RBFKernel = SqExponentialKernel
const GaussianKernel = SqExponentialKernel

"""
`ExponentialKernel([ρ=1.0])`
The exponential kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-ρ‖x-y‖)
```
"""
struct ExponentialKernel{Tr} <: Kernel{Tr}
    transform::Tr
end

@inline kappa(κ::ExponentialKernel, d::Real) = exp(-d)
@inline iskroncompatible(::ExponentialKernel) = true
metric(::ExponentialKernel) = Euclidean()

## Alias ##
const LaplacianKernel = ExponentialKernel

"""
`GammaExponentialKernel([ρ=1.0, [γ=2.0]])`
The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-ρ^(2γ)‖x-y‖^(2γ))
```
"""
struct GammaExponentialKernel{Tr, Tγ<:Real} <: Kernel{Tr}
    transform::Tr
    γ::Tγ
    function GammaExponentialKernel{Tr,Tγ}(t::Tr, γ::Tγ) where {Tr<:Transform,Tγ<:Real}
        @check_args(GammaExponentialKernel, γ, γ >= zero(Tγ), "γ > 0")
        return new{Tr, Tγ}(t, γ)
    end
end

params(k::GammaExponentialKernel) = (params(transform),γ)
opt_params(k::GammaExponentialKernel) = (opt_params(transform),γ)

function GammaExponentialKernel(ρ::Real=1.0, γ::Real=2.0)
    GammaExponentialKernel(ScaleTransform(ρ), γ)
end

function GammaExponentialKernel(ρ::AbstractVector{<:Real}, γ::Real=2.0)
    GammaExponentialKernel(ARDTransform(ρ), γ)
end

function GammaExponentialKernel(t::Tr, γ::Tγ=2.0) where {Tr<:Transform, Tγ<:Real}
    GammaExponentialKernel{Tr, Tγ}(t, γ)
end

@inline kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^κ.γ)
@inline iskroncompatible(::GammaExponentialKernel) = true
metric(::GammaExponentialKernel) = SqEuclidean()


text:
"""
We demonstrate core capabilities of EquivariantOperators.jl: particle mesh, finite differences, and machine learning. In a particle mesh context, we place 2 +1 point charges in a grid creating a scalar field. We then do a finite difference calculation of the electric field and potential. Switching gears to machine learning, we train our equivariant neural network to learn this transformation from the charge distribution (scalar field) to the electric potential (scalar field) and field (vector field), essentially learning the Green's function solution to Poisson's Equation. Finally, we treat our point charges as proton nuclei and train another neural network to predict the ground state electronic density of H2.
"""

repo = ".."
include("$repo/src/neural_networks.jl")
# include("../src/plotutils.jl")
# include("$repo/src/EquivariantOperators.jl")
# using .EquivariantOperators

using FileIO
using Plots
using Random
Random.seed!(1)

name = "h2"
case = load("..\\density-prediction\\data\\$name.jld2", "cases")[1]
@show positions = case.positions'
@show charges = case.charges
natoms = length(charges)
@show r = abs(positions[1, 1] - positions[1, 2]) / 2
@show center_charge = vec(sum(positions, dims = 2) / 2)

# grid params
@show dx = case.resolution
@show origin = case.origin
electronic_density = case.density
sz = size(electronic_density)
grid = Grid(dx, sz; origin)
@unpack sz, dV = grid

chargesum = sum(electronic_density) * dV
@assert chargesum ≈ sum(charges)

# place dipole point charges 2r apart in grid as input
rank = 0
dims = 3
ρp = Field(; grid, rank)
put_point_source!(
    ρp,
    grid,
    positions,
    reshape(charges, (1, length(charges))),
)
# plot(ρp(1))

# tensor convolution params

# generate data with Green's functions for Poisson's Equation, Gauss's law
rmax = 3.0
rank_max = 1
in_ = Props((1,), grid) # input scalar field
x = cat(ρp, dims = 4)
y = cat(electronic_density, dims = 4)
# define tensor field convolution layer
L = EquivConv(Props((1,), grid), Props((1,), grid), rmax)
Q = EquivProd(L.out; spectra = true)
D1 = EquivDense(Q.out, Props((1,), grid))

function normalize_density(density, total)
    # return density
    s = sum(density)
    if s == 0
        return density
    end
    density * total / s
end

function calcforces(rvecs,charges,ρ,grid;rmax=grid.dx*norm(grid.sz))
    name = :inverse_squared_field
    rmin=2grid.dx
    efield_op = LinearOperator(name; dx,rmin, rmax)
    efield =efield_op(ρ)
    vcat([charge*get(efield,grid,rvec) for (charge,rvec) in zip(charges,eachcol(rvecs))]...)
end
a=rand(length(Q.out.ranks))
ps = Flux.params(L, a)
# ps = Flux.params(L, D1)
f0(x)=abs.(sum([a[i]*Q(L(x))[:,:,:,i:i] for i=eachindex(a)]))
# f0(x) =sum( abs.((Q ∘ L)(x)))
# f0(x) = abs.((D1 ∘ Q ∘ L)(x))
# f0(x)=abs.(L(x))
# f0(x)=normalize_density(abs.((D1∘Q∘L)(x)),chargesum,dV)
# f(x) = normalize_density(f0(x) - f0(zeros(sz..., 1)),chargesum,dV)
f(x) = normalize_density(f0(x) - f0(zeros(sz..., 1)),chargesum/dV)

sumy=sum(y)
function loss()
    yhat = f(x)
    @show l = nae(yhat,y;sumy)
    l
end
loss()

data = [()]
opt = ADAM(0.1)

for i = 1:20
    Flux.train!(loss, ps, data, opt)
end
ρ=electronic_density+ρp
@show forces=calcforces(positions,charges,ρ,grid)
electronic_density_hat=yhat = f(x)
ρhat=electronic_density_hat+ρp
@show forces_hat=calcforces(positions,charges,ρ_hat,grid)

plot(x[:, :, :, 1])
plot(yhat[:, :, :, 1])
plot(y[:, :, :, 1])


text:
# hand made using grep
# these could probably be weeded out more

const  ADAMS_Q_MAX                =  12
const  BDF_Q_MAX                  =  5
const  CVDIAG_ILL_INPUT           =  -3
const  CVDIAG_INV_FAIL            =  -5
const  CVDIAG_LMEM_NULL           =  -2
const  CVDIAG_MEM_FAIL            =  -4
const  CVDIAG_MEM_NULL            =  -1
const  CVDIAG_NO_ADJ              =  -101
const  CVDIAG_RHSFUNC_RECVR       =  -7
const  CVDIAG_RHSFUNC_UNRECVR     =  -6
const  CVDIAG_SUCCESS             =  0
const  CVDLS_ILL_INPUT            =  -3
const  CVDLS_JACFUNC_RECVR        =  -6
const  CVDLS_JACFUNC_UNRECVR      =  -5
const  CVDLS_LMEMB_NULL           =  -102
const  CVDLS_LMEM_NULL            =  -2
const  CVDLS_MEM_FAIL             =  -4
const  CVDLS_MEM_NULL             =  -1
const  CVDLS_NO_ADJ               =  -101
const  CVDLS_SUCCESS              =  0
const  CVSPILS_DGMAX              =  0.2
const  CVSPILS_EPLIN              =  0.05
const  CVSPILS_ILL_INPUT          =  -3
const  CVSPILS_LMEMB_NULL         =  -102
const  CVSPILS_LMEM_NULL          =  -2
const  CVSPILS_MAXL               =  5
const  CVSPILS_MEM_FAIL           =  -4
const  CVSPILS_MEM_NULL           =  -1
const  CVSPILS_MSBPRE             =  50
const  CVSPILS_NO_ADJ             =  -101
const  CVSPILS_PMEM_NULL          =  -5
const  CVSPILS_SUCCESS            =  0
const  CV_ADAMS                   =  1
const  CV_BAD_DKY                 =  -26
const  CV_BAD_IS                  =  -45
const  CV_BAD_K                   =  -24
const  CV_BAD_T                   =  -25
const  CV_BAD_TB0                 =  -104
const  CV_BDF                     =  2
const  CV_CENTERED                =  1
const  CV_CONV_FAILURE            =  -4
const  CV_ERR_FAILURE             =  -3
const  CV_FAIL_BAD_J              =  1
const  CV_FAIL_OTHER              =  2
const  CV_FIRST_QRHSFUNC_ERR      =  -32
const  CV_FIRST_QSRHSFUNC_ERR     =  -52
const  CV_FIRST_RHSFUNC_ERR       =  -9
const  CV_FIRST_SRHSFUNC_ERR      =  -42
const  CV_FORWARD                 =  2
const  CV_FUNCTIONAL              =  1
const  CV_FWD_FAIL                =  -106
const  CV_GETY_BADT               =  -107
const  CV_HERMITE                 =  1
const  CV_ILL_INPUT               =  -22
const  CV_LINIT_FAIL              =  -5
const  CV_LSETUP_FAIL             =  -6
const  CV_LSOLVE_FAIL             =  -7
const  CV_MEM_FAIL                =  -20
const  CV_MEM_NULL                =  -21
const  CV_NEWTON                  =  2
const  CV_NORMAL                  =  1
const  CV_NO_ADJ                  =  -101
const  CV_NO_BCK                  =  -103
const  CV_NO_FAILURES             =  0
const  CV_NO_FWD                  =  -102
const  CV_NO_MALLOC               =  -23
const  CV_NO_QUAD                 =  -30
const  CV_NO_QUADSENS             =  -50
const  CV_NO_SENS                 =  -40
const  CV_ONE_STEP                =  2
const  CV_POLYNOMIAL              =  2
const  CV_QRHSFUNC_FAIL           =  -31
const  CV_QSRHSFUNC_FAIL          =  -51
const  CV_REIFWD_FAIL             =  -105
const  CV_REPTD_QRHSFUNC_ERR      =  -33
const  CV_REPTD_QSRHSFUNC_ERR     =  -53
const  CV_REPTD_RHSFUNC_ERR       =  -10
const  CV_REPTD_SRHSFUNC_ERR      =  -43
const  CV_RHSFUNC_FAIL            =  -8
const  CV_ROOT_RETURN             =  2
const  CV_RTFUNC_FAIL             =  -12
const  CV_SIMULTANEOUS            =  1
const  CV_SRHSFUNC_FAIL           =  -41
const  CV_STAGGERED               =  2
const  CV_STAGGERED1              =  3
const  CV_SUCCESS                 =  0
const  CV_TOO_CLOSE               =  -27
const  CV_TOO_MUCH_ACC            =  -2
const  CV_TOO_MUCH_WORK           =  -1
const  CV_TSTOP_RETURN            =  1
const  CV_UNREC_QRHSFUNC_ERR      =  -34
const  CV_UNREC_QSRHSFUNC_ERR     =  -54
const  CV_UNREC_RHSFUNC_ERR       =  -11
const  CV_UNREC_SRHSFUNC_ERR      =  -44
const  CV_WARNING                 =  99
const  FCMIX_CVODE                =  1
const  FCMIX_IDA                  =  2
const  FCMIX_KINSOL               =  3
const  HMAX_INV_DEFAULT           =  0.0
const  HMIN_DEFAULT               =  0.0
const  IDADLS_ILL_INPUT           =  -3
const  IDADLS_JACFUNC_RECVR       =  -6
const  IDADLS_JACFUNC_UNRECVR     =  -5
const  IDADLS_LMEMB_NULL          =  -102
const  IDADLS_LMEM_NULL           =  -2
const  IDADLS_MEM_FAIL            =  -4
const  IDADLS_MEM_NULL            =  -1
const  IDADLS_NO_ADJ              =  -101
const  IDADLS_SUCCESS             =  0
const  IDASPILS_ILL_INPUT         =  -3
const  IDASPILS_LMEMB_NULL        =  -102
const  IDASPILS_LMEM_NULL         =  -2
const  IDASPILS_MEM_FAIL          =  -4
const  IDASPILS_MEM_NULL          =  -1
const  IDASPILS_NO_ADJ            =  -101
const  IDASPILS_PMEM_NULL         =  -5
const  IDASPILS_SUCCESS           =  0
const  IDA_BAD_DKY                =  -27
const  IDA_BAD_EWT                =  -24
const  IDA_BAD_IS                 =  -43
const  IDA_BAD_K                  =  -25
const  IDA_BAD_T                  =  -26
const  IDA_BAD_TB0                =  -104
const  IDA_CENTERED               =  1
const  IDA_CONSTR_FAIL            =  -11
const  IDA_CONV_FAIL              =  -4
const  IDA_EE                     =  4
const  IDA_ERR_FAIL               =  -3
const  IDA_FIRST_QRHS_ERR         =  -32
const  IDA_FIRST_QSRHS_ERR        =  -52
const  IDA_FIRST_RES_FAIL         =  -12
const  IDA_FORWARD                =  2
const  IDA_FWD_FAIL               =  -106
const  IDA_GETY_BADT              =  -107
const  IDA_HERMITE                =  1
const  IDA_ILL_INPUT              =  -22
const  IDA_LINESEARCH_FAIL        =  -13
const  IDA_LINIT_FAIL             =  -5
const  IDA_LSETUP_FAIL            =  -6
const  IDA_LSOLVE_FAIL            =  -7
const  IDA_MEM_FAIL               =  -21
const  IDA_MEM_NULL               =  -20
const  IDA_NN                     =  0
const  IDA_NORMAL                 =  1
const  IDA_NO_ADJ                 =  -101
const  IDA_NO_BCK                 =  -103
const  IDA_NO_FWD                 =  -102
const  IDA_NO_MALLOC              =  -23
const  IDA_NO_QUAD                =  -30
const  IDA_NO_QUADSENS            =  -50
const  IDA_NO_RECOVERY            =  -14
const  IDA_NO_SENS                =  -40
const  IDA_ONE_STEP               =  2
const  IDA_POLYNOMIAL             =  2
const  IDA_QRHS_FAIL              =  -31
const  IDA_QSRHS_FAIL             =  -51
const  IDA_REIFWD_FAIL            =  -105
const  IDA_REP_QRHS_ERR           =  -33
const  IDA_REP_QSRHS_ERR          =  -53
const  IDA_REP_RES_ERR            =  -9
const  IDA_REP_SRES_ERR           =  -42
const  IDA_RES_FAIL               =  -8
const  IDA_ROOT_RETURN            =  2
const  IDA_RTFUNC_FAIL            =  -10
const  IDA_SIMULTANEOUS           =  1
const  IDA_SRES_FAIL              =  -41
const  IDA_SS                     =  1
const  IDA_STAGGERED              =  2
const  IDA_SUCCESS                =  0
const  IDA_SV                     =  2
const  IDA_TOO_MUCH_ACC           =  -2
const  IDA_TOO_MUCH_WORK          =  -1
const  IDA_TSTOP_RETURN           =  1
const  IDA_WARNING                =  99
const  IDA_WF                     =  3
const  IDA_YA_YDP_INIT            =  1
const  IDA_Y_INIT                 =  2
const  KINBBDPRE_FUNC_UNRECVR     =  -12
const  KINBBDPRE_PDATA_NULL       =  -11
const  KINBBDPRE_SUCCESS          =  0
const  KINDLS_ILL_INPUT           =  -3
const  KINDLS_JACFUNC_RECVR       =  -6
const  KINDLS_JACFUNC_UNRECVR     =  -5
const  KINDLS_LMEM_NULL           =  -2
const  KINDLS_MEM_FAIL            =  -4
const  KINDLS_MEM_NULL            =  -1
const  KINDLS_SUCCESS             =  0
const  KINSPILS_ILL_INPUT         =  -3
const  KINSPILS_LMEM_NULL         =  -2
const  KINSPILS_MAXL              =  10
const  KINSPILS_MEM_FAIL          =  -4
const  KINSPILS_MEM_NULL          =  -1
const  KINSPILS_PMEM_NULL         =  -5
const  KINSPILS_SUCCESS           =  0
const  KIN_ETACHOICE1             =  1
const  KIN_ETACHOICE2             =  2
const  KIN_ETACONSTANT            =  3
const  KIN_FIRST_SYSFUNC_ERR      =  -14
const  KIN_ILL_INPUT              =  -2
const  KIN_INITIAL_GUESS_OK       =  1
const  KIN_LINESEARCH             =  1
const  KIN_LINESEARCH_BCFAIL      =  -8
const  KIN_LINESEARCH_NONCONV     =  -5
const  KIN_LINIT_FAIL             =  -10
const  KIN_LINSOLV_NO_RECOVERY    =  -9
const  KIN_LSETUP_FAIL            =  -11
const  KIN_LSOLVE_FAIL            =  -12
const  KIN_MAXITER_REACHED        =  -6
const  KIN_MEM_FAIL               =  -4
const  KIN_MEM_NULL               =  -1
const  KIN_MXNEWT_5X_EXCEEDED     =  -7
const  KIN_NONE                   =  0
const  KIN_NO_MALLOC              =  -3
const  KIN_REPTD_SYSFUNC_ERR      =  -15
const  KIN_STEP_LT_STPTOL         =  2
const  KIN_SUCCESS                =  0
const  KIN_SYSFUNC_FAIL           =  -13
const  KIN_WARNING                =  99
const  MSBSET_DEFAULT             =  10
const  MSBSET_SUB_DEFAULT         =  5
const  MXHNIL_DEFAULT             =  10
const  MXITER_DEFAULT             =  200
const  MXNBCF_DEFAULT             =  10
const  MXORDP1                    =  6
const  MXSTEP_DEFAULT             =  500
const  NUM_TESTS                  =  5
const  OMEGA_MAX                  =  0.9
const  OMEGA_MIN                  =  0.00001
const  PRINTFL_DEFAULT            =  0
const  Q_MAX                      =  ADAMS_Q_MAX
const  SPBCG_ATIMES_FAIL_REC      =  4
const  SPBCG_ATIMES_FAIL_UNREC    =  -2
const  SPBCG_CONV_FAIL            =  2
const  SPBCG_MEM_NULL             =  -1
const  SPBCG_PSET_FAIL_REC        =  5
const  SPBCG_PSET_FAIL_UNREC      =  -4
const  SPBCG_PSOLVE_FAIL_REC      =  3
const  SPBCG_PSOLVE_FAIL_UNREC    =  -3
const  SPBCG_RES_REDUCED          =  1
const  SPBCG_SUCCESS              =  0
const  SPGMR_ATIMES_FAIL_REC      =  5
const  SPGMR_ATIMES_FAIL_UNREC    =  -2
const  SPGMR_CONV_FAIL            =  2
const  SPGMR_GS_FAIL              =  -4
const  SPGMR_MEM_NULL             =  -1
const  SPGMR_PSET_FAIL_REC        =  6
const  SPGMR_PSET_FAIL_UNREC      =  -6
const  SPGMR_PSOLVE_FAIL_REC      =  4
const  SPGMR_PSOLVE_FAIL_UNREC    =  -3
const  SPGMR_QRFACT_FAIL          =  3
const  SPGMR_QRSOL_FAIL           =  -5
const  SPGMR_RES_REDUCED          =  1
const  SPGMR_SUCCESS              =  0
const  SPTFQMR_ATIMES_FAIL_REC    =  4
const  SPTFQMR_ATIMES_FAIL_UNREC  =  -2
const  SPTFQMR_CONV_FAIL          =  2
const  SPTFQMR_MEM_NULL           =  -1
const  SPTFQMR_PSET_FAIL_REC      =  5
const  SPTFQMR_PSET_FAIL_UNREC    =  -4
const  SPTFQMR_PSOLVE_FAIL_REC    =  3
const  SPTFQMR_PSOLVE_FAIL_UNREC  =  -3
const  SPTFQMR_RES_REDUCED        =  1
const  SPTFQMR_SUCCESS            =  0


text:
module proj
using LinearAlgebra
# ! indicates that the variable is modified

function infnorm(x,n,γ)
# compute ∞-norm of x
# ---
# x : input vector 				(real n-vector)
# n : length of x
# γ : temp storage variable		(real scalar)    !

	γ = 0.0
	for j=1:n
		γ = max(γ,abs(x[j]))
	end
	return γ
end

function hlfspace!(x,z,u,ζ,n,β)
# projection onto halfspace { x | u⋅x ≤ ζ }
# ---
# x : projected result					(real n-vector) !
# z : vector to be projected			(real n-vector)
# u : normal to supporting hyperplane 	(real n-vector)
# ζ : scalar on RHS 					(real scalar)
# β : temporary variable                (real scalar)   !

	@assert BLAS.nrm2(n,u,1)≈1 "Normal vector u should have unit norm."

	β = BLAS.dot(n,z,1,u,1) - ζ

	if β > 0
		BLAS.blascopy!(n,z,1,x,1)
		axpy!(-β,u,x)
	else
		BLAS.blascopy!(n,z,1,x,1)
	end

end

function hlfspace2!(x,z,u1,u2,ζ1,ζ2,n,β,α1,α2)
# projection onto intersection of two halfspaces { x | u1⋅x ≤ ζ1 } ∩ { x | u2⋅x ≤ ζ2 }
# ---
# x  : projected result					(real n-vector) !
# z  : vector to be projected			(real n-vector)
# u1 : normal to supporting hyperplane 	(real n-vector)
# u2 : normal to supporting hyperplane 	(real n-vector)
# ζ1 : scalar on RHS 					(real scalar)
# ζ1 : scalar on RHS 					(real scalar)
# β  : temporary variable               (real scalar)   !
# α1 : temporary variable               (real scalar)   !
# α2 : temporary variable               (real scalar)   !

	@assert BLAS.nrm2(n,u1,1)≈1 "Normal vector u1 should have unit norm."
	@assert BLAS.nrm2(n,u2,1)≈1 "Normal vector u2 should have unit norm."

	β = BLAS.dot(n,u1,1,u2,1)
	α1 = BLAS.dot(n,z,1,u1,1)
	α2 = BLAS.dot(n,z,1,u2,1)

	if β ≈ 1
		β = min(ζ1,ζ2)
		if α1 > β
			BLAS.blascopy!(n,z,1,x,1)
			axpy!(β-α1,u1,x)
		else
			BLAS.blascopy!(n,z,1,x,1)
		end
	elseif β ≈ -1
		@assert ζ1+ζ2>0 "The two halfspaces do not intersect."
		if α1 < -ζ2
			BLAS.blascopy!(n,z,1,x,1)
			axpy!(-α1-ζ2,u1,x)
		elseif α1 > ζ1
			BLAS.blascopy!(n,z,1,x,1)
			axpy!(-α1+ζ1,u1,x)
		else
			BLAS.blascopy!(n,z,1,x,1)
		end
	else
		if α1-ζ1+β*(ζ2-α2)>0 && α2-ζ2+β*(ζ1-α1)>0
			BLAS.blascopy!(n,z,1,x,1)
			axpy!((α1-ζ1+β*(ζ2-α2))/(β*β-1),u1,x)
			axpy!((α2-ζ2+β*(ζ1-α1))/(β*β-1),u2,x)
		elseif α2>ζ2 && α1-ζ1+β*(ζ2-α2)≤0
			BLAS.blascopy!(n,z,1,x,1)
			axpy!(ζ2-α2,u2,x)
		elseif α1>ζ1 || α2>ζ2
			BLAS.blascopy!(n,z,1,x,1)
			axpy!(ζ1-α1,u1,x)
		else
			BLAS.blascopy!(n,z,1,x,1)
		end
	end

end

function box!(x,z,l,u,n)
# projection onto l∞-norm ball
# ---
# x : projected result			(real n-vector) !
# z : vector to be projected	(real n-vector)
# l : lower bound			 	(real scalar)
# u : upper bound 				(real scalar)
# n : length of z 				(integer)

	@assert u>l "Upper bound should be greater than lower bound."

	for j = 1:n
		x[j] = min(u,max(z[j],l))
	end

end

function ball!(x,z,r,n)
# projection onto l2-norm ball
# ---
# x : projected result			(real n-vector) !
# z : vector to be projected	(real n-vector)
# r : radius of l2-norm ball 	(real scalar)
# n : length of z 				(integer)

	@assert r>0 "Radius should be positive."

	BLAS.blascopy!(n,z,1,x,1)
	BLAS.scal!(n,r/max(BLAS.nrm2(n,x,1),r),x,1)

end

function soc!(x,z,α,n,β)
# projection onto second-order cone (aka soc, ice-cream, lorentz cone)
# type 1 : cone axis along last coordinate and arctan(α) is the cone angle
# ---
# x : projected result 		     (real n-vector) !
# z : vector to be projected	 (real n-vector)
# α : tan(cone angle)			 (real scalar)
# n : length of z (integer)		 (integer)
# β : temp storage variable  	 (real scalar)   !

	@assert α>0 "Cone angle arctan(α) must be positive."

	β = BLAS.nrm2(n-1,z,1)

	if β ≤ α*z[n]
		BLAS.blascopy!(n,z,1,x,1)
	elseif α*β + z[n] ≤ 0
		fill!(x,0.0)
	else
		BLAS.blascopy!(n-1,z,1,x,1)
		x[n] = (β*α+z[n])/(α*α+1)
		BLAS.scal!(n-1,α*x[n]/β,x,1)
	end

end

function soc2!(x,z,θ1,θ2,u,n,β,γ,κ)
# projection onto second-order cone (aka soc, ice-cream, lorentz cone)
# type 2 : cone axis along unit vector u and arccos(θ) is the cone angle
# ---
# x  : projected result 			(real n-vector) !
# z  : vector to be projected		(real n-vector)
# θ1 : cos(cone angle) 				(real scalar)
# θ2 : sin(cone angle) 				(real scalar)
# u  : unit vector along cone axis	(real n-vector)
# n  : length of z 					(integer)
# β  : temp storage variable 		(real scalar)   !
# γ  : temp storage variable 		(real scalar)   !
# κ  : temp storage variable 		(real n-vector) !

	@assert min(θ1,1-θ1)>0 && min(θ2,1-θ2)>0 "θ1,θ2 ∈ [0,1]"
	@assert BLAS.nrm2(n,u,1)≈1 "u must be a unit vector"

	β = BLAS.nrm2(n,z,1)
	γ = BLAS.dot(n,z,1,u,1)

	if θ1*β ≤ γ
		BLAS.blascopy!(n,z,1,x,1)
	elseif β*θ2 + γ ≤ 0
		fill!(x,0.0)
	else
		BLAS.blascopy!(n,u,1,x,1)
		BLAS.blascopy!(n,z,1,κ,1)
		axpy!(-γ,u,κ)
		axpby!(θ2/BLAS.nrm2(n,κ,1),κ,θ1,x)
		BLAS.scal!(n,BLAS.dot(n,z,1,x,1),x,1)
	end

end

function soc_ball!(x,z,α,r,n,β)
# projection onto intersection of soc cone and l2-norm ball
# ---
# x : projected result 			(real n-vector) !
# z : vector to be projected	(real n-vector)
# α : tan(cone angle) 			(real scalar)
# r : radius of l2-norm ball 	(real scalar)
# n : length of z 				(integer)
# β : temp storage variable 	(real scalar)   !

	soc!(x,z,α,n,β)
	ball!(x,x,r,n)

end

function soc2_ball!(x,z,θ1,θ2,u,r,n,β,γ,κ)
# projection onto intersection of soc cone and l2-norm ball
# ---
# x  : projected result 			(real n-vector) !
# z  : vector to be projected		(real n-vector)
# θ1 : cos(cone angle) 				(real scalar)
# θ2 : cos(cone angle) 				(real scalar)
# u  : unit vector along cone axis	(real n-vector)
# r  : radius of l2-norm ball 		(real scalar)
# n  : length of z 					(integer)
# β  : temp storage variable 		(real scalar)   !
# γ  : temp storage variable 		(real scalar)   !
# κ  : temp storage variable 		(real n-vector) !

	soc2!(x,z,θ1,θ2,u,n,β,γ,κ)
	ball!(x,x,r,n)

end

function admm!(x,y,u,z,f1,f2,α,β,ϵ,n,γ1,γ2,κ1,κ2)
# projection onto intersection of two sets with respective projections defined by f1 and f2
# ---
# x  : projected result (copy 1)	(real n-vector) !
# y  : projected result (copy 2)	(real n-vector) !
# u  : dual variable 				(real n-vector) !
# z  : vector to be projected		(real n-vector)
# f1 : projection function			(function: f!(θ,ψ) where θ is modified)
# f2 : projection function 			(function: f!(θ,ψ) where θ is modified)
# α  : ADMM step size				(real scalar)
# β  : 1/(1+α) 						(real scalar)
# ϵ  : ADMM termination criteria  	(real scalar)
# n  : length of z					(integer)
# γ  : temp storage variable        (real scalar)   !
# γ2 : temp storage variable 		(real scalar)   !
# κ1 : temp storage variable 		(real n-vector) !
# κ2 : temp storage variable 		(real n-vector) !

# x, y and u may hold values from previous call to admm! (warm-starting)

	f1(κ1,z)
	f2(κ2,z)
	κ1 .= κ1 .- z
	κ2 .= κ2 .- z
	γ1 = BLAS.nrm2(n,κ2,1)
	γ1 = max(γ1,BLAS.nrm2(n,κ1,1))
	if γ1 ≤ ϵ
	    BLAS.blascopy!(n,z,1,x,1)
		BLAS.blascopy!(n,z,1,y,1)
		# κ2[1] will be zero fortunately
		κ2[1] = 0.0 # inserted for type-stability
	else
		κ1 .= z .- u
		axpy!(α,y,κ1)
		BLAS.scal!(n,β,κ1,1)
		f1(x,κ1)
		κ2 .= z .+ u
		axpy!(α,x,κ2)
		BLAS.scal!(n,β,κ2,1)
		f2(y,κ2)
		κ1 .= x .- y
		axpy!(α,κ1,u)

		κ2[1] = infnorm(κ1,n,γ2)
		γ1 = 1.0
		while κ2[1]≥ϵ
			κ1 .= z .- u
			axpy!(α,y,κ1)
			BLAS.scal!(n,β,κ1,1)
			f1(x,κ1)
			κ2 .= z .+ u
			axpy!(α,x,κ2)
			BLAS.scal!(n,β,κ2,1)
			f2(y,κ2)
			κ1 .= x .- y
			axpy!(α,κ1,u)

			γ1 = γ1 + 1.0 # iteration counter

			κ2[1] = infnorm(κ1,n,γ2)
		end
		κ2[1] = γ1 # store the iteration count
	end

end

function multi_admm!(x,y,u,z,f,α,β,ϵ,n,N,γ1,γ2,ψ,κ)
# projection onto intersection of N sets with projections defined by f[i], i=1,...,N
# ---
# x  : projected results (copy 1)	(N-vector of real n-vectors) !
# y  : projected results (copy 2)	(N-vector of real n-vectors) !
# u  : dual variable 				(N-vector of real n-vectors) !
# z  : vector to be projected		(real n-vector)
# f  : N projection functions   	(N-vector of functions: f[i]!(θ,ψ) where θ is modified, i=1,...,N)
# α  : ADMM step-size				(positive scalar)
# β  : 1/(2α(N-1)+1)				(positive scalar)
# n  : length of z 					(integer)
# N  : no. of sets					(integer)
# γ1 : temp storage variable  		(real scalar)                !
# γ2 : temp storage variable  		(real scalar)                !
# ψ  : temp storage variable  		(real scalar)                !
# κ  : temp storage variable        (real n-vector)              !

	@assert length(f)==N "No. of projection functions is not N."
	# @assert N ≥ 3 "Use admm! for fewer than 3 sets."

	# γ1 = 100.0
	# ψ = 0.0 # iteration counter

	for i = 1:N
		κ .= x[i] # BLAS.blascopy!(n,x[i],1,κ,1)
		lmul!(N-2,κ)
		for j = 1:N
			κ .= κ .+ x[j]
		end
		lmul!(α,κ)
		κ .= κ .- u[i] .+ z
		lmul!(β,κ)
		f[i](y[i],κ)
	end

	γ1 = 0.0 # holds termination criteria value
	x[1] .= y[1]
	for i = 2:N
		x[i] .= y[i]
		κ .= x[1] .- x[i]

		γ1 = max(γ1,infnorm(κ,n,γ2))
	end

	for i = 1:N
		κ .= x[i]
		lmul!(N,κ)
		for j = 1:N
			κ .= κ .- x[j]
		end
		axpy!(α,κ,u[i])
	end

	ψ = 1.0

	while γ1≥ϵ
		for i = 1:N
			κ .= x[i] # BLAS.blascopy!(n,x[i],1,κ,1)
			lmul!(N-2,κ)
			for j = 1:N
				κ .= κ .+ x[j]
			end
			lmul!(α,κ)
			κ .= κ .- u[i] .+ z
			lmul!(β,κ)
			f[i](y[i],κ)
		end

		γ1 = 0.0 # holds termination criteria value
		x[1] .= y[1]
		for i = 2:N
			x[i] .= y[i]
			κ .= x[1] .- x[i]

			γ1 = max(γ1,infnorm(κ,n,γ2))
		end

		for i = 1:N
			κ .= x[i]
			lmul!(N,κ)
			for j = 1:N
				κ .= κ .- x[j]
			end
			axpy!(α,κ,u[i])
		end

		ψ = ψ + 1.0
	end
	κ[1] = ψ

end

end

max_stars_repo_name: purnanandelango/pipg-demo

max_issues_repo_name: UW-ACL/pipg-demo

max_forks_repo_name: UW-ACL/pipg-demo


text:
export swish

doc"""
    Swish

Swish activation function.

```math
f(x) = x \cdot \sigma (\beta x)
```
where ``\beta`` is a leanable parameter.

# References
* Ramachandran et al. ["Searching for Activation Functions"](https://arxiv.org/abs/1710.05941), arXiv 2017.
"""
struct Swish
    beta::Var
end
Swish(::Type{T}) where T = Swish(zerograd(ones(T,1)))
(f::Swish)(x) = swish(x, f.beta)

swish(x::Var, beta::Var) = Var(swish.(x.data,beta.data), (swish,x,beta))
swish(x::T, beta::T) where T = x * sigmoid(beta*x)

function addgrad!(y::Var, ::typeof(swish), x::Var, beta::Var)
    isvoid(x.grad) || ∇swish_x!(y.data, y.grad, x.data, x.grad, beta.data)
    isvoid(beta.grad) || ∇swish_beta!(y.data, y.grad, x.data, beta.data, beta.grad)
end

function ∇swish_x!(y::Array{T}, gy::Array{T}, x::Array{T}, gx::Array{T}, beta::Vector{T}) where T
    @inbounds for i = 1:length(x)
        s = sigmoid(beta[1] * x[i])
        gx[i] += gy[i] * (s + beta[1]*x[i] * s * (1-s))
    end
end

function ∇swish_beta!(y::Array{T}, gy::Array{T}, x::Array{T}, beta::Vector{T}, gbeta::Vector{T}) where T
    @inbounds for i = 1:length(gy)
        s = sigmoid(beta[1] * x[i])
        gbeta[1] += gy[i] * x[i] * x[i] * s * (1-s)
    end
end


text:
using PowerSystems
using JuMP
using PowerSimulations
const PS = PowerSimulations

# ED Testing
base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir,"data/data_5bus.jl"))
sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing, 100.0);
#simple_reserve = PowerSystems.StaticReserve("test_reserve",sys5.generators.thermal,60.0,[gen.tech for gen in sys5.generators.thermal])

# ED with thermal gen, static load, copper plate
@test try
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    #JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
true finally end

# ED with thermal and curtailable renewable gen, static load, copper plate
@test try
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    #JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
true finally end

# ED with thermal and fixed renewable gen, interruptable load, copper plate
@test try
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation = PS.ThermalDispatch),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            [(device = ElectricLoad, formulation = PS.InterruptibleLoad)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    #JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
true finally end

# ED with thermal gen, copper plate, and reserve
@test try
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            [(service = reserve5, formulation = PS.RampLimitedReserve)],
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    #JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
    #ED.model.moi_backend.model.optimizer.termination_status
true finally end

# UC Testing
base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir,"data/data_5bus_uc.jl"))
sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing, 100.0);

# UC with thermal gen, static load, copper plate
@test try
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    #JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
true finally end

# UC with thermal and curtailable renewable gen, static load, copper plate
@test try
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    #JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
true finally end

# UC with thermal and fixUC renewable gen, interruptable load, copper plate
@test try
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation = PS.StandardThermalCommitment),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            [(device = ElectricLoad, formulation = PS.InterruptibleLoad)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    #JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
true finally end

# UC with thermal gen, copper plate, and reserve
@test try
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            [(service = reserve5, formulation = PS.RampLimitedReserve)],
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    #JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
    #UC.model.moi_backend.model.optimizer.termination_status
true finally end


text:
using Bio3DView
using Test