text:
using Libdl, Dates, Printf, DecFP, WeakRefStrings

# Link to ODBC Driver Manager (system-dependent)
const odbc_dm = let
    if !isdefined(@__MODULE__, :odbc_dm)
        Sys.islinux()   && (lib_choices = ["libodbc", "libodbc.so", "libodbc.so.1", "libodbc.so.2", "libodbc.so.3"])
        Sys.iswindows() && (lib_choices = ["odbc32"])
        Sys.isapple()   && (lib_choices = ["libodbc.2.dylib","libodbc.dylib","libiodbc","libiodbc.dylib","libiodbc.1.dylib","libiodbc.2.dylib","libiodbc.3.dylib"])
        Libdl.find_library(lib_choices)
    end
end

function setODBC(x)
    global odbc_dm
    odbc_dm = x
end

# Translation of sqltypes.h; C typealiases for SQL functions
# http://msdn.microsoft.com/en-us/library/windows/desktop/ms716298(v=vs.85).aspx
# http://msdn.microsoft.com/en-us/library/windows/desktop/aa383751(v=vs.85).aspx
const SQLCHAR =       UInt8
const SQLSMALLINT =   Cshort
const SQLUSMALLINT =  Cushort

const SQLSCHAR =      Cchar
const SQLDATE =       Cuchar
const SQLDECIMAL =    Cuchar
const SQLDOUBLE =     Cdouble
const SQLFLOAT =      Cdouble

const SQLVARCHAR =    Cuchar
const SQLNUMERIC =    Cuchar
const SQLREAL =       Cfloat
const SQLTIME =       Cuchar
const SQLTIMESTAMP =  Cuchar

if occursin("iodbc", odbc_dm)
    const SQLWCHAR = UInt32
else
    # correct for windows + unixODBC
    const SQLWCHAR = Cushort
end

# ODBC API	64-bit platform	32-bit platform
# SQLINTEGER	32 bits	32 bits
# SQLUINTEGER	32 bits	32 bits
# SQLLEN	64 bits	32 bits
# SQLULEN	64 bits	32 bits
# SQLSETPOSIROW	64 bits	16 bits
# SQL_C_BOOKMARK	64 bits	32 bits
# BOOKMARK	64 bits	32 bits

const SQLINTEGER =  Cint
const SQLUINTEGER = Cuint
const SQLLEN = Int
const SQLULEN = UInt

# if WORD_SIZE == 64
#     const SQLINTEGER =    Cint
#     const SQLUINTEGER =   Cuint
# else
#     const SQLINTEGER =    Clong
#     const SQLUINTEGER =   Culong
# end
#
# const SQLLEN =        SQLINTEGER
# const SQLULEN =       SQLUINTEGER
const SQLSETPOSIROW = SQLUSMALLINT

const SQLROWCOUNT =   SQLULEN
const SQLROWSETSIZE = SQLULEN
const SQLTRANSID =    SQLULEN
const SQLROWOFFSET =  SQLLEN
const SQLPOINTER =    Ptr{Cvoid}
const SQLRETURN =     SQLSMALLINT
const SQLHANDLE =     Ptr{Cvoid}
const SQLHENV =       SQLHANDLE
const SQLHDBC =       SQLHANDLE
const SQLHSTMT =      SQLHANDLE
const SQLHDESC =      SQLHANDLE
const ULONG =         Cuint
const PULONG =        Ptr{ULONG}
const USHORT =        Cushort
const PUSHORT =       Ptr{USHORT}
const UCHAR =         Cuchar
const PUCHAR =        Ptr{Cuchar}
const PSZ =           Ptr{Cchar}
const SCHAR =         Cchar
const SDWORD =        Cint
const SWORD =         Cshort
const UDWORD =        Cuint
const UWORD =         Cushort
const SLONG =         Cint
const SSHORT =        Cshort
const SDOUBLE =       Cdouble
const LDOUBLE =       Cdouble
const SFLOAT =        Cfloat
const PTR =           Ptr{Cvoid}
const HENV =          Ptr{Cvoid}
const HDBC =          Ptr{Cvoid}
const HSTMT =         Ptr{Cvoid}
const RETCODE =       Cshort
const SQLHWND =       Ptr{Cvoid}

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

# provide lowercase conversion functions for all types
# e.g., sqlchar(x) = convert(SQLCHAR, x)

for t in [:SQLCHAR, :SQLSCHAR, :SQLWCHAR, :SQLDATE, :SQLDECIMAL,
          :SQLDOUBLE, :SQLFLOAT, :SQLINTEGER, :SQLUINTEGER,
          :SQLSMALLINT, :SQLUSMALLINT, :SQLLEN, :SQLULEN,
          :SQLSETPOSIROW, :SQLROWCOUNT, :SQLROWSETSIZE, :SQLTRANSID,
          :SQLROWOFFSET, :SQLNUMERIC, :SQLPOINTER, :SQLREAL, :SQLTIME,
          :SQLTIMESTAMP, :SQLVARCHAR, :SQLRETURN, :SQLHANDLE,
          :SQLHENV, :SQLHDBC, :SQLHSTMT, :SQLHDESC, :ULONG, :PULONG,
          :USHORT, :PUSHORT, :UCHAR, :PUCHAR, :PSZ, :SCHAR, :SDWORD,
          :SWORD, :UDWORD, :UWORD, :SLONG, :SSHORT, :SDOUBLE,
          :LDOUBLE, :SFLOAT, :PTR, :HENV, :HDBC, :HSTMT, :RETCODE,
          :SQLHWND]
    fn = Symbol(lowercase(string(t)))
    @eval $fn(x) = convert($t, x)
end

# Data Type Mappings
# SQL data types are returned in resultset metadata calls (ODBCMetadata)
# C data types are used in SQLBindCols (ODBCFetch) to allocate column memory; the driver manager converts from the SQL type to this C type in memory
# Julia types indicate how julia should read the returned C data type memory from the previous step

# SQL Data Type                     C Data Type                         Julia Type
# ---------------------------------------------------------------------------------
# SQL_CHAR                          SQL_C_CHAR                          UInt8
# SQL_VARCHAR                       SQL_C_CHAR                          UInt8
# SQL_LONGVARCHAR                   SQL_C_CHAR                          UInt8
# SQL_WCHAR                         SQL_C_WCHAR                         Cwchar_t
# SQL_WVARCHAR                      SQL_C_WCHAR                         Cwchar_t
# SQL_WLONGVARCHAR                  SQL_C_WCHAR                         Cwchar_t
# SQL_DECIMAL                       SQL_C_DOUBLE                        SQLNumeric
# SQL_NUMERIC                       SQL_C_DOUBLE                        SQLNumeric
# SQL_SMALLINT                      SQL_C_SHORT                         Int16
# SQL_INTEGER                       SQL_C_LONG                          Int32
# SQL_REAL                          SQL_C_FLOAT                         Float64
# SQL_FLOAT                         SQL_C_DOUBLE                        Float64
# SQL_DOUBLE                        SQL_C_DOUBLE                        Float64
# SQL_BIT                           SQL_C_BIT                           Int8
# SQL_TINYINT                       SQL_C_TINYINT                       Int8
# SQL_BIGINT                        SQL_C_BIGINT                        Int64
# SQL_BINARY                        SQL_C_BINARY                        UInt8
# SQL_VARBINARY                     SQL_C_BINARY                        UInt8
# SQL_LONGVARBINARY                 SQL_C_BINARY                        UInt8
# SQL_TYPE_DATE                     SQL_C_TYPE_DATE                     SQLDate
# SQL_TYPE_TIME                     SQL_C_TYPE_TIME                     SQLTime
# SQL_TYPE_TIMESTAMP                SQL_C_TYPE_TIMESTAMP                SQLTimestamp
# SQL_INTERVAL_MONTH                SQL_C_INTERVAL_MONTH                UInt8
# SQL_INTERVAL_YEAR                 SQL_C_INTERVAL_YEAR                 UInt8
# SQL_INTERVAL_YEAR_TO_MONTH        SQL_C_INTERVAL_YEAR_TO_MONTH        UInt8
# SQL_INTERVAL_DAY                  SQL_C_INTERVAL_DAY                  UInt8
# SQL_INTERVAL_HOUR                 SQL_C_INTERVAL_HOUR                 UInt8
# SQL_INTERVAL_MINUTE               SQL_C_INTERVAL_MINUTE               UInt8
# SQL_INTERVAL_SECOND               SQL_C_INTERVAL_SECOND               UInt8
# SQL_INTERVAL_DAY_TO_HOUR          SQL_C_INTERVAL_DAY_TO_HOUR          UInt8
# SQL_INTERVAL_DAY_TO_MINUTE        SQL_C_INTERVAL_DAY_TO_MINUTE        UInt8
# SQL_INTERVAL_DAY_TO_SECOND        SQL_C_INTERVAL_DAY_TO_SECOND        UInt8
# SQL_INTERVAL_HOUR_TO_MINUTE       SQL_C_INTERVAL_HOUR_TO_MINUTE       UInt8
# SQL_INTERVAL_HOUR_TO_SECOND       SQL_C_INTERVAL_HOUR_TO_SECOND       UInt8
# SQL_INTERVAL_MINUTE_TO_SECOND     SQL_C_INTERVAL_MINUTE_TO_SECOND     UInt8
# SQL_GUID                          SQL_C_GUID                          SQLGUID

# SQL Data Type Definitions
const SQL_NULL_DATA     = -1
const SQL_NO_TOTAL      = -4
const SQL_CHAR          = Int16(  1) # Character string of fixed string length n.
const SQL_VARCHAR       = Int16( 12) # Variable-length character string with a maximum string length n.
const SQL_LONGVARCHAR   = Int16( -1) # Variable length character data. Maximum length is data source–dependent.
const SQL_WCHAR         = Int16( -8) # Unicode character string of fixed string length n
const SQL_WVARCHAR      = Int16( -9) # Unicode variable-length character string with a maximum string length n
const SQL_WLONGVARCHAR  = Int16(-10) # Unicode variable-length character data. Maximum length is data source–dependent
const SQL_DECIMAL       = Int16(  3)
const SQL_NUMERIC       = Int16(  2)
const SQL_SMALLINT      = Int16(  5) # Exact numeric value with precision 5 and scale 0 (signed: –32,768 <= n <= 32,767, unsigned: 0 <= n <= 65,535)
const SQL_INTEGER       = Int16(  4) # Exact numeric value with precision 10 and scale 0 (signed: –2[31] <= n <= 2[31] – 1, unsigned: 0 <= n <= 2[32] – 1)
const SQL_REAL          = Int16(  7) # Signed, approximate, numeric value with a binary precision 24 (zero or absolute value 10[–38] to 10[38]).
const SQL_FLOAT         = Int16(  6) # Signed, approximate, numeric value with a binary precision of at least p. (The maximum precision is driver-defined.)
const SQL_DOUBLE        = Int16(  8) # Signed, approximate, numeric value with a binary precision 53 (zero or absolute value 10[–308] to 10[308]).
const SQL_BIT           = Int16( -7) # Single bit binary data.
const SQL_TINYINT       = Int16( -6) # Exact numeric value with precision 3 and scale 0 (signed: –128 <= n <= 127, unsigned: 0 <= n <= 255)
const SQL_BIGINT        = Int16( -5) # Exact numeric value with precision 19 (if signed) or 20 (if unsigned) and scale 0 (signed: –2[63] <= n <= 2[63] – 1, unsigned: 0 <= n <= 2[64] – 1)
const SQL_BINARY        = Int16( -2) # Binary data of fixed length n.
const SQL_VARBINARY     = Int16( -3) # Variable length binary data of maximum length n. The maximum is set by the user.
const SQL_LONGVARBINARY = Int16( -4) # Variable length binary data. Maximum length is data source–dependent.
const SQL_TYPE_DATE     = Int16( 91) # Year, month, and day fields, conforming to the rules of the Gregorian calendar.
const SQL_TYPE_TIME     = Int16( 92) # Hour, minute, and second fields, with valid values for hours of 00 to 23,
                                     # valid values for minutes of 00 to 59, and valid values for seconds of 00 to 61. Precision p indicates the seconds precision.
const SQL_TYPE_TIMESTAMP = Int16( 93) # Year, month, day, hour, minute, and second fields, with valid values as defined for the DATE and TIME data types.
# SQL Server specific
const SQL_SS_TIME2       = Int16(-154)
const SQL_SS_TIMESTAMPOFFSET = Int16(-155)

#const SQL_INTERVAL_MONTH            = Int16(102)
#const SQL_INTERVAL_YEAR             = Int16(101)
#const SQL_INTERVAL_YEAR_TO_MONTH    = Int16(107)
#const SQL_INTERVAL_DAY              = Int16(103)
#const SQL_INTERVAL_HOUR             = Int16(104)
#const SQL_INTERVAL_MINUTE           = Int16(105)
#const SQL_INTERVAL_SECOND           = Int16(106)
#const SQL_INTERVAL_DAY_TO_HOUR      = Int16(108)
#const SQL_INTERVAL_DAY_TO_MINUTE    = Int16(109)
#const SQL_INTERVAL_DAY_TO_SECOND    = Int16(110)
#const SQL_INTERVAL_HOUR_TO_MINUTE   = Int16(111)
#const SQL_INTERVAL_HOUR_TO_SECOND   = Int16(112)
#const SQL_INTERVAL_MINUTE_TO_SECOND = Int16(113)
const SQL_GUID                      = Int16(-11) # Fixed length GUID.

# C Data Types
const SQL_C_CHAR      = Int16(  1)
const SQL_C_WCHAR     = Int16( -8)
const SQL_C_DOUBLE    = Int16(  8)
const SQL_C_SHORT     = Int16(  5)
const SQL_C_LONG      = Int16(  4)
const SQL_C_FLOAT     = Int16(  7)
const SQL_C_NUMERIC   = Int16(  2)
const SQL_C_BIT       = Int16( -7)
const SQL_C_TINYINT   = Int16( -6)
const SQL_C_BIGINT    = Int16(-27)
const SQL_C_BINARY    = Int16( -2)
const SQL_C_TYPE_DATE = Int16( 91)
const SQL_C_TYPE_TIME = Int16( 92)
const SQL_C_TYPE_TIMESTAMP = Int16( 93)

#const SQL_C_INTERVAL_MONTH            = Int16(102)
#const SQL_C_INTERVAL_YEAR             = Int16(101)
#const SQL_C_INTERVAL_YEAR_TO_MONTH    = Int16(107)
#const SQL_C_INTERVAL_DAY              = Int16(103)
#const SQL_C_INTERVAL_HOUR             = Int16(104)
#const SQL_C_INTERVAL_MINUTE           = Int16(105)
#const SQL_C_INTERVAL_SECOND           = Int16(106)
#const SQL_C_INTERVAL_DAY_TO_HOUR      = Int16(108)
#const SQL_C_INTERVAL_DAY_TO_MINUTE    = Int16(109)
#const SQL_C_INTERVAL_DAY_TO_SECOND    = Int16(110)
#const SQL_C_INTERVAL_HOUR_TO_MINUTE   = Int16(111)
#const SQL_C_INTERVAL_HOUR_TO_SECOND   = Int16(112)
#const SQL_C_INTERVAL_MINUTE_TO_SECOND = Int16(113)
const SQL_C_GUID                      = Int16(-11)

# Julia mapping C structs
struct SQLDate <: Dates.AbstractTime
    year::Int16
    month::Int16
    day::Int16
end

Base.show(io::IO,x::SQLDate) = print(io,"$(x.year)-$(lpad(x.month,2,'0'))-$(lpad(x.day,2,'0'))")
SQLDate(x::Dates.Date) = SQLDate(Dates.yearmonthday(x)...)
SQLDate() = SQLDate(0,0,0)
import Base: ==
==(x::SQLDate, y::Dates.Date) = x.year == Dates.year(y) && x.month == Dates.month(y) && x.day == Dates.day(y)
==(y::Dates.Date, x::SQLDate) = x.year == Dates.year(y) && x.month == Dates.month(y) && x.day == Dates.day(y)
Dates.Date(x::SQLDate) = Dates.Date(x.year, x.month, x.day)

struct SQLTime <: Dates.AbstractTime
    hour::Int16
    minute::Int16
    second::Int16
end

Base.show(io::IO,x::SQLTime) = print(io,"$(lpad(x.hour,2,'0')):$(lpad(x.minute,2,'0')):$(lpad(x.second,2,'0'))")
SQLTime() = SQLTime(0,0,0)

struct SQLTimestamp <: Dates.AbstractTime
    year::Int16
    month::Int16
    day::Int16
    hour::Int16
    minute::Int16
    second::Int16
    fraction::Int32 #nanoseconds
end

Base.show(io::IO,x::SQLTimestamp) = print(io,"$(x.year)-$(lpad(x.month,2,'0'))-$(lpad(x.day,2,'0'))T$(lpad(x.hour,2,'0')):$(lpad(x.minute,2,'0')):$(lpad(x.second,2,'0'))$(x.fraction == 0 ? "" : strip(Printf.@sprintf("%.9f",x.fraction/1e+9),'0'))")
function SQLTimestamp(x::Dates.DateTime)
    y, m, d = Dates.yearmonthday(x)
    h, mm, s = Dates.hour(x), Dates.minute(x), Dates.second(x)
    frac = Dates.millisecond(x) * 1_000_000
    return SQLTimestamp(y, m, d, h, mm, s, frac)
end
SQLTimestamp() = SQLTimestamp(0,0,0,0,0,0,0)
==(x::SQLTimestamp, y::Dates.DateTime) = x.year == Dates.year(y) && x.month == Dates.month(y) && x.day == Dates.day(y) &&
                                   x.hour == Dates.hour(y) && x.minute == Dates.minute(y) && x.second == Dates.second(y)
==(y::Dates.DateTime, x::SQLTimestamp) = x.year == Dates.year(y) && x.month == Dates.month(y) && x.day == Dates.day(y) &&
                               x.hour == Dates.hour(y) && x.minute == Dates.minute(y) && x.second == Dates.second(y)
Dates.DateTime(x::SQLTimestamp) = Dates.DateTime(x.year, x.month, x.day, x.hour, x.minute, x.second, x.fraction ÷ 1_000_000)

const SQL_MAX_NUMERIC_LEN = 16
struct SQLNumeric
    precision::SQLCHAR
    scale::SQLSCHAR
    sign::SQLCHAR
    val::NTuple{SQL_MAX_NUMERIC_LEN,SQLCHAR}
end

Base.show(io::IO,x::SQLNumeric) = print(io,"SQLNumeric($(x.sign == 1 ? '+' : '-') precision: $(x.precision) scale: $(x.scale) val: $(x.val))")
SQLNumeric() = SQLNumeric(0,0,0,(0,))

# typedef struct  tagSQLGUID
# {
#     DWORD Data1;
#     WORD Data2;
#     WORD Data3;
#     BYTE Data4[ 8 ];
# } SQLGUID;
struct SQLGUID
    Data1::Cuint
    Data2::Cushort
    Data3::Cushort
    Data4::NTuple{8,Cuchar}
end

# for representing SQL LONG types
struct Long{T} end

"""
Dict for mapping SQL types to C types.
When executing an SQL query that returns results, the DBMS will return the SQL types of the resultset;
The application then tells the ODBC Driver Manager how to actually return the data
(by specifying the equivalent C type or something else entirely).
"""
const SQL2C = Dict(
    SQL_CHAR           => SQL_C_CHAR,
    SQL_VARCHAR        => SQL_C_CHAR,
    SQL_LONGVARCHAR    => SQL_C_CHAR,
    SQL_WCHAR          => SQL_C_WCHAR,
    SQL_WVARCHAR       => SQL_C_WCHAR,
    SQL_WLONGVARCHAR   => SQL_C_WCHAR,
    SQL_DECIMAL        => SQL_C_CHAR,
    SQL_NUMERIC        => SQL_C_CHAR,
    SQL_SMALLINT       => SQL_C_SHORT,
    SQL_INTEGER        => SQL_C_LONG,
    SQL_REAL           => SQL_C_FLOAT,
    SQL_FLOAT          => SQL_C_DOUBLE,
    SQL_DOUBLE         => SQL_C_DOUBLE,
    SQL_BIT            => SQL_C_BIT,
    SQL_TINYINT        => SQL_C_TINYINT,
    SQL_BIGINT         => SQL_C_BIGINT,
    SQL_BINARY         => SQL_C_BINARY,
    SQL_VARBINARY      => SQL_C_BINARY,
    SQL_LONGVARBINARY  => SQL_C_BINARY,
    SQL_TYPE_DATE      => SQL_C_TYPE_DATE,
    SQL_TYPE_TIME      => SQL_C_TYPE_TIME,
    SQL_TYPE_TIMESTAMP => SQL_C_TYPE_TIMESTAMP,
    SQL_SS_TIME2       => SQL_C_TYPE_TIME,
    SQL_SS_TIMESTAMPOFFSET => SQL_C_TYPE_TIMESTAMP,
    SQL_GUID           => SQL_C_GUID)

const julia2SQL = Dict(
    String                => SQL_CHAR,
    WeakRefString{UInt8}  => SQL_CHAR,
    WeakRefString{UInt16} => SQL_WCHAR,
    WeakRefString{UInt32} => SQL_WCHAR,
    Dec64                 => SQL_DOUBLE,
    Float16               => SQL_FLOAT,
    Float32               => SQL_FLOAT,
    Float64               => SQL_DOUBLE,
    Int8                  => SQL_TINYINT,
    Int16                 => SQL_SMALLINT,
    Int32                 => SQL_INTEGER,
    Int64                 => SQL_BIGINT,
    Bool                  => SQL_BIT,
    Vector{UInt8}         => SQL_BINARY,
    Dates.Date            => SQL_TYPE_DATE,
    Dates.DateTime        => SQL_TYPE_TIMESTAMP,
    SQLDate               => SQL_TYPE_DATE,
    SQLTime               => SQL_TYPE_TIME,
    SQLTimestamp          => SQL_TYPE_TIMESTAMP
)

const julia2C = Dict(
    String                => SQL_C_CHAR,
    WeakRefString{UInt8}  => SQL_C_CHAR,
    WeakRefString{UInt16} => SQL_C_WCHAR,
    WeakRefString{UInt32} => SQL_C_WCHAR,
    Dec64                 => SQL_C_DOUBLE,
    Float16               => SQL_C_FLOAT,
    Float32               => SQL_C_FLOAT,
    Float64               => SQL_C_DOUBLE,
    Int8                  => SQL_C_TINYINT,
    Int16                 => SQL_C_SHORT,
    Int32                 => SQL_C_LONG,
    Int64                 => SQL_C_BIGINT,
    Bool                  => SQL_C_BIT,
    Vector{UInt8}         => SQL_C_BINARY,
    Dates.Date            => SQL_C_TYPE_DATE,
    Dates.DateTime        => SQL_C_TYPE_TIMESTAMP,
    SQLDate               => SQL_C_TYPE_DATE,
    SQLTime               => SQL_C_TYPE_TIME,
    SQLTimestamp          => SQL_C_TYPE_TIMESTAMP
)

"""
maps SQL types from the ODBC manager to Julia types;
in particular, it returns a Tuple{A,B,Bool}, where `A` is the Julia type
used to allocate and return data from the ODBC manager, while `B`
represents the final column type of the data; conversion from `A` to `B` happens
through the `cast(::Type{B})` function in utils.jl.
The 3rd `Bool` value indicates whether the column is a LONGTEXT or LONGBINARY SQL types, since these
tend to require special result-handling rules.
"""
const SQL2Julia = Dict(
    SQL_CHAR           => (SQLCHAR, Union{String, Missing}, false),
    SQL_VARCHAR        => (SQLCHAR, Union{String, Missing}, false),
    SQL_LONGVARCHAR    => (SQLCHAR, Union{String, Missing}, true),
    SQL_WCHAR          => (SQLWCHAR, Union{String, Missing}, false),
    SQL_WVARCHAR       => (SQLWCHAR, Union{String, Missing}, false),
    SQL_WLONGVARCHAR   => (SQLWCHAR, Union{String, Missing}, true),
    SQL_DECIMAL        => (SQLCHAR, Union{Dec64, Missing}, false),
    SQL_NUMERIC        => (SQLCHAR, Union{Dec64, Missing}, false),
    SQL_SMALLINT       => (SQLSMALLINT, Union{SQLSMALLINT, Missing}, false),
    SQL_INTEGER        => (SQLINTEGER,  Union{SQLINTEGER, Missing}, false),
    SQL_REAL           => (SQLREAL,   Union{SQLREAL, Missing}, false),
    SQL_FLOAT          => (SQLFLOAT,  Union{SQLFLOAT, Missing}, false),
    SQL_DOUBLE         => (SQLDOUBLE, Union{SQLDOUBLE, Missing}, false),
    SQL_BIT            => (Int8,  Union{Int8, Missing}, false),
    SQL_TINYINT        => (Int8,  Union{Int8, Missing}, false),
    SQL_BIGINT         => (Int64, Union{Int64, Missing}, false),
    SQL_BINARY         => (UInt8, Union{Vector{UInt8}, Missing}, false),
    SQL_VARBINARY      => (UInt8, Union{Vector{UInt8}, Missing}, false),
    SQL_LONGVARBINARY  => (UInt8, Union{Vector{UInt8}, Missing}, true),
    SQL_TYPE_DATE      => (SQLDate, Union{SQLDate, Missing}, false),
    SQL_TYPE_TIME      => (SQLTime, Union{SQLTime, Missing}, false),
    SQL_TYPE_TIMESTAMP => (SQLTimestamp, Union{SQLTimestamp, Missing}, false),
    SQL_SS_TIME2       => (SQLTime, Union{SQLTime, Missing}, false),
    SQL_SS_TIMESTAMPOFFSET => (SQLTimestamp, Union{SQLTimestamp, Missing}, false),
    SQL_GUID           => (SQLGUID, Union{SQLGUID, Missing}, false))

"Convenience mapping of SQL types to their string representation"
const SQL_TYPES = Dict(
      1 => "SQL_CHAR",
     12 => "SQL_VARCHAR",
     -1 => "SQL_LONGVARCHAR",
     -8 => "SQL_WCHAR",
     -9 => "SQL_WVARCHAR",
    -10 => "SQL_WLONGVARCHAR",
      3 => "SQL_DECIMAL",
      2 => "SQL_NUMERIC",
      5 => "SQL_SMALLINT",
      4 => "SQL_INTEGER",
      7 => "SQL_REAL",
      6 => "SQL_FLOAT",
      8 => "SQL_DOUBLE",
     -7 => "SQL_BIT",
     -6 => "SQL_TINYINT",
     -5 => "SQL_BIGINT",
     -2 => "SQL_BINARY",
     -3 => "SQL_VARBINARY",
     -4 => "SQL_LONGVARBINARY",
     91 => "SQL_TYPE_DATE",
     92 => "SQL_TYPE_TIME",
     93 => "SQL_TYPE_TIMESTAMP",
   -154 => "SQL_SS_TIME2",
   -155 => "SQL_SS_TIMESTAMPOFFSET",
    102 => "SQL_INTERVAL_MONTH",
    101 => "SQL_INTERVAL_YEAR",
    107 => "SQL_INTERVAL_YEAR_TO_MONTH",
    103 => "SQL_INTERVAL_DAY",
    104 => "SQL_INTERVAL_HOUR",
    105 => "SQL_INTERVAL_MINUTE",
    106 => "SQL_INTERVAL_SECOND",
    108 => "SQL_INTERVAL_DAY_TO_HOUR",
    109 => "SQL_INTERVAL_DAY_TO_MINUTE",
    110 => "SQL_INTERVAL_DAY_TO_SECOND",
    111 => "SQL_INTERVAL_HOUR_TO_MINUTE",
    112 => "SQL_INTERVAL_HOUR_TO_SECOND",
    113 => "SQL_INTERVAL_MINUTE_TO_SECOND",
    -11 => "SQL_GUID")

"Convenience mapping of SQL types to their C-type equivalent as a string"
const C_TYPES = Dict(
    SQL_CHAR           => "SQL_C_CHAR",
    SQL_VARCHAR        => "SQL_C_CHAR",
    SQL_LONGVARCHAR    => "SQL_C_CHAR",
    SQL_WCHAR          => "SQL_C_WCHAR",
    SQL_WVARCHAR       => "SQL_C_WCHAR",
    SQL_WLONGVARCHAR   => "SQL_C_WCHAR",
    SQL_DECIMAL        => "SQL_C_NUMERIC",
    SQL_NUMERIC        => "SQL_C_NUMERIC",
    SQL_SMALLINT       => "SQL_C_SHORT",
    SQL_INTEGER        => "SQL_C_LONG",
    SQL_REAL           => "SQL_C_FLOAT",
    SQL_FLOAT          => "SQL_C_DOUBLE",
    SQL_DOUBLE         => "SQL_C_DOUBLE",
    SQL_BIT            => "SQL_C_BIT",
    SQL_TINYINT        => "SQL_C_TINYINT",
    SQL_BIGINT         => "SQL_C_BIGINT",
    SQL_BINARY         => "SQL_C_BINARY",
    SQL_VARBINARY      => "SQL_C_BINARY",
    SQL_LONGVARBINARY  => "SQL_C_BINARY",
    SQL_TYPE_DATE      => "SQL_C_TYPE_DATE",
    SQL_TYPE_TIME      => "SQL_C_TYPE_TIME",
    SQL_TYPE_TIMESTAMP => "SQL_C_TYPE_TIMESTAMP",
    SQL_C_GUID         => "SQL_C_GUID")


text:
use_cuda = false

using Mocha

srand(12345678)

############################################################
# Prepare Random Data
############################################################
N = 5    # works with arbitrary minibatch size as long as
         # N == batch_size in MemoryDataLayer so it cycles through
         # and gets the same data during forward()
M = 10
P = 4

X = rand(M, N)
W = rand(M, P)
B = rand(P, 1)

Y = (W'*X .+ B)
Y = Y + 0.01*randn(size(Y))

############################################################
# Define network
############################################################
backend = CPUBackend()
init(backend)

data_layer = MemoryDataLayer(batch_size=N, data=Array[X,Y])

w1 = InnerProductLayer(neuron=Neurons.Sigmoid(), name="ip1",output_dim=20, tops=[:a], bottoms=[:data])
w2 = InnerProductLayer(neuron=Neurons.Identity(), name="ip2",output_dim=4, tops=[:b], bottoms=[:a])
loss_layer = SquareLossLayer(name="loss", bottoms=[:b, :label] )

net = Net("TEST", backend, [w1,w2, loss_layer, data_layer])

# epsilon:     milage may vary 1e-4 - 1e-8
# digit:       compare this many digits to check for 'identity'
# visualize:   prints out correct and failed positions of gradients
# return value: false if any gradients fails, true if all passes
test_gradients(net, epsilon=1e-8, digit=6, visual=true )

shutdown(backend)

ext: jl


text:
abstract type AbstractConsolidation end
abstract type AbstractConsolidator end

# No consolidation
# ------------------------------------------------------------
"""
    NoRegularization

Empty functor object for running the L-shaped algorithm without consolidation.

"""
struct NoConsolidation <: AbstractConsolidation end

function consolidate!(::AbstractLShapedSolver, ::NoConsolidation)
    return nothing
end

function add_cut!(::AbstractLShapedSolver, ::NoConsolidation, ::AbstractHyperPlane)
    return nothing
end

function add_cut!(::AbstractLShapedSolver, ::NoConsolidation, ::Integer, ::AbstractHyperPlane)
    return nothing
end
# Consolidation
# ------------------------------------------------------------
"""
    Consolidation

Functor object for using consolidation in an L-shaped algorithm. Create by supplying a [`Consolidate`](@ref) object through `consolidate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Algorithm parameters
- `tresh::T` = 0.95: Relative amount of redundant cuts in a former iteration required to consider the iteration redundant
- `at::Int = 5.0`: Number of times an iteration can be redundant before consolidation is triggered
- `rebuild::Function = at_tolerance()`: Function deciding when the master model should be rebuilt according to performed consolidations
...
"""
struct Consolidation{T <: AbstractFloat} <: AbstractConsolidation
    cuts::Vector{Vector{AnySparseOptimalityCut{T}}}
    feasibility_cuts::Vector{Vector{SparseFeasibilityCut{T}}}
    consolidated::Vector{Bool}
    redundance_count::Vector{Int}
    redundance_treshold::T
    consolidation_trigger::Int
    rebuild::Function

   function Consolidation(::Type{T}, tresh::T, at::Int, rebuild::Function) where T <: AbstractFloat
        cuts = Vector{Vector{AnySparseOptimalityCut{T}}}()
        push!(cuts, Vector{AnySparseOptimalityCut{T}}())
        feasibility_cuts = Vector{Vector{SparseFeasibilityCut{T}}}()
        push!(feasibility_cuts, Vector{SparseFeasibilityCut{T}}())
        return new{T}(cuts, feasibility_cuts, [false], [0], tresh, at, rebuild)
    end
end

function consolidate!(lshaped::AbstractLShapedSolver, consolidation::Consolidation{T}) where T <: AbstractFloat
    for i in findall(map(!,consolidation.consolidated))
        nc = sum(nsubproblems.(consolidation.cuts[i]))
        nc == 0 && continue
        count = 0
        for cut in consolidation.cuts[i]
            # Check for redundance
            if !active(lshaped, cut)
                count += nsubproblems(cut)
            end
        end
        if count/nc >= consolidation.redundance_treshold
            consolidation.redundance_count[i] += 1
        end
        if consolidation.redundance_count[i] == consolidation.consolidation_trigger
            consolidated_cut = aggregate(consolidation.cuts[i])
            empty!(consolidation.cuts[i])
            push!(consolidation.cuts[i], consolidated_cut)
            consolidation.consolidated[i] = true
        end
    end
    for i in findall(c->length(c)>0, consolidation.feasibility_cuts)
        inactive = Vector{Int}()
        for (j,cut) in enumerate(consolidation.feasibility_cuts[i])
            # Check for redundance
            if !active(lshaped, cut)
                push!(inactive, j)
            end
        end
        deleteat!(consolidation.feasibility_cuts[i], inactive)
    end
    if minimum_requirements(lshaped,consolidation) && consolidation.rebuild(lshaped, consolidation)
        # Rebuild sparingly
        rebuild_master!(lshaped, consolidation)
        lshaped.data.consolidations += 1
    end
    # Prepare memory for next iteration
    push!(consolidation.cuts, Vector{AnySparseOptimalityCut{T}}())
    push!(consolidation.feasibility_cuts, Vector{SparseFeasibilityCut{T}}())
    push!(consolidation.consolidated, false)
    push!(consolidation.redundance_count, 0)
    return nothing
end

function add_cut!(::AbstractLShapedSolver, ::Consolidation, ::AbstractHyperPlane)
    return nothing
end

function add_cut!(lshaped::AbstractLShapedSolver, consolidation::Consolidation{T}, cut::HyperPlane{FeasibilityCut}) where T <: AbstractFloat
    push!(consolidation.feasibility_cuts[timestamp(lshaped)], cut)
    return nothing
end

function add_cut!(lshaped::AbstractLShapedSolver, consolidation::Consolidation{T}, cut::AnySparseOptimalityCut) where T <: AbstractFloat
    push!(consolidation.cuts[timestamp(lshaped)], cut)
    return nothing
end

function rebuild_master!(lshaped::AbstractLShapedSolver, consolidation::Consolidation)
    ncuts = ncutconstraints(lshaped)
    cut_indices = first_stage_nconstraints(lshaped.stochasticprogram)+1:MPB.numconstr(lshaped.mastersolver.lqmodel)
    MPB.delconstrs!(lshaped.mastersolver.lqmodel, collect(cut_indices))
    if lshaped.regularization isa LevelSet
        # Reset projection problem if using level sets
        lv = lshaped.regularization
        MPB.delconstrs!(lv.projectionsolver.lqmodel, [lv.data.levelindex])
        cut_indices = first_stage_nconstraints(lshaped.stochasticprogram)+1:MPB.numconstr(lv.projectionsolver.lqmodel)
        MPB.delconstrs!(lv.projectionsolver.lqmodel, collect(cut_indices))
        lv.data.levelindex = -1
        if :index ∈ fieldnames(typeof(lv.penaltyterm))
            lv.penaltyterm.index = -1
        end
    elseif lshaped.regularization isa RegularizedDecomposition
        rd = lshaped.regularization
        if :index ∈ fieldnames(typeof(rd.penaltyterm))
            rd.penaltyterm.index = -1
        end
    end
    lshaped.data.ncuts -= ncuts
    readd_cuts!(lshaped, consolidation)
    return nothing
end

noptimalitycuts(consolidation::Consolidation) = sum(length.(consolidation.cuts))
nfeasibilitycuts(consolidation::Consolidation) = sum(length.(consolidation.feasibility_cuts))
function ncutconstraints(consolidation::Consolidation)
    return noptimalitycuts(consolidation) + nfeasibilitycuts(consolidation)
end

function ncutconstraints(lshaped::AbstractLShapedSolver)
    return MPB.numconstr(lshaped.mastersolver.lqmodel)-first_stage_nconstraints(lshaped.stochasticprogram)
end
# Rebuild functions
# ------------------------------------------------------------
function minimum_requirements(lshaped::AbstractLShapedSolver, consolidation::Consolidation)
    return (noptimalitycuts(consolidation) > 0 &&
        count(consolidation.consolidated) > 0) ||
        (nfeasibilitycuts(consolidation) > 0 &&
         nfeasibilitycuts(consolidation)/ncutconstraints(lshaped) <= 0.1)
end
"""
    at_tolerance(τ = 0.4, miniter = 0)

Rebuild master when at least nconsolidations*`miniter` iterations has passed and the ratio of number of cuts in the consolidated collection and the number of cuts in the master model has decreased below `τ`.

"""
at_tolerance() = at_tolerance(0.4, 0)
function at_tolerance(τ, miniter)
    return (lshaped, consolidation) -> begin
        return lshaped.data.iterations >= (lshaped.data.consolidations+1)*miniter && ncutconstraints(consolidation)/ncutconstraints(lshaped) <= τ
    end
end

# API
# ------------------------------------------------------------
"""
    DontConsolidate

Factory object for [`NoConsolidation`](@ref). Passed by default to `consolidate` in the `LShapedSolver` factory function.

"""
struct DontConsolidate <: AbstractConsolidator end

function (::DontConsolidate)(::Type{T} where T <: AbstractFloat)
    return NoConsolidation()
end

"""
    Consolidate

Factory object for [`Consolidation`](@ref). Pass to `consolidate` in the `LShapedSolver` factory function. See ?Consolidation for parameter descriptions.

"""
struct Consolidate <: AbstractConsolidator
    redundance_treshold::Float64
    consolidation_trigger::Int
    rebuild::Function
end
Consolidate(; tresh::Float64 = 0.95, at::Int = 5, rebuild = at_tolerance()) = Consolidate(tresh, at, rebuild)

function (consolidator::Consolidate)(::Type{T}) where T <: AbstractFloat
    return Consolidation(T,
                         consolidator.redundance_treshold,
                         consolidator.consolidation_trigger,
                         consolidator.rebuild)
end


text:
module Gap

using JuMP
using Cbc

include("Computation.jl")
using .Computation

include("Format.jl")
using .Format

include("FC.jl")
using .FC

include("Half.jl")
using .Half

export gap, vgap

# Determines upper bound alpha with the Gap Method
# Author: Jason Liu
function gap(m::Int64, s::Int64; output::Int64 = 2)
    output > 0 && printHeader(center("GAP METHOD"))

    if m <= 0 || s <= 0
        if output > 0
            printf("Gap Method does not apply", line = true)
            printEnd()
        end
        return 1 // 1
    elseif m < s
        alpha = m // s * gap(s, m, output = 0)
        if output > 1
            printfT(
                "Duality Theorem",
                "muffins($m,$s) = $m/$s · muffins($s,$m)",
                "Bounding muffins($s,$m) instead",
            )
            gap(s, m, output = output)
        elseif output > 0
            printfT("Gap Method", "Upper bound α of muffins($m,$s) is $(formatFrac(alpha))")
            printEnd()
        end
        return alpha
    elseif m % s == 0
        if output > 0
            printfT("Gap Method", "Since $m % $s = 0, muffins($m,$s) = 1")
            printEnd()
        end
        return 1 // 1
    end

    # Compute alpha candidates
    alphas = []
    upper = min(fc(m, s, output = 0), half(m, s, output = 0))
    for b = 1:m
        for n = Int64(ceil(b * s / 3)):Int64(floor(b * s * upper))
            append!(alphas, n // (b * s))
        end
    end
    append!(alphas, 1 // 3)
    alphas = sort(unique(alphas))

    # Test alpha candidates with binary search
    potAlphas = []
    (left, right) = (1, length(alphas))
    while left <= right
        mid = Int64(floor((left + right) / 2))
        append!(potAlphas, alphas[mid])
        if vgap(m, s, alphas[mid], output = 0)
            right = mid - 1
        else
            left = mid + 1
        end
    end

    for alpha in reverse(potAlphas)
        if vgap(m, s, alpha, output = 0)
            vgap(m, s, alpha, output = output)
            return alpha
        end
    end

    if output > 0
        printf("All α candidates were inconclusive, Gap Method inconclusive", line = true)
        printEnd()
    end
    1 // 1
end

# Helper function for gap -- verifies whether gap(m, s, alpha) is conclusive
function vgap(m::Int64, s::Int64, alpha::Rational{Int64}; output::Int64 = 2)
    if !(m > s > 0) || m % s == 0
        output > 0 && printf("VGap does not apply", line = true)
        false
    elseif alpha < 1 / 3
        output > 0 && printfT("No piece size < 1/3", "For m ≥ s, α must be ≥ 1/3")
        false
    elseif alpha > 1
        output > 0 && printfT("No piece size > 1", "α must be ≤ 1")
        false
    else
        (V, W, sV, sW) = sv(m, s)
        numV = (V)sV
        numW = (W)sW

        # Initialize Gap Method proof
        if output > 1
            # Define and format variables for proof
            alphaF = formatFrac(alpha)
            size = formatFrac(m // s)
            a = formatFrac(m // (V + 1)s)
            b = max(1 - m // (W - 1)s, 0)
            b1 = formatFrac(1 - b)
            b = formatFrac(b)

            # Establish assumptions and premises
            printHeader("CLAIM")
            printfT(
                "Claim",
                "For all ($m, $s) procedures where $s students each receive $size muffins, the smallest piece size is bounded above by α = $alphaF",
                "Proven by contradicting the assumption",
            )

            printHeader("ASSUMPTION")
            printfT("Assumption", "Assume muffins($m,$s) > α ≥ 1/3")
            printfT(
                "Property of Buddies",
                "Since muffins($m,$s) > 1/3, each muffin must be cut into 2 shs, totaling $(2m) shs",
                "",
                "Since each muffin is cut into 2 shs that are buddies, each sh of size x must imply the existence of a sh of size 1-x",
            )

            # Begin casework
            printHeader("CASEWORK")

            printfT(
                "Case 1",
                "Alice has ≥ $(V+1) shs",
                "One of them is ≤ ($size × 1/$(V+1)) = $a",
                "",
                "Contradicts assumption if α ≥ $a",
            )

            if W - 1 <= 1 && m / s > 1
                printfT(
                    "Case 2",
                    "Bob cannot have ≤ $(W-1) shs since $size > 1, so this case is impossible",
                )
            else
                printfT(
                    "Case 2",
                    "Bob has ≤ $(W-1) shs",
                    "One of them is ≥ ($size × 1/$(W-1)) = $b1",
                    "Its buddy is ≤ (1 - $b1) = $b",
                    "",
                    "Contradicts assumption if α ≥ $b",
                )
            end
        end

        # Check if V-Conjecture applies
        if m // s * 1 // (V + 1) > alpha || 1 - m // s * 1 // (W - 1) > alpha
            output > 0 && printf("V-Conjecture does not apply, VGap failed", line = true)
            return false
        end

        ((_, x), (y, _)) = findend(m, s, alpha, V)
        xOriginal = m // s - alpha * (V - 1)
        yOriginal = m // s - (1 - alpha) * (V - 2)

        # Check if FindEnd works
        if x != xOriginal && y != yOriginal
            output > 0 && printf("FindEnd inconclusive, VGap failed", line = true)
            return false
        end

        # Define and format variables for proof
        cd = lcm(s, denominator(x), denominator(y), 2)
        alphaF = formatFrac(alpha, cd)
        alpha1 = formatFrac(1 - alpha, cd)
        size = formatFrac(m // s)
        a = formatFrac(m // s - x, cd)
        b = formatFrac(m // s - y, cd)
        xF = formatFrac(x, cd)
        x1 = formatFrac(1 - x, cd)
        yF = formatFrac(y, cd)
        y1 = formatFrac(1 - y, cd)
        diff = abs(numV - numW)

        # Continue proof
        if output > 1
            # Continue casework
            printfT(
                "Note",
                "The remaining cases deal with everyone having either $W or $V shs, so:",
                "",
                "$(W)·s_$W + $(V)·s_$V = $(2m)  (total shs)",
                "s_$W + s_$V = $s  (total students)",
                "where s_N = # of students with N shs",
                "",
                "The solution to the system is:",
                "s_$W = $sW, s_$V = $sV",
                "So there are $numW $W-shs and $numV $V-shs",
            )

            if x == xOriginal
                printfT(
                    "Case 3",
                    "Alice has a $V-sh ≥ $xF",
                    "Her other $(V-1) $V-shs sum to ≤ ($size - $xF) = $a",
                    "One of them is ≤ ($a × 1/$(V-1)) = $alphaF",
                    "",
                    "Contradicts assumption if α ≥ $alphaF",
                )
            else
                printfT("Case 3", "FindEnd did not produce a conclusive bound for $V-shs")
            end

            if y == yOriginal
                printfT(
                    "Case 4",
                    "Bob has a $W-sh ≤ $yF",
                    "His other $(W-1) $W-shs sum to ≥ ($size - $yF) = $b",
                    "One of them is ≥ ($b × 1/$(W-1)) = $alpha1",
                    "Its buddy is ≤ (1 - $alpha1) = $alphaF",
                    "",
                    "Contradicts assumption if α ≥ $alphaF",
                )
            else
                printfT("Case 4", "FindEnd did not produce a conclusive bound for $W-shs")
            end

            printHeader("CASE 5: INTERVAL ANALYSIS")
            println()
            printf("The following intervals capture the negation of the previous cases:")
            if alpha < x <= y < 1 - alpha
                println(
                    "\n",
                    interval(
                        ["(", alphaF],
                        [")[", xF],
                        ["](", yF],
                        [")", alpha1],
                        labels = ["$numV $V-shs", "0", "$numW $W-shs"],
                    ),
                )
                printLine()
            else
                if xF != alphaF
                    println(
                        "\n",
                        interval(
                            ["(", alphaF],
                            [")[", xF],
                            [")", alpha1],
                            labels = ["$numV $V-shs", "0 $V-shs"],
                        ),
                    )
                else
                    println("\n", "No negation interval for $V-shs")
                end

                if yF != alpha1
                    println(
                        "\n",
                        interval(
                            ["(", alphaF],
                            ["](", yF],
                            [")", alpha1],
                            labels = ["0 $W-shs", "$numW $W-shs"],
                        ),
                    )
                else
                    println("\n", "No negation interval for $W-shs")
                end

                printfT(
                    "Case 5",
                    "The Gap Method is inconclusive on these intervals, VGap failed",
                )
            end
        end

        # Fail if interval inconclusive
        if x > y || x == alpha || y == 1 - alpha
            output == 1 &&
                printf("Could not generate disjoint intervals, VGap failed", line = true)
            return false
        end

        # Define and format variables for proof
        (vMin, vMax, sMin, sMax) = (V, W, sV, sW)
        (i, j, k, l) = (xF, yF, y1, x1)
        (rngMin, rngMax, rngS, rngL) = ((alphaF, xF), (x1, alpha1), (yF, y1), (x1, alpha1))
        if numV > numW
            (vMin, vMax, sMin, sMax) = (W, V, sW, sV)
            (i, j, k, l) = (y1, x1, xF, yF)
            (rngMin, rngMax, rngS, rngL) =
                ((yF, alpha1), (alphaF, y1), (alphaF, y1), (x1, xF))
        end
        numMin = (vMin)sMin

        # Conclude interval case with case by case analysis
        if output > 1
            printfT(
                "Property of Buddies",
                "Because 0 shs lie in [$xF,$yF], 0 shs must also lie in [$y1,$x1]",
                "",
                "Similary, the $numMin $vMin-shs in ($(rngMin[1]),$(rngMin[2])) implies $numMin $vMax-shs in ($(rngMax[1]),$(rngMax[2]))",
            )
            println()
            printf("The following intervals capture the previous statements:")
            println(
                "\n",
                interval(
                    ["(", alphaF],
                    [")[", i],
                    ["](", j],
                    [")[", k],
                    ["](", l],
                    [")", alpha1],
                    labels = [
                        "$numMin $V-shs",
                        "0",
                        "$diff $vMax-shs",
                        "0",
                        "$numMin $W-shs",
                    ],
                ),
            )
            printLine()

            printfT(
                "Note",
                "Let the $vMax-shs that lie in ($(rngS[1]),$(rngS[2])) be small $vMax-shs",
                "Let the $vMax-shs that lie in ($(rngL[1]),$(rngL[2])) be large $vMax-shs",
            )
        end

        # Determine possible small-sh, large-sh combinations
        posCombo = []
        for k = 0:vMax
            upperB = (vMax - k) * toFrac(rngL[2]) + (k) * toFrac(rngS[2])
            lowerB = (vMax - k) * toFrac(rngL[1]) + (k) * toFrac(rngS[1])
            if upperB > m // s > lowerB
                append!(posCombo, k)
            end
        end

        # Diverge from Interval Method
        if length(posCombo) > 0
            mid = formatFrac(1 // 2, cd)
            if output > 1
                combos =
                    ["$k small $vMax-shs and $(vMax-k) large $vMax-shs" for k in posCombo]
                printfT(
                    "Case 5",
                    "All combinations of small and large $vMax-shs are impossible except for:",
                    "",
                    combos...,
                    "",
                    "because the range of muffin amounts for the others does not include $size",
                )

                printHeader("CASE 5: MIDPOINT ANALYSIS")
                printfT(
                    "Property of Buddies",
                    "The # of $vMax-shs in ($j, $mid) is equal to the # of $vMax-shs in ($mid, $k)",
                )
            end

            bounds = [[toFrac(j), 1 // 2], [1 // 2, toFrac(k)], [toFrac(l), toFrac(alpha1)]]
            smallBounds = [toFrac(j), toFrac(k)]
            solid = [toFrac(l), toFrac(alpha1)]
            intervals = [["(", j], ["|", mid], [")[", k], ["](", l], [")", alpha1]]
            labels = ["z $W-shs", "z $W-shs", "0", "$numV $W-shs"]
            intervalDefs = [
                "A = ($j,$mid)",
                "B = ($mid,$k) (|A| = |B| = z)",
                "C = ($l,$alpha1) (|C| = $numV)",
            ]
            if numW < numV
                bounds = [[alpha, toFrac(i)], [toFrac(j), 1 // 2], [1 // 2, toFrac(k)]]
                smallBounds = [alpha, toFrac(i)]
                solid = [alpha, toFrac(i)]
                intervals = [["(", alphaF], [")[", i], ["](", j], ["|", mid], [")", k]]
                labels = ["$numW $V-shs", "0", "z $V-shs", "z $V-shs"]
                intervalDefs = [
                    "A = ($alphaF, $i) (|A| = $numW)",
                    "B = ($j, $mid)",
                    "C = ($mid, $k) (|B| = |C| = z)",
                ]
            end

            if output > 1
                println()
                printf("The following intervals capture what we know about the $V-shs:")
                println("\n", interval(intervals..., labels = labels))
                printLine()
                printfT("Note", "We define the following intervals:", intervalDefs...)
            end

            # Determine possible interval combinations
            posInt = intCombo(m, s, vMax, bounds, smallBounds, posCombo)

            if length(posInt) > 0
                if output > 1
                    printfT(
                        "Note",
                        "Let the notation (a, b, c, ...) denote a combination of $vMax-shs where there are:",
                        "a shs in interval A",
                        "b shs in interval B",
                        "c shs in interval C",
                        "...",
                        "",
                        "A (a, b, c, ...)-student denotes a student with that combination of $vMax-shs",
                    )
                    printfT(
                        "Midpoint Analysis I",
                        "Given the possible small and large $vMax-sh combinations,",
                        "all combinations of $vMax-shs from intervals A, B, and C are impossible except for:",
                        "",
                        posInt...,
                        "",
                        "because the range of muffin amounts for the others does not include $size",
                    )
                end

                # Check if posInt is too large
                if length(posInt) > 5
                    output > 1 && printf(
                        "The # of possible interval combinations is too large, skipping Midpoint Analysis",
                        line = true,
                    )
                    solution = [0, 0]
                else
                    # Set up and solve system of equations
                    iter = 1:length(posInt)
                    (a, b, c) =
                        (getindex.(posInt, 1), getindex.(posInt, 2), getindex.(posInt, 3))
                    (A, B, C) = ("A", "B", "C")
                    if numW < numV
                        (a, b, c) = (b, c, a)
                        (A, B, C) = (B, C, A)
                    end

                    equations = [
                        join(["$(a[i])·x_$i" for i in iter], " + ") *
                        " = " *
                        join(["$(b[i])·x_$i" for i in iter], " + "),
                        join(["$(c[i])·x_$i" for i in iter], " + ") * " = |$C| = $numMin",
                        join(["x_$i" for i in iter], " + ") * " = s_$vMax = $sMax",
                        "",
                        ["where x_$i is the # of $(posInt[i])-students" for i in iter]...,
                    ]

                    # Test all possible non-negative integer solutions with system
                    model = Model(Cbc.Optimizer)
                    set_optimizer_attribute(model, "logLevel", 0)
                    set_optimizer_attribute(model, "cuts", "off")

                    @variable(model, 0 <= x[i = 1:length(iter)], Int)

                    @constraint(model, sum(x .* a) == sum(x .* b))
                    @constraint(model, sum(x .* c) == numMin)
                    @constraint(model, sum(x) == sMax)

                    @objective(model, Max, x[1])
                    optimize!(model)
                    if termination_status(model) == MOI.OPTIMAL
                        x = convert.(Int64, round.(value.(x)))
                        solution = [
                            "One solution is " *
                            join(["x_$i = $(x[i])" for i in iter], ", "),
                            "The system has a non-negative integer solution, so Case 5 is still possible",
                        ]
                    else
                        solution =
                            ["The system has no non-negative integer solutions, so Case 5 is impossible"]
                    end
                end

                if length(solution) > 1
                    if output > 1
                        printf("Midpoint Analysis inconclusive, continuing to Gap Analysis", line=true)
                        printHeader("CASE 5: GAP ANALYSIS")
                    end

                    upChr = join('A':'Z')
                    (loChr, loOrd) = (join('z':-1:'a'), 1)
                    while true
                        # Find gaps
                        gaps = []
                        equations = []
                        for i = 1:length(bounds)
                            ints = []
                            for tup in posInt
                                if tup[i] > 0
                                    upper =
                                        m // s - sum(getindex.(bounds, 1) .* tup) +
                                        bounds[i][1]
                                    lower =
                                        m // s - sum(getindex.(bounds, 2) .* tup) +
                                        bounds[i][2]
                                    disjoint = true
                                    merged = []
                                    for j = 1:length(ints)
                                        if upper <= ints[j][1] || ints[j][2] <= lower
                                            continue
                                        else
                                            disjoint = false
                                            ints[j][1] = min(lower, ints[j][1])
                                            ints[j][2] = max(upper, ints[j][2])
                                            append!(merged, j)
                                        end
                                    end
                                    if length(merged) > 1
                                        append!(
                                            ints,
                                            [[
                                                minimum(ints[j][1] for j in merged),
                                                maximum(ints[j][2] for j in merged),
                                            ]],
                                        )
                                        deleteat!(ints, merged)
                                    end
                                    if disjoint
                                        append!(ints, [[lower, upper]])
                                    end
                                    append!(
                                        equations,
                                        ["For $tup-student: all $vMax-shs in interval $(upChr[i]) must lie in ($(formatFrac(lower, cd)), $(formatFrac(upper, cd)))"],
                                    )
                                end
                            end
                            sort!(ints)
                            append!(equations, [""])
                            for j = 1:length(ints)-1
                                (a, b, c, d) = (
                                    ints[j][2],
                                    ints[j+1][1],
                                    1 - ints[j+1][1],
                                    1 - ints[j][2],
                                )
                                append!(gaps, [[a, b], [c, d]])
                                append!(
                                    equations,
                                    [
                                        "Gap found: [$(formatFrac(a, cd)), $(formatFrac(b, cd))]",
                                        "Its buddy: [$(formatFrac(c, cd)), $(formatFrac(d, cd))]",
                                    ],
                                )
                            end
                            if length(ints) < 2
                                append!(
                                    equations,
                                    ["No gaps found in interval $(upChr[i])"],
                                )
                            end
                            append!(equations, [""])
                        end

                        # Break if no gaps found
                        if length(gaps) == 0
                            break
                        end

                        # Sort and merge gaps
                        sort!(gaps)
                        mergedGaps = []
                        i = 1
                        while i < length(gaps)
                            if gaps[i][2] < gaps[i+1][1]
                                append!(mergedGaps, [gaps[i]])
                            else
                                append!(mergedGaps, [[gaps[i][1], gaps[i+1][2]]])
                                i += 1
                            end
                            i += 1
                        end
                        if i == length(gaps)
                            append!(mergedGaps, [gaps[i]])
                        end
                        gaps = mergedGaps

                        if output > 1
                            printfT("New Gaps", equations...)
                        end

                        # Insert gaps
                        for gap in gaps
                            (l, u) = (gap[1], gap[2])
                            i = 1
                            while i <= length(bounds)
                                (bl, bu) = (bounds[i][1], bounds[i][2])
                                if bl < l <= u < bu
                                    # Update bounds
                                    insert!(bounds, i + 1, [u, bu])
                                    bounds[i][2] = l
                                    # Update intervals
                                    for j = 2:length(intervals)
                                        if toFrac(intervals[j-1][2]) == bl &&
                                           toFrac(intervals[j][2]) == bu
                                            insert!(intervals, j, ["](", formatFrac(u, cd)])
                                            insert!(intervals, j, [")[", formatFrac(l, cd)])
                                            break
                                        end
                                    end
                                end
                                i += 1
                            end
                        end

                        # Fail if gap finding excessive
                        if length(bounds) > 26
                            output > 0 &&
                                printf("Excessive gap finding, VGap failed", line = true)
                            return false
                        end

                        # Update posInt
                        posInt = intCombo(m, s, vMax, bounds, smallBounds, posCombo)
                    end

                    if length(bounds) == 3
                        output > 0 && printf(
                            "No further gaps could be generated, VGap failed",
                            line = true,
                        )
                        return false
                    end

                    # Update labels and intervalDefs
                    labels = repeat(["0"], length(intervals) - 1)
                    intervalDefs = repeat([""], length(bounds))
                    equivs = []
                    solidInt = 0
                    visited = []
                    iter = 1:length(intervals)-1
                    for i = 1:length(bounds)
                        if !(i in visited)
                            int = bounds[i]
                            chr = upChr[i]
                            for k in iter
                                if toFrac(intervals[k][2]) == int[1] &&
                                   toFrac(intervals[k+1][2]) == int[2]
                                    (lower, upper) = (intervals[k][2], intervals[k+1][2])
                                    if int == solid
                                        labels[k] = "$numMin $vMax-shs"
                                        intervalDefs[i] = "$chr = ($lower, $upper) (|$chr| = $numMin)"
                                        solidInt = i
                                    else
                                        chrL = loChr[loOrd]
                                        labels[k] = "$chrL $vMax-shs"
                                        intervalDefs[i] = "$chr = ($lower, $upper)"
                                        equiv = [i]
                                        for j = i:length(bounds)
                                            intOth = bounds[j]
                                            if intOth == [1 - int[2], 1 - int[1]] || (
                                                vMax == 2 &&
                                                intOth ==
                                                [m // s - int[2], m // s - int[1]]
                                            )
                                                for l in iter
                                                    if toFrac(intervals[l][2]) ==
                                                       intOth[1] &&
                                                       toFrac(intervals[l+1][2]) ==
                                                       intOth[2]
                                                        (lowerOth, upperOth) = (
                                                            intervals[l][2],
                                                            intervals[l+1][2],
                                                        )
                                                        labels[l] = "$chrL $vMax-shs"
                                                        intervalDefs[j] = "$(upChr[j]) = ($lowerOth, $upperOth)"
                                                        append!(equivs, [(i, j)])
                                                        append!(equiv, j)
                                                        append!(visited, j)
                                                    end
                                                end
                                            end
                                        end
                                        intervalDefs[equiv[length(
                                            equiv,
                                        )]] *= " ($(join(["|$(upChr[j])|" for j in equiv], " = ")) = $chrL)"
                                        loOrd += 1
                                    end
                                    break
                                end
                            end
                        end
                    end

                    if output > 1
                        println()
                        printf("After finding gaps and buddying, this is the interval of the $vMax-shs:")
                        println("\n", interval(intervals..., labels = labels))
                        printLine()
                        printfT(
                            "Note",
                            "We define the following intervals:",
                            intervalDefs...,
                        )
                    end

                    if length(posInt) > 0
                        if output > 1
                            printfT(
                                "Gap Analysis I",
                                "Given the possible small and large $vMax-sh combinations,",
                                "all combinations of $vMax-shs from the defined intervals are impossible except for:",
                                "",
                                posInt...,
                                "",
                                "because the range of muffin amounts for the others does not include $size",
                            )
                        end

                        # Set up and solve system of equations
                        iter = 1:length(posInt)
                        consts = [getindex.(posInt, i) for i = 1:length(bounds)]
                        l = solidInt

                        equations = []
                        for (j, k) in equivs
                            append!(
                                equations,
                                [
                                    join(["$(consts[j][i])·x_$i" for i in iter], " + ") *
                                    " = " *
                                    join(["$(consts[k][i])·x_$i" for i in iter], " + "),
                                ],
                            )
                        end
                        if l > 0
                            append!(
                                equations,
                                [
                                    join(["$(consts[l][i])·x_$i" for i in iter], " + ") * " = |$(upChr[l])| = $numMin",
                                ],
                            )
                        end
                        append!(
                            equations,
                            [
                                join(["x_$i" for i in iter], " + ") * " = s_$vMax = $sMax",
                                "",
                            ],
                        )

                        append!(
                            equations,
                            [
                                "where x_$i is the # of ($(join(posInt[i], ",")))-students"
                                for i in iter
                            ],
                        )

                        # Test all possible non-negative integer solutions with system
                        model = Model(Cbc.Optimizer)
                        set_optimizer_attribute(model, "logLevel", 0)
                        set_optimizer_attribute(model, "cuts", "off")

                        @variable(model, 0 <= x[i = 1:length(iter)], Int)

                        for (j, k) in equivs
                            @constraint(model, sum(x .* consts[j]) == sum(x .* consts[k]))
                        end
                        if l > 0
                            @constraint(model, sum(x .* consts[l]) == numMin)
                        end
                        @constraint(model, sum(x) == sMax)

                        @objective(model, Max, x[1])
                        optimize!(model)
                        if termination_status(model) == MOI.OPTIMAL
                            x = convert.(Int64, round.(value.(x)))
                            solution = [
                                "One solution is " *
                                join(["x_$i = $(x[i])" for i in iter], ", "),
                                "The system has a non-negative integer solution, so Case 5 is still possible",
                            ]
                        else
                            solution =
                                ["The system has no non-negative integer solutions, so Case 5 is impossible"]
                        end

                        if output > 1
                            printfT(
                                "Gap Analysis II",
                                "Given the new interval equalities, we can set up the following system of equations:",
                                "",
                                equations...,
                                "",
                                solution...,
                            )
                        end

                        if length(solution) > 1
                            output == 1 && printf(
                                "Could not generate conclusive system of equations, VGap failed",
                                line = true,
                            )
                            return false
                        end
                    else
                        if output > 1
                            printfT(
                                "Gap Analysis",
                                "Given the possible small and large $vMax-sh combinations,",
                                "all combinations of $vMax-shs from the intervals are impossible since their total amount can not be $size",
                            )
                        end
                    end

                else
                    if output > 1
                        printfT(
                            "Midpoint Analysis II",
                            "Since |$A| = |$B|, we can set up the following system of equations:",
                            "",
                            equations...,
                            "",
                            solution...,
                        )
                    end
                end
            else
                if output > 1
                    printfT(
                        "Midpoint Analysis",
                        "Given the possible small and large $vMax-sh combinations,",
                        "all combinations of $vMax-shs from intervals A, B, and C are impossible since their total amount can not be $size",
                    )
                end
            end

        else
            if output > 1
                printfT(
                    "Case 5",
                    "All combinations of small and large $vMax-shs are impossible since their total amount can not be $size",
                )
            end
        end

        # Conclude proof
        alpha = formatFrac(alpha)
        if output > 1
            # Conclude with alpha's value
            printHeader("CONCLUSION")
            printfT(
                "Conclusion",
                "All possible cases contradict the assumption iff. α ≥ $alpha",
                "muffins($m,$s) ≤ α, ∀ α ≥ $alpha",
                "",
                "muffins($m,$s) ≤ $alpha",
            )
            printEnd()
        elseif output > 0
            printfT("Gap Method", "Upper bound α of muffins($m,$s) is $alpha")
            printEnd()
        end

        true
    end
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkSecurityGroupResult <: SwaggerModel
    securityRuleAccessResult::Any # spec type: Union{ Nothing, SecurityRuleAccess } # spec name: securityRuleAccessResult
    evaluatedNetworkSecurityGroups::Any # spec type: Union{ Nothing, Vector{EvaluatedNetworkSecurityGroup} } # spec name: evaluatedNetworkSecurityGroups

    function NetworkSecurityGroupResult(;securityRuleAccessResult=nothing, evaluatedNetworkSecurityGroups=nothing)
        o = new()
        validate_property(NetworkSecurityGroupResult, Symbol("securityRuleAccessResult"), securityRuleAccessResult)
        setfield!(o, Symbol("securityRuleAccessResult"), securityRuleAccessResult)
        validate_property(NetworkSecurityGroupResult, Symbol("evaluatedNetworkSecurityGroups"), evaluatedNetworkSecurityGroups)
        setfield!(o, Symbol("evaluatedNetworkSecurityGroups"), evaluatedNetworkSecurityGroups)
        o
    end
end # type NetworkSecurityGroupResult

const _property_map_NetworkSecurityGroupResult = Dict{Symbol,Symbol}(Symbol("securityRuleAccessResult")=>Symbol("securityRuleAccessResult"), Symbol("evaluatedNetworkSecurityGroups")=>Symbol("evaluatedNetworkSecurityGroups"))
const _property_types_NetworkSecurityGroupResult = Dict{Symbol,String}(Symbol("securityRuleAccessResult")=>"SecurityRuleAccess", Symbol("evaluatedNetworkSecurityGroups")=>"Vector{EvaluatedNetworkSecurityGroup}")
Base.propertynames(::Type{ NetworkSecurityGroupResult }) = collect(keys(_property_map_NetworkSecurityGroupResult))
Swagger.property_type(::Type{ NetworkSecurityGroupResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkSecurityGroupResult[name]))}
Swagger.field_name(::Type{ NetworkSecurityGroupResult }, property_name::Symbol) =  _property_map_NetworkSecurityGroupResult[property_name]

function check_required(o::NetworkSecurityGroupResult)
    true
end

function validate_property(::Type{ NetworkSecurityGroupResult }, name::Symbol, val)
end


text:
export list, image, link, abbr

@api list => (List <: Tile) begin
    doc("Stylize a list of tiles as an itemized list.")
    curry(tiles::TileList , doc="A tile or a vector of tiles.")
    kwarg(ordered::Bool=false, doc=md"If set to `true`, numbering will be used.")
end

render(l::List, state) =
    Elem(l.ordered ? :ol : :ul,
         map(x -> Elem(:li, render(x, state)), l.tiles.tiles))

@api image => (Image <: Tile) begin
    doc(md"""Show an image from a `url`. To read an image and display it, use
             [`Images.imread`](https://github.com/timholy/Images.jl#readme).""")
    arg(url::AbstractString, doc="The url of the image.")
    kwarg(alt::AbstractString="", doc="Text to display if image does not load.")
end

render(i::Image, state) =
    Elem(:img, src=i.url, alt=i.alt, style=@d("width"=>"auto", "height"=>"auto", "display" => "block"))

@api link => (Hyperlink <: Tile) begin
    doc("A hyperlink.")
    arg(url::AbstractString, doc="The destination of the link.")
    curry(
        tiles::TileList,
        doc="""A tile or a vector of tiles. These tiles link to the url."""
    )
end

render(a::Hyperlink, state) =
    Elem(:a, render(a.tiles, state), href=a.url)

@api abbr => (Abbr <: Tile) begin
    doc(md"""An abbreviation. When you hover over an abbreviation, the `title` is
          shown in a tooltip.""")
    arg(title::AbstractString, doc="The title to show.")
    curry(
        tiles::TileList,
        doc="A tile or a vector of tiles. Hovering over these tiles triggers the tooltip."
    )
end

render(a::Abbr, state) =
    Elem(:abbr, render(a.tiles, state), title=a.title)


text:
@testset "Solver: Random QP w/ equalities (QDLDL)" begin
    """
        minimize   x' P x + q' x
        subject to    Ax = b
    """

    #####
    n = 10
    m = 3

    nw = n + m
    nθ = 2n + m * n + m

    function obj(z, θ)
        x = z[1:n]

        P = Diagonal(θ[1:n])
        p = θ[n .+ (1:n)]
        return transpose(x) * P * x + transpose(p) * x
    end

    function constraints(z, θ)
        x = z[1:n]
        A = reshape(θ[2n .+ (1:(m * n))], m, n)
        b = θ[2n + m * n .+ (1:m)]

        A * x - b
    end

    function lagrangian(w, θ)
        z = w[1:n]
        y = w[n .+ (1:m)]

        L = 0.0
        L += obj(z, θ)
        L += dot(y, constraints(z, θ))

        return L
    end

    # residual
    function _r(w, θ, κ)
        L = lagrangian(w, θ)
        Lz = Symbolics.gradient(L, w[1:n])
        c = constraints(w[1:n], θ)

        [
            Lz;
            c;
        ]
    end

    @variables z_sym[1:nw]
    @variables θ_sym[1:nθ]
    @variables κ_sym[1:1]

    r_sym = _r(z_sym, θ_sym, κ_sym)
    rf! = eval(Symbolics.build_function(r_sym, z_sym, θ_sym, κ_sym)[2])
    rz_exp = Symbolics.jacobian(r_sym, z_sym)
    rθ_exp = Symbolics.jacobian(r_sym, θ_sym)
    rz_sp = similar(rz_exp, Float64)
    rθ_sp = similar(rθ_exp, Float64)
    rzf! = eval(Symbolics.build_function(rz_exp, z_sym, θ_sym)[2])
    rθf! = eval(Symbolics.build_function(rθ_exp, z_sym, θ_sym)[2])

    x0 = randn(n)
    A = rand(m, n)
    b = A * x0
    z = [randn(n); zeros(m)]
    θ = [ones(n); zeros(n); vec(A); b]
    κ = [1.0]

    r = zeros(nw)
    rz = zeros(nw, nw)
    rθ = zeros(nw, nθ)

    opts = CALIPSO.InteriorPointOptions(diff_sol=true,
        undercut=10.0,
        max_ls=25,
        max_iter=100,
        verbose=false,
        solver=:lu_solver)

    idx = CALIPSO.IndicesOptimization(
        nw, nw,
        [collect(1:0), collect(1:0)],
        [collect(1:0), collect(1:0)],
        Vector{Vector{Int}}[],
        Vector{Vector{Int}}[],
        collect(1:(n + m)),
        Vector{Int}(),
        Vector{Int}(),
        Vector{Int}[],
        Vector{Int}())

    # solver
    ip = CALIPSO.interior_point(z, θ,
        idx=idx,
        r! = rf!, rz! = rzf!, rθ! = rθf!,
        rz = rz,
        rθ = rθ,
        opts = opts)

    status = CALIPSO.interior_point_solve!(ip)
    # @benchmark CALIPSO.interior_point_solve!($ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test norm(A * ip.z[1:n] - b, Inf) < 1.0e-6
    @test norm(ip.δz, 1) != 0.0
end


text:
module StructuralIdentifiability

# General purpose packages
using Base.Iterators
using Dates
using IterTools
using LinearAlgebra
using Logging
using MacroTools
using Primes
using DataStructures

# Algebra packages
using AbstractAlgebra
using Nemo
using GroebnerBasis
using Singular

# For testing (TODO: move to the test-specific dependencies)
using Test
using TestSetExtensions
using IterTools
using MacroTools

using ModelingToolkit

# defining a model
export ODE, @ODEmodel, PreprocessODE

# assessing identifiability
export assess_local_identifiability, assess_global_identifiability, assess_identifiability

# auxuliary function
export set_parameter_values

# extra functionality
export find_ioequations, find_identifiable_functions

# exporting to other formats
export print_for_maple, print_for_DAISY, print_for_COMBOS, print_for_GenSSI

# function for creating linear compartment models
export linear_compartment_model

# structure for storing the io-equations for doing reductions further
export PBRepresentation, diffreduce, io_switch!, pseudodivision

# would be great to merge with the Julia logger
_runtime_logger = Dict()

include("util.jl")
include("power_series_utils.jl")
include("ODE.jl")
include("ODEexport.jl")
include("local_identifiability.jl")
include("wronskian.jl")
include("elimination.jl")
include("primality_check.jl")
include("io_equation.jl")
include("global_identifiability.jl")
include("lincomp.jl")
include("pb_representation.jl")

"""
    assess_identifiability(ode::ODE{P}, p::Float64=0.99) where P <: MPolyElem{fmpq}
Input:
- `ode` - the ODE model
- `p` - probability of correctness.

Assesses identifiability (both local and global) of a given ODE model (parameters detected automatically). The result is guaranteed to be correct with the probability
at least `p`.

"""
function assess_identifiability(ode::ODE{P}, p::Float64 = 0.99) where {P<:MPolyElem{fmpq}}
    result_list = assess_identifiability(ode, ode.parameters, p)
    return Dict(param => res for (param, res) in zip(ode.parameters, result_list))
end

"""
    assess_identifiability(ode, [funcs_to_check, p=0.99])

Input:
- `ode` - the ODE model
- `p` - probability of correctness.

Assesses identifiability of a given ODE model. The result is guaranteed to be correct with the probability
at least `p`.

If `funcs_to_check` are given, then the function will assess the identifiability of the provided functions
and return a list of the same length with each element being one of `:nonidentifiable`, `:locally`, `:globally`.

If `funcs_to_check` are not given, the function will assess identifiability of the parameters, and the result will
be a dictionary from the parameters to their identifiability properties (again, one of `:nonidentifiable`, `:locally`, `:globally`).
"""
function assess_identifiability(ode::ODE{P}, funcs_to_check::Array{<:RingElem,1}, p::Float64 = 0.99) where {P<:MPolyElem{fmpq}}
    p_glob = 1 - (1 - p) * 0.9
    p_loc = 1 - (1 - p) * 0.1

    @info "Assessing local identifiability"
    runtime = @elapsed local_result, bound = assess_local_identifiability(ode, funcs_to_check, p_loc, :ME)
    @info "Local identifiability assessed in $runtime seconds"
    _runtime_logger[:loc_time] = runtime

    if bound > 1
        @warn "For this model single-experiment identifiable functions are not the same as multi-experiment identifiable. The analysis performed by this program is right now for multi-experiment only. If you still  would like to assess single-experiment identifiability, we recommend using SIAN (https://github.com/alexeyovchinnikov/SIAN-Julia)"
        @debug "Bound: $bound"
    end

    locally_identifiable = Array{Any,1}()
    for (loc, f) in zip(local_result, funcs_to_check)
        if loc
            push!(locally_identifiable, f)
        end
    end

    @info "Assessing global identifiability"
    runtime = @elapsed global_result = assess_global_identifiability(ode, locally_identifiable, p_glob)
    @info "Global identifiability assessed in $runtime seconds"
    _runtime_logger[:glob_time] = runtime

    result = Array{Symbol,1}()
    glob_ind = 1
    for i in 1:length(funcs_to_check)
        if !local_result[i]
            push!(result, :nonidentifiable)
        else
            if global_result[glob_ind]
                push!(result, :globally)
            else
                push!(result, :locally)
            end
            glob_ind += 1
        end
    end

    return result
end

"""
    assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], p = 0.99)
Input:
- `ode` - the ModelingToolkit.ODESystem object that defines the model
- `measured_quantities` - the output functions of the model
- `funcs_to_check` - functions of parameters for which to check the identifiability
- `p` - probability of correctness.

Assesses identifiability (both local and global) of a given ODE model (parameters detected automatically). The result is guaranteed to be correct with the probability
at least `p`.

"""
function assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], p = 0.99)
    if length(measured_quantities)==0
        if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode))
            @info "Measured quantities are not provided, trying to find the outputs in input ODE."
            measured_quantities = filter(eq->(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode))
        else
            throw(error("Measured quantities (output functions) were not provided and no outputs were found."))
        end
    end
    if length(funcs_to_check) == 0
        funcs_to_check = ModelingToolkit.parameters(ode)
    end
    ode, syms, gens_ = PreprocessODE(ode, measured_quantities)
    out_dict = Dict{Num,Symbol}()
    funcs_to_check_ = [eval_at_nemo(each, Dict(syms .=> gens_)) for each in funcs_to_check]
    tmp = Dict(param => res for (param, res) in zip(funcs_to_check_, assess_identifiability(ode, funcs_to_check_, p)))
    nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check)
    for (func, res) in pairs(tmp)
        out_dict[nemo2mtk[func]] = res
    end
    return out_dict
end

end


text:
const DEFAULT_URI = "https://management.azure.com"
const AUTH_URI = "https://login.microsoftonline.com"
const DEFAULT_TOKEN_EXPIRE = 5*60

"""
Auth provider provides a way to retrieve credentials for authentication / re-authentication after token expiry.
"""
abstract type AzureAuthProvider end

"""
A simple auth provider that stores secrets in memory.
"""
mutable struct AzureCredentials <: AzureAuthProvider
    tenant_id::String
    client_id::String
    client_secret::String
end

tenant_id(creds::AzureCredentials) = creds.tenant_id
client_id(creds::AzureCredentials) = creds.client_id
client_secret(creds::AzureCredentials) = creds.client_secret

"""
Holds a swagger client to make API calls with, an auth provider to authenticate to Azure with, and an authenticated token.
"""
mutable struct AzureContext
    client::Swagger.Client
    auth_provider::AzureAuthProvider
    token::Dict{String,Any}
    expires::Float64

    function AzureContext(auth::AzureAuthProvider)
        new(Swagger.Client(DEFAULT_URI), auth, Dict{Symbol,Any}(), 0)
    end
end

"""
Authenticates the context if not authenticated yet.
Re-authenticates if token has expired.
"""
function authenticate(ctx::AzureContext)
    (time() < ctx.expires) && return

    @info("authenticating...")
    data = Dict{String,String}(
        "client_id" => client_id(ctx.auth_provider),
        "client_secret" => client_secret(ctx.auth_provider),
        "resource" => DEFAULT_URI*"/",  # must end with /
        "grant_type" => "client_credentials"
    )

    headers = Dict{String,String}(
        "Cookie" => "flight-uxoptin=true; stsservicecookie=ests; x-ms-gateway-slice=productionb; stsservicecookie=ests",
        "Content-Type" => "application/x-www-form-urlencoded",
        "Connection" => "close"
    )

    tenant = tenant_id(ctx.auth_provider)
    auth_url = AUTH_URI * "/" * tenant * "/oauth2/token"

    input = PipeBuffer()
    write(input, URIs.escapeuri(data))
    output = IOBuffer()
    resp = Downloads.request(auth_url; method="POST", input=input, output=output, headers=headers)
    (isa(resp, Downloads.Response) && (200 <= resp.status <= 206)) || throw(Swagger.ApiException(resp))

    ctx.token = JSON.parse(String(take!(output)))
    if "expires_in" in keys(ctx.token)
        ctx.expires = time() + parse(Int, ctx.token["expires_in"])
    elseif "expires_on" in keys(ctx.token)
        ctx.expires = parse(Float64, ctx.token["expires_on"])
    elseif "not_before" in keys(ctx.token)
        ctx.expires = parse(Float64, ctx.token["not_before"])
    else
        ctx.expires = time() + DEFAULT_TOKEN_EXPIRE
    end

    token_type = get(ctx.token, "token_type", "Bearer")
    ctx.client.headers["Authorization"] = token_type * " " * ctx.token["access_token"]

    @info("authenticated")
    nothing
end

"""
Ensures context is authenticated, returns the specified API context to use in an API call
"""
function api(ctx::AzureContext, apicls)
    authenticate(ctx)
    apicls(ctx.client)
end

apiver(mod::Module) = _module_versions[mod]
apiver(api::DataType) = _api_versions[api]

export api, apiver, AzureContext, AzureAuthProvider, AzureCredentials


text:
# this file contains some information theoretic functionals of probability distributions
# should ultimately spin this off into separate package
# functional benchmarked just as fast as direct methods!!!

# TODO make this work for functions with ℵ₀ support
function functional(d::Distribution{T, Discrete}, f::Function, st::Function=(x -> true)) where T
    o = 0.0
    for k ∈ support(d)
        if st(k)
            o += f(k)
        end
    end
    o
end

function functional(d::SampleDistribution, f::Function, st::Function=(x -> true))
    o = 0.0
    for i ∈ 1:ncategories(d)
        if st(getsample(d, i))
            o += f(pdfbyindex(d, i))
        end
    end
    o
end
export functional

Hfunc(p::Real) = -p*log(p)
H(d::SampleDistribution, st::Function=(x -> true)) = functional(d, Hfunc, st)
export H


text:
function interpolate_output!(x_out::Vector{Float64}, rr::RadauIntegrator{T_object, N, n_stage_max}, h::Float64, is_extrapolate_prev::Bool=true) where {n_stage_max, N<:Int64, T_object}
    rr.dense.is_has_X || error("no dense output to interpolate")
    dense = rr.dense
    theta = h / dense.h + is_extrapolate_prev
    theta_prod = theta
    n_stage = dense.s
    theta_vec = ones(n_stage)
    theta_vec = cumprod(theta_vec)
    table = rr.table[n_stage]
    bi = table.bi
    weights = bi * theta_vec
    LinearAlgebra.BLAS.blascopy!(N, dense.X_stage[1], 1, x_out, 1)
    weight_1 = weights[1]
    LinearAlgebra.BLAS.scal!(N, weight_1, x_out, 1)
    for k = 2:n_stage
        weight_k = weights[k]
        LinearAlgebra.BLAS.axpy!(weight_k, dense.X_stage[k], x_out)
    end
    return nothing
end

function initialize_X_with_interp!(rr::RadauIntegrator{T_object, N, n_stage_max}, table::RadauTable{n_stage}) where {n_stage, n_stage_max, N, T_object}
    h = rr.step.h
    c = table.c
    for k = 1:n_stage
        interpolate_output!(rr.ct.X_stage[k], rr, h * c[k], true)
    end
    return nothing
end

@inline function update_dense_unsuccessful!(rr::RadauIntegrator{T_object, N, n_stage_max}) where {n_stage_max, N, T_object}
    rr.dense.is_has_X = false
    return nothing
end

function update_dense_successful!(rr::RadauIntegrator{T_object, N, n_stage_max}, table::RadauTable{n_stage}) where {n_stage_max, N, T_object, n_stage}
    rr.dense.is_has_X = true
    rr.dense.h = rr.step.h
    rr.dense.s = n_stage
    for i = 1:n_stage
        LinearAlgebra.BLAS.blascopy!(N, rr.ct.X_stage[i], 1, rr.dense.X_stage[i], 1)
    end
    return nothing
end


text:
using Flux, MLDataPattern, Mill, JsonGrinder, JSON, IterTools, Statistics, BenchmarkTools, ThreadTools
import JsonGrinder: suggestextractor, ExtractCategorical, ExtractBranch, ExtractString, MultipleRepresentation, extractscalar
import Mill: mapdata, sparsify, reflectinmodel

###############################################################
# start by loading all samples
###############################################################
samples = open("examples/recipes.json","r") do fid
	Vector{Dict}(JSON.parse(read(fid, String)))
end
JSON.print(samples[1],2)

###############################################################
# create schema of the JSON
###############################################################
sch = JsonGrinder.schema(samples)

###############################################################
# create extractor and split it into one for loading targets and
# one for loading data, using custom function to set conditions for using n-gram representation
###############################################################
delete!(sch.childs,"id")

limituse(d::Dict{T,Int}, limit) where {T<:AbstractString} = String.(limituse(d, limit))
limituse(d::Dict{T,Int}, limit) where {T} = collect(filter(k -> d[k] >= limit, keys(d)))

function custom_scalar_extractor()
	[(e -> promote_type(unique(typeof.(keys(e.counts)))...) <: String,
		e -> MultipleRepresentation((ExtractCategorical(limituse(e.counts, 10)), ExtractString(String)))),
	 (e -> (length(keys(e.counts)) / e.updated < 0.1  && length(keys(e.counts)) <= 10000),
		e -> ExtractCategorical(collect(keys(e.counts)))),
	 (e -> true,
		e -> extractscalar(promote_type(unique(typeof.(keys(e.counts)))...))),]
end

extractor = suggestextractor(sch, (scalar_extractors=custom_scalar_extractor(), mincount=100,))

extract_data = ExtractBranch(nothing,deepcopy(extractor.other))
extract_target = ExtractBranch(nothing,deepcopy(extractor.other))
delete!(extract_target.other,"ingredients")
delete!(extract_data.other,"cuisine")
extract_target.other["cuisine"] = JsonGrinder.ExtractCategorical(keys(sch["cuisine"]))

extract_data(JsonGrinder.sample_synthetic(sch))
###############################################################
# we convert JSONs to Datasets
###############################################################
# advised to use all the samples, this is just speedup to demonstrate functionality
data = tmap(extract_data, samples[1:5_000])
data = reduce(catobs, data)
target = tmap(extract_target, samples[1:5_000])
target = reduce(catobs, target).data

e = sch["cuisine"]

###############################################################
# 	create the model according to the data
###############################################################
m = reflectinmodel(sch, extract_data,
	k -> Dense(k,20,relu),
	d -> SegmentedMeanMax(d),
	b = Dict("" => k -> Dense(k, size(target, 1))),
)

###############################################################
#  train
###############################################################
opt = Flux.Optimise.ADAM()
loss = (x,y) -> Flux.logitcrossentropy(m(x).data, y)
valdata = data[1:1_000],target[:,1:1_000]
data, target = data[1_001:5_000], target[:,1001:5_000]
# for less recourse-chungry training, we use only part of data for trainng, but it is advised to used all, as i following line:
# data, target = data[1001:nobs(data)], target[:,1001:size(target,2)]

cb = () -> println("accuracy = ",mean(Flux.onecold(m(valdata[1]).data) .== Flux.onecold(valdata[2])))
batch_size = 100
#iterations = 1000
iterations = 100
ps = Flux.params(m)
mean(Flux.onecold(m(data).data) .== Flux.onecold(target))

@info "testing the gradient"
gs = gradient(() -> loss(data, target), ps)
Flux.Optimise.update!(opt, ps, gs)

Flux.Optimise.train!(loss, ps, repeatedly(() -> (data, target), 20), opt, cb = Flux.throttle(cb, 2))
# feel free to train for longer period of time, this example is learns only 20 iterations, so it runs fast
# Flux.Optimise.train!(loss, ps, repeatedly(() -> (data, target), 500), opt, cb = Flux.throttle(cb, 10))

#calculate the accuracy
mean(Flux.onecold(m(data).data) .== Flux.onecold(target))

# samples = open("recipes_test.json","r") do fid
# 	Array{Dict}(JSON.parse(readstring(fid)))
# end
# ids = map(s -> s["id"],samples)
# tstdata = @>> samples map(extract_data);
# tstdata = cat(tstdata...);
# tstdata = mapdata(sentence2ngrams,tstdata)
# tstdata = mapdata(i -> sparsify(Float32.(i),0.05),tstdata)
# names = extract_target.other["cuisine"].items
# y = map(i -> names[i],Flux.argmax(m(tstdata)));

# using DataFrames
# CSV.write("cuisine.csv",DataFrame(id = ids,cuisine = y ),delim=',')

# schema can also be created in parallel for better performance, compare:

# single threaded
# @btime JsonGrinder.schema(samples)
# multi threaded
# @btime merge(tmap(x->JsonGrinder.schema(collect(x)), Iterators.partition(samples, 5_000))...)


text:
#    Problem : GROUP A
#    *********
#	 A multi-dimensional variant of HUMPS, a two dimensional function
#    with a lot of humps. The density of humps increases with the
#    parameter ZETA, making the problem more difficult.
#
#	 The problem is nonconvex.
#
#    Origonal SIF Source:
#    Ph. Toint, private communication, 1997
#
#    GENHUMP.SIF classification OUR2-AN-V-0
#
#    Number of variables is variable
#
# Daniel Henderson, 08/2021
f = (x) -> begin
	fx = 0.0
    ζ = 20.0
    for i in firstindex(x):lastindex(x)-1
        fx += sin(ζ*x[i])^2*sin(ζ*x[i+1])^2+0.05*(x[i]^2+x[i+1]^2)
    end
    return fx
end

g! = (g, x) -> begin
    ζ = 20.0
    n = lastindex(x)

    g[1] = ζ*sin(2ζ*x[1])*sin(ζ*x[2])^2+0.1x[1]
    for i in 1:n-2
        g[i+1] = 0.2x[i+1]+40cos(20x[i+1])sin(20x[i])^2*sin(20x[i+1])+40cos(20x[i+1])sin(20x[i+1])sin(20x[i+2])^2
    end
    g[n] = 2ζ*cos(ζ*x[n])*sin(ζ*x[n-1])^2*sin(ζ*x[n])+0.1x[n]
    return g
end

fg! = (g, x) -> begin
    ζ = 20.0
    n = lastindex(x)
    fx = 0.0
    g[1] = ζ*sin(2ζ*x[1])*sin(ζ*x[2])^2+0.1x[1]
    for i in 1:n-2
        fx += sin(ζ*x[i])^2*sin(ζ*x[i+1])^2+0.05*(x[i]^2+x[i+1]^2)
        g[i+1] = 0.2x[i+1] + 40cos(20x[i+1])sin(20x[i])^2*sin(20x[i+1])+40cos(20x[i+1])sin(20x[i+1])sin(20x[i+2])^2
    end
    fx += sin(ζ*x[n-1])^2*sin(ζ*x[n])^2+0.05*(x[n-1]^2+x[n]^2)
    g[n] = 2ζ*cos(ζ*x[n])*sin(ζ*x[n-1])^2*sin(ζ*x[n])+0.1x[n]
    return g
end

init = (n::Int=5000) -> begin
	n < 2 && @warn("GENHUMPS: number of variables must be ≥ 2")
	n = max(n, 2)
    x0 = -506.2*ones(n)
    x0[1] = -506
    return n, x0
end

TestSet["GENHUMPS"] = UncProgram("GENHUMPS", f, g!, fg!, init)


text:
using HTMLParser
using Test

@testset "HTMLParser.jl" begin
    # Write your tests here.
end


text:
##############################################################################
##
## Least Squares via Cholesky Factorization
##
##############################################################################

type CholfactFixedEffectProblem <: FixedEffectProblem
    fes::Vector{FixedEffect}
    m::SparseMatrixCSC{Float64,Int}
    chol::Base.SparseArrays.CHOLMOD.Factor{Float64}
    x::Vector{Float64}
end

function FixedEffectProblem(fes::Vector{FixedEffect}, ::Type{Val{:cholesky}})
    m = sparse(fes)
    chol = cholfact(At_mul_B(m, m))
    total_len = reduce(+, map(fe -> sum(fe.scale .!= 0), fes))
    CholfactFixedEffectProblem(fes, m, chol, Array(Float64, total_len))
end

function solve!(fep::CholfactFixedEffectProblem, r::AbstractVector{Float64} ; kwargs...)
    fep.chol \ At_mul_B!(fep.x, fep.m, r)
end

##############################################################################
##
## Least Squares via QR Factorization
##
##############################################################################

type QRfactFixedEffectProblem <: FixedEffectProblem
    fes::Vector{FixedEffect}
    m::SparseMatrixCSC{Float64,Int}
    qr::Base.SparseArrays.SPQR.Factorization{Float64}
    b::Vector{Float64}
end

function FixedEffectProblem(fes::Vector{FixedEffect}, ::Type{Val{:qr}})
    m = sparse(fes)
    qr = qrfact(m)
    b = Array(Float64, length(fes[1].refs))
    QRfactFixedEffectProblem(fes, m, qr, b)
end

function solve!(fep::QRfactFixedEffectProblem, r::AbstractVector{Float64} ; kwargs...)
    # since \ needs a vector
    copy!(fep.b, r)
    fep.qr \ fep.b
end

##############################################################################
##
## Methods used by all matrix factorization
##
##############################################################################

# construct the sparse matrix of fixed effects A in  A'Ax = A'r
function Base.sparse(fes::Vector{FixedEffect})
    # construct model matrix A constituted by fixed effects
    nobs = length(fes[1].refs)
    N = length(fes) * nobs
    I = Array(Int, N)
    J = similar(I)
    V = Array(Float64, N)
    start = 0
    idx = 0
    for fe in fes
       for i in 1:length(fe.refs)
           idx += 1
           I[idx] = i
           J[idx] = start + fe.refs[i]
           V[idx] = fe.interaction[i] * fe.sqrtw[i]
       end
       start += sum(fe.scale .!= 0)
    end
    sparse(I, J, V)
end

# updates r as the residual of the projection of r on A
function solve_residuals!(fep::Union{CholfactFixedEffectProblem, QRfactFixedEffectProblem}, r::AbstractVector{Float64}; kwargs...)
    x = solve!(fep, r; kwargs...)
    A_mul_B!(-1.0, fep.m, x, 1.0, r)
    return r, 1, true
end

# solves A'Ax = A'r
# transform x from Vector{Float64} (stacked vector of coefficients)
# to Vector{Vector{Float64}} (vector of coefficients for each categorical variable)
function solve_coefficients!(fep::Union{CholfactFixedEffectProblem, QRfactFixedEffectProblem}, r::AbstractVector{Float64}; kwargs...)
    x = solve!(fep, r; kwargs...)
    out = Vector{Float64}[]
    iend = 0
    for fe in get_fes(fep)
        istart = iend + 1
        iend = istart + sum(fe.scale .!= 0) - 1
        push!(out, x[istart:iend])
    end
    return out, 1, true
end

get_fes(fep::Union{CholfactFixedEffectProblem, QRfactFixedEffectProblem}) = fep.fes


text:
# Exercise 5-1
# Infinite recursion will terminate when stack maximum is reached.
function recurse()
    recurse()
end

# recurse()

# Exercise 5-2
function daytime()
    reftime = round(Int64, time())
    minutes, seconds = divrem(reftime, 60)
    hours, minutes = divrem(minutes, 60)
    days, hours = divrem(hours, 24)
    println("Number of days since epoch: $days")
    println("GMT Time now: $hours:$minutes:$seconds")
end

# daytime()

# Exercise 5-3
function checkfermat(a, b, c, n)
    if n > 2 && a^n + b^n == c^n
        print("Holy smokes, Fermat was wrong!")
    else
        print("No that doesn't work.")
    end
end

# checkfermat(2, 3, 4, 10)

function checkfermat()
    println("Input value for a:")
    a = parse(Int64, readline())
    println("Input value for b:")
    b = parse(Int64, readline())
    println("Input value for c:")
    c = parse(Int64, readline())
    println("Input value for n:")
    n = parse(Int64, readline())
    checkfermat(a, b, c, n)
end

# checkfermat()

# Exercise 5-4
function istriangle(a, b, c)
    return a + b > c && a + c > b && b + c > a
end

# println(istriangle(5, 4, 3))
# println(istriangle(5, 7, 6))
# println(istriangle(11, 1, 3))

function istriangle()
    println("Input value for 1st side:")
    a = parse(Int64, readline())
    println("Input value for 2nd side:")
    b = parse(Int64, readline())
    println("Input value for 3rd side:")
    c = parse(Int64, readline())
    println("istriangle(a, b, c): $(istriangle(a, b, c))")
end

# istriangle()

# Exercise 5-5
function recurse(n, s)
    if n == 0
        println(s)
    else
        recurse(n-1, n+s)
    end
end

# This prints the sum of all numbers from 1 to n.
# recurse(10, 0)
# This results in a infinite recursion.
# recurse(-1, 0)

# Exercise 5-6
using ThinkJulia

function draw(t, length, n)
    if n == 0
        return
    end
    angle = 50
    forward(t, length*n)
    turn(t, -angle)
    draw(t, length, n-1)
    turn(t, 2*angle)
    draw(t, length, n-1)
    turn(t, -angle)
    forward(t, -length*n)
end

t = Turtle()
@svg begin
    draw(t, 10, 5)
end


text:
@inline function ExaTronKernels.dprsrch(n::Int, x::CuDeviceArray{Float64,1},
                         xl::CuDeviceArray{Float64,1},
                         xu::CuDeviceArray{Float64,1},
                         A::CuDeviceArray{Float64,2},
                         g::CuDeviceArray{Float64,1},
                         w::CuDeviceArray{Float64,1},
                         wa1::CuDeviceArray{Float64,1},
                         wa2::CuDeviceArray{Float64,1})
    one = 1.0
    p5 = 0.5

    # Constant that defines sufficient decrease.
    mu0 = 0.01

    # Interpolation factor.

    interpf = 0.5

    # Set the initial alpha = 1 because the quadratic function is
    # decreasing in the ray x + alpha*w for 0 <= alpha <= 1.

    alpha = one
    nsteps = 0

    # Find the smallest break-point on the ray x + alpha*w.
    nbrpt,brptmin,brptmax = dbreakpt(n,x,xl,xu,w)

    search = true
    while (search && alpha > brptmin)

        # Calculate P[x + alpha*w] - x and check the sufficient
        # decrease condition.

        nsteps = nsteps + 1
        dgpstep(n,x,xl,xu,alpha,w,wa1)
        dssyax(n,A,wa1,wa2)
        gts = ddot(n,g,1,wa1,1)
        q = p5*ddot(n,wa1,1,wa2,1) + gts
        if q <= mu0*gts
            search = false
        else

            # This is a crude interpolation procedure that
            # will be replaced in future versions of the code.

            alpha = interpf*alpha
        end
    end

    # Force at least one more constraint to be added to the active
    # set if alpha < brptmin and the full step is not successful.
    # There is sufficient decrease because the quadratic function
    # is decreasing in the ray x + alpha*w for 0 <= alpha <= 1.

    if (alpha < one && alpha < brptmin)
        alpha = brptmin
    end

    # Compute the final iterate and step.

    dgpstep(n,x,xl,xu,alpha,w,wa1)
    daxpy(n,alpha,w,1,x,1)
    dmid(n,x,xl,xu)
    dcopy(n,wa1,1,w,1)

    return
end