text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# To propose an update to the Clp binaries please follow these steps:
# 1) Fork this ClpBuilder repository and create a new branch
# 2) Modify the build script (e.g., changing the portion marked by START/END-EASY-CHANGE-BLOCK) and test that it compiles
#    for all architectures by running `julia --color=yes build_tarballs.jl --verbose`. (Note that this build make take a
#    long time. Also, adding the --debug argument to the command will drop you into the build environment for debugging)
# 3) Create a github release in your fork of ClpBuilder with the generated tarballs and build_ClpBuilder.vX.X.X.jl
#    (https://github.com/tcnksm/ghr is usefull to automate this release creation)
# 4) Fork the Clp.jl repository (https://github.com/JuliaOpt/Clp.jl) and create a new branch
# 5) Update the portion marked by ## START-VERSION-UPDATE-BLOCK / ## END-VERSION-UPDATE-BLOCK in Clp.jl/deps/build.jl
#    on your branch of Clp.jl with the corresponding lines in the build_ClpBuilder.vX.X.X.jl file you generated in steps 2)/3)
# 6) Check that tests pass on your branch of Clp (e.g. `]test Clp`)
# 7) Create a PR for ClpBuilder.jl (No need to create a PR for Clp.jl)

# Ideally, any update would need changes only on the EASY-CHANGE-BLOCK below
# Note that if two sources have the same version then they need to have
# different extension so that the corresponding filenames are different.
# For instance, consider the following case where COINMumps and COINLapack
# both use versio 1.6.0:
#    COINMumps_version = v"1.6.0"
#    COINMumps_extension = "tar.gz"
#    COINLapack_version = v"1.6.0"
#    COINLapack_extension = "zip"
# In this case, we need one of the libraries to use the "tar.gz" extension and
# the other to use the "zip" extension. The reason for this is that the file names
# for github releases for COINMumps and COINLapack for the configuration above will
# be respectively `1.6.0.tar.gz` and `1.6.0.zip`. If we instead had used the same
# extensions, one of the files would have been overwritten.

##START-EASY-CHANGE-BLOCK
Clp_version = v"1.16.11"
Clp_extension = "tar.gz"
Clp_hash = "ac42c00ba95e1e034ae75ba0e3a5ff03b452191e0c9b2f5e2d5e65bf652fb0a1"
Osi_version = v"0.107.9"
Osi_extension = "tar.gz"
Osi_hash = "e2c8a0ee4a2a0abe7475d67f7f98230e8bfbbcb6e74487877e757c996bfd6d30"
CoinUtils_version = v"2.10.14"
CoinUtils_extension = "tar.gz"
CoinUtils_hash = "929b6eae0aaf62cf4467e506f24dfab1df7ab8d2e5a1ea71e9bab5480e872d84"
COINMumps_version = v"1.6.0"
COINMumps_extension = "tar.gz"
COINMumps_hash = "3f2bb7d13333e85a29cd2dadc78a38bbf469bc3920c4c0933a90b7d8b8dc798a"
COINMetis_version = v"1.3.5"
COINMetis_extension = "tar.gz"
COINMetis_hash = "98a6110d5d004a16ad42ee26cfac508477f44aa6fe296b90a6413fe0273ebe24"
COINLapack_version = v"1.5.6"
COINLapack_extension = "tar.gz"
COINLapack_hash = "c625dbb227e54e496430ffa708ddf23df5dbf173a0fcf570e1c249e13e411ba1"
COINBLAS_version = v"1.4.6"
COINBLAS_extension = "tar.gz"
COINBLAS_hash = "f9601efb98f04fdba220d49d5bda98d2a5a5e2ed7564df339bc7149b0c303f0c"
ASL_version = v"3.1.0"
ASL_extension = "tar.gz"
ASL_hash = "587c1a88f4c8f57bef95b58a8586956145417c8039f59b1758365ccc5a309ae9"

# List of symbols that will be externally visible in the Clp library
# The true symbol names are matched with a * pre- and post-appended (e.g. Clp matches *Clp*)
PRESERVE_SYMBOLS = ["Clp","maximumIterations"]
##END-EASY-CHANGE-BLOCK

name = "ClpBuilder"
# Collection of sources required to build ClpBuilder
sources = [
    "./bundled",
    "https://github.com/coin-or/Clp/archive/releases/$(Clp_version).$(Clp_extension)" =>
    Clp_hash,
    "https://github.com/coin-or/Osi/archive/releases/$(Osi_version).$(Osi_extension)" =>
    Osi_hash,
    "https://github.com/coin-or/CoinUtils/archive/releases/$(CoinUtils_version).$(CoinUtils_extension)" =>
    CoinUtils_hash,
    "https://github.com/coin-or-tools/ThirdParty-Mumps/archive/releases/$(COINMumps_version).$(COINMumps_extension)" =>
    COINMumps_hash,
    "https://github.com/coin-or-tools/ThirdParty-Metis/archive/releases/$(COINMetis_version).$(COINMetis_extension)" =>
    COINMetis_hash,
    "https://github.com/coin-or-tools/ThirdParty-Lapack/archive/releases/$(COINLapack_version).$(COINLapack_extension)" =>
    COINLapack_hash,
    "https://github.com/coin-or-tools/ThirdParty-Blas/archive/releases/$(COINBLAS_version).$(COINBLAS_extension)"  =>
    COINBLAS_hash,
    "https://github.com/ampl/mp/archive/$(ASL_version).tar.gz"  =>
    ASL_hash,
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
set -e
### Preliminaries
# Standard BB updates
update_configure_scripts
# Fix some paths
for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done
# Osi, Cgl, Clp, CoinUtils, mumps, metis, lapack, blas, asl
if [ $target = "x86_64-apple-darwin14" ]; then
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi
### Compile ASL
# Use staticfloat's cross-compile trick for ASL https://github.com/ampl/mp/issues/115
cd $WORKSPACE/srcdir/mp-*
rm -rf thirdparty/benchmark
patch -p1 < $WORKSPACE/srcdir/asl-extra/no_benchmark.patch
# Build ASL
mkdir build
cd build
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DRUN_HAVE_STD_REGEX=0 -DRUN_HAVE_STEADY_CLOCK=0 -DHAVE_ACCESS_DRIVER_EXITCODE=0 -DHAVE_EXCEL_DRIVER_EXITCODE=0 -DHAVE_ODBC_TEXT_DRIVER_EXITCODE=0    ../
else
   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=$prefix  -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain       -DRUN_HAVE_STD_REGEX=0       -DRUN_HAVE_STEADY_CLOCK=0       ../
fi
# Copy over pregenerated files after building arithchk, so as to fake out cmake,
# because cmake will delete our arith.h
## If this fails it is ok, we just want to prevend cmake deleting arith.h
set +e
make arith-h VERBOSE=1
set -e
mkdir -p src/asl
cp -v $WORKSPACE/srcdir/asl-extra/expr-info.cc ../src/expr-info.cc
cp -v $WORKSPACE/srcdir/asl-extra/arith.h.${target} src/asl/arith.h
# Build and install ASL
make -j${nproc} VERBOSE=1
make install VERBOSE=1
### Compile Lapack
cd $WORKSPACE/srcdir
cd ThirdParty-Blas-releases-*/
./get.Blas
mkdir build
cd build/
# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
make -j${nproc}
make install
### Compile Lapack
cd $WORKSPACE/srcdir
cd ThirdParty-Lapack-releases-*/
./get.Lapack
mkdir build
cd build/
# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config  --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-blas-lib="-L$prefix/lib -lcoinblas"
make -j${nproc}
make install
### Compile Metis
cd $WORKSPACE/srcdir
cd ThirdParty-Metis-releases-*/
./get.Metis
mkdir build
cd build/
# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
make -j${nproc}
make install
### Compile Mumps
cd $WORKSPACE/srcdir
cd ThirdParty-Mumps-releases-*/
./get.Mumps
patch -p1 < $WORKSPACE/srcdir/mumps-extra/quiet.diff
mkdir build
cd build/
# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
--with-blas-lib="-L$prefix/lib -lcoinblas"
make -j${nproc}
make install
### Compile CoinUtils
cd $WORKSPACE/srcdir
cd CoinUtils-releases-*/
mkdir build
cd build/
# Set env vars: CPPFLAGS
if [ $target = "aarch64-linux-gnu" ] || [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
fi
# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config  --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-blas-lib="-L${prefix}/lib -lcoinblas" \
--with-lapack-lib="-L${prefix}/lib -lcoinlapack"
make -j${nproc}
make install
# Clear env vars: CPPFLAGS
export CPPFLAGS=""
### Compile Osi
cd $WORKSPACE/srcdir
cd Osi-releases-*/
mkdir build
cd build/
# Set env vars: CXXFLAGS and CPPFLAGS
export CXXFLAGS="-std=c++11"
if [ $target = "aarch64-linux-gnu" ] || [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-std=c++11"
fi
# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
--with-lapack-lib="-L${prefix}/lib -lcoinlapack" \
--with-blas-lib="-L${prefix}/lib -lcoinblas"
make -j${nproc}
make install
# Clear env vars: CXXFLAGS and CPPFLAGS
export CXXFLAGS=""
export CPPFLAGS=""

## Compile Clp

cd $WORKSPACE/srcdir
cd Clp-releases-*/
mkdir build
cd build/

# Set env vars: CPPFLAGS
if [ $target = "aarch64-linux-gnu" ]; then
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H -D__arm__ -std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
elif [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H -D__arm__ -mfpu=neon -std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
else
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H"
fi

# Staticly link all dependencies and export only Clp symbols
# force only exporting symbols related to Clp
sed -i~ -e 's/LT_LDFLAGS="-no-undefined"/LT_LDFLAGS="-no-undefined -export-symbols-regex \\\\"BB_PRESERVE_SYMBOLS\\\\""/g' ../configure
sed -i~ -e 's/LT_LDFLAGS="-no-undefined"/LT_LDFLAGS="-no-undefined -export-symbols-regex \\\\"BB_PRESERVE_SYMBOLS\\\\""/g' ../Clp/configure

# configure, make and install
if [ $target = "x86_64-apple-darwin14" ]; then
  # Ignore the "# Don't fix this by using the ld -exported_symbols_list flag, it doesn't exist in older darwin lds"
  # seems to work for the current version and otherwise a long list of non-Clp symbols are exported
  sed -i~ -e "s|~nmedit -s \$output_objdir/\${libname}-symbols.expsym \${lib}| -exported_symbols_list \$output_objdir/\${libname}-symbols.expsym|g" ../configure
  # fix linking issue
  export OSICLPLIB_LIBS=" -lbz2 -lz ${prefix}/lib/libcoinlapack.a ${prefix}/lib/libCoinUtils.a ${prefix}/lib/libOsi.a ${prefix}/lib/libcoinblas.a"
  ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --enable-static  \
  --with-asl-lib="${prefix}/lib/libasl.a" --with-asl-incdir="$prefix/include/asl" \
  --with-blas-lib="${prefix}/lib/libcoinblas.a -lgfortran" \
  --with-lapack-lib="${prefix}/lib/libcoinlapack.a" \
  --with-metis-lib="${prefix}/lib/libcoinmetis.a" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
  --with-mumps-lib="${prefix}/lib/libcoinmetis.a ${prefix}/lib/libcoinblas.a ${prefix}/lib/libcoinmumps.a -L/opt/x86_64-apple-darwin14/x86_64-apple-darwin14/lib -lgfortran" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
  --with-coinutils-lib="${prefix}/lib/libcoinlapack.a ${prefix}/lib/libcoinblas.a ${prefix}/lib/libCoinUtils.a -lbz2 -lz" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="${prefix}/lib/libOsi.a" --with-osi-incdir="$prefix/include/coin" \
  lt_cv_deplibs_check_method=pass_all
elif [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  # fix linking issue
  export OSICLPLIB_LIBS="${prefix}/lib/libOsi.a   ${prefix}/lib/libCoinUtils.a ${prefix}/lib/libcoinlapack.a  ${prefix}/lib/libcoinblas.a  -lgfortran"
  ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --enable-static \
  --with-asl-lib="${prefix}/lib/libasl.a" --with-asl-incdir="$prefix/include/asl" \
  --with-lapack-lib="${prefix}/lib/libcoinlapack.a" \
  --with-mumps-lib="${prefix}/lib/libcoinmumps.a -lgfortran ${prefix}/lib/libcoinmetis.a" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
  --with-metis-lib="${prefix}/lib/libcoinmetis.a" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
  --with-coinutils-lib="${prefix}/lib/libCoinUtils.a ${prefix}/lib/libcoinblas.a ${prefix}/lib/libcoinlapack.a" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="${prefix}/lib/libOsi.a" --with-osi-incdir="$prefix/include/coin" \
  --with-blas-lib="${prefix}/lib/libcoinblas.a -lgfortran" \
  lt_cv_deplibs_check_method=pass_all
  # fix linking issue
  sed -i~ -e 's|libClpSolver_la_LIBADD = \$(CLPLIB_LIBS) libClp\.la|libClpSolver_la_LIBADD = $(CLPLIB_LIBS) libClp.la ${prefix}/lib/libcoinblas.a -lgfortran|g' Clp/src/Makefile
else
 ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --enable-static \
 --with-asl-lib="${prefix}/lib/libasl.a" --with-asl-incdir="$prefix/include/asl" \
 --with-lapack-lib="${prefix}/lib/libcoinlapack.a" \
 --with-mumps-lib="${prefix}/lib/libcoinmumps.a -lgfortran ${prefix}/lib/libcoinmetis.a" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
 --with-metis-lib="${prefix}/lib/libcoinmetis.a" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
 --with-coinutils-lib="${prefix}/lib/libCoinUtils.a ${prefix}/lib/libcoinblas.a ${prefix}/lib/libcoinlapack.a" --with-coinutils-incdir="$prefix/include/coin" \
 --with-osi-lib="${prefix}/lib/libOsi.a" --with-osi-incdir="$prefix/include/coin" \
 --with-blas-lib="${prefix}/lib/libcoinblas.a -lgfortran" \
 lt_cv_deplibs_check_method=pass_all LDFLAGS=-ldl;
fi

# COIN makefiles uses -retain-symbols-file in linux for symbol filtering, which does not seem to filter the dynamic table so
# we switch to -version-script
if [ $target = "x86_64-linux-gnu" ] || [ $target = "i686-linux-gnu" ]; then
  echo "{ global:" > $WORKSPACE/srcdir/names.ver
  echo "*BB_PRESERVE_SYMBOLS*;" | sed -e "s/|/*;*/g" >> $WORKSPACE/srcdir/names.ver
  echo "local: *; };" >> $WORKSPACE/srcdir/names.ver
  sed -i~ -e 's/archive_expsym_cmds=.*CC.*/archive_expsym_cmds="\\$CC -shared -nostdlib \\$predep_objects \\$libobjs \\$deplibs \\$postdep_objects \\$compiler_flags \\${wl}-soname \\$wl\\$soname \\${wl}-version-script \\${wl}\\$WORKSPACE\/srcdir\/names.ver -o \\$lib"/g' libtool
fi

make -j${nproc}

# Clean-up bin directory before installing
rm ${prefix}/bin/*

# Install
make install

# Clean-up lib directory
rm ${prefix}/lib/*.a
"""
BB_PRESERVE_SYMBOLS = join(PRESERVE_SYMBOLS, raw"""|""")
script = replace(script, "BB_PRESERVE_SYMBOLS" => BB_PRESERVE_SYMBOLS)

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)
# To fix gcc4 bug in Windows
platforms = setdiff(platforms, [Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)), Windows(:i686, compiler_abi=CompilerABI(:gcc4))])
push!(platforms, Windows(:i686,compiler_abi=CompilerABI(:gcc6)))
push!(platforms, Windows(:x86_64,compiler_abi=CompilerABI(:gcc6)))

# It seems Clp in ARM now needs C++11, so gcc4 is not working.
# A possible fix is to switch to gcc6 as for windows, but this will require updating all dependencies below too.
platforms = setdiff(platforms, [Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)), Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4))])

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libClp", :libClp),
    LibraryProduct(prefix, "libClpSolver", :libClpSolver)
]

# Dependencies that must be installed before this package can be built
dependencies = [
   # "https://github.com/JuliaOpt/OsiBuilder/releases/download/v0.107.9-1-static/build_OsiBuilder.v0.107.9.jl",
   # "https://github.com/JuliaOpt/CoinUtilsBuilder/releases/download/v2.10.14-1-static/build_CoinUtilsBuilder.v2.10.14.jl",
   # "https://github.com/JuliaOpt/COINMumpsBuilder/releases/download/v1.6.0-1-static/build_COINMumpsBuilder.v1.6.0.jl",
   # "https://github.com/JuliaOpt/COINMetisBuilder/releases/download/v1.3.5-1-static/build_COINMetisBuilder.v1.3.5.jl",
   # "https://github.com/JuliaOpt/COINLapackBuilder/releases/download/v1.5.6-1-static/build_COINLapackBuilder.v1.5.6.jl",
   # "https://github.com/JuliaOpt/COINBLASBuilder/releases/download/v1.4.6-1-static/build_COINBLASBuilder.v1.4.6.jl",
   # "https://github.com/JuliaOpt/ASLBuilder/releases/download/v3.1.0-1-static/build_ASLBuilder.v3.1.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, Clp_version, sources, script, platforms, products, dependencies)


text:
# Julia GARCH package
# Copyright 2013 Andrey Kolev
# Distributed under MIT license (see LICENSE.md)

type GarchFit
  data::Vector
  params::Vector
  llh::Float64
  status::Symbol
  converged::Bool
  sigma::Vector
  hessian::Array{Float64,2}
  cvar::Array{Float64,2}
  secoef::Vector
  tval::Vector
end

function Base.show(io::IO ,fit::GarchFit)
  pnorm(x) = 0.5*(1+erf(x/sqrt(2)))
  prt(x) = 2*(1-pnorm(abs(x)))
  @printf io "Fitted garch model \n"
  @printf io " * Coefficient(s): \tomega \t\talpha \t\tbeta\n"
  @printf io "   \t\t\t%f\t%f\t%f\n" fit.params[1] fit.params[2] fit.params[3]
  @printf io " * Log Likelihood: %f\n" fit.llh
  @printf io " * Converged: %s\n" fit.converged
  @printf io " * Solver status: %s\n\n" fit.status
  println(io," * Standardised Residuals Tests:")
  println(io,"   \t\t\t\tStatistic\tp-Value")
  jbstat,jbp = jbtest(fit.data./fit.sigma);
  @printf io "   Jarque-Bera Test\t\U1D6D8\u00B2\t%.6f\t%.6f\n\n" jbstat jbp
  println(io," * Error Analysis:")
  println(io,"   \t\tEstimate\t\Std.Error\tt value \tPr(>|t|)")
  @printf io "   omega\t%f\t%f\t%f\t%f\n" fit.params[1] fit.secoef[1] fit.tval[1] prt(fit.tval[1])
  @printf io "   alpha\t%f\t%f\t%f\t%f\n" fit.params[2] fit.secoef[2] fit.tval[2] prt(fit.tval[2])
  @printf io "   beta \t%f\t%f\t%f\t%f\n"  fit.params[3] fit.secoef[3] fit.tval[3] prt(fit.tval[3])
end

function cdHessian(par,LLH)
  eps = 1e-4 * par
  n = length(par)
  H = zeros(n,n)
  for(i = 1:n)
    for(j = 1:n)
      x1 = copy(par)
      x1[i] += eps[i]
      x1[j] += eps[j]
      x2 = copy(par)
      x2[i] += eps[i]
      x2[j] -= eps[j]
      x3 = copy(par)
      x3[i] -= eps[i]
      x3[j] += eps[j]
      x4 = copy(par)
      x4[i] -= eps[i]
      x4[j] -= eps[j]
      H[i,j] = (LLH(x1)-LLH(x2)-LLH(x3)+LLH(x4)) / (4.*eps[i]*eps[j])
    end
  end
  H
end

function garchLLH(rets::Vector,x::Vector)
  rets2   = rets.^2;
  T = length(rets);
  ht = zeros(T);
  omega,alpha,beta = x;
  ht[1] = sum(rets2)/T;
  for i=2:T
    ht[i] = omega + alpha*rets2[i-1] + beta * ht[i-1];
  end
  -0.5*(T-1)*log(2*pi)-0.5*sum( log(ht) + (rets./sqrt(ht)).^2 );
end

function predict(fit::GarchFit)
 omega, alpha, beta = fit.params;
 rets = fit.data
 rets2   = rets.^2;
 T = length(rets);
 ht    = zeros(T);
 ht[1] = sum(rets2)/T;
 for i=2:T
    ht[i] = omega + alpha*rets2[i-1] + beta * ht[i-1];
 end
 sqrt(omega + alpha*rets2[end] + beta*ht[end]);
end

function garchFit(data::Vector)
  rets = data
  rets2   = rets.^2;
  T = length(rets);
  ht = zeros(T);
  function garchLike(x::Vector, grad::Vector)
    omega,alpha,beta = x;
    ht[1] = sum(rets2)/T;
    for i=2:T
      ht[i] = omega + alpha*rets2[i-1] + beta * ht[i-1];
    end
    sum( log(ht) + (rets./sqrt(ht)).^2 );
  end
  opt = Opt(:LN_SBPLX,3)
  lower_bounds!(opt,[1e-10, 0.0, 0.0])
  upper_bounds!(opt,[1; 0.3; 0.99])
  min_objective!(opt, garchLike)
  (minf,minx,ret) = NLopt.optimize(opt, [1e-5, 0.09, 0.89])
  converged = minx[1]>0 && all(minx[2:3].>=0) && sum(minx[2:3])<1.0
  H = cdHessian(minx,x->garchLLH(rets,x))
  cvar = -inv(H)
  secoef = sqrt(diag(cvar))
  tval = minx./secoef
  out = GarchFit(data, minx, -0.5*(T-1)*log(2*pi)-0.5*minf, ret, converged, sqrt(ht),H,cvar,secoef,tval)
end

# function garchPkgTest()
#   println("Running GARCH package test...")
#   try
#     include(Pkg.dir("GARCH", "test","GARCHtest.jl"))
#     println("All tests passed!")
#   catch err
#     throw(err)
#   end
# end

ext: jl


text:
# Julia wrapper for header: /opt/ros/crystal/include/rmw/allocators.h
# Automatically generated using Clang.jl wrap_c

function rmw_allocate(size)
    ccall((:rmw_allocate, librmw), Ptr{Cvoid}, (Csize_t,), size)
end

function rmw_free(pointer)
    ccall((:rmw_free, librmw), Cvoid, (Ptr{Cvoid},), pointer)
end

function rmw_node_allocate()
    ccall((:rmw_node_allocate, librmw), Ptr{rmw_node_t}, ())
end

function rmw_node_free(node)
    ccall((:rmw_node_free, librmw), Cvoid, (Ptr{rmw_node_t},), node)
end

function rmw_publisher_allocate()
    ccall((:rmw_publisher_allocate, librmw), Ptr{rmw_publisher_t}, ())
end

function rmw_publisher_free(publisher)
    ccall((:rmw_publisher_free, librmw), Cvoid, (Ptr{rmw_publisher_t},), publisher)
end

function rmw_subscription_allocate()
    ccall((:rmw_subscription_allocate, librmw), Ptr{rmw_subscription_t}, ())
end

function rmw_subscription_free(subscription)
    ccall((:rmw_subscription_free, librmw), Cvoid, (Ptr{rmw_subscription_t},), subscription)
end

function rmw_guard_condition_allocate()
    ccall((:rmw_guard_condition_allocate, librmw), Ptr{rmw_guard_condition_t}, ())
end

function rmw_guard_condition_free(guard_condition)
    ccall((:rmw_guard_condition_free, librmw), Cvoid, (Ptr{rmw_guard_condition_t},), guard_condition)
end

function rmw_client_allocate()
    ccall((:rmw_client_allocate, librmw), Ptr{rmw_client_t}, ())
end

function rmw_client_free(client)
    ccall((:rmw_client_free, librmw), Cvoid, (Ptr{rmw_client_t},), client)
end

function rmw_service_allocate()
    ccall((:rmw_service_allocate, librmw), Ptr{rmw_service_t}, ())
end

function rmw_service_free(service)
    ccall((:rmw_service_free, librmw), Cvoid, (Ptr{rmw_service_t},), service)
end

function rmw_wait_set_allocate()
    ccall((:rmw_wait_set_allocate, librmw), Ptr{rmw_wait_set_t}, ())
end

function rmw_wait_set_free(wait_set)
    ccall((:rmw_wait_set_free, librmw), Cvoid, (Ptr{rmw_wait_set_t},), wait_set)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/convert_rcutils_ret_to_rmw_ret.h
# Automatically generated using Clang.jl wrap_c

function rmw_convert_rcutils_ret_to_rmw_ret(rcutils_ret)
    ccall((:rmw_convert_rcutils_ret_to_rmw_ret, librmw), rmw_ret_t, (rcutils_ret_t,), rcutils_ret)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/error_handling.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/get_node_info_and_types.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_subscriber_names_and_types_by_node(node, allocator, node_name, node_namespace, demangle, topics_names_and_types)
    ccall((:rmw_get_subscriber_names_and_types_by_node, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rcutils_allocator_t}, Cstring, Cstring, Bool, Ptr{rmw_names_and_types_t}), node, allocator, node_name, node_namespace, demangle, topics_names_and_types)
end

function rmw_get_publisher_names_and_types_by_node(node, allocator, node_name, node_namespace, demangle, topic_names_and_types)
    ccall((:rmw_get_publisher_names_and_types_by_node, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rcutils_allocator_t}, Cstring, Cstring, Bool, Ptr{rmw_names_and_types_t}), node, allocator, node_name, node_namespace, demangle, topic_names_and_types)
end

function rmw_get_service_names_and_types_by_node(node, allocator, node_name, node_namespace, service_names_and_types)
    ccall((:rmw_get_service_names_and_types_by_node, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rcutils_allocator_t}, Cstring, Cstring, Ptr{rmw_names_and_types_t}), node, allocator, node_name, node_namespace, service_names_and_types)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/get_service_names_and_types.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_service_names_and_types(node, allocator, service_names_and_types)
    ccall((:rmw_get_service_names_and_types, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rcutils_allocator_t}, Ptr{rmw_names_and_types_t}), node, allocator, service_names_and_types)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/get_topic_names_and_types.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_topic_names_and_types(node, allocator, no_demangle, topic_names_and_types)
    ccall((:rmw_get_topic_names_and_types, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rcutils_allocator_t}, Bool, Ptr{rmw_names_and_types_t}), node, allocator, no_demangle, topic_names_and_types)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/init.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_zero_initialized_context()
    ccall((:rmw_get_zero_initialized_context, librmw), rmw_context_t, ())
end

function rmw_init(options, context)
    ccall((:rmw_init, librmw), rmw_ret_t, (Ptr{rmw_init_options_t}, Ptr{rmw_context_t}), options, context)
end

function rmw_shutdown(context)
    ccall((:rmw_shutdown, librmw), rmw_ret_t, (Ptr{rmw_context_t},), context)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/init_options.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_zero_initialized_init_options()
    ccall((:rmw_get_zero_initialized_init_options, librmw), rmw_init_options_t, ())
end

function rmw_init_options_init(init_options, allocator)
    ccall((:rmw_init_options_init, librmw), rmw_ret_t, (Ptr{rmw_init_options_t}, rcutils_allocator_t), init_options, allocator)
end

function rmw_init_options_copy(src, dst)
    ccall((:rmw_init_options_copy, librmw), rmw_ret_t, (Ptr{rmw_init_options_t}, Ptr{rmw_init_options_t}), src, dst)
end

function rmw_init_options_fini(init_options)
    ccall((:rmw_init_options_fini, librmw), rmw_ret_t, (Ptr{rmw_init_options_t},), init_options)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/macros.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/names_and_types.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_zero_initialized_names_and_types()
    ccall((:rmw_get_zero_initialized_names_and_types, librmw), rmw_names_and_types_t, ())
end

function rmw_names_and_types_check_zero(names_and_types)
    ccall((:rmw_names_and_types_check_zero, librmw), rmw_ret_t, (Ptr{rmw_names_and_types_t},), names_and_types)
end

function rmw_names_and_types_init(names_and_types, size, allocator)
    ccall((:rmw_names_and_types_init, librmw), rmw_ret_t, (Ptr{rmw_names_and_types_t}, Csize_t, Ptr{rcutils_allocator_t}), names_and_types, size, allocator)
end

function rmw_names_and_types_fini(names_and_types)
    ccall((:rmw_names_and_types_fini, librmw), rmw_ret_t, (Ptr{rmw_names_and_types_t},), names_and_types)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/node_security_options.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_zero_initialized_node_security_options()
    ccall((:rmw_get_zero_initialized_node_security_options, librmw), rmw_node_security_options_t, ())
end

function rmw_get_default_node_security_options()
    ccall((:rmw_get_default_node_security_options, librmw), rmw_node_security_options_t, ())
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/qos_profiles.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/ret_types.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/rmw.h
# Automatically generated using Clang.jl wrap_c

function rmw_get_implementation_identifier()
    ccall((:rmw_get_implementation_identifier, librmw), Cstring, ())
end

function rmw_get_serialization_format()
    ccall((:rmw_get_serialization_format, librmw), Cstring, ())
end

function rmw_create_node(context, name, namespace_, domain_id, security_options)
    ccall((:rmw_create_node, librmw), Ptr{rmw_node_t}, (Ptr{rmw_context_t}, Cstring, Cstring, Csize_t, Ptr{rmw_node_security_options_t}), context, name, namespace_, domain_id, security_options)
end

function rmw_destroy_node(node)
    ccall((:rmw_destroy_node, librmw), rmw_ret_t, (Ptr{rmw_node_t},), node)
end

function rmw_node_get_graph_guard_condition(node)
    ccall((:rmw_node_get_graph_guard_condition, librmw), Ptr{rmw_guard_condition_t}, (Ptr{rmw_node_t},), node)
end

function rmw_create_publisher(node, type_support, topic_name, qos_policies)
    ccall((:rmw_create_publisher, librmw), Ptr{rmw_publisher_t}, (Ptr{rmw_node_t}, Ptr{rosidl_message_type_support_t}, Cstring, Ptr{rmw_qos_profile_t}), node, type_support, topic_name, qos_policies)
end

function rmw_destroy_publisher(node, publisher)
    ccall((:rmw_destroy_publisher, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rmw_publisher_t}), node, publisher)
end

function rmw_publish(publisher, ros_message)
    ccall((:rmw_publish, librmw), rmw_ret_t, (Ptr{rmw_publisher_t}, Ptr{Cvoid}), publisher, ros_message)
end

function rmw_publisher_count_matched_subscriptions(publisher, subscription_count)
    ccall((:rmw_publisher_count_matched_subscriptions, librmw), rmw_ret_t, (Ptr{rmw_publisher_t}, Ptr{Csize_t}), publisher, subscription_count)
end

function rmw_publish_serialized_message(publisher, serialized_message)
    ccall((:rmw_publish_serialized_message, librmw), rmw_ret_t, (Ptr{rmw_publisher_t}, Ptr{rmw_serialized_message_t}), publisher, serialized_message)
end

function rmw_serialize(ros_message, type_support, serialized_message)
    ccall((:rmw_serialize, librmw), rmw_ret_t, (Ptr{Cvoid}, Ptr{rosidl_message_type_support_t}, Ptr{rmw_serialized_message_t}), ros_message, type_support, serialized_message)
end

function rmw_deserialize(serialized_message, type_support, ros_message)
    ccall((:rmw_deserialize, librmw), rmw_ret_t, (Ptr{rmw_serialized_message_t}, Ptr{rosidl_message_type_support_t}, Ptr{Cvoid}), serialized_message, type_support, ros_message)
end

function rmw_create_subscription(node, type_support, topic_name, qos_policies, ignore_local_publications)
    ccall((:rmw_create_subscription, librmw), Ptr{rmw_subscription_t}, (Ptr{rmw_node_t}, Ptr{rosidl_message_type_support_t}, Cstring, Ptr{rmw_qos_profile_t}, Bool), node, type_support, topic_name, qos_policies, ignore_local_publications)
end

function rmw_destroy_subscription(node, subscription)
    ccall((:rmw_destroy_subscription, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rmw_subscription_t}), node, subscription)
end

function rmw_subscription_count_matched_publishers(subscription, publisher_count)
    ccall((:rmw_subscription_count_matched_publishers, librmw), rmw_ret_t, (Ptr{rmw_subscription_t}, Ptr{Csize_t}), subscription, publisher_count)
end

function rmw_take(subscription, ros_message, taken)
    ccall((:rmw_take, librmw), rmw_ret_t, (Ptr{rmw_subscription_t}, Ptr{Cvoid}, Ptr{Bool}), subscription, ros_message, taken)
end

function rmw_take_with_info(subscription, ros_message, taken, message_info)
    ccall((:rmw_take_with_info, librmw), rmw_ret_t, (Ptr{rmw_subscription_t}, Ptr{Cvoid}, Ptr{Bool}, Ptr{rmw_message_info_t}), subscription, ros_message, taken, message_info)
end

function rmw_take_serialized_message(subscription, serialized_message, taken)
    ccall((:rmw_take_serialized_message, librmw), rmw_ret_t, (Ptr{rmw_subscription_t}, Ptr{rmw_serialized_message_t}, Ptr{Bool}), subscription, serialized_message, taken)
end

function rmw_take_serialized_message_with_info(subscription, serialized_message, taken, message_info)
    ccall((:rmw_take_serialized_message_with_info, librmw), rmw_ret_t, (Ptr{rmw_subscription_t}, Ptr{rmw_serialized_message_t}, Ptr{Bool}, Ptr{rmw_message_info_t}), subscription, serialized_message, taken, message_info)
end

function rmw_create_client(node, type_support, service_name, qos_policies)
    ccall((:rmw_create_client, librmw), Ptr{rmw_client_t}, (Ptr{rmw_node_t}, Ptr{rosidl_service_type_support_t}, Cstring, Ptr{rmw_qos_profile_t}), node, type_support, service_name, qos_policies)
end

function rmw_destroy_client(node, client)
    ccall((:rmw_destroy_client, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rmw_client_t}), node, client)
end

function rmw_send_request(client, ros_request, sequence_id)
    ccall((:rmw_send_request, librmw), rmw_ret_t, (Ptr{rmw_client_t}, Ptr{Cvoid}, Ptr{Int64}), client, ros_request, sequence_id)
end

function rmw_take_response(client, request_header, ros_response, taken)
    ccall((:rmw_take_response, librmw), rmw_ret_t, (Ptr{rmw_client_t}, Ptr{rmw_request_id_t}, Ptr{Cvoid}, Ptr{Bool}), client, request_header, ros_response, taken)
end

function rmw_create_service(node, type_support, service_name, qos_policies)
    ccall((:rmw_create_service, librmw), Ptr{rmw_service_t}, (Ptr{rmw_node_t}, Ptr{rosidl_service_type_support_t}, Cstring, Ptr{rmw_qos_profile_t}), node, type_support, service_name, qos_policies)
end

function rmw_destroy_service(node, service)
    ccall((:rmw_destroy_service, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rmw_service_t}), node, service)
end

function rmw_take_request(service, request_header, ros_request, taken)
    ccall((:rmw_take_request, librmw), rmw_ret_t, (Ptr{rmw_service_t}, Ptr{rmw_request_id_t}, Ptr{Cvoid}, Ptr{Bool}), service, request_header, ros_request, taken)
end

function rmw_send_response(service, request_header, ros_response)
    ccall((:rmw_send_response, librmw), rmw_ret_t, (Ptr{rmw_service_t}, Ptr{rmw_request_id_t}, Ptr{Cvoid}), service, request_header, ros_response)
end

function rmw_create_guard_condition(context)
    ccall((:rmw_create_guard_condition, librmw), Ptr{rmw_guard_condition_t}, (Ptr{rmw_context_t},), context)
end

function rmw_destroy_guard_condition(guard_condition)
    ccall((:rmw_destroy_guard_condition, librmw), rmw_ret_t, (Ptr{rmw_guard_condition_t},), guard_condition)
end

function rmw_trigger_guard_condition(guard_condition)
    ccall((:rmw_trigger_guard_condition, librmw), rmw_ret_t, (Ptr{rmw_guard_condition_t},), guard_condition)
end

function rmw_create_wait_set(max_conditions)
    ccall((:rmw_create_wait_set, librmw), Ptr{rmw_wait_set_t}, (Csize_t,), max_conditions)
end

function rmw_destroy_wait_set(wait_set)
    ccall((:rmw_destroy_wait_set, librmw), rmw_ret_t, (Ptr{rmw_wait_set_t},), wait_set)
end

function rmw_wait(subscriptions, guard_conditions, services, clients, wait_set, wait_timeout)
    ccall((:rmw_wait, librmw), rmw_ret_t, (Ptr{rmw_subscriptions_t}, Ptr{rmw_guard_conditions_t}, Ptr{rmw_services_t}, Ptr{rmw_clients_t}, Ptr{rmw_wait_set_t}, Ptr{rmw_time_t}), subscriptions, guard_conditions, services, clients, wait_set, wait_timeout)
end

function rmw_get_node_names(node, node_names, node_namespaces)
    ccall((:rmw_get_node_names, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rcutils_string_array_t}, Ptr{rcutils_string_array_t}), node, node_names, node_namespaces)
end

function rmw_count_publishers(node, topic_name, count)
    ccall((:rmw_count_publishers, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Cstring, Ptr{Csize_t}), node, topic_name, count)
end

function rmw_count_subscribers(node, topic_name, count)
    ccall((:rmw_count_subscribers, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Cstring, Ptr{Csize_t}), node, topic_name, count)
end

function rmw_get_gid_for_publisher(publisher, gid)
    ccall((:rmw_get_gid_for_publisher, librmw), rmw_ret_t, (Ptr{rmw_publisher_t}, Ptr{rmw_gid_t}), publisher, gid)
end

function rmw_compare_gids_equal(gid1, gid2, result)
    ccall((:rmw_compare_gids_equal, librmw), rmw_ret_t, (Ptr{rmw_gid_t}, Ptr{rmw_gid_t}, Ptr{Bool}), gid1, gid2, result)
end

function rmw_service_server_is_available(node, client, is_available)
    ccall((:rmw_service_server_is_available, librmw), rmw_ret_t, (Ptr{rmw_node_t}, Ptr{rmw_client_t}, Ptr{Bool}), node, client, is_available)
end

function rmw_set_log_severity(severity)
    ccall((:rmw_set_log_severity, librmw), rmw_ret_t, (rmw_log_severity_t,), severity)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/sanity_checks.h
# Automatically generated using Clang.jl wrap_c

function rmw_check_zero_rmw_string_array(array)
    ccall((:rmw_check_zero_rmw_string_array, librmw), rmw_ret_t, (Ptr{rcutils_string_array_t},), array)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/serialized_message.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/types.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/validate_full_topic_name.h
# Automatically generated using Clang.jl wrap_c

function rmw_validate_full_topic_name(topic_name, validation_result, invalid_index)
    ccall((:rmw_validate_full_topic_name, librmw), rmw_ret_t, (Cstring, Ptr{Cint}, Ptr{Csize_t}), topic_name, validation_result, invalid_index)
end

function rmw_validate_full_topic_name_with_size(topic_name, topic_name_length, validation_result, invalid_index)
    ccall((:rmw_validate_full_topic_name_with_size, librmw), rmw_ret_t, (Cstring, Csize_t, Ptr{Cint}, Ptr{Csize_t}), topic_name, topic_name_length, validation_result, invalid_index)
end

function rmw_full_topic_name_validation_result_string(validation_result)
    ccall((:rmw_full_topic_name_validation_result_string, librmw), Cstring, (Cint,), validation_result)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/validate_namespace.h
# Automatically generated using Clang.jl wrap_c

function rmw_validate_namespace(namespace_, validation_result, invalid_index)
    ccall((:rmw_validate_namespace, librmw), rmw_ret_t, (Cstring, Ptr{Cint}, Ptr{Csize_t}), namespace_, validation_result, invalid_index)
end

function rmw_validate_namespace_with_size(namespace_, namespace_length, validation_result, invalid_index)
    ccall((:rmw_validate_namespace_with_size, librmw), rmw_ret_t, (Cstring, Csize_t, Ptr{Cint}, Ptr{Csize_t}), namespace_, namespace_length, validation_result, invalid_index)
end

function rmw_namespace_validation_result_string(validation_result)
    ccall((:rmw_namespace_validation_result_string, librmw), Cstring, (Cint,), validation_result)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/validate_node_name.h
# Automatically generated using Clang.jl wrap_c

function rmw_validate_node_name(node_name, validation_result, invalid_index)
    ccall((:rmw_validate_node_name, librmw), rmw_ret_t, (Cstring, Ptr{Cint}, Ptr{Csize_t}), node_name, validation_result, invalid_index)
end

function rmw_validate_node_name_with_size(node_name, node_name_length, validation_result, invalid_index)
    ccall((:rmw_validate_node_name_with_size, librmw), rmw_ret_t, (Cstring, Csize_t, Ptr{Cint}, Ptr{Csize_t}), node_name, node_name_length, validation_result, invalid_index)
end

function rmw_node_name_validation_result_string(validation_result)
    ccall((:rmw_node_name_validation_result_string, librmw), Cstring, (Cint,), validation_result)
end
# Julia wrapper for header: /opt/ros/crystal/include/rmw/visibility_control.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/impl/config.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/impl/cpp/demangle.hpp
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/impl/cpp/key_value.hpp
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rmw/impl/cpp/macros.hpp
# Automatically generated using Clang.jl wrap_c


text:
"""
not finished experiments on using WMMA tensor cores
"""
# """
# calculates Mahalanobis distance between two segmentations
# 1)so easist is to think about those segmentations as sparse matricies where we store in 3 columns x,y,z coordinates of all entries where we have true
# 	-we need mean from those entries only we ignore zeros
# 	- hence first pass can be iterating over both arrays - do the reduction understood in two ways
# 		- we need 3 sums and count of the x , y znad z components - in order to achieve this we can apply classic reduction techniques
# 		all of blocks will accumulate sum of elements  and count, it will also add point coordinates (x,y,z) to shared memory queue - if shared memory queue will
# 		get filled to its capacity we will update the main global sparse matrix representation using atomics to update the counter
# 		- summation will be done as ussual using first accumulation in private variables than summation across
# 2) so we should have 3 sums (of x,y and z coordinates) and their count and a 3 column table with all x,y,z coordinates of the non empty entries
# 	- we will have this data for both  arrays - representing gold standard and other segmentation
# 	now we need to calculate covariance matriscies of each point using global means and accumulate (add) all of those matricies
# 	We will use WMMA to speed things up - and in order to utilize WMMA constraints we will write calculation of covariance matrix in
# 	matrix multiplication notation look into https://docs.google.com/document/d/1G5FqjRj7WrDs4LWtBH667_orA8HzGa_JtTQ5Yb8vu8E/edit
# 	- in order to get those matricies we just need to load values prepared in step 1 - probably best way would be to keep means in shared memory

# https://github.com/JuliaGPU/GemmKernels.jl/blob/40c1dacb2ff2d24d7d392fc784c389e5e7fd8307/src/kernel.jl

# example of WMMA https://github.com/JuliaGPU/CUDA.jl/blob/800b7b89c4c19b9b98a7d150a813a0e3d0e18be5/examples/wmma/high-level.jl
# look into https://math.stackexchange.com/questions/4240707/mahalanobis-distance-between-two-3-dimensional-boolean-arrays

# taken from https://github.com/JuliaGPU/GemmKernels.jl/blob/40c1dacb2ff2d24d7d392fc784c389e5e7fd8307/src/kernel.jl
# - probably how to load element by element to fragment

# """
# module Mahalanobis
# export WMMAkernel
# using CUDA,..IterationUtils, ..ReductionUtils, ..MemoryUtils
# """

# First we upload means meanx, mean y and mean z to shared memory
# We defined shared memory for fragmentA , fragmentB and C - A and B we will define on every iteration
# 	fragment C will accumulate all
# Before we start we need to set all shared memory to 0

# We want to iterate over our arrays - first the gold than other one
# we will check is in given spot of the array the predicate of equality is met - if yes we will write the output into shared memory fragments and update local counter
# counters and fragments will be warp owned and managed
# every time the  true will be found the counter will be updated and on this basis value (old value from counter) will be put into the fragA and fragB
# x loop
# 	- every time we will finish  the iteration of x loop we will check the warp counter and on the basis of it we will either continue  or not
# 	in the simplest version if warp private counter will be greater  than (224-33) we will update fragmentA and B  execute WMMA
# 	- prepareing the frag A and B will be mainly related to fill with as many (xses - meanx) as we have first  14 rows for frag A and first 14 columns in B
# 		last columns/rows need to be filled by current (y - mean),(z-mean)
# 	-when we will fill the matrix before pushing it into WMMA we need to add some modificatiions
# 		- in the main diagonal we have variances of x entries (apart from last 2 in basic algorithm) - here we do not need to achange anything
# 		- in last two entries of main diagonal we will have 16 y variances and 16 z variances - this number needs to adjusted to be equal to local counter
# 			so we divide by 16 and multiply by local counter value
# 		- in entries [15,16] and [16,15] we will have yz covariance
# 		- covariances of xses with x and y will be in last two columns apart from right bottom corner
# 	-after each WMMA we need to set all entries of fragments and counter to 0
# y loop
#    - we check wheather result counter evaluated to more than 0 if yes we fire procedure described above if not we get to next y we are intrested in
# z loop
# 	we check weather our results are not empty by checking the value of [1,1] of frag A if it is grater than 0
#     - in case we are intreseted in slicewise results we save the covariance data into  global memory slice wise covariance data
# 	- in this case we reset
# after we got through all iterations we add data about covariances to global memory (atomically)

# goldArr,segmArr  - arrays we analyze
# numberToLooFor - number we are intrested in in the array
# arrDims - dimensions of main array
# loopXdim, loopYdim,  loopZdim - indicates how many times we need to loop over dimensions to cover all
# """
# function WMMAkernel(goldArr,segmArr
#     ,numberToLooFor
# 	,arrDims, loopXdim, loopYdim,  loopZdim
# 	 )

# 	warpNum = 1
# 	#store matricies we will multiply
# 	warpFragA = @cuStaticSharedMem(Float16, (16,16))
# 	warpFragB = @cuStaticSharedMem(Float16, (16,16))
# 	#matrix with accumulated results
# 	warpFragC = @cuStaticSharedMem(Float16, (16,16))
# 	#keeps track how many x true entries we have
# 	warpXCounter = @cuStaticSharedMem(UInt16, (1))
# 	#resetting all
# 	@ifXY warpNum 1 warpXCounter=UInt16(0)
# 	#resetting using multiple warps
# 	@ifY warpNum begin clearFrag(warpFragA) ;  clearFrag(warpFragB)  ;clearFrag(warpFragC)  end

# 	#iterating over data looking for places that meet our predicate
# 	@iter3dAdditionalxyzActs(arrDims,loopXdim,loopYdim,loopZdim
# 	#expression in very center
# 	,begin    end
# 	#additionalActionAfterX
# 	,begin    end
# 	#additionalActionAfterY
# 	,begin    end
# 	#additionalActionAfterZ
# 	,begin    end
# 	)

# 	 conf = WMMA.Config{16, 16, 16, Float16}
# 	 a_frag = WMMA.load_a(pointer(fragA), 16, WMMA.ColMajor, conf)
# 	 b_frag = WMMA.load_b(pointer(fragB), 16, WMMA.ColMajor, conf)
# 	 c_frag = WMMA.load_c(pointer(ones), 16, WMMA.ColMajor, conf)
# 	d_frag = WMMA.mma(a_frag, b_frag, c_frag, conf)
# 	WMMA.store_d(pointer(d_out), d_frag, 16, WMMA.ColMajor, conf)

# 		 return
# 	end

# end#module

# """
# we have a 16 by 16 fragment and we need to fill it with 0s
# """
# function clearFrag(frag)
# 	@unroll for i in 1:8
# 		#indexing so for example on first iteration so i=1 when threadIdx<=16 we will cover in that wayy all first row ; and if it is bigger than 16 we will covering second row
# 		frag[((i-1) & (2^4 - 1))+1,(((threadIdxX()-1)>>4 )+1)*i]= Float16(0)
# 	end#for
# end#clearFrag

# # function WMMAkernel(dataShmem,ones,d_out,fragA,fragB)
# # 	#  fragA = @cuStaticSharedMem(Float16, (4,4))
# # 	#  fragB = @cuStaticSharedMem(Float16, (4,4))
# # 	# fragC = @cuStaticSharedMem(Float16, (4,4))

# # 	 if(threadIdxX()<17)
# # 		a =  ((threadIdxX()-1) & (3))+1
# # 		b = ~((threadIdxX()+3)>>2 & 1)+3
# # 		c = ((threadIdxX()-1)>>2 )+1
# # 		fragA[a,c] = dataShmem[b,a]*( ((threadIdxX()>4 && threadIdxX()<13)*-2)+1 )
# # 	else
# # 		a = ((threadIdxX()-1) & (3))+1
# # 		c = (((threadIdxX()-16)-1)>>2 )+1
# # 		d = ((threadIdxX()-17) >>3)+1
# # 		fragB[ c,a] = dataShmem[ d,a ]
# # 	end
# # 	 conf = WMMA.Config{16, 16, 16, Float16}
# # 	 a_frag = WMMA.load_a(pointer(fragA), 16, WMMA.ColMajor, conf)
# # 	 b_frag = WMMA.load_b(pointer(fragB), 16, WMMA.ColMajor, conf)
# # 	 c_frag = WMMA.load_c(pointer(ones), 16, WMMA.ColMajor, conf)
# # 	d_frag = WMMA.mma(a_frag, b_frag, c_frag, conf)
# # 	WMMA.store_d(pointer(d_out), d_frag, 16, WMMA.ColMajor, conf)

# # 		 return
# # 	end


text:
# DOCTODO  add docstring for TechRenewable
struct TechRenewable <: TechnicalParams
    installedcapacity::Float64 # [MW]
    reactivepowerlimits::Union{NamedTuple{(:min, :max),Tuple{Float64,Float64}},Nothing} # [MVar]
    powerfactor::Union{Float64,Nothing} # [-1. -1]
end

# DOCTODO document this constructor for TechRenewable
TechRenewable(; InstalledCapacity = 0,
                reactivepowerlimits = nothing,
                powerfactor = nothing
              ) = TechRenewable(InstalledCapacity, reactivepowerlimits,
                                powerfactor)

"""
    TechThermal(activepower::Float64,
            activepowerlimits::NamedTuple{(:min, :max),Tuple{Float64,Float64}},
            reactivepower::Union{Float64,Nothing},
            reactivepowerlimits::Union{(min::Float64,max::Float64),Nothing},
            ramplimits::Union{NamedTuple{(:up, :down),Tuple{Float64,Float64}},Nothing},
            timelimits::Union{NamedTuple{(:min, :max),Tuple{Float64,Float64}},Nothing})

Data Structure for the economical parameters of thermal generation technologies.
    The data structure can be called calling all the fields directly or using named fields.
    Two examples [DOCTODO only one example exists below] are provided one with minimal data definition and a more comprenhensive one

    # Examples

    ```jldoctest

    julia> Tech = TechThermal(activepower = 100.0, activepowerlimits = (min = 50.0, max = 200.0))
   [ Info: 'Reactive Power' limits defined as nothing
    TechThermal:
       activepower: 100.0
       activepowerlimits: (min = 50.0, max = 200.0)
       reactivepower: nothing
       reactivepowerlimits: nothing
       ramplimits: nothing
       timelimits: nothing

"""
struct TechThermal <: TechnicalParams
    activepower::Float64 # [MW]
    activepowerlimits::NamedTuple{(:min, :max),Tuple{Float64,Float64}} # [MW]
    reactivepower::Union{Float64,Nothing} # [MVAr]
    reactivepowerlimits::Union{NamedTuple{(:min, :max),Tuple{Float64,Float64}},Nothing} # [MVAr]
    ramplimits::Union{NamedTuple{(:up, :down),Tuple{Float64,Float64}},Nothing} #MW/Hr
    timelimits::Union{NamedTuple{(:up, :down),Tuple{Float64,Float64}},Nothing} #Hr
    function TechThermal(activepower, activepowerlimits, reactivepower, reactivepowerlimits, ramplimits, timelimits)

        new(activepower, PowerSystems.orderedlimits(activepowerlimits, "Real Power"), reactivepower, PowerSystems.orderedlimits(reactivepowerlimits, "Reactive Power"), ramplimits, timelimits)

    end
end

# DOCTODO document this constructor
TechThermal(; activepower = 0.0,
            activepowerlimits = (min = 0.0, max = 0.0),
            reactivepower = nothing,
            reactivepowerlimits = nothing,
            ramplimits = nothing,
            timelimits = nothing
            ) = TechThermal(activepower, activepowerlimits, reactivepower,
                            reactivepowerlimits, ramplimits, timelimits)


text:
@testset "Test Stopping Remote Control" begin
    src = StopRemoteControl()
    src2 = StopRemoteControl(user_check = false)
    cheap_src = cheap_stop_remote_control()
    cheap_src2 = cheap_stop_remote_control(user_check = false)

    for k in fieldnames(StopRemoteControl)
        if k == :cheap_check
            @test !getfield(src, k)
        else
            @test getfield(src, k)
        end
    end
    for k in fieldnames(StopRemoteControl)
        if k ∈ (:cheap_check, :user_check)
            @test !getfield(src2, k)
        else
            @test getfield(src2, k)
        end
    end
    for k in fieldnames(StopRemoteControl)
        if k ∈ (:unbounded_and_domain_x_check, :domain_check, :unbounded_problem_check, :main_pb_check)
            @test !getfield(cheap_src, k)
        else
            @test getfield(cheap_src, k)
        end
    end
    for k in fieldnames(StopRemoteControl)
        if k ∈ (:unbounded_and_domain_x_check, :domain_check, :unbounded_problem_check, :main_pb_check, :user_check)
            @test !getfield(cheap_src2, k)
        else
            @test getfield(cheap_src2, k)
        end
    end

    io = IOBuffer()
    show(io, src)

end


text:
# Iterative decomposition of a system into uncoupled surfaces.

"""
The iteration has encountered a degeneracy that it cannot handle.
"""
struct DegeneracyException <: Exception
    n::Int
end

Base.show(io::IO, ex::DegeneracyException) = show(io, "Cannot handle $(ex.n)-fold degeneracy.")

"""
Iterative decomposition of a system into uncoupled surfaces.
"""
struct IterativeDecomposition{S,S_,M}
    "System which was decomposed."
    sys::System{S_,M}

    "Degeneracies at iteration steps."
    degens::Vector{Int}

    "Transformation matrix (S_, S)."
    vs::Matrix{Float64}

    "Coefficients for the energy (S)."
    ws_energy::Vector{Float64}
    "Coefficients for the linear coupling (M, S)."
    ws_lin::Matrix{Float64}
    "Coefficients for the quadratic coupling (M, M, S)."
    ws_quad::Array{Float64,3}
end

"""
    IterativeDecomposition(sys::System, max_iter::Int; sigma_cutoff::Float64=0.01)

Generate an iterative decomposition of `sys` using no more than `max_iter`
iterations.

Iterations are stopped when the largest singular value drops below the initial
singular value multiplied by `sigma_cutoff`.
"""
function IterativeDecomposition(sys::System{S,M}, max_iter::Int; sigma_cutoff::Float64=0.01) where {S,M}
    all(sys.freq .== sys.freq[:, 1]) || throw(DomainError(sys.freq, "All surfaces must have the same frequencies."))

    # Remaining tensor values.
    x_energy = copy(sys.energy)
    x_lin = copy(sys.lin)
    x_quad = copy(sys.quad)

    # Generated vectors and coefficients.
    vs = Matrix{Float64}(undef, S, 0)
    ws_energy = Vector{Float64}(undef, 0)
    ws_lin = Matrix{Float64}(undef, M, 0)
    ws_quad = Array{Float64,3}(undef, M, M, 0)

    degens = Int[]

    first_sigma = nothing
    for _ in 1:max_iter
        # Construct and SVD the matrix of coefficients.
        x = zeros(S, S + S*M + S*M*M)
        for s1 in 1:S
            row = s1
            for s2 in 1:S
                col0 = s2
                x[row, col0] = x_energy[s1, s2]
                for m1 in 1:M
                    col1 = S + (s2-1)*M + m1
                    x[row, col1] = x_lin[m1, s1, s2]
                    for m2 in 1:M
                        col2 = S + S*M + (s2-1)*M*M + (m1-1)*M + m2
                        x[row, col2] = x_quad[m1, m2, s1, s2]
                    end
                end
            end
        end
        F = svd(x)
        sigma = F.S[1]

        # Stopping criterion.
        if first_sigma === nothing
            first_sigma = sigma
        elseif sigma < sigma_cutoff*first_sigma
            break
        end

        # Determine degeneracy.
        num_degen = 0
        for s in 1:S
            F.S[s] != sigma && break
            num_degen += 1
        end
        push!(degens, num_degen)

        # Compute the vectors.
        pre_vs = Tuple{Vector{Float64},Float64,Vector{Float64},Vector{Float64}}[]
        if num_degen == 1
            v = F.U[:, 1]
            push!(pre_vs, (v, +1, v, v))
        elseif num_degen == 2
            v1 = F.U[:, 1]
            push!(pre_vs, (v1, +1, v1, v1))
            v2 = F.U[:, 2]
            push!(pre_vs, (v2, +1, v2, v2))
            vp = (v1 .+ v2) ./ sqrt(2)
            push!(pre_vs, (vp, +1, v1, v2))
            vm = (v1 .- v2) ./ sqrt(2)
            push!(pre_vs, (vm, -1, v1, v2))
            # Check symmetry.
            v1' * x_energy * v2 == v2' * x_energy * v1 || error("Asymmetric x_energy")
            for m1 in 1:M
                v1' * x_lin[m1, :, :] * v2 == v2' * x_lin[m1, :, :] * v1 || error("Asymmetric x_lin")
                for m2 in 1:M
                    v1' * x_quad[m2, m1, :, :] * v2 == v2' * x_quad[m2, m1, :, :] * v1 || error("Asymmetric x_quad")
                end
            end
        else
            throw(DegeneracyException(num_degen))
        end

        # Save vectors and coefficients.
        for (v, s, vL, vR) in pre_vs
            w_energy = s * vL' * x_energy * vR
            w_lin = zeros(M)
            w_quad = zeros(M, M)
            for m1 in 1:M
                w_lin[m1] = s * vL' * x_lin[m1, :, :] * vR
                for m2 in 1:M
                    w_quad[m2, m1] = s * vL' * x_quad[m2, m1, :, :] * vR
                end
            end
            @cat(vs, v)
            @cat(ws_energy, w_energy)
            @cat(ws_lin, w_lin)
            @cat(ws_quad, w_quad)
        end

        # Update arrays.
        for n in 1:length(pre_vs)
            v = pre_vs[n][1]
            idx = size(vs, 2) - length(pre_vs) + n
            x_energy .-= ws_energy[idx] * v * v'
            for m1 in 1:M
                x_lin[m1, :, :] .-= ws_lin[m1, idx] * v * v'
                for m2 in 1:M
                    x_quad[m2, m1, :, :] .-= ws_quad[m2, m1, idx] * v * v'
                end
            end
        end
    end

    # Number of generated surfaces.
    S_iter = length(ws_energy)

    # Check that we recover the original tensors.
    recovered_energy = x_energy + sum([ws_energy[n]*vs[:, n]*vs[:, n]' for n in 1:S_iter])
    err = maximum(abs.((sys.energy - recovered_energy)))
    err < 1e-14 || error("Maximum deviation: $(err)")
    for m1 in 1:M
        recovered_lin = x_lin[m1, :, :] + sum([ws_lin[m1, n]*vs[:, n]*vs[:, n]' for n in 1:S_iter])
        err = maximum(abs.((sys.lin[m1, :, :] - recovered_lin)))
        err < 1e-14 || error("Maximum deviation: $(err)")
        for m2 in 1:M
            recovered_quad = x_quad[m2, m1, :, :] + sum([ws_quad[m2, m1, n]*vs[:, n]*vs[:, n]' for n in 1:S_iter])
            err = maximum(abs.((sys.quad[m2, m1, :, :] - recovered_quad)))
            err < 1e-14 || error("Maximum deviation: $(err)")
        end
    end

    IterativeDecomposition{S_iter,S,M}(sys, degens, vs, ws_energy, ws_lin, ws_quad)
end

"""
    DiagonalSystem(decomp::IterativeDecomposition)

Generate a `DiagonalSystem` from `decomp`.
"""
function DiagonalSystem(decomp::IterativeDecomposition{S,S_,M}) where {S,S_,M}
    energy = zeros(S, S)
    freq = zeros(M, S)
    lin = zeros(M, S, S)
    quad = zeros(M, M, S, S)
    for s in 1:S
        energy[s, s] = decomp.ws_energy[s]
        for m1 in 1:M
            freq[m1, s] = decomp.sys.freq[m1, 1]
            lin[m1, s, s] = decomp.ws_lin[m1, s]
            for m2 in 1:M
                quad[m1, m2, s, s] = decomp.ws_quad[m1, m2, s]
            end
        end
    end
    DiagonalSystem(energy, freq, lin, quad)
end


text:
module TestUtilities

using Test

export @testset_accelerated, include_everywhere

using ComputationalResources
using ComputationalResources: CPUProcesses

macro testset_accelerated(name::String, var, ex)
    testset_accelerated(name, var, ex)
end
macro testset_accelerated(name::String, var, opts::Expr, ex)
    testset_accelerated(name, var, ex; eval(opts)...)
end
function testset_accelerated(name::String, var, ex; exclude=[])
    final_ex = quote
        $var = CPU1()
        @testset $name $ex
    end
    resources = Any[CPUProcesses()]
    @static if VERSION >= v"1.3.0-DEV.573"
        push!(resources, CPUThreads())
    end
    for res in resources
        if any(x->typeof(res)<:x, exclude)
            push!(final_ex.args, quote
               $var = $res
               @testset $(name*" ($(typeof(res).name))") begin
                   @test_broken false
               end
            end)
        else
            push!(final_ex.args, quote
               $var = $res
               @testset $(name*" ($(typeof(res).name))") $ex
            end)
        end
    end
    return esc(final_ex)
end

end


text:
"""
  IsingJ

  ver 0.2.2 Released Sep 2020
  This code is distributed under the constitution
  of MIT LICENSE.

  module Initial

  2019/04/24  Released by PearCandy                 """

module Initial
    export random_start,cold_start

    #random start : initialization of spin
    function random_start(N)
        spin=ones(N,N)
        for i in 1:N
            for j in 1:N
                if(rand()<=0.5)
                    spin[i,j]=-1
                else
                    spin[i,j]=1
                end
            end
        end
        return spin
    end

    #cold start : initialization of spin
    function cold_start(N)
        spin=ones(N,N)
        return spin
    end

end # module


text:
using ITensors
using ITensorChemistry
using OptimKit
using Random
using Zygote

# General definitions
ITensors.inner(ψ, Hs::Vector, ϕ; kwargs...) = sum([inner(ψ, H, ϕ; kwargs...) for H in Hs])
⊗(x...) = kron(reverse(x)...)

molecule = "H₂O"
basis = "sto-3g"

@show molecule
@show basis

# Split the Hamiltonian into `nparts` sub-Hamiltonians
nparts = 2

println("\nRunning Hartree-Fock")
(; hamiltonian, state, hartree_fock_energy) = @time molecular_orbital_hamiltonian(nparts; molecule, basis)
println("Hartree-Fock complete")

println("Basis set size = ", length(state))

s = siteinds("Electron", length(state); conserve_qns=true)

println("\nConstruct MPO")

H = Vector{MPO}(undef, nparts)
@time for n in 1:nparts
  H[n] = splitblocks(linkinds, MPO(hamiltonian[n], s))
end
println("MPO constructed")

@show maxlinkdim.(H)

ψhf = MPS(s, state)

@show inner(ψhf, H, ψhf)
@show hartree_fock_energy

sweeps = Sweeps(4)
setmaxdim!(sweeps, 100, 200)
setcutoff!(sweeps, 1e-6)
setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)

println("\nRunning DMRG")
@show sweeps

e, ψ = dmrg(H, ψhf, sweeps)
println("DMRG complete")

function ITensors.op(::OpName"ucc_1e", st::SiteType"Electron"; f)
  T = zeros(4^2, 4^2)
  for σ in ["↑", "↓"]
    T += op("c†$σ", st) ⊗ op("c$σ", st)
    T -= op("c$σ", st) ⊗ op("c†$σ", st)
  end
  return f(T)
end

function ITensors.op(::OpName"ucc_2e", st::SiteType"Electron"; f)
  T = zeros(4^4, 4^4)
  for σ1 in ["↑", "↓"], σ2 in ["↑", "↓"]
    T += op("c†$σ1", st) ⊗ op("c†$σ2", st) ⊗ op("c$σ2", st) ⊗ op("c$σ1", st)
    T -= op("c$σ1", st) ⊗ op("c$σ2", st) ⊗ op("c†$σ2", st) ⊗ op("c†$σ1", st)
  end
  return f(T)
end

# Get the number of parameters
nα = length(s)
sites_1e = [(i, j) for i in 1:nα, j in 1:nα if i < j]
sites_2e = [(i, j, k, l) for i in 1:nα, j in 1:nα, k in 1:nα, l in 1:nα if i < j < k < l]

# TODO: implement a simpler, more local circuit, like one from:
# https://arxiv.org/abs/1805.04340
function ucc_circuit(s, t)
  nt_1e = length(sites_1e)
  nt_2e = length(sites_2e)
  t_1e = t[1:nt_1e]
  t_2e = t[(nt_1e + 1):end]

  U_1e = [
    (
      f(x) = exp(t_1e[i] * x);
      n⃗ = sites_1e[i];
      s⃗ = [s[n] for n in n⃗];
      op("ucc_1e", s⃗...; f=f)
    )
    for i in eachindex(sites_1e)
  ]
  U_2e = [
    (
      f(x) = exp(t_2e[i] * x);
      n⃗ = sites_2e[i];
      s⃗ = [s[n] for n in n⃗];
      op("ucc_2e", s⃗...; f=f)
    )
    for i in eachindex(sites_2e)
  ]
  return [U_1e; U_2e]
end

# Directly perform VQE.
#function loss(t)
#  U = ucc_circuit(s, t)
#  Uψhf = apply(U, ψhf; cutoff=1e-2)
#  return inner(Uψhf', H, Uψhf; cutoff=1e-2)
#end

# Quantum state preparation using DMRG
# as the target state.
function loss(t)
  U = ucc_circuit(s, t)
  Uψhf = apply(U, ψhf; cutoff=1e-6)
  return -abs2(inner(ψ, Uψhf))
end

nt_1e = length(sites_1e)
nt_2e = length(sites_2e)

# Start at HF GS
# When trying to optimize starting from this state,
# for some reason the paramaters go to NaN.
# Need to debug.
#t0_1e = zeros(nt_1e)
#t0_2e = zeros(nt_2e)

Random.seed!(1234)
random_scale = 0.1
t0_1e = random_scale * randn(nt_1e)
t0_2e = random_scale * randn(nt_2e)

t0 = [t0_1e; t0_2e]

U0 = ucc_circuit(s, t0)

U0ψhf = apply(U0, ψhf)
@show inner(U0ψhf, H, U0ψhf)

@show loss(t0)

loss_∇loss(t) = (loss(t), convert(Vector, loss'(t)))
algorithm = LBFGS(; gradtol=1e-3, verbosity=2)
tₒₚₜ, lossₒₚₜ, ∇lossₒₚₜ, numfg, normgradhistory = optimize(loss_∇loss, t0, algorithm)

@show loss(tₒₚₜ)

Ut = ucc_circuit(s, tₒₚₜ)
Utψhf = apply(Ut, ψhf; cutoff=1e-6)
@show inner(Utψhf', H, Utψhf)


text:
"""
    FullyObservablePOMDP(mdp)

Turn `MDP` `mdp` into a `POMDP` where the observations are the states of the MDP.
"""
struct FullyObservablePOMDP{M,S,A} <: POMDP{S,A,S}
    mdp::M
end

function FullyObservablePOMDP(m::MDP)
    return FullyObservablePOMDP{typeof(m), statetype(m), actiontype(m)}(m)
end

mdptype(::Type{FullyObservablePOMDP{M,S,A}}) where {M,S,A} = M

POMDPs.observations(pomdp::FullyObservablePOMDP) = states(pomdp.mdp)
POMDPs.obsindex(pomdp::FullyObservablePOMDP{S, A}, o::S) where {S, A} = stateindex(pomdp.mdp, o)

POMDPs.convert_o(T::Type{V}, o, pomdp::FullyObservablePOMDP) where {V<:AbstractArray} = convert_s(T, o, pomdp.mdp)
POMDPs.convert_o(T::Type{S}, vec::V, pomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, pomdp.mdp)

function POMDPs.observation(pomdp::FullyObservablePOMDP, a, sp)
    return Deterministic(sp)
end

function POMDPs.observation(pomdp::FullyObservablePOMDP, s, a, sp)
    return Deterministic(sp)
end

# inherit other function from the MDP type

POMDPs.states(pomdp::FullyObservablePOMDP) = states(pomdp.mdp)
POMDPs.actions(pomdp::FullyObservablePOMDP) = actions(pomdp.mdp)
POMDPs.transition(pomdp::FullyObservablePOMDP, s, a) = transition(pomdp.mdp, s, a)
POMDPs.isterminal(pomdp::FullyObservablePOMDP, s) = isterminal(pomdp.mdp, s)
POMDPs.discount(pomdp::FullyObservablePOMDP) = discount(pomdp.mdp)
POMDPs.stateindex(pomdp::FullyObservablePOMDP, s) = stateindex(pomdp.mdp, s)
POMDPs.actionindex(pomdp::FullyObservablePOMDP, a) = actionindex(pomdp.mdp, a)
POMDPs.convert_s(T::Type{V}, s, pomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_s(T, s, pomdp.mdp)
POMDPs.convert_s(T::Type{S}, vec::V, pomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, pomdp.mdp)
POMDPs.convert_a(T::Type{V}, a, pomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_a(T, a, pomdp.mdp)
POMDPs.convert_a(T::Type{A}, vec::V, pomdp::FullyObservablePOMDP) where {A,V<:AbstractArray} = convert_a(T, vec, pomdp.mdp)
POMDPs.reward(pomdp::FullyObservablePOMDP, s, a) = reward(pomdp.mdp, s, a)
POMDPs.initialstate(m::FullyObservablePOMDP) = initialstate(m.mdp)
POMDPs.initialobs(m::FullyObservablePOMDP, s) = Deterministic(s)

max_stars_repo_name: JuliaPOMDP/POMDPModelTools

max_issues_repo_name: JuliaPOMDP/POMDPModelTools

max_forks_repo_name: JuliaPOMDP/POMDPModelTools


text:
print("The_Hedge_Fund_Industry.jl here")

greet() = print("The Hedge Fund Industry");


text:
module VerySimpleRasters

using HDF5
using RecipesBase
using Mmap
using Dates

import Base: getindex, setindex!

include("datatypes.jl")
include("importASCII.jl")
include("grd_file.jl")
include("interface.jl")
include("operations.jl")
include("plotrecipes.jl")

export VerySimpleRaster, importASCII, write_grd
export extract, aggregate, crop, mask
export writeraster
export eachcoordinate, coordinates, coordinate_to_index, index_to_coordinate
end # module


text:
@cache struct EulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

@cache struct SplitEulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  SplitEulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct SplitEulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SplitEulerConstantCache()

function alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  EulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct EulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = EulerConstantCache()

@cache struct HeunCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
end

@cache struct RalstonCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  HeunCache(u,uprev,similar(u),similar(u, uEltypeNoUnits),zero(rate_prototype),
            zero(rate_prototype))
end

function alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  RalstonCache(u,uprev,similar(u),similar(u, uEltypeNoUnits),zero(rate_prototype),
               zero(rate_prototype))
end

struct HeunConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = HeunConstantCache()

struct RalstonConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RalstonConstantCache()

@cache struct MidpointCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
end

struct MidpointConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  MidpointCache(u,uprev,k,tmp,atmp,fsalfirst)
end

alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = MidpointConstantCache()

@cache struct RK4Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k₂::rateType
  k₃::rateType
  k₄::rateType
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
end

struct RK4ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  k₄ = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  RK4Cache(u,uprev,k₁,k₂,k₃,k₄,k,tmp,atmp)
end

alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RK4ConstantCache()

@cache struct CarpenterKennedy2N54Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct CarpenterKennedy2N54ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  A2::T
  A3::T
  A4::T
  A5::T
  B1::T
  B2::T
  B3::T
  B4::T
  B5::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2

  function CarpenterKennedy2N54ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    A2 = T(-567301805773/1357537059087)
    A3 = T(-2404267990393/2016746695238)
    A4 = T(-3550918686646/2091501179385)
    A5 = T(-1275806237668/842570457699)
    B1 = T(1432997174477/9575080441755)
    B2 = T(5161836677717/13612068292357)
    B3 = T(1720146321549/2090206949498)
    B4 = T(3134564353537/4481467310338)
    B5 = T(2277821191437/14882151754819)
    c2 = T2(1432997174477/9575080441755)
    c3 = T2(2526269341429/6820363962896)
    c4 = T2(2006345519317/3224310063776)
    c5 = T2(2802321613138/2924317926251)
    new{T,T2}(A2, A3, A4, A5, B1, B2, B3, B4, B5, c2, c3, c4, c5)
  end
end

function alg_cache(alg::CarpenterKennedy2N54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CarpenterKennedy2N54ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  CarpenterKennedy2N54Cache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::CarpenterKennedy2N54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  CarpenterKennedy2N54ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
end

@cache struct BS3Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen3Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen4Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen4Cache(u,uprev,k1,k2,k3,k4,k5,k6,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct BS5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Tsit5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

@cache struct RK46NLCache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct RK46NLConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α2::T
  α3::T
  α4::T
  α5::T
  α6::T
  β1::T
  β2::T
  β3::T
  β4::T
  β5::T
  β6::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2

  function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α2 = T(-0.737101392796)
    α3 = T(-1.634740794343)
    α4 = T(-0.744739003780)
    α5 = T(-1.469897351522)
    α6 = T(-2.813971388035)
    β1 = T(0.032918605146)
    β2 = T(0.823256998200)
    β3 = T(0.381530948900)
    β4 = T(0.200092213184)
    β5 = T(1.718581042715)
    β6 = T(0.27)
    c2 = T2(0.032918605146)
    c3 = T2(0.249351723343)
    c4 = T2(0.466911705055)
    c5 = T2(0.582030414044)
    c6 = T2(0.847252983783)
    new{T,T2}(α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6)
  end
end

function alg_cache(alg::RK46NL,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = RK46NLConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  RK46NLCache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::RK46NL,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  RK46NLConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct DP5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = k2
  k7 = zero(rate_prototype) # This is FSAL'd to k1
  dense_tmp3 = k2
  dense_tmp4 = k5
  bspl = k3

  tmp = similar(u) # has to be separate for FSAL
  utilde = tmp

  if eltype(u) != uEltypeNoUnits || calck
    update = zero(rate_prototype)
    atmp = similar(u,uEltypeNoUnits)
  else
    update = k7
    atmp = k3
  end

  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  cache = DP5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
  cache
end

alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct DP5ThreadedCache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::DP5Threaded,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  dense_tmp3 = zero(rate_prototype)
  dense_tmp4 = zero(rate_prototype)
  update = zero(rate_prototype)
  bspl = zero(rate_prototype)
  utilde = similar(u)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits)
  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  DP5ThreadedCache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
end

@cache struct Anas5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Anas5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct CFRLDDRK64Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct CFRLDDRK64ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α1::T
  α2::T
  α3::T
  α4::T
  α5::T
  β1::T
  β2::T
  β3::T
  β4::T
  β5::T
  β6::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2

  function CFRLDDRK64ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α1 = T(0.17985400977138)
    α2 = T(0.14081893152111)
    α3 = T(0.08255631629428)
    α4 = T(0.65804425034331)
    α5 = T(0.31862993413251)
    β1 = T(0.10893125722541)
    β2 = T(0.13201701492152)
    β3 = T(0.38911623225517)
    β4 = T(-0.59203884581148)
    β5 = T(0.47385028714844)
    β6 = T(0.48812405426094)
    c2 = T2(0.28878526699679)
    c3 = T2(0.38176720366804)
    c4 = T2(0.71262082069639)
    c5 = T2(0.69606990893393)
    c6 = T2(0.83050587987157)
    new{T,T2}(α1, α2, α3, α4, α5, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6)
  end
end

function alg_cache(alg::CFRLDDRK64,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  CFRLDDRK64Cache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::CFRLDDRK64,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  CFRLDDRK64ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct NDBLSRK124Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct NDBLSRK124ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α2::T
  α3::T
  α4::T
  α5::T
  α6::T
  α7::T
  α8::T
  α9::T
  α10::T
  α11::T
  α12::T
  β1::T
  β2::T
  β3::T
  β4::T
  β5::T
  β6::T
  β7::T
  β8::T
  β9::T
  β10::T
  β11::T
  β12::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2
  c8::T2
  c9::T2
  c10::T2
  c11::T2
  c12::T2

  function NDBLSRK124ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α2  = T(-0.0923311242368072)
    α3  = T(-0.9441056581158819)
    α4  = T(-4.3271273247576394)
    α5  = T(-2.1557771329026072)
    α6  = T(-0.9770727190189062)
    α7  = T(-0.7581835342571139)
    α8  = T(-1.7977525470825499)
    α9  = T(-2.6915667972700770)
    α10 = T(-4.6466798960268143)
    α11 = T(-0.1539613783825189)
    α12 = T(-0.5943293901830616)
    β1  = T(0.0650008435125904)
    β2  = T(0.0161459902249842)
    β3  = T(0.5758627178358159)
    β4  = T(0.1649758848361671)
    β5  = T(0.3934619494248182)
    β6  = T(0.0443509641602719)
    β7  = T(0.2074504268408778)
    β8  = T(0.6914247433015102)
    β9  = T(0.3766646883450449)
    β10 = T(0.0757190350155483)
    β11 = T(0.2027862031054088)
    β12 = T(0.2167029365631842)
    c2  = T2(0.0650008435125904)
    c3  = T2(0.0796560563081853)
    c4  = T2(0.1620416710085376)
    c5  = T2(0.2248877362907778)
    c6  = T2(0.2952293985641261)
    c7  = T2(0.3318332506149405)
    c8  = T2(0.4094724050198658)
    c9  = T2(0.6356954475753369)
    c10 = T2(0.6806551557645497)
    c11 = T2(0.7143773712418350)
    c12 = T2(0.9032588871651854)
    new{T,T2}(α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12)
  end
end

function alg_cache(alg::NDBLSRK124,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = NDBLSRK124ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  NDBLSRK124Cache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::NDBLSRK124,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  NDBLSRK124ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct NDBLSRK134Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
    u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct NDBLSRK134ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α2::T
  α3::T
  α4::T
  α5::T
  α6::T
  α7::T
  α8::T
  α9::T
  α10::T
  α11::T
  α12::T
  α13::T
  β1::T
  β2::T
  β3::T
  β4::T
  β5::T
  β6::T
  β7::T
  β8::T
  β9::T
  β10::T
  β11::T
  β12::T
  β13::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2
  c8::T2
  c9::T2
  c10::T2
  c11::T2
  c12::T2
  c13::T2

  function NDBLSRK134ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α2  = T(-0.6160178650170565)
    α3  = T(-0.4449487060774118)
    α4  = T(-1.0952033345276178)
    α5  = T(-1.2256030785959187)
    α6  = T(-0.2740182222332805)
    α7  = T(-0.0411952089052647)
    α8  = T(-0.1797084899153560)
    α9  = T(-1.1771530652064288)
    α10 = T(-0.4078831463120878)
    α11 = T(-0.8295636426191777)
    α12 = T(-4.7895970584252288)
    α13 = T(-0.6606671432964504)
    β1  = T(0.0271990297818803)
    β2  = T(0.1772488819905108)
    β3  = T(0.0378528418949694)
    β4  = T(0.6086431830142991)
    β5  = T(0.2154313974316100)
    β6  = T(0.2066152563885843)
    β7  = T(0.0415864076069797)
    β8  = T(0.0219891884310925)
    β9  = T(0.9893081222650993)
    β10 = T(0.0063199019859826)
    β11 = T(0.3749640721105318)
    β12 = T(1.6080235151003195)
    β13 = T(0.0961209123818189)
    c2  = T2(0.0271990297818803)
    c3  = T2(0.0952594339119365)
    c4  = T2(0.1266450286591127)
    c5  = T2(0.1825883045699772)
    c6  = T2(0.3737511439063931)
    c7  = T2(0.5301279418422206)
    c8  = T2(0.5704177433952291)
    c9  = T2(0.5885784947099155)
    c10 = T2(0.6160769826246714)
    c11 = T2(0.6223252334314046)
    c12 = T2(0.6897593128753419)
    c13 = T2(0.9126827615920843)
    new{T,T2}(α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13)
  end
end

function alg_cache(alg::NDBLSRK134,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = NDBLSRK134ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  NDBLSRK134Cache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::NDBLSRK134,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  NDBLSRK134ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct NDBLSRK144Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct NDBLSRK144ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α2::T
  α3::T
  α4::T
  α5::T
  α6::T
  α7::T
  α8::T
  α9::T
  α10::T
  α11::T
  α12::T
  α13::T
  α14::T
  β1::T
  β2::T
  β3::T
  β4::T
  β5::T
  β6::T
  β7::T
  β8::T
  β9::T
  β10::T
  β11::T
  β12::T
  β13::T
  β14::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2
  c8::T2
  c9::T2
  c10::T2
  c11::T2
  c12::T2
  c13::T2
  c14::T2

  function NDBLSRK144ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α2  = T(-0.7188012108672410)
    α3  = T(-0.7785331173421570)
    α4  = T(-0.0053282796654044)
    α5  = T(-0.8552979934029281)
    α6  = T(-3.9564138245774565)
    α7  = T(-1.5780575380587385)
    α8  = T(-2.0837094552574054)
    α9  = T(-0.7483334182761610)
    α10 = T(-0.7032861106563359)
    α11 = T(0.0013917096117681)
    α12 = T(-0.0932075369637460)
    α13 = T(-0.9514200470875948)
    α14 = T(-7.1151571693922548)
    β1  = T(0.0367762454319673)
    β2  = T(0.3136296607553959)
    β3  = T(0.1531848691869027)
    β4  = T(0.0030097086818182)
    β5  = T(0.3326293790646110)
    β6  = T(0.2440251405350864)
    β7  = T(0.3718879239592277)
    β8  = T(0.6204126221582444)
    β9  = T(0.1524043173028741)
    β10 = T(0.0760894927419266)
    β11 = T(0.0077604214040978)
    β12 = T(0.0024647284755382)
    β13 = T(0.0780348340049386)
    β14 = T(5.5059777270269628)
    c2  = T2(0.0367762454319673)
    c3  = T2(0.1249685262725025)
    c4  = T2(0.2446177702277698)
    c5  = T2(0.2476149531070420)
    c6  = T2(0.2969311120382472)
    c7  = T2(0.3978149645802642)
    c8  = T2(0.5270854589440328)
    c9  = T2(0.6981269994175695)
    c10 = T2(0.8190890835352128)
    c11 = T2(0.8527059887098624)
    c12 = T2(0.8604711817462826)
    c13 = T2(0.8627060376969976)
    c14 = T2(0.8734213127600976)
    new{T,T2}(α2, α3, α4, α5, α6, α7, α8, α9, α10, α11, α12, α13, α14, β1, β2, β3, β4, β5, β6, β7, β8, β9, β10, β11, β12, β13, β14, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14)
  end
end

function alg_cache(alg::NDBLSRK144,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = NDBLSRK144ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  NDBLSRK144Cache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::NDBLSRK144,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  NDBLSRK144ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end


text:
using Test
using Transpiler

@testset "Examples" begin
    include("SCAMPy.jl")
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(γ::EmpiricalVariogram; showbins=true)
  # get the data
  x, y, n = values(γ)

  binsize = x[2] - x[1]

  # discard empty bins
  x = x[n .> 0]; y = y[n .> 0]; n = n[n .> 0]

  # draw bin counts as a measure of confidence
  if showbins
    @series begin
      # plot a "frequency" instead of raw counts
      f = n*(maximum(y) / maximum(n)) / 10.

      seriestype := :bar
      seriescolor := :black
      fillalpha := .1
      label := "bin counts"

      x, f
    end
  end

  seriestype --> :scatterpath
  linecolor --> :black
  xlims --> (0, maximum(x) + binsize/2)
  xguide --> "h"
  yguide --> "γ(h)"
  label --> "variogram"

  x, y
end


text:
using DelimitedFiles
using StaticArrays, Parameters
using OrdinaryDiffEq, Plots

include("eom.jl")
include("parameters.jl")
include("specifieds.jl")
include("functions.jl")
include("contact_force.jl")
include("../plotting.jl")
include("setup.jl")

# set up model
p, u₀ = setup()

# time span
tspan = (0.0, 0.484)

# set up problem
prob = ODEProblem(eom, u₀, tspan, p)

# solve
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat=0.001)

# animate
animate_model(sol)


text:
# This file is a part of Pukeko.jl.
# License is MIT. https://github.com/IainNZ/Pukeko.jl

import Pukeko

module PukekoTests

    using Pukeko

    parametric_by_function_called = 0
    function parametric_by_function(value)
        global parametric_by_function_called += 1
        @test value in ["foo", "bar"]
        @test parametric_by_function_called >= 1
    end
    Pukeko.parametric(PukekoTests, parametric_by_function, ["foo", "bar"])

    parametric_by_macro_called = 0
    function parametric_by_macro(value)
        global parametric_by_macro_called += 1
        @test value in ["foo", "bar"]
        @test parametric_by_macro_called >= 1
    end
    Pukeko.@parametric parametric_by_macro ["foo", "bar"]

    parametric_many_called = 0
    function parametric_many(value1, value2)
        global parametric_many_called += 1
        @test value1 + 1 == value2
        @test parametric_many_called >= 1
    end
    Pukeko.@parametric parametric_many ((1, 2), (3, 4), (5, 6))

    function __test_equal()
        @test 1 == 1
    end
end

Pukeko.run_tests(PukekoTests)
@assert PukekoTests.parametric_by_function_called == 2
@assert PukekoTests.parametric_by_macro_called == 2
@assert PukekoTests.parametric_many_called == 3

module TestThrowsTests

    using Pukeko

    function test_throw_stuff()
        @test_throws DomainError log(-1)
    end

    function test_throw_message()
        err = ErrorException("Err Msg")

        function i_throw()
            throw(ErrorException("Err Msg"))
        end

        @test_throws err i_throw()
    end
end
Pukeko.run_tests(TestThrowsTests)

module FailFastTests
    using Pukeko
    function test_fail()
        @test 1 == 2
    end
    should_be_zero = 0
    function test_not_run()
        global should_be_zero += 1
    end
end

println("Message about failing test should appear next")
should_not_run = true
try
    Pukeko.run_tests(FailFastTests, fail_fast=true)
    global should_not_run = false
catch
    @assert should_not_run
end
@assert should_not_run
@assert FailFastTests.should_be_zero == 0

module FailSlowTests
    using Pukeko
    function test_fail()
        @test 1 == 2
    end
    should_be_one = 0
    function test_run()
        global should_be_one += 1
    end
end

should_not_run = true
try
    Pukeko.run_tests(FailSlowTests)
    global should_not_run = false
catch
    @assert should_not_run
end
@assert should_not_run
@assert FailSlowTests.should_be_one == 1

module Empty
end
Pukeko.run_tests(Empty)

module OnlyCallFunctions
    test_variable = 2
    module test_stuff
    end
end
Pukeko.run_tests(OnlyCallFunctions)

include("base_int.jl")


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using DifferentialEquations

function main()

# # Model Specifications
# - SH num of high risk susceptible
# - IH num of high risk infected
# - SL num of low risk susceptible
# - IL num of low risk infected
# # Parameters
# - beta represents the determines the contact and transmission rates
# - gamma the rate at which treatment is sought

function sis_ode(du,u,p,t)
        SH,IH,SL,IL = u
        betaHH,betaHL,betaLH,betaLL,gamma=p
        du[1]=-(betaHH*IH+betaHL*IL)*SH+gamma*IH
        du[2]=+(betaHH*IH+betaHL*IL)*SH-gamma*IH
        du[3]=-(betaLH*IH+betaLL*IL)*SL+gamma*IL
        du[4]=+(betaLH*IH+betaLL*IL)*SL-gamma*IL
end

parms =[10,0.1,0.1,1,1]
init=[0.19999,0.00001,0.799,0.001]
tspan=tspan = (0.0,15.0)

sis_prob = ODEProblem(sis_ode,init,tspan,parms)
sis_sol = solve(sis_prob,saveat=0.1);

end