text:
# # rewrites the expressions generated by clang

# using DataStructures

# # node: the input obuff has the type annotations in the function signature and the types in the ccall as references to the same Expr object.  A deepcopy and overwrite is needed to fix hthis

# # dictionary to map pointer types to desired type
# # isdefined(modulename, :name) -> Bool
# # use wrap_contex.commonbuf for dictionary of names that have been defined

# ############################################################################### Only these parameters need to be modified for the different
# # version of PETSc

# # used to modify function signatures
# type_dict = Dict{Any, Any} (
#   :PetscScalar => :Float32,
#   :PetscReal => :Float32,
#   :PetscInt => :Int64,
# )

# const petsc_libname = :petscRealSingle

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

# val_tmp = type_dict[:PetscScalar]
# type_dict_single = Dict{Any, Any} (
# :(Ptr{UInt8}) => Union{Cstring, String, Symbol, Array{UInt8}, Ptr{UInt8}}
# )

# # used to convert typealiases to immutable type definionts
# # currently, if the key exists, it is converted
# # if value == 1, create new immutable type
# # otherwise replace key with value
# # also used for function signatures
# new_type_dict = Dict{Any, Any} (
#   :Vec => :(Vec{$val_tmp}),
#   :Mat => :(Mat{$val_tmp}),
#   :KSP => :(KSP{$val_tmp}),
#   :PC => :(PC{$val_tmp}),
#   :PetscViewer => :(PetscViewer{$val_tmp}),
#   :PetscOption => :(PetscOption{$val_tmp}),
#   :IS => :(IS{$val_tmp}),
#   :ISLocalToGlobalMapping => :(ISLocalToGlobalMapping{$val_tmp}),
#   :ISColoring => :(ISColoring{$val_tmp}),
#   :PetscLayout => :(PetscLayout{$val_tmp}),
#   :VecScatter => :(VecScatter{$val_tmp}),
#   :AO => :(AO{$val_tmp}),
#   :TS => :(TS{$val_tmp})
# )

# # definitions that will be provided, but don't come from Petsc
# # mostly MPI stuff
# const_defs = Dict{Any, Any} (
#   :MPI_COMM_SELF => 1,
#   :MPI_Comm => 1,
#   :comm_type => 1,
# )

# # dictionary to hold names of types
# # that are aliased to Symbol
# symbol_type_dict = Dict{Any, Any} ()

# # create a string array mirroring a Symbol array
# function send_symbol(x::AbstractArray)
#   string_arr = similar(x, ASCIIString)
#   for i=1:length(x)
#     if isdefined(x[i])
#       string_arr[i] = copy(x[i])
#     end
#   end
#   return string_arr
# end

# # get a string array and turn it into a Symbol array
# function return_symbol(string_array::AbstractArray, x::AbstractArray)
#   for i=1:length(x)
#     x[i] = string(string_array[i])
#   end
# end

# # things to be recursively replaced in function signatures
# sig_rec_dict = Dict{Any, Any} ()

# for i in keys(type_dict)
#   get!(sig_rec_dict, i, type_dict[i])
# end

# for i in keys(new_type_dict)
#   get!(sig_rec_dict, i, new_type_dict[i])
# end

# # things to be replaced in function signatures only if they
# # are top level (ie. this does a non recursive replace)
# sig_single_dict = Dict{Any, Any} (
#   :(Ptr{UInt8}) => :(Union{String, Cstring, Symbol, Array{UInt8}, Ptr{UInt8}}),
#   :Int32 => :Integer,
#   :Int64 => :Integer,
#   :Cint => :Integer,
#   :PetscInt => :Integer,
# )

# # list of symbols to search for
# # if none are found, then a dummy argument is added
# sig_dummyarg_dict = Dict{Any, Any}(
#   :PetscScalar => 1,
# )

# # add keys from new_type_dict because dummy arg check
# # is done after type replacement
# for i in values(new_type_dict)
#   get!(sig_dummyarg_dict, i, 1)
# end

# # things to be recursively replaced in the ccall argument list
# ccall_rec_dict = Dict{Any, Any} (
# :MPI_Comm => :comm_type
# )

# for i in keys(type_dict)
#   get!(ccall_rec_dict, i, type_dict[i])
# end

# for i in keys(new_type_dict)
#   get!(ccall_rec_dict, i, new_type_dict[i])
# end

# # things to be replace in the ccall argument list only if
# # they are top level (ie. this does a non recursive replace)
# ccall_single_dict = Dict{Any, Any} (
#   :(Ptr{UInt8}) => :Cstring,
# )

# # things to be replaced recurisvely in typealias rhs
# typealias_rec_dict = Dict{Any, Any} ()

# for i in keys(type_dict)
#   get!(typealias_rec_dict, i, type_dict[i])
# end

# for i in keys(new_type_dict)
#   get!(typealias_rec_dict, i, new_type_dict[i])
# end

# # dictionary of typealiases to exclude based on the lhs argument
# typealias_lhs_dict = Dict{Any, Any} ()

# for i in keys(type_dict)
#   get!(typealias_lhs_dict, i, i)
# end

# # things to be replaced in typealias rhs, non recursive
# # this creates a potential loophole for string handling
# # because a Ptr{typealias} == Ptr{Ptr{UInt8}} will be
# # handled incorrectly
# typealias_single_dict = Dict{Any, Any} (
# :(Ptr{UInt8}) => :Symbol
# )

# # key = typealias rhs values to exclude
# # if values == 1, then create an immutable type
# typealias_exclude_dict = Dict{Any, Any} (
# )

# for i in keys(new_type_dict)
#   get!(typealias_exclude_dict, i, 1)
# end

# # things to be recurisvely replace in constant declaration rhs
# const_rec_dict = Dict{Any, Any} (
# :NULL => :C_NULL
# )

# for i in keys(type_dict)
#   get!(const_rec_dict, i, type_dict[i])
# end

# for i in keys(new_type_dict)
#   get!(const_rec_dict, i, new_type_dict[i])
# end

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

# function petsc_rewriter(obuf)
#   for i=1:length(obuf)
#     ex_i = obuf[i]  # obuf is an array of expressions
#     println("rewriting expression ", ex_i)

#     if typeof(ex_i) == Expr
#       head_i = ex_i.head  # function
#       # each ex_i has 2 args, a function signature and a function body

#       # figure out what kind of expression it is, do the right modification
#       if head_i == :function  # function declaration
#         obuf[i] = process_func(ex_i)
#       elseif ex_i.head == :const  # constant declaration
#         obuf[i] = process_const(ex_i)
#       elseif ex_i.head == :typealias  # typealias definition
#         obuf[i] = fix_typealias(ex_i)
#       elseif ex_i.head == :type
#         obuf[i] = process_type(ex_i)

#       else  # some other kind of expression
#         println("not processing expression", ex_i)

#         # convert to concrete types
#         for j in keys(type_dict)
#           obuf[i] = replace_symbol(ex_i, j, type_dict[j])
#         end

#         # purge anything unknown
#         # this will always omit the expression because a newly declared name
#         # will be unknown
#         tmp = are_syms_defined(ex_i)
#         if tmp != 0
#           obuf[i] = "#= skipping undefined expression $ex_i =#"
#         end

#         # for types at least, arg[2] is the type name
#         obj_name = ex_i.args[2]
#         if haskey(wc.common_buf, obj_name)
#           delete!(wc.common_buf, obj_name)
#         end

#       end

#     elseif typeof(ex_i) == ASCIIString
#       obuf[i] = process_string(ex_i)
#     else
#       println("not processing ", typeof(ex_i))
#     end  # end if Expr

#     println("final expression = ", obuf[i])
#   end

#   return obuf  # return modified obuf
# end

# ##### functions to rewrite function signature #####
# function process_func(ex)
#   @assert ex.head == :function  # this is a function declaration

#   println("at beginning of process_func, processing function ", ex.args[1])
#   println("with body ", ex.args[2])

#   ex.args[1] = rewrite_sig(ex.args[1])  # function signature

#   println("after rewrite_sig, processing function ", ex.args[1])
#   println("with body ", ex.args[2])

#   ex.args[2] = rewrite_body(ex.args[2])  # function body

#   println("processing function ", ex.args[1])
#   println("with body ", ex.args[2])
#   # now check if any undefined type annotations remain
#   sum = 0
#   ex_sig = ex.args[1]
#   ex_body = ex.args[2]

#   for i=2:length(ex_sig.args)  # loop over all arguments to the function
#     # get the type annotation Expr or Symbol of argument i-1
#     ex_sig_i = ex_sig.args[i].args[2]
#     sum += are_syms_defined(ex_sig_i)
#   end

#   # now check ccall argument types
#   ex_ccall = ex_body.args[1]
#   ex_ccall_types = ex_ccall.args[3]  # the tuple of argument types
#   sum += are_syms_defined(ex_ccall_types)

#   if sum != 0
#     return "#= skipping function with undefined symbols: \n $ex \n=#"
#   end

#   # now add any extra expression needed
#   ex = add_body(ex)
#   return ex
# end

# # add any needed expressions to the body of the function
# # this is post rewrite_sig, rewrite_body
# function add_body(ex)

#   ex_sig = ex.args[1]
#   ex_body = ex.args[2]

#   fname = ex_sig.args[1]
#   @assert typeof(fname) == Symbol
#   fname_str = string(fname)
#   # check if a type annotation == Union(Ptr{symbol_type}...)
#   println("checking function ", fname_str, " for symbols-string conversions")

#   # check number of ccall args vs. function signature args
#   numargs_ccall = length(ex_body.args[1].args) - 3
#   numargs_func = length(ex_sig.args) - 1
#   offset = numargs_func - numargs_ccall  # skip any function args in excess of ccall args

#   for i in keys(symbol_type_dict)
#     for j=(2 + offset):length(ex_sig.args)  # loop over arguments to function
# #      println("j = ", j)
#       type_annot_j = ex_sig.args[j].args[2]  # get the teyp annotation
#       argname_j = ex_sig.args[j].args[1]
#       # check for arrays of symbols that need to be copied into a string array
#       if type_annot_j == :(Union{Ptr{$i}, StridedArray{$i}, Ptr{$i}, Ref{$i}})
#         if contains(fname_str, "Get")  # array is to be populated
#           # add calls to function body
#           resize!(ex_body.args, length(ex_body.args) + 2)
#           for i=2:(length(ex_body.args) - 1)
#             ex_body.args[i] = deepcopy(ex_body.args[i-1])  # shift ccall to 2nd arg
#           end
#           ex_body.args[1] = nothing
#           ex_body.args[end] = nothing

#           # figure out where the ccall is
#           ccall_index = get_ccall_index(ex_body)

#           # construct the before function call
#           call_ex =  Expr(:call, :symbol_get_before, argname_j)
#           new_argname = Symbol(string( argname_j, "_"))  # add underscore
#           ex_body.args[1] = Expr(:(=), :($new_argname, tmp), deepcopy(call_ex))

#           # construct the after call
#           ex_body.args[length(ex_body.args)] = Expr(:call, :symbol_get_after, new_argname, argname_j)

#           # modify the ccall argument name
#           println("ccall_index = ", ccall_index)
#           ccall_ex = ex_body.args[ccall_index]

#           @assert ccall_ex.head == :ccall
#           ccall_ex.args[j + 3 - 1 - offset] = new_argname

#         elseif contains(fname_str, "Set")  # array is already populated
#           # add calls to function body
#           resize!(ex_body.args, length(ex_body.args) + 1)
#           for i=2:(length(ex_body.args))
#             ex_body.args[i] = deepcopy(ex_body.args[i-1])  # shift ccall to 2nd arg
#           end
#           ex_body.args[1] = nothing

#           # figure out where the ccall is
#           ccall_index = get_ccall_index(ex_body)

#           # construct the before function call
#           call_ex =  Expr(:call, :symbol_set_before, argname_j)
#           new_argname = Symbol(string( argname_j, "_"))  # add underscore
#           ex_body.args[1] = Expr(:(=), new_argname, deepcopy(call_ex))

#           # modify the ccall argument name
#           ccall_ex = ex_body.args[ccall_index]
#           @assert ccall_ex.head == :ccall
#           ccall_ex.args[j + 3 - 1 - offset] = new_argname

#         else
#           println(stderr, "Warning, Symbol type conversion not handled in function ", fname_str)

#         end  # end if contains(get)
#       end  # end if type_annot_j == ...
#     end  # end for j
#   end # end loop over symbol_type_dict

#   # assign ccall error code to a variable, then return the variable
#   # this should be the last modification made, because get_ccall_index
#   # won't work after this
#   ccall_index = get_ccall_index(ex_body)
#   ccall_ex = deepcopy(ex_body.args[ccall_index])
#   ex_body.args[ccall_index] = :(err = $ccall_ex)
#   resize!(ex_body.args, length(ex_body.args) + 1)
#   ex_body.args[end] = :(return err)

#   # convert array of symbols to array of strings
#   # check signature for
#   return deepcopy(ex)
# end

# function get_ccall_index(ex)
#   # returns the index of the first ccall in the arguments of ex

#   index = 0

#   for i=1:length(ex.args)
#     if typeof(ex.args[i]) == Expr
#       if ex.args[i].head  == :ccall
#         index = i
#         break
#       end
#     end
#   end

#   if index == 0
#     println("Warning, ccall not found")
#   end

#   return index
# end

# function rewrite_sig(ex)  # rewrite the function signature

#   @assert ex.head == :call  # verify this is a function signature expression

#   # ex.args[1] = function name, as a Symbol

#   # check if function contains a PetscScalar

#   # if yes, add paramterization, do search and replace PetscScalar -> S
#   # also add macro

#   println("rewrite_sig ex = ", ex)

#   # Process all arguments of the function
#   # Replace all pointers with Unions(Ptr, AbstractArray...)
#   # do any other transformations in ptr_dict
#   for i=2:length(ex.args)
#     # each of ex.args is an expression containing arg name, argtype
#     println("typeof(ex.args[$i]) = ", typeof(ex.args[i]))
#     @assert typeof(ex.args[i]) == Expr  || typeof(ex.args[i]) == Symbol  # verify these are all expressions
#     ex.args[i] = process_sig_arg(ex.args[i])  # process each expression
#   end

#   # check for any Symbol that will uniquely identify which
#   # version of petsc to call
#   println("checking for uniqueness of signature")
#   println("ex = ", ex)

#   val = contains_symbol(ex, :PetscScalar)
#   for i in keys(sig_dummyarg_dict)
#     val += contains_symbol(ex, i)
#   end

#   println("contains_symbol = ", val)

#   if val == 0  # if no arguments will make the function signature unique

#     println("adding dummy arg")
#     ex = add_dummy_arg(ex)
#   end

#   #=
#   println("replacing typealiases")
#   println("ex = ", ex)
#   # do second pass to replace Petsc typealiases with a specific type
#   for i=2:length(ex.args)
#   for j in keys(type_dict)  # check for all types

#   println("replacing ", j, ", with ", type_dict[j])
#   ex.args[i] = replace_symbol(ex.args[i], j, type_dict[j])
# end
# end

# println("after modification rewrite_sig ex = ", ex)
# =#
# return ex
# end

# function process_sig_arg(ex)
#   # take the expression that is an argument to the function and rewrite it

#   @assert ex.head == :(::)

#   # modify args here
#   #   arg_name = ex.args[1]  # Symbol
#   #   arg_type = ex.args[2]  # Expr contianing type tag
#   println("ex.args[2] = ", ex.args[2])
#   println("type = ", typeof(ex.args[2]))

#   #   if typeof(ex.args[1]) == Symbol

#   # get only typetag expression, modify it
#   ex.args[2] = modify_typetag(ex.args[2])

#   return ex
# end

# function modify_typetag(ex)
#   # do non recursive replace first
#   ex = deepcopy(get(sig_single_dict, ex, ex))

#   # now do recursive search and replace
#   for i in keys(sig_rec_dict)
#     ex = deepcopy(replace_symbol(ex, i, sig_rec_dict[i]))
#   end

#   # should make sure there are no nested pointers?
#   #  @assert ex.head == :curly || ex.head == :Symbol # verify this is a typetag
#   if typeof(ex) == Expr
#     # replace pointer with Union of ptr, array, c_null
#     if ex.head == :curly && ex.args[1] == :Ptr
#       ptr_type = deepcopy(ex.args[2])
#       ex =  :(Union{Ptr{$ptr_type}, StridedArray{$ptr_type}, Ptr{$ptr_type}, Ref{$ptr_type}})
#     end
#   end

#   return ex
# end

# function add_param(ex)
#   # add the {S <: PetscScalar} to a function declaration
#   @assert typeof(ex) == Symbol  # make sure we don't already have a paramaterization

#   # could do more extensive operations here
#   return :($ex{S <: PetscScalars})
# end

# function add_dummy_arg(ex)
#   # add a dummy argument to the function argument list
#   # ex is the function name + args (ie. the entire function signature)

#   println("adding dummy argument to ", ex)
#   @assert ex.head == :call
#   len = length(ex.args)
#   resize!(ex.args, len + 1)

#   println("length(ex.args) = ", length(ex.args))
#   println("ex.args = ", ex.args)
#   # shift  args up by one
#   # first arg is function name, so don't modify it

#   if length(ex.args) >= 3 # if there are any args to shift
#     println("shifting arguments, length(ex.args) = ", length(ex.args))
#     for i=length(ex.args):(-1):3
#       ex.args[i] = ex.args[i-1]
#       println("ex.args[$i] = ", ex.args[i])
#     end
#   end

#   # insert dummy argument into first position
#   # the PetscScalar will get rewritten to
#   # appropriate type later
#   val = type_dict[:PetscScalar]
#   ex.args[2] = :(arg0::Type{$val})

#   println("finished adding dummy arg, ex = ", ex)

#   return ex
# end

# #####  function to  rewrite the body of the function #####
# function rewrite_body(ex)  # rewrite body of a function
#   @assert ex.head == :block
#   # ex has only one argument, the ccall
#   # could insert other statements arond the ccall here?
#   ex.args[1] = process_ccall(ex.args[1])
#   return ex
# end

# function process_ccall(ex)
#   @assert ex.head == :ccall  # verify this is the ccall statement
#   # args[1] = Expr, tuple of fname, libname
#   # args[2] = return type, Symbol
#   # args[3] = Expr, tuple of types to ccall
#   # args[4] = Symbol, first argument name,
#   # ...
#   println("processing ccall ", ex)
#   ex.args[1] = modify_libname(ex.args[1])  # change the library name
#   ex.args[2] = modify_rettype(ex.args[2]) # change return type
#   ex.args[3] = modify_types(ex.args[3])
#   println("changing argument names")
#   ex3 = ex.args[3]  # get expression of argument types
#   #  ex4 = ex.args[4]  # get expression argument names

#   @assert ex3.head == :tuple

#   return ex
# end

# function modify_rettype(ex)

#   println("modifying return type")
#   println("ex = ", ex)
#   return deepcopy(ex)
# end

# function modify_types(ex)
#   @assert ex.head == :tuple

#   # do non recursive replace first
#   for i=1:length(ex.args)
#     println("checking if ", ex.args[i], " has entry in ccall_single_dict")
#     println("ex.args[$i] before = ", ex.args[i])
#     ex.args[i] = get(ccall_single_dict, ex.args[i], ex.args[i])
#     println("  after = ", ex.args[i])
#   end

#   println("  after all non-recursive replacements = ", ex.args)

#   println("  ccall_rec_dict = ", ccall_rec_dict)
#   # do recursive replacement
#   for i in keys(ccall_rec_dict)
#     for j=1:length(ex.args)
#       ex.args[j] = replace_symbol(ex.args[j], i, ccall_rec_dict[i])
#     end
#   end

#   println("  after recursive replacement = ", ex.args)

#   return deepcopy(ex)
# end

# function modify_libname(ex)
#   @assert ex.head == :tuple
#   ex.args[2] = petsc_libname  # replace libname with a new one
#   return ex
# end

# ##### functions to make consts into global consts #####

# function process_const(ex)
#   @assert ex.head == :const

#   # do any replacements
#   for j in keys(const_rec_dict)
#     replace_symbol(ex, j, const_rec_dict[j])
#   end

#   ex2 = ex.args[1]  # get the assignment

#   rhs = ex2.args[2]
#   tmp = are_syms_defined(rhs)
#   if tmp != 0  # something is undefined
#     lhs = ex2.args[1]
#     delete!(wc.common_buf, lhs)  # remove lh from list of defined symbols
#     return "#skipping undefined $ex"  # replace expression with constant
#   end

#   # turn strings into symbols
#   if typeof(rhs) <: AbstractString
#     ex2.args[2] = :(Symbol($rhs))
#   end

#   return deepcopy(ex)
# end

# function fix_typealias(ex)
#   @assert ex.head == :typealias
#   new_type = ex.args[1]

#   println("fixing typealias ", ex)

#   if haskey(typealias_exclude_dict, new_type)
#     if typealias_exclude_dict[new_type] == 1
#       # construct immutable type definition
#       fields = Expr(:(::), :pobj, :(Ptr{Void}))
#       body = Expr(:block, fields)
#       typename = Expr(:curly, new_type, :T)  # add static parameter T
#       ex_new = Expr(:type, false, typename, body)
#       return ex_new
#     else
#       if haskey(wc.common_buf, new_type)
#         delete!(wc.common_buf, new_type)  # record that Symbol is now undefined
#       end
#       return "# excluding $ex"
#     end
#   end

#   # if we didn't create immutable type, transform the rhs according to the
#   # dictionaries
#   lhs = deepcopy(ex.args[1])
#   rhs = deepcopy(ex.args[2])

#   # check lhs for exclusion criteria
#   if haskey(typealias_lhs_dict, lhs)
#     # record that Symbol is now undefined
#     delete!(wc.common_buf, lhs)
#     return "# excluding lhs of  $ex"
#   end

#   # do non recursive replacement
#   rhs = get(typealias_single_dict, rhs, rhs)

#   if rhs == :Symbol
#     get!(ccall_rec_dict, lhs, :(Ptr{UInt8}))  # make the ccall argument a UInt8 recursive
#     get!(ccall_single_dict, lhs, :Cstring)  # make it a cstring for single level replacement
#     get!(symbol_type_dict, lhs, 1)  # record that this type is a Symbol
#   end

#   # do recursive replacement
#   for i in keys(typealias_rec_dict)
#     rhs = replace_symbol(rhs, i, typealias_rec_dict[i])
#   end

#   # check for undefined symbols

#   tmp = are_syms_defined(rhs)
#   if tmp != 0
#     delete!(wc.common_buf, lhs)  # record that the lhs Symbol is now undefined
#     return "# skipping undefined typealias $ex"
#   end

#   # if no conditions met
#   # construct a new typealias

#   return :(typealias $lhs $rhs)
# end

# # string processing
# function process_string(ex)
#   return "#= $ex =#"
# end

# ##### Type declaration functions #####
# function process_type(ex)

#   println("processing type declaration ", ex)

#   # make all types (structs) immutable
#   ex.args[1] = false

#   ex_body = ex.args[3]  # body of type declaration

#   # check that every type annotation is fully defined
#   for i=1:length(ex_body.args)
#     # if this is a type annotation expression
#     if typeof(ex_body.args[i]) == Expr && ex_body.args[i].head == :(::)
#       type_sym = ex_body.args[i].args[2]
#       tmp = are_syms_defined(type_sym)
#       if tmp != 0
#         new_name = ex.args[2]  # get the name of the type being declared
#         delete!(wc.common_buf, new_name)

#         return "#= skipping type declaration with undefined symbols:\n$ex \n=#"
#       end
#     end
#   end

#   return ex
# end

# ##### Misc. functions ####
# function contains_symbol(ex, sym)
#   # recursively check Symbol in expression ex to see if sym is present
#   sum = 0
#   if ex == sym
#     return 1
#   elseif typeof(ex) == Expr
#     #check the arguments of the exprssion
#     for i=1:length(ex.args)
#       sum += contains_symbol(ex.args[i], sym)
#     end
#   else  # something unknown
#     return 0  # assuming we did not find the Symbol
#   end

#   return sum
# end

# function check_annot(ex, sym::Symbol)
#   # check a type annotation
#   for i=1:length(ex.args)
#     if ex.args[i] == sym
#       return true
#     end
#   end
#   return false
# end

# function replace_symbol(ex, sym_old, sym_new)
#   if ex == sym_old
#     return deepcopy(sym_new)
#   elseif typeof(ex) == Expr  # recurse the arguments
#     for i=1:length(ex.args)
#       ex.args[i] = replace_symbol(ex.args[i], sym_old,sym_new)
#     end
#   else  # don't know/care what this is
#     return deepcopy(ex)
#   end

#   return deepcopy(ex)
# end

# function count_depth(seed::Integer, ex)
#   # primitive attempt to count number of nodes on tree
#   if typeof(ex) == Expr
#     for i=1:length(ex.args)

#       seed += count_depth(seed, ex.args[i])
#     end
#   else
#     println("ex = ", ex, " returning 1")
#     return 1
#   end
#   return seed
# end

# function are_syms_defined(ex)
#   # check all symbols, see if they are defined by julia or by
#   # the wrap contex
#   undefined_syms = 0  # approximate counter of undefined symbols found

#   if typeof(ex)  == Expr  # keep recursing
#     for i=1:length(ex.args)
#       undefined_syms +=  are_syms_defined(ex.args[i])
#     end  # end loop over args
#   elseif typeof(ex) == Symbol  # if this is a Symbol
#     if isdefined(ex) || haskey(wc.common_buf, ex) || haskey(const_defs, ex)
#       return 0  # Symbol is defined
#     else
#       return 1
#     end
#   else # we don't know/care what this expression is
#     println("  not counting unknown expression ", ex)
#     return undefined_syms
#   end  # end if ... elseifa

#   return undefined_syms
# end


text:
using Flux.Tracker

function loss(x, y, func)
    ŷ = func(x)
    sum(@. (y - ŷ)^2)
end

function main()
    x, y = rand(5), rand(2)
    W, b = rand(2,5), rand(2)
    affine(x) = W*x + b
    @show loss(x, y, affine)
    W = param(W)
    b = param(b)
    learning_rate=0.0001
    for i in 1:500
        ∇loss = Tracker.gradient(()->loss(x, y, affine), Params([W, b]))
        Tracker.update!(W,-learning_rate * ∇loss[W])
        Tracker.update!(b,-learning_rate * ∇loss[b])
        @show i, loss(x, y)
    end
end

main()

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function under(robj1, robj2, direction, gap=0f0)
    bb1 = robj1.boundingbox.value;
    bb2 = robj2.boundingbox.value;
    position1 = minimum(bb1)
    position2 = minimum(bb2)
    position2_upper = maximum(bb2)
    to_move = -(position2_upper-position1)
    robj2[:model] = lift(*, translationmatrix(to_move.*direction + Vec3f0(gap).*direction),  robj2[:model])
end
under_x(robj1, robj2, gap=0f0) = under(robj1, robj2, unit(Vec3f0, 1), gap)
under_y(robj1, robj2, gap=0f0) = under(robj1, robj2, unit(Vec3f0, 2), gap)
under_z(robj1, robj2, gap=0f0) = under(robj1, robj2, unit(Vec3f0, 3), gap)

function align_left(robj1, robj2, gap=0f0)
    bb1 = robj1.boundingbox.value;
    bb2 = robj2.boundingbox.value;
    position1 = minimum(bb1)
    position2 = minimum(bb2)
    to_move = (position2-position1)
    robj2[:model] = lift(*, translationmatrix(Vec3f0(-to_move[1],0,0)),  robj2[:model])
end


text:
using Test
using LazyDistances, Distances

@testset "LazyDistances.jl" begin
    include("pairwise.jl")
end


text:
module Cocoa

export window

using ObjectiveC

ObjectiveC.framework("AppKit")
Libdl.dlopen(Pkg.dir("ObjectiveC", "deps", "cocoa"))
include("constants.jl")

@classes NSWindow, NSApplication, NSSound, NSTimer

# Struct creation

nsrect(x, y, w, h) =
  ccall(:nsmakerect, Ptr{Nothing}, (Cfloat, Cfloat, Cfloat, Cfloat),
        x, y, w, h)

NSApp() = cglobal(:NSApp, Ptr{Nothing}) |> unsafe_load |> ObjectiveC.Object

function initapp()
  @objc begin
    [NSApplication sharedApplication]
#     [NSApp() setActivationPolicy:NSApplicationActivationPolicyRegular]
    [NSApp() setActivationPolicy:NSApplicationActivationPolicyAccessory]
  end
end

function window(title = "Julia", width = 600, height = 400)
  @objc begin
    win = [[NSWindow alloc] initWithContentRectRef:nsrect(0, 0, width, height)
                                            styleMask:(NSTitledWindowMask |
                                                       NSClosableWindowMask |
                                                       NSMiniaturizableWindowMask |
                                                       NSResizableWindowMask |
                                                       NSUnifiedTitleAndToolbarWindowMask)
                                              backing:NSBackingStoreBuffered
                                                defer:NO]
    [win setTitle:title]
    [win center]
    [win makeKeyAndOrderFront:nil]
    return win
  end
end

@class type Yielder
  @+ (Nothing) tick:timer begin
#     pop()
    yield()
  end
end

function eventloop()
  @async @objc begin
    try
      [NSTimer scheduledTimerWithTimeInterval:1/1000
                                      target:Yielder
                                    selector:sel"tick:"
                                    userInfo:nil
                                     repeats:YES]
      [NSApp() run]
    catch e
      showerror(STDERR, e)
      rethrow()
    end
  end
end

function init()
  initapp()
  eventloop()
end

# Sound

export play, pop

function play(name::String)
  @objc begin
    sound = [NSSound soundNamed:name]
    if [sound isPlaying] |> bool
      [sound stop]
    end
    [sound play]
  end
end

pop() = play("Pop")

# Alert window

@classes NSAlert
export NSAlert

export alert

function alert(header, body)
  @objc begin
    msg = [[NSAlert alloc] init]
    [msg setMessageText: header]
    [msg setInformativeText: body]
    [msg runModal]
    [msg release]
  end
end

end


text:
using SequentialMonteCarlo
using RNGPool
import SMCExamples.Netto.defaultNettoModel
using Test

setRNGs(0)

N = 8192
n = 10
nt = Threads.nthreads()

model, theta, ys = defaultNettoModel(n)

smcio = SMCIO{model.particle, model.pScratch}(N, n, nt, false)
smc!(model, smcio)


text:
function sesam3C_revMM(;name, k_N=60.0)
    # @parameters t [unit = uT]
    # D = Differential(t)

    # @parameters ϵ ϵ_tvr κ_E a_E [unit = uT^-1] m [unit = uT^-1] τ [unit = uT^-1]
    # @parameters k_L [unit = uT^-1] k_R [unit = uT^-1] k_mN [unit = uQ*uT^-1]
    # @parameters α_R0

    # # state variables
    # @variables B(t) [unit = uQ] L(t) [unit = uQ]  R(t) [unit = uQ] α_R(t) cumresp(t) [unit = uQ]

    # @variables α_L(t) tmpNoUnit(t)
    # # fluxes
    # @variables (begin
    #     syn_Enz(t),[unit = uQ/uT], r_M(t),[unit = uQ/uT], tvr_B(t),[unit = uQ/uT],
    #     dec_LPot(t),[unit = uQ/uT], dec_L(t),[unit = uQ/uT], dec_RPot(t),[unit = uQ/uT],
    #     dec_R(t),[unit = uQ/uT], u_C(t),[unit = uQ/uT],
    #     C_synBCt(t),[unit = uQ/uT], C_synBC(t),[unit = uQ/uT], r_G(t),[unit = uQ/uT],
    #     r_tvr(t),[unit = uQ/uT], syn_B(t),[unit = uQ/uT],
    #     r_B(t),[unit = uQ/uT], r_GEnz(t),[unit = uQ/uT], r_O(t),[unit = uQ/uT]
    # end)

    # #------------- N ---------------------
    # @variables L_N(t) [unit = uQ] R_N(t) [unit = uQ] I_N(t) [unit = uQ]
    # @variables (begin
    #     Φ_N(t),[unit = uQ/uT], Φ_Nu(t),[unit = uQ/uT], Φ_NB(t),[unit = uQ/uT],
    #     u_PlantN(t),[unit = uQ/uT], u_NOM(t),[unit = uQ/uT], u_immNPot(t),[unit = uQ/uT],
    #     u_NPot(t),[unit = uQ/uT], N_synBN(t),[unit = uQ/uT], M_ImbN(t),[unit = uQ/uT],
    #     β_NL(t), β_NR(t)
    # end)
    # @parameters β_NEnz β_NB l_N [unit = uT^-1] ν_N i_BN [unit = uT^-1]

    # # need to be specified by coupled system:
    # @variables (begin
    #     i_L(t),[unit = uQ/uT], β_Ni(t), i_IN(t),[unit = uQ/uT],
    #     u_PlantNmax(t),[unit = uQ/uT], k_PlantN(t),[unit = uT^-1]
    # end)

    # replace patttern to remove units: ,?\[unit[^\]]+\]
    @parameters t
    D = Differential(t)

    @parameters ϵ ϵ_tvr κ_E a_E  m  τ
    @parameters k_L  k_R  k_mN_L k_mN_R k_N=k_N
    @parameters α_R0

    @variables (begin
        B(t),  L(t),   R(t),  cumresp(t),
        dB(t), dL(t), dR(t), r_tot(t),
        syn_Enz(t), tvr_Enz(t), r_M(t), tvr_B(t),
        E_L(t), E_R(t),
        dec_LPot(t), dec_L(t), dec_RPot(t),
        dec_R(t), u_C(t),
        C_synBCt(t), C_synBC(t), r_G(t),
        r_tvr(t),
        r_B(t), r_GEnz(t), r_O(t),
        # need to be defined by component across all elements
        α_L(t), α_R(t),
        # need to be specified by coupled system:
        i_L(t), syn_B(t)
    end)

    eqs = [
        D(B) ~ dB, dB ~ syn_B - tvr_B,
        D(L) ~ dL, dL ~ -dec_L + i_L,
        D(R) ~ dR, dR ~ -dec_R + ϵ_tvr*tvr_B + (1-κ_E)*syn_Enz,
        syn_Enz ~ a_E*B, tvr_Enz ~ syn_Enz,
        r_M ~ m*B,
        tvr_B ~ τ*B,
        dec_LPot ~ k_L * L,
        dec_L ~ dec_LPot*(α_L * syn_Enz)/(k_mN_L + α_L*syn_Enz),
        dec_RPot ~ k_R * R,
        dec_R ~ dec_RPot*(α_R * syn_Enz)/(k_mN_R + α_R*syn_Enz),
        u_C ~ dec_L + dec_R + κ_E*syn_Enz,
        C_synBCt ~ u_C - syn_Enz/ϵ - r_M,
        #C_synBC ~ IfElse.ifelse(C_synBCt > 0.0, ϵ*C_synBCt, C_synBCt),
        C_synBC ~ C_synBCt - ((C_synBCt > 0.0) * (1-ϵ)*C_synBCt),
        #syn_B ~ C_synBC # define in all-elements model
        r_O ~ u_C - (syn_Enz/ϵ +  syn_B + r_G + r_M),
        #D(cumresp) ~ r_B + r_tvr,
        r_tvr ~ (1-ϵ_tvr)*tvr_B,
        r_B ~ r_GEnz + r_G + r_M + r_O,
        r_GEnz ~ (1-ϵ)/ϵ * syn_Enz,
        #r_G ~ IfElse.ifelse(syn_B.val > 0.0, (1-ϵ)/ϵ * syn_B, 0.0),
        r_G ~ (syn_B.val > 0.0) * (1-ϵ)/ϵ * syn_B,
        r_tot ~ r_B + r_tvr,
        E_L ~ (α_L * syn_Enz)/k_N,
        E_R ~ (α_R * syn_Enz)/k_N,
        ]
    ODESystem(eqs; name)
end

function sesam3N_revMM(;name, sC = sesam3C_revMM(name=:sC))
    @parameters t
    D = Differential(t)

    sts = @variables (begin
        L_N(t),  R_N(t),  I_N(t),
        dL_N(t),  dR_N(t),  dI_N(t),
        Φ_N(t), Φ_Nu(t), Φ_NB(t), Φ_Ntvr(t),
        u_PlantN(t), u_NOM(t), u_immNPot(t),
        u_NPot(t), N_synBN(t), M_ImbN(t),
        β_NL(t), β_NR(t),
        leach_N(t),
        # need to be specified by coupled system:
        β_Ni(t), i_IN(t),
        u_PlantNmax(t), k_PlantN(t)
    end)
    ps = @parameters β_NEnz β_NB l_N  ν_N i_BN

    @unpack L, R, dec_L, dec_R, i_L, ϵ_tvr, tvr_B, syn_B, syn_Enz, tvr_Enz, r_tvr, κ_E = sC

    eqs = [
        β_NL ~ L/L_N, β_NR ~ R/R_N,
        D(L_N) ~ dL_N, dL_N ~ -dec_L/β_NL + i_L/β_Ni,
        D(R_N) ~ dR_N, dR_N ~ -dec_R/β_NR + ϵ_tvr*tvr_B/β_NB + (1-κ_E)*tvr_Enz/β_NEnz,
        D(I_N) ~ dI_N,
        u_PlantN ~ min(u_PlantNmax, k_PlantN*I_N),
        dI_N ~ i_IN - u_PlantN - leach_N + Φ_N,
        leach_N ~ l_N*I_N,
        Φ_N ~ Φ_Nu + Φ_NB + Φ_Ntvr,
        Φ_Ntvr ~ r_tvr/β_NB,
        Φ_Nu ~ (1-ν_N) * u_NOM,
        u_NPot ~ ν_N * u_NOM + u_immNPot,
        u_NOM ~ dec_L/β_NL + dec_R/β_NR + κ_E*tvr_Enz/β_NEnz,
        u_immNPot ~ i_BN * I_N,
        N_synBN ~ u_NPot - syn_Enz/β_NEnz,
        M_ImbN ~ u_NPot - (syn_B/β_NB + syn_Enz/β_NEnz),
        Φ_NB ~ M_ImbN - u_immNPot,
        ]
    extend(ODESystem(eqs, t, sts, ps; name), sC)
end

function sesam3CN_revMM(;name, δ=40.0, max_w=12, use_seam_revenue=false, sN=sesam3N_revMM(name=:sN))
    @parameters t
    D = Differential(t)
    @unpack α_L, α_R, syn_B, B, C_synBC, β_NB, N_synBN, tvr_B, τ = sN
    sts = @variables (begin
        C_synBN(t),
        C_synBmC(t), C_synBmN(t),
        dα_L(t), dα_R(t),
        w_C(t), w_N(t), lim_C(t), lim_N(t)
    end)
    ps = @parameters δ=δ
    eqs_rev, sts_rev = use_seam_revenue ?
        get_revenue_eq_seam(sN) : get_revenue_eq_sesam3CN(sN)
    @variables α_LT(t) α_RT(t)
    eqs = [
        C_synBN ~ β_NB*N_synBN,
        syn_B ~ min(C_synBC, C_synBN),
        C_synBmC ~ min(C_synBN),
        C_synBmN ~ min(C_synBC),
        # need minimum, otherwise danger of Inf and nan -> instability
        w_C ~ exp(min(max_w, -δ/tvr_B*(C_synBC - syn_B))),
        w_N ~ exp(min(max_w, -δ/tvr_B*(C_synBN - syn_B))),
        # w_C ~ min(max_w, exp(δ/tvr_B*(C_synBmC - C_synBC))),
        # w_N ~ min(max_w, exp(δ/tvr_B*(C_synBmN - C_synBN))),
        # w_C ~ exp(δ/tvr_B*(C_synBmC - C_synBC)),
        # w_N ~ exp(δ/tvr_B*(C_synBmN - C_synBN)),
        lim_C ~ w_C/(w_C + w_N), lim_N ~ w_N/(w_C + w_N), # normalized for plot
        # α_LT, α_RT by get_revenue_eq_X
        D(α_L) ~ dα_L, dα_L ~ (α_LT - α_L)*(τ + abs(syn_B)/B),
        D(α_R) ~ dα_R, dα_R ~ (α_RT - α_R)*(τ + abs(syn_B)/B),
        ]
    extend(ODESystem(vcat(eqs,eqs_rev), t, vcat(sts, sts_rev), ps; name), sN)
end

sesam3_revMM(args...;kwargs...) = sesam3CN_revMM(args...;kwargs...)

function get_revenue_eq_seam(sN)
    # needs dec_RPot, which is specific to revMM decomposition formulation
    @parameters t
    @unpack dec_LPot, dec_RPot, k_mN_L, k_mN_R, syn_Enz, α_L, α_R, β_NL, β_NR = sN
    sts = @variables (begin
        α_LT(t), α_RT(t),
        rev_LC(t), rev_RC(t), rev_LN(t), rev_RN(t),
        α_RC(t), α_RN(t)
    end)
    # need to be defined in coupled component:
    @variables w_C(t), w_N(t)
    eqs = [
        rev_LC ~ dec_LPot/(k_mN_L + α_L*syn_Enz),
        rev_RC ~ dec_RPot/(k_mN_R + α_R*syn_Enz),
        rev_LN ~ rev_LC/β_NL,
        rev_RN ~ rev_RC/β_NR,
        α_RC ~ rev_RC/(rev_RC + rev_LC),
        α_RN ~ rev_RN/(rev_RN + rev_LN),
        α_RT ~ (w_C * α_RC + w_N * α_RN)/(w_C + w_N),
        α_LT ~ 1.0 - α_RT,
    ]
    (;eqs, sts)
end


text:
using Base.Test, DiffEqApproxFun

tic()
@time @testset "Direct ApproxFun Test" begin include("direct_approxfun.jl") end
@time @testset "Indirect ApproxFun Test" begin include("indirect_approxfun.jl") end
toc()


text:
## This file contains direct Julia implementations of depwthwise convolutions

"""
    depthwiseconv_direct!(y, x, w, cdims; alpha=1, beta=0)

Direct depthwise convolution implementation; used for debugging, tests, and mixing/
matching of strange datatypes within a single convolution.  Uses naive nested for loop
implementation and does not attempt to optimize performance.  Rather, this implementation
is intended to be maximally understandable and debuggable, to aid in testing other, more
performant implementations.  We also explicitly support mixing and matching of strange
datatypes, so that if the user really wants to convolve an image of `UInt8`'s with a
`Float16` kernel, storing the result in a `Float32` output, there is at least a function
call for that madness.

One subtlety about depthwise convolutions; the shape of a depthwise convolutional kernel
is `(spatial_dims..., C_mult, C_in)`, so the axis that must match with the number of
channels in `x` is the last, not the second-to-last, as in a normal dense convolution.

See the docstring for `conv_direct!()` for more on the optional parameters.
"""
function depthwiseconv_direct!(y::AbstractArray{yT,5}, x::AbstractArray{xT,5},
                      w::AbstractArray{wT,5}, cdims::DepthwiseConvDims;
                      alpha::yT=yT(1), beta=false) where {yT, xT, wT}
    check_dims(size(x), size(w), size(y), cdims)

    width, height, depth = input_size(cdims)
    kernel_w, kernel_h, kernel_d = kernel_size(cdims)
    out_c = channels_out(cdims)
    pad_w_lo, pad_w_hi, pad_h_lo, pad_h_hi, pad_d_lo, pad_d_hi = padding(cdims)
    dil_w, dil_h, dil_d = dilation(cdims)
    stride_w, stride_h, stride_d = stride(cdims)
    out_width, out_height, out_depth = output_size(cdims)

    # Create a method that, at compile-time, determines how we're going to index into `w`
    kproj(k, M, cdims::DepthwiseConvDims{N,K,C_mult,C_in,S,P,D,true}) where {N, K, C_mult, C_in, S, P, D} = k
    kproj(k, M, cdims::DepthwiseConvDims{N,K,C_mult,C_in,S,P,D,false}) where {N, K, C_mult, C_in, S, P, D} = M - k + 1

    # A helper function to project from output (w, h) to input (input_w, input_h)
    project(idx, stride, pad) = (idx - 1)*stride - pad + 1

    # Use `calc_padding_regions` to determine where we do or don't need to worry about padding
    padded_regions, central_region = calc_padding_regions(cdims)

    # Start with the central region
    w_region, h_region, d_region = central_region
    @inbounds for batch in 1:size(x)[end],
        c_mult in 1:channel_multiplier(cdims),
        c_in in 1:channels_in(cdims),
        d_idx in d_region,
        h_idx in h_region,
        w_idx in w_region

        # Since we're in the central region, we don't need to worry about clamping
        dotprod = yT(0)
        c_out = (c_in - 1)*channel_multiplier(cdims) + c_mult
        for kd in 1:kernel_d,
            kh in 1:kernel_h,
            kw in 1:kernel_w

            # Hoist me, you coward.
            x_d = project(d_idx, stride_d, pad_d_lo) + (kd - 1)*dil_d
            x_h = project(h_idx, stride_h, pad_h_lo) + (kh - 1)*dil_h
            x_w = project(w_idx, stride_w, pad_w_lo) + (kw - 1)*dil_w

            x_val = x[x_w, x_h, x_d, c_in, batch]
            w_val = w[kproj(kw, kernel_w, cdims),
                      kproj(kh, kernel_h, cdims),
                      kproj(kd, kernel_d, cdims),
                      c_mult, c_in]
            dotprod = muladd(x_val, w_val, dotprod)
        end
        y[w_idx, h_idx, d_idx, c_out, batch] = alpha*dotprod + beta*y[w_idx, h_idx, d_idx, c_out, batch]
    end

    # Next, do potentially-padded regions:
    @inbounds for (w_region, h_region, d_region) in padded_regions,
        batch in 1:size(x)[end],
        c_mult in 1:channel_multiplier(cdims),
        c_in in 1:channels_in(cdims),
        d_idx in d_region,
        h_idx in h_region,
        w_idx in w_region

        # Probe for out-of-bounds accesses on `x` and `continue` if we hit one
        dotprod = yT(0)
        c_out = (c_in - 1)*channel_multiplier(cdims) + c_mult
        for c_in in 1:channels_in(cdims),
            kd in 1:kernel_d

            x_d = project(d_idx, stride_d, pad_d_lo) + (kd - 1)*dil_d
            if x_d <= 0 || x_d > depth
                continue
            end

            for kh in 1:kernel_h
                x_h = project(h_idx, stride_h, pad_h_lo) + (kh - 1)*dil_h
                if x_h <= 0 || x_h > height
                    continue
                end

                for kw in 1:kernel_w
                    x_w = project(w_idx, stride_w, pad_w_lo) + (kw - 1)*dil_w
                    if x_w <= 0 || x_w > width
                        continue
                    end

                    x_val = x[x_w, x_h, x_d, c_in, batch]
                    w_val = w[kproj(kw, kernel_w, cdims),
                              kproj(kh, kernel_h, cdims),
                              kproj(kd, kernel_d, cdims),
                              c_mult, c_in]
                    dotprod = muladd(x_val, w_val, dotprod)
                end
            end
        end

        y[w_idx, h_idx, d_idx, c_out, batch] = alpha*dotprod + beta*y[w_idx, h_idx, d_idx, c_out, batch]
    end

    return y
end

"""
    ∇depthwiseconv_data_direct!(dx, dy, w, cdims; alpha=1, beta=0)

Calculate the gradient imposed upon `x` in the depthwise convolution `y = x * w`.
We make use of the fact that a depthwise convolution is equivalent to `C_in` separate
normal convolutions between that channel of `x` and the `C_mult` different kernels that
get applied to it.  The output of such a convolution is the gradient imposed upon that
particular channel of `x`, and so we simply walk through `x`, calculating the gradient
for each batch and channel independently.
"""
∇depthwiseconv_data_direct!

function ∇depthwiseconv_data_direct!(
                dx::AbstractArray{xT,5}, dy::AbstractArray{yT,5},
                w::AbstractArray{wT,5}, cdims::DepthwiseConvDims;
                alpha::xT=xT(1), beta=false) where {xT, yT, wT}
    # We do a separate convolution for each channel in x
    @inbounds for cidx in 1:channels_in(cdims)
        # For this batch and in-channel, we have a normal transposed convolution
        # between this slice of `x` and the corresponding slices of `w` and `dy`:
        dx_slice = view(dx, :, :, :, cidx:cidx, :)
        C_mult = channel_multiplier(cdims)
        dy_slice = view(dy, :, :, :, ((cidx-1)*C_mult + 1):cidx*C_mult, :)
        w_slice = permutedims(view(w, :, :, :, :, cidx:cidx), (1, 2, 3, 5, 4))

        # Adapt a DenseConvDims out of this DepthwiseConvDims, setting the in/out
        # channels appropriately for this one convolution.
        cdims_slice = DenseConvDims(cdims;
            C_in=1,
            C_out=channel_multiplier(cdims),
        )

        ∇conv_data_direct!(dx_slice, dy_slice, w_slice, cdims_slice;
                                               alpha=alpha, beta=beta)
    end
    return dx
end

"""
    ∇depthwiseconv_filter_direct!(dw, x, dy, cdims; alpha=1, beta=0)

Calculate the gradient imposed upon `w` in the depthwise convolution `y = x * w`.
"""
∇depthwiseconv_filter_direct!

function ∇depthwiseconv_filter_direct!(
                dw::AbstractArray{wT,5}, x::AbstractArray{xT,5},
                dy::AbstractArray{yT,5}, cdims::DepthwiseConvDims;
                alpha::wT=wT(1),beta=false) where {xT, yT, wT}
    # We do a separate convolution for each channel in x
    @inbounds for cidx in 1:channels_in(cdims)
        # For this batch and in-channel, we have a normal transposed convolution
        # between this slice of `x` and the corresponding slices of `w` and `dy`:
        x_slice = view(x, :, :, :, cidx:cidx, :)
        C_mult = channel_multiplier(cdims)
        dy_slice = view(dy, :, :, :, ((cidx-1)*C_mult + 1):cidx*C_mult, :)
        dw_slice = permutedims(view(dw, :, :, :, :, cidx:cidx), (1, 2, 3, 5, 4))

        # Adapt a DenseConvDims out of this DepthwiseConvDims, setting the in/out
        # channels appropriately for this one convolution.
        cdims_slice = DenseConvDims(cdims;
            C_in=1,
            C_out=channel_multiplier(cdims),
        )

        ∇conv_filter_direct!(dw_slice, x_slice, dy_slice, cdims_slice;
                                                alpha=alpha, beta=beta)
        dw[:, :, :, :, cidx:cidx] .= permutedims(dw_slice, (1, 2, 3, 5, 4))
    end
    return dw
end


text:
!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
using DSP, Test, FilterTestHelpers, Polynomials
using Polynomials.PolyCompat

@testset "convert to SOS" begin
    # Test conversion to SOS against MATLAB
    # Test poles/zeros generated with:
    #=
    mag = rand(10)
    arg = rand(10)
    z = mag.*cospi(arg) + im*mag.*sinpi(arg)
    z = [z, conj(z)]
    mag = rand(10)
    arg = rand(10)
    p = mag.*cospi(arg) + im*mag.*sinpi(arg)
    p = [p, conj(p)]
    k = real(prod(p))/real(prod(z))
    =#
    z = [0.07674942314081117 + 0.5605934468331276im,
        -0.10637764563083413 + 0.00938830970435945im,
        0.19723688182100613 + 0.20500254264692958im,
        0.07019769777809694 + 0.00040337356287483566im,
        -0.029225401438832663 + 0.35500551559734517im,
        -0.35980697033738923 + 0.21595798755003964im,
        0.3712083544916094 + 0.41767928564585416im,
        -0.17585090520154495 + 0.32300027988596314im,
        -0.2308322620393009 + 0.3539248310694154im,
        -0.008498685920569043 + 0.028356975487226484im,
        0.07674942314081117 - 0.5605934468331276im,
        -0.10637764563083413 - 0.00938830970435945im,
        0.19723688182100613 - 0.20500254264692958im,
        0.07019769777809694 - 0.00040337356287483566im,
        -0.029225401438832663 - 0.35500551559734517im,
        -0.35980697033738923 - 0.21595798755003964im,
        0.3712083544916094 - 0.41767928564585416im,
        -0.17585090520154495 - 0.32300027988596314im,
        -0.2308322620393009 - 0.3539248310694154im,
        -0.008498685920569043 - 0.028356975487226484im]
    p = [-0.946179900464128 + 0.23187351222922784im,
        0.05191136844411161 + 0.8713748123758278im,
        -0.05109307872385063 + 0.7440280322767342im,
        -0.032892467440199066 + 0.024218620496805687im,
        -0.1266287653888985 + 0.33150956246323654im,
        -0.15051989876024102 + 0.9373469058653078im,
        -0.6672740143157547 + 0.5034138963468052im,
        -0.790336466064852 + 0.109069102200402im,
        -0.009317017740249631 + 0.013158400271660778im,
        0.5075897927503011 + 0.02902816288546107im,
        -0.946179900464128 - 0.23187351222922784im,
        0.05191136844411161 - 0.8713748123758278im,
        -0.05109307872385063 - 0.7440280322767342im,
        -0.032892467440199066 - 0.024218620496805687im,
        -0.1266287653888985 - 0.33150956246323654im,
        -0.15051989876024102 - 0.9373469058653078im,
        -0.6672740143157547 - 0.5034138963468052im,
        -0.790336466064852 - 0.109069102200402im,
        -0.009317017740249631 - 0.013158400271660778im,
        0.5075897927503011 - 0.02902816288546107im]
    k = 10586.686805294861

    # MATLAB:
    #=
    [sos, g] = zp2sos(z, p, k)
    save('...', 'sos', '-ascii', '-double')
    =#
    m_sos_full = [
    1.0000000000000000e+00  -1.4039539555619387e-01   4.9278794835762620e-03   1.0000000000000000e+00   1.8634035480499262e-02   2.5995031728136877e-04
    1.0000000000000000e+00   1.6997371841138087e-02   8.7634572115964215e-04   1.0000000000000000e+00   6.5784934880398133e-02   1.6684559930728524e-03
    1.0000000000000000e+00   2.1275529126166826e-01   1.1404343849064294e-02   1.0000000000000000e+00   2.5325753077779700e-01   1.2593343422848324e-01
    1.0000000000000000e+00  -3.9447376364201225e-01   8.0928430042179728e-02   1.0000000000000000e+00  -1.0151795855006023e+00   2.5849003194479853e-01
    1.0000000000000000e+00  -7.4241670898321876e-01   3.1225162810199936e-01   1.0000000000000000e+00   1.0218615744770126e-01   5.5618821550707065e-01
    1.0000000000000000e+00   3.5170181040308990e-01   1.3525272166661328e-01   1.0000000000000000e+00   1.5806729321297039e+00   6.3652779864668074e-01
    1.0000000000000000e+00   4.6166452407860181e-01   1.7854631924569470e-01   1.0000000000000000e+00   1.3345480286315095e+00   6.9868016121613397e-01
    1.0000000000000000e+00   5.8450802877665325e-02   1.2688304019379779e-01   1.0000000000000000e+00  -1.0382273688822322e-01   7.6198885381674941e-01
    1.0000000000000000e+00  -1.5349884628162233e-01   3.2015548658469395e-01   1.0000000000000000e+00   3.0103979752048204e-01   9.0127546185805940e-01
    1.0000000000000000e+00   7.1961394067477846e-01   1.7609890829003397e-01   1.0000000000000000e+00   1.8923598009282561e+00   9.4902172971582510e-01
    ]
    @test m_sos_full ≈ sosfilter_to_matrix(convert(SecondOrderSections, ZeroPoleGain(z, p, k)))

    # And with half of the zeros removed
    # MATLAB:
    # [sos, g] = zp2sos(z([1:5, 11:15]), p, k)
    zp = z[[1:5; 11:15]]
    m_sos_half = [
    0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00   1.0000000000000000e+00   1.8634035480499262e-02   2.5995031728136877e-04
    0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00   1.0000000000000000e+00   6.5784934880398133e-02   1.6684559930728524e-03
    0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00   1.0000000000000000e+00   2.5325753077779700e-01   1.2593343422848324e-01
    0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00   1.0000000000000000e+00  -1.0151795855006023e+00   2.5849003194479853e-01
    0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00   1.0000000000000000e+00   1.0218615744770126e-01   5.5618821550707065e-01
    1.0000000000000000e+00  -3.9447376364201225e-01   8.0928430042179728e-02   1.0000000000000000e+00   1.5806729321297039e+00   6.3652779864668074e-01
    1.0000000000000000e+00  -1.4039539555619387e-01   4.9278794835762620e-03   1.0000000000000000e+00   1.3345480286315095e+00   6.9868016121613397e-01
    1.0000000000000000e+00   5.8450802877665325e-02   1.2688304019379779e-01   1.0000000000000000e+00  -1.0382273688822322e-01   7.6198885381674941e-01
    1.0000000000000000e+00  -1.5349884628162233e-01   3.2015548658469395e-01   1.0000000000000000e+00   3.0103979752048204e-01   9.0127546185805940e-01
    1.0000000000000000e+00   2.1275529126166826e-01   1.1404343849064294e-02   1.0000000000000000e+00   1.8923598009282561e+00   9.4902172971582510e-01
    ]
    @test m_sos_half ≈ sosfilter_to_matrix(convert(SecondOrderSections, ZeroPoleGain(zp, p, k)))

    # And with an extra real pole
    pp = [p; 0.7]
    # MATLAB:
    #=
    pp = [p'; 0.7]
    [sos, g] = zp2sos(z, pp, k)
    =#
    m_sos_extra = [
    0.0000000000000000e+00   1.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00  -6.9999999999999996e-01   0.0000000000000000e+00
    1.0000000000000000e+00  -1.4039539555619387e-01   4.9278794835762620e-03   1.0000000000000000e+00   1.8634035480499262e-02   2.5995031728136877e-04
    1.0000000000000000e+00   1.6997371841138087e-02   8.7634572115964215e-04   1.0000000000000000e+00   6.5784934880398133e-02   1.6684559930728524e-03
    1.0000000000000000e+00   2.1275529126166826e-01   1.1404343849064294e-02   1.0000000000000000e+00   2.5325753077779700e-01   1.2593343422848324e-01
    1.0000000000000000e+00  -3.9447376364201225e-01   8.0928430042179728e-02   1.0000000000000000e+00  -1.0151795855006023e+00   2.5849003194479853e-01
    1.0000000000000000e+00  -7.4241670898321876e-01   3.1225162810199936e-01   1.0000000000000000e+00   1.0218615744770126e-01   5.5618821550707065e-01
    1.0000000000000000e+00   3.5170181040308990e-01   1.3525272166661328e-01   1.0000000000000000e+00   1.5806729321297039e+00   6.3652779864668074e-01
    1.0000000000000000e+00   4.6166452407860181e-01   1.7854631924569470e-01   1.0000000000000000e+00   1.3345480286315095e+00   6.9868016121613397e-01
    1.0000000000000000e+00   5.8450802877665325e-02   1.2688304019379779e-01   1.0000000000000000e+00  -1.0382273688822322e-01   7.6198885381674941e-01
    1.0000000000000000e+00  -1.5349884628162233e-01   3.2015548658469395e-01   1.0000000000000000e+00   3.0103979752048204e-01   9.0127546185805940e-01
    1.0000000000000000e+00   7.1961394067477846e-01   1.7609890829003397e-01   1.0000000000000000e+00   1.8923598009282561e+00   9.4902172971582510e-01
    ]
    @test m_sos_extra ≈ sosfilter_to_matrix(convert(SecondOrderSections, ZeroPoleGain(z, pp, k)))

    # And with only poles (no zeros)
    m_sos_only_poles = copy(m_sos_full)
    m_sos_only_poles[:, 1:2] .= 0
    m_sos_only_poles[:, 3] .= 1
    @test m_sos_only_poles ≈ sosfilter_to_matrix(convert(SecondOrderSections, ZeroPoleGain(Float64[], p, k)))

    # Test that a filter with repeated zeros is handled properly
    # MATLAB:
    #=
    [z,p,k] = butter(2, [49.5 50.5]/500, 'stop')
    [sos,g] = zp2sos(z, p, k)
    =#
    m_sos_butterworth_bs = [
    1.0000000000000000e+00  -1.9021224191804869e+00   1.0000000000000000e+00   1.0000000000000000e+00  -1.8964983429993663e+00   9.9553672990017417e-01
    1.0000000000000000e+00  -1.9021224191804869e+00   1.0000000000000000e+00   1.0000000000000000e+00  -1.8992956433548462e+00   9.9559721515078736e-01
    ]
    f = convert(SecondOrderSections, digitalfilter(Bandstop(49.5, 50.5; fs=1000), DSP.Butterworth(2)))
    @test m_sos_butterworth_bs ≈ sosfilter_to_matrix(f)
    @test f.g ≈ 0.995566972017647

    # Test that a numerically challenging filter (high order, clustered
    # roots) has acceptable errors in its coefficients after conversion to
    # SOS
    f = ZeroPoleGain(ones(100), 0.99*ones(100), 1)
    g = convert(SecondOrderSections, f)
    tffilter_eq(convert(PolynomialRatio, f), convert(PolynomialRatio, g))

    H = ZeroPoleGain{:z}([1+im, 1-im, 0.5+im, 0.5-im], [1, 0.0, 0.0, 0.0], 1.0)
    H′ = ZeroPoleGain(SecondOrderSections(H))
    @test sort(H.p, by=z->(real(z), imag(z))) ≈ sort(H′.p, by=z->(real(z), imag(z)))
    @test sort(H.z, by=z->(real(z), imag(z))) ≈ sort(H′.z, by=z->(real(z), imag(z)))
    @test H.k ≈ H′.k
end

@testset "conversions" begin
    for f in (digitalfilter(Lowpass(0.5), Butterworth(1)), digitalfilter(Lowpass(0.5), Butterworth(2)),
              digitalfilter(Bandpass(0.25, 0.75), Butterworth(1)))
        for ftype1 in (ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections)
            f2 = convert(ftype1, f)
            for ftype2 in (ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections)
                f3 = convert(ftype2, f)
                try
                    zpkfilter_eq(convert(ZeroPoleGain, f), convert(ZeroPoleGain, f3), sqrt(eps()))
                catch e
                    println("Conversion from $ftype1 to $ftype2 failed:")
                    rethrow(e)
                end
            end
        end
    end

    for proto in (Butterworth(3), Chebyshev1(3, 1), Chebyshev2(3, 1))
        f = digitalfilter(Lowpass(0.5), proto)
        for ftype1 in (ZeroPoleGain, PolynomialRatio, SecondOrderSections)
            f2 = convert(ftype1, f)
            for ftype2 in (ZeroPoleGain, PolynomialRatio, SecondOrderSections)
                f3 = convert(ftype2, f2)
                zpkfilter_eq(convert(ZeroPoleGain, f), convert(ZeroPoleGain, f3), 2e-5)
            end
        end
    end
end

@testset "filter gain and composition" begin
    # Test setting filter gain
    x = randn(100)
    f1 = digitalfilter(Lowpass(0.3), Butterworth(2))
    y = filt(f1, x)
    for ty in (ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections)
        @test filt(3*convert(ty, f1), x) ≈ 3*y
        @test filt(convert(ty, f1)*3, x) ≈ 3*y
    end

    # Test composing filters
    f2 = digitalfilter(Highpass(0.5), Butterworth(1))
    y = filt(f2, y)
    for ty in (ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections)
        @test filt(convert(ty, f1)*convert(ty, f2), x) ≈ y
    end

    f3 = digitalfilter(Bandstop(0.35, 0.4), Butterworth(1))
    y = filt(f3, y)
    for ty in (ZeroPoleGain, PolynomialRatio, Biquad, SecondOrderSections)
        @test filt(convert(ty, f1)*convert(ty, f2)*convert(ty, f3), x) ≈ y
    end
    @test filt(convert(Biquad, f1)*(convert(Biquad, f2)*convert(Biquad, f3)), x) ≈ y
    @test filt((convert(Biquad, f1)*convert(Biquad, f2))*convert(Biquad, f3), x) ≈ y
end

@testset "filter inversion and exponentiation" begin
    # A filter is stable if all poles within unit circle and causal if not more zeros
    # than poles. For the inverse to also be stable and causal, it follows that the
    # zeros also have to be within the unit circle and the number of poles and zeros
    # must be equal. But poles and zeros at zero lead to degenerate polynomials in z⁻¹,
    # so ensure these are tested as well. Further, we want both complex conjugate and
    # real poles and zeros.
    for Npr ∈ 0:2, Npc ∈ 0:2, Nzr ∈ 0:2, Nzc ∈ 0:2
        z = rand(ComplexF64, Nzc) .- (0.5 + 0.5im);
        z = [z; conj(z); rand(Nzr) .- 0.5; zeros(max(2Npc+Npr-2Nzc-Nzr, 0))]
        p = rand(ComplexF64, Npc) .- (0.5 + 0.5im);
        p = [p; conj(p); rand(Npr) .- 0.5; zeros(max(2Nzc+Nzr-2Npc-Npr, 0))]
        H′ = ZeroPoleGain(z, p,
            (rand() + 0.5) * rand([-1, 1]), # non-zero gain with random sign
        )
        maybe_biquad = length(z) ≤ 2 && length(p) ≤ 2 ? (Biquad,) : ()
        for T ∈ (PolynomialRatio, ZeroPoleGain, SecondOrderSections, maybe_biquad...)
            H = T(H′)
            H⁻¹ = inv(H)
            x = rand(100)
            @test filt(H⁻¹, filt(H, x)) ≈ x
            for p ∈ 1:3
                Hᵖ = H^p
                H⁻ᵖ = H^-p
                @test filt(Hᵖ, x) ≈ reduce(∘, fill(x -> filt(H, x), p))(x)
                @test filt(H⁻ᵖ, filt(Hᵖ, x)) ≈ x rtol=5e-8
            end
            @test filt(H^0, x) ≈ x
        end
    end
end

@testset "types" begin
    # normalizes and promotes to result type of division
    @test @inferred(PolynomialRatio{:z}([1, 2], [3, 4])) isa PolynomialRatio{:z,Float64}
    # does not normalize, keeps type
    @test @inferred(PolynomialRatio{:z,Int}([1, 2], [1, 4])) isa PolynomialRatio{:z,Int}
    # normalizes, keeps type
    @test @inferred(PolynomialRatio{:z,Int}([2, 4], [2, 8])) isa PolynomialRatio{:z,Int}
    # throws because normalization is impossible within Int
    @test_throws InexactError PolynomialRatio{:z,Int}([1, 2], [3, 4])
    # throws because normalization is impossible for a0 = 0
    @test_throws ArgumentError PolynomialRatio{:z,Float64}([1.0, 2.0], [0.0, 4.0])
    # does not normalize, uses type from input
    @test @inferred(PolynomialRatio{:s}([1, 2], [3, 4])) isa PolynomialRatio{:s,Int}
    # throws because denominator must not be zero
    @test_throws ArgumentError PolynomialRatio{:s}([1.0, 2.0], [0.0])
    @test_throws ArgumentError PolynomialRatio{:s}([1.0, 2.0], Float64[])
end

@testset "misc" begin
    f1 = digitalfilter(Lowpass(0.3), Butterworth(2))
    f2 = digitalfilter(Highpass(0.5), Butterworth(1))
    # Test some otherwise untested code paths
    @test promote_type(ZeroPoleGain{:z,ComplexF32,ComplexF64,Float32}, ZeroPoleGain{:z,ComplexF64,ComplexF32,Float64}) == ZeroPoleGain{:z,ComplexF64,ComplexF64,Float64}
    @test convert(ZeroPoleGain{:z,ComplexF64,ComplexF64,Float64}, f1) === f1
    f1f = convert(ZeroPoleGain{:z,ComplexF32,ComplexF32,Float32}, f1)
    @test f1f.z == convert(Vector{ComplexF32}, f1.z)
    @test f1f.p == convert(Vector{ComplexF32}, f1.p)
    @test f1f.k == convert(Float32, f1.k)

    @test_throws ArgumentError PolynomialRatio(Float64[], Float64[])
    @test promote_type(PolynomialRatio{:z,Float32}, PolynomialRatio{:z,Int}) == PolynomialRatio{:z,Float32}
    f1f = convert(PolynomialRatio{:z,Float32}, f1)
    f1p = convert(PolynomialRatio, f1)
    @test coefb(f1f) == convert(Array{Float32}, coefb(f1p))
    @test coefa(f1f) == convert(Array{Float32}, coefa(f1p))

    b = Biquad(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 2)
    @test b.b0 === 0.0
    @test b.b1 === 2*1/3
    @test b.b2 === 2*2/3
    @test b.a1 === 4/3
    @test b.a2 === 5/3
    @test promote_type(Biquad{:z,Float32}, Biquad{:z,Int}) == Biquad{:z,Float32}
    @test convert(Biquad{:z,Float32}, b) == Biquad{:z,Float32}(0.0, 2*1/3, 2*2/3, 4/3, 5/3)
    zpkfilter_eq(convert(ZeroPoleGain{:z,ComplexF64,ComplexF64,Float64}, b), convert(ZeroPoleGain, b), 0.0)
    f = convert(PolynomialRatio, Biquad(2.0, 0.0, 0.0, 0.0, 0.0))
    @test coefb(f) == [2.0]
    @test coefa(f) == [1.0]
    @test convert(Biquad, PolynomialRatio([4.0], [2.0])) == Biquad(2.0, 0.0, 0.0, 0.0, 0.0)
    @test Biquad(2.0, 0.0, 0.0, 0.0, 0.0)*2 == Biquad(4.0, 0.0, 0.0, 0.0, 0.0)
    @test convert(Biquad{:z,Float64}, f1) == convert(Biquad, f1)
    f = PolynomialRatio(Float64[1.0], Float64[1.0])

    @test_throws ArgumentError convert(SecondOrderSections, ZeroPoleGain([0.5 + 0.5im, 0.5 + 0.5im], [0.5 + 0.5im, 0.5 - 0.5im], 1))
    @test_throws ArgumentError convert(SecondOrderSections, ZeroPoleGain([0.5 + 0.5im, 0.5 - 0.5im], [0.5 + 0.5im, 0.5 + 0.5im], 1))
    @test_throws ArgumentError convert(SecondOrderSections, ZeroPoleGain([1+im, 1+im, 1-im], [1, 0, 0], 1))
    @test_throws ArgumentError convert(SecondOrderSections, ZeroPoleGain([1+im, 1-im, 1-im], [1, 0, 0], 1))

    @test promote_type(SecondOrderSections{:z,Float64,Float32}, SecondOrderSections{:z,Float32,Float64}) == SecondOrderSections{:z,Float64,Float64}
    @test convert(SecondOrderSections{:z,Float32,Float32}, convert(SecondOrderSections, b)).biquads == convert(SecondOrderSections, convert(Biquad{:z,Float32}, b)).biquads
    @test convert(Biquad, convert(SecondOrderSections, b)) == b
    @test_throws ArgumentError convert(Biquad, convert(SecondOrderSections, f1*f2))
    f1p = convert(PolynomialRatio{:z,Float32}, convert(PolynomialRatio, convert(SecondOrderSections, f1)))
    f1f = convert(PolynomialRatio{:z,Float32}, convert(SecondOrderSections, f1))
    @test f1p.b == f1f.b
    @test f1p.a == f1f.a
end

@testset "coef() testing" begin
    # relies on conversion for non-polynomial ratio filter coeff objects
    # so no need to test conversion correctness

    @testset "Biquad" begin
        bs = [1, 2, 3, 4, 5]
        B = Biquad(bs...)
        @test coefa(B) == [1, 4, 5]
        @test coefb(B) == [1, 2, 3]

        bs = Float64[20, 16, 13, 31, 33]
        B = Biquad(bs...)
        @test coefa(B) == [1.0, 31, 33]
        @test coefb(B) == [20.0, 16, 13]
    end

    @testset "SecondOrderSections" begin
        bs = [2.0, 0, 0, 0, 0]
        B = SecondOrderSections(repeat([Biquad(bs...)], 2), 0.25)
        @test coefb(B) == [1.0]
        @test coefa(B) == [1.0]

        bs = [0, 1, 0, 0, 0]
        B = SecondOrderSections(repeat([Biquad(bs...)], 2), 1)
        @test coefb(B) == [0, 0, 1]
        @test coefa(B) == [1]
    end

    @testset "ZeroPoleGain" begin
        f = ZeroPoleGain([0], [-1, 1], 1)
        @test coefa(f) == [1, 0, -1]
        @test coefb(f) == [0, 1]

        f = ZeroPoleGain(Int[], [-0.25, 0.25], 1)
        @test coefa(f) == [1.0, 0, -1/16]
        @test coefb(f) == [0.0, 0.0, 1.0]
    end
end

@testset "issue #432" begin
    bp1 = 0.75
    bp2 = 10.0
    Fsamp = 180
    responsetype = Bandpass(bp1, bp2; fs = Fsamp)
    designmethod = Elliptic(11, 0.25, 40)
    H = digitalfilter(responsetype, designmethod)
    H′ = ZeroPoleGain(SecondOrderSections(H))
    @test sort(H.p, by=z->(real(z), imag(z))) ≈ sort(H′.p, by=z->(real(z), imag(z)))
    @test sort(H.z, by=z->(real(z), imag(z))) ≈ sort(H′.z, by=z->(real(z), imag(z)))
    @test H.k ≈ H′.k
end


text:
# module ImmerseCompose

#import Base: start, next, done

# import Compose
using Measures: Vec, Measure, AbsoluteBox, resolve
using Compose: Backend
using Compose: Container, Context, Table
using Compose: Form, Line, Circle
using Compose: Property, Stroke, Fill, LineWidth, Visible, SVGClass
using Compose: List, ListNode, ListNull
using Compose: Transform, IdentityTransform, MatrixTransform, UnitBox
using Compose: absolute_native_units

start(l::List) = l
next(::List, l::List) = (l.head, l.tail)
done(::List, l::List) = typeof(l) <: ListNull
cons(value, l::List{T}) where T = ListNode{T}(value, l)

# using Compat, Colors #, GtkUtilities
# import Gtk

# export
#     find_object,
#     find_tagged,
#     find_panelforms,
#     find_inmask,
#     bareobj,
#     getproperty,
#     setproperty!,
#     nearest,
#     hitcenter,
#     absolute_to_data

const ContainersWithChildren = Table
const Iterables = Union{ContainersWithChildren, AbstractArray}

iterable(cnt::ContainersWithChildren) = cnt.children
iterable(a::AbstractArray) = a

# Testing utilities
#
# function test_by_drawing(c, obj, transform, units, box)
#     backend = Image{CairoBackend}(Main.Gtk.cairo_surface(c))
#     form = transform_form(transform, units, box, obj)
#     backend.stroke = RGBA(1.0,0,0,1.0)
#     Compose.draw(backend, form.primitives[1])
# end

# function transform_form(t::Transform, units::UnitBox,
#                         box::AbsoluteBox, form::Form)
#     Form([Compose.absolute_units(primitive, t, units, box)
#                                for primitive in form.primitives])
# end

# For now only Forms have tags, but this may change

# Find a tagged object
function find_object(cnt::Iterables, tag)
    for item in iterable(cnt)
        if has_tag(item, tag)
            return item
        end
        ret = find_object(item, tag)
        if ret != nothing
            return ret
        end
    end
    nothing
end

function find_object(ctx::Context, tag)
    for item in ctx.form_children
        if has_tag(item, tag)
            return item
        end
    end
    for item in ctx.container_children
        ret = find_object(item, tag)
        if ret != nothing
            return ret
        end
    end
    nothing
end

find_object(obj, tag) = nothing

# Find a tagged object's parent Context
function find_parent(root, tag)
    found, ret = _find_parent(root, tag)
    found == false && error(tag, " not found")
    ret
end

function _find_parent(obj::Iterables, tag)
    for item in iterable(obj)
        if has_tag(item, tag)
            return true, obj
        end
        found, ret = _find_parent(item, tag)
        found && return found, ret
    end
    false, obj
end

function _find_parent(obj::Context, tag)
    for item in obj.form_children
        if has_tag(item, tag)
            return true, obj
        end
    end
    for item in obj.container_children
        found, ret = _find_parent(item, tag)
        found && return found, ret
    end
    false, obj
end

_find_parent(obj, tag) = false, obj

# Find the enclosing Context (parent) of all tagged Forms
function find_tagged(root)
    handles = Dict{Symbol,Context}()
    find_tagged!(handles, root)
end

function find_tagged!(handles, obj::Iterables)
    for item in iterable(obj)
        if has_tag(item)
            handles[item.tag] = obj
        else
            find_tagged!(handles, item)
        end
    end
    handles
end

function find_tagged!(handles, obj::Context)
    for item in obj.form_children
        if has_tag(item)
            handles[item.tag] = obj
        end
    end
    for item in obj.container_children
        find_tagged!(handles, item)
    end
    handles
end

find_tagged!(handles, obj) = obj

# Find the nested series of containers leading to a tagged object
function find_path(root::Context, tag)
    ret = Any[]
    find_path!(ret, root, tag)
    ret
end

function find_path!(ret, cnt::Iterables, tag)
    for item in iterable(cnt)
        if find_path!(ret, item, tag)
            isa(item, Array) || push!(ret, item)
            return true
        end
    end
    false
end

function find_path!(ret, ctx::Context, tag)
    for iter in (ctx.form_children, ctx.container_children)
        for item in iter
            if find_path!(ret, item, tag)
                push!(ret, item)
                return true
            end
        end
    end
    false
end

find_path!(ret, form::Form, tag) = has_tag(form, tag)

find_path!(ret, obj, tag) = false

has_tag(form::Form, tag) = form.tag == tag
has_tag(form::Form)      = form.tag != Compose.empty_tag

has_tag(obj, tag) = false
has_tag(obj)      = false

# Finding all Forms in a plotpanel
# This is for working with untagged objects
find_panelforms(cnt) = find_panelforms!(Any[], cnt, false)

function find_panelforms!(forms, ctx::Context, inpanel::Bool)
    if !inpanel
        for child in ctx.property_children
            if isa(child, SVGClass) && length(child.primitives) == 1 && child.primitives[1].value == "plotpanel"
                inpanel = true
                break
            end
        end
        for child in ctx.container_children
            find_panelforms!(forms, child, inpanel)
        end
    else
        for child in ctx.form_children
            push!(forms, child)
        end
        for child in ctx.container_children
            find_panelforms!(forms, child, true)
        end
    end
    forms
end

function find_panelforms!(forms, cnt::Iterables, inpanel::Bool)
    for child in iterable(cnt)
        find_panelforms!(forms, child, inpanel)
    end
    forms
end

find_panelforms!(forms, obj, inpanel::Bool) = forms

function bareobj(ctx::Context)
    for c in ctx.form_children
        return c
    end
    error("No bareobj found in ", string_compact(ctx))
end
bareobj(f::Form) = f

# Modifying objects

sym2proptype(sym::Symbol) =
    sym == :stroke    ? Stroke :
    sym == :fill      ? Fill :
    sym == :linewidth ? LineWidth :
    sym == :visible   ? Visible :
    error(sym, " not a recognized property")

proptype2sym(::Type{Stroke})    = :stroke
proptype2sym(::Type{Fill})      = :fill
proptype2sym(::Type{LineWidth}) = :linewidth
proptype2sym(::Type{Visible})   = :visible

getproperty(ctx::Context, sym::Symbol) = getproperty(ctx, sym2proptype(sym))

function getproperty(ctx::Context, ::Type{P}, default=nothing) where {P<:Property}
    for c in ctx.property_children
        isa(c, P) && return _getvalue(c)
    end
    if default==nothing
        error(proptype2sym(P), " not found in ", string_compact(ctx))
    end
    return default
end
getproperty(ctx::Context, ::Type{Visible}) = getproperty(ctx, Visible, true)

_getvalue(p::Stroke)    = [prim.color for prim in p.primitives]
_getvalue(p::Fill)      = [prim.color for prim in p.primitives]
_getvalue(p::LineWidth) = [prim.value for prim in p.primitives]
_getvalue(p::Visible)   = [prim.value for prim in p.primitives]

setproperty!(ctx::Context, val, sym::Symbol) = setproperty!(ctx, val, sym2proptype(sym))

setproperty!(ctx::Context, val::Union{Colorant,AbstractString,AbstractArray}, ::Type{P}) where {P<:Stroke} =
    setproperty!(ctx, Compose.stroke(val))

setproperty!(ctx::Context, val::Union{Colorant,AbstractString,AbstractArray}, ::Type{P}) where {P<:Fill} =
    setproperty!(ctx, Compose.fill(val))

setproperty!(ctx::Context, val::Union{Measure,Number}, ::Type{P}) where {P<:LineWidth} =
    setproperty!(ctx, Compose.linewidth(val))

setproperty!(ctx::Context, val::Bool, ::Type{P}) where {P<:Visible} =
    setproperty!(ctx, Compose.visible(val))

function setproperty!(ctx::Context, val::P) where {P<:Property}
    ctx.property_children = _setproperty!(ctx.property_children, val)
    ctx
end

# Substitutes or adds a new property node in the list of children
function _setproperty!(iter, val::P, i) where {P}
    if done(iter, i)
        # add a node
        return ListNode{Property}(val, i)
    end
    item, inew = next(iter, i)
    if isa(item, P)
        # replace the node
        return ListNode{Property}(val, inew)
    end
    ListNode{Property}(item, _setproperty!(iter, val, inew))
end
_setproperty!(iter, val) = _setproperty!(iter, val, start(iter))

function string_compact(obj)
    io = IOBuffer()
    showcompact(io, obj)
    takebuf_string(io)
end

# Absolute-to-relative coordinate transformations
function absolute_to_data(x, y, transform, unit_box::UnitBox, parent_box::AbsoluteBox)
    xt, yt = invert_transform(transform, x, y)
    (unit_box.x0 + unit_box.width *(xt-parent_box.x0[1])/Measures.width(parent_box),
     unit_box.y0 + unit_box.height*(yt-parent_box.x0[2])/Measures.height(parent_box))
end

invert_transform(::IdentityTransform, x, y) = x, y

function invert_transform(t::MatrixTransform, x, y)
    @assert t.M[3,1] == t.M[3,2] == 0
    xyt = t.M\[x, y, 1.0]
    xyt[1], xyt[2]
end

# Iterators/functions that report object specs in native (screen) units
# These are useful for hit-testing

abstract type FormIterator end

struct FormNativeIterator{F<:Form,B<:Backend,C} <: FormIterator
    form::F
    backend::B
    coords::C
end

start(iter::FormNativeIterator) = 1  # some will specialize this
done(iter::FormNativeIterator, state::Int) = state > length(iter.form.primitives)
done(iter::FormNativeIterator, state::Tuple{Int,Int}) = state[1] > length(iter.form.primitives)

native(form::F, backend::B, coords::C) where {F,B,C} = FormNativeIterator{F,B,C}(form, backend, coords)

@inline inc(state::Tuple{Int,Int}, len::Int) =
    ifelse(state[2]+1 > len, (state[1]+1, 1), (state[1], state[2]+1))

# Scalars and points
function native(m::Measure, backend::Backend, coords)
    box, units, transform = coords
    absolute_native_units(backend, resolve(box, units, transform, m))
end

function native(m::Tuple{Measure,Measure}, backend::Backend, coords)
    box, units, transform = coords
    absolute_native_units(backend, resolve(box, units, transform, m))
end

# Iteration for specific Forms
start(::FormNativeIterator{F}) where {F<:Line} = (1,1)
function next(iter::FormNativeIterator{F}, state) where {F<:Line}
    nxy = native(iter.form.primitives[state[1]].points[state[2]], iter.backend, iter.coords)
    nxy, inc(state, length(iter.form.primitives[state[1]].points))
end

@inline function next(iter::FormNativeIterator{F}, state) where {F<:Circle}
    prim = iter.form.primitives[state]
    nx, ny = native(prim.center, iter.backend, iter.coords)
    nr = native(prim.radius, iter.backend, iter.coords)
    (nx, ny, nr), state+1
end

# Hit testing
function nearest(backend::Backend, coords, form::Circle, x, y)
    box, units, transform = coords
    mindist = Inf
    mini = 0
    for (i,nxyr) in enumerate(native(form, backend, coords))
        nx, ny, nr = nxyr
        dist = Float64((nx-x)^2 + (ny-y)^2 - nr^2)
        if dist < mindist
            mindist = max(0.0, dist)
            mini = i
        end
    end
    sqrt(mindist), mini
end

# Here we don't use the iterator because we need to break it out by segments
# Return index is (iline, isegment, fraction) where 0.0 <= fraction <= 1.0
# measures the fraction of the distance along the isegment-th segment
# of the nearest point to the click.
function nearest(backend::Backend, coords, form::Line, x, y)
    box, units, transform = coords
    mindist = Inf
    mini = (0,0,0.0)
    for iline = 1:length(form.primitives)
        thisline = form.primitives[iline]
        length(thisline.points) < 2 && continue
        nxold, nyold = native(thisline.points[1], backend, coords)
        skipping = false
        for isegment = 1:length(thisline.points)-1
            nx, ny = native(thisline.points[isegment+1], backend, coords)
            len = Float64(sqrt((nx-nxold)^2 + (ny-nyold)^2)) # length of segment
            if len == 0
                # The two points defining the segment are the same,
                # we'll report that it's either the incoming or outgoing
                # segment unless there are no other segments
                if !skipping || length(thisline.points)>isegment+1
                    skipping = true
                    continue
                end
                skipping = false
                dist = Float64(sqrt((nx-x)^2 + (ny-y)^2))
                if dist < mindist
                    mindist = dist
                    mini = (iline, isegment, 0.0)
                end
            else
                skipping = false
                # Compute the projection of (x,y) onto the segment
                vx, vy = (nx-nxold)/len, (ny-nyold)/len  # unit vector tanget
                l = (x-nxold)*vx + (y-nyold)*vy
                frac = l/len
                px, py = l*vx+nxold, l*vy+nyold
                nxold, nyold = nx, ny
                0.0 <= frac <= 1.0 || continue
                # Compute the distance
                dist = Float64(sqrt((px-x)^2 + (py-y)^2))
                if dist < mindist
                    mindist = dist
                    mini = (iline, isegment, frac)
                end
            end

        end
    end
    mindist, mini
end

function hitcenter(backend::Backend, coords, form::Line, index)
    prim = form.primitives[index[1]]
    frac = index[3]
    idx = index[2] + (frac >= 0.5 ? 1 : 0)
    pt = prim.points[idx]
    native(pt, backend, coords)
end

function hitcenter(backend, coords, form::Circle, index)
    circ = form.primitives[index]
    native(circ.center, backend, coords)
end

# selection by mask
function find_inmask(backend::Backend, coords, form::Circle, mask)
    index = Int[]
    for (i, nxyr) in enumerate(native(form, backend, coords))
        nx, ny, nr = nxyr
        inx, iny = round(Int, nx), round(Int, ny)
        1 <= inx <= size(mask,1) && 1 <= iny <= size(mask,2) || continue
        if mask[inx,iny]
            push!(index, i)
        end
    end
    index
end

function find_inmask(backend::Backend, coords, form::Line, mask)
    index = Tuple{Int,Int}[]
    iter = native(form, backend, coords)
    i = start(iter)
    while !done(iter, i)
        nxy, i = next(iter, i)
        nx, ny = nxy
        inx, iny = round(Int, nx), round(Int, ny)
        1 <= inx <= size(mask,1) && 1 <= iny <= size(mask,2) || continue
        if mask[inx,iny]
            push!(index, i)
        end
    end
    index
end

# end # module


text:
### X9.62 spec talks about the way optimal basis is selected as for polynomial and for gaussian normal basis. Currently, selection rule is implemented only for gaussian normal basis

function _isprime(p)

    for i in 2:(p - 1)
        mod(p, i) == 0 && return false
    end

    return true
end

function gn_basis_exist(m, T)

    @assert mod(m, 8) != 0 # May as well return false?

    p = T*m + 1

    !_isprime(p) && return false

    k = compute_integer_order(2, p)

    h = div(T*m, k)

    d = gcd(h, m)

    if d == 1
        return true
    else
        return false
    end

end

function gn_basis_representation_rule(m)

    @assert mod(m, 8) != 0

    if gn_basis_exist(m, 2)
        return 2
    elseif gn_basis_exist(m, 1)
        return 1
    else
        T = 3

        while true
            if gn_basis_exist(m, T)
                return T
            end
            T += 1
        end
    end
end

function specialize(::Type{F2GNB}, N::Int)

    @assert div(N, 8) != 0 "Out of X9.62 spec"
    T = gn_basis_representation_rule(m)

    return F2GNB{N, T}
end

function F2GNB(x::BitVector)

    N = length(x)
    F = specialize(F2GNB, N)

    return F(x)
end

max_stars_repo_name: PeaceFounder/DiffieHellmanGroups

max_issues_repo_name: PeaceFounder/DiffieHellmanGroups

max_forks_repo_name: PeaceFounder/DiffieHellmanGroups


text:
using aoc2019.star_chart: total_orbits, total_transfers
using DelimitedFiles

input = joinpath(@__DIR__, "input")
orbits = readdlm(input, String)

p1 = total_orbits(orbits)
@assert p1 == 301100

p2 = total_transfers(orbits)
@assert p2 == 547

print("-----------------------------------------------------------------------\n")
print("universal orbit map -- part one\n    total orbits    : $p1\n")
print("universal orbit map -- part two\n    total transfers : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
using Test
using DataStructure101
const DS = DataStructure101
using Algorithm101
const ALGO = Algorithm101

@testset "MST" begin
    #= number shows on the edge represents weight:

        v[4]---3---v[1]---2---v[6]
                    | \      / |
                    |  5    7  |
                    |   \  /   |
                    6    v[5]  10
                    |   /  \   |
                    |  1    8  |
                    | /      \ |
        v[3]---4---v[2]---9---v[7]
    =#
    g = DS.Graph(7, :w_list, T=Int64)
    DS.relate!(g, 2, 5, 1)
    DS.relate!(g, 1, 6, 2)
    DS.relate!(g, 1, 4, 3)
    DS.relate!(g, 2, 3, 4)
    DS.relate!(g, 1, 5, 5)
    DS.relate!(g, 1, 2, 6)
    DS.relate!(g, 5, 6, 7)
    DS.relate!(g, 5, 7, 8)
    DS.relate!(g, 2, 7, 9)
    DS.relate!(g, 6, 7, 10)

    @testset "Prim's" begin
        #= number shows on the edge represents order:

            v[4]---2---v[1]---1---v[6]
                          \
                           3
                            \
                             v[5]
                            /  \
                           4    6
                          /      \
            v[3]---5---v[2]      v[7]
        =#
        g_ans = DS.Graph(7, :w_list, T=Int64)
        DS.relate!(g_ans, 1, 6, 2)
        DS.relate!(g_ans, 1, 4, 3)
        DS.relate!(g_ans, 1, 5, 5)
        DS.relate!(g_ans, 2, 5, 1)
        DS.relate!(g_ans, 2, 3, 4)
        DS.relate!(g_ans, 5, 7, 8)

        g_prims = ALGO.mst!(deepcopy(g), algo=:prims)
        @test g_prims.relation == g_ans.relation
        @test g_prims.weight == g_ans.weight
    end

    @testset "Kruskal's" begin
        #= number shows on the edge represents order:

            v[4]---3---v[1]---2---v[6]
                          \
                           5
                            \
                             v[5]
                            /  \
                           1    6
                          /      \
            v[3]---4---v[2]      v[7]
        =#
        g_ans = DS.Graph(7, :w_list, T=Int64)
        DS.relate!(g_ans, 2, 5, 1)
        DS.relate!(g_ans, 1, 6, 2)
        DS.relate!(g_ans, 1, 4, 3)
        DS.relate!(g_ans, 2, 3, 4)
        DS.relate!(g_ans, 1, 5, 5)
        DS.relate!(g_ans, 5, 7, 8)

        g_kruskals = ALGO.mst!(deepcopy(g), algo=:kruskals)
        @test g_kruskals.relation == g_ans.relation
        @test g_kruskals.weight == g_ans.weight
    end
end


text:
"""
This file istarts the ClimateGNSS module by creating the root directory and by
specifying whether the data is to be downloaded or analyzed.  Functionalities
include:
    - Creation of root directory

"""

function gnssroot(dataID::AbstractString)

    path = joinpath("$(homedir())","research","data","GNSS-$(dataID)");
    @warn "$(Dates.now()) - No directory path was given.  Setting to default path: $(path) for ClimateGNSS data downloads."

    if isdir(path)
        @info "$(Dates.now()) - The default path $(path) exists and therefore can be used as a directory for ClimateGNSS data download and analysis from $(dataID)."
    else
        if dataID == "EOS-SG"
            error("$(Dates.now()) - The path $(path) does not exist.  Since you are using EOS-SG GNSS data, please move the data to this directory first before proceeding further.")
        else
            @warn "$(Dates.now()) - The path $(path) does not exist.  A new directory will be created here.  Therefore if you have already uploaded the data, make sure that $(path) is the correct location."
            @info "$(Dates.now()) - Creating path $(path) to hold GNSS data."
            mkpath(path)
        end
    end

    return Dict("raw"=>joinpath(path,"raw"),
                "data"=>joinpath(path,"data"),
                "ana"=>joinpath(path,"ana"));

end

function gnssroot(path::AbstractString,dataID::AbstractString)

    path = joinpath(path,"GNSS-$(dataID)");

    if isdir(path)
        @info "$(Dates.now()) - The directory $(path) exists and therefore can be used for ClimateGNSS data download and analysis from $(dataID)."
    else
        if dataID == "EOS-SG"
            error("$(Dates.now()) - The path $(path) does not exist.  Since you are using EOS-SG GNSS data, please move the data to this directory first before proceeding further.")
        else
            @warn "$(Dates.now()) - The path $(path) does not exist.  A new directory will be created here.  Therefore if you have already uploaded the data, make sure that $(path) is the correct location."
            @info "$(Dates.now()) - Creating path $(path) to hold GNSS data."
            mkpath(path)
        end
    end

    return Dict("raw"=>joinpath(path,"raw"),
                "data"=>joinpath(path,"data"),
                "ana"=>joinpath(path,"ana"));

end


text:
mutable struct AccumulatorLogPsiGPU{N,Ac,Ag,B} <: AbstractMachineValAccumulator
    bnet::N         # A batched version of the cached neural network

    in_buf_cpu::Ac  # the matrix of Nsites x batchsz used as input on cpu
    in_buf_gpu::Ag  # the matrix of Nsites x batchsz used as input on gpu
    out_buf::B      # The row vector of outputs

    buf_n::Int      # Counter for elements in buffer
    batch_sz::Int   # batch size
end

function AccumulatorLogPsi(net::NeuralNetwork, v::Union{gpuAStateBatch, gpuADoubleStateBatch},
                           batch_sz)
    #bnet    = cached(net, batch_sz)
    bnet    = net

    in_buf_cpu  = state_collect(v)
    in_buf_gpu  = deepcopy(v)
    out_buf     = out_similar(net)

    return AccumulatorLogPsiGPU(
        bnet,
        in_buf_cpu, in_buf_gpu, out_buf,
        0, batch_sz)
end

Base.length(a::AccumulatorLogPsiGPU) = a.batch_sz
Base.count(a::AccumulatorLogPsiGPU)  = a.buf_n
isfull(a::AccumulatorLogPsiGPU) = count(a) == length(a)

"""
    init!(c::ScalarBatchAccumulator)

Resets the internal counter of the accumulator, deleting
all previously accumulated (but not computed) values.
"""
init!(c::AccumulatorLogPsiGPU) = c.buf_n = 0
function init!(c::AccumulatorLogPsiGPU, σ)
    init!(c)

    # Reset the state so that later is faster to apply the changes
    state_copy!(c.in_buf_cpu, σ)
end

function (c::AccumulatorLogPsiGPU)(v::Union{AState, ADoubleState})
    @assert !isfull(c) "Pushed data to a full accumulator."

    # Increase the step in our internal buffer
    # this should be guaranteed to always be < max_capacity
    c.buf_n = c.buf_n + 1

    if v isa AState
        buf_i = unsafe_get_batch(c.in_buf_cpu, c.buf_n)
        buf_i .= v
    elseif v isa ADoubleState
        buf_r_i = unsafe_get_batch(row(c.in_buf_cpu), c.buf_n)
        buf_c_i = unsafe_get_batch(col(c.in_buf_cpu), c.buf_n)
        buf_r_i .= row(v)
        buf_c_i .= col(v)
    else
        throw("Not implemented")
    end

    return nothing
end

function (c::AccumulatorLogPsiGPU)(cngs)
    @assert !isfull(c) "Pushed data to a full accumulator."
    c.buf_n = c.buf_n + 1

    σp = unsafe_get_batch(c.in_buf_cpu, c.buf_n)
    apply!(σp, cngs)
end

function (c::AccumulatorLogPsiGPU)(cngs_l, cngs_r)
    @assert !isfull(c) "Pushed data to a full accumulator."
    c.buf_n = c.buf_n + 1

    σp = unsafe_get_batch(c.in_buf_cpu, c.buf_n)
    apply!(σp, cngs_l, cngs_r)
end

function process_accumulator!(c::AccumulatorLogPsiGPU)
    out = c.out_buf

    # copuy to gpu
    copy!(c.in_buf_gpu, c.in_buf_cpu)
    # Compute the batch of logψ neural networks
    logψ!(out, c.bnet, c.in_buf_gpu)

    # Reset the counter of the batch accumulator
    init!(c)

    return out
end


text:
# create pde structure and mesh structure
pde[2], mesh[2] = Preprocessing.initializeexasim(version);

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde[2].model = "ModelD";            # ModelC, ModelD, ModelW
pde[2].modelfile = "pdemodel2";     # name of a file defining the PDE model
include(pde[2].modelfile * ".jl");  # include the PDE model file
pde[2].modelnumber = 2;

# Set discretization parameters, physical parameters, and solver parameters
pde[2].porder = 3;                  # polynomial degree
pde[2].physicsparam = [1.0 0.0];    # thermal conductivity and boundary value
pde[2].tau = [1.0];                 # DG stabilization parameter

# Choose computing platform and set number of processors
#pde[2].platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde[2].mpiprocs = 1;                # number of MPI processors

# set indices to obtain v from the solutions of the other PDE models
# first column : model index
# second column: solution index
pde[2].vindx = [1 1];

# create a linear mesh for a square domain
mesh[2].p, mesh[2].t = Mesh.SquareMesh(8,8,1); # a mesh of 8 by 8 quadrilaterals
# expressions for disjoint boundaries
mesh[2].boundaryexpr = [p -> (p[2,:] .< 1e-3), p -> (p[1,:] .> 1-1e-3), p -> (p[2,:] .> 1-1e-3), p -> (p[1,:] .< 1e-3)];
mesh[2].boundarycondition = [1 1 1 1]; # Set boundary condition for each disjoint boundary

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim