text:
#EasyPlotGrace: Render EasyPlot-plots with Grace/xmgrace (through GracePlot)
#-------------------------------------------------------------------------------
module EasyPlotGrace

#Hoping that environment of including module points to *correct* CMDimData:
import CMDimData
#Ensure *correct* modules are imported:
import CMDimData.EasyPlot #Avoid name collisions
using CMDimData.MDDatasets
using CMDimData.Colors
using GracePlot

import CMDimData.EasyPlot: Optional, Target, DS
import GracePlot: LineAttributes, GlyphAttributes

include("base.jl")
include("defaults.jl")
include("builder.jl")

#==Initialization
===============================================================================#
function __init__()
	global defaults
	_initialize(defaults)

	EasyPlot.register(:EasyPlotGrace)
	return
end

end

#Last line


text:
#=

  Curious set of integers in Julia ConstraintSolver.jl

  Martin Gardner (February 1967):
  """
  The integers 1,3,8, and 120 form a set with a remarkable property:
  the product of any two integers is one less than a perfect square.
  Find a fifth number that can be added to the set without destroying
  this property.
  """

  Solution: The number is 0.

  There are however other sets of five numbers with this property.
  Here are the one in the range of 0..10000:

  [0, 1, 3, 8, 120]
  [0, 1, 3, 120, 1680]
  [0, 1, 8, 15, 528]
  [0, 1, 8, 120, 4095]
  [0, 1, 15, 24, 1520]
  [0, 1, 24, 35, 3480]
  [0, 1, 35, 48, 6888]
  [0, 2, 4, 12, 420]
  [0, 2, 12, 24, 2380]
  [0, 2, 24, 40, 7812]
  [0, 3, 5, 16, 1008]
  [0, 3, 8, 21, 2080]
  [0, 3, 16, 33, 6440]
  [0, 4, 6, 20, 1980]
  [0, 4, 12, 30, 5852]
  [0, 5, 7, 24, 3432]
  [0, 6, 8, 28, 5460]
  [0, 7, 9, 32, 8160]

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function curious_set_of_integers(print_solutions=true,all_solutions=false,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 5
    dom_max = 120
    @variable(model, 0 <= x[1:n] <= dom_max, Int)

    @constraint(model, x in CS.AllDifferent())
    increasing_strict(model, x)

    # It's faster without these hints
    # @variable(model, 1 <= ixs[1:n] <= n, Int) # indices in X
    # @constraint(model, ixs in CS.AllDifferent())
    # my_element(model,ixs[1],x,1)
    # my_element(model,ixs[2],x,3)
    # my_element(model,ixs[3],x,8)
    # my_element(model,ixs[4],x,120)

    table = resize_matrix([ [i,j,i * j] for i in 0:dom_max, j in 0:dom_max])
    ps = []
    for i in 1:n, j in 1:i-1
        p = @variable(model, integer=true, lower_bound=0, upper_bound=dom_max)
        # @constraint(model, p*p -1 == x[i]*x[j]) # non-linear constraints!

        p_sq = @variable(model, integer=true, lower_bound=0, upper_bound=dom_max^2)
        @constraint(model, [p, p, p_sq] in CS.TableSet(table))

        xij_sq = @variable(model, integer=true, lower_bound=0, upper_bound=dom_max^2)
        @constraint(model, [x[i], x[j], xij_sq] in CS.TableSet(table))

        @constraint(model,p_sq-1 == xij_sq)

        push!(ps,p)
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time curious_set_of_integers(true,false)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
isdefined(Base, :__precompile__) && __precompile__()

module FredData

using Requests
using DataFrames
using TimeZones
import Requests: get
import JSON

export
       # Fred object
       Fred, get_api_url, set_api_url!, get_api_key,

       # FredSeries object
       FredSeries, id, title, units_short, units, seas_adj_short, seas_adj, freq_short,
       freq, realtime_start, realtime_end, last_updated, notes, trans_short, df,

       # Download data
       get_data

const MAX_ATTEMPTS       = 3
const FIRST_REALTIME     = Date(1776,07,04)
const LAST_REALTIME      = Date(9999,12,31)
const FRED_TIME_ZONE     = TimeZone("America/Chicago")
const EARLY_VINTAGE_DATE = "1991-01-01"
const DEFAULT_API_URL    = "https://api.stlouisfed.org/fred/"
const API_KEY_LENGTH     = 32
const KEY_ENV_NAME       = "FRED_API_KEY"
const KEY_FILE_NAME      = ".freddatarc"

# Fred connection type
"""
A connection to the Fred API.

Constructors
------------
- `Fred()`: Key detected automatically. First, looks for the environment variable
    FRED_API_KEY, then looks for the file ~/.freddatarc.
- `Fred(key::AbstractString)`: User specifies key

Arguments
---------
- `key`: Registration key provided by the Fred.

Notes
-----
- Set the API url with `set_api_url!(f::Fred, url::AbstractString)`
"""
type Fred
    key::AbstractString
    url::AbstractString
end
Fred(key) =  Fred(key, DEFAULT_API_URL)
function Fred()
    key = ""
    if KEY_ENV_NAME in keys(ENV)
        key = ENV[KEY_ENV_NAME]
    elseif isfile(joinpath(homedir(), KEY_FILE_NAME))
        open(joinpath(homedir(), KEY_FILE_NAME), "r") do file
            key = readall(file)
        end
        key = rstrip(key)
    else
        error("FRED API Key not detected.")
    end

    @printf "API key loaded.\n"

    # Key validation
    if length(key) > API_KEY_LENGTH
        key = key[1:API_KEY_LENGTH]
        warn("Key too long. First ", API_KEY_LENGTH, " chars used.")
    elseif length(key) < API_KEY_LENGTH
        error("Invalid FRED API key: ", key, ". Key too short.")
    end
    if !isxdigit(key)
        error("Invalid FRED API key: ", key, ". Invalid characters.")
    end

    return Fred(key)
end
get_api_key(f::Fred) = f.key
get_api_url(f::Fred) = f.url
set_api_url!(f::Fred, url::AbstractString) = setfield!(f, :url, url)

function Base.show(io::IO, f::Fred)
    @printf io "FRED API Connection\n"
    @printf io "\turl: %s\n" get_api_url(f)
    @printf io "\tkey: %s\n" get_api_key(f)
end

"""
```
FredSeries(...)
```

Represent a single data series, and all associated metadata, return from Fred.

### Field access
- `id(f)`
- `title(f)`
- `units_short(f)`
- `units(f)`
- `seas_adj_short(f)`
- `seas_adj(f)`
- `freq_short(f)`
- `freq(f)`
- `realtime_start(f)`
- `realtime_end(f)`
- `last_updated(f)`
- `notes(f)`
- `trans_short(f)`
- `df(f)`

"""
immutable FredSeries
    # From series query
    id::AbstractString
    title::AbstractString
    units_short::AbstractString
    units::AbstractString
    seas_adj_short::AbstractString
    seas_adj::AbstractString
    freq_short::AbstractString
    freq::AbstractString
    realtime_start::AbstractString
    realtime_end::AbstractString
    last_updated::DateTime
    notes::AbstractString

    # From series/observations query
    trans_short::AbstractString # "units"
    df::DataFrames.DataFrame
end
id(f::FredSeries) = f.id
title(f::FredSeries) = f.title
units_short(f::FredSeries) = f.units_short
units(f::FredSeries) = f.units
seas_adj_short(f::FredSeries) = f.seas_adj_short
seas_adj(f::FredSeries) = f.seas_adj
freq_short(f::FredSeries) = f.freq_short
freq(f::FredSeries) = f.freq
realtime_start(f::FredSeries) = f.realtime_start
realtime_end(f::FredSeries) = f.realtime_end
last_updated(f::FredSeries) = f.last_updated
function notes(f::FredSeries)
    str = f.notes
    str = replace(str, r"[\r\n]", " ")
    str = replace(str, r" +", " ")
    str = strip(str)
    return str
end
trans_short(f::FredSeries) = f.trans_short
df(f::FredSeries) = f.df

function Base.show(io::IO, s::FredSeries)
    @printf io "FredSeries\n"
    @printf io "\tid: %s\n"                id(s)
    @printf io "\ttitle: %s\n"             title(s)
    @printf io "\tunits: %s\n"             units(s)
    @printf io "\tseas_adj (native): %s\n" seas_adj(s)
    @printf io "\tfreq (native): %s\n"     freq(s)
    @printf io "\trealtime_start: %s\n"    realtime_start(s)
    @printf io "\trealtime_end: %s\n"      realtime_end(s)
    @printf io "\tlast_updated: %s\n"      last_updated(s)
    @printf io "\tnotes: %s\n"             notes(s)
    @printf io "\ttrans_short: %s\n"       trans_short(s)
    @printf io "\tdf: %dx%d DataFrame with columns %s\n" size(df(s))... names(df(s))
end

include("get_data.jl")

end # module


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
theme(:wong)
gr(aspect_ratio=1,legend=nothing,axis=nothing,border=nothing,ticks=nothing,
   markerstrokewidth=0,markersize=1)
r1D, w1D = gauss_quad(0,0,2)
Nfp = length(r1D) # number of points per face
e = ones(Nfp) # vector of all ones
z = zeros(Nfp) # vector of all zeros
rf = [r1D; -r1D; -e];
sf = [-e; r1D; -r1D];

r_tri,s_tri,_ = Basis2DTri.quad_nodes_tri(4)
# Front triangle
# plot([-1;-1],[-1;1],seriescolor=:black,lw=3,label = nothing)
# plot!([-1;1],[1;-1],seriescolor=:black,lw=3,label = nothing)
# plot!([-1;1],[-1;-1],seriescolor=:black,lw=3,label = nothing)
# plot!([1;4],[-1,1],seriescolor=:black,lw=3,label = nothing)
# plot!([-1;2],[-1,1],seriescolor=:black,lw=3,ls=:dash,label = nothing)
# plot!([-1;2],[1,3],seriescolor=:black,lw=3,label = nothing)

plot([-1;-1].+3/4,[-1;1].+2/4,seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+3/4,[1;-1].+2/4,seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+3/4,[-1;-1].+2/4,seriescolor=:black,lw=3,label = nothing)
scatter!(r_tri.+3/4,s_tri.+2/4,markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(rf.+3/4,sf.+2/4,markersize=7,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)

# plot!([-1;-1].+3/4*2,[-1;1].+2/4*2,seriescolor=:black,lw=3,ls=:dash,label = nothing)
# plot!([-1;1].+3/4*2,[1;-1].+2/4*2,seriescolor=:black,lw=3,ls=:dash,label = nothing)
# plot!([-1;1].+3/4*2,[-1;-1].+2/4*2,seriescolor=:black,lw=3,ls=:dash,label = nothing)
# scatter!(r_tri.+3/4*2,s_tri.+2/4*2,markersize=5,markercolor=:white,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(rf.+3/4*2,sf.+2/4*2,markersize=3,markershape=:rect,markercolor=:white,markerstrokecolor=:black,markerstrokewidth=2)

# plot!([-1;-1].+3/4*3,[-1;1].+2/4*3,seriescolor=:black,lw=3,ls=:dash,label = nothing)
# plot!([-1;1].+3/4*3,[1;-1].+2/4*3,seriescolor=:black,lw=3,ls=:dash,label = nothing)
# plot!([-1;1].+3/4*3,[-1;-1].+2/4*3,seriescolor=:black,lw=3,ls=:dash,label = nothing)
# scatter!(r_tri.+3/4*3,s_tri.+2/4*3,markersize=5,markercolor=:white,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(rf.+3/4*3,sf.+2/4*3,markersize=3,markershape=:rect,markercolor=:white,markerstrokecolor=:black,markerstrokewidth=2)

# plot!([-1;-1].+3,[-1;1].+2,seriescolor=:black,lw=3,label = nothing)
# plot!([-1;1].+3,[1;-1].+2,seriescolor=:black,lw=3,label = nothing)
# plot!([-1;1].+3,[-1;-1].+2,seriescolor=:black,lw=3,label = nothing)

# # savefig("reference_tri.png")
# scatter!(r_tri.+3,s_tri.+2,markersize=5,markercolor=:white,markerstrokecolor=:black,markerstrokewidth=2)
# scatter!(rf.+3,sf.+2,markersize=3,markershape=:rect,markercolor=:white,markerstrokecolor=:black,markerstrokewidth=2)

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_circles(fname)
    pagewidth, pageheight = 1200, 1400

    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setopacity(0.5)
    gsave()
    t = Tiler(1000, 1000, 20, 20)
    for (pos, n) in t
        randomhue()
        box(pos, rand(5:50), rand(5:50), rand(5:55), :fillpreserve)
        randomhue()
        rotate(2pi * rand())
        strokepath()
        # p = pathtopoly()
        # for pl in p
        #     prettypoly(pl, :fill, () -> circle(O, 1, :fill))
        # end
    end

    grestore()

    sethue("black") # hide
    setline(4)
    # round corners with different radii
    box(O, 200, 150, collect(range(5, 20, step=4)), :stroke)

    for (pos, n) in t
        randomhue()

        polysmooth(box(O, 200, 150, vertices=true), 10, :stroke)

        box(pos, rand(5:50), rand(5:50), rand(5:15), :path)
        randomhue()
        strokepath()
        p = pathtopoly()
        for pl in p
            prettypoly(pl, :fill, () -> circle(O, 5, :stroke))
        end
    end

    @test finish() == true
    println("...finished boxtest, saved in $(fname)")
end

fname = "box-test.png"
test_circles(fname)


text:
# boundary functional evaluation and derivatives

include("boundary_entropy_diss.jl")

@doc """
### EulerEquationMod.evalFunctional

High level function that evaluates all the functionals specified over
various edges. This function is agnostic to the type of the functional being
computed and calls a mid level functional-type specific function for the actual
evaluation.

**Arguments**

*  `mesh` :  Abstract mesh object
*  `sbp`  : Summation-By-Parts operator
*  `eqn`  : Euler equation object
*  `opts` : Options dictionary
*  `func` : Object of type AbstractBoundaryFunctional. This is type is associated
                      with the functional being computed and holds all the
                      relevant data.

**Outputs**

 * val: functional value
"""->
function _evalFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol}, opts,
            func::AbstractBoundaryFunctional) where {Tmsh, Tsol}

  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  # Calculate functional over edges
  val = calcBndryFunctional(mesh, sbp, eqn, opts, func)

  return val
end

@doc """
### EulerEquationMod.evalFunctional_revm

Reverse mode of EulerEquationMod.evalFunctional, It takes in functional value
and return `mesh.nrm_bndry_bar`. Different functionals will need to be added
to the if statement to further extend this function.

**Arguments**

*  `mesh` :  Abstract mesh object
*  `sbp`  : Summation-By-Parts operator
*  `eqn`  : Euler equation object
*  `opts` : Options dictionary
*  `func` : Object of type AbstractBoundaryFunctional. This is type is associated
                      with the functional being computed and holds all the
                      relevant data.
*  `functionalName` : Name of the functional being evaluated.

"""->
function _evalFunctionalDeriv_m(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::AbstractSolutionData{Tsol}, opts,
                           func::AbstractBoundaryFunctional,
                           val_bar::Number=1,
                           ) where {Tmsh, Tsol}

  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  calcBndryFunctional_revm(mesh, sbp, eqn, opts, func, val_bar)
  return nothing
end

"""
  This is the general case function,  It calls
  [`calcBoundaryFunctionalIntegrand`](@ref) to get the functional integrand at
  each boundary node and computes the integral.
"""
function calcBndryFunctional(mesh::AbstractDGMesh{Tmsh},
                      sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                      opts, func::AbstractBoundaryFunctional,
                      ) where {Tmsh, Tsol, Tres, Tdim}

  functional_val = zeros(Tres, 1)
  node_info = zeros(Int, 3)

  # loop over boundary conditions that have this functional
  for itr = 1:length(func.bcnums)
    bcnum = func.bcnums[itr]

    start_index = mesh.bndry_offsets[bcnum]
    end_index = mesh.bndry_offsets[bcnum+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
    integrand = zeros(Tsol, 1, mesh.sbpface.numnodes, nfaces)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.sbpface.numnodes
        q = ro_sview(eqn.q_bndry, :, j, global_facenum)
#        convertToConservative(eqn.params, q, q2)
        aux_vars = ro_sview(eqn.aux_vars_bndry, :, j, global_facenum)
        x = ro_sview(mesh.coords_bndry, :, j, global_facenum)
        phys_nrm = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
        node_info[1] = itr; node_info[2] = j; node_info[3] = i

        integrand[1, j, i] = calcBoundaryFunctionalIntegrand(eqn.params, q,
                                           aux_vars, phys_nrm, node_info, func)
      end  # End for j = 1:mesh.sbpface.numnodes
    end    # End for i = 1:nfaces

    bndry_facenums_arr = sview(mesh.bndryfaces, idx_range)
    integratefunctional!(mesh.sbpface, bndry_facenums_arr, integrand, functional_val)
  end  # end loop itr

  # global sum
  val = MPI.allreduce(functional_val, MPI.SUM, eqn.comm)[1]

  return val
end  # function calcBndryFunctional

function calcBndryFunctional_revm(mesh::AbstractDGMesh{Tmsh},
     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
     opts, func::AbstractBoundaryFunctional, _val_bar::Number=1,
     ) where {Tmsh, Tsol, Tres, Tdim}

  functional_val = zeros(Tres, 1)
  node_info = zeros(Int, 3)

  val_bar = Array{Tres}(1)
  val_bar[1] = _val_bar
#  val_bar = Tres[_val_bar]  # value is implicitly known, no need to reverse allreduce

  # loop over boundary conditions that have this functional
  for itr = 1:length(func.bcnums)
    bcnum = func.bcnums[itr]

    start_index = mesh.bndry_offsets[bcnum]
    end_index = mesh.bndry_offsets[bcnum+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
    integrand_bar = zeros(Tsol, 1, mesh.sbpface.numnodes, nfaces)

    bndry_facenums_arr = sview(mesh.bndryfaces, idx_range)
    integratefunctional_rev!(mesh.sbpface, bndry_facenums_arr, integrand_bar, val_bar)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.sbpface.numnodes
        q = ro_sview(eqn.q_bndry, :, j, global_facenum)
#        convertToConservative(eqn.params, q, q2)
        aux_vars = ro_sview(eqn.aux_vars_bndry, :, j, global_facenum)
        x = ro_sview(mesh.coords_bndry, :, j, global_facenum)
        phys_nrm = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
        nrm_bar = sview(mesh.nrm_bndry_bar, :, j, global_facenum)
        node_info[1] = itr; node_info[2] = j; node_info[3] = i

        calcBoundaryFunctionalIntegrand_revm(eqn.params, q,  aux_vars, phys_nrm,
                              node_info, func, nrm_bar, integrand_bar[1, j, i])
      end  # End for j = 1:mesh.sbpface.numnodes
    end    # End for i = 1:nfaces

  end  # end loop itr

  return
end  # function calcBndryFunctional

function calcBndryFunctional(mesh::AbstractDGMesh{Tmsh},
     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
     opts, func::Union{LiftCoefficient, DragCoefficient}
    ) where {Tmsh, Tsol, Tres, Tdim}

  val = calcBndryFunctional(mesh, sbp, eqn, opts, func.func)
  fac = 0.5*eqn.params.rho_free*eqn.params.Ma*eqn.params.Ma

  val = val/fac

  return val
end

function calcBndryFunctional(mesh::AbstractDGMesh{Tmsh},
     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
     opts, func::SolutionDeviation) where {Tmsh, Tsol, Tres, Tdim}

  u_bar = zeros(Tsol, mesh.numDofPerNode)
  val = zero(Tres)
  for i=1:mesh.numEl

    # compute u_bar
    vol = zero(Tmsh)  # area/volume of the element
    fill!(u_bar, 0)
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      vol += w_j
      for k=1:mesh.numDofPerNode
        u_bar[k] += eqn.q[k, j, i]*w_j
      end
    end

    for k=1:mesh.numDofPerNode
      u_bar[k] /= vol
    end

    # compute \int || u - u_bar ||
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      for k=1:mesh.numDofPerNode
        delta_u = eqn.q[k, j, i] - u_bar[k]
        val += delta_u*w_j*delta_u
      end
    end
  end  # end i

  return val
end

function calcBndryFunctional_revm(mesh::AbstractDGMesh{Tmsh},
                       sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                       opts, func::AeroCoefficients, val_bar::Number=1,
                       ) where {Tmsh, Tsol, Tres, Tdim}

  fac = 0.5*eqn.params.rho_free*eqn.params.Ma*eqn.params.Ma
  val_bar = val_bar/fac

  calcBndryFunctional_revm(mesh, sbp, eqn, opts, func.func, val_bar)

  return nothing
end

#TODO: reverse mode for revq

function calcBndryFunctional_revm(mesh::AbstractDGMesh{Tmsh},
                       sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                       opts, func::SolutionDeviation, val_bar::Number=1,
                       ) where {Tmsh, Tsol, Tres, Tdim}

  u_bar = zeros(Tsol, mesh.numDofPerNode)
  u_bar_bar = zeros(Tsol, mesh.numDofPerNode)
  for i=1:mesh.numEl

    # compute u_bar
    vol = zero(Tmsh)  # area/volume of the element
    fill!(u_bar, 0)
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      vol += w_j
      for k=1:mesh.numDofPerNode
        u_bar[k] += eqn.q[k, j, i]*w_j
      end
    end

    for k=1:mesh.numDofPerNode
      u_bar[k] /= vol
    end
#=
    # compute \int || u - u_bar ||
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      for k=1:mesh.numDofPerNode
        delta_u = eqn.q[k, j, i] - u_bar[k]
        val += delta_u*w_j*delta_u
      end
    end
=#
    #-------------------------------
    # reverse sweep
    fill!(u_bar_bar, 0)
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      w_j_bar = zero(Tmsh)
      for k=1:mesh.numDofPerNode
        delta_u = eqn.q[k, j, i] - u_bar[k]

        u_bar_bar[k] += 2*delta_u*w_j*val_bar
        w_j_bar      += delta_u*delta_u*val_bar
      end
      mesh.jac_bar[j, i] += -sbp.w[j]*w_j_bar/(mesh.jac[j, i]*mesh.jac[j, i])
    end

    vol_bar = zero(Tmsh)
    for k=1:mesh.numDofPerNode
      u_bar[k] *= vol  # restore primal value
      vol_bar      += -u_bar[k]*u_bar_bar[k]/(vol*vol)
      u_bar_bar[k] /= vol
    end

    for j=1:mesh.numNodesPerElement
      #w_j = sbp.w[j]/mesh.jac[j, i]
      w_j_bar = zero(Tmsh)
      for k=1:mesh.numDofPerNode
        w_j_bar += eqn.q[k, j, i]*u_bar_bar[k]
      end
      w_j_bar += vol_bar
      mesh.jac_bar[j, i] += -sbp.w[j]*w_j_bar/(mesh.jac[j, i]*mesh.jac[j, i])
    end

  end  # end i

  return nothing
end

@doc """
### EulerEquationMod.calcBoundaryFunctionalIntegrand

Computes the integrand for boundary functional at a surface SBP node. Every
functional of a different type may need a corresponding method to compute the
integrand. The type of the functional object, which is a subtype of
`AbstractBoundaryFunctional`.

**Arguments**

*  `params` : eqn.params object
*  `q` : Nodal solution
*  `aux_vars` : Auxiliary variables
*  `nrm` : Face normal vector in the physical space
*  `node_info` : Information about the SBP node
*  `objective` : Functional data type
*  `val` : Function output value

"""->
function calcBoundaryFunctionalIntegrand(params::ParamType{2},
                       q::AbstractArray{Tsol,1},
                       aux_vars::AbstractArray{Tres, 1},
                       nrm::AbstractArray{Tmsh},
                       node_info::AbstractArray{Int},
                       objective::BoundaryForceData,
                       ) where {Tsol, Tres, Tmsh}

  # Compute the numerical flux for the euler equation and extract the X & Y
  # momentum values. The normal vector supplied has already been converted
  # to the physical space from the parametric space.

  euler_flux = objective.euler_flux # Reuse existing memory

  fac = 1.0/(sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2]))
  # normalize normal vector
  nx = nrm[1]*fac
  ny = nrm[2]*fac

  normal_momentum = nx*q[2] + ny*q[3]

  qg = objective.qg
  for i=1:length(q)
    qg[i] = q[i]
  end
  qg[2] -= nx*normal_momentum
  qg[3] -= ny*normal_momentum

  calcEulerFlux(params, qg, aux_vars, nrm, euler_flux)

  return euler_flux[2]*objective.facx + euler_flux[3]*objective.facy
end # End calcBoundaryFunctionalIntegrand 2D

function calcBoundaryFunctionalIntegrand(params::ParamType{3},
                       q::AbstractArray{Tsol,1},
                       aux_vars::AbstractArray{Tres, 1},
                       nrm::AbstractArray{Tmsh},
                       node_info::AbstractArray{Int},
                       objective::BoundaryForceData,
                       ) where {Tsol, Tres, Tmsh}

  fac = 1.0/(sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3]))
  # normalize normal vector
  nx = nrm[1]*fac
  ny = nrm[2]*fac
  nz = nrm[3]*fac

  normal_momentum = nx*q[2] + ny*q[3] + nz*q[4]
  qg = objective.qg
  for i=1:length(q)
    qg[i] = q[i]
  end
  qg[2] -= nx*normal_momentum
  qg[3] -= ny*normal_momentum
  qg[4] -= nz*normal_momentum

  euler_flux = objective.euler_flux # Reuse existing memory
  calcEulerFlux(params, qg, aux_vars, nrm, euler_flux)

  val = euler_flux[2]*objective.facx + euler_flux[3]*objective.facy +
        euler_flux[4]*objective.facz
  return val
end # End calcBoundaryFunctionalIntegrand 3D

function calcBoundaryFunctionalIntegrand(params::ParamType{2},
                       q::AbstractArray{Tsol,1},
                       aux_vars::AbstractArray{Tres, 1},
                       nrm::AbstractArray{Tmsh},
                       node_info::AbstractArray{Int},
                       objective::MassFlowData,
                       ) where {Tsol, Tres, Tmsh}

  val = q[2]*nrm[1] + q[3]*nrm[2]

  return val
end

function calcBoundaryFunctionalIntegrand(params::ParamType{Tdim},
                       q::AbstractArray{Tsol,1},
                       aux_vars::AbstractArray{Tres, 1},
                       nrm::AbstractArray{Tmsh},
                       node_info::AbstractArray{Int},
                       objective::EntropyFluxData,
                       ) where {Tsol, Tres, Tmsh, Tdim}

  # compute u dot n
  unet = zero(Tres)
  for i=1:Tdim
    unet += q[i+1]*nrm[i]
  end
  unet/q[1]

  p = calcPressure(params, q)
  s = log(p) - params.gamma*log(q[1])
  U = -q[1]*s/params.gamma_1

  return unet*U
end

@doc """
### EulerEquationMod. calcBoundaryFunctionalIntegrand_revm

Reverse mode for boundary functional integrand w.r.t. nrm. Takes in input
val_bar and return nrm_bar for further reverse propagation.

**Arguments**

*  `params` : eqn.params object
*  `q` : Nodal solution
*  `aux_vars` : Auxiliary variables
*  `nrm` : Face normal vector in the physical space
*  `node_info` : Information about the SBP node
*  `objective` : Functional data type
*  `nrm_bar` : Resulting vector
*  `val_bar` : Nodal portion of the seeding vector

"""->

function calcBoundaryFunctionalIntegrand_revm(params::ParamType{2},
                                         q::AbstractArray{Tsol,1},
                                         aux_vars::AbstractArray{Tres, 1},
                                         nrm::AbstractArray{Tmsh},
                                         node_info::AbstractArray{Int},
                                         objective::BoundaryForceData,
                                         nrm_bar::AbstractArray{Tmsh,1},
                                         val_bar::Number) where {Tsol, Tres, Tmsh}

  #---- Forward sweep
  fac = 1.0/(sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2]))
  nx = nrm[1]*fac # Normalized unit vectors
  ny = nrm[2]*fac #
  normal_momentum = nx*q[2] + ny*q[3]
  qg = objective.qg
  for i=1:length(q)
    qg[i] = q[i]
  end
  qg[2] -= nx*normal_momentum
  qg[3] -= ny*normal_momentum

  #---- Reverse Sweep
  euler_flux_bar = zeros(Tsol, 4) # For 2D
  qg_bar = zeros(Tsol, 4)
  q_bar = zeros(Tsol,4)

  # Reverse diff val = euler_flux[2]*objective.facx + euler_flux[3]*objective.facy
  euler_flux_bar[2] = val_bar*objective.facx
  euler_flux_bar[3] = val_bar*objective.facy

  # Reverse diff calcEulerFlux
  calcEulerFlux_revm(params, qg, aux_vars, nrm, nrm_bar, euler_flux_bar)
  calcEulerFlux_revq(params, qg, qg_bar, aux_vars, nrm, euler_flux_bar)
  ny_bar = zero(Tres)               # Initialize
  nx_bar = zero(Tres)               #
  normal_momentum_bar = zero(Tsol)  #

  # Reverse diff qg[3] -= ny*normal_momentum
  ny_bar -= qg_bar[3]*normal_momentum
  normal_momentum_bar -= qg_bar[3]*ny
  qg_bar[3] += qg_bar[3]

  # Reverse diff qg[2] -= nx*normal_momentum
  nx_bar -= qg_bar[2]*normal_momentum
  normal_momentum_bar -= qg_bar[2]*nx
  qg_bar[2] += qg_bar[2]

  # Reverse diff qg[:] = q[:]
  for i=1:length(q_bar)
    q_bar[i] += qg_bar[i]
  end

  # Reverse diff normal_momentum = nx*q[2] + ny*q[3]
  nx_bar += normal_momentum_bar*q[2]
  ny_bar += normal_momentum_bar*q[3]
  q_bar[2] += normal_momentum_bar*nx
  q_bar[3] += normal_momentum_bar*ny

  # Reverse diff ny = nrm[2]*fac
  fac_bar = zero(Tsol)
  nrm_bar[2] += ny_bar*fac
  fac_bar += ny_bar*nrm[2]

  # Reverse diff nx = nrm[1]*fac
  nrm_bar[1] += nx_bar*fac
  fac_bar += nx_bar*nrm[1]

  # Reverse diff fac = 1.0/(sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2]))
  nrm_bar[1] -= fac_bar*((nrm[1]*nrm[1] + nrm[2]*nrm[2])^(-1.5))*nrm[1]
  nrm_bar[2] -= fac_bar*((nrm[1]*nrm[1] + nrm[2]*nrm[2])^(-1.5))*nrm[2]

  return nothing
end # End calcBoundaryFunctionalIntegrand_revm 2D

function calcBoundaryFunctionalIntegrand_revm(params::ParamType{3},
                                         q::AbstractArray{Tsol,1},
                                         aux_vars::AbstractArray{Tres, 1},
                                         nrm::AbstractArray{Tmsh},
                                         node_info::AbstractArray{Int},
                                         objective::BoundaryForceData,
                                         nrm_bar::AbstractArray{Tmsh,1},
                                         val_bar::Number) where {Tsol, Tres, Tmsh}

  # Forward Sweep
  fac = 1.0/(sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3]))
  nx = nrm[1]*fac
  ny = nrm[2]*fac
  nz = nrm[3]*fac

  normal_momentum = nx*q[2] + ny*q[3] + nz*q[4]
  qg = objective.qg
  for i=1:length(q)
    qg[i] = q[i]
  end
  qg[2] -= nx*normal_momentum
  qg[3] -= ny*normal_momentum
  qg[4] -= nz*normal_momentum

  # Reverse Sweep
  #TODO: preallocate these
  euler_flux_bar = zeros(Tsol, 5) # For 2D
  qg_bar = zeros(Tsol, 5)
  q_bar = zeros(Tsol,5)

  # Reverse diff val[:] = euler_flux[2:4]
  euler_flux_bar[2] = objective.facx*val_bar
  euler_flux_bar[3] = objective.facy*val_bar
  euler_flux_bar[4] = objective.facz*val_bar

  # Reverse diff calcEulerFlux
  calcEulerFlux_revm(params, qg, aux_vars, nrm, nrm_bar, euler_flux_bar)
  calcEulerFlux_revq(params, qg, qg_bar, aux_vars, nrm, euler_flux_bar)
  nz_bar = zero(Tsol)               #
  ny_bar = zero(Tsol)               # Initialize
  nx_bar = zero(Tsol)               #
  normal_momentum_bar = zero(Tsol)  #

  # qg[4] -= nz*normal_momentum
  nz_bar -= qg_bar[4]*normal_momentum
  normal_momentum_bar -= qg_bar[4]*nz
  qg_bar[4] += qg_bar[4]

  # Reverse diff qg[3] -= ny*normal_momentum
  ny_bar -= qg_bar[3]*normal_momentum
  normal_momentum_bar -= qg_bar[3]*ny
  qg_bar[3] += qg_bar[3]

  # Reverse diff qg[2] -= nx*normal_momentum
  nx_bar -= qg_bar[2]*normal_momentum
  normal_momentum_bar -= qg_bar[2]*nx
  qg_bar[2] += qg_bar[2]

  # Reverse diff qg[:] = q[:]a
  for i=1:length(q_bar)
    q_bar[i] += qg_bar[i]
  end

  # normal_momentum = nx*q[2] + ny*q[3] + nz*q[4]
  nx_bar += normal_momentum_bar*q[2]
  ny_bar += normal_momentum_bar*q[3]
  nz_bar += normal_momentum_bar*q[4]

  # nz = nrm[3]*fac
  nrm_bar[3] += nz_bar*fac
  fac_bar = nz_bar*nrm[3]

  # Reverse diff ny = nrm[2]*fac
  nrm_bar[2] += ny_bar*fac
  fac_bar += ny_bar*nrm[2]

  # Reverse diff nx = nrm[1]*fac
  nrm_bar[1] += nx_bar*fac
  fac_bar += nx_bar*nrm[1]

  # fac = 1.0/(sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3]))
  nrm_bar[1] -= fac_bar*((nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3])^(-1.5))*nrm[1]
  nrm_bar[2] -= fac_bar*((nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3])^(-1.5))*nrm[2]
  nrm_bar[3] -= fac_bar*((nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3])^(-1.5))*nrm[3]

  return nothing
end # End calcBoundaryFunctionalIntegrand_revm 3D


text:
# Specify contrasts for coding categorical data in model matrix. Contrasts types
# are a subtype of AbstractContrasts. ContrastsMatrix types hold a contrast
# matrix, levels, and term names and provide the interface for creating model
# matrix columns and coefficient names.
#
# Contrasts types themselves can be instantiated to provide containers for
# contrast settings (currently, just the base level).
#
# ModelFrame will hold a Dict{Symbol, ContrastsMatrix} that maps column
# names to contrasts.
#
# ModelMatrix will check this dict when evaluating terms, falling back to a
# default for any categorical data without a specified contrast.

"""
Interface to describe contrast coding schemes for categorical variables.

Concrete subtypes of `AbstractContrasts` describe a particular way of converting a
categorical data vector into numeric columns in a `ModelMatrix`. Each
instantiation optionally includes the levels to generate columns for and the base
level. If not specified these will be taken from the data when a `ContrastsMatrix` is
generated (during `ModelFrame` construction).

# Constructors

For `C <: AbstractContrast`:

```julia
C()                                     # levels are inferred later
C(levels = ::Vector{Any})               # levels checked against data later
C(base = ::Any)                         # specify base level
C(levels = ::Vector{Any}, base = ::Any) # specify levels and base
```

If specified, levels will be checked against data when generating a
`ContrastsMatrix`. Any mismatch will result in an error, because missing data
levels would lead to empty columns in the model matrix, and missing contrast
levels would lead to empty or undefined rows.

You can also specify the base level of the contrasts. The actual interpretation
of this depends on the particular contrast type, but in general it can be
thought of as a "reference" level.  It defaults to the first level.

# Concrete types

* `DummyCoding` - Code each non-base level as a 0-1 indicator column.
* `EffectsCoding` - Code each non-base level as 1, and base as -1.
* `HelmertCoding` - Code each non-base level as the difference from the mean of
  the lower levels
* `ContrastsCoding` - Manually specify contrasts matrix

The last coding type, `ContrastsCoding`, provides a way to manually specify a
contrasts matrix. For a variable `x` with k levels, a contrasts matrix `M` is a
k by k-1 matrix, that maps the k levels onto k-1 model matrix columns.
Specifically, let ``X^*`` be the full-rank indicator matrix for `x`, where
``X^*_{i,j} = 1`` if `x[i]` is level `j`, and 0 otherwise. Then the model matrix
columns generated by the contrasts matrix `M` are ``X = X^* M``.

To implement your own `AbstractContrasts` type, implement a constructor, a
`contrasts_matrix` method for constructing the actual contrasts matrix that maps
from levels to `ModelMatrix` column values, and (optionally) a `termnames`
method:

```julia
type MyCoding <: AbstractContrasts
    ...
end

contrasts_matrix(C::MyCoding, baseind, n) = ...
termnames(C::MyCoding, levels, baseind) = ...
```

"""
@compat abstract type AbstractContrasts end

# Contrasts + Levels (usually from data) = ContrastsMatrix
type ContrastsMatrix{C <: AbstractContrasts, T}
    matrix::Matrix{Float64}
    termnames::Vector{T}
    levels::Vector{T}
    contrasts::C
end

"""
    ContrastsMatrix{C<:AbstractContrasts}(contrasts::C, levels::AbstractVector)

Compute contrasts matrix for given data levels.

If levels are specified in the `AbstractContrasts`, those will be used, and likewise
for the base level (which defaults to the first level).
"""
function ContrastsMatrix{C <: AbstractContrasts}(contrasts::C, levels::AbstractVector)

    # if levels are defined on contrasts, use those, validating that they line up.
    # what does that mean? either:
    #
    # 1. contrasts.levels == levels (best case)
    # 2. data levels missing from contrast: would generate empty/undefined rows.
    #    better to filter data frame first
    # 3. contrast levels missing from data: would have empty columns, generate a
    #    rank-deficient model matrix.
    c_levels = get(contrasts.levels, levels)
    if eltype(c_levels) != eltype(levels)
        throw(ArgumentError("mismatching levels types: got $(eltype(levels)), expected " *
                            "$(eltype(c_levels)) based on contrasts levels."))
    end
    mismatched_levels = symdiff(c_levels, levels)
    if !isempty(mismatched_levels)
        throw(ArgumentError("contrasts levels not found in data or vice-versa: " *
                            "$mismatched_levels." *
                            "\n  Data levels: $levels." *
                            "\n  Contrast levels: $c_levels"))
    end

    n = length(c_levels)
    if n == 0
        throw(ArgumentError("empty set of levels found (need at least two to compute " *
                            "contrasts)."))
    elseif n == 1
        throw(ArgumentError("only one level found: $(c_levels[1]) (need at least two to " *
                            "compute contrasts)."))
    end

    # find index of base level. use contrasts.base, then default (1).
    baseind = isnull(contrasts.base) ?
              1 :
              findfirst(c_levels, get(contrasts.base))
    if baseind < 1
        throw(ArgumentError("base level $(get(contrasts.base)) not found in levels " *
                            "$c_levels."))
    end

    tnames = termnames(contrasts, c_levels, baseind)

    mat = contrasts_matrix(contrasts, baseind, n)

    ContrastsMatrix(mat, tnames, c_levels, contrasts)
end

# Methods for constructing ContrastsMatrix from data. These are called in
# ModelFrame constructor and setcontrasts!.
# TODO: add methods for new categorical types

ContrastsMatrix(C::AbstractContrasts, v::PooledDataArray) =
    ContrastsMatrix(C, levels(v))
ContrastsMatrix{C <: AbstractContrasts}(c::Type{C}, col::PooledDataArray) =
    throw(ArgumentError("contrast types must be instantiated (use $c() instead of $c)"))
# given an existing ContrastsMatrix, check that all of the levels present in the
# data are present in the contrasts. Note that this behavior is different from the
# ContrastsMatrix constructor, which requires that the levels be exactly the same.
# This method exists to support things like `predict` that can operate on new data
# which may contain only a subset of the original data's levels. Checking here
# (instead of in `modelmat_cols`) allows an informative error message.
function ContrastsMatrix(c::ContrastsMatrix, col::PooledDataArray)
    if !isempty(setdiff(levels(col), c.levels))
        throw(ArgumentError("there are levels in data that are not in ContrastsMatrix: " *
                            "$(setdiff(levels(col), c.levels))" *
                            "\n  Data levels: $(levels(col))" *
                            "\n  Contrast levels: $(c.levels)"))
    end
    return c
end

function termnames(C::AbstractContrasts, levels::AbstractVector, baseind::Integer)
    not_base = [1:(baseind-1); (baseind+1):length(levels)]
    levels[not_base]
end

nullify(x::Nullable) = x
nullify(x) = Nullable(x)

# Making a contrast type T only requires that there be a method for
# contrasts_matrix(T, v::PooledDataArray). The rest is boilerplate.
for contrastType in [:DummyCoding, :EffectsCoding, :HelmertCoding]
    @eval begin
        type $contrastType <: AbstractContrasts
            base::Nullable{Any}
            levels::Nullable{Vector}
        end
        ## constructor with optional keyword arguments, defaulting to Nullables
        $contrastType(;
                      base=Nullable{Any}(),
                      levels=Nullable{Vector}()) =
                          $contrastType(nullify(base),
                                        nullify(levels))
    end
end

"""
    FullDummyCoding()

Coding that generates one indicator (1 or 0) column for each level,
__including__ the base level.

Needed internally when a term is non-redundant with lower-order terms (e.g., in
`~0+x` vs. `~1+x`, or in the interactions terms in `~1+x+x&y` vs. `~1+x+y+x&y`. In the
non-redundant cases, we can expand x into `length(levels(x))` columns
without creating a non-identifiable model matrix (unless the user has done
something foolish in specifying the model, which we can't do much about anyway).
"""
type FullDummyCoding <: AbstractContrasts
# Dummy contrasts have no base level (since all levels produce a column)
end

ContrastsMatrix{T}(C::FullDummyCoding, lvls::Vector{T}) =
    ContrastsMatrix(eye(Float64, length(lvls)), lvls, lvls, C)

"Promote contrasts matrix to full rank version"
Base.convert(::Type{ContrastsMatrix{FullDummyCoding}}, C::ContrastsMatrix) =
    ContrastsMatrix(FullDummyCoding(), C.levels)

"""
    DummyCoding([base[, levels]])

Contrast coding that generates one indicator column (1 or 0) for each non-base level.

Columns have non-zero mean and are collinear with an intercept column (and
lower-order columns for interactions) but are orthogonal to each other. In a
regression model, dummy coding leads to an intercept that is the mean of the
dependent variable for base level.

Also known as "treatment coding" (`contr.treatment` in R) or "one-hot encoding".
"""
DummyCoding

contrasts_matrix(C::DummyCoding, baseind, n) = eye(n)[:, [1:(baseind-1); (baseind+1):n]]

"""
    EffectsCoding([base[, levels]])

Contrast coding that generates columns that code each non-base level as the
deviation from the base level.  For each non-base level `x` of `variable`, a
column is generated with 1 where `variable .== x` and -1 where `col .== base`.

`EffectsCoding` is like `DummyCoding`, but using -1 for the base level instead
of 0.

When all levels are equally frequent, effects coding generates model matrix
columns that are mean centered (have mean 0).  For more than two levels the
generated columns are not orthogonal.  In a regression model with an
effects-coded variable, the intercept corresponds to the grand mean.

Also known as "sum coding" (`contr.sum` in R) or "simple coding" (SPSS). Note
though that the default in R and SPSS is to use the _last_ level as the base.
Here we use the _first_ level as the base, for consistency with other coding
schemes.
"""
EffectsCoding

function contrasts_matrix(C::EffectsCoding, baseind, n)
    not_base = [1:(baseind-1); (baseind+1):n]
    mat = eye(n)[:, not_base]
    mat[baseind, :] = -1
    return mat
end

"""
    HelmertCoding([base[, levels]])

Contrasts that code each level as the difference from the average of the lower
levels.

For each non-base level, Helmert coding generates a columns with -1 for each of
n levels below, n for that level, and 0 above.

# Examples

```julia
julia> ContrastsMatrix(HelmertCoding(), collect(1:4)).matrix
4x3 Array{Float64,2}:
 -1.0  -1.0  -1.0
  1.0  -1.0  -1.0
  0.0   2.0  -1.0
  0.0   0.0   3.0
```

When all levels are equally frequent, Helmert coding generates columns that are
mean-centered (mean 0) and orthogonal.
"""
HelmertCoding

function contrasts_matrix(C::HelmertCoding, baseind, n)
    mat = zeros(n, n-1)
    for i in 1:n-1
        mat[1:i, i] = -1
        mat[i+1, i] = i
    end

    # re-shuffle the rows such that base is the all -1.0 row (currently first)
    mat = mat[[baseind; 1:(baseind-1); (baseind+1):end], :]
    return mat
end

"""
    ContrastsCoding(mat::Matrix[, base[, levels]])

Coding by manual specification of contrasts matrix. For k levels, the contrasts
must be a k by k-1 Matrix.
"""
type ContrastsCoding <: AbstractContrasts
    mat::Matrix
    base::Nullable{Any}
    levels::Nullable{Vector}

    function ContrastsCoding(mat, base, levels)
        if !isnull(levels)
            check_contrasts_size(mat, length(get(levels)))
        end
        new(mat, base, levels)
    end
end

check_contrasts_size(mat::Matrix, n_lev) =
    size(mat) == (n_lev, n_lev-1) ||
    throw(ArgumentError("contrasts matrix wrong size for $n_lev levels. " *
                        "Expected $((n_lev, n_lev-1)), got $(size(mat))"))

## constructor with optional keyword arguments, defaulting to Nullables
ContrastsCoding(mat::Matrix; base=Nullable{Any}(), levels=Nullable{Vector}()) =
    ContrastsCoding(mat, nullify(base), nullify(levels))

function contrasts_matrix(C::ContrastsCoding, baseind, n)
    check_contrasts_size(C.mat, n)
    C.mat
end


text:
## Real-to-complex and complex-to-real transforms.

"""
    RFFT()

Real-to-complex FFT.

See also
[`AbstractFFTs.rfft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.rfft).
"""
struct RFFT <: AbstractTransform end

"""
    BRFFT(d::Integer)
    BRFFT((d1, d2, ..., dN))

Unnormalised inverse of [`RFFT`](@ref).

To obtain the inverse transform, divide the output by the length of the
transformed dimension (of the real output array).

As described in the [AbstractFFTs docs](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.irfft),
the length of the output cannot be fully inferred from the input length.
For this reason, the `BRFFT` constructor accepts an optional `d` argument
indicating the output length.

For multidimensional datasets, a tuple of dimensions
`(d1, d2, ..., dN)` may also be passed.
This is equivalent to passing just `dN`.
In this case, the **last** dimension (`dN`) is the one that changes size between
the input and output.
Note that this is the opposite of `FFTW.brfft`.
The reason is that, in PencilFFTs, the **last** dimension is the one along which
a complex-to-real transform is performed.

See also
[`AbstractFFTs.brfft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.brfft).
"""
struct BRFFT <: AbstractTransform
    even_output :: Bool
end

_show_extra_info(io::IO, tr::BRFFT) = print(io, tr.even_output ? "{even}" : "{odd}")

BRFFT(d::Integer) = BRFFT(iseven(d))
BRFFT(ts::Tuple) = BRFFT(last(ts))  # c2r transform is applied along the **last** dimension (opposite of FFTW)

is_inplace(::Union{RFFT, BRFFT}) = false

length_output(::RFFT, length_in::Integer) = div(length_in, 2) + 1
length_output(tr::BRFFT, length_in::Integer) = 2 * length_in - 1 - tr.even_output

eltype_output(::RFFT, ::Type{T}) where {T <: FFTReal} = Complex{T}
eltype_output(::BRFFT, ::Type{Complex{T}}) where {T <: FFTReal} = T

eltype_input(::RFFT, ::Type{T}) where {T <: FFTReal} = T
eltype_input(::BRFFT, ::Type{T}) where {T <: FFTReal} = Complex{T}

plan(::RFFT, args...; kwargs...) = FFTW.plan_rfft(args...; kwargs...)

# NOTE: unlike most FFTW plans, this function also requires the length `d` of
# the transform output along the first transformed dimension.
function plan(tr::BRFFT, A, dims; kwargs...)
    Nin = size(A, first(dims))  # input length along first dimension
    d = length_output(tr, Nin)
    FFTW.plan_brfft(A, d, dims; kwargs...)
end

binv(::RFFT, d) = BRFFT(d)
binv(::BRFFT, d) = RFFT()

function scale_factor(tr::BRFFT, A::ComplexArray, dims)
    prod(dims; init = one(Int)) do i
        n = size(A, i)
        i == last(dims) ? length_output(tr, n) : n
    end
end

scale_factor(::RFFT, A::RealArray, dims) = _prod_dims(A, dims)

# r2c along the first dimension, then c2c for the other dimensions.
expand_dims(tr::RFFT, ::Val{N}) where {N} =
    N === 0 ? () : (tr, expand_dims(FFT(), Val(N - 1))...)

expand_dims(tr::BRFFT, ::Val{N}) where {N} = (BFFT(), expand_dims(tr, Val(N - 1))...)
expand_dims(tr::BRFFT, ::Val{1}) = (tr, )
expand_dims(tr::BRFFT, ::Val{0}) = ()


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
np = pyimport("numpy")

stepsize = 1
if length(ARGS)==1
  global stepsize = parse(Int64, ARGS[1])
end
@info stepsize

## alpha-scheme
β = 1/4; γ = 1/2
a = b = 0.1    # damping

mmesh = Mesh(10, 10, 1/10, degree=1)
NT = 500
Δt = 2.0/NT
ηmax = 1
ηmin = 0.5

bdedge = bcedge((x1, y1, x2, y2)->(x1>1-1e-5 && x2>1-1e-5), mmesh)
bdnode = bcnode((x1, y1)->(y1<1e-5), mmesh)

invη = 50 * constant(ones(get_ngauss(mmesh)))
μ = 1.0
λ = 1.0

fn_G = invη->begin
  G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
    -μ/3*invη 1/Δt+2/3*μ*invη 0.0
    0.0 0.0 1/Δt+μ*invη])
  invG = inv(G)
end
invG = map(fn_G, invη)
S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
    λ/Δt 2μ/Δt+λ/Δt 0.0
    0.0 0.0 μ/Δt])
H = invG*S

M = compute_fem_mass_matrix1(mmesh)
Zero = spzero(mmesh.ndof, mmesh.ndof)
M = [M Zero;Zero M]

K = compute_fem_stiffness_matrix(H, mmesh)
C = a*M + b*K # damping matrix
L = M + γ*Δt*C + β*Δt^2*K
bddof = [bdnode; bdnode .+ mmesh.ndof]
L, _ = impose_Dirichlet_boundary_conditions(L, zeros(2mmesh.ndof), bddof, zeros(length(bddof)))
L = factorize(L)

a = TensorArray(NT+1); a = write(a, 1, zeros(2mmesh.ndof))
v = TensorArray(NT+1); v = write(v, 1, zeros(2mmesh.ndof))
d = TensorArray(NT+1); d = write(d, 1, zeros(2mmesh.ndof))
U = TensorArray(NT+1); U = write(U, 1, zeros(2mmesh.ndof))
Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(get_ngauss(mmesh), 3))
Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1, zeros(get_ngauss(mmesh), 3))

Forces = zeros(NT, 2mmesh.ndof)
T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, mmesh)
T = [-T T]
rhs = compute_fem_traction_term(T, bdedge, mmesh)
for i = 1:NT
  Forces[i, :] = rhs
end
Forces = constant(Forces)

function condition(i, tas...)
  i <= NT
end

function body(i, tas...)
  a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
  a = read(a_, i)
  v = read(v_, i)
  d = read(d_, i)
  U = read(U_, i)
  Sigma = read(Sigma_, i)
  Varepsilon = read(Varepsilon_, i)

  res = batch_matmul(invG/Δt, Sigma)
  F = compute_strain_energy_term(res, mmesh) - K * U
  rhs = Forces[i] - F

  td = d + Δt*v + Δt^2/2*(1-2β)*a
  tv = v + (1-γ)*Δt*a
  rhs = rhs - C*tv - K*td
  rhs = scatter_update(rhs, bddof, zeros(length(bddof)))

  ## alpha-scheme
  a = L\rhs # bottleneck
  d = td + β*Δt^2*a
  v = tv + γ*Δt*a
  U_new = d

  Varepsilon_new = eval_strain_on_gauss_pts(U_new, mmesh)
  Sigma_new = update_stress_viscosity(Varepsilon_new, Varepsilon, Sigma, invη, μ*ones(get_ngauss(mmesh)), λ*ones(get_ngauss(mmesh)), Δt)

  i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
        write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
end

i = constant(1, dtype=Int32)
_, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                  [i, a, v, d, U, Sigma, Varepsilon])

U = stack(u)
Sigma = stack(sigma)
Varepsilon = stack(varepsilon)

sess = Session(); init(sess)
Uval,Sigmaval, Varepsilonval = run(sess, [U, Sigma, Varepsilon])

# p = visualize_von_mises_stress(Sigmaval[1:5:end,:,:], m, n, h); saveanim(p, "space_s.gif")
p = visualize_displacement(Uval[1:5:end,:], mmesh)
saveanim(p, "space_u.gif")

# visualize_von_mises_stress(Sigmaval[end,:,:], m, n, h); savefig("space_s.pdf")
# visualize_displacement(Uval[end,:], m, n, h); savefig("space_u.pdf")

# visualize_inv_eta(run(sess, invη), "true")


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct SoundAnalysisPreprocessing_Vggish <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SoundAnalysisPreprocessing_Vggish(; kwargs...)
        obj = new(meta(SoundAnalysisPreprocessing_Vggish), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SoundAnalysisPreprocessing_Vggish
const __meta_SoundAnalysisPreprocessing_Vggish = Ref{ProtoMeta}()
function meta(::Type{SoundAnalysisPreprocessing_Vggish})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SoundAnalysisPreprocessing_Vggish)
            __meta_SoundAnalysisPreprocessing_Vggish[] = target = ProtoMeta(SoundAnalysisPreprocessing_Vggish)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, SoundAnalysisPreprocessing_Vggish, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SoundAnalysisPreprocessing_Vggish[]
    end
end

mutable struct SoundAnalysisPreprocessing <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SoundAnalysisPreprocessing(; kwargs...)
        obj = new(meta(SoundAnalysisPreprocessing), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SoundAnalysisPreprocessing
const __meta_SoundAnalysisPreprocessing = Ref{ProtoMeta}()
function meta(::Type{SoundAnalysisPreprocessing})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SoundAnalysisPreprocessing)
            __meta_SoundAnalysisPreprocessing[] = target = ProtoMeta(SoundAnalysisPreprocessing)
            fnum = Int[20]
            allflds = Pair{Symbol,Union{Type,String}}[:vggish => SoundAnalysisPreprocessing_Vggish]
            oneofs = Int[1]
            oneof_names = Symbol[Symbol("SoundAnalysisPreprocessingType")]
            meta(target, SoundAnalysisPreprocessing, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_SoundAnalysisPreprocessing[]
    end
end
function Base.getproperty(obj::SoundAnalysisPreprocessing, name::Symbol)
    if name === :vggish
        return (obj.__protobuf_jl_internal_values[name])::SoundAnalysisPreprocessing_Vggish
    else
        getfield(obj, name)
    end
end

export SoundAnalysisPreprocessing_Vggish, SoundAnalysisPreprocessing


text:
module Clustering

    include("knn_clustering.jl")

    export knn_clustering

end # module

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO


text:
function Base.read(fpath::String, ::Type{SearchTrace}; accept::Symbol = :accept)
    json_string = read(fpath, String)
    stfile = JSON3.read(json_string)

    df = DataFrame(jsontable(stfile[:trace]))
    df[!, accept] = Bool.(df[:, accept])

    rep = nothing
    if :s in keys(stfile[:initial])
        rep = stfile[:initial][:s]
    end

    st = SearchTrace(
        rep,
        stfile[:initial][:v],
        df
    )
end

function Base.write(fpath::String, st::SearchTrace; accept::Symbol= :accept)
    open(fpath, "w") do io
        write(io, serialise(st, accept))
    end
end

function serialise(st::SearchTrace, accept::Symbol)
    st.df[!, accept] = Int8.(st.df[:, accept])
    str = string(
        "{\n" ,
            "\"initial\": { \n\t\"s\": " , serialise(st.s_init) ,
                                ", \n\t\t\"v\": " , serialise(st.v_init) , "}" ,
            ", \n\"trace\": " , serialise(st.df) ,
        "\n}"
    )
    return str
end

function serialise(st::SearchTrace{Nothing, V}, accept) where {V}
    st.df[!, accept] = Int8.(st.df[:, accept])
    str = string(
        "{\n" ,
            "\"initial\": { \n\t\"v\": " , serialise(st.v_init) , "}" ,
            ", \n\"trace\": " , serialise(st.df) ,
        "\n}"
    )
    return str
end

function serialise(df::DataFrame)
    df.op = string.(df.op)
    return arraytable(df)
end

function serialise(t::Any)
    return JSON3.write(t)
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz

max_forks_repo_name: Michiel-VL/PolyaViz


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

@testset "zhang ip" begin
	numStations = 2
	stationSlots = [1, 1, 2, 2] # each station can hold two ambulances
	benefitSlots = [4.0, 1.0, 3.0, 2.0] # benefitSlots[i] gives value of placing ambulance at stationSlots[i]
	@assert(all(j -> issorted(benefitSlots[stationSlots .== j], lt=<=, rev=true), numStations)) # value of placing kth amb at station j > value of k+1th amb

	function solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, numMovableAmbs, numAtHospitalAmbs)
		if typeof(ambToStationCosts) <: Vector
			ambToStationCosts = convert(Array{Float,2}, ambToStationCosts')
		end
		solution1 = JEMSS.solveZhangIp(stationSlots, benefitSlots, ambToStationCosts, numMovableAmbs, numAtHospitalAmbs, true)
		solution2 = JEMSS.solveZhangIpAssignmentProblem(stationSlots, benefitSlots, ambToStationCosts, numMovableAmbs, numAtHospitalAmbs)
		return solution1 == solution2 ? solution1 : (solution1, solution2)
	end

	# case: no ambs
	ambToStationCosts = zeros(Float, 0, numStations) # ambToStationCosts[i,j] is cost of assigning amb i to station j (excludes benefit)
	expectedSolution = (Int[], Int[])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 0, 0)

	## 1 amb cases

	# case: 1 amb at station (not moved)
	ambToStationCosts = Float[1.5 0] # amb at station 2; cost to move is greater than benefit increase
	expectedSolution = ([2], Int[])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 1, 0)

	# case: 1 amb at station (moved)
	ambToStationCosts = Float[0.5 0] # amb at station 2; cost to move is less than benefit increase
	expectedSolution = ([1], Int[]) # moved to station 1
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 1, 0)

	# case: 1 amb at station (assigned to station even with negative objective value)
	ambToStationCosts = Float[5 5]
	expectedSolution = ([1], Int[])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 1, 0)

	# case: one amb at hospital (assigned to station 1)
	ambToStationCosts = Float[1.5 1] # value for station 1: 4 - 1.5 = 3.5; for station 2: 3 - 1 = 2
	expectedSolution = (Int[], [1])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 0, 1)

	# case: one amb at hospital (not assigned to a station)
	ambToStationCosts = Float[5 5] # costs outweight benefits, will not redeploy amb
	expectedSolution = (Int[], [nullIndex])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 0, 1)

	## 2 amb cases

	# case: 2 ambs at stations (not moved)
	ambToStationCosts = Float[0 1; 1 0]
	expectedSolution = ([1,2], Int[])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 2, 0)

	# case: 2 ambs at 1 station (not moved)
	ambToStationCosts = Float[0 4; 0 4]
	expectedSolution = ([1,1], Int[])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 2, 0)

	# case: 2 ambs at 1 station (amb 2 moved)
	ambToStationCosts = Float[0 2; 0 1]
	expectedSolution = ([1,2], Int[])
	@test expectedSolution == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 2, 0)

	# case: 1 amb at station, 1 at hospital; amb at hospital keeps amb from moving to better station
	ambToStationCosts = Float[0.5 0; 1 1] # amb 1 would move to station 1 if second amb wasn't at hospital
	# note: ambToStationCosts[1,1] + ambToStationCosts[2,2] > ambToStationCosts[1,2] + ambToStationCosts[2,1]
	@test ([1], Int[]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts[1,:], 1, 0)
	@test ([2], [1]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 1, 1)

	# case: 1 amb at station, 1 at hospital; amb at hospital causes amb to move to worse station
	ambToStationCosts = Float[0 0.5; 1 2]
	# note: ambToStationCosts[1,1] + ambToStationCosts[2,2] > ambToStationCosts[1,2] + ambToStationCosts[2,1]
	@test ([1], Int[]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts[1,:], 1, 0) # without amb at hospital
	@test ([2], [1]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 1, 1)

	## 4 amb cases

	# case: 4 ambs at stations
	ambToStationCosts = Float[0 1; 0 1; 1 0; 1 0]
	@test ([1,1,2,2], Int[]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 4, 0)

	# case: 3 ambs at stations, 1 at hospital; amb at hospital causes amb 3 to move from station 2 to 1
	ambToStationCosts = Float[0 2; 1 0; 0.5 0; 1 0.1]
	# note: ambToStationCosts[3,1] + ambToStationCosts[4,2] < ambToStationCosts[3,2] + ambToStationCosts[4,1]
	@test ([1,2,2], Int[]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts[1:3,:], 3, 0) # without amb at hospital
	@test ([1,2,1], [2]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 3, 1) # with amb at hospital

	# case: 2 ambs at station 1, 2 at hospital; ambs at hospital cause ambs at station 1 to move to station 2
	ambToStationCosts = Float[0 2.2; 0 2.1; -2 4; -2 4] # ambs at station 1 would not move if not for ambs at hospitals
	@test ([1,1], Int[]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts[1:2,:], 2, 0) # no ambs at hospital
	@test ([1,2], [1]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts[1:3,:], 2, 1) # 1 amb at hospital
	@test ([2,2], [1,1]) == solveZhangIpTest(stationSlots, benefitSlots, ambToStationCosts, 2, 2) # 2 ambs at hospital

end


text:
using Documenter, SimpleDifferentialOperators

# figure out what to paste instead of TAG_GOES_HERE
replacement = get(ENV, "DOCS_TAG", "dev")

tmp = @__DIR__
fileset = readdir(tmp * "/src/")
for file in fileset
    if occursin(".md", file)
        txt = read(tmp * "/src/" * file, String)
        open(tmp * "/src/" * file, "w") do f
           write(f, replace(txt, "TAG_GOES_HERE" => replacement))
        end
    end
end

# Compile the online documentation.
makedocs(sitename = "SimpleDifferentialOperators.jl",
	pages = [
        "Home" => "index.md",
        "Gallery" => [
            "Examples" => "examples.md",
            "Notebooks" => "notebooks.md"
        ],
        "api.md"
    ],
	doctest = :fix)

# Deploy to git
deploydocs(
    repo = "github.com/QuantEcon/SimpleDifferentialOperators.jl.git",
)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random, Serialization

# Check that serializer hasn't gone out-of-frame
@test Serialization.sertag(Symbol) == 1
@test Serialization.sertag(()) == 68
@test Serialization.sertag(false) == 76

function create_serialization_stream(f::Function)
    s = IOBuffer()
    f(s)
    close(s)
end

# Tags
create_serialization_stream() do s
    Serialization.writetag(s, Serialization.sertag(Bool))
    @test take!(s)[end] == UInt8(Serialization.sertag(Bool))
end

create_serialization_stream() do s
    Serialization.write_as_tag(s, Serialization.sertag(Bool))
    @test take!(s)[end] == UInt8(Serialization.sertag(Bool))
end

create_serialization_stream() do s
    Serialization.write_as_tag(s, Serialization.sertag(Symbol))
    data = take!(s)
    @test data[end-1] == 0x00
    @test data[end] == UInt8(Serialization.sertag(Symbol))
end

# SubString

create_serialization_stream() do s
    ss = Any[SubString("12345", 2, 4),
             SubString(GenericString("12345"), 2, 4),
             SubString(s"12345", 2, 4),
             SubString(GenericString(s"12345"), 2, 4),]
    for x in ss
        serialize(s, x)
    end
    seek(s, 0)
    for x in ss
        y = deserialize(s)
        @test x == y
        @test typeof(x) == typeof(y)
    end
end

# Boolean & Empty & Nothing
create_serialization_stream() do s
    serialize(s, true)
    serialize(s, ())
    serialize(s, nothing)
    seek(s, 0)
    @test deserialize(s) === true
    @test deserialize(s) === ()
    @test deserialize(s) === nothing
end

# Ptr
create_serialization_stream() do s
    serialize(s, C_NULL)
    seekstart(s)
    @test deserialize(s) === C_NULL
end

# Integer
create_serialization_stream() do s
    serialize(s, 0x1)
    serialize(s, 32)
    serialize(s, 33)
    seek(s,0)
    @test deserialize(s) === 0x01
    @test deserialize(s) === 32
    @test deserialize(s) === 33
end

# Tuple
create_serialization_stream() do s
    tpl = (0x1,'a')
    serialize(s, tpl)

    len = 257
    lt = ntuple(i->0x1, len)
    serialize(s, lt)

    serialize(s, Tuple{})

    seek(s, 0)
    @test deserialize(s) === tpl
    @test deserialize(s) === lt
    @test deserialize(s) === Tuple{}
end

# Dict
create_serialization_stream() do s
    dct = Dict("A"=>1, "B"=>2)
    serialize(s, dct)
    seek(s, 0)
    @test deserialize(s) == dct
end

# Symbol
create_serialization_stream() do s
    gensym(len) = Symbol(repeat("A", len))

    smbl = gensym(1)
    serialize(s, smbl)

    smbl2 = gensym(10)
    serialize(s, smbl2)

    smbl3 = gensym(257)
    serialize(s, smbl3)

    seek(s, 0)
    @test deserialize(s) === smbl
    @test deserialize(s) === smbl2
    @test deserialize(s) === smbl3
end

# Module
create_serialization_stream() do s # user-defined module
    mod = b"SomeModule"
    modstring = String(mod)
    eval(Meta.parse("module $(modstring); end"))
    modtype = eval(Meta.parse("$(modstring)"))
    serialize(s, modtype)
    seek(s, 0)
    @test deserialize(s) === modtype
end

# DataType
create_serialization_stream() do s # standard types
    typ1 = UInt8
    serialize(s, typ1)
    typ2 = Bool
    serialize(s, typ2)

    seek(s, 0)
    @test deserialize(s) == typ1
    @test deserialize(s) == typ2
end

create_serialization_stream() do s # user-defined type
    usertype = "SerializeSomeType"
    eval(Meta.parse("abstract type $(usertype) end"))
    utype = eval(Meta.parse("$(usertype)"))
    serialize(s, utype)
    seek(s, 0)
    @test deserialize(s) === utype
end

create_serialization_stream() do s # user-defined type
    usertype = "SerializeSomeType1"
    eval(Meta.parse("mutable struct $(usertype); end"))
    utype = eval(Meta.parse("$(usertype)"))
    serialize(s, utype)
    seek(s, 0)
    @test deserialize(s) === utype
end

create_serialization_stream() do s # user-defined type
    usertype = "SerializeSomeType2"
    eval(Meta.parse("abstract type $(usertype){T} end"))
    utype = eval(Meta.parse("$(usertype)"))
    serialize(s, utype)
    seek(s, 0)
    @test deserialize(s) == utype
end

create_serialization_stream() do s # immutable struct with 1 field
    usertype = "SerializeSomeType3"
    eval(Meta.parse("struct $(usertype){T}; a::T; end"))
    utype = eval(Meta.parse("$(usertype)"))
    serialize(s, utype)
    seek(s, 0)
    @test deserialize(s) == utype
end

create_serialization_stream() do s # immutable struct with 2 field
    usertype = "SerializeSomeType4"
    eval(Meta.parse("struct $(usertype){T}; a::T; b::T; end"))
    utval = eval(Meta.parse("$(usertype)(1,2)"))
    serialize(s, utval)
    seek(s, 0)
    @test deserialize(s) === utval
end

create_serialization_stream() do s # immutable struct with 3 field
    usertype = "SerializeSomeType5"
    eval(Meta.parse("struct $(usertype){T}; a::T; b::T; c::T; end"))
    utval = eval(Meta.parse("$(usertype)(1,2,3)"))
    serialize(s, utval)
    seek(s, 0)
    @test deserialize(s) === utval
end

create_serialization_stream() do s # immutable struct with 4 field
    usertype = "SerializeSomeType6"
    eval(Meta.parse("struct $(usertype){T}; a::T; b::T; c::T; d::T; end"))
    utval = eval(Meta.parse("$(usertype)(1,2,3,4)"))
    serialize(s, utval)
    seek(s, 0)
    @test deserialize(s) === utval
end

# Expression
create_serialization_stream() do s
    expr = Meta.parse("a = 1")
    serialize(s, expr)

    expr2 = Meta.parse(repeat("a = 1;", 300))
    serialize(s, expr2)

    seek(s, 0)
    @test deserialize(s) == expr
    @test deserialize(s) == expr2
end

# Array
mutable struct TA1
    v::UInt8
end
create_serialization_stream() do s # small 1d array
    arr1 = fill(0x01, 10)    # small 1d array
    serialize(s, arr1)
    arr2 = fill(0x01, 3, 3)  # small 2d array
    serialize(s, arr2)
    arr3 = fill(0x0001, 50)  # large 1d array
    serialize(s, arr3)
    arr4 = reshape([true, false, false, false, true, false, false, false, true], 3, 3)
    serialize(s, arr4)       # boolean array

    arr5 = Vector{TA1}(uninitialized, 3)
    arr5[2] = TA1(0x01)
    serialize(s, arr5)

    seek(s, 0)
    @test deserialize(s) == arr1
    @test deserialize(s) == arr2
    @test deserialize(s) == arr3
    @test deserialize(s) == arr4

    result = deserialize(s)
    @test !isassigned(result,1)
    @test !isassigned(result,3)
    @test result[2].v == arr5[2].v
end

# SubArray
create_serialization_stream() do s # slices
    slc1 = view(UInt8[1,1,1,1], 2:3)
    serialize(s, slc1)
    slc2 = view(repeat(UInt8[1,2,3,4], 1, 4), 1, 2:4)
    serialize(s, slc2)

    seek(s, 0)
    @test deserialize(s) == slc1
    @test deserialize(s) == slc2
end

# Objects that have a SubArray as a type in a type-parameter list
module ArrayWrappers

struct ArrayWrapper{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    data::A
end
ArrayWrapper(data::AbstractArray{T,N}) where {T,N} = ArrayWrapper{T,N,typeof(data)}(data)
Base.size(A::ArrayWrapper) = size(A.data)
Base.size(A::ArrayWrapper, d) = size(A.data, d)
Base.getindex(A::ArrayWrapper, i::Real...) = getindex(A.data, i...)

end

let A = rand(3,4)
    for B in (view(A, :, 2:4), view(A, 2, 1:3))
        C = ArrayWrappers.ArrayWrapper(B)
        io = IOBuffer()
        serialize(io, C)
        seek(io, 0)
        Cd = deserialize(io)
        @test size(Cd) == size(C)
        @test Cd == B
    end
end

# Function
serialize_test_function() = 1
serialize_test_function2 = ()->1
create_serialization_stream() do s # Base generic function
    serialize(s, sin)
    serialize(s, typeof)
    serialize(s, serialize_test_function)
    serialize(s, serialize_test_function2)

    seek(s, 0)
    @test deserialize(s) === sin
    @test deserialize(s) === typeof
    @test deserialize(s)() === 1
    @test deserialize(s)() === 1
end

# Anonymous Functions
main_ex = quote
    using Serialization
    $create_serialization_stream() do s
        local g() = :magic_token_anon_fun_test
        serialize(s, g)
        serialize(s, g)

        seekstart(s)
        ds = Serializer(s)
        local g2 = deserialize(ds)
        $Test.@test g2 !== g
        $Test.@test g2() == :magic_token_anon_fun_test
        $Test.@test g2() == :magic_token_anon_fun_test
        $Test.@test deserialize(ds) === g2

        # issue #21793
        y = x -> (() -> x)
        seekstart(s)
        serialize(s, y)
        seekstart(s)
        y2 = deserialize(s)
        x2 = y2(2)
        $Test.@test x2() == 2
    end
end
# This needs to be run on `Main` since the serializer treats it differently.
eval(Main, main_ex)

# Task
create_serialization_stream() do s # user-defined type array
    f = () -> begin task_local_storage(:v, 2); return 1+1 end
    t = Task(f)
    Base._wait(schedule(t))
    serialize(s, t)
    seek(s, 0)
    r = deserialize(s)
    @test r.storage[:v] == 2
    @test r.state == :done
    @test r.exception === nothing
end

struct MyErrorTypeTest <: Exception end
create_serialization_stream() do s # user-defined type array
    t = Task(()->throw(MyErrorTypeTest()))
    @test_throws MyErrorTypeTest Base._wait(schedule(t))
    serialize(s, t)
    seek(s, 0)
    r = deserialize(s)
    @test r.state == :failed
    @test isa(t.exception, MyErrorTypeTest)
end

# corner case: undefined inside immutable struct
struct MyNullable{T}
    hasvalue::Bool
    value::T

    MyNullable{T}() where {T} = new(false)
end
create_serialization_stream() do s
    serialize(s, MyNullable{Any}())
    seekstart(s)
    n = deserialize(s)
    @test isa(n, MyNullable)
    @test !isdefined(n, :value)
end

# cycles
module CycleFoo
    echo(x)=x
end
create_serialization_stream() do s
    echo(x) = x
    afunc = (x)->x
    A = Any[1,2,3,abs,abs,afunc,afunc,echo,echo,CycleFoo.echo,CycleFoo.echo,4,5]
    A[3] = A
    serialize(s, A)
    seekstart(s)
    b = deserialize(s)
    @test b[3] === b
    @test b[1] == 1
    @test b[4] === abs
    @test b[5] === b[4]
    @test b[5](-1) == 1

    @test b[6] === b[7]
    @test b[6]("Hello") == "Hello"

    @test b[8] === b[9]
    @test b[8]("World") == "World"

    @test b[10] === b[11]
    @test b[10]("foobar") == "foobar"

    @test b[end] == 5
    @test length(b) == length(A)
    @test isa(b,Vector{Any})
end

# shared references
create_serialization_stream() do s
    A = [1,2]
    B = [A,A]
    serialize(s, B)
    seekstart(s)
    C = deserialize(s)
    @test C == B
    @test C[1] === C[2]
end

# Regex
create_serialization_stream() do s
    r1 = r"a?b.*"
    serialize(s, r1)
    seekstart(s)
    r2 = deserialize(s)
    @test r1 == r2
end

# issue #13452
module Test13452
using Test, Serialization

module Shell
export foo
foo(x) = error("Instances must implement foo")
end

module Instance1
using ..Shell
Shell.foo(x::Int) = "This is an Int"
end

using .Shell, .Instance1
io = IOBuffer()
serialize(io, foo)
str = String(take!(io))
@test !contains(str, "Instance1")
@test contains(str, "Shell")

end  # module Test13452

# issue #15163
mutable struct B15163{T}
    x::Array{T}
end
let b = IOBuffer()
    serialize(b,B15163([1]))
    seekstart(b)
    c = deserialize(b)
    @test isa(c,B15163) && c.x == [1]
end
# related issue #20066
let b = IOBuffer()
    serialize(b, Dict{Vector, Vector}())
    seekstart(b)
    c = deserialize(b)
    @test isa(c, Dict{Vector, Vector})
end

# issue #15849
let b = IOBuffer()
    vt = Tuple[(1,)]
    serialize(b, vt)
    seekstart(b)
    @test deserialize(b) == vt
end

# issue #1770
let a = ['T', 'e', 's', 't'],
    f = IOBuffer()

    # issue #1770
    serialize(f, a)
    seek(f, 0)
    @test deserialize(f) == a
    f = IOBuffer()
    serialize(f, a)
    seek(f, 0)
    @test deserialize(f) == a

    # issue #4414
    seek(f,0)
    serialize(f, :β)
    seek(f,0)
    @test deserialize(f) === :β
end

# issue #20324
struct T20324{T}
    x::T
end
let x = T20324[T20324(1) for i = 1:2]
    b = IOBuffer()
    serialize(b, x)
    seekstart(b)
    y = deserialize(b)
    @test isa(y,Vector{T20324})
    @test y == x
end

# serializer header
let io = IOBuffer()
    serialize(io, ())
    seekstart(io)
    b = read(io)
    @test b[1] == Serialization.HEADER_TAG
    @test b[2:3] == b"JL"
    @test b[4] == Serialization.ser_version
    @test (b[5] & 0x3) == (ENDIAN_BOM == 0x01020304)
    @test ((b[5] & 0xc)>>2) == (sizeof(Int) == 8)
    @test (b[5] & 0xf0) == 0
    @test all(b[6:8] .== 0)
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: lex_runtime_service
using AWS.Compat
using AWS.UUIDs

"""
    delete_session(bot_alias, bot_name, user_id)
    delete_session(bot_alias, bot_name, user_id, params::Dict{String,<:Any})

Removes session information for a specified bot, alias, and user ID.

# Arguments
- `bot_alias`: The alias in use for the bot that contains the session data.
- `bot_name`: The name of the bot that contains the session data.
- `user_id`: The identifier of the user associated with the session data.

"""
delete_session(botAlias, botName, userId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("DELETE", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/session"; aws_config=aws_config)
delete_session(botAlias, botName, userId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("DELETE", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/session", params; aws_config=aws_config)

"""
    get_session(bot_alias, bot_name, user_id)
    get_session(bot_alias, bot_name, user_id, params::Dict{String,<:Any})

Returns session information for a specified bot, alias, and user ID.

# Arguments
- `bot_alias`: The alias in use for the bot that contains the session data.
- `bot_name`: The name of the bot that contains the session data.
- `user_id`: The ID of the client application user. Amazon Lex uses this to identify a
  user's conversation with your bot.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"checkpointLabelFilter"`: A string used to filter the intents returned in the
  recentIntentSummaryView structure.  When you specify a filter, only intents with their
  checkpointLabel field set to that string are returned.
"""
get_session(botAlias, botName, userId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("GET", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/session/"; aws_config=aws_config)
get_session(botAlias, botName, userId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("GET", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/session/", params; aws_config=aws_config)

"""
    post_content(content-_type, bot_alias, bot_name, input_stream, user_id)
    post_content(content-_type, bot_alias, bot_name, input_stream, user_id, params::Dict{String,<:Any})

 Sends user input (text or speech) to Amazon Lex. Clients use this API to send text and
audio requests to Amazon Lex at runtime. Amazon Lex interprets the user input using the
machine learning model that it built for the bot.  The PostContent operation supports audio
input at 8kHz and 16kHz. You can use 8kHz audio to achieve higher speech recognition
accuracy in telephone audio applications.   In response, Amazon Lex returns the next
message to convey to the user. Consider the following example messages:     For a user
input \"I would like a pizza,\" Amazon Lex might return a response with a message eliciting
slot data (for example, PizzaSize): \"What size pizza would you like?\".     After the user
provides all of the pizza order information, Amazon Lex might return a response with a
message to get user confirmation: \"Order the pizza?\".     After the user replies \"Yes\"
to the confirmation prompt, Amazon Lex might return a conclusion statement: \"Thank you,
your cheese pizza has been ordered.\".     Not all Amazon Lex messages require a response
from the user. For example, conclusion statements do not require a response. Some messages
require only a yes or no response. In addition to the message, Amazon Lex provides
additional context about the message in the response that you can use to enhance client
behavior, such as displaying the appropriate client user interface. Consider the following
examples:     If the message is to elicit slot data, Amazon Lex returns the following
context information:     x-amz-lex-dialog-state header set to ElicitSlot
x-amz-lex-intent-name header set to the intent name in the current context
x-amz-lex-slot-to-elicit header set to the slot name for which the message is eliciting
information     x-amz-lex-slots header set to a map of slots configured for the intent with
their current values       If the message is a confirmation prompt, the
x-amz-lex-dialog-state header is set to Confirmation and the x-amz-lex-slot-to-elicit
header is omitted.     If the message is a clarification prompt configured for the intent,
indicating that the user intent is not understood, the x-amz-dialog-state header is set to
ElicitIntent and the x-amz-slot-to-elicit header is omitted.     In addition, Amazon Lex
also returns your application-specific sessionAttributes. For more information, see
Managing Conversation Context.

# Arguments
- `content-_type`:  You pass this value as the Content-Type HTTP header.   Indicates the
  audio format or text. The header value must start with one of the following prefixes:
  PCM format, audio data must be in little-endian byte order.   audio/l16; rate=16000;
  channels=1   audio/x-l16; sample-rate=16000; channel-count=1   audio/lpcm;
  sample-rate=8000; sample-size-bits=16; channel-count=1; is-big-endian=false      Opus
  format   audio/x-cbr-opus-with-preamble; preamble-size=0; bit-rate=256000;
  frame-size-milliseconds=4     Text format   text/plain; charset=utf-8
- `bot_alias`: Alias of the Amazon Lex bot.
- `bot_name`: Name of the Amazon Lex bot.
- `input_stream`:  User input in PCM or Opus audio format or text format as described in
  the Content-Type HTTP header.  You can stream audio data to Amazon Lex or you can create a
  local buffer that captures all of the audio data before sending. In general, you get better
  performance if you stream audio data rather than buffering the data locally.
- `user_id`: The ID of the client application user. Amazon Lex uses this to identify a
  user's conversation with your bot. At runtime, each request must contain the userID field.
  To decide the user ID to use for your application, consider the following factors.   The
  userID field must not contain any personally identifiable information of the user, for
  example, name, personal identification numbers, or other end user personal information.
  If you want a user to start a conversation on one device and continue on another device,
  use a user-specific identifier.   If you want the same user to be able to have two
  independent conversations on two different devices, choose a device-specific identifier.
  A user can't have two independent conversations with two different versions of the same
  bot. For example, a user can't have a conversation with the PROD and BETA versions of the
  same bot. If you anticipate that a user will need to have conversation with two different
  versions, for example, while testing, include the bot alias in the user ID to separate the
  two conversations.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Accept"`:  You pass this value as the Accept HTTP header.   The message Amazon Lex
  returns in the response can be either text or speech based on the Accept HTTP header value
  in the request.     If the value is text/plain; charset=utf-8, Amazon Lex returns text in
  the response.     If the value begins with audio/, Amazon Lex returns speech in the
  response. Amazon Lex uses Amazon Polly to generate the speech (using the configuration you
  specified in the Accept header). For example, if you specify audio/mpeg as the value,
  Amazon Lex returns speech in the MPEG format.   If the value is audio/pcm, the speech
  returned is audio/pcm in 16-bit, little endian format.    The following are the accepted
  values:   audio/mpeg   audio/ogg   audio/pcm   text/plain; charset=utf-8   audio/*
  (defaults to mpeg)
- `"x-amz-lex-active-contexts"`: A list of contexts active for the request. A context can
  be activated when a previous intent is fulfilled, or by including the context in the
  request, If you don't specify a list of contexts, Amazon Lex will use the current list of
  contexts for the session. If you specify an empty list, all contexts for the session are
  cleared.
- `"x-amz-lex-request-attributes"`: You pass this value as the x-amz-lex-request-attributes
  HTTP header. Request-specific information passed between Amazon Lex and a client
  application. The value must be a JSON serialized and base64 encoded map with string keys
  and values. The total size of the requestAttributes and sessionAttributes headers is
  limited to 12 KB. The namespace x-amz-lex: is reserved for special attributes. Don't create
  any request attributes with the prefix x-amz-lex:. For more information, see Setting
  Request Attributes.
- `"x-amz-lex-session-attributes"`: You pass this value as the x-amz-lex-session-attributes
  HTTP header. Application-specific information passed between Amazon Lex and a client
  application. The value must be a JSON serialized and base64 encoded map with string keys
  and values. The total size of the sessionAttributes and requestAttributes headers is
  limited to 12 KB. For more information, see Setting Session Attributes.
"""
post_content(Content_Type, botAlias, botName, inputStream, userId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("POST", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/content", Dict{String, Any}("inputStream"=>inputStream, "headers"=>Dict{String, Any}("Content-Type"=>Content_Type)); aws_config=aws_config)
post_content(Content_Type, botAlias, botName, inputStream, userId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("POST", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/content", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("inputStream"=>inputStream, "headers"=>Dict{String, Any}("Content-Type"=>Content_Type)), params)); aws_config=aws_config)

"""
    post_text(bot_alias, bot_name, input_text, user_id)
    post_text(bot_alias, bot_name, input_text, user_id, params::Dict{String,<:Any})

Sends user input to Amazon Lex. Client applications can use this API to send requests to
Amazon Lex at runtime. Amazon Lex then interprets the user input using the machine learning
model it built for the bot.   In response, Amazon Lex returns the next message to convey to
the user an optional responseCard to display. Consider the following example messages:
For a user input \"I would like a pizza\", Amazon Lex might return a response with a
message eliciting slot data (for example, PizzaSize): \"What size pizza would you like?\"
  After the user provides all of the pizza order information, Amazon Lex might return a
response with a message to obtain user confirmation \"Proceed with the pizza order?\".
After the user replies to a confirmation prompt with a \"yes\", Amazon Lex might return a
conclusion statement: \"Thank you, your cheese pizza has been ordered.\".     Not all
Amazon Lex messages require a user response. For example, a conclusion statement does not
require a response. Some messages require only a \"yes\" or \"no\" user response. In
addition to the message, Amazon Lex provides additional context about the message in the
response that you might use to enhance client behavior, for example, to display the
appropriate client user interface. These are the slotToElicit, dialogState, intentName, and
slots fields in the response. Consider the following examples:    If the message is to
elicit slot data, Amazon Lex returns the following context information:    dialogState set
to ElicitSlot     intentName set to the intent name in the current context     slotToElicit
set to the slot name for which the message is eliciting information     slots set to a map
of slots, configured for the intent, with currently known values       If the message is a
confirmation prompt, the dialogState is set to ConfirmIntent and SlotToElicit is set to
null.    If the message is a clarification prompt (configured for the intent) that
indicates that user intent is not understood, the dialogState is set to ElicitIntent and
slotToElicit is set to null.     In addition, Amazon Lex also returns your
application-specific sessionAttributes. For more information, see Managing Conversation
Context.

# Arguments
- `bot_alias`: The alias of the Amazon Lex bot.
- `bot_name`: The name of the Amazon Lex bot.
- `input_text`: The text that the user entered (Amazon Lex interprets this text).
- `user_id`: The ID of the client application user. Amazon Lex uses this to identify a
  user's conversation with your bot. At runtime, each request must contain the userID field.
  To decide the user ID to use for your application, consider the following factors.   The
  userID field must not contain any personally identifiable information of the user, for
  example, name, personal identification numbers, or other end user personal information.
  If you want a user to start a conversation on one device and continue on another device,
  use a user-specific identifier.   If you want the same user to be able to have two
  independent conversations on two different devices, choose a device-specific identifier.
  A user can't have two independent conversations with two different versions of the same
  bot. For example, a user can't have a conversation with the PROD and BETA versions of the
  same bot. If you anticipate that a user will need to have conversation with two different
  versions, for example, while testing, include the bot alias in the user ID to separate the
  two conversations.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"activeContexts"`: A list of contexts active for the request. A context can be activated
  when a previous intent is fulfilled, or by including the context in the request, If you
  don't specify a list of contexts, Amazon Lex will use the current list of contexts for the
  session. If you specify an empty list, all contexts for the session are cleared.
- `"requestAttributes"`: Request-specific information passed between Amazon Lex and a
  client application. The namespace x-amz-lex: is reserved for special attributes. Don't
  create any request attributes with the prefix x-amz-lex:. For more information, see Setting
  Request Attributes.
- `"sessionAttributes"`: Application-specific information passed between Amazon Lex and a
  client application. For more information, see Setting Session Attributes.
"""
post_text(botAlias, botName, inputText, userId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("POST", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/text", Dict{String, Any}("inputText"=>inputText); aws_config=aws_config)
post_text(botAlias, botName, inputText, userId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("POST", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/text", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("inputText"=>inputText), params)); aws_config=aws_config)

"""
    put_session(bot_alias, bot_name, user_id)
    put_session(bot_alias, bot_name, user_id, params::Dict{String,<:Any})

Creates a new session or modifies an existing session with an Amazon Lex bot. Use this
operation to enable your application to set the state of the bot. For more information, see
Managing Sessions.

# Arguments
- `bot_alias`: The alias in use for the bot that contains the session data.
- `bot_name`: The name of the bot that contains the session data.
- `user_id`: The ID of the client application user. Amazon Lex uses this to identify a
  user's conversation with your bot.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Accept"`: The message that Amazon Lex returns in the response can be either text or
  speech based depending on the value of this field.   If the value is text/plain;
  charset=utf-8, Amazon Lex returns text in the response.   If the value begins with audio/,
  Amazon Lex returns speech in the response. Amazon Lex uses Amazon Polly to generate the
  speech in the configuration that you specify. For example, if you specify audio/mpeg as the
  value, Amazon Lex returns speech in the MPEG format.   If the value is audio/pcm, the
  speech is returned as audio/pcm in 16-bit, little endian format.   The following are the
  accepted values:    audio/mpeg     audio/ogg     audio/pcm     audio/* (defaults to mpeg)
   text/plain; charset=utf-8
- `"activeContexts"`: A list of contexts active for the request. A context can be activated
  when a previous intent is fulfilled, or by including the context in the request, If you
  don't specify a list of contexts, Amazon Lex will use the current list of contexts for the
  session. If you specify an empty list, all contexts for the session are cleared.
- `"dialogAction"`: Sets the next action that the bot should take to fulfill the
  conversation.
- `"recentIntentSummaryView"`: A summary of the recent intents for the bot. You can use the
  intent summary view to set a checkpoint label on an intent and modify attributes of
  intents. You can also use it to remove or add intent summary objects to the list. An intent
  that you modify or add to the list must make sense for the bot. For example, the intent
  name must be valid for the bot. You must provide valid values for:    intentName    slot
  names    slotToElict    If you send the recentIntentSummaryView parameter in a PutSession
  request, the contents of the new summary view replaces the old summary view. For example,
  if a GetSession request returns three intents in the summary view and you call PutSession
  with one intent in the summary view, the next call to GetSession will only return one
  intent.
- `"sessionAttributes"`: Map of key/value pairs representing the session-specific context
  information. It contains application information passed between Amazon Lex and a client
  application.
"""
put_session(botAlias, botName, userId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("POST", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/session"; aws_config=aws_config)
put_session(botAlias, botName, userId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_service("POST", "/bot/$(botName)/alias/$(botAlias)/user/$(userId)/session", params; aws_config=aws_config)


text:
"""
We map a charge distribution (scalar field) to electric potential (scalar field) and electric field (vector field).
In fluid mechanics the equivalent mapping is from source/sink distribution to flow potential and velocity field.
"""
dir = ".."
include("$dir/src/neural_networks.jl")
include("$dir/src/plotutils.jl")
Random.seed!(1)

# grid params
dx = 0.1
sz = fill(11,3) #
grid = Grid(dx, sz)
@show grid.dV, grid.origin

rank = 0
charges = zeros(sz...,1)
# @show charges.components

# place dipole point charges 2r apart in grid as input
pos = zeros(3)
vals = [1.0]
put_point_source!(charges,grid, pos, vals)
# field_plot(charges)

# tensor convolution params

# generate data with Green's functions for Poisson's Equation, Gauss's law
rmax = 1.0
name = :inverse_squared_field
field_op = LinearOperator(name; dx, rmax)
name = :potential
potential_op = LinearOperator(name; dx, rmax)

# output data
potential = potential_op(charges,grid)
efield = field_op(charges,grid)
# field_plot(potential)
# field_plot(field)

# check
r=.2
rvec = [r, 0, 0]
# automatic interpolation at rvec
# @assert get_tensor(rvec,potential,grid) ≈ [1 / (4π * r) - 1 / (4π * 3r)]
# @assert get_tensor(rvec,efield,grid) ≈ [1 / (4π * r^2) - 1 / (4π * (3r)^2), 0, 0]
@assert get(potential,grid,rvec) ≈ [1 / (4π * r)]
@assert get(efield,grid,rvec) ≈ [1 / (4π * r^2), 0, 0]

rank_max = 1
in_ =Props((1,),grid) # input scalar field
out =Props((1,1),grid) # input scalar field
X = charges
Y = cat(potential, efield,dims=4)
# define tensor field convolution layer
L = EquivConv(in_, out, rmax)

function loss()
    Yhat= L(X)
    l = nae(Yhat,Y)
    println(l)
    l
end
loss()
##
ps = Flux.params(L)
data = [()]
opt = ADAM(0.1)

for i = 1:30
    Flux.train!(loss, ps, data, opt)
end

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
using Spirit
using Test
using Sockets

httptest1 = include("./data/http/test1.jl")
httptest2 = include("./data/http/test2.jl")
httptest3 = include("./data/http/test3.jl")

# write your own tests here
@testset "Headers" begin
    include("test_headers.jl")
end

@testset "Continuation" begin
    include("test_continuation.jl")
end

@testset "Pipeline" begin
    include("test_pipeline.jl")
end

@testset "Util" begin
    include("test_util.jl")
end

@testset "Connection" begin
    include("test_connection.jl")
end

@testset "ChunkedTransferDecodeStream" begin
    include("test_chunkedtransferdecodestream.jl")
end

@testset "SizedStream" begin
    include("test_sizedstream.jl")
end

@testset "URI" begin
    include("test_uri.jl")
end

@testset "Request" begin
    include("test_request.jl")
end