text:
module OpenVR
  using Libdl
  using Cxx
  import Cxx.CxxCore.CppEnum
  import Cxx.CxxCore.CppMFptr
  import Cxx.CxxCore.CppValue
  import Cxx.CxxCore.CppPtr
  import Cxx.CxxCore.CxxQualType
  import Cxx.CxxCore.CppBaseType
  import Cxx.CxxCore.cppcall_member
  import Cxx.CxxCore.cppcall
  import Cxx.CxxCore.CppNNS
  import Cxx.CxxCore.CppTemplate
  import Cxx.CxxCore.CppRef

  reinterpret4cxx(::Type{T}     , x::T     ) where {T}                           = x
  reinterpret4cxx(::Type{Ptr{T}}, x::Ptr{U}) where {T,U}                         = convert(Ptr{T},x)
  reinterpret4cxx(::Type{Ptr{T}}, x::Ref{U}) where {T,U}                         = convert(Ptr{T},pointer_from_objref(x))
  reinterpret4cxx(::Type{T}     , x::U     ) where {T <: Cxx.CxxCore.CppEnum, U} = T(x)
  reinterpret4cxx(::Type{T}     , x::T     ) where {T <: Cxx.CxxCore.CppEnum}    = x

  const openvr_version = "1.9.16"
  const openvr_dirname = "openvr-$openvr_version"
  const module_dir = dirname(dirname(@__FILE__))
  const src_dir = joinpath(module_dir,"src")
  const openvr_dir = joinpath(module_dir,"deps",openvr_dirname)
  const openvr_bin_relpath =
    if     Sys.islinux()   && Sys.ARCH == :x86_32; "linux32/libopenvr_api.so"
    elseif Sys.islinux()   && Sys.ARCH == :x86_64; "linux64/libopenvr_api.so"
    elseif Sys.iswindows() && Sys.ARCH == :x86_32; "win32/openvr_api.dll"
    elseif Sys.iswindows() && Sys.ARCH == :x86_64; "win64/openvr_api.dll"
    elseif Sys.isapple()   && Sys.ARCH == :x86_32; "osx32/libopenvr_api.dylib"
    elseif Sys.isapple()   && Sys.ARCH == :x86_64; "osx64/OpenVR.framework/OpenVR"
    else error("system not detected or unsupported")
    end

  const openvr_binpath = joinpath(openvr_dir,"bin",openvr_bin_relpath)
  const openvr_incpath = joinpath(openvr_dir,"headers")

  function __init__()
    # /Cxx.cpp:1:1: note: '/usr/lib64/gcc/x86_64-pc-linux-gnu/9.2.0/../../../../include/openvr.h' included multiple times, additional include site here
    # __current_compiler__ = Cxx.new_clang_instance()
    addHeaderDir(openvr_incpath, kind=C_System)
    Libdl.dlopen(openvr_binpath, Libdl.RTLD_GLOBAL)
    try
      cxxinclude("openvr.h")
    catch e
      # TODO: this throws when openvr is included multiple times, e.g. first during build and second at `using OpenVR`
    end
  end

  const include_openvr_enums               = joinpath(@__DIR__,"openvr_enums.jl")
  const include_openvr_consts              = joinpath(@__DIR__,"openvr_consts.jl")
  const include_openvr_functions           = joinpath(@__DIR__,"openvr_functions.jl")
  const include_openvr_methods             = joinpath(@__DIR__,"openvr_methods.jl")
  const include_openvr_enums_generated     = joinpath(@__DIR__,"openvr_enums_generated.jl")
  const include_openvr_consts_generated    = joinpath(@__DIR__,"openvr_consts_generated.jl")
  const include_openvr_functions_generated = joinpath(@__DIR__,"openvr_functions_generated.jl")
  const include_openvr_methods_generated   = joinpath(@__DIR__,"openvr_methods_generated.jl")

  cur_output_file = stdout

  # const GENERATE_BINDINGS          = true
  const GENERATE_BINDINGS          = !isfile(include_openvr_methods_generated)
  const INCLUDE_GENERATED_BINDINGS = true
  const DIRECT_DEFINITIONS         = false

  if GENERATE_BINDINGS || DIRECT_DEFINITIONS
    println("creating OpenVR bindings...")
    addHeaderDir(openvr_incpath, kind=C_System)
    Libdl.dlopen(openvr_binpath, Libdl.RTLD_GLOBAL)
    try
      cxxinclude("openvr.h")
    catch e
      # TODO: this throws when openvr is included multiple times, e.g. first during build and second at `using OpenVR`
    end

    macro mkCppConst(cname)
      cxx_quot = :( vr::$cname )
      println(cxx_quot)
      cppconst = __module__.eval(:( @cxx $cxx_quot ))
      if cppconst isa Cxx.CxxCore.CppPtr{UInt8,(true, false, false)}
        printval = "\"" * unsafe_string(reinterpret(Ptr{Int8},cppconst)) * "\""
      elseif cppconst isa Ptr{UInt8}
        printval = "\"" * unsafe_string(convert(Ptr{Int8},cppconst)) * "\""
      elseif cppconst isa CppEnum
        constructor = string(typeof(cppconst).parameters[1])[length("vr::")+1:end]
        printval = "$constructor($(cppconst.val))"
      elseif cppconst isa DataType
        printval = replace(string(cppconst), "Cxx.CxxCore." => "")
      else
        constructor = typeof(cppconst)
        printval = "$constructor($cppconst)"
      end
      if cppconst isa DataType
        GENERATE_BINDINGS && println(cur_output_file,"export $(string(cname))")
        GENERATE_BINDINGS && println(cur_output_file,"const $(string(cname)) = $printval")
        DIRECT_DEFINITIONS || return :()
        return quote
          export $(esc(cname))
          const $(esc(cname)) = $cppconst
        end
      else
        GENERATE_BINDINGS && println(cur_output_file,"const $(string(cname)) = $printval")
        DIRECT_DEFINITIONS || return :()
        return quote
          const $(esc(cname)) = $cppconst
        end
      end
    end
    quote2string(q) = strip(replace(replace(repr(q)[length("quote\n")+1:end-length("\nend")], r" *#=.*=# *\n" => s""), r"\n *" => s"; "))
    joinargstrindent(arr::Vector{String},del,i;forcenl = false) =
      ( length(arr) == 0
      ? ""
      : length(arr) == 1 && !forcenl
      ? arr[1]
      : let istr = repeat(" ",i)
            "\n" * istr * repeat(" ",length(del)) * join(arr,"\n" * istr * del) * "\n" * istr
        end
      )

    macro cppFun(fname)
      cppfun = __module__.eval(quote
        typeof(@cxx vr::&$fname).parameters[1]
      end)
    end
    # reinterpret4cxx(::Type{T},x::Ref{U}) where {T,U} = reinterpret(T,convert(Ptr{U},pointer_from_objref(x)))
    macro mkCppFun(fname)
      cppfun      = __module__.eval(:( typeof($(__module__).@cxx vr::&$fname).parameters[1] ))
      rtype       = cppfun.parameters[1]
      args        = cppfun.parameters[2].parameters
      argstyp     = [:(                      $(esc(Symbol("arg$n")))::$T  ) for (n,T) in enumerate(args)]
      arglist     = [:(                      $(esc(Symbol("arg$n")))      ) for (n,T) in enumerate(args)]
      argsipt     = [:( $reinterpret4cxx($T, $(esc(Symbol("arg$n")))    ) ) for (n,T) in enumerate(args)]
      argstyp_str = [                       "arg$n :: $T" for (n,T) in enumerate(args)]
      arglist_str = [                       "arg$n"       for (n,T) in enumerate(args)]
      argsipt_str = [   "reinterpret4cxx($T, arg$n)"      for (n,T) in enumerate(args)]
      q1 = quote
        function $(esc(fname))( $(argstyp...) )::$rtype
          # $Cxx.@cxx vr::$fname( $(arglist...) )
          # ($(__module__.Cxx.CxxCore.cppcall))(
          #      $(__module__.__current_compiler__)
          #   , ($(__module__.Cxx.CxxCore.CppNNS)){(Tuple){:vr, $(QuoteNode(Symbol(fname)))}}()
          #     , $(arglist...) )
          $cppcall( $__current_compiler__
                  , $CppNNS{(Tuple){:vr, $(QuoteNode(Symbol(fname)))}}()
                  , $(arglist...) )
        end
      end
      s1 = """
      $(string(fname))($(joinargstrindent(argstyp_str,", ",2))) $(length(arglist) < 2 ? "" : "     "):: $rtype =
        cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :$(string(fname))}}()$(length(arglist_str) > 0 ? ", " : "")$(join(arglist_str,", ")))"""
      s1 = replace(s1, "Cxx.CxxCore." => "")
      GENERATE_BINDINGS && println(cur_output_file,"export $fname")
      GENERATE_BINDINGS && println(cur_output_file,s1)
      if isempty(arglist)
        GENERATE_BINDINGS && println(cur_output_file)
        DIRECT_DEFINITIONS || return :()
        return quote
            export $(esc(fname))
            $q1
          end
      else
        q2 = quote
               function $(esc(fname))( $(arglist...) )::$rtype
                 $(esc(fname))( $(argsipt...) )
               end
             end
        s2 = """
        $(string(fname))($(join(arglist_str,", "))) =
          $(string(fname))($(joinargstrindent(argsipt_str,", ",4)))"""
        GENERATE_BINDINGS && println(cur_output_file,s2)
        GENERATE_BINDINGS && println(cur_output_file)
        DIRECT_DEFINITIONS || return :()
        return quote
            export $(esc(fname))
            $q1
            $q2
          end
      end
    end

    macro mkCppMethodCast(cname,fname,rtype,argtypes)
      argtypes_str =
        if argtypes isa Expr && argtypes.head == :tuple
          argtypes_str = join(string.(argtypes.args), ",")
        elseif argtypes isa String
          if startswith(argtypes,"(") && endswith(argtypes,")")
            argtypes[2:end-1]
          else
            argtypes
          end
        else
          string(argtypes)
        end
      # .method("SetString", static_cast<void(vr::CVRSettingHelper::*)(const char *, const char *, const char *, vr::EVRSettingsError *)>(&vr::CVRSettingHelper::SetString))
      cname_str = string(cname)
      fname_str = string(fname)
      rtype_str = string(rtype)
      icxx_str = "static_cast<$rtype_str(vr::$cname_str::*)($argtypes_str)>(&vr::$cname_str::$fname_str);"
      println(icxx_str)
      cppmet = __module__.eval(:( typeof(@icxx_str $icxx_str)) )
      mkCppMethod_expr(cppmet,cname,fname; createConvenienceWrapper=false)
    end

    macro mkCppMethod(cname,fname)
      icxx_str = "&vr::$cname::$fname;"
      println(icxx_str)
      cppmet  = __module__.eval(:( typeof(@icxx_str $icxx_str) ))
      mkCppMethod_expr(cppmet,cname,fname)
    end

    macro mkCppMethodInline(cname,fname)
      cxx_quot = :( vr::$cname::&$fname )
      println(cxx_quot)
      cppmet  = __module__.eval(:( typeof(@cxx $cxx_quot) ))
      mkCppMethod_expr(cppmet,cname,fname)
    end

    function mkCppMethod_expr(cppmet,cname,fname; createConvenienceWrapper = true)
      cppmet <: Cxx.CxxCore.CppMFptr || return :( error("vr::$($cname)::$($fname) is not a method reference") )
      this = ( cppmet.parameters[1] <: (Cxx.CxxCore.CppValue{U,N} where {U,N})
             ? Cxx.CxxCore.CppPtr{cppmet.parameters[1].body.parameters[1], N} where N
             : return :( error("method call on *this :: $($(cppmet.parameters[1]))") )
             )
      cppfun  = cppmet.parameters[2]
      rtype   = cppfun.parameters[1]
      args    = cppfun.parameters[2].parameters
      argstyp = [:(                     $(esc(Symbol("arg$n")))::$T  ) for (n,T) in enumerate(args)]
      arglist = [:(                     $(esc(Symbol("arg$n")))      ) for (n,T) in enumerate(args)]
      argsipt = [:( $reinterpret4cxx($T, $(esc(Symbol("arg$n")))    ) ) for (n,T) in enumerate(args)]
      argstyp_str = [                       "arg$n :: $T" for (n,T) in enumerate(args)]
      arglist_str = [                       "arg$n"       for (n,T) in enumerate(args)]
      argsipt_str = [   "reinterpret4cxx($T, arg$n)"      for (n,T) in enumerate(args)]
      pushfirst!(argstyp, :( $(esc(Symbol("this")))::$this ))
      pushfirst!(arglist, :( $(esc(Symbol("this"))) ))
      pushfirst!(argsipt, :( $(esc(Symbol("this"))) ))
      pushfirst!(argstyp_str, "this :: $this" )
      pushfirst!(arglist_str, "this" )
      pushfirst!(argsipt_str, "this" )
      q1 = quote
        function $(esc(fname))( $(argstyp...) )::$rtype
          # $Cxx.@cxx vr::$fname( $(arglist...) )
          ($(Cxx.CxxCore.cppcall_member))(
               $(__current_compiler__)
            , ($(Cxx.CxxCore.CppNNS)){(Tuple){$fname}}()
            , $(arglist...)
            )
        end
      end
      s1 = """
      $(string(fname))($(joinargstrindent(argstyp_str,", ",2; forcenl=true)))      :: $rtype =
        cppcall_member(__current_compiler__, CppNNS{(Tuple){:$(string(fname))}}()$(length(arglist_str) > 0 ? ", " : "")$(join(arglist_str,", ")))"""
      s1 = replace(s1, "Cxx.CxxCore." => "")
      GENERATE_BINDINGS && println(cur_output_file,"export $fname")
      GENERATE_BINDINGS && println(cur_output_file,s1)
      if length(arglist) == 1 || !createConvenienceWrapper
        GENERATE_BINDINGS && println(cur_output_file)
        DIRECT_DEFINITIONS || return :()
        return quote
            export $(esc(fname))
            $q1
          end
      else
        q2 = quote
               function $(esc(fname))( $(vcat(argstyp[1:1],arglist[2:end])...) )::$rtype
                 $(esc(fname))( $(argsipt...) )
               end
             end
        s2 = """
        $(string(fname))($(join(vcat(argstyp_str[1:1],arglist_str[2:end]),", "))) =
          $(string(fname))($(joinargstrindent(argsipt_str,", ",4)))"""
        s2 = replace(s2, "Cxx.CxxCore." => "")
        GENERATE_BINDINGS && println(cur_output_file,s2)
        GENERATE_BINDINGS && println(cur_output_file)
        DIRECT_DEFINITIONS || return :()
        return quote
            export $(esc(fname))
            $q1
            $q2
          end
      end
    end

    println("creating $include_openvr_enums_generated...")
    cur_output_file = open(include_openvr_enums_generated, "w")
      include(include_openvr_enums)
    close(cur_output_file)
    println("creating $include_openvr_consts_generated...")
    cur_output_file = open(include_openvr_consts_generated, "w")
      include(include_openvr_consts)
    close(cur_output_file)
    println("creating $include_openvr_functions_generated...")
    cur_output_file = open(include_openvr_functions_generated, "w")
      include(include_openvr_functions)
    close(cur_output_file)
    println("creating $include_openvr_methods_generated...")
    cur_output_file = open(include_openvr_methods_generated, "w")
      include(include_openvr_methods)
    close(cur_output_file)
    cur_output_file = stdout
  end

  if INCLUDE_GENERATED_BINDINGS
    include(include_openvr_enums_generated)
    include(include_openvr_consts_generated)
    include(include_openvr_functions_generated)
    include(include_openvr_methods_generated)
  end

  # julia> cxx"struct A{ int x; };"
  # true
  # julia> @cxxnew A()
  # (struct A *) @0x000055f2048a1eb0
  # julia> typeof(@cxxnew A())
  # Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{:A},(false, false, false)},N} where N,(false, false, false)}

  # ivrsystem
  # (class vr::IVRSystem *) @0x00007f543a2f00e0
  # julia> typeof(ivrsystem)
  # Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("vr::IVRSystem")},(false, false, false)},(false, false, false)}

  # @cxx vr::IVRSystem::&GetEyeToHeadTransform
  # Cxx.CxxCore.CppMFptr{
  #   Cxx.CxxCore.CppValue{
  #     Cxx.CxxCore.CxxQualType{
  #       Cxx.CxxCore.CppBaseType{Symbol("vr::IVRSystem")},(false, false, false)}
  #     , N
  #     } where N
  #   , Cxx.CxxCore.CppFunc{
  #       Cxx.CxxCore.CppValue{
  #         Cxx.CxxCore.CxxQualType{
  #           Cxx.CxxCore.CppBaseType{Symbol("vr::HmdMatrix34_t")},(false, false, false)},N} where N,Tuple{Cxx.CxxCore.CppEnum{Symbol("vr::EVREye"),UInt32}}}}(0x0000000000000021, 0x0000000000000000)

  # @mkCppConst IVRSystem
  # @mkCppConst IVRApplications
  # @mkCppConst IVRSettings
  # @mkCppConst IVRChaperone
  # @mkCppConst IVRChaperoneSetup
  # @mkCppConst IVRCompositor
  # @mkCppConst IVRNotifications
  # @mkCppConst IVROverlay
  # @mkCppConst IVRRenderModels
  # @mkCppConst IVRExtendedDisplay
  # @mkCppConst IVRTrackedCamera
  # @mkCppConst IVRScreenshots
  # @mkCppConst IVRResources
  # @mkCppConst IVRDriverManager
  # @mkCppConst IVRInput
  # @mkCppConst IVRIOBuffer
  # @mkCppConst IVRSpatialAnchors
  # @mkCppConst COpenVRContext

  const IVRSystem_ptr          = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSystem")},(false, false, false)},N} where N
  const IVRApplications_ptr    = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRApplications")},(false, false, false)},N} where N
  const IVRSettings_ptr        = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSettings")},(false, false, false)},N} where N
  const IVRChaperone_ptr       = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRChaperone")},(false, false, false)},N} where N
  const IVRChaperoneSetup_ptr  = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRChaperoneSetup")},(false, false, false)},N} where N
  const IVRCompositor_ptr      = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRCompositor")},(false, false, false)},N} where N
  const IVRNotifications_ptr   = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRNotifications")},(false, false, false)},N} where N
  const IVROverlay_ptr         = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVROverlay")},(false, false, false)},N} where N
  const IVRRenderModels_ptr    = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRRenderModels")},(false, false, false)},N} where N
  const IVRExtendedDisplay_ptr = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRExtendedDisplay")},(false, false, false)},N} where N
  const IVRTrackedCamera_ptr   = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRTrackedCamera")},(false, false, false)},N} where N
  const IVRScreenshots_ptr     = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRScreenshots")},(false, false, false)},N} where N
  const IVRResources_ptr       = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRResources")},(false, false, false)},N} where N
  const IVRDriverManager_ptr   = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRDriverManager")},(false, false, false)},N} where N
  const IVRInput_ptr           = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRInput")},(false, false, false)},N} where N
  const IVRIOBuffer_ptr        = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRIOBuffer")},(false, false, false)},N} where N
  const IVRSpatialAnchors_ptr  = CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSpatialAnchors")},(false, false, false)},N} where N
  const COpenVRContext_ptr     = CppPtr{CxxQualType{CppBaseType{Symbol("vr::COpenVRContext")},(false, false, false)},N} where N

  #  NOTE: it also might help to set the vr*-process' niceness to 19 and move the vrcompmositor windows "out of the visual viewport"
  function killall()
    println("VR_Shutdown...")
    VR_Shutdown()
    println("waiting for VR to shut down")
    sleep(2.0)
    println("sending term signal...")
    try; run(`pkill --signal TERM -f vrserver`);     catch e; end
    try; run(`pkill --signal TERM -f vrdashboard`);  catch e; end
    try; run(`pkill --signal TERM -f vrwebhelper`);  catch e; end
    try; run(`pkill --signal TERM -f vrcompositor`); catch e; end
    sleep(2.0)
    println("sending kill signal...")
    try; run(`pkill --signal KILL -f vrserver`);     catch e; end
    try; run(`pkill --signal KILL -f vrdashboard`);  catch e; end
    try; run(`pkill --signal KILL -f vrwebhelper`);  catch e; end
    try; run(`pkill --signal KILL -f vrcompositor`); catch e; end
    nothing
  end
end


text:
include("types.jl")
include("shared.jl")
include("util.jl")
include("kernel.jl")
include("adv_kernel_cg.jl")

## set
solver = :mosek
# solver = :gurobi

log = 0
psdtol = 1e-6
perturb = 1e-12
obj_reltol_cv = 0.0
obj_reltol_test = 0.0
verbose = false

smaller_cv = false   # 3 instead of 5

### prepare data

dname = "glass"
D_all = readcsv("data-example/" * dname * ".csv")
id_train = readcsv("data-example/" * dname * ".train")
id_test = readcsv("data-example/" * dname * ".test")

id_train = round(Int64, id_train)
id_test = round(Int64, id_test)

println(dname)

### Cross Validation, using first split
## First stage

id_tr = vec(id_train[1,:])
id_ts = vec(id_test[1,:])
X_train = D_all[id_tr,1:end-1]
y_train = round(Int, D_all[id_tr, end])

X_test = D_all[id_ts,1:end-1]
y_test = round(Int, D_all[id_ts, end])

X_train, mean_vector, std_vector = standardize(X_train)
X_test = standardize(X_test, mean_vector, std_vector)

if smaller_cv
  Cs =  [2.0^i for i=0:4:8]
  Gs =  [2.0^(-8+i) for i=0:4:8]
else
  Cs =  [2.0^i for i=0:3:12]
  Gs =  [2.0^(-12+i) for i=0:3:12]
end

ncs = length(Cs)

Pars = [ Tuple{Float64,Float64}((Cs[i], Gs[j])) for i=1:ncs, j=1:ncs ]
Pars = vec(Pars)
npar = length(Pars)

# fold
n_train = size(X_train, 1)
n_test = size(X_test, 1)
kf = 5

# k folds
folds = k_fold(n_train, kf)

loss_list = zeros(npar)
loss01_list = zeros(npar)

# The first stage of CV
idx = randperm(n_train)
X_train = X_train[idx,:]
y_train = y_train[idx]

println("First CV")

for i = 1:npar

  println(i, " | Adversarial | C = ", Pars[i][1], ", Gamma = ", Pars[i][2])

  losses = zeros(n_train)
  losses1 = zeros(n_train)
  # k fold
  for j = 1:kf
    # prepare training and validation
    id_tr = vcat(folds[[1:j-1; j+1:end]]...)
    id_val = folds[j]

    X_tr = X_train[id_tr, :];  y_tr = y_train[id_tr]
    X_val = X_train[id_val, :];  y_val = y_train[id_val]

    C = Pars[i][1]
    gamma = Pars[i][2]

    print("    ",j, "-th fold : ")
    @time model = train_adv_kernel_cg(X_tr, y_tr, C, :gaussian, [gamma], perturb=perturb, obj_reltol=obj_reltol_cv, solver=solver, log=log, psdtol=psdtol, verbose=verbose)

    _, ls, _, ls1, _, _ = test_adv_kernel(model, X_val, y_val, X_tr, y_tr)

    losses[id_val] = ls
    losses1[id_val] = ls1

  end

  loss_list[i] = mean(losses)
  loss01_list[i] = mean(losses1)
  # println("loss : ", string(mean(losses)))
  println("  => loss01 : ", string(mean(losses1)))

end

ind_max= indmin(loss01_list)
C0 = Pars[ind_max][1]
G0 = Pars[ind_max][2]
if smaller_cv
  Cs =  [C0*2.0^(i-2) for i=1:3]
  Gs =  [G0*2.0^(i-2) for i=1:3]
else
  Cs =  [C0*2.0^(i-3) for i=1:5]
  Gs =  [G0*2.0^(i-3) for i=1:5]
end
ncs = length(Cs)

Pars = [ Tuple{Float64,Float64}((Cs[i], Gs[j])) for i=1:ncs, j=1:ncs ]
Pars = vec(Pars)
npar = length(Pars)

## Second stage
idx = randperm(n_train)
X_train = X_train[idx,:]
y_train = y_train[idx]

println("Second CV")
for i = 1:npar

  println(i, " | Adversarial | C = ", Pars[i][1], ", Gamma = ", Pars[i][2])

  losses = zeros(n_train)
  losses1 = zeros(n_train)
  # k fold
  for j = 1:kf
    # prepare training and validation
    id_tr = vcat(folds[[1:j-1; j+1:end]]...)
    id_val = folds[j]

    X_tr = X_train[id_tr, :];  y_tr = y_train[id_tr]
    X_val = X_train[id_val, :];  y_val = y_train[id_val]

    C = Pars[i][1]
    gamma = Pars[i][2]

    print("    ",j, "-th fold : ")
    @time model = train_adv_kernel_cg(X_tr, y_tr, C, :gaussian, [gamma], perturb=perturb, obj_reltol=obj_reltol_cv, solver=solver, log=log, psdtol=psdtol, verbose=verbose)

    _, ls, _, ls1, _, _ = test_adv_kernel(model, X_val, y_val, X_tr, y_tr)

    losses[id_val] = ls
    losses1[id_val] = ls1
    #println("loss : ", string(ls))
    #println("loss01 : ", string(ls))

  end

  loss_list[i] = mean(losses)
  loss01_list[i] = mean(losses1)
  #println("loss : ", string(mean(losses)))
  println("  => loss01 : ", string(mean(losses1)))

end

ind_max= indmin(loss01_list)
C_best = Pars[ind_max][1]
G_best = Pars[ind_max][2]

### Evaluation

n_split = size(id_train, 1)

v_model = Vector{ClassificationModel}()
v_result = Vector{Tuple}()
v_acc = zeros(n_split)
v_cs_result = zeros(n_split, 5)

println("Evaluation")
for i = 1:n_split
  # standardize
  id_tr = vec(id_train[i,:])
  id_ts = vec(id_test[i,:])
  X_train = D_all[id_tr,1:end-1]
  y_train = round(Int, D_all[id_tr, end])

  X_test = D_all[id_ts,1:end-1]
  y_test = round(Int, D_all[id_ts, end])

  X_train, mean_vector, std_vector = standardize(X_train)
  X_test = standardize(X_test, mean_vector, std_vector)

  #train and test
  @time model = train_adv_kernel_cg(X_train, y_train, C_best, :gaussian, [G_best], perturb=perturb, obj_reltol=obj_reltol_test, solver=solver, log=log, psdtol=psdtol, verbose=verbose)

  result = test_adv_kernel(model, X_test, y_test, X_train, y_train)
  loss01 = result[3]
  acc = 1.0 - loss01
  cs_result = count_constraints(model, y_train)

  println("accuracy : ", acc)

  push!(v_model, model)
  push!(v_result, result)
  v_acc[i] = acc
  v_cs_result[i, :] = collect(cs_result)
end

println(dname)
println("mean accuracy : ", mean(v_acc))
println("std accuracy : ", std(v_acc))

ext: jl

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
using CompactBasisFunctions
using Test

include("vandermonde_tests.jl")
include("basis_tests.jl")
include("bernstein_tests.jl")
include("chebyshev_tests.jl")
include("lagrange_tests.jl")
include("legendre_tests.jl")


text:
module TestIntroduction
include("../../../examples/introduction.jl")
end  # module


text:
module StipplePlotlyExport

import Genie, Stipple, StipplePlotly, PlotlyBase, Base64

function save(data::Union{StipplePlotly.Charts.PlotData,Vector{StipplePlotly.Charts.PlotData}},
              layout::StipplePlotly.Charts.PlotLayout = StipplePlotly.Charts.PlotLayout();
              width::Int = 700,
              height::Int = 500,
              scale::Real = 1,
              format::String = "png")

  format in PlotlyBase.ALL_FORMATS || error("Unknown format $format. Expected one of $(PlotlyBase.ALL_FORMATS)")

  payload = Dict(
    :data => Dict(
      :data => Stipple.render(data),
      :layout => Stipple.render(layout)
    ),
    :width => width,
    :height => height,
    :scale => scale,
    :format => format
  )

  PlotlyBase._start_kaleido_process()
  PlotlyBase._ensure_kaleido_running() || PlotlyBase._restart_kaleido_process()

  # Integrally copied from PlotlyBase.jl -- big thank you to the contributors!

  bytes = transcode(UInt8, Stipple.JSONParser.json(payload))
  write(PlotlyBase.P.stdin, bytes)
  write(PlotlyBase.P.stdin, transcode(UInt8, "\n"))
  flush(PlotlyBase.P.stdin)

  # read stdout and parse to json
  res = readline(PlotlyBase.P.stdout)
  js = Stipple.JSONParser.parse(res)

  # check error code
  code = get(js, "code", 0)
  if code != 0
      msg = get(js, "message", nothing)
      error("Transform failed with error code $code: $msg")
  end

  # get raw image
  img = String(js["result"])

  # base64 decode if needed, otherwise transcode to vector of byte
  if format in PlotlyBase.TEXT_FORMATS
      return transcode(UInt8, img)
  else
      return Base64.base64decode(img)
  end
end

function save(io::Union{String,IO},
              data::Union{StipplePlotly.Charts.PlotData,Vector{StipplePlotly.Charts.PlotData}},
              layout::StipplePlotly.Charts.PlotLayout = StipplePlotly.Charts.PlotLayout();
              width::Int = 700,
              height::Int = 500,
              scale::Real = 1,
              format::String = "png")
  write(io, save(data, layout; width=width, height=height, scale=scale, format=format))
end

end


text:
module Keys

import Base: getindex, haskey, merge, convert, &, |, @pure, Bool, tail, Generator
import Base.Meta: quot
import MacroTools: @capture
import Base.Iterators: flatten
import Base: Generator

include("typed_bools.jl")
include("recur_unroll.jl")

export Key
"""
    struct Key{K}

A typed key. See [`@__str`](@ref) for an easy way to create keys. Use to create
[`Keyed`](@ref) values.
"""
struct Key{K} end

inner_value(k::Key{K}) where K = K

export @__str
"""
    @__str

make a key

```jldoctest
julia> using Keys

julia> _"a"
.a
```
"""
macro __str(s::String)
    esc(:($Key{$(quot(Symbol(s)))}()))
end

const SomeKeys = NTuple{N, Key} where N
const PairOfKeys = Pair{T1, T2} where {T1 <: Key, T2 <: Key}

function Base.show(io::IO, key::Key{K}) where K
    print(io, :.)
    print(io, K)
end

export Keyed
"""
    struct Keyed{K, V}

an alias for a [`Key`](@ref)-value pair. a tuple of `Keyed` values is aliased
as a [`KeyedTuple`](@ref).
"""
const Keyed{K, V} = Pair{Key{K}, V} where {K, V}

export key
"""
    key(keyed::Keyed)

get the key of a [`Keyed`](@ref) value.

```jldoctest
julia> using Keys

julia> key.((_"a" => 1, _"b" => 2))
(.a, .b)
```
"""
key(::Keyed{K}) where K = Key{K}()

export value
"""
    value(key::Keyed)

get the value of a [`Keyed`](@ref) value.

```jldoctest
julia> using Keys

julia> value.((_"a" => 1, _"b" => 2))
(1, 2)
```
"""
value(keyed_tuple::Keyed) = keyed_tuple.second

export KeyedTuple
"""
    const KeyedTuple

A tuple with only [`Keyed`](@ref) values. You can index them with [`Key`](@ref)s or
access them with dots. Duplicated keys are allowed; will return the
first match.

```jldoctest
julia> using Keys

julia> keyed_tuple = (_"a" => 1, _"b" => 2)
(.a => 1, .b => 2)

julia> keyed_tuple.b
2

julia> keyed_tuple[(_"a", _"b")]
(.a => 1, .b => 2)

julia> keyed_tuple.c
ERROR: Key .c not found
[...]

julia> haskey(keyed_tuple, _"b")
True()

julia> merge(keyed_tuple, (_"a" => 4, _"c" => 3))
(.b => 2, .a => 4, .c => 3)
```
"""
const KeyedTuple = Tuple{Keyed, Vararg{Keyed}}

match_key(::Keyed{K}, ::Key{K}) where K = True()
match_key(::Keyed, ::Key) = False()

function match_key(keyed::Keyed, keys::SomeKeys)
    reduce_unrolled(|, map(
        let keyed = keyed
            key -> match_key(keyed, key)
        end,
        keys
    ))
end

first_error(::Tuple{}, key::Key) = error("Key $key not found")
first_error(keyed_tuple::KeyedTuple, key::Key) = value(keyed_tuple[1])

which_key(keyed_tuple::KeyedTuple, key::Union{Key, SomeKeys}) = map(
    let key = key
        keyed -> match_key(keyed, key)
    end,
    keyed_tuple
)

_getindex(keyed_tuple, keys) =
    getindex_unrolled(keyed_tuple, which_key(keyed_tuple, keys))

getindex(keyed_tuple::KeyedTuple, key::Key) =
    first_error(_getindex(keyed_tuple, key), key)

getindex(keyed_tuple::KeyedTuple, keys::SomeKeys) =
    _getindex(keyed_tuple, keys)

haskey(keyed_tuple::KeyedTuple, key::Key) =
    reduce_unrolled(|, which_key(keyed_tuple, key))

export delete
"""
    delete(keyed_tuple::KeyedTuple, keys::Key...)

delete all [`Keyed`](@ref) values matching [`Key`](@ref)s in a [`KeyedTuple`](@ref).

```jldoctest
julia> using Keys

julia> delete((_"a" => 1, _"b" => 2), _"a")
(.b => 2,)
```
"""
delete(keyed_tuple::KeyedTuple, keys::Key...) =
    getindex_unrolled(keyed_tuple, map(
        not,
        which_key(keyed_tuple, keys)))

export push
"""
    push(keyed_tuple::KeyedTuple, pairs::Keyed...)

push the [`Keyed`](@ref) values in `pairs` into the
[`KeyedTuple`](@ref), replacing common [`Key`](@ref)s.

```jldoctest
julia> using Keys

julia> push((_"a" => 1, _"b" => 2), _"b" => 4, _"c" => 3)
(.a => 1, .b => 4, .c => 3)
```
"""
push(keyed_tuple::KeyedTuple, pairs::Keyed...) =
    delete(keyed_tuple, key.(pairs)...)..., pairs...

@inline Base.getproperty(key::KeyedTuple, s::Symbol) = getindex(key, Key{s}())

export map_values
"""
    map_values(f, keyed_tuple::KeyedTuple)

map `f` over the values of a [`KeyedTuple`](@ref).

```jldoctest
julia> using Keys

julia> map_values(x -> x + 1, (_"a" => 1, _"b" => 2))
(.a => 2, .b => 3)
```
"""
map_values(f, keyed_tuple::KeyedTuple) = map(
    let f = f
        keyed -> key(keyed) => f(value(keyed))
    end,
    keyed_tuple
)

rename_one(pair_of_keys::Keyed{New, Key{Old}}, old_keyed::Keyed{Old}) where {Old, New} =
    pair_of_keys.first => old_keyed.second
rename_one(pair_of_keys::PairOfKeys, old_keyed::Keyed) = old_keyed

rename_single(pair_of_keys::PairOfKeys, ::Tuple{}) = ()
rename_single(pair_of_keys::PairOfKeys, keyed_tuple::KeyedTuple) =
    rename_one(pair_of_keys, keyed_tuple[1]),
    rename_single(pair_of_keys, tail(keyed_tuple))...

rename(keyed_tuple::KeyedTuple) = keyed_tuple

export rename
"""
    rename(keyed_tuple::KeyedTuple, pairs_of_keys::PairOfKeys...)

for each pair of [`Key`](@ref)s, where the first key matches in
[`KeyedTuple`](@ref), it will be replaced by the second.

```jldoctest
julia> using Keys

julia> rename((_"a" => 1, _"b" => 2), _"c" => _"a")
(.c => 1, .b => 2)
```
"""
rename(keyed_tuple::KeyedTuple, pairs_of_keys::PairOfKeys...) =
    rename(rename_single(pairs_of_keys[1], keyed_tuple), tail(pairs_of_keys)...)

common_keys(x::KeyedTuple, y::KeyedTuple) =
    first.(filter_unrolled(pair -> pair[1] === pair[2], product_unrolled(key.(x), key.(y))))

merge(a::KeyedTuple, b::KeyedTuple) =
    (delete(a, common_keys(a, b)...)..., b...)

#=
@require DataFrames begin

    import DataFrames: DataFrame

    (::Type{KeyedTuple})(d::DataFrame) = (map(
        (name, column) -> Key{name}() => (column),
        names(d),
        d.columns
    )...)

    DataFrame(k::KeyedTuple) =
        DataFrame(map(x -> inner_value(x.first) => x.second, k)...)

end
=#

substitute_underscores!(dictionary, body) = body
substitute_underscores!(dictionary, body::Symbol) =
    if all(isequal('_'), string(body))
        if !haskey(dictionary, body)
            dictionary[body] = gensym("argument")
        end
        dictionary[body]
    else
        body
    end
substitute_underscores!(dictionary, body::Expr) =
    if body.head == :quote
        body
    elseif @capture body @_ args__
        body
    else
        Expr(body.head,
            map(body -> substitute_underscores!(dictionary, body), body.args)
        ...)
    end

string_length(something) = something |> String |> length

function anonymize(body, line, file)
    dictionary = Dict{Symbol, Symbol}()
    new_body = substitute_underscores!(dictionary, body)
    sorted_dictionary = sort(
        lt = (pair1, pair2) ->
            isless(string_length(pair1.first), string_length(pair2.first)),
        collect(dictionary)
    )
    Expr(:->,
        Expr(:tuple, Generator(pair -> pair.second, sorted_dictionary)...),
        Expr(:block, LineNumberNode(line, file), new_body)
    )
end

export @_
"""
    macro _(body::Expr)

Another syntax for anonymous functions. The arguments are inside the body; the
first arguments is `_`, the second argument is `__`, etc.

```jldoctest
julia> using Keys

julia> 1 |> (@_ _ + 1)
2

julia> map((@_ __ - _), (1, 2), (2, 1))
(1, -1)
```
"""
macro _(body::Expr)
    anonymize(body, @__LINE__, @__FILE__) |> esc
end

build_call(afunction, arguments, parity, line, file) =
    if length(arguments) >= parity
        anonymous_arguments = ((anonymize(argument, line, file), quot(argument))
            for argument in arguments[parity+1:end])
        Expr(:call,
            afunction,
            arguments[1:parity]...,
            flatten(anonymous_arguments)...
        )
    else
        error("Expecting at least $parity argument(s)")
    end

anonymize_arguments(atail, line, file) =
    if @capture atail numberedfunction_(arguments__)
        string_function = string(numberedfunction)
        parity = tryparse(Int, string(string_function[end]))
        if parity == nothing
            atail
        else
            build_call(Symbol(chop(string_function)), arguments, parity, line, file)
        end
    else
        atail
    end

query(body, line, file)  =
    if @capture body head_ |> atail_
        Expr(:call, anonymize(anonymize_arguments(atail, line, file), line, file),
            query(head, line, file)
        )
    else
        body
    end

export @query
"""
    macro query(body::Expr)

Query your code. If body is a chain `head_ |> tail_`, recur on
head. If tail is a function call, and the function ends
with a number (the parity), anonymize and quote arguments past that parity.
Either way, anonymize the whole tail, then call it on head.

```jldoctest
julia> using Keys

julia> call(source1, source2, anonymous, quoted) = anonymous(source1, source2);

julia> @query 1 |> (_ - 2) |> abs(_) |> call2(_, 2, _ + __)
3

julia> @query 1 |> call2(_)
ERROR: LoadError: Expecting at least 2 argument(s)
[...]
```
"""
macro query(body)
    line = @__LINE__
    file = @__FILE__
    query(body, line, file) |> esc
end

end


text:
@testset "MRI tests" begin

phasefile = joinpath("data", "small", "Phase.nii")
magfile = joinpath("data", "small", "Mag.nii")
echo = 3
echo2 = 2
phase4D = niread(phasefile).raw
mag4D = niread(magfile).raw
phase = phase4D[:,:,:,echo]
mag = mag4D[:,:,:,echo]
phase2 = niread(phasefile).raw[:,:,:,echo2]
TEs = [echo, echo2]

## accuracy tests
function unwrap_test(wrapped; keyargs...)
    unwrapped = unwrap(wrapped; keyargs...)
    # test that unwrapped is not a copy of phase
    @test unwrapped != wrapped
    # test that all resulting values are only 2π different
    @test all(isapprox.(rem2pi.(unwrapped - wrapped, RoundNearest), 0; atol=1e-5))
    unwrapped
end

t = []
push!(t, unwrap_test(phase))
push!(t, unwrap_test(phase; mag=mag))
#push!(t, unwrap_test(phase4D)) #TODO use test data set with noise to see difference
push!(t, unwrap_test(phase4D; mag=mag4D, TEs=TEs))
push!(t, unwrap_individual(phase4D; mag=mag4D, TEs=TEs))
push!(t, unwrap_test(phase; weights=:bestpath))
#push!(t, unwrap_test(phase; weights=:romeo, mask=robustmask(mag)))
push!(t, unwrap_test(phase; weights=:romeo2, mag=mag, TEs=TEs, phase2=phase2))
push!(t, unwrap_test(phase; weights=:romeo3, mag=mag, TEs=TEs, phase2=phase2))
push!(t, unwrap_test(phase; weights=:romeo4, mag=mag, TEs=TEs, phase2=phase2))
push!(t, unwrap_test(phase; mag=mag, maxseeds=50))
push!(t, unwrap_test(phase4D; mag=mag4D, TEs=TEs, temporal_uncertain_unwrapping=true))

# all results should be different
for i in 1:length(t), j in 1:(i-1)
    @test t[i] != t[j]
end

t = []
push!(t, unwrap_test(phase; mag=mag, maxseeds=50))
push!(t, unwrap_test(phase; mag=mag, maxseeds=50, merge_regions=true, correct_regions=true))
#TODO correct_regions does not affect outcome
#push!(t, unwrap_test(phase; mag=mag, maxseeds=50, merge_regions=true))
for i in 1:length(t), j in 1:(i-1)
    @test t[i] != t[j]
end

@test unwrap_individual(phase4D; mag=mag4D, TEs=TEs) == unwrap(phase4D; mag=mag4D, TEs=TEs, individual=true)

## performance tests (at end to avoid first run overhead)
@test (@timed unwrap(phase))[5].poolalloc < 6e3
@test (@timed unwrap(phase; mag=mag))[5].poolalloc < 7e3
@test (@timed unwrap(phase; weights=:bestpath))[5].poolalloc < 3e4

## NaN tests
nanphase = copy(phase)
nanphase[1,:,:] .= NaN
nan_unwrapped = unwrap_test(phase)
nan_unwrapped[1,:,:] .= NaN
@test nan_test(unwrap(nanphase), nan_unwrapped)

nanmag = copy(mag)
nanmag[1,:,:] .= NaN
@test nan_test(unwrap(phase; mag=nanmag)[2:end,:,:], unwrap_test(phase; mag=mag)[2:end,:,:])

end


text:
"""
Create a tensor product of the supplied arguments.

The function `tensorproduct` applies some simplifications and does not necessarily
return a Product type.

A `tensorproduct(a)` with just a single element returns `a`.

For integer `n`, `tensorproduct(a, n)` becomes `tensorproduct(a, a, ..., a)`.
A type-safe variant is `tensorproduct(a, Val{N})`.
"""
tensorproduct() = nothing

# Don't create a tensor product of just one element
tensorproduct(a::Tuple) = tensorproduct(a...)
tensorproduct(a::AbstractArray) = tensorproduct(a...)
tensorproduct(a) = a

# Create a tensor product with n times the same element
tensorproduct(a, n::Int) = tensorproduct(ntuple(t->a, n)...)

tensorproduct(a, ::Type{Val{N}}) where {N} = tensorproduct(ntuple(t->a, Val{N})...)

# Use \otimes as notation for tensor product.
⊗ = tensorproduct

# Flatten a sequence of elements that may be recursively composite
# For example: a ProductDomain of ProductDomains will yield a list of each of the
# individual domains, like the leafs of a tree structure.
function flatten(::Type{T}, elements::Array, BaseType = Any) where {T}
    flattened = BaseType[]
    for element in elements
        append_flattened!(T, flattened, element)
    end
    flattened
end

flatten(T, elements...) = tuple(flatten(T, [el for el in elements])...)

function append_flattened!(::Type{T}, flattened::Vector, element::T) where {T}
    for el in components(element)
        append_flattened!(T, flattened, el)
    end
end

function append_flattened!(::Type{T}, flattened::Vector, element) where {T}
    append!(flattened, [element])
end

# All tensor products are created using the generic 'tensorproduct' function.
# This function calls a suitable constructor for the tensor product.

for (BaseType,TPType) in [(:DictionaryOperator,:TensorProductOperator),
           (:Dictionary,:TensorProductDict)]
    # In order to avoid strange nested structures, we flatten the arguments
    @eval tensorproduct(args::$BaseType...) = $TPType(flatten($TPType, args...)...)
    @eval tensorproduct(arg::$BaseType, n::Int) = tensorproduct([arg for i in 1:n]...)
    # Disallow tensor products with just one argument
    @eval tensorproduct(arg::$BaseType) = arg

end

function ishomogeneous(tp)
    T = typeof(component(tp, 1))
    mapreduce(el->typeof(el)==T, &, components(tp))
end

basetype(tp) = promote_type(map(typeof, components(tp))...)


text:
module ConstrainedControl

using ConstrainedDynamics
using ConstrainedDynamics: svcat, szeros
using LinearAlgebra
using StaticArrays
using Rotations
using Rotations: rotation_error

export PID,
    LQR,
    TrackingLQR,

    ilqr

include(joinpath("util","util.jl"))
include(joinpath("control","pid.jl"))
include(joinpath("control","lqr.jl"))
include(joinpath("control","lqr_tracking.jl"))
# include(joinpath("control","ilqr.jl"))

end


text:
#!/usr/local/bin/julia
# last tested Julia version: 1.6.1
# Fith-Hugh Nagumo model on a 2D lattice
# FvW 03/2018

using NPZ
using PyCall
using PyPlot
using Statistics
using VideoIO
@pyimport matplotlib.animation as anim

function fhn2d(N, T, t0, dt, s, D, a, b, c, I0, stim, blocks)
    c1 = 1/c
    # initialize Fitz-Hugh Nagumo system
    v = zeros(Float64,N,N)
    w = zeros(Float64,N,N)
    dv = zeros(Float64,N,N)
    dw = zeros(Float64,N,N)
    s_sqrt_dt = s*sqrt(dt)
    X = zeros(Float64,T,N,N)
    # stimulation protocol
    I = zeros(Float64,t0+T,N,N)
    for st in stim
        t_on, t_off = st[1]
        x0, x1 = st[2]
        y0, y1 = st[3]
        I[t0+t_on:t0+t_off, x0:x1, y0:y1] .= I0
    end
    # iterate
    for t in range(2, stop=t0+T, step=1)
        (t%100 == 0) && print("    t = ", t, "\r")
        # FHN equations
        dv = c1*(v - 1/3*v.*v.*v - w + I[t,:,:]) + D*L(v)
        dw = c*(v - a*w .+ b)
        # Ito stochastic integration
        v += (dv*dt + s_sqrt_dt*randn(N,N))
        w += (dw*dt)
        # dead block(s):
        for bl in blocks
            #x0, x1 = bl[1]
            #y0, y1 = bl[2]
            v[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
            w[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
        end
        (t > t0) && (X[t-t0,:,:] = v)
    end
    println("\n")
    return X
end

function animate_pyplot(fname, data)
    """
    Animate 3D array as .mp4 using PyPlot, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    vmin = minimum(data)
    vmax = maximum(data)
    # setup animation image
    fig = figure(figsize=(6,6))
    axis("off")
    t = imshow(data[1,:,:], origin="lower", cmap=ColorMap("gray"),
			   vmin=vmin, vmax=vmax)
    tight_layout()
    # frame generator
    println("[+] animate")
    function animate(i)
        (i%100 == 0) && print("    t = ", i, "/", nt, "\r")
        t.set_data(data[i+1,:,:])
    end
    # create animation
    ani = anim.FuncAnimation(fig, animate, frames=nt, interval=10)
    println("\n")
    # save animation
    ani[:save](fname, bitrate=-1,
               extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    show()
end

function animate_video(fname, data)
    """
    Animate 3D array as .mp4 using VideoIO, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    # BW
    y = UInt8.(round.(255*(data .- minimum(data))/(maximum(data)-minimum(data))))
    # BW inverted
    #y = UInt8.(round.(255 .- 255*(data .- minimum(data))/(maximum(data)-minimum(data))))
    encoder_options = (color_range=2, crf=0, preset="medium")
    framerate=30
    T = size(data,1)
    open_video_out(fname, y[1,end:-1:1,:], framerate=framerate,
                   encoder_options=encoder_options) do writer
        for i in range(2,stop=T,step=1)
            write(writer, y[i,end:-1:1,:])
        end
    end
end

function L(x)
    # Laplace operator
    # periodic boundary conditions
    xU = circshift(x, [-1 0])
    xD = circshift(x, [1 0])
    xL = circshift(x, [0 -1])
    xR = circshift(x, [0 1])
    Lx = xU + xD + xL + xR - 4x
    # non-periodic boundary conditions
    Lx[1,:] .= 0.0
    Lx[end,:] .= 0.0
    Lx[:,1] .= 0.0
    Lx[:,end] .= 0.0
    return Lx
end

function main()
    println("FitzHugh-Nagumo (FHN) lattice model\n")
    N = 128
    T = 1000
    t0 = 0
    dt = 0.1
    s = 0.02 # 0.02 # 0.10
    D = 1.0
    # FitzHugh-Nagumo parameters
    a = 0.5
    b = 0.7
    c = 0.3
    I = 1.0 # 1.0, 0.5
    println("[+] Lattice size N: ", N)
    println("[+] Time steps T: ", T)
    println("[+] Warm-up steps t0: ", t0)
    println("[+] Integration time step dt: ", dt)
    println("[+] Noise std. dev.: ", s)
    println("[+] Diffusion coefficient D: ", D)
    println("[+] FHN parameter a: ", a)
    println("[+] FHN parameter b: ", b)
    println("[+] FHN parameter c: ", c)
    println("[+] Stimulation current I: ", I)

    # stim protocol, array of elements [[t0,t1], [x0,x1], [y0,y1]]
    n_2 = Int(round(N/2)) # 1/2 lattice size
    n_4 = Int(round(N/4)) # 1/4 lattice size
    n_5 = Int(round(N/5)) # 1/5 lattice size
    stim = [ [[25,50], [1,N], [3,8]],
             [[130,150], [n_2-2,n_2+2], [10,25]] ]
    #stim = []

    # dead blocks, array of elementy [[x0,x1], [y0,y1]]
    blocks = [ [[2*n_4,3*n_4], [15,20]], [[2*n_4+10,3*n_4+10], [40,45]] ]
    #blocks = []

    # run simulation
    data = fhn2d(N, T, t0, dt, s, D, a, b, c, I, stim, blocks)
    println("[+] Data dimensions: ", size(data))

    # plot mean voltage
    m = mean(reshape(data, (T,N*N)), dims=2)
    plot(m, "-sk"); show()

    # save data
    I_str = rpad(I, 4, '0') # stim. current amplitude as 4-char string
    s_str = rpad(s, 4, '0') # noise as 4-char string
    D_str = rpad(D, 4, '0') # diffusion coefficient as 4-char string
    fname1 = string("fhn2d_I_", I_str, "_s_", s_str, "_D_", D_str, ".npy")
    #npzwrite(data_filename, data)
    #println("[+] Data saved as: ", data_filename)

    # video
    fname2 = string("fhn2d_I_", I_str, "_s_", s_str, "_D_", D_str, ".mp4")
    #animate_pyplot(fname2, data) # slow
    animate_video(fname2, data) # fast
    println("[+] Data saved as: ", fname2)
end

main()


text:
struct SDEDiffusionTermWrapper{iip,G,H}
  g::G
  h::H
end
function SDEDiffusionTermWrapper(g::G,h::H) where {G,H}
  SDEDiffusionTermWrapper{isinplace(g,5),G,H}(g,h)
end
(g::SDEDiffusionTermWrapper{true})(du, u, p, t) = g.g(du, u, g.h, p, t)
(g::SDEDiffusionTermWrapper{false})(u, p, t) = g.g(u, g.h, p, t)

struct SDEFunctionWrapper{iip,F,G,H,TMM,Ta,Tt,TJ,JP,TW,TWt,TPJ,S,TCV} <: DiffEqBase.AbstractRODEFunction{iip}
  f::F
  g::G
  h::H
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jac_prototype::JP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

(f::SDEFunctionWrapper{true})(du, u, p, t) = f.f(du, u, f.h, p, t)
(f::SDEFunctionWrapper{false})(u, p, t) = f.f(u, f.h, p, t)

function wrap_functions_and_history(f::SDDEFunction, g, h)
  gwh = SDEDiffusionTermWrapper{isinplace(g,5),typeof(g),typeof(h)}(g,h)

  if f.jac === nothing
    jac = nothing
  else
    if isinplace(f)
      jac = let f_jac = f.jac, h = h
        (J, u, p, t) -> f_jac(J, u, h, p, t)
      end
    else
      jac = let f_jac = f.jac, h = h
        (u, p, t) -> f_jac(u, h, p, t)
      end
    end
  end

  SDEFunctionWrapper{isinplace(f),typeof(f.f),typeof(gwh),typeof(h),typeof(f.mass_matrix),
                     typeof(f.analytic),typeof(f.tgrad),typeof(jac),
                     typeof(f.jac_prototype),typeof(f.Wfact),typeof(f.Wfact_t),
                     typeof(f.paramjac),typeof(f.syms),typeof(f.colorvec)}(
                       f.f, gwh, h, f.mass_matrix, f.analytic, f.tgrad, jac,
                       f.jac_prototype, f.Wfact, f.Wfact_t, f.paramjac, f.syms,
                       f.colorvec), gwh
end


text:
function randomPointsNearVertices(V, point, NumPoints, radioIncidencia)

    numVertices = size(V,1);

    V_rand = poly2D.randomPointInPolygon(V, NumPoints)
    flagMat = zeros(NumPoints,numVertices)
    dist = sqrt.(sum( (V_rand - repeat(point,NumPoints,1)).^2, dims=2))
    flagIn = dist .<= radioIncidencia

    V_out = V_rand[flagIn[:] .== true, :]

    return V_out

end


text:
using Test
using ReversibleSeismic
using NiLang

"""
the reversible loss
"""
@i function i_loss!(out::T, param, src, srcv::AbstractVector{T}, c::AbstractMatrix{T},
          tu::AbstractArray{T,3}, tφ::AbstractArray{T,3}, tψ::AbstractArray{T,3}) where T
     i_solve!(param, src, srcv, c, tu, tφ, tψ)
     for i=1:size(tu, 1)
        for j=1:size(tu, 2)
            out += tu[i,j,end] ^ 2
        end
     end
end

"""
obtain gradients with NiLang.AD
"""
function getgrad(c::AbstractMatrix{T}; nstep::Int) where T
     param = AcousticPropagatorParams(nx=size(c,1)-2, ny=size(c,2)-2,
          Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep)

     c = copy(c)
     tu = zeros(T, size(c)..., nstep+1)
     tφ = zeros(T, size(c)..., nstep+1)
     tψ = zeros(T, size(c)..., nstep+1)

     src = size(c) .÷ 2 .- 1
     srcv = Ricker(param, 100.0, 500.0)
     NiLang.AD.gradient(Val(1), i_loss!, (0.0, param, src, srcv, c, tu, tφ, tψ))[end-3]
end

"""
obtain gradients numerically, for gradient checking.
"""
function getngrad(c::AbstractMatrix{T}, i, j; nstep::Int, δ=1e-4) where T
     param = AcousticPropagatorParams(nx=size(c,1)-2, ny=size(c,2)-2,
          Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep)

     c = copy(c)
     tu = zeros(T, size(c)..., nstep+1)
     tφ = zeros(T, size(c)..., nstep+1)
     tψ = zeros(T, size(c)..., nstep+1)

     src = size(c) .÷ 2 .- 1
     srcv = Ricker(param, 100.0, 500.0)
     c[i,j] += δ
     fpos = i_loss!(0.0, param, src, srcv, copy(c), copy(tu), copy(tφ), copy(tψ))[1]
     c[i,j] -= 2δ
     fneg = i_loss!(0.0, param, src, srcv, copy(c), copy(tu), copy(tφ), copy(tψ))[1]
     return (fpos - fneg)/2δ
end

@testset "loss" begin
     nx = ny = 100
     nstep = 2000
     param = AcousticPropagatorParams(nx=nx, ny=ny,
          Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep)

     tu = zeros(nx+2, ny+2, nstep+1)
     tφ = zeros(nx+2, ny+2, nstep+1)
     tψ = zeros(nx+2, ny+2, nstep+1)

     c = 1000*ones(nx+2, ny+2)
     src = (nx ÷ 2, ny ÷ 2)
     srcv = Ricker(param, 100.0, 500.0)

     loss = i_loss!(0.0, param, src, srcv, c, tu, tφ, tψ)[1]
     @test check_inv(i_loss!, (0.0, param, src, srcv, c, tu, tφ, tψ); atol=1e-6)
     #@test loss ≈ 10.931466822080788
     @test loss ≈ 6.234873084873294e-5
end

@testset "gradient" begin
     nx = ny = 100
     nstep = 2000
     c = 1000*ones(nx+2, ny+2)
     g4545 = getgrad(c; nstep=nstep)[45,45]
     ng4545 = getngrad(c, 45, 45; nstep=nstep, δ=1e-4)
     @test isapprox(g4545, ng4545; rtol=1e-2)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Solids (Modia3D/Solids/_module.jl)
#
#
# Content:
#
# Solid geometries and operations on the geometries:
#    volume(geo)                      - Returns volume of solid
#    centroid(geo)                    - Return position vector from solid reference frame to centroid of solid
#    inertiaMatrix(geo,mass)          - Return inertia matrix of solid with respect to reference frame
#    supportPoint(geo,r_abs,R_abs,e)    - Return support point in direction of unit vector e
#    boundingBox!(geo,r_abs,R_abs,AABB) - Update vertex positions of smallest AABB box that contains geo
#
# Useful web pages:
#    https://en.wikipedia.org/wiki/List_of_moments_of_inertia
#    https://en.wikipedia.org/wiki/List_of_second_moments_of_area

using LinearAlgebra
EYE3() = Matrix(1.0I,3,3)

const InertiaMatrix  = SMatrix{3,3,Float64,9}

#------------------------------------ Sphere ----------------------------------
struct SolidSphere <: Modia3D.AbstractSolidGeometry
  Dx::Float64
  function SolidSphere(Dx::Float64; rsmall=0.001)    # rsmall is ignored
    @assert(Dx >= 0.0)
    new(Dx)
  end
end

#------------------------------------ Ellipsoid -------------------------------
struct SolidEllipsoid <: Modia3D.AbstractSolidGeometry
  Dx::Float64
  Dy::Float64
  Dz::Float64
  function SolidEllipsoid(Dx::Float64,Dy::Float64,Dz::Float64; rsmall=0.001)   # rsmall is ignored
    @assert(Dx >= 0.0)
    @assert(Dy  >= 0.0)
    @assert(Dz >= 0.0)
    new(Dx,Dy,Dz)
  end
end

#------------------------------------ Box ------------------------------------
struct SolidBox <: Modia3D.AbstractSolidGeometry
  Lx::Float64
  Ly::Float64
  Lz::Float64
  rsmall::Float64   # Radius of small sphere to smooth surface for collision handling
  function SolidBox(Lx::Float64,Ly::Float64,Lz::Float64; rsmall=0.001)
    @assert(Lx >= 0.0)
    @assert(Ly >= 0.0)
    @assert(Lz >= 0.0)
    @assert(rsmall >= 0.0)

    rsmall2 = min(rsmall, 0.1*min(Lx, Ly, Lz)) # at most 10% of the smallest edge length
    if rsmall2 < rsmall
        println("rsmall in SolidBox changed from $rsmall to $rsmall2")
    end
    @assert(Lx >= 2*rsmall2)
    @assert(Ly >= 2*rsmall2)
    @assert(Lz >= 2*rsmall2)

    new(Lx,Ly,Lz,rsmall2)
  end
end

#------------------------------------ Cylinder --------------------------------
struct SolidCylinder <: Modia3D.AbstractSolidGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  rsmall::Float64   # Radius of small sphere to smooth surface for collision handling
  function SolidCylinder(Dx::Float64, Lz::Float64; Dy=Dx, rsmall=0.001)
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(rsmall >= 0.0)
    rsmall2 = min(rsmall, 0.1*min(Dx, Dy, Lz)) # at most 10% of the smallest edge length
    new(Dx,Dy,Lz,rsmall2)
  end
end

#------------------------------------ Capsule --------------------------------
struct SolidCapsule <: Modia3D.AbstractSolidGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  rsmall::Float64   # Radius of small sphere to smooth surface for collision handling
  function SolidCapsule(Dx::Float64, Lz::Float64; Dy=Dx, rsmall=0.001)
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(rsmall >= 0.0)
    rsmall2 = min(rsmall, 0.1*min(Dx, Dy, Lz)) # at most 10% of the smallest edge length
    new(Dx,Dy,Lz,rsmall2)
  end
end

#------------------------------------ Beam -----------------------------------
struct SolidBeam <: Modia3D.AbstractSolidGeometry
  Lx::Float64
  Dy::Float64
  Lz::Float64
  rsmall::Float64   # Radius of small sphere to smooth surface for collision handling
  function SolidBeam(Lx::Float64,Dy::Float64,Lz::Float64; rsmall=0.001)
    #@assert(length >= height)
    @assert(Lx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    rsmall2 = min(rsmall, 0.1*min(Lx, Dy, Lz)) # at most 10% of the smallest edge length
    new(Lx,Dy,Lz,rsmall2)
  end
end

#------------------------------------ Cone -----------------------------------
struct SolidCone <: Modia3D.AbstractSolidGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  extras::MVector{3,Float64}
  rsmall::Float64   # Radius of small sphere to smooth surface for collision handling
  function SolidCone(Dx::Float64, Lz::Float64; Dy=Dx,relativeTipDiameter=0.5,relativeInnerDiameter=0.0, rsmall=0.001)
    # Tip diameter >= 0.0, Tip diameter is relative to base diameter (0 = real cone, 1 = cylinder),
    # Inner diameter any real number, Inner diameter of the cone is proportional to tip size (0..1)
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(relativeTipDiameter>=0.0 && relativeTipDiameter<=1.0)
    @assert(relativeInnerDiameter>=0.0 && relativeInnerDiameter<1.0)
    @assert(rsmall >= 0.0)
    rsmall2 = min(rsmall, 0.1*min(Dx, Dy, Lz)) # at most 10% of the smallest edge length
    extras=[relativeTipDiameter,relativeInnerDiameter,0.0]
    new(Dx,Dy,Lz,extras,rsmall2)
  end
end

#------------------------------------ Pipe -----------------------------------
struct SolidPipe <: Modia3D.AbstractSolidGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  extras::MVector{3,Float64}
  rsmall::Float64   # Radius of small sphere to smooth surface for collision handling
  function SolidPipe(Dx::Float64, Lz::Float64; Dy=Dx, relativeInnerDiameter=0.9, rsmall=0.001)
    #relativeInnerDiameter ... Inner diameter of pipe (0 = rigid, 1 = fully hollow), every real number is allowed
    #relativeTipDiameter ... Tip diameter of pipe relative to base diameter (0 = real cone, 1 = cylinder), every real number >= 0 is possible
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(relativeInnerDiameter >=0.0 && relativeInnerDiameter<1.0)
    @assert(rsmall >= 0.0)
    relativeTipDiameter=1.0
    rsmall2 = min(rsmall, 0.1*min(Dx, Dy, Lz)) # at most 10% of the smallest edge length
    extras=[0.0,relativeInnerDiameter,relativeTipDiameter]
    new(Dx,Dy,Lz,extras,rsmall2)
  end
end

#------------------------------------ SolidFileMesh -----------------------------------
struct SolidFileMesh <: Modia3D.AbstractSolidGeometry
  filename::AbstractString
  scaleFactor::MVector{3,Float64}
  useGraphicsMaterialColor::Bool
  smoothNormals::Bool
  centroid::MVector{3,Float64}
  longestEdge::Float64
  objPoints::Vector{MVector{3,Float64}}
  facesIndizes::Vector{AbstractVector{}}
  volume::Float64
  centroidAlgo::SVector{3,Float64}
  inertia::SMatrix{3,3,Float64,9}
  function SolidFileMesh(filename::AbstractString; scaleFactor=Basics.onesMVector()  , useGraphicsMaterialColor::Bool=false, smoothNormals::Bool=false)
     if !isfile(filename)
        error("SolidFileMesh(\"$filename\",...): file not found.")
     end
     # @assert(isfile(filename))
     @assert(scaleFactor[1] >= 0.0)
     @assert(scaleFactor[2] >= 0.0)
     @assert(scaleFactor[3] >= 0.0)
     (centroid, longestEdge, objPoints, facesIndizes) = getObjInfos(filename, scaleFactor)
     (volume, centroidAlgo, inertia) = computeMassProperties(objPoints, facesIndizes; bodyCoords=false)
     #println("volume = ", volume)
     fileMesh = new(filename, scaleFactor, useGraphicsMaterialColor, smoothNormals, centroid, longestEdge, objPoints, facesIndizes, volume, centroidAlgo, inertia)
     return fileMesh
  end
end
SolidFileMesh(filename::AbstractString, scale::Number; useGraphicsMaterialColor::Bool=false, smoothNormals::Bool=false) =
        SolidFileMesh(filename; scaleFactor=[scale,scale,scale], useGraphicsMaterialColor=useGraphicsMaterialColor, smoothNormals=smoothNormals)

#------------ bottom area ----------------------------------
bottomArea(geo::SolidBox)       = geo.Lx*geo.Ly
bottomArea(geo::SolidCylinder)  = pi*geo.Dx/2*geo.Dy/2
function bottomArea(geo::SolidBeam)
  w = geo.Dy
  h = geo.Lx
  # Area = rectangle + circle
  return h*w + pi*(w/2)^2
end
function bottomArea(geo::SolidCone)
  A = geo.Dx/2
  B = geo.Dy/2
  # maybe for inner hole (circle)
  a_inner = geo.Dx/2 * geo.extras[1] * geo.extras[2]
  b_inner = geo.Dy/2 * geo.extras[1] * geo.extras[2]
  if geo.extras[1] == 0.0   # total cone
    return pi*A*B
  else  # frustum of a cone, maybe with a hole (inner circle)
    return pi*A*B - pi*a_inner*b_inner
  end
end
function bottomArea(geo::SolidPipe)
  A = geo.Dx/2
  B = geo.Dy/2
  a = geo.Dx/2*geo.extras[2]
  b = geo.Dy/2*geo.extras[2]
  return pi*(A*B - a*b)
end
bottomArea(geo::Modia3D.AbstractSolidGeometry)  = error(typeof(geo), ": has no bottom or top area!")

#------------top area ----------------------------------
function topArea(geo::SolidCone)
  # top area
  a = geo.Dx/2 * geo.extras[1]
  b = geo.Dy/2 * geo.extras[1]
  # for cylinder (also top area)
  a_inner = geo.Dx/2 * geo.extras[1] * geo.extras[2]
  b_inner = geo.Dy/2 * geo.extras[1] * geo.extras[2]
  if geo.extras[1] > 0.0   # total cone
    return pi*a*b - pi*a_inner*b_inner
  end
end
topArea(geo::Modia3D.AbstractSolidGeometry) = bottomArea(geo)

#------------------volume-----------------------------------------------
"""
    V = volume(geo)

Return the volume of the solid geometry `geo::Modia3D.AbstractSolidGeometry` in [m^3].
"""
volume(geo::Modia3D.AbstractSolidGeometry) = bottomArea(geo)*lengthGeo(geo)
volume(geo::SolidSphere)                   = 4/3*pi*(geo.Dx/2)^3
function volume(geo::SolidEllipsoid)
  a = geo.Dx/2
  b = geo.Dy/2
  c = geo.Dz/2
  return 4/3*pi*a*b*c
end
function volume(geo::SolidCapsule)
  a = geo.Dx/2
  b = geo.Dy/2
  h = geo.Lz
  # volume cylinder: pi*h*a*b
  # volume ellipsoid: only changing in x and y direction: 4/3*pi*b^2*a
  return pi*h*a*b + 4/3*pi*b^2*a
end
function volume(geo::SolidCone)
  h = geo.Lz
  # bottom area
  A = geo.Dx/2
  B = geo.Dy/2
  # top area
  a = geo.Dx/2 * geo.extras[1]
  b = geo.Dy/2 * geo.extras[1]
  # for cylinder (also top area)
  a_inner = geo.Dx/2 * geo.extras[1] * geo.extras[2]
  b_inner = geo.Dy/2 * geo.extras[1] * geo.extras[2]
  if geo.extras[1] == 0.0   # total cone
    return (h*pi/3)*A*B
  else    # frustum of a cone (with elliptic ground area): (h*pi/3) * (A*B^2 - a*b^2)/(B - b)
    # solid frustum of a cone - cylinder
    return (h*pi/3) * (A*B^2 - a*b^2)/(B - b) - pi*h*a_inner*b_inner
  end
end
function volume(geo::SolidFileMesh)
    if !isempty(geo.facesIndizes)
    #  println("geo.volume = ", geo.volume)
      return geo.volume
    else
      println("SolidFileMesh: ", geo.filename, ". The surface areas must be triangular, and each triangle should be specified in right-handed/counter-clockwise order. Otherwise it is not possible to compute a volume.")
      return nothing
    end
end

#--------------------longest edge (maximum of all directions)---------------------------------
longestEdge(geo::SolidSphere)    = geo.Dx
longestEdge(geo::SolidEllipsoid) = max(geo.Dx,geo.Dy,geo.Dz)
longestEdge(geo::SolidBox)       = max(geo.Lx,geo.Ly,geo.Lz)
longestEdge(geo::SolidBeam)      = max(geo.Lx,geo.Dy,geo.Lz)
longestEdge(geo::SolidFileMesh)  = geo.longestEdge
longestEdge(geo::Modia3D.AbstractSolidGeometry) = max(geo.Dx,geo.Dy,geo.Lz)

#---------------------length of geometry (length in z direction)----------------------------------------------------------
# length of geometries
lengthGeo(geo::Modia3D.AbstractSolidGeometry) = geo.Lz
lengthGeo(geo::SolidSphere)                   = geo.Dx
lengthGeo(geo::SolidEllipsoid)                = geo.Dz # error("There is no length for an ellipsoid!")
lengthGeo(geo::SolidFileMesh)                 = error("lengthGeo(SolidFileMesh): It is not implemented yet!")

#------------------------centroid of geometries----------------------------------------
"""
    r = centroid(geo)

Return position vector from solid reference frame to [centroid](https://en.wikipedia.org/wiki/Centroid)
of solid `geo::Modia3D.AbstractSolidGeometry` in [m]. If the solid has a uniform density,
the centroid is identical to the *center of mass*.
"""
centroid(geo::Modia3D.AbstractSolidGeometry) = ModiaMath.ZeroVector3D
centroid(geo::SolidCone)                     = SVector{3,Float64}([0.0,0.0,geo.Lz/4])
function centroid(geo::SolidFileMesh)
    if !isempty(geo.facesIndizes)
        return geo.centroidAlgo
    else
        return geo.centroid
    end
end

#------------------------inertia matrix of geometries-----------------------------------
"""
   I = inertiaMatrix(geo, mass)

Return [inertia matrix](https://en.wikipedia.org/wiki/Moment_of_inertia) `I` of solid
`geo::Modia3D.AbstractSolidGeometry` with respect to the
solid reference frame in [kg*m^2] as `SMatrix{3,3,Float64,9}`. Hereby it is assumed that
`geo` has uniform density and `mass` is the mass of `geo` in [kg].
"""
inertiaMatrix(geo::SolidSphere, mass::Number) = InertiaMatrix(mass/10*geo.Dx^2*EYE3())
inertiaMatrix(geo::SolidEllipsoid, mass::Number) = InertiaMatrix(mass/20*Diagonal([geo.Dy^2 + geo.Dz^2,
                                                                                geo.Dx^2 + geo.Dz^2,
                                                                                geo.Dx^2 + geo.Dy^2]))
inertiaMatrix(geo::SolidBox, mass::Number) = InertiaMatrix(1/12*mass*Diagonal([geo.Ly^2 + geo.Lz^2,
                                                                               geo.Lx^2 + geo.Lz^2,
                                                                               geo.Lx^2 + geo.Ly^2]))
function inertiaMatrix(geo::SolidCylinder, mass::Number)
    #println("mass = ", mass)
    return InertiaMatrix(mass*Diagonal([1/4*(geo.Dy/2)^2 + 1/3*geo.Lz^2,
                                        1/4*(geo.Dx/2)^2 + 1/3*geo.Lz^2,
                                        1/4*((geo.Dx/2)^2 + (geo.Dy/2)^2)]))
end
function inertiaMatrix(geo::SolidCapsule, massGeo::Number)
    @warn "from Modia3D.inertiaMatrix(SolidCapsule): inertia matrix is not fully tested yet!"
    # mass = rho * volume
    #=
    volGeo = volume(geo)
    rho = massGeo/vollGeo
    volBox = geo.Lx*geo.Lz*geo.Dy
    massBox = volBox*rho
    =#

    volGeo = volume(geo)
    rho = massGeo/vollGeo

    volCyl = geo.Lz*(geo.Dx/2)*(geo.Dy/2)*pi
    massCyl = volCyl*rho

    volHalfEll = (4/3*pi*(geo.Dy/2)^2*(geo.Dx/2))/2
    massHalfEll = volHalfEll*rho

    IxxCyl = massCyl*(1/4*(geo.Dy/2)^2 + 1/3*(geo.Lz^2))
    IyyCyl = massCyl*(1/4*(geo.Dx/2)^2 + 1/3*(geo.Lz^2))
    IzzCyl = massCyl*(1/4*((geo.Dx/2)^2 + (geo.Dy/2)^2))

    # SolidEllipsoid: with Dx, Dy, Dz = Dy
    IxxEll = massHalfEll*1/20*((geo.Dy/2)^2 + (geo.Dy/2)^2)
    IyyEll = massHalfEll*1/20*((geo.Dx/2)^2 + (geo.Dy/2)^2)
    IzzEll = massHalfEll*1/20*((geo.Dx/2)^2 + (geo.Dy/2)^2)

    IxxSteiner = IxxEll + (geo.Lz/2)^2*massHalfEll
    IyySteiner = IyyEll + (geo.Lz/2)^2*massHalfEll

    Ixx    = IxxCyl + 2*IxxSteiner
    Iyy    = IyyCyl + 2*IyySteiner
    Izz    = IzzCyl + 2*IzzEll
    return InertiaMatrix(Diagonal([Ixx,Iyy,Izz]))
end
function inertiaMatrix(geo::SolidBeam, massGeo::Number)
    # mass = rho * volume
    volGeo = volume(geo)
    rho = massGeo/volGeo
    volBox = geo.Lx*geo.Lz*geo.Dy
    massBox = volBox*rho

    volHalfCyl  = geo.Lz*(geo.Dy/2)^2*pi/2
    massHalfCyl = volHalfCyl*rho

    IxxBox = 1/12*massBox * (geo.Dy^2 + geo.Lz^2)
    IyyBox = 1/12*massBox * (geo.Lx^2 + geo.Lz^2)
    IzzBox = 1/12*massBox * (geo.Lx^2 + geo.Dy^2)

    # http://www.efunda.com/math/solids/solids_display.cfm?SolidName=HalfCircularCylinder
    IxxCylHalf = massHalfCyl*(1/4*(geo.Dy/2)^2 + 1/3*(geo.Lz^2))
    IyyCylHalf = massHalfCyl*(1/4*(geo.Dy/2)^2 + 1/3*(geo.Lz^2))
    IzzCylHalf = massHalfCyl*(1/2*(geo.Dy/2)^2)

    IyyCylHalfSteiner = IyyCylHalf + massHalfCyl*(geo.Lx/2)^2
    IzzCylHalfSteiner = IzzCylHalf + massHalfCyl*(geo.Lx/2)^2

    Ixx = IxxBox + 2*IxxCylHalf
    Iyy = IyyBox + 2*IyyCylHalfSteiner
    Izz = IzzBox + 2*IzzCylHalfSteiner
    return InertiaMatrix(Diagonal([Ixx,Iyy,Izz]))
end
function inertiaMatrix(geo::SolidCone, mass::Number)
    if geo.Dx == geo.Dy && geo.extras[2] == 0.0
        r = geo.Dx/2
        h = geo.Lz
        if geo.extras[1] == 0.0   # total cone
            # https://en.wikipedia.org/wiki/List_of_moments_of_inertia
            return InertiaMatrix(3*mass*Diagonal([1/20*r^2 + 1/5*h^2,
                                                  1/20*r^2 + 1/5*h^2,
                                                  1/10*r^2]))
        else
            # https://de.wikipedia.org/wiki/Tr%C3%A4gheitsmoment 18.7.18
            # https://matheplanet.com/default3.html?call=viewtopic.php?topic=202081&ref=https%3A%2F%2Fwww.bing.com%2F 18.7.18
            h2 = geo.Lz*geo.extras[1]
            r2 = geo.Dx/2
            r1 = geo.Dx/2*geo.extras[1]
            xks = h2/4*((r2^2+2*r1*r2+3*r1^2)/(r2^2+r1*r2+r1^2))
            PHI = 1/(r2^3-r1^3)*( 3/80*(4+h2^2/(r2-r1)^2)*(r2^5-r1^5) + r2^3*(h2/4*r2/(r2-r1)-xks)^2 - r1^3*(h2/4*(4*r2-3*r1)/(r2-r1)-xks)^2 )
            return InertiaMatrix(mass*Diagonal([PHI,
                                                PHI,
                                                3/10*(r2^5-r1^5)/(r2^3-r1^3)]))
        end
    else
        error("SolidCone: at the moment inertiaMatrix is defined for Dx=Dy and without an inner hole.")
    end
end
function inertiaMatrix(geo::SolidPipe, mass::Number)
    if geo.Dx == geo.Dy
        r1 = geo.Dx/2*geo.extras[2]
        r2 = geo.Dx/2
        h  = geo.Lz
        return InertiaMatrix(mass*Diagonal([1/12*(3*(r1^2+r2^2)^2 + h^2),
                                            1/12*(3*(r1^2+r2^2)^2 + h^2),
                                            1/2*(r1^2+r2^2)^2]))
    else
        r1x = geo.Dx/2*geo.extras[2]
        r2x = geo.Dx/2
        r1y = geo.Dy/2*geo.extras[2]
        r2y = geo.Dy/2
        h  = geo.Lz
        return InertiaMatrix(mass*Diagonal([1/12*(3*(r1y^2+r2y^2)^2 + h^2),
                                            1/12*(3*(r1x^2+r2x^2)^2 + h^2),
                                            1/4*(r1x^2+r2x^2)^2 + 1/4*(r1y^2+r2y^2)^2]))
    end
end
function inertiaMatrix(geo::SolidFileMesh, mass::Number)
  #println("mass = ", mass)
    if !isempty(geo.facesIndizes)
        return geo.inertia.*mass
    else
        println("SolidFileMesh: ", geo.filename, ". The surface areas must be triangular, and each triangle should be specified in right-handed/counter-clockwise order. Otherwise it is not possible to compute an inertia tensor.")
        return nothing
    end
end

include("computePropertiesFileMes.jl")

include("boundingBoxes.jl")

include("utilities.jl")


text:
println("Number of processes running: ",nprocs())

#Define some quantities and a very simple TargetModel for all tests below
rinitseed1 = 879
rnburnin1 = 2
rniter1 = 4
smhnprops1 = 1 #number of test proposals for standard mh
gmhnprops1 = 2 #number of test proposals for generalized mh
rtime1 = linspace(0.0,10.0,100)
rlower1 = [2.5]
rupper1 = [3.5]
rdefault1 = [3.0]
rcov1 = [0.01]
rparas1 = parameters([:a],rlower1,rupper1,rdefault1)
rdata1 = data(:array,rtime1,sin(3*rtime1))
rnoise1 = noise(:gaussian,rcov1)
rmodel1 = model(:target,rparas1,rdata1,rnoise1,(t,p)->sin(p[1]*t);name="RunnersTestModel")
rnparas1 = numparas(rmodel1)
rsampler1 = sampler(:mh,:normal,0.1,eye(1))
rtuner1 = tuner(:monitor,1)
rtuner2 = tuner(:scale,2,0.5,:logistic)

function testfrom(state,vals,ll)
    @test state.from.values == vals
    @test state.from.loglikelihood == ll
end

function testproposals(state,vals,ll)
    @test state.proposals.values == vals
    @test state.proposals.loglikelihood == ll
end

function teststore(chain,vals,ll,i,a,p)
    @test chain.values[:,i] == vals
    @test chain.loglikelihood[i] == ll
    @test chain.accepted == a
    @test chain.proposed == p
end

function setfrom(state,vals,loglikelihood)
    copy!(state.from.values,vals)
    copy!(state.from.loglikelihood,loglikelihood)
end

function setproposals(state,vals,loglikelihood)
    copy!(state.proposals.values,vals)
    copy!(state.proposals.loglikelihood,loglikelihood)
end

function setindicator(indicator,probs,samples)
    copy!(indicator.stationary,probs)
    copy!(indicator.samples,samples)
end

function setsegment(rsegs,from,ll1,prop,ll2,j)
    l = length(prop)
    p,s = GeneralizedMetropolisHastings._prop2seg(rsegs,j)
    r = rsegs.remote[s]
    @sync begin
        @spawnat r.where copy!(fetch(r).samplerstate.from.values,from)
        @spawnat r.where copy!(fetch(r).samplerstate.from.loglikelihood,ll1)
        @spawnat r.where copy!(fetch(r).samplerstate.proposals.values,l*(p-1)+1,prop,1,l)
        @spawnat r.where copy!(fetch(r).samplerstate.proposals.loglikelihood,p,ll2,1,1)
    end
    rsegs.prop2collected[j] = (p,s)
end

function testtransitionprobability(states,indicator)
    acc = Float64[]
    for x in states
        acc = vcat(acc,acceptanceratio(x))
    end
    acc = vcat(acc,[0.0])
    acc = exp(acc - maximum(acc))
    acc /= sum(acc)
    @test_approx_eq indicator.stationary acc
end

function testupdatefrom(indicatorstate,states,indicator)
    indend = indicatorsamples(indicator)[end]
    if indend != numproposals(indicator) + 1
        s,t = ind2sub(GeneralizedMetropolisHastings._ind2subtuple(states),indend)
        @test from(indicatorstate).values == proposals(states[t]).values[:,s]
        @test from(indicatorstate).loglikelihood[1] == proposals(states[t]).loglikelihood[s]
        @test from(indicatorstate).logprior[1] == proposals(states[t]).logprior[s]
    end
end

function dumpiterate1(samplerstates,indicator)
    dump(samplerstates[1])
    println("===============================================")
    dump(indicator1)
    println("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
end

function runiterateandstore1!(runner,model,samplerstates,indicator,chain)
    prevacc = chain.accepted
    prevprop = chain.proposed
    iterateandstore!(runner,model,samplerstates,indicator,chain)
    teststore(samplerstates[1],samplerstates,indicator,chain,prevacc,prevprop)
    testupdatefrom(samplerstates1[1],samplerstates1,indicator1)
    dumpiterate1(samplerstates,indicator)
end

function dumpiterate2(indicatorstate,samplerstates,indicator)
    dump(indicatorstate)
    println("000000000000000000000000000000000000000000000000000000000000000000000000000")
    dump(samplerstates[1])
    println("111111111111111111111111111111111111111111111111111111111111111111111111111")
    dump(samplerstates[2])
    println("222222222222222222222222222222222222222222222222222222222222222222222222222")
    dump(samplerstates[3])
    println("333333333333333333333333333333333333333333333333333333333333333333333333333")
    dump(indicator)
    println("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
end

function runiterateandstore2!(runner,model,indicatorstate,samplerstates,indicator,chain)
    prevacc = chain.accepted
    prevprop = chain.proposed
    iterateandstore!(runner,model,indicatorstate,samplerstates,indicator,chain)
    teststore(indicatorstate,samplerstates,indicator,chain,prevacc,prevprop)
    testupdatefrom(indicatorstate,samplerstates,indicator)
    dumpiterate2(indicatorstate,samplerstates,indicator)
end


text:
using MLDatasets: MNIST
using Flux: onehotbatch
using ProtoGrad: Conv, Linear, Dropout, Compose, maxpool, relu, softmax
using ProtoGrad: SupervisedObjective, cross_entropy, class_error
using ProtoGrad: Adam, forever, within_training_mode
using ProgressMeter: @showprogress
using StatsBase: sample
using Serialization: serialize, deserialize

train_x, train_y = MNIST.traindata(Float32, dir=joinpath(".", "datasets", "MNIST"))
test_x, test_y = MNIST.testdata(Float32, dir=joinpath(".", "datasets", "MNIST"))

train_x_with_channel = reshape(train_x, size(train_x)[1:2]..., 1, size(train_x)[3])
train_y_onehot = onehotbatch(train_y, 0:9)

input_model_filename = "./serialized_conv_model"
println("loading model from $(input_model_filename)")
m_original = deserialize(input_model_filename)

function get_complete_model(m_top)
    return Compose(
        m_original.layers[1:end-3]...,
        m_top.layers...
    )
end

m_top = Compose(m_original.layers[end-2:end]...)

batch_size = 128

training_batches = forever() do
    idx = sample(1:size(train_x_with_channel)[end], batch_size, replace = false)
    return (train_x_with_channel[:, :, :, idx], train_y_onehot[:, idx])
end

f = SupervisedObjective(cross_entropy, training_batches) ∘ get_complete_model

optimizer = Adam()
iterations = Iterators.Stateful(optimizer(m_top, f))

num_epochs = 5
num_batches_per_epoch = 500

within_training_mode() do
    @showprogress for epoch in 1:num_epochs
        for output in Iterators.take(iterations, num_batches_per_epoch)
            global m_top = output.solution
        end
    end
end

m = get_complete_model(m_top)

model_filename = "./serialized_conv_finetuned_model"
println("saving model to $(model_filename)")
serialize(model_filename, m)

m_deserialized = deserialize(model_filename)

test_x_with_channel = reshape(test_x, size(test_x)[1:2]..., 1, size(test_x)[3])
test_y_onehot = onehotbatch(test_y, 0:9)

err = class_error(m_deserialized(test_x_with_channel), test_y_onehot)

println("test error = $(err * 100)%")


text:
###############################################################################
# PackageEvaluator
###############################################################################
# output.jl
# Output functions for test results
# Exports:
#  featuresToJSON
#  getDetailsString
###############################################################################

# keyToJSON
# Helper function that takes a key and value and writes them in the official
# JSON format.
function keyToJSON(key, value, last=false)
    return "  \"$key\": \"$value\"$(!last?",":"")\n"
end

# featuresToJSON
# Takes test results and formats them as a JSON string that can be handled
# by status.julialang.org
export featuresToJSON
function featuresToJSON(pkg_name, features)
    json_str = "{\n"
    json_str = json_str * keyToJSON("name",     pkg_name)
    json_str = json_str * keyToJSON("url",      features[:URL])
    json_str = json_str * keyToJSON("version",  features[:VERSION])
    json_str = json_str * keyToJSON("gitsha",   chomp(features[:GITSHA]))
    json_str = json_str * keyToJSON("license",  features[:LICENSE])
    json_str = json_str * keyToJSON("licfile",  features[:LICENSE_FILE])
    json_str = json_str * keyToJSON("status",   features[:TEST_STATUS])
    json_str = json_str * keyToJSON("possible", features[:TEST_POSSIBLE] ? "true" : "false")
    json_str = json_str * keyToJSON("details",  getDetailsString(pkg_name, features))
    json_str = json_str * keyToJSON("pkgreq",   features[:REQUIRE_VERSION] ? "true" : "false")
    json_str = json_str * keyToJSON("metareq",  features[:REQUIRES_OK] ? "true" : "false")
    json_str = json_str * keyToJSON("travis",   features[:TEST_TRAVIS] ? "true" : "false", true)
    json_str = json_str * "}"
    return json_str
end

# getDetailsString
# Take full results and builds a human-readable string that summarizes the
# testing results (TEST_EXIST, TEST_STATUS, TEST_MASTERFILE)
export getDetailsString
function getDetailsString(pkg_name, features)
  t_exist     = features[:TEST_EXIST]
  t_status    = features[:TEST_STATUS]
  t_master    = features[:TEST_MASTERFILE]
  t_possible  = features[:TEST_POSSIBLE]

  if !(t_possible)
    return "This package can't be automatically tested - see package's README."
  end

  details = ""
  if t_exist && t_master == ""
    details = "Tests exist, no master file to run, tried 'using $pkg_name'"
    if t_status == "using_fail"
      details = string(details, ", failed!")
    else
      details = string(details, ", no errors.")
    end
  elseif t_exist && t_master != ""
    details = "Tests exist, ran 'julia $(splitdir(t_master)[2])'"
    if t_status == "full_fail"
      details = string(details, ", failed!")
    else
      details = string(details, ", passed!")
    end
  else
    details = "No tests, tried 'using $pkg_name'"
    if t_status == "using_fail"
      details = string(details, ", failed!")
    else
      details = string(details, ", no errors.")
    end
  end

  return details
end


text:
module PencilIO

using ..PencilArrays
import ..PencilArrays: MaybePencilArrayCollection, collection_size, collection

using MPI
using Requires: @require
using StaticArrays: SVector
using TimerOutputs

"""
    ParallelIODriver

Abstract type specifying a parallel I/O driver.
"""
abstract type ParallelIODriver end

"""
    open([f::Function], driver::ParallelIODriver, filename, comm::MPI.Comm; keywords...)

Open parallel file using the chosen driver.

## Keyword arguments

Supported keyword arguments include:

- open mode arguments: `read`, `write`, `create`, `append` and `truncate`.
  These have the same behaviour and defaults as `Base.open`.
  Some of them may be ignored by the chosen driver (see driver-specific docs).

- as in [`MPI.File.open`](https://juliaparallel.github.io/MPI.jl/latest/io/#MPI.File.open),
  other arguments are passed via an `MPI.Info` object.

Note that driver-specific options (such as HDF5 property lists) must be passed
to each driver's constructor.

## See also

- [`open(::MPIIODriver)`](@ref) for MPI-IO specific options
- [`open(::PHDF5Driver)`](@ref) for HDF5 specific options

"""
function Base.open(::ParallelIODriver) end

function Base.open(f::Function, driver::ParallelIODriver, args...; kw...)
    fid = open(driver, args...; kw...)
    try
        f(fid)
    finally
        close(fid)
    end
end

# Metadata to be attached to each dataset (as HDF5 attributes or in an external
# metadata file).
function metadata(x::MaybePencilArrayCollection)
    pen = pencil(x)
    topo = topology(x)
    edims = extra_dims(x)  # this may be an empty tuple, with no type information
    (
        permutation = Tuple(permutation(x)),
        extra_dims = SVector{length(edims),Int}(edims),
        decomposed_dims = decomposition(pen),
        process_dims = size(topo),
    )
end

function keywords_to_open(; read=nothing, write=nothing, create=nothing,
                          truncate=nothing, append=nothing, other_kws...)
    flags = Base.open_flags(read=read, write=write, create=create,
                            truncate=truncate, append=append)
    flags, other_kws
end

include("mpi_io.jl")

function __init__()
    @require HDF5="f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" @eval include("hdf5.jl")
end

end


text:
update(::S, m::AbstractLedger) where {S<:System}= error("No update method implemented for $S")

requested_components(::System) = ()

const Stage = Pair{Symbol, Vector{System}}

Base.push!(s::Stage, sys) = push!(last(s), sys)

Base.insert!(s::Stage, i::Integer, sys) = insert!(last(s), i, sys)

function requested_components(stage::Stage)
    comps = Type{<:ComponentData}[]
	for s in last(stage)
		for c in requested_components(s)
			push!(comps, c)
		end
	end
	return comps
end

function prepare(s::Stage, m::AbstractLedger)
    for sys in last(s)
        prepare(sys, m)
    end
end

prepare(::System, ::AbstractLedger) = nothing


text:
import Statistics: mean, std, stdm, var, varm

mean(a::KnetArray; o...) = (b = sum(a; o...); b .* convert(eltype(b),(length(b)/length(a))))
mean(f::Base.Callable, a::KnetArray) = sum(f, a) / length(a)
std(x::KnetArray, args...; kws...) = sqrt.(var(x, args...; kws...))
stdm(x::KnetArray, args...; kws...) = sqrt.(varm(x, args...; kws...))
var(x::KnetArray; corrected::Bool=true, mean=nothing, dims=:)=_varm(x, mean; corrected=corrected, dims=dims)
varm(x::KnetArray, m; corrected::Bool=true, dims=:)=_varm(x, m; corrected=corrected, dims=dims)

function _varm(x, m; corrected::Bool=true, dims=:)
    if m === nothing; m = mean(x,dims=dims); end
    s = sum(abs2, x .- m; dims=dims)
    r = length(x) ÷ length(s) - Int(corrected)
    s ./ r
end

if first(methods(varm, (AutoGrad.Value,Any))).module !== AutoGrad # TODO: delete after AutoGrad 1.1.4
    stdm(x::AutoGrad.Value, args...; kws...) = sqrt.(varm(x, args...; kws...))
    varm(x::AutoGrad.Value, m; corrected::Bool=true, dims=:)=_varm(x, m; corrected=corrected, dims=dims)
end

# TODO: stdm, varm to be deprecated after Julia 1.x


text:
export ProductMeasure

using MappedArrays
using Base: @propagate_inbounds
import Base
using FillArrays

abstract type AbstractProductMeasure <: AbstractMeasure end

struct ProductMeasure{F,S,I} <: AbstractProductMeasure
    f::Kernel{F,S}
    pars::I
end

# TODO: Test for equality without traversal, probably by first converting to a
# canonical form
function Base.:(==)(a::ProductMeasure, b::ProductMeasure)
    all(zip(a.pars, b.pars)) do (aᵢ, bᵢ)
        a.f(aᵢ) == b.f(bᵢ)
    end
end

Base.size(μ::ProductMeasure) = size(marginals(μ))

Base.length(m::ProductMeasure{T}) where {T} = length(marginals(μ))

basemeasure(d::ProductMeasure) = productmeasure(basekernel(d.f), d.pars)

# TODO: Do we need these methods?
# basemeasure(d::ProductMeasure) = ProductMeasure(basemeasure ∘ d.f, d.pars)
# basemeasure(d::ProductMeasure{typeof(identity)}) = ProductMeasure(identity, map(basemeasure, d.pars))
# basemeasure(d::ProductMeasure{typeof(identity), <:FillArrays.Fill}) = ProductMeasure(identity, map(basemeasure, d.pars))

export marginals

function marginals(d::ProductMeasure{F,S,I}) where {F,S,I}
    _marginals(d, isiterable(I))
end

function _marginals(d::ProductMeasure, ::Iterable)
    return (d.f(i) for i in d.pars)
end

function _marginals(d::ProductMeasure{F,S,I}, ::NonIterable) where {F,S,I}
    error("Type $I is not iterable. Add an `iterate` or `marginals` method to fix.")
end

testvalue(d::ProductMeasure) = map(testvalue, marginals(d))

function Pretty.tile(μ::ProductMeasure{F,S,NamedTuple{N,T}}) where {F,S,N,T}
    result = Pretty.literal("Product(")
    result *= Pretty.pair_layout(μ.f, μ.pars; sep = ", ")
    result *= Pretty.literal(")")
end

function Base.rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure) where {T}
    _rand(rng, T, d, marginals(d))
end

function _rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure, mar) where {T}
    (rand(rng, T, m) for m in mar)
end

###############################################################################
# I <: Tuple

struct TupleProductMeasure{T} <: AbstractProductMeasure
    pars::T
end

export ⊗
⊗(μs::AbstractMeasure...) = productmeasure(μs)

marginals(d::TupleProductMeasure{T}) where {F,T<:Tuple} = d.pars

function Pretty.tile(μ::TupleProductMeasure{T}) where {F,T<:Tuple}
    mar = marginals(μ)
    Pretty.list_layout(Pretty.Layout[Pretty.tile.(mar)...]; sep = " ⊗ ")
end

@inline function logdensity(d::TupleProductMeasure, x::Tuple) where {T<:Tuple}
    mapreduce(logdensity, +, d.pars, x)
end

function Base.rand(rng::AbstractRNG, ::Type{T}, d::TupleProductMeasure) where {T}
    rand.(d.pars)
end

###############################################################################
# I <: AbstractArray

marginals(d::ProductMeasure{F,S,A}) where {F,S,A<:AbstractArray} = mappedarray(d.f, d.pars)

function marginals(d::ProductMeasure{<:Returns,S,A}) where {F,S,A<:AbstractArray}
    Fill(d.f.f.value, size(d.pars))
    # mappedarray(d.f.f, d.pars)
end

function logdensity(d::ProductMeasure, x)
    mapreduce(logdensity, +, marginals(d), x)
end

function logdensity(d::ProductMeasure{<:Returns}, x)
    sum(x -> logdensity(d.f.f.value, x), x)
end

function Pretty.tile(d::ProductMeasure{F,S,A}) where {F,S,A}
    result = Pretty.literal("For(")
    result *= Pretty.pair_layout(Pretty.tile(d.f), Pretty.tile(d.pars); sep = ", ")
    result *= Pretty.literal(")")
end

###############################################################################
# I <: CartesianIndices

function Pretty.tile(d::ProductMeasure{F,S,I}) where {F,S,I<:CartesianIndices}
    result = Pretty.literal("For(")
    result *= Pretty.pair_layout(Pretty.tile(d.f), Pretty.tile(size(d.pars)); sep = ", ")
    result *= Pretty.literal(")")
end

# function Base.rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure{F,S,I}) where {T,F,I<:CartesianIndices}

# end

###############################################################################
# I <: Base.Generator

export rand!
using Random: rand!, GLOBAL_RNG, AbstractRNG

function logdensity(d::ProductMeasure{F,S,I}, x) where {F,S,I<:Base.Generator}
    sum((logdensity(dj, xj) for (dj, xj) in zip(marginals(d), x)))
end

@propagate_inbounds function Random.rand!(
    rng::AbstractRNG,
    d::ProductMeasure,
    x::AbstractArray
)
    # TODO: Generalize this
    T = Float64
    for (j, m) in zip(eachindex(x), marginals(d))
        @inbounds x[j] = rand(rng, T, m)
    end
    return x
end

export rand!
using Random: rand!, GLOBAL_RNG, AbstractRNG

function _rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure, mar::AbstractArray) where {T}
    elT = typeof(rand(rng, T, first(mar)))

    sz = size(mar)
    x = Array{elT,length(sz)}(undef, sz)
    rand!(rng, d, x)
end

# TODO:
# function Base.rand(rng::AbstractRNG, d::ProductMeasure)
#     return rand(rng, sampletype(d), d)
# end

# function Base.rand(T::Type, d::ProductMeasure)
#     return rand(Random.GLOBAL_RNG, T, d)
# end

# function Base.rand(d::ProductMeasure)
#     T = sampletype(d)
#     return rand(Random.GLOBAL_RNG, T, d)
# end

function sampletype(d::ProductMeasure{A}) where {T,N,A<:AbstractArray{T,N}}
    S = @inbounds sampletype(marginals(d)[1])
    Array{S,N}
end

function sampletype(d::ProductMeasure{<:Tuple})
    Tuple{sampletype.(marginals(d))...}
end

# function logdensity(μ::ProductMeasure{Aμ}, x::Ax) where {Aμ <: MappedArray, Ax <: AbstractArray}
#     μ.data
# end

function ConstructionBase.constructorof(::Type{P}) where {F,S,I,P<:ProductMeasure{F,S,I}}
    p -> productmeasure(d.f, p)
end

# function Accessors.set(d::ProductMeasure{N}, ::typeof(params), p) where {N}
#     setproperties(d, NamedTuple{N}(p...))
# end

# function Accessors.set(d::ProductMeasure{F,T}, ::typeof(params), p::Tuple) where {F, T<:Tuple}
#     set.(marginals(d), params, p)
# end

# function logdensity(μ::ProductMeasure, ν::ProductMeasure, x)
#     sum(zip(marginals(μ), marginals(ν), x)) do μ_ν_x
#         logdensity(μ_ν_x...)
#     end
# end

function kernelfactor(μ::ProductMeasure{F,S,<:Fill}) where {F,S}
    k = kernel(first(marginals(μ)))
    (p -> k.f(p)^size(μ), k.ops)
end

function kernelfactor(μ::ProductMeasure{F,S,A}) where {F,S,A<:AbstractArray}
    (p -> set.(marginals(μ), params, p), μ.pars)
end


text:
using pigpiod_if2_jll
include("./constants.jl")

# SERIAL
#  serial_open                  Opens a serial device
#  int serial_open(int pi, char *ser_tty, unsigned baud, unsigned ser_flags)
#  This function opens a serial device at a specified baud rate with specified flags. The device name must start with /dev/tty or /dev/serial.
function serial_open(p::Int32, ser_tty::String, baud::Int, ser_flags::Int)
    return ccall((:serial_open, libpigpiod_if2), Int32, (Int32, Cstring, UInt32, UInt32), p, ser_tty, baud, ser_flags)
end

#  serial_close                 Closes a serial device
#  int serial_close(int pi, unsigned handle)
#  Closes serial device associated with handle.
function serial_close(p::Int32, handle::Int32)
    return ccall((:serial_close, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

#  serial_read_byte             Reads a byte from a serial device
#  int serial_read_byte(int pi, unsigned handle)
#  This function reads a byte from the serial port associated with handle.
function serial_read_byte(p::Int32, handle::Int32)
    return ccall((:serial_read_byte, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

#  serial_write_byte            Writes a byte to a serial device
#  int serial_write_byte(int pi, unsigned handle, unsigned bVal)
#  This function writes bVal to the serial port associated with handle.
function serial_write_byte(p::Int32, handle::Int32, bVal::Int)
    return ccall((:serial_write_byte, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, handle, bVal)
end

#  serial_read	                Reads bytes from a serial device
#  int serial_read(int pi, unsigned handle, char *buf, unsigned count)
#  This function reads up to count bytes from the the serial port associated with handle and writes them to buf.
function serial_read(p::Int32, handle::Int32, buf::Vector{Int}, count::Int)
    return ccall((:serial_read, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32), p, handle, buf, count)
end

#  serial_write	                Writes bytes to a serial device
#  int serial_write(int pi, unsigned handle, char *buf, unsigned count)
#  This function writes count bytes from buf to the the serial port associated with handle.
function serial_write(p::Int32, handle::Int32, buf::Vector{Int}, count::Int)
    return ccall((:serial_write, libpigpiod_if2), Int32, (Int32, UInt32, Ref{UInt8}, UInt32), p, handle, buf, count)
end

#  serial_data_available	Returns number of bytes ready to be read
#  int serial_data_available(int pi, unsigned handle)
#  Returns the number of bytes available to be read from the device associated with handle.
function serial_data_available(p::Int32, handle::Int32)
    return ccall((:serial_data_available, libpigpiod_if2), Int32, (Int32, UInt32), p, handle)
end

# SERIAL BIT BANG (read only)
#  bb_serial_read_open	    Opens a GPIO for bit bang serial reads
#  int bb_serial_read_open(int pi, unsigned user_gpio, unsigned baud, unsigned data_bits)
#  This function opens a GPIO for bit bang reading of serial data.
function bb_serial_read_open(p::Int32, gpio::Int, baud::Int, bits::Int)
    return ccall((:bb_serial_read_open, libpigpiod_if2), Int32, (Int32, UInt32, UInt32, UInt32), p, gpio, baud, bits)
end

#  bb_serial_read_close	    Closes a GPIO for bit bang serial reads
#  int bb_serial_read_close(int pi, unsigned user_gpio)
#  This function closes a GPIO for bit bang reading of serial data.
function bb_serial_read_close(p::Int32, gpio::Int)
    return ccall((:bb_serial_read_close, libpigpiod_if2), Int32, (Int32, UInt32), p, gpio)
end

#  bb_serial_invert	    Invert serial logic (1 invert, 0 normal)
#  int bb_serial_invert(int pi, unsigned user_gpio, unsigned invert)
#  This function inverts serial logic for big bang serial reads.
function bb_serial_invert(p::Int32, gpio::Int, invert::Int)
    return ccall((:bb_serial_invert, libpigpiod_if2), Int32, (Int32, UInt32, UInt32), p, gpio, invert)
end

#  bb_serial_read	    Reads bit bang serial data from a GPIO
#  int bb_serial_read(int pi, unsigned user_gpio, void *buf, size_t bufSize)
#  This function copies up to bufSize bytes of data read from the bit bang serial cyclic buffer to the buffer starting at buf.
#  VERIFY THIS FUNCTION SINCE IT USES A VOID POINTER
function bb_serial_read(p::Int32, gpio::Int, buf, bufSize::Int)
    return ccall((:bb_serial_read, libpigpiod_if2), Int32, (Int32, UInt32, Ref{Any}, UInt32), p, gpio, buf, bufSize)
end


text:
using HTSEngineAPI
using Base.Test

using WAV

### Data ###

mei_htsvoice = joinpath(dirname(@__FILE__), "data", "mei_normal.htsvoice")
labelpath = joinpath(dirname(@__FILE__), "data","nitech_jp_atr503_m001_a01.lab")

### Tests ###

function test_hts_engine_setup()
    engine = HTS_Engine()
    load(engine, mei_htsvoice)
    @test_throws Exception load(engine, "doesnt_exists.htsvoice")
    engine = HTS_Engine(mei_htsvoice)
end

function test_htscall()
    engine = HTS_Engine(mei_htsvoice)
    synthesize_from_fn(engine, labelpath)

    # Get number of generated samples directly from GStream
    nsamples = signed(@htscall(:HTS_GStreamSet_get_total_nsamples, Csize_t,
                               (Ptr{HTS_GStreamSet},), &engine.gss))
    @test get_nsamples(engine) == nsamples
end

function test_hts_engine_properties()
    engine = HTS_Engine()

    for (name, val) in [
                        (:sampling_frequency, 48000),
                        (:fperiod, 240),
                        (:audio_buff_size, 10),
                        (:stop_flag, false),
                        (:volume, 2.0),
                        (:alpha, 0.42),
                        (:beta, 0.01)
                       ]
        setter = eval(Symbol(:set_, name))
        getter = eval(Symbol(:get_, name))
        println("test setter/getter property: $name")

        setter(engine, val)
        @test getter(engine) ≈ val
    end

    # no throws
    set_speed(engine, 0.7)
    set_phoneme_alignment_flag(engine, false)
    add_half_tone(engine, 0.5)

    engine = HTS_Engine(mei_htsvoice)

    for stream_index in 1:get_nstream(engine)
        for (name, val) in [
                            (:msd_threshold, 0.5),
                            (:gv_weight, 0.5)
                            ]
            setter = eval(Symbol(:set_, name))
            getter = eval(Symbol(:get_, name))
            println("test setter/getter property: $name")

            setter(engine, stream_index, val)
            @test getter(engine, stream_index) ≈ val
        end
    end

    for voice_index in 1:get_nvoices(engine)
        for (name, val) in [
                            (:duration_interpolation_weight, 0.5)
                            ]
            setter = eval(Symbol(:set_, name))
            getter = eval(Symbol(:get_, name))
            println("test setter/getter property: $name")

            setter(engine, voice_index, val)
            @test getter(engine, voice_index) ≈ val
        end
    end

    for voice_index in 1:get_nvoices(engine)
        for stream_index in 1:get_nstream(engine)
            for (name, val) in [
                                (:parameter_interpolation_weight, 0.5),
                                (:gv_interpolation_weight, 0.5)
                                ]
                setter = eval(Symbol(:set_, name))
                getter = eval(Symbol(:get_, name))
                println("test setter/getter property: $name")

                setter(engine, voice_index, stream_index, val)
                @test getter(engine, voice_index, stream_index) ≈ val
            end
        end
    end

    open(labelpath) do f
        synthesize_from_strings(engine, readlines(f))
    end

    @test get_total_state(engine) > 0
    for state_index in 1:get_nstate(engine)
        @test get_state_duration(engine, state_index) > 0
    end
    @test get_total_frame(engine) > 0
    @test get_nsamples(engine) > 0

    clear(engine)
end

function test_hts_engine_synthesis_funcs()
    engine = HTS_Engine(mei_htsvoice)

    # From filepath
    synthesize_from_fn(engine, labelpath)
    @test_throws Exception synthesize_from_fn(engine, "doesnt_exists.lab")

    synthesized_fn = get_generated_speech(engine)
    refresh(engine)

    # From file contents (lines)
    open(labelpath) do f
        synthesize_from_strings(engine, readlines(f))
    end
    synthesized_strings = get_generated_speech(engine)

    # should be same
    @test synthesized_fn == synthesized_strings
    @test length(synthesized_fn) == get_nsamples(engine)
end

function test_hts_engine_fullcontext_label()
    engine = HTS_Engine(mei_htsvoice)
    @test get_fullcontext_label_format(engine) == "HTS_TTS_JPN"
    @test get_fullcontext_label_version(engine) == "1.0"
end

function test_hts_engine_save()
    engine = HTS_Engine(mei_htsvoice)
    synthesize_from_fn(engine, labelpath)

    for savefunc in [
                     save_information,
                     save_label,
                     save_generated_speech,
                     save_riff
                     ]
        file = tempname()
        savefunc(engine, file)
        @test isfile(file)
        rm(file)
    end

    for stream_index = 1:get_nstream(engine)
        file = tempname()
        save_generated_parameter(engine, stream_index, file)
        @test isfile(file)
        rm(file)
    end

    clear(engine)
end

function test_hts_engine_synthesis_and_save_wav()
    # Load htsvoice
    engine = HTS_Engine(mei_htsvoice)
    @test get_sampling_frequency(engine) == 48000
    @test get_fperiod(engine) == 240

    # Synthesize waveform from label file
    synthesize_from_fn(engine, labelpath)

    # Save as wav file
    file = tempname()
    save_riff(engine, file)
    @test isfile(file)

    # Load wavfile
    x, fs = wavread(file, format="native")
    @test fs == get_sampling_frequency(engine)
    rm(file)
    @test length(x) == get_nsamples(engine)

    # get generated waveform
    synthesized = get_generated_speech(engine)
    synthesized = map(Int, map(trunc, synthesized))
    synthesized = map(x -> clamp(x, typemin(Int16), typemax(Int16)), synthesized)

    # check sample-by-sample correctess
    @test all(x .== synthesized)
end

### Run tests ###

info("test: hts_engine setup")
test_hts_engine_setup()

info("test: @htscall")
test_htscall()

info("test: hts_engine properties")
test_hts_engine_properties()

info("test: hts_engine synthesis functions")
test_hts_engine_synthesis_funcs()

info("test: hts_engine fullcontext label")
test_hts_engine_fullcontext_label()

info("test: hts_engine save")
test_hts_engine_save()

info("test: hts_engine syntheis and save to wav")
test_hts_engine_synthesis_and_save_wav()