text:
#   Generating exports
export curl_strequal
export curl_strnequal
export curl_formadd
export curl_formget
export curl_formfree
export curl_getenv
export curl_version
export curl_easy_escape
export curl_escape
export curl_easy_unescape
export curl_unescape
export curl_free
export curl_global_init
export curl_global_init_mem
export curl_global_cleanup
export curl_slist_append
export curl_slist_free_all
export curl_getdate
export curl_share_init
export curl_share_setopt
export curl_share_cleanup
export curl_version_info
export curl_easy_strerror
export curl_share_strerror
export curl_easy_pause
export curl_easy_init
export curl_easy_setopt
export curl_easy_perform
export curl_easy_cleanup
export curl_easy_getinfo
export curl_easy_duphandle
export curl_easy_reset
export curl_easy_recv
export curl_easy_send
export curl_multi_init
export curl_multi_add_handle
export curl_multi_remove_handle
export curl_multi_fdset
export curl_multi_perform
export curl_multi_cleanup
export curl_multi_info_read
export curl_multi_strerror
export curl_multi_socket
export curl_multi_socket_action
export curl_multi_socket_all
export curl_multi_timeout
export curl_multi_setopt
export curl_multi_assign
export CURL
export curl_socket_t
export curl_progress_callback
export curl_write_callback
export CURLFILETYPE_FILE
export CURLFILETYPE_DIRECTORY
export CURLFILETYPE_SYMLINK
export CURLFILETYPE_DEVICE_BLOCK
export CURLFILETYPE_DEVICE_CHAR
export CURLFILETYPE_NAMEDPIPE
export CURLFILETYPE_SOCKET
export CURLFILETYPE_DOOR
export CURLFILETYPE_UNKNOWN
export curlfiletype
export curl_chunk_bgn_callback
export curl_chunk_end_callback
export curl_fnmatch_callback
export curl_seek_callback
export curl_read_callback
export CURLSOCKTYPE_IPCXN
export CURLSOCKTYPE_LAST
export curlsocktype
export curl_sockopt_callback
export curl_opensocket_callback
export curl_closesocket_callback
export CURLIOE_OK
export CURLIOE_UNKNOWNCMD
export CURLIOE_FAILRESTART
export CURLIOE_LAST
export curlioerr
export CURLIOCMD_NOP
export CURLIOCMD_RESTARTREAD
export CURLIOCMD_LAST
export curliocmd
export curl_ioctl_callback
export curl_malloc_callback
export curl_free_callback
export curl_realloc_callback
export curl_strdup_callback
export curl_calloc_callback
export CURLINFO_TEXT
export CURLINFO_HEADER_IN
export CURLINFO_HEADER_OUT
export CURLINFO_DATA_IN
export CURLINFO_DATA_OUT
export CURLINFO_SSL_DATA_IN
export CURLINFO_SSL_DATA_OUT
export CURLINFO_END
export curl_infotype
export curl_debug_callback
export CURLE_OK
export CURLE_UNSUPPORTED_PROTOCOL
export CURLE_FAILED_INIT
export CURLE_URL_MALFORMAT
export CURLE_NOT_BUILT_IN
export CURLE_COULDNT_RESOLVE_PROXY
export CURLE_COULDNT_RESOLVE_HOST
export CURLE_COULDNT_CONNECT
export CURLE_FTP_WEIRD_SERVER_REPLY
export CURLE_REMOTE_ACCESS_DENIED
export CURLE_FTP_ACCEPT_FAILED
export CURLE_FTP_WEIRD_PASS_REPLY
export CURLE_FTP_ACCEPT_TIMEOUT
export CURLE_FTP_WEIRD_PASV_REPLY
export CURLE_FTP_WEIRD_227_FORMAT
export CURLE_FTP_CANT_GET_HOST
export CURLE_OBSOLETE16
export CURLE_FTP_COULDNT_SET_TYPE
export CURLE_PARTIAL_FILE
export CURLE_FTP_COULDNT_RETR_FILE
export CURLE_OBSOLETE20
export CURLE_QUOTE_ERROR
export CURLE_HTTP_RETURNED_ERROR
export CURLE_WRITE_ERROR
export CURLE_OBSOLETE24
export CURLE_UPLOAD_FAILED
export CURLE_READ_ERROR
export CURLE_OUT_OF_MEMORY
export CURLE_OPERATION_TIMEDOUT
export CURLE_OBSOLETE29
export CURLE_FTP_PORT_FAILED
export CURLE_FTP_COULDNT_USE_REST
export CURLE_OBSOLETE32
export CURLE_RANGE_ERROR
export CURLE_HTTP_POST_ERROR
export CURLE_SSL_CONNECT_ERROR
export CURLE_BAD_DOWNLOAD_RESUME
export CURLE_FILE_COULDNT_READ_FILE
export CURLE_LDAP_CANNOT_BIND
export CURLE_LDAP_SEARCH_FAILED
export CURLE_OBSOLETE40
export CURLE_FUNCTION_NOT_FOUND
export CURLE_ABORTED_BY_CALLBACK
export CURLE_BAD_FUNCTION_ARGUMENT
export CURLE_OBSOLETE44
export CURLE_INTERFACE_FAILED
export CURLE_OBSOLETE46
export CURLE_TOO_MANY_REDIRECTS
export CURLE_UNKNOWN_OPTION
export CURLE_TELNET_OPTION_SYNTAX
export CURLE_OBSOLETE50
export CURLE_PEER_FAILED_VERIFICATION
export CURLE_GOT_NOTHING
export CURLE_SSL_ENGINE_NOTFOUND
export CURLE_SSL_ENGINE_SETFAILED
export CURLE_SEND_ERROR
export CURLE_RECV_ERROR
export CURLE_OBSOLETE57
export CURLE_SSL_CERTPROBLEM
export CURLE_SSL_CIPHER
export CURLE_SSL_CACERT
export CURLE_BAD_CONTENT_ENCODING
export CURLE_LDAP_INVALID_URL
export CURLE_FILESIZE_EXCEEDED
export CURLE_USE_SSL_FAILED
export CURLE_SEND_FAIL_REWIND
export CURLE_SSL_ENGINE_INITFAILED
export CURLE_LOGIN_DENIED
export CURLE_TFTP_NOTFOUND
export CURLE_TFTP_PERM
export CURLE_REMOTE_DISK_FULL
export CURLE_TFTP_ILLEGAL
export CURLE_TFTP_UNKNOWNID
export CURLE_REMOTE_FILE_EXISTS
export CURLE_TFTP_NOSUCHUSER
export CURLE_CONV_FAILED
export CURLE_CONV_REQD
export CURLE_SSL_CACERT_BADFILE
export CURLE_REMOTE_FILE_NOT_FOUND
export CURLE_SSH
export CURLE_SSL_SHUTDOWN_FAILED
export CURLE_AGAIN
export CURLE_SSL_CRL_BADFILE
export CURLE_SSL_ISSUER_ERROR
export CURLE_FTP_PRET_FAILED
export CURLE_RTSP_CSEQ_ERROR
export CURLE_RTSP_SESSION_ERROR
export CURLE_FTP_BAD_FILE_LIST
export CURLE_CHUNK_FAILED
export CURL_LAST
export CURLcode
export curl_conv_callback
export curl_ssl_ctx_callback
export CURLPROXY_HTTP
export CURLPROXY_HTTP_1_0
export CURLPROXY_SOCKS4
export CURLPROXY_SOCKS5
export CURLPROXY_SOCKS4A
export CURLPROXY_SOCKS5_HOSTNAME
export curl_proxytype
export CURLKHSTAT_FINE_ADD_TO_FILE
export CURLKHSTAT_FINE
export CURLKHSTAT_REJECT
export CURLKHSTAT_DEFER
export CURLKHSTAT_LAST
export CURLKHMATCH_OK
export CURLKHMATCH_MISMATCH
export CURLKHMATCH_MISSING
export CURLKHMATCH_LAST
export curl_sshkeycallback
export CURLUSESSL_NONE
export CURLUSESSL_TRY
export CURLUSESSL_CONTROL
export CURLUSESSL_ALL
export CURLUSESSL_LAST
export curl_usessl
export CURLFTPSSL_CCC_NONE
export CURLFTPSSL_CCC_PASSIVE
export CURLFTPSSL_CCC_ACTIVE
export CURLFTPSSL_CCC_LAST
export curl_ftpccc
export CURLFTPAUTH_DEFAULT
export CURLFTPAUTH_SSL
export CURLFTPAUTH_TLS
export CURLFTPAUTH_LAST
export curl_ftpauth
export CURLFTP_CREATE_DIR_NONE
export CURLFTP_CREATE_DIR
export CURLFTP_CREATE_DIR_RETRY
export CURLFTP_CREATE_DIR_LAST
export curl_ftpcreatedir
export CURLFTPMETHOD_DEFAULT
export CURLFTPMETHOD_MULTICWD
export CURLFTPMETHOD_NOCWD
export CURLFTPMETHOD_SINGLECWD
export CURLFTPMETHOD_LAST
export curl_ftpmethod
export CURLOPT_FILE
export CURLOPT_URL
export CURLOPT_PORT
export CURLOPT_PROXY
export CURLOPT_USERPWD
export CURLOPT_PROXYUSERPWD
export CURLOPT_RANGE
export CURLOPT_INFILE
export CURLOPT_ERRORBUFFER
export CURLOPT_WRITEFUNCTION
export CURLOPT_READFUNCTION
export CURLOPT_TIMEOUT
export CURLOPT_INFILESIZE
export CURLOPT_POSTFIELDS
export CURLOPT_REFERER
export CURLOPT_FTPPORT
export CURLOPT_USERAGENT
export CURLOPT_LOW_SPEED_LIMIT
export CURLOPT_LOW_SPEED_TIME
export CURLOPT_RESUME_FROM
export CURLOPT_COOKIE
export CURLOPT_HTTPHEADER
export CURLOPT_HTTPPOST
export CURLOPT_SSLCERT
export CURLOPT_KEYPASSWD
export CURLOPT_CRLF
export CURLOPT_QUOTE
export CURLOPT_WRITEHEADER
export CURLOPT_COOKIEFILE
export CURLOPT_SSLVERSION
export CURLOPT_TIMECONDITION
export CURLOPT_TIMEVALUE
export CURLOPT_CUSTOMREQUEST
export CURLOPT_STDERR
export CURLOPT_POSTQUOTE
export CURLOPT_WRITEINFO
export CURLOPT_VERBOSE
export CURLOPT_HEADER
export CURLOPT_NOPROGRESS
export CURLOPT_NOBODY
export CURLOPT_FAILONERROR
export CURLOPT_UPLOAD
export CURLOPT_POST
export CURLOPT_DIRLISTONLY
export CURLOPT_APPEND
export CURLOPT_NETRC
export CURLOPT_FOLLOWLOCATION
export CURLOPT_TRANSFERTEXT
export CURLOPT_PUT
export CURLOPT_PROGRESSFUNCTION
export CURLOPT_PROGRESSDATA
export CURLOPT_AUTOREFERER
export CURLOPT_PROXYPORT
export CURLOPT_POSTFIELDSIZE
export CURLOPT_HTTPPROXYTUNNEL
export CURLOPT_INTERFACE
export CURLOPT_KRBLEVEL
export CURLOPT_SSL_VERIFYPEER
export CURLOPT_CAINFO
export CURLOPT_MAXREDIRS
export CURLOPT_FILETIME
export CURLOPT_TELNETOPTIONS
export CURLOPT_MAXCONNECTS
export CURLOPT_CLOSEPOLICY
export CURLOPT_FRESH_CONNECT
export CURLOPT_FORBID_REUSE
export CURLOPT_RANDOM_FILE
export CURLOPT_EGDSOCKET
export CURLOPT_CONNECTTIMEOUT
export CURLOPT_HEADERFUNCTION
export CURLOPT_HTTPGET
export CURLOPT_SSL_VERIFYHOST
export CURLOPT_COOKIEJAR
export CURLOPT_SSL_CIPHER_LIST
export CURLOPT_HTTP_VERSION
export CURLOPT_FTP_USE_EPSV
export CURLOPT_SSLCERTTYPE
export CURLOPT_SSLKEY
export CURLOPT_SSLKEYTYPE
export CURLOPT_SSLENGINE
export CURLOPT_SSLENGINE_DEFAULT
export CURLOPT_DNS_USE_GLOBAL_CACHE
export CURLOPT_DNS_CACHE_TIMEOUT
export CURLOPT_PREQUOTE
export CURLOPT_DEBUGFUNCTION
export CURLOPT_DEBUGDATA
export CURLOPT_COOKIESESSION
export CURLOPT_CAPATH
export CURLOPT_BUFFERSIZE
export CURLOPT_NOSIGNAL
export CURLOPT_SHARE
export CURLOPT_PROXYTYPE
export CURLOPT_ACCEPT_ENCODING
export CURLOPT_PRIVATE
export CURLOPT_HTTP200ALIASES
export CURLOPT_UNRESTRICTED_AUTH
export CURLOPT_FTP_USE_EPRT
export CURLOPT_HTTPAUTH
export CURLOPT_SSL_CTX_FUNCTION
export CURLOPT_SSL_CTX_DATA
export CURLOPT_FTP_CREATE_MISSING_DIRS
export CURLOPT_PROXYAUTH
export CURLOPT_FTP_RESPONSE_TIMEOUT
export CURLOPT_IPRESOLVE
export CURLOPT_MAXFILESIZE
export CURLOPT_INFILESIZE_LARGE
export CURLOPT_RESUME_FROM_LARGE
export CURLOPT_MAXFILESIZE_LARGE
export CURLOPT_NETRC_FILE
export CURLOPT_USE_SSL
export CURLOPT_POSTFIELDSIZE_LARGE
export CURLOPT_TCP_NODELAY
export CURLOPT_FTPSSLAUTH
export CURLOPT_IOCTLFUNCTION
export CURLOPT_IOCTLDATA
export CURLOPT_FTP_ACCOUNT
export CURLOPT_COOKIELIST
export CURLOPT_IGNORE_CONTENT_LENGTH
export CURLOPT_FTP_SKIP_PASV_IP
export CURLOPT_FTP_FILEMETHOD
export CURLOPT_LOCALPORT
export CURLOPT_LOCALPORTRANGE
export CURLOPT_CONNECT_ONLY
export CURLOPT_CONV_FROM_NETWORK_FUNCTION
export CURLOPT_CONV_TO_NETWORK_FUNCTION
export CURLOPT_CONV_FROM_UTF8_FUNCTION
export CURLOPT_MAX_SEND_SPEED_LARGE
export CURLOPT_MAX_RECV_SPEED_LARGE
export CURLOPT_FTP_ALTERNATIVE_TO_USER
export CURLOPT_SOCKOPTFUNCTION
export CURLOPT_SOCKOPTDATA
export CURLOPT_SSL_SESSIONID_CACHE
export CURLOPT_SSH_AUTH_TYPES
export CURLOPT_SSH_PUBLIC_KEYFILE
export CURLOPT_SSH_PRIVATE_KEYFILE
export CURLOPT_FTP_SSL_CCC
export CURLOPT_TIMEOUT_MS
export CURLOPT_CONNECTTIMEOUT_MS
export CURLOPT_HTTP_TRANSFER_DECODING
export CURLOPT_HTTP_CONTENT_DECODING
export CURLOPT_NEW_FILE_PERMS
export CURLOPT_NEW_DIRECTORY_PERMS
export CURLOPT_POSTREDIR
export CURLOPT_SSH_HOST_PUBLIC_KEY_MD5
export CURLOPT_OPENSOCKETFUNCTION
export CURLOPT_OPENSOCKETDATA
export CURLOPT_COPYPOSTFIELDS
export CURLOPT_PROXY_TRANSFER_MODE
export CURLOPT_SEEKFUNCTION
export CURLOPT_SEEKDATA
export CURLOPT_CRLFILE
export CURLOPT_ISSUERCERT
export CURLOPT_ADDRESS_SCOPE
export CURLOPT_CERTINFO
export CURLOPT_USERNAME
export CURLOPT_PASSWORD
export CURLOPT_PROXYUSERNAME
export CURLOPT_PROXYPASSWORD
export CURLOPT_NOPROXY
export CURLOPT_TFTP_BLKSIZE
export CURLOPT_SOCKS5_GSSAPI_SERVICE
export CURLOPT_SOCKS5_GSSAPI_NEC
export CURLOPT_PROTOCOLS
export CURLOPT_REDIR_PROTOCOLS
export CURLOPT_SSH_KNOWNHOSTS
export CURLOPT_SSH_KEYFUNCTION
export CURLOPT_SSH_KEYDATA
export CURLOPT_MAIL_FROM
export CURLOPT_MAIL_RCPT
export CURLOPT_FTP_USE_PRET
export CURLOPT_RTSP_REQUEST
export CURLOPT_RTSP_SESSION_ID
export CURLOPT_RTSP_STREAM_URI
export CURLOPT_RTSP_TRANSPORT
export CURLOPT_RTSP_CLIENT_CSEQ
export CURLOPT_RTSP_SERVER_CSEQ
export CURLOPT_INTERLEAVEDATA
export CURLOPT_INTERLEAVEFUNCTION
export CURLOPT_WILDCARDMATCH
export CURLOPT_CHUNK_BGN_FUNCTION
export CURLOPT_CHUNK_END_FUNCTION
export CURLOPT_FNMATCH_FUNCTION
export CURLOPT_CHUNK_DATA
export CURLOPT_FNMATCH_DATA
export CURLOPT_RESOLVE
export CURLOPT_TLSAUTH_USERNAME
export CURLOPT_TLSAUTH_PASSWORD
export CURLOPT_TLSAUTH_TYPE
export CURLOPT_TRANSFER_ENCODING
export CURLOPT_CLOSESOCKETFUNCTION
export CURLOPT_CLOSESOCKETDATA
export CURLOPT_GSSAPI_DELEGATION
export CURLOPT_DNS_SERVERS
export CURLOPT_ACCEPTTIMEOUT_MS
export CURLOPT_TCP_KEEPALIVE
export CURLOPT_TCP_KEEPIDLE
export CURLOPT_TCP_KEEPINTVL
export CURLOPT_SSL_OPTIONS
export CURLOPT_MAIL_AUTH
export CURLOPT_LASTENTRY
export CURLoption
export CURL_HTTP_VERSION_NONE
export CURL_HTTP_VERSION_1_0
export CURL_HTTP_VERSION_1_1
export CURL_HTTP_VERSION_LAST
export CURL_NETRC_IGNORED
export CURL_NETRC_OPTIONAL
export CURL_NETRC_REQUIRED
export CURL_NETRC_LAST
export CURL_TLSAUTH_NONE
export CURL_TLSAUTH_SRP
export CURL_TLSAUTH_LAST
export CURL_TIMECOND_NONE
export CURL_TIMECOND_IFMODSINCE
export CURL_TIMECOND_IFUNMODSINCE
export CURL_TIMECOND_LASTMOD
export CURL_TIMECOND_LAST
export curl_TimeCond
export CURLFORM_NOTHING
export CURLFORM_COPYNAME
export CURLFORM_PTRNAME
export CURLFORM_NAMELENGTH
export CURLFORM_COPYCONTENTS
export CURLFORM_PTRCONTENTS
export CURLFORM_CONTENTSLENGTH
export CURLFORM_FILECONTENT
export CURLFORM_ARRAY
export CURLFORM_OBSOLETE
export CURLFORM_FILE
export CURLFORM_BUFFER
export CURLFORM_BUFFERPTR
export CURLFORM_BUFFERLENGTH
export CURLFORM_CONTENTTYPE
export CURLFORM_CONTENTHEADER
export CURLFORM_FILENAME
export CURLFORM_END
export CURLFORM_OBSOLETE2
export CURLFORM_STREAM
export CURLFORM_LASTENTRY
export CURLformoption
export CURL_FORMADD_OK
export CURL_FORMADD_MEMORY
export CURL_FORMADD_OPTION_TWICE
export CURL_FORMADD_NULL
export CURL_FORMADD_UNKNOWN_OPTION
export CURL_FORMADD_INCOMPLETE
export CURL_FORMADD_ILLEGAL_ARRAY
export CURL_FORMADD_DISABLED
export CURL_FORMADD_LAST
export CURLFORMcode
export curl_formget_callback
export CURLINFO_NONE
export CURLINFO_EFFECTIVE_URL
export CURLINFO_RESPONSE_CODE
export CURLINFO_TOTAL_TIME
export CURLINFO_NAMELOOKUP_TIME
export CURLINFO_CONNECT_TIME
export CURLINFO_PRETRANSFER_TIME
export CURLINFO_SIZE_UPLOAD
export CURLINFO_SIZE_DOWNLOAD
export CURLINFO_SPEED_DOWNLOAD
export CURLINFO_SPEED_UPLOAD
export CURLINFO_HEADER_SIZE
export CURLINFO_REQUEST_SIZE
export CURLINFO_SSL_VERIFYRESULT
export CURLINFO_FILETIME
export CURLINFO_CONTENT_LENGTH_DOWNLOAD
export CURLINFO_CONTENT_LENGTH_UPLOAD
export CURLINFO_STARTTRANSFER_TIME
export CURLINFO_CONTENT_TYPE
export CURLINFO_REDIRECT_TIME
export CURLINFO_REDIRECT_COUNT
export CURLINFO_PRIVATE
export CURLINFO_HTTP_CONNECTCODE
export CURLINFO_HTTPAUTH_AVAIL
export CURLINFO_PROXYAUTH_AVAIL
export CURLINFO_OS_ERRNO
export CURLINFO_NUM_CONNECTS
export CURLINFO_SSL_ENGINES
export CURLINFO_COOKIELIST
export CURLINFO_LASTSOCKET
export CURLINFO_FTP_ENTRY_PATH
export CURLINFO_REDIRECT_URL
export CURLINFO_PRIMARY_IP
export CURLINFO_APPCONNECT_TIME
export CURLINFO_CERTINFO
export CURLINFO_CONDITION_UNMET
export CURLINFO_RTSP_SESSION_ID
export CURLINFO_RTSP_CLIENT_CSEQ
export CURLINFO_RTSP_SERVER_CSEQ
export CURLINFO_RTSP_CSEQ_RECV
export CURLINFO_PRIMARY_PORT
export CURLINFO_LOCAL_IP
export CURLINFO_LOCAL_PORT
export CURLINFO_LASTONE
export CURLINFO
export CURLCLOSEPOLICY_NONE
export CURLCLOSEPOLICY_OLDEST
export CURLCLOSEPOLICY_LEAST_RECENTLY_USED
export CURLCLOSEPOLICY_LEAST_TRAFFIC
export CURLCLOSEPOLICY_SLOWEST
export CURLCLOSEPOLICY_CALLBACK
export CURLCLOSEPOLICY_LAST
export curl_closepolicy
export CURL_LOCK_DATA_NONE
export CURL_LOCK_DATA_SHARE
export CURL_LOCK_DATA_COOKIE
export CURL_LOCK_DATA_DNS
export CURL_LOCK_DATA_SSL_SESSION
export CURL_LOCK_DATA_CONNECT
export CURL_LOCK_DATA_LAST
export curl_lock_data
export CURL_LOCK_ACCESS_NONE
export CURL_LOCK_ACCESS_SHARED
export CURL_LOCK_ACCESS_SINGLE
export CURL_LOCK_ACCESS_LAST
export curl_lock_access
export curl_lock_function
export curl_unlock_function
export CURLSH
export CURLSHE_OK
export CURLSHE_BAD_OPTION
export CURLSHE_IN_USE
export CURLSHE_INVALID
export CURLSHE_NOMEM
export CURLSHE_NOT_BUILT_IN
export CURLSHE_LAST
export CURLSHcode
export CURLSHOPT_NONE
export CURLSHOPT_SHARE
export CURLSHOPT_UNSHARE
export CURLSHOPT_LOCKFUNC
export CURLSHOPT_UNLOCKFUNC
export CURLSHOPT_USERDATA
export CURLSHOPT_LAST
export CURLSHoption
export CURLVERSION_FIRST
export CURLVERSION_SECOND
export CURLVERSION_THIRD
export CURLVERSION_FOURTH
export CURLVERSION_LAST
export CURLversion
export CURLM
export CURLM_CALL_MULTI_PERFORM
export CURLM_OK
export CURLM_BAD_HANDLE
export CURLM_BAD_EASY_HANDLE
export CURLM_OUT_OF_MEMORY
export CURLM_INTERNAL_ERROR
export CURLM_BAD_SOCKET
export CURLM_UNKNOWN_OPTION
export CURLM_LAST
export CURLMcode
export CURLMSG_NONE
export CURLMSG_DONE
export CURLMSG_LAST
export CURLMSG
export curl_socket_callback
export curl_multi_timer_callback
export CURLMOPT_SOCKETFUNCTION
export CURLMOPT_SOCKETDATA
export CURLMOPT_PIPELINING
export CURLMOPT_TIMERFUNCTION
export CURLMOPT_TIMERDATA
export CURLMOPT_MAXCONNECTS
export CURLMOPT_LASTENTRY
export CURLMoption
export CURL_SOCKET_BAD
export CURL_MAX_WRITE_SIZE
export CURL_MAX_HTTP_HEADER
export CURL_WRITEFUNC_PAUSE
export CURLFINFOFLAG_KNOWN_FILENAME
export CURLFINFOFLAG_KNOWN_FILETYPE
export CURLFINFOFLAG_KNOWN_TIME
export CURLFINFOFLAG_KNOWN_PERM
export CURLFINFOFLAG_KNOWN_UID
export CURLFINFOFLAG_KNOWN_GID
export CURLFINFOFLAG_KNOWN_SIZE
export CURLFINFOFLAG_KNOWN_HLINKCOUNT
export CURL_CHUNK_BGN_FUNC_OK
export CURL_CHUNK_BGN_FUNC_FAIL
export CURL_CHUNK_BGN_FUNC_SKIP
export CURL_CHUNK_END_FUNC_OK
export CURL_CHUNK_END_FUNC_FAIL
export CURL_FNMATCHFUNC_MATCH
export CURL_FNMATCHFUNC_NOMATCH
export CURL_FNMATCHFUNC_FAIL
export CURL_SEEKFUNC_OK
export CURL_SEEKFUNC_FAIL
export CURL_SEEKFUNC_CANTSEEK
export CURL_READFUNC_ABORT
export CURL_READFUNC_PAUSE
export CURL_SOCKOPT_OK
export CURL_SOCKOPT_ERROR
export CURL_SOCKOPT_ALREADY_CONNECTED
export CURLE_OBSOLETE10
export CURLE_OBSOLETE12
export CURLOPT_ENCODING
export CURLE_UNKNOWN_TELNET_OPTION
export CURLE_SSL_PEER_CERTIFICATE
export CURLE_OBSOLETE
export CURLE_BAD_PASSWORD_ENTERED
export CURLE_BAD_CALLING_ORDER
export CURLE_FTP_USER_PASSWORD_INCORRECT
export CURLE_FTP_CANT_RECONNECT
export CURLE_FTP_COULDNT_GET_SIZE
export CURLE_FTP_COULDNT_SET_ASCII
export CURLE_FTP_WEIRD_USER_REPLY
export CURLE_FTP_WRITE_ERROR
export CURLE_LIBRARY_NOT_FOUND
export CURLE_MALFORMAT_USER
export CURLE_SHARE_IN_USE
export CURLE_URL_MALFORMAT_USER
export CURLE_FTP_ACCESS_DENIED
export CURLE_FTP_COULDNT_SET_BINARY
export CURLE_FTP_QUOTE_ERROR
export CURLE_TFTP_DISKFULL
export CURLE_TFTP_EXISTS
export CURLE_HTTP_RANGE_ERROR
export CURLE_FTP_SSL_FAILED
export CURLE_OPERATION_TIMEOUTED
export CURLE_HTTP_NOT_FOUND
export CURLE_HTTP_PORT_FAILED
export CURLE_FTP_COULDNT_STOR_FILE
export CURLE_FTP_PARTIAL_FILE
export CURLE_FTP_BAD_DOWNLOAD_RESUME
export CURLE_ALREADY_COMPLETE
export CURLAUTH_NONE
export CURLAUTH_BASIC
export CURLAUTH_DIGEST
export CURLAUTH_GSSNEGOTIATE
export CURLAUTH_NTLM
export CURLAUTH_DIGEST_IE
export CURLAUTH_NTLM_WB
export CURLAUTH_ONLY
export CURLAUTH_ANY
export CURLAUTH_ANYSAFE
export CURLSSH_AUTH_ANY
export CURLSSH_AUTH_NONE
export CURLSSH_AUTH_PUBLICKEY
export CURLSSH_AUTH_PASSWORD
export CURLSSH_AUTH_HOST
export CURLSSH_AUTH_KEYBOARD
export CURLSSH_AUTH_DEFAULT
export CURLGSSAPI_DELEGATION_NONE
export CURLGSSAPI_DELEGATION_POLICY_FLAG
export CURLGSSAPI_DELEGATION_FLAG
export CURL_ERROR_SIZE
export CURLSSLOPT_ALLOW_BEAST
export CURLFTPSSL_NONE
export CURLFTPSSL_TRY
export CURLFTPSSL_CONTROL
export CURLFTPSSL_ALL
export CURLFTPSSL_LAST
export CURLPROTO_HTTP
export CURLPROTO_HTTPS
export CURLPROTO_FTP
export CURLPROTO_FTPS
export CURLPROTO_SCP
export CURLPROTO_SFTP
export CURLPROTO_TELNET
export CURLPROTO_LDAP
export CURLPROTO_LDAPS
export CURLPROTO_DICT
export CURLPROTO_FILE
export CURLPROTO_TFTP
export CURLPROTO_IMAP
export CURLPROTO_IMAPS
export CURLPROTO_POP3
export CURLPROTO_POP3S
export CURLPROTO_SMTP
export CURLPROTO_SMTPS
export CURLPROTO_RTSP
export CURLPROTO_RTMP
export CURLPROTO_RTMPT
export CURLPROTO_RTMPE
export CURLPROTO_RTMPTE
export CURLPROTO_RTMPS
export CURLPROTO_RTMPTS
export CURLPROTO_GOPHER
export CURLPROTO_ALL
export CURLOPTTYPE_LONG
export CURLOPTTYPE_OBJECTPOINT
export CURLOPTTYPE_FUNCTIONPOINT
export CURLOPTTYPE_OFF_T
export CURLOPT_SERVER_RESPONSE_TIMEOUT
export CURLOPT_POST301
export CURLOPT_SSLKEYPASSWD
export CURLOPT_FTPAPPEND
export CURLOPT_FTPLISTONLY
export CURLOPT_FTP_SSL
export CURLOPT_SSLCERTPASSWD
export CURLOPT_KRB4LEVEL
export CURL_IPRESOLVE_WHATEVER
export CURL_IPRESOLVE_V4
export CURL_IPRESOLVE_V6
export CURLOPT_WRITEDATA
export CURLOPT_READDATA
export CURLOPT_HEADERDATA
export CURLOPT_RTSPHEADER
export CURL_REDIR_GET_ALL
export CURL_REDIR_POST_301
export CURL_REDIR_POST_302
export CURL_REDIR_POST_303
export CURL_REDIR_POST_ALL
export CURLINFO_STRING
export CURLINFO_LONG
export CURLINFO_DOUBLE
export CURLINFO_SLIST
export CURLINFO_MASK
export CURLINFO_TYPEMASK
export CURLINFO_HTTP_CODE
export CURL_GLOBAL_SSL
export CURL_GLOBAL_WIN32
export CURL_GLOBAL_ALL
export CURL_GLOBAL_NOTHING
export CURL_GLOBAL_DEFAULT
export CURLVERSION_NOW
export CURL_VERSION_IPV6
export CURL_VERSION_KERBEROS4
export CURL_VERSION_SSL
export CURL_VERSION_LIBZ
export CURL_VERSION_NTLM
export CURL_VERSION_GSSNEGOTIATE
export CURL_VERSION_DEBUG
export CURL_VERSION_ASYNCHDNS
export CURL_VERSION_SPNEGO
export CURL_VERSION_LARGEFILE
export CURL_VERSION_IDN
export CURL_VERSION_SSPI
export CURL_VERSION_CONV
export CURL_VERSION_CURLDEBUG
export CURL_VERSION_TLSAUTH_SRP
export CURL_VERSION_NTLM_WB
export CURLPAUSE_RECV
export CURLPAUSE_RECV_CONT
export CURLPAUSE_SEND
export CURLPAUSE_SEND_CONT
export CURLPAUSE_ALL
export CURLPAUSE_CONT
export CURLM_CALL_MULTI_SOCKET
export CURL_POLL_NONE
export CURL_POLL_IN
export CURL_POLL_OUT
export CURL_POLL_INOUT
export CURL_POLL_REMOVE
export CURL_SOCKET_TIMEOUT
export CURL_CSELECT_IN
export CURL_CSELECT_OUT
export CURL_CSELECT_ERR


text:
"""
    MinTimeZones <: AbstractSequencingAlgorithm

Algorithm that decides the sequence based on minimizing the time from the current
stop to the next stop.
"""
struct MinTimeZones <: AbstractSequencingAlgorithm end

function get_waiting_list(route, seq)
    zones_list = Vector{String}(undef, length(route.stops))
    i = 1
    for S_seq in seq.stops
        zones_list[i] = zone(route.stops[findfirst(x->x.id==S_seq.id, route.stops)])
        i += 1
    end
    zones_list_small = filter(x->x!="0", zones_list)
    zones_array = [(i,count(x->x==i,zones_list_small)) for i in unique(zones_list_small)]
    return zones_array
end

function filter_travel_times(route, T, stop, waiting_list, zone_num)
    stops_per_zone = name.(filter(x->zone(x)==waiting_list[zone_num][1], route.stops))
    println("STOPSPERZONE :"*string(stops_per_zone))
    println("STOP: "*string(stop))
    #is_stop_involved = (findfirst((x,_)->x==stop,stops_per_zone)!=nothing)
    is_stop_involved = (findfirst(x->x==stop,stops_per_zone)!=nothing)
    println("STOPINVOLVED: "*string(is_stop_involved))
    if(is_stop_involved)
        travel_times_out = stops_per_zone
    else
        current_stop_name = [stop]
        travel_times_out = [current_stop_name; stops_per_zone]
    end
    return T[travel_times_out, stops_per_zone]
end

function solve(alg::MinTimeZones, route::AbstractRoute, times::TravelTimes; seq::Sequence, initial_stop=station(route))
    T = times.travel_times
    return solve_zones(alg, route, T, seq, initial_stop=initial_stop)
end

function solve_zones(alg::MinTimeZones, route::AbstractRoute, T::AbstractMatrix, Q::Sequence;
               initial_stop=station(route))

    S = stops(route)
    num_stops = length(S)

    # salida
    route_id = name(route)
    seq = Vector{BareStop}()
    out = Sequence(route_id, seq)

    # orden de zonas
    waiting_list = get_waiting_list(route, Q)
    #println(waiting_list)
    stops_per_zone = 0
    zone_num = 1

    # encontrar la base y guardarla en el vector de salida
    x = initial_stop |> name
    push!(out, BareStop(x))

    isFirstStop = true

    T_mod = T

    # para cada parada, encontrar la siguiente con el menor tiempo
    while length(out) < num_stops
        if(isFirstStop)
            s = sort(view(T, x, :))
            i = 2
        else
            s = sort(view(T_mod, x, :))
            i = 1
        end

        y = s.dicts[1].keys[i]
        while (y ∈ out) # iterar hasta encontrar uno que ya visite
            i += 1
            i == length(s.dicts[1]) && return out
            if(i<=length(s.dicts[1].keys))
                y = s.dicts[1].keys[i]
            else
                break
            end
        end
        push!(out, BareStop(y))
        x = y
        stops_per_zone += 1
        if(stops_per_zone == waiting_list[zone_num][2])
            print(stops_per_zone)
            print(waiting_list[zone_num])
            zone_num += 1
            T_mod = filter_travel_times(route, T, x, waiting_list, zone_num)
            stops_per_zone = 0
            isFirstStop = false
        end
    end
    println("OUT:")
    println(out)
    return out
end

ext: jl


text:
# test Petsc matrix functions
function make_mat()
   A = PetscMat(PETSC_DECIDE, PETSC_DECIDE,  "mpiaij", comm, mlocal=sys_size_local, nlocal=sys_size_local)
  SetUp(A)
  low, high = MatGetOwnershipRange(A)
  global_indices = collect(low:PetscInt(high-1))

  MatSetValues(A, global_indices, global_indices, A_julia, INSERT_VALUES)
  MatAssemblyBegin(A)
  MatAssemblyEnd(A)

  return A
end

function make_vec()
  x = PetscVec(PETSC_DECIDE, VECMPI,  comm, mlocal=sys_size_local)

  low, high = VecGetOwnershipRange(x)
  global_indices = collect(low:PetscInt(high - 1))

  for i=1:sys_size_local
    idxm = [global_indices[i]]   # index
    val = [ rhs[i] ]  # value
    VecSetValues(x, idxm, val, INSERT_VALUES)
  end

  VecAssemblyBegin(x)
  VecAssemblyEnd(x)

  return x
end

function test_indexing()
  @testset "----- Testing Matrix indexing -----" begin

    A = make_mat()
    B = make_mat()

    x = make_vec()
    low, high = VecGetOwnershipRange(x)
    global_indices = collect(low:PetscInt(high - 1))
    xvec_julia = copy(rhs)

    low, high = MatGetOwnershipRange(A)
    global_indices = collect(low:PetscInt(high - 1))

    B_julia = A_julia + PetscScalar(1)

    # test set_values1!
    global_indices1 = global_indices + PetscInt(1)
    global_indices2 = global_indices + PetscInt(1)
    vals = rand(PetscScalar, sys_size_local, sys_size_local)

    set_values1!(A, global_indices1, global_indices2, vals)
    MatAssemblyBegin(A)
    MatAssemblyEnd(A)
    vals2 = zeros(vals)
    MatGetValues(A, global_indices, global_indices, vals2)

    @test isapprox( vals, vals2) atol=1e-13
    fill!(vals2, 0.0)
    get_values1!(A, global_indices1, global_indices2, vals2)
    @test isapprox( vals, vals2) atol=1e-13

    A2s = zeros(3, 3)
    idxm = collect(PetscInt, 1:3)
    idxn = collect(PetscInt, 1:3)
    vals = rand(3,3)
    set_values1!(A2s, idxm, idxn, vals)

    @test isapprox( A2s[idxm,  idxn], vals) atol=1e-13

    set_values1!(A2s, idxm, idxn, vals, ADD_VALUES)

    @test isapprox( A2s[idxm, idxn], 2*vals) atol=1e-13
    vals2 = zeros(vals)
    get_values1!(A2s, idxm, idxn, vals2)
    @test isapprox( vals2, 2*vals) atol=1e-13

    # test the sparse optimized version of +=
    A2ss = sparse(A2s)
    set_values1!(A2ss, idxm, idxn, vals, ADD_VALUES)
    fill!(vals2, 0.0)
    get_values1!(A2ss, idxm, idxn, vals2)
    @test isapprox( vals2, 3*vals) atol=1e-13

    for i=1:sys_size_local
      for j = 1:sys_size_local
        idxm = [ global_indices[i] ] # row index
        idxn = [ global_indices[j] ] # column index
        MatSetValues(A,idxm, idxn, [A_julia[i,j]],INSERT_VALUES);
        MatSetValues(B,idxm, idxn, [A_julia[i,j] + PetscScalar(1)],INSERT_VALUES);
      end
    end

    MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    PetscDestroy(A)
    PetscDestroy(B)
    PetscDestroy(x)

  end  # end facts block

  return nothing
end  # end function

function test_transpose()
  @testset "----- Testig Matrix Transpose -----" begin

    A = make_mat()
    # test CreateMatTransposed
    At = MatCreateTranspose(A)

    b2t = make_vec()
    low, high = VecGetOwnershipRange(b2t)
    global_indices = collect(low:PetscInt(high - 1))

    x2 = make_vec()

    MatMult(At, x2, b2t)

    b2t_julia = A_julia*rhs

    vals = zeros(PetscScalar, sys_size_local)
    VecGetValues(b2t, global_indices, vals)

    for i=1:sys_size_local
      @test isapprox( b2t_julia[i], vals[i]) atol=1e-12
    end

    PetscDestroy(At)
    PetscDestroy(b2t)
    PetscDestroy(x2)

    # test MatTranpose
    At = MatTranspose(A, inplace=false)
    vals = zeros(PetscScalar, sys_size_local, sys_size_local)
    valst = zeros(vals)
    MatGetValues(A, global_indices, global_indices, vals)
    MatGetValues(At, global_indices, global_indices, valst)

    for i=1:sys_size_local
      for j=1:sys_size_local
        @test isapprox( vals[i, j], valst[j, i]) atol=1e-12
      end
    end

    PetscDestroy(At)

    At = MatTranspose(A, inplace=true)
  #  fill!(vals, 0.0)
    fill!(valst, 0.0)
  #  MatGetValues(A, global_indices, global_indices, vals)
    MatGetValues(At, global_indices, global_indices, valst)

    for i=1:sys_size_local
      for j=1:sys_size_local
        @test isapprox( vals[i, j], valst[j, i]) atol=1e-12
      end
    end

    # transpose A back to normal
    MatTranspose(A, inplace=true)
  #  PetscDestroy(At)

    PetscDestroy(A)
    PetscDestroy(b2t)
    PetscDestroy(x2)

  end

  return nothing
end

#=
  for i=1:sys_size_local
    for j=1:sys_size_local
      idxm = [global_indices[i] ]  # row index
      idxn = [ global_indices[j] ] # column index
      v = zeros(PetscScalar, 1,1)
      MatGetValues(A, idxm, idxn, v)
#      println("i = ", i, " j = ", j, " v = ", v[1,1], " A[i,j] = ", A_julia[i,j])
      @test isapprox( v[1,1], A_julia[i,j])
      @test isapprox( A[i,j], A_julia[i,j])
    end
  end

  for i=1:sys_size_local*sys_size_local
    A[i] --> roughly(A_julia[i])
  end
=#

function test_blas()
@testset "----- Testing Matrix BLAS -----" begin
  A = make_mat()
  B = make_mat()
  C = make_mat()

  # transpose everything to be column major
  MatTranspose(A, inplace=true)
  MatTranspose(B, inplace=true)
  MatTranspose(C, inplace=true)

  B_julia = copy(A_julia)
  B_copy = zeros(PetscScalar, sys_size_local, sys_size_local)

  y = make_vec()
  y_julia = copy(rhs)

  x = make_vec()
  x_julia = copy(rhs)
  xvec_julia = copy(rhs)

  low, high = VecGetOwnershipRange(x)
  global_indices = collect(low:PetscInt(high - 1))

  @test ( size(A) )== (sys_size_local, sys_size_local)
  @test ( size(A, 1) )== sys_size_local
  @test ( size(A, 2) )== sys_size_local

  @test ( MatGetSize(A) )== (comm_size*sys_size_local, comm_size*sys_size_local)
  @test ( MatGetLocalSize(A) )== (sys_size_local, sys_size_local)
  # testing non zero exist code is all we can really do here
  @test ( PetscView(A) )== 0

  alpha = PetscScalar(2.3)

  MatAXPY(B, alpha, A, DIFFERENT_NONZERO_PATTERN)
  B_julia = alpha*A_julia + B_julia
  MatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size_local
    for j=1:sys_size_local
      @test isapprox( B_julia[j, i], B_copy[i, j])
    end
  end

  MatAYPX(B, alpha, A, DIFFERENT_NONZERO_PATTERN)
  B_julia = alpha*B_julia + A_julia
  MatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size_local
    for j=1:sys_size_local
      @test isapprox( B_julia[j, i], B_copy[i, j])
    end
  end

  MatScale(B, alpha)
  B_julia = alpha*B_julia
  MatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size_local
    for j=1:sys_size_local
      @test isapprox( B_julia[j, i], B_copy[i, j])
    end
  end

  fac = 2.0
  C_julia = rand(3,3)
  C_julia2 = fac*C_julia
  scale!(C_julia, 2.0)
  @test isapprox( C_julia, C_julia2)

  C_julia = sprand(10, 10, 0.1)
  C_julia2 = fac*C_julia
  scale!(C_julia, fac)
  @test isapprox( C_julia, C_julia2)

  MatShift(B, alpha)
  B_julia += alpha*eye(PetscScalar, sys_size_local)
  MatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size_local
    for j=1:sys_size_local
      @test isapprox( B_julia[j, i], B_copy[i, j])
    end
  end

  y_copy = zeros(PetscScalar, sys_size_local)
  MatMult(B, x, y)
  y_julia = B_julia*xvec_julia

  VecGetValues(y, sys_size_local, global_indices, y_copy)
  for i=1:sys_size_local
      @test isapprox( y_julia[i], y_copy[i])
  end

  MatMultAdd(B, x, y, y)
  y_julia = y_julia + B_julia*xvec_julia
  VecGetValues(y, sys_size_local, global_indices, y_copy)
  for i=1:sys_size_local
      @test isapprox( y_julia[i], y_copy[i])
  end

  MatMultTranspose(A, x, y)
  y_julia = A_julia.'*xvec_julia
  VecGetValues(y, sys_size_local, global_indices, y_copy)
  for i=1:sys_size_local
      @test isapprox( y_julia[i], y_copy[i])
  end

  MatMultHermitianTranspose(A, x, y)
  y_julia = A_julia'*xvec_julia
  VecGetValues(y, sys_size_local, global_indices, y_copy)
  for i=1:sys_size_local
      @test isapprox( y_julia[i], y_copy[i])
  end

  D = PetscMat(comm)
  MatMatMult(A, B, MAT_INITIAL_MATRIX, PetscReal(1.0), D)
  D_julia = A_julia*B_julia
  D_copy = zeros(PetscScalar, sys_size_local, sys_size_local)
  MatGetValues(D, global_indices, global_indices, D_copy)

  for i=1:sys_size_local
    for j=1:sys_size_local
      @test isapprox( D_copy[j,i], D_julia[i,j])
    end
  end

  # test matrix norms
  fnorm = MatNorm(D, NORM_FROBENIUS)
  infnorm = MatNorm(D, NORM_INFINITY)

  @test isapprox( fnorm, sqrt(comm_size)*vecnorm(D_julia)) atol=1e-2
  @test isapprox( infnorm, norm(D_julia, Inf))

  PetscDestroy(A)
  PetscDestroy(B)
  PetscDestroy(C)
  PetscDestroy(D)
  PetscDestroy(x)
  PetscDestroy(y)

end
return nothing
end

function test_prealloc()
@testset "----- Testing Matrix Preallocation -----" begin

  println("testing preallocation")
  nb = PetscInt(100)  # number of times/blocks to insert
  C = PetscMat(comm)
  MatSetType(C, "mpiaij")
# SetFromOptions(C)

  MatSetSizes(C, nb*sys_size_local, nb*sys_size_local, PetscInt(nb*comm_size*sys_size_local),PetscInt(nb*comm_size*sys_size_local));

  # preallocation parameters
  bs = PetscInt(1)
  dnnz = PetscInt(3)*ones(PetscInt, nb*sys_size_local)  # on diagonal (row + column owned by this process)
  onnz = zeros(PetscInt, nb*sys_size_local)  # no off diagonal (column not owned by this process)
  dnnzu = Array{PetscInt}(0)  # this is not a symmetric matrix, so unused
  onnzu = Array{PetscInt}(0)  # this is not a symmetric matrix, so unused

  MatXAIJSetPreallocation(C, bs, dnnz, onnz, dnnzu, onnzu)

#  SetUp(C);

  A_julia_t = A_julia.'  # transpose because C is row major

  low, high = MatGetOwnershipRange(C)
  idi = zeros(PetscInt, sys_size_local)  # row indices
  idj = zeros(PetscInt, sys_size_local)  # column indices

  for i=1:nb

    # get indices
    pos = 1
    for j=1:sys_size_local # insert block on diagonal
	idi[pos] = low + sys_size_local*(i - 1) + j - 1
	idj[pos] = low + sys_size_local*(i - 1) + j - 1
	pos += 1
    end

    # add a random component
    for i=1:sys_size_local
      for j=1:sys_size_local
	A_julia_t[i,j] += rand()
      end
    end

#    println("idi = ", idi)
#    println("idj = ", idj)
    MatSetValues(C, idi, idj, A_julia_t, INSERT_VALUES)
  end

  MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)
  MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)

  matinfo = MatGetInfo(C, MAT_LOCAL)

  @test isapprox( matinfo.mallocs, 0.0)

#  PetscView(C, 0)

  println("finished testing preallocation")

  MatZeroEntries(C)
  cnorm = MatNorm(C, NORM_FROBENIUS)
  @test isapprox( cnorm, 0.0) atol=1e-3# norm is zero iff matrix is zero

  @test ( PetscDestroy(C) )== 0

end
return nothing
end

test_indexing()
test_transpose()
test_blas()
test_prealloc()


text:
mutable struct iSplitUnit
    Channel::Int           # Channel Index
    RecordDate::String     # "180903"
    RecordTime::String     # "18-09-03 18:33:20"
    Comment::Dict{String, String} # Comments
    PhysicalUnit::Float32  # mV = PhysicalUnit * Value
    SamplingFreq::Float32   # sampling rate in hertz
    Value::Vector{Int16}   # raw value, shoule be integer.
end

mutable struct iSplitContainer
    Name::String             # file name or other label
    Comment::String          # Comments
    Data::Vector{iSplitUnit} # data list
end

# function chunk_iSplit(isplitdata::iSplitContainer, markers::Vector{T}, mbias::Float32, roi=Tuple{T, T}; flatten=true, selector::Function=x->true) where {T <: Number}
#     _targets = filter(selector, isplitdata.Data)
# end

"""load JLD file"""
function loadisplit(filename::String)
    _raw = load(filename)
    _data = [_raw[item] for item in keys(_raw) if !(item in ["Name", "Comment"])]

    iSplitContainer(_raw["Name"], _raw["Comment"], _data)
end


text:
module TestRaggedArrays

using FactCheck
using RaggedArrays

Adata = [1.0,2.0,
         3.0,4.0,5.0,
         6.0,
         7.0,8.0,9.0,10.0];
Aoffs = [0,2,5,6,10];
A = RaggedArray(Adata,Aoffs);

facts("Testing length(A)") do
    @fact length(A) => 4;
end

facts("Testing size(A)") do
    (sizes,n) = size(A);
    @fact n => 4;
    @fact sizes => [2,3,1,4];
end

facts("Testing index getters") do
    @fact A[2,1] => 2.0
    @fact_throws A[3,1];
    @fact_throws A[5,4];
    @fact_throws A[1,10];

end

facts("Testing index setters") do
    A[2,1] = 2.5;
    @fact A[2,1] => 2.5;
    A[2,1]=2.0;

    @fact_throws A[3,1]=1.0;
    @fact_throws A[5,4]=1.0;
    @fact_throws A[1,10]=1.0;
end

facts("Testing subarray lengths") do
    @fact getsubarraysize(A,1) => 2
    @fact getsubarraysize(A,2) => 3
    @fact getsubarraysize(A,3) => 1
    @fact getsubarraysize(A,4) => 4

    @fact_throws getsubarraysize(A,5)
end

facts("Testing subarray getter") do
    @fact getsubarray(A,1)=>[1.0,2.0];
    @fact getsubarray(A,2)=>[3.0,4.0,5.0];
    @fact getsubarray(A,3)=>[6.0];
    @fact getsubarray(A,4)=>[7.0,8.0,9.0,10.0];

    @fact_throws getsubarraysize(A,5);

end

facts("Testing sizes array constructor") do
    B=RaggedArray(Float64,[2,3,1,4]);
    @fact B.offs=>A.offs
end

facts("Testing regular array constructor") do
    nrows = 53;
    ncols = 101;
    A = rand(nrows,ncols);
    B = RaggedArray(A);
    b::Bool = true;
    for r = 1 : nrows
        for c = 1 : ncols
            b = b && (A[r,c]==B[r,c]);
        end
    end

    @fact b => true;
end

end


text:
# Run tests for the deflation

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra
using SparseArrays;
using BenchmarkTools;

#
#@show "nlevp_native"
#nep=nep_gallery("nlevp_native_gun");
#n=size(nep,1);
#(λ,v)=augnewton(nep,v=ones(n),λ=200^2,tol=1e-11)
#
#dnep1=deflate_eigpair(nep,λ,v,mode=:Generic);
#dnep2=deflate_eigpair(nep,λ,v,mode=:SPMF);
#dnep3=deflate_eigpair(nep,λ,v,mode=:MM);
#
#@btime augnewton(dnep1,v=ones(size(dnep1,1)),λ=250^2,
#                      tol=1e-11,maxit=300,armijo_factor=0.5)
#@btime augnewton(dnep2,v=ones(size(dnep1,1)),λ=250^2,
#                      tol=1e-11,maxit=300,armijo_factor=0.5)
#@btime augnewton(dnep3,v=ones(size(dnep1,1)),λ=complex(250^2),
#                      tol=1e-11,maxit=300,armijo_factor=0.5)
#asd()

nep=nep_gallery("dep0_sparse",100);
n=size(nep,1);
local λ,v;
(λ,v)=quasinewton(nep,v=ones(n),λ=2,tol=1e-11,logger=1,maxit=200)
dnep1=deflate_eigpair(nep,λ,v,mode=:Generic);
@btime mslp(dnep1,λ=1im,tol=1e-11)
dnep2=deflate_eigpair(nep,λ,v,mode=:SPMF);
@btime mslp(dnep2,λ=1im,tol=1e-11)
dnep3=deflate_eigpair(nep,λ,v,mode=:MM);
@btime mslp(dnep3,λ=1im,tol=1e-11)

asd()
begin
  nep=nep_gallery("dep0");
  n=size(nep,1);
  local λ,v;
  (λ,v)=quasinewton(nep,v=ones(n),λ=0,tol=1e-11)
  local dnep
  dnep=deflate_eigpair(nep,λ,v,mode=:Generic);
  for k=1:4
      (λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=-1+0.1im,
                      tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

      (λ2,V2)=get_deflated_eigpairs(dnep,λ,v)
      @show norm(compute_Mlincomb(nep,λ2[end],V2[:,end]))
      dnep=deflate_eigpair(dnep,λ,v);
  end
  (λ,V)=get_deflated_eigpairs(dnep)
  @show norm(compute_Mlincomb(nep,λ[end],V[:,end]))
  @show norm(compute_Mlincomb(nep,λ[1],V[:,1]))
end

begin
  nep=nep_gallery("nlevp_native_gun");
  n=size(nep,1);
  local λ,v;
  (λ,v)=augnewton(nep,v=ones(n),λ=200^2,tol=1e-11)
  local dnep
  dnep=deflate_eigpair(nep,λ,v,mode=:Generic);
  for k=1:2
      (λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=200^2,
                      tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

      (λ2,V2)=get_deflated_eigpairs(dnep,λ,v)
      @show norm(compute_Mlincomb(nep,λ2[end],V2[:,end]))
      dnep=deflate_eigpair(dnep,λ,v);
  end
  (λ,V)=get_deflated_eigpairs(dnep)
  @show norm(compute_Mlincomb(nep,λ[end],V[:,end]))
  @show norm(compute_Mlincomb(nep,λ[1],V[:,1]))
end

asd()

begin
  nep=nep_gallery("dep0");
  n=size(nep,1);
  local λ,v;
  (λ,v)=quasinewton(nep,v=ones(n),λ=0,tol=1e-11)
  local dnep
  dnep=deflate_eigpair(nep,λ,v,mode=:MM);
  for k=1:4
      (λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=-1+0.1im,
                      tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

      (λ2,V2)=get_deflated_eigpairs(dnep,λ,v)
      @show norm(compute_Mlincomb(nep,λ2[end],V2[:,end]))
      dnep=deflate_eigpair(dnep,λ,v);
  end
  (λ,V)=get_deflated_eigpairs(dnep)
  @show norm(compute_Mlincomb(nep,λ[end],V[:,end]))
  @show norm(compute_Mlincomb(nep,λ[1],V[:,1]))
end

#(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=0,
#                tol=1e-11,displaylevel=1,maxit=300,armijo_factor=0.5)
#dnep=deflate_eigpair(dnep,λ,v);
#(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=0.1im	,
#                tol=1e-11,displaylevel=1,maxit=300,armijo_factor=0.5)
#dnep=deflate_eigpair(dnep,λ,v);
#(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=0.1im	,
#                tol=1e-11,displaylevel=1,maxit=300,armijo_factor=0.5)
#dnep=deflate_eigpair(dnep,λ,v);

asd()

#nep=nep_gallery("nlevp_native_gun");
n=size(nep,1);
(λ,v)=quasinewton(nep,v=ones(n),λ=150^2,tol=1e-11)
dnep=deflate_eigpair(nep,λ,v,mode=:Generic);

asd()
(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=300^2,
                tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

dnep=deflate_eigpair(dnep,λ,v);
(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=300^2,
                tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)
dnep=deflate_eigpair(dnep,λ,v);
(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=300^2,
                tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

(λ,V)=get_deflated_eigpairs(dnep)

asd()

nep=nep_gallery("nlevp_native_gun");
n=size(nep,1);
(λ,v)=quasinewton(nep,v=ones(n),λ=150^2,tol=1e-11)

dnep=deflate(dnep,λ,v)

(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=300^2,
                tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)
dnep=deflate(dnep,λ,v)

(λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=300^2,
                tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

asd()

# Wanted usage:
# (λ,v)=quasinewton(nep)
# dnep=deflate(nep,λ,v)
# # dnep=deflate(nep,λ,v,mode=:SPMF); # Default
# # dnep=deflate(nep,λ,v,mode=:NaiveSPMF)
# # dnep=deflate(nep,λ,v,mode=:Generic)
# (λ,v)=quasinewton(dnep)
# dnep=deflate(dnep,λ,v)
# # dnep=deflate_schur(nep,S,V)
#
# # Extraction usage 1:
# (D,V)=deflated_eigpairs(dnep)
# # Extraction usage 2:
# (D,V)=deflated_eigpairs(dnep,λ,v)

# # Extraction Schur pairs
# (S,X)=deflated_schurfact(dnep)
# (S,X)=deflated_schurfact(dnep,λ,v)
#

nep=nep_gallery("nlevp_native_gun");
n=size(nep,1);
(λ,v)=quasinewton(nep,v=ones(n),λ=150^2,tol=1e-11)

#nep=nep_gallery("dep0_sparse");
#n=size(nep,1);
#(λ,v)=augnewton(nep,v=ones(n),λ=-0.4+0.3im,tol=1e-11)
v=v/norm(v);
S0=reshape([λ],1,1);
V0=reshape(v,n,1);
begin

    local S1=reshape([λ],1,1);
    local V1=reshape(v,n,1);

    @show 1
    for k=1:3
        normalize_schur_pair!(S1,V1);
        @show norm(compute_MM(nep,S1,V1))
        dnep=create_spmf_dnep(nep,S1,V1)
        #λ,V=get_deflated_eigpairs(dnep);
        @show norm(compute_Mlincomb(nep,λ[end],V[:,end]))
        (λ,v)=augnewton(dnep,v=ones(size(dnep,1)),λ=300^2,
                        tol=1e-11,logger=1,maxit=300,armijo_factor=0.5)

        V=V1;
        S=S1;
        V1=zeros(eltype(v),n,size(V,2)+1);
        S1=zeros(eltype(v),size(S,1)+1,size(S,2)+1);
        V1[1:n,1:end-1]=V[1:n,:];
        V1[1:n,end]=v[1:n];
        S1[1:end-1,1:end-1]=S;
        S1[1:end,end]=[v[n+1:end];λ];

        normalize_schur_pair!(S1,V1);

    end
    @show diag(S1)
    @show diag(V1'*V1)
end

asd()

# dnep1_eff=effenberger_deflation(nep,S0,V0)

dnep1_new=create_spmf_dnep(nep,S0,V0)

#W=randn(size(nep,1)+1,3);
##W[3:end].=0;
#Z=randn(3,3);
#@show norm(compute_MM(dnep,Z,W)-compute_MM(dnep2,Z,W))

(λ2,v2)=resinv(dnep1_new,λ=290^2,armijo_factor=0.9,logger=1,maxit=100,v=ones(n+1),tol=1e-12)
v2=v2/norm(v2);

#Z=randn(1,1); X=randn(n+1,1);
#compute_MM(dnep,Z,X)-compute_MM(dnep2,Z,X)

# Create the new invariant pair
V=V0;
S=S0;

V1=zeros(eltype(v2),n,size(V,2)+1);
S1=zeros(eltype(v2),size(S,1)+1,size(S,2)+1);
V1[1:n,1:end-1]=V[1:n,:];
V1[1:n,end]=v2[1:n];
S1[1:end-1,1:end-1]=S;
S1[1:end,end]=[v2[n+1:end];λ2];

normalize_schur_pair!(S1,V1);

@show norm(compute_MM(nep,S1,V1))

dnep2_new=create_spmf_dnep(nep,S1,V1)
#dnep2_eff=effenberger_deflation(nep,S1,V1);
#Z=randn(1,1); X=randn(n+2,1);
#compute_MM(dnep3b,Z,X)-compute_MM(dnep3,Z,X)

(λ3,v3)=quasinewton(dnep2_new,λ=300^2,armijo_factor=0.9,
                    logger=1,maxit=100,v=ones(n+2),tol=1e-12)

V=V1;
S=S1;

V1=zeros(ComplexF64,n,size(V,2)+1);
S1=zeros(ComplexF64,size(S,1)+1,size(S,2)+1);
V1[1:n,1:end-1]=V[1:n,:];
V1[1:n,end]=v3[1:n];
S1[1:end-1,1:end-1]=S;
S1[1:end,end]=[v3[(n+1):end];λ3];

@show norm(compute_MM(nep,S1,V1))

dnep3_new=create_spmf_dnep(nep,S1,V1)
#dnep2_eff=effenberger_deflation(nep,S1,V1);
#Z=randn(1,1); X=randn(n+2,1);
#compute_MM(dnep3b,Z,X)-compute_MM(dnep3,Z,X)

(λ4,v4)=quasinewton(dnep3_new,λ=295^2,armijo_factor=0.9,
                    logger=1,maxit=200,v=ones(n+3),tol=1e-12)

normalize!(v4)

V=V1;
S=S1;

V1=zeros(ComplexF64,n,size(V,2)+1);
S1=zeros(ComplexF64,size(S,1)+1,size(S,2)+1);
V1[1:n,1:end-1]=V[1:n,:];
V1[1:n,end]=v4[1:n];
S1[1:end-1,1:end-1]=S;
S1[1:end,end]=[v4[(n+1):end];λ4];

normalize_schur_pair!(S1,V1);

@show norm(compute_MM(nep,S1,V1))

dnep4_new=create_spmf_dnep(nep,S1,V1)

(λ5,v5)=quasinewton(dnep4_new,λ=295^2,armijo_factor=0.9,
                    logger=1,maxit=200,v=ones(n+4),tol=1e-12)

V=V1;
S=S1;

V1=zeros(ComplexF64,n,size(V,2)+1);
S1=zeros(ComplexF64,size(S,1)+1,size(S,2)+1);
V1[1:n,1:end-1]=V[1:n,:];
V1[1:n,end]=v5[1:n];
S1[1:end-1,1:end-1]=S;
S1[1:end,end]=[v5[(n+1):end];λ5];

normalize_schur_pair!(S1,V1);
#
#(QQ,RR)=qr(V1);
#V1=Matrix(QQ);
#S1=(RR*S1)/RR;
#

@show norm(compute_MM(nep,S1,V1))


text:
## DeltaF/F
function deltaFF(fs::AxisArray,Fo::Array{Float64,1},B=0.0)
    length(Fo) != size(fs)[2] ? error("Fo vector should be the same length as the number of ROIs in the series.") :
    Fo = permutedims(Fo)
    newF = (fs .- Fo)./(Fo .-B)
    AxisArray(newF,AxisArrays.axes(fs)...)
end


text:
include("../src/Transit.jl")

using Base.Test
import Transit

function square_trip(inval)
  io = IOBuffer()
  JSON.print(io, inval)
  seek(io, 0)
  Transit.parse(io)
end

@test square_trip(["~~foo"]) == ["~foo"]
@test square_trip(["~#'","~~foo"]) == "~foo"

@test square_trip(["~#'",1]) == 1
@test square_trip(Dict{Any,Any}("~#'" => 1)) == 1

@test square_trip([1, 2, 3, 4]) == [1, 2, 3, 4]
@test square_trip(["some", "funny", "words"]) == ["some", "funny", "words"]
@test square_trip([1, "and", 2, "we", Dict{Any,Any}("mix"=>"it up")]) ==
                 [1, "and", 2, "we", Dict{Any,Any}("mix"=>"it up")]
@test square_trip(Dict{Any,Any}("a" => 1, "b" => 2)) == Dict{Any,Any}("a" => 1, "b" => 2)
@test square_trip(Dict{Any,Any}("a" => 1)) == Dict{Any,Any}("a" => 1)
@test square_trip(Dict{Any,Any}("a" => Dict{Any,Any}("b" => 2))) ==
                 Dict{Any,Any}("a" => Dict{Any,Any}("b" => 2))
@test square_trip(["~:aaaa", "~:bbbb"]) == [:aaaa, :bbbb]
@test square_trip([Dict{Any,Any}("~:b" => "~i3"), "~i2"]) == [Dict{Any,Any}(:b => 3), 2]
@test square_trip(["~zINF", "~z-INF"]) == [Inf, -Inf]
@test isnan(square_trip(["~zNaN"])[1])
@test square_trip(["~f3.14"])[1] == Decimals.decimal("3.14")
@test square_trip(["~ude305d54-75b4-431b-adb2-eb6b9e546014"]) == [Base.Random.UUID("de305d54-75b4-431b-adb2-eb6b9e546014")]
@test square_trip(["~'ok"]) == ["ok"]
@test square_trip(["~_"]) == [nothing]
# fix for list change
#@test square_trip(Any["~#list",Any[1,2,"five"]]) == list(1, 2, "five")
@test square_trip(Any["~#set", Any[1, 2, 3]]) == Transit.TSet([1, 2, 3])
@test square_trip(Any["~#cmap",Any[Any[2,2],"two",Any[1,1],"one"]]) == Dict{Any, Any}(Any[2,2] => "two", Any[1,1] => "one")
# fix for URI change
#@test square_trip(Any["~rhttp://example.com","~rftp://example.com"]) == Any[URIParser.URI("http://example.com"), URIParser.URI("ftp://example.com")]
#@test square_trip(["~m-6106017600000","~m0","~m946728000000","~m1396909037000"]) == TBD - what do you typically want?


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

### Open Database

### Work database
"""
    workdb()

Return a [`FameDatabase`](@ref) instance for the work database.

Since the work database can be open only once, we create an instance the first
time and return the existing instance each subsequent time.
"""
function workdb end
export workdb

function cfmopwk()
    dbkey = Ref{Cint}(-1)
    @cfm_call_check(cfmopwk, (Ref{Cint},), dbkey)
    return dbkey[]
end

function workdb()
    global chli
    if !isopen(chli.workdb[])
        chli.workdb[] = FameDatabase(cfmopwk(), "WORK", :update)
    end
    return chli.workdb[]
end

function workdb(F::Function)
    db = workdb()
    try
        F(db)
    finally
        closedb!(db)
    end
end

### Open database

"""
    opendb(dbname, [mode])

Open a FAME database and return an instance of [`FameDatabase`](@ref) for it.

`dbname` can be a path to a .db file or a string specifying a database over a
remote connection in the following format.

    "[<tcp_port>@]<host> [<username> [<password>] ] <db>"

`mode` can be an integer (consult the CHLI help), or a `Symbol`. Valid modes
include `:readonly`, `:create`, `:overwrite`, `:update`, `:shared`, `:write`,
`:direct_write`.
"""
function opendb end
export opendb

function cfmopdb(dbname::String, mode::Int32)
    dbkey = Ref{Cint}(-1)
    @cfm_call_check(cfmopdb, (Ref{Cint}, Cstring, Cint), dbkey, dbname, mode)
    return dbkey[]
end

opendb(dbname::String, mode=:readonly) =
    FameDatabase(cfmopdb(dbname, val_to_int(mode, access_mode)), dbname, mode)

function opendb(F::Function, args...)
    db = opendb(args...)
    try
        F(db)
    finally
        closedb!(db)
    end
end

### Post database

export postdb
"""
    postdb(db::FameDatabase)

Post the given FAME database. If you've made any updates to the database you
must post it before closing, otherwise all your changes will be lost.
"""
postdb(db::FameDatabase) = @cfm_call_check(cfmpodb, (Cint,), db.key)

### Close database

export closedb!

cfmcldb(dbkey::Int32) = @cfm_call_check(cfmcldb, (Cint,), dbkey)

"""
    closedb!(db::FameDatabase)

Close the given FAME database.
"""
function closedb!(db::FameDatabase)
    cfmcldb(db.key)
    ind = indexin([db], chli.opendb)[1]
    if ind !== nothing
        deleteat!(chli.opendb, ind)
    else
        @warn "Closing untracked database. " (key = db.key, name = db.name, mode = db.mode)
    end
    db.key = -1
    return db
end


text:
struct SEChain{N} <: MorphologySE{N}
    data::Vector{<:AbstractArray{Bool}}
end

struct SEChainArray{N,AT<:AbstractArray{Bool,N}} <: MorphologySEArray{N}
    data::Vector{<:AbstractArray{Bool}}
    _data::AT # pre-calculated chained data as mask
    function SEChainArray{N}(data::Vector{<:AbstractArray{Bool}}) where {N}
        # TODO(johnnychen94): defer the calculation of chained data until it is needed since
        # this is mainly used for visualization purposes
        _data = _strel_chain(data)
        return new{N,typeof(_data)}(data, _data)
    end
end
function SEChainArray(data::Vector{<:AbstractArray{Bool}})
    data = convert(Vector, data)
    N = maximum(ndims.(data))
    return SEChainArray{N}(data)
end
SEChainArray(data::Tuple) = SEChainArray(collect(data))
strel_type(A::SEChainArray{N}) where {N} = SEChain{N}(A.data)
strel_size(A::SEChainArray) = size(A._data)

Base.axes(A::SEChainArray) = axes(A._data)
Base.size(A::SEChainArray) = size(A._data)
Base.@propagate_inbounds Base.getindex(A::SEChainArray, inds::Int...) = getindex(A._data, inds...)

"""
    strel_chain(A, B, ...)
    strel_chain(As)

For structuring elements of the same dimensions, chain them together to build a bigger one.

The output dimension is the same as the inputs dimensions. See also [`strel_product`](@ref)
that cartesian producting each SE.

!!! note "structuring element decomposition"
    For some morphological operations `f` such as dilation and erosion, if `se` can be
    decomposed into smaller ones `se1, se2, ..., seN`, then `f(img, se)` is equivalent to
    `f(...f(f(img, se1), se2), ..., seN)`. Because applying `f` to smaller SEs is more
    efficient than to the original big one, this trick is used widely in image morphology.

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> img = rand(512, 512);

julia> se1, se2 = [centered(rand(Bool, 3, 3)) for _ in 1:2];

julia> se = strel_chain(se1, se2);

julia> out_se = dilate(img, se);

julia> out_pipe = dilate(dilate(img, se1), se2);

julia> out_se[2:end-1, 2:end-1] == out_pipe[2:end-1, 2:end-1] # except for the boundary
true
```
"""
strel_chain(se, se_rest...) = strel_chain([se, se_rest...])
strel_chain(se_list::Vector{<:AbstractArray{T,N}}) where {T,N} = SEChainArray{N}(se_list)
strel_chain(se_list::Tuple) = SEChainArray(se_list)
strel_chain(se) = se

function _strel_chain(data::AbstractVector{<:AbstractArray})
    isempty(data) && throw(ArgumentError("data cannot be empty"))
    data = strel.(CartesianIndex, data)
    out = first(data)
    for i in axes(data, 1)[2:end]
        # TODO: preallocating the output can reduce a few more
        out = _simple_dilate(out, data[i])
    end
    out = strel(Bool, strel(CartesianIndex, out))
    return out
end

# a simple dilate function that automatically extends the boundary and dimension
function _simple_dilate(A::AbstractArray{T,N}, B::AbstractArray{T,N}) where {T,N}
    r = strel_size(B) .÷ 2
    sz = max.(strel_size(A), strel_size(B))
    out_sz = @. sz + 2r
    out = centered(falses(out_sz))

    R = strel(CartesianIndex, A)
    offsets = strel(CartesianIndex, B)
    i = zero(eltype(R))
    out[i] = true
    for o in offsets
        out[i + o] = true
    end
    for i in R
        out[i] = true
        for o in offsets
            out[i + o] = true
        end
    end

    # remove unnecessary zero boundaries
    return strel(CartesianIndex, out)
end

"""
    strel_product(A, B, ...)
    strel_product(se_list)

Cartesian product of multiple structuring elements; the output dimension `ndims(out) ==
sum(ndims, se_list)`.

See also [`strel_chain`](@ref) that chains SEs in the same dimension.

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> strel_product(strel_diamond((5, 5)), centered(Bool[1, 1, 1]))
5×5×3 SEChainArray{3, OffsetArrays.OffsetArray{Bool, 3, BitArray{3}}} with indices -2:2×-2:2×-1:1:
[:, :, -1] =
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0

[:, :, 0] =
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0

[:, :, 1] =
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0
```
"""
strel_product(se, se_rest...) = strel_product([se, se_rest...])

# TODO(johnnychen94): fix the type instability if it really matters in practice
function strel_product(se_list)
    N = sum(ndims, se_list)
    new_se_list = Array{Bool,N}[]
    ni = 0
    for se in se_list
        pre_ones = ntuple(_ -> one(Int), max(0, ni))
        post_ones = ntuple(_ -> one(Int), max(0, N - ndims(se) - ni))
        new_se = reshape(se, pre_ones..., size(se)..., post_ones...)
        push!(new_se_list, convert(Array{Bool,N}, new_se))
        ni += ndims(se)
    end
    return SEChainArray{N}(map(centered, new_se_list))
end


text:
module TestFunctions

using Bio.Seq

export get_bio_fmt_specimens,
    random_array,
    random_seq,
    random_dna,
    random_rna,
    random_aa,
    intempdir

function get_bio_fmt_specimens()
    path = Pkg.dir("Bio", "test", "BioFmtSpecimens")
    if !isdir(path)
        run(`git clone --depth 1 https://github.com/BioJulia/BioFmtSpecimens.git $(path)`)
    end
end

# The generation of random test cases...

function random_array(n::Integer, elements, probs)
    cumprobs = cumsum(probs)
    x = Array(eltype(elements), n)
    for i in 1:n
        x[i] = elements[searchsorted(cumprobs, rand()).start]
    end
    return x
end

# Return a random DNA/RNA sequence of the given length.
function random_seq(n::Integer, nts, probs)
    cumprobs = cumsum(probs)
    x = Array(Char, n)
    for i in 1:n
        x[i] = nts[searchsorted(cumprobs, rand()).start]
    end
    return convert(AbstractString, x)
end

function random_dna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end

function random_rna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end

function intempdir(fn::Function, parent=tempdir())
    dirname = mktempdir(parent)
    try
        cd(fn, dirname)
    finally
        rm(dirname, recursive=true)
    end
end

end


text:
#######################################################################
# Example: approximate density given by mixture model with a Gaussian #
#######################################################################

using PyPlot

# Define means for three-component Gaussian mixture model
# All components are implicitly equally weighted and have unit covariance
μ = [zeros(2), [2.9; 0.0], [2.5; 1.2]]

# Define log-likelihood
logp(θ) = log(exp(-0.5*sum((μ[1].-θ).^2)) + exp(-0.5*sum((μ[2].-θ).^2)) + exp(-0.5*sum((μ[3].-θ).^2)))

# Plot target density
figure()
θval = collect(LinRange(-5.0, 5.0, 100))
title("target unnormalised posterior")
contourf( repeat(θval,1,length(θval)),  repeat(θval',length(θval),1), map(θ -> exp(logp(θ)), [[x;y] for x in θval, y in θval]), cmap=plt.cm.binary)
axis("equal")

# Approximate mixture model with Gaussian density
posteriorfull,   =         VI(logp, randn(2); S = 150, iterations = 100)
posteriordiag,   =     VIdiag(logp, randn(2); S = 150, iterations = 100)
posteriorfixC,   = VIfixedcov(logp, randn(2), 1.0*Matrix(I, 2, 2); S = 150, iterations = 100)
posteriorlap,    =    laplace(logp, randn(2); iterations = 100)
posteriormvi,    =        MVI(logp, posteriorlap; S = 150, iterations = 100)
# posteriormvidiag,  = MVI(logp, MvNormal(zeros(2),Matrix(I,2,2)); S = 150, iterations = 100)

# Plot our approximations
ApproximateVI.plot_ellipse(posteriorfull,   "b", "full")
ApproximateVI.plot_ellipse(posteriordiag,   "r", "diag")
ApproximateVI.plot_ellipse(posteriorfixC,   "c", "sphere-fix")
ApproximateVI.plot_ellipse(posteriormvi,    "g", "mvi")
ApproximateVI.plot_ellipse(posteriorlap,    "m", "laplace")
# ApproximateVI.plot_ellipse(posteriormvidiag,  "k", "diag by mvi")

legend();

@info("Note how the covariance matrices of Laplace and MVI are rotated the same way")


text:
using JAMSDWriter, Compat, JuMP, MathProgBase
using Compat.Test, Compat.LinearAlgebra, Compat.SparseArrays

# This doesn't appears to be in Compat?
@static if VERSION >= v"0.7.0"
    function isdef(symbol)
        return isdefined(@__MODULE__, symbol)
    end
else
    function isdef(symbol)
        return isdefined(symbol)
    end
end

@testset "Basic Functionality" begin
    include("nl_convert.jl")
    include("nl_linearity.jl")
    include("nl_write.jl")
end

#solver = JuMP.UnsetSolver()
solvers = Any[]
push!(solvers, JAMSDSolver())

@testset "Examples" begin
    examples_path = joinpath(dirname(dirname(@__FILE__)), "examples")
    for example in ["jump_nltrig.jl",
                    "jump_nlexpr.jl",
                    "jump_minlp.jl",
                    "jump_nonlinearbinary.jl",
                    "jump_no_obj.jl",
                    "jump_const_obj.jl",
                    "jump_pruning.jl"
                    ]
        include(joinpath(examples_path, example))
    end
end
# need to support n-ary max and min operators
#for example in ["jump_maxmin.jl"]
#    include(joinpath(examples_path, example))
#end

# The following list is obtained though
#    grep "solvers = Any" JuMP/test/solvers.jl
lp_solvers = Any[]
ip_solvers = Any[]
ip_dual_solvers = Any[]
semi_solvers = Any[]
sos_solvers = Any[]
lazy_solvers, lazy_soc_solvers, lazylocal_solvers, cut_solvers, cutlocal_solvers, heur_solvers, info_solvers = Any[], Any[], Any[], Any[], Any[], Any[], Any[]
quad_solvers = Any[]
quad_soc_solvers = Any[]
rsoc_solvers = Any[]
nlp_solvers = Any[]
convex_nlp_solvers = Any[]
minlp_solvers = Any[]
sdp_solvers = Any[]

push!(lp_solvers, JAMSDSolver())
push!(lp_solvers, JAMSDSolver("cplex"))
#push!(lp_solvers, JAMSDSolver("gurobi"))
push!(lp_solvers, JAMSDSolver("mosek"))
push!(lp_solvers, JAMSDSolver("xpress"))

push!(quad_solvers, JAMSDSolver())
push!(quad_solvers, JAMSDSolver("cplex"))
#push!(quad_solvers, JAMSDSolver("gurobi"))
push!(quad_solvers, JAMSDSolver("mosek"))
push!(lp_solvers, JAMSDSolver("xpress"))

quad_mip_solvers = copy(quad_solvers)

# Investigate if we can solve those
#soc_solvers = copy(quad_solvers)
soc_solvers = Any[]

push!(convex_nlp_solvers, JAMSDSolver())
push!(convex_nlp_solvers, JAMSDSolver("conopt"))
# buggy on my system
#push!(convex_nlp_solvers, JAMSDSolver("ipopt"))
push!(convex_nlp_solvers, JAMSDSolver("minos"))
push!(convex_nlp_solvers, JAMSDSolver("snopt"))
push!(convex_nlp_solvers, JAMSDSolver("mosek"))

# push!(nlp_solvers, JAMSDSolver("baron"))
push!(nlp_solvers, JAMSDSolver("knitro"))

append!(convex_nlp_solvers, nlp_solvers)

push!(minlp_solvers, JAMSDSolver())

if VERSION < v"0.7"
    include(Pkg.dir("JuMP","test","qcqpmodel.jl"))
else
    include(joinpath(dirname(pathof(JuMP)), "..", "test", "qcqpmodel.jl"))
end

#
ipt = false

include("jump_test_nonlinear.jl")


text:
using LightXML
using Colors

docstr = read("basic.itermcolors", String)
xdoc = parse_string(docstr)
xroot = root(xdoc)
nodes = child_nodes(xroot)
top = collect(nodes)[2]

key = nothing
reorder = Dict()
for keydictnode in child_nodes(top)
    if name(keydictnode) == "key"
        global key
        key = content(keydictnode)
    elseif name(keydictnode) == "dict"
        keyname = nothing
        colors = []
        for dictnode in child_nodes(keydictnode)
            if name(dictnode) == "key"
                keyname = content(dictnode)
            elseif name(dictnode) == "real"
                if endswith(keyname, "Component")
                    color = parse(Float64, content(dictnode))
                    push!(colors, color)
                end
            end
        end
        if startswith(key, "Ansi")
            n = parse(Int, split(key, ' ')[2])
            color = hex(BGR(colors[2:end]...))
            reorder[n] = color
        end
    end
end

for n in sort(collect(keys(reorder)))
     color = reorder[n]
     println(".color", rpad(n, 2), " {fill: #", color, ";}")
end


text:
using LimitedLDLFactorizations
using Base.Test

# this matrix possesses an LDL factorization without pivoting
A = [ 1.7     0     0     0     0     0     0     0   .13     0
        0    1.     0     0   .02     0     0     0     0   .01
        0     0   1.5     0     0     0     0     0     0     0
        0     0     0   1.1     0     0     0     0     0     0
        0   .02     0     0   2.6     0   .16   .09   .52   .53
        0     0     0     0     0   1.2     0     0     0     0
        0     0     0     0   .16     0   1.3     0     0   .56
        0     0     0     0   .09     0     0   1.6   .11     0
      .13     0     0     0   .52     0     0   .11   1.4     0
        0   .01     0     0   .53     0   .56     0     0   3.1 ]
A = sparse(A)

L, d, α = lldl(A, memory=0)
nnzl0 = nnz(L)
@test nnzl0 == nnz(tril(A, -1))
@test α == 0

L, d, α = lldl(A, memory=5)
nnzl5 = nnz(L)
@test nnzl5 ≥ nnzl0
@test α == 0

L, d, α = lldl(A, memory=10)  # should be the exact factorization
@test nnz(L) ≥ nnzl5
@test α == 0
L = L + speye(10)
@test vecnorm(L * diagm(d) * L' - A) ≤ sqrt(eps()) * vecnorm(A)

# this matrix requires a shift
A = [ 1.  1.
      1.  0. ]
L, d, α = lldl(A)
@test α ≥ 1.0e-3
@test d[1] > 0
@test d[2] < 0

# specify our own shift
L, d, α = lldl(A, α=1.0e-2)
@test α ≥ 1.0e-2
@test d[1] > 0
@test d[2] < 0


text:
# shorter display of MersenneTwister:
Base.show(stream::IO, t::Random.MersenneTwister) =
    print(stream, "MersenneTwister($(t.seed))")

## WEIGHTED ENSEMBLES OF FITRESULTS

# Atom is atomic model type, eg, DecisionTree
# R will be the tightest type of the atom fit-results.
using StatsBase
mutable struct WrappedEnsemble{R,Atom <: Supervised} <: MLJType
    atom::Atom
    ensemble::Vector{R}
end

# Corner case: To address very special case of wrapped ensembles of
# categorical elements (eg, ensembles of fitresults for
# ConstantClassifier). These appear because doing comprehension with
# categorical elements gives CategoricalVector instead of Vector, but
# Vector is required in above struct definition.
function WrappedEnsemble(atom, ensemble::AbstractVector{L}) where L
    ensemble_vec = Vector{L}(undef, length(ensemble))
    for k in eachindex(ensemble)
        ensemble_vec[k] = ensemble[k]
    end
    return WrappedEnsemble(atom, ensemble_vec)
end

# to enable trait-based dispatch of predict:
predict(wens::WrappedEnsemble{R,Atom}, weights, Xnew) where {R,Atom<:Deterministic} =
    predict(wens, weights, Xnew, Deterministic, target_scitype(Atom))
predict(wens::WrappedEnsemble{R,Atom}, weights, Xnew) where {R,Atom<:Probabilistic} =
    predict(wens, weights, Xnew, Probabilistic, target_scitype(Atom))

function predict(wens::WrappedEnsemble,
                 weights,
                 Xnew,
                 ::Type{Deterministic}, ::Type{<:AbstractVector{<:Finite}})

    # weights ignored in this case

    ensemble = wens.ensemble
    atom = wens.atom

    n_atoms = length(ensemble)

    n_atoms > 0  || @error "Empty ensemble cannot make predictions."

    # TODO: make this more memory efficient but note that the type of
    # Xnew is unknown (ie, model dependent)
    predictions =
        reduce(hcat, [predict(atom, fitresult, Xnew) for fitresult in ensemble])
    classes = levels(predictions)
    n = size(predictions, 1)
    prediction =
        categorical(vcat([mode(predictions[i,:]) for i in 1:n], classes))[1:n]
    return prediction
end

function predict(wens::WrappedEnsemble, weights, Xnew, ::Type{Deterministic}, ::Type{<:AbstractVector{<:Continuous}})
    ensemble = wens.ensemble

    atom = wens.atom

    n_atoms = length(ensemble)

    n_atoms > 0  || @error "Empty ensemble cannot make predictions."

    # TODO: make more memory efficient:
    predictions = reduce(hcat, [weights[k]*predict(atom, ensemble[k], Xnew) for k in 1:n_atoms])
    prediction =  [sum(predictions[i,:]) for i in 1:size(predictions, 1)]

    return prediction
end

function predict(wens::WrappedEnsemble, weights, Xnew, ::Type{Probabilistic}, ::Type{<:AbstractVector{<:Finite}})

    ensemble = wens.ensemble

    atom = wens.atom

    n_atoms = length(ensemble)

    n_atoms > 0  || @error "Empty ensemble cannot make predictions."

    # TODO: make this more memory efficient but note that the type of
    # Xnew is unknown (ie, model dependent):

    # a matrix of probability distributions:
    predictions = reduce(hcat, [predict(atom, fitresult, Xnew) for fitresult in ensemble])
    n_rows = size(predictions, 1)

    # the weighted averages over the ensemble of the discrete pdf's:
    predictions  = [MLJBase.average([predictions[i,k] for k in 1:n_atoms], weights=weights) for i in 1:n_rows]

    return predictions
end

function predict(wens::WrappedEnsemble, weights, Xnew, ::Type{Probabilistic}, ::Type{<:AbstractVector{<:Continuous}})

    ensemble = wens.ensemble

    atom = wens.atom

    n_atoms = length(ensemble)

    n_atoms > 0  || @error "Empty ensemble cannot make predictions."

    # TODO: make this more memory efficient but note that the type of
    # Xnew is unknown (ie, model dependent):

    # a matrix of probability distributions:
    predictions = reduce(hcat, [predict(atom, fitresult, Xnew) for fitresult in ensemble])

    # n_rows = size(predictions, 1)
    # # the weighted average over the ensemble of the pdf means and pdf variances:
    # μs  = [sum([weights[k]*mean(predictions[i,k]) for k in 1:n_atoms]) for i in 1:n_rows]
    # σ2s = [sum([weights[k]*var(predictions[i,k]) for k in 1:n_atoms]) for i in 1:n_rows]

    # # a vector of normal probability distributions:
    # prediction = [Distributions.Normal(μs[i], sqrt(σ2s[i])) for i in 1:n_rows]

    prediction = [Distributions.MixtureModel(predictions[i,:], weights) for i in 1:size(predictions, 1)]

    return prediction

end

## CORE ENSEMBLE-BUILDING FUNCTIONS

# for when out-of-bag performance estimates are requested:
function get_ensemble_and_indices(atom::Supervised, verbosity, X, y, n, n_patterns,
                      n_train, rng, progress_meter)

    ensemble_indices =
        [StatsBase.sample(rng, 1:n_patterns, n_train, replace=false)
         for i in 1:n]
    ensemble = map(ensemble_indices) do train_rows
        verbosity == 1 && next!(progress_meter)
        verbosity < 2 ||  print("#")
        atom_fitresult, atom_cache, atom_report =
            fit(atom, verbosity - 1, selectrows(X, train_rows), selectrows(y, train_rows))
        atom_fitresult
    end
    verbosity < 1 || println()

    return (ensemble, ensemble_indices)

end

# for when out-of-bag performance estimates are not requested:
function get_ensemble(atom::Supervised, verbosity, X, y, n, n_patterns,
                      n_train, rng, progress_meter)

    # define generator of training rows:
    ensemble_indices = (StatsBase.sample(rng, 1:n_patterns, n_train, replace=false)
                        for i in 1:n)
    ensemble = map(ensemble_indices) do train_rows
        verbosity == 1 && next!(progress_meter)
        verbosity < 2 ||  print("#")
        atom_fitresult, atom_cache, atom_report =
            fit(atom, verbosity - 1, selectrows(X, train_rows), selectrows(y, train_rows))
        atom_fitresult
    end
    verbosity < 1 || println()

    return ensemble

end

pair_vcat(p, q) = (vcat(p[1], q[1]), vcat(p[2], q[2]))

## ENSEMBLE MODEL FOR DETERMINISTIC MODELS

mutable struct DeterministicEnsembleModel{Atom<:Deterministic} <: Deterministic
    atom::Atom
    weights::Vector{Float64}
    bagging_fraction::Float64
    rng::Union{Int,AbstractRNG}
    n::Int
    acceleration::AbstractResource
    out_of_bag_measure # TODO: type this
end

function clean!(model::DeterministicEnsembleModel)

    target_scitype(model.atom) <: Union{AbstractVector{<:Finite}, AbstractVector{<:Continuous}} ||
        error("`atom` has unsupported target_scitype "*
              "`$(target_scitype(model.atom))`. ")

    message = ""

    if model.bagging_fraction > 1 || model.bagging_fraction <= 0
        message = message*"`bagging_fraction` should be "*
        "in the range (0,1]. Reset to 1. "
        model.bagging_fraction = 1.0
    end

    if target_scitype(model.atom) <: AbstractVector{<:Finite} && !isempty(model.weights)
        message = message*"weights will be ignored to form predictions. "
    elseif !isempty(model.weights)
        total = sum(model.weights)
        if !(total ≈ 1.0)
            message = message*"weights should sum to one and are being automatically normalized. "
            model.weights = model.weights/total
        end
    end

    return message

end

# constructor to infer type automatically:
DeterministicEnsembleModel(atom::Atom, weights,
                           bagging_fraction, rng, n, acceleration, out_of_bag_measure) where Atom<:Deterministic =
                               DeterministicEnsembleModel{Atom}(atom, weights,
                                                                   bagging_fraction, rng, n, acceleration, out_of_bag_measure)

# lazy keyword constructors:
function DeterministicEnsembleModel(;atom=DeterministicConstantClassifier(),
                                    weights=Float64[],
                                    bagging_fraction=0.8,
                                    rng=Random.GLOBAL_RNG,
                                    n::Int=100,
                                    acceleration=DEFAULT_RESOURCE[],
                                    out_of_bag_measure=[])

    model = DeterministicEnsembleModel(atom, weights, bagging_fraction, rng,
                                       n, acceleration, out_of_bag_measure)

    message = clean!(model)
    isempty(message) || @warn message

    return model
end

## ENSEMBLE MODEL FOR PROBABILISTIC MODELS

mutable struct ProbabilisticEnsembleModel{Atom<:Probabilistic} <: Probabilistic
    atom::Atom
    weights::Vector{Float64}
    bagging_fraction::Float64
    rng::Union{Int, AbstractRNG}
    n::Int
    acceleration::AbstractResource
    out_of_bag_measure
end

function clean!(model::ProbabilisticEnsembleModel)

    message = ""

    if model.bagging_fraction > 1 || model.bagging_fraction <= 0
        message = message*"`bagging_fraction` should be "*
        "in the range (0,1]. Reset to 1. "
        model.bagging_fraction = 1.0
    end

    if !isempty(model.weights)
        total = sum(model.weights)
        if !(total ≈ 1.0)
            message = message*"Weights should sum to one and are being automatically normalized. "
            model.weights = model.weights/total
        end
    end

    return message

end

# constructor to infer type automatically:
ProbabilisticEnsembleModel(atom::Atom, weights, bagging_fraction, rng, n, acceleration, out_of_bag_measure) where Atom<:Probabilistic =
                               ProbabilisticEnsembleModel{Atom}(atom, weights, bagging_fraction, rng, n, acceleration, out_of_bag_measure)

# lazy keyword constructor:
function ProbabilisticEnsembleModel(;atom=ConstantProbabilisticClassifier(),
                                    weights=Float64[],
                                    bagging_fraction=0.8,
                                    rng=Random.GLOBAL_RNG,
                                    n::Int=100,
                                    acceleration=DEFAULT_RESOURCE[],
                                    out_of_bag_measure=[])

    model = ProbabilisticEnsembleModel(atom, weights, bagging_fraction, rng, n, acceleration, out_of_bag_measure)

    message = clean!(model)
    isempty(message) || @warn message

    return model
end

## COMMON CONSTRUCTOR

"""
    EnsembleModel(atom=nothing,
                  weights=Float64[],
                  bagging_fraction=0.8,
                  n=100,
                  rng=GLOBAL_RNG,
                  acceleration=DEFAULT_RESOURCE[],
                  out_of_bag_measure=[])

Create a model for training an ensemble of `n` learners, with optional
bagging, each with associated model `atom`. Ensembling is useful if
`fit!(machine(atom, data...))` does not create identical models on
repeated calls (ie, is a stochastic model, such as a decision tree
with randomized node selection criteria), or if `bagging_fraction` is
set to a value less than 1.0, or both. The constructor fails if no
`atom` is specified.

Only atomic models supporting targets with scitype
`AbstractVector{<:Finite}` (univariate classifiers) or
`AbstractVector{<:Continuous}` (univariate regressors) are supported.

If `rng` is an integer, then `MersenneTwister(rng)` is the random
number generator used for bagging. Otherwise some `AbstractRNG` object
is expected.

The atomic predictions are weighted according to the vector `weights` (to allow
for external optimization) except in the case that `atom` is a
`Deterministic` classifier. Uniform weights are used if `weight` has
zero length.

The ensemble model is `Deterministic` or `Probabilistic`, according to
the corresponding supertype of `atom`. In the case of deterministic
classifiers (`target_scitype(atom) <: Abstract{<:Finite}`), the
predictions are majority votes, and for regressors
(`target_scitype(atom)<: AbstractVector{<:Continuous}`) they are
ordinary averages.  Probabilistic predictions are obtained by
averaging the atomic probability distribution/mass functions; in
particular, for regressors, the ensemble prediction on each input
pattern has the type `MixtureModel{VF,VS,D}` from the Distributions.jl
package, where `D` is the type of predicted distribution for `atom`.

If a single measure or non-empty vector of measures is specified by
`out_of_bag_measure`, then out-of-bag estimates of performance are
reported.

The `acceleration` keyword argument is used to specify the compute resource (a
subtype of `ComputationalResources.AbstractResource`) that will be used to
accelerate/parallelize ensemble fitting.

"""
function EnsembleModel(; args...)
    d = Dict(args)
    :atom in keys(d) || error("No atomic model specified. Use EnsembleModel(atom=...)")
    if d[:atom] isa Deterministic
        return DeterministicEnsembleModel(; d...)
    elseif d[:atom] isa Probabilistic
        return ProbabilisticEnsembleModel(; d...)
    end
    error("$(d[:atom]) does not appear to be a Supervised model.")
end

## THE COMMON FIT AND PREDICT METHODS

const EitherEnsembleModel{Atom} = Union{DeterministicEnsembleModel{Atom}, ProbabilisticEnsembleModel{Atom}}

MLJBase.is_wrapper(::Type{<:EitherEnsembleModel}) = true

function _fit(res::CPU1, func, verbosity, args)
    atom, X, y, n, n_patterns, n_train, rng, progress_meter = args
    verbosity < 2 ||  @info "One hash per new atom trained: "
    return func(atom, verbosity, X, y,
                                    n, n_patterns, n_train, rng,
                                    progress_meter)
end
function _fit(res::CPUProcesses, func, verbosity, args)
    atom, X, y, n, n_patterns, n_train, rng, progress_meter = args
    if verbosity > 0
        println("Ensemble-building in parallel on $(nworkers()) processors.")
    end
    chunk_size = div(n, nworkers())
    left_over = mod(n, nworkers())
    return @distributed (pair_vcat) for i = 1:nworkers()
        if i != nworkers()
            func(atom, 0, X, y, chunk_size, n_patterns, n_train,
                         rng, progress_meter)
        else
            func(atom, 0, X, y, chunk_size + left_over, n_patterns, n_train,
                         rng, progress_meter)
        end
    end
end
@static if VERSION >= v"1.3.0-DEV.573"
    function _fit(res::CPUThreads, func, verbosity, args)
        atom, X, y, n, n_patterns, n_train, rng, progress_meter = args
        if verbosity > 0
            println("Ensemble-building in parallel on $(Threads.nthreads()) threads.")
        end
        nthreads = Threads.nthreads()
        chunk_size = div(n, nthreads)
        left_over = mod(n, nthreads)
        resvec = Vector(undef, nthreads) # FIXME: Make this type-stable?
        Threads.@threads for i = 1:nthreads
            resvec[i] = if i != nworkers()
                func(atom, 0, X, y, chunk_size, n_patterns, n_train,
                             rng, progress_meter)
            else
                func(atom, 0, X, y, chunk_size + left_over, n_patterns, n_train,
                             rng, progress_meter)
            end
        end
        return reduce(pair_vcat, resvec)
    end
end

function fit(model::EitherEnsembleModel{Atom}, verbosity::Int, X, y) where Atom<:Supervised

    acceleration = model.acceleration
    if acceleration isa CPUProcesses && nworkers() == 1
        acceleration = DEFAULT_RESOURCE[]
    end

    if model.out_of_bag_measure isa Vector
        out_of_bag_measure = model.out_of_bag_measure
    else
        out_of_bag_measure = [model.out_of_bag_measure,]
    end

    if model.rng isa Integer
        rng = MersenneTwister(model.rng)
    else
        rng = model.rng
    end

    atom = model.atom
    n = model.n
    n_patterns = nrows(X)
    n_train = round(Int, floor(model.bagging_fraction*n_patterns))

    progress_meter = Progress(n, dt=0.5, desc="Training ensemble: ",
                              barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

    if !isempty(out_of_bag_measure)

        args = atom, X, y, n, n_patterns, n_train, rng, progress_meter
        ensemble, ensemble_indices = _fit(acceleration, get_ensemble_and_indices, verbosity, args)

    else

        args = atom, X, y, n, n_patterns, n_train, rng, progress_meter
        ensemble = _fit(acceleration, get_ensemble, verbosity, args)

    end

    fitresult = WrappedEnsemble(model.atom, ensemble)

    if !isempty(out_of_bag_measure)

        metrics=zeros(length(ensemble),length(out_of_bag_measure))
        for i= 1:length(ensemble)
            #oob indices
            ooB_indices=  setdiff(1:n_patterns, ensemble_indices[i])
            if isempty(ooB_indices)
                error("Empty out-of-bag sample. "*
                      "Data size too small or "*
                      "bagging_fraction too close to 1.0. ")
            end
            predictions = predict(atom, ensemble[i], selectrows(X,ooB_indices))

            for k in eachindex(out_of_bag_measure)
                metrics[i,k] = out_of_bag_measure[k](predictions,selectrows(y, ooB_indices))
            end

        end
        metrics=mean(metrics, dims=1)

        names = Symbol.(string.(out_of_bag_measure))
        oob_estimates=NamedTuple{Tuple(names)}(Tuple(vec(metrics)))

    else
        oob_estimates=NamedTuple()
    end

    report=(oob_estimates=oob_estimates,)
    cache = deepcopy(model)

    return fitresult, cache, report

end

# if n is only parameter that changes, we just append to the existing
# ensemble, or truncate it:
function update(model::EitherEnsembleModel, verbosity::Int, fitresult, old_model, X, y)

    n = model.n

    if model.atom == old_model.atom &&
        model.bagging_fraction == old_model.bagging_fraction
        if n > old_model.n
            verbosity < 1 || @info "Building on existing ensemble of length $(old_model.n)"
            model.n = n - old_model.n # temporarily mutate the model
            wens, model_copy, report = fit(model, verbosity, X, y)
            append!(fitresult.ensemble, wens.ensemble)
            model.n = n         # restore model
            model_copy.n = n    # new copy of the model
        else
            verbosity < 1 || @info "Truncating existing ensemble."
            fitresult.ensemble = fitresult.ensemble[1:n]
            model_copy = deepcopy(model)
        end
        cache, report = model_copy, NamedTuple()
        return fitresult, cache, report
    else
        return fit(model, verbosity, X, y)
    end

end

function predict(model::EitherEnsembleModel, fitresult, Xnew)

    n = model.n
    if isempty(model.weights)
        weights = fill(1/n, n)
    else
        length(model.weights) == n || error("Ensemble size and number of weights not the same.")
        weights = model.weights
    end
    predict(fitresult, weights, Xnew)
end

## METADATA

# Note: input and target traits are inherited from atom

MLJBase.load_path(::Type{<:DeterministicEnsembleModel}) = "MLJ.DeterministicEnsembleModel"
MLJBase.package_name(::Type{<:DeterministicEnsembleModel}) = "MLJ"
MLJBase.package_uuid(::Type{<:DeterministicEnsembleModel}) = ""
MLJBase.package_url(::Type{<:DeterministicEnsembleModel}) = "https://github.com/alan-turing-institute/MLJ.jl"
MLJBase.is_pure_julia(::Type{<:DeterministicEnsembleModel{Atom}}) where Atom = MLJBase.is_pure_julia(Atom)
MLJBase.input_scitype(::Type{<:DeterministicEnsembleModel{Atom}}) where Atom = MLJBase.input_scitype(Atom)
MLJBase.target_scitype(::Type{<:DeterministicEnsembleModel{Atom}}) where Atom = MLJBase.target_scitype(Atom)

MLJBase.load_path(::Type{<:ProbabilisticEnsembleModel}) = "MLJ.ProbabilisticEnsembleModel"
MLJBase.package_name(::Type{<:ProbabilisticEnsembleModel}) = "MLJ"
MLJBase.package_uuid(::Type{<:ProbabilisticEnsembleModel}) = ""
MLJBase.package_url(::Type{<:ProbabilisticEnsembleModel}) = "https://github.com/alan-turing-institute/MLJ.jl"
MLJBase.is_pure_julia(::Type{<:ProbabilisticEnsembleModel{Atom}}) where Atom = MLJBase.is_pure_julia(Atom)
MLJBase.input_scitype(::Type{<:ProbabilisticEnsembleModel{Atom}}) where Atom = MLJBase.input_scitype(Atom)
MLJBase.target_scitype(::Type{<:ProbabilisticEnsembleModel{Atom}}) where Atom = MLJBase.target_scitype(Atom)

### old KoalaEnsembles code for optimizing the weights in the deterministic regressor case:

#     # Optimize weights:

#     n = length(ensemble)

#     if model.weight_regularization == 1
#         weights = ones(n)/n
#         verbosity < 1 || @info "Weighting atoms uniformly."
#     else
#         verbosity < 1 || print("\nOptimizing weights...")
#         Y = Array{Float64}(undef, n, n_patterns)
#         for k in 1:n
#             Y[k,:] = predict(model.atom, ensemble[k], X, false, false)
#         end

#         # If I rescale all predictions by the same amount it makes no
#         # difference to the values of the optimal weights:
#         ybar = mean(abs.(y))
#         Y = Y/ybar

#         A = Y*Y'
#         b = Y*(y/ybar)

#         scale = abs(det(A))^(1/n)

#         if scale < eps(Float64)

#             verbosity < 0 || @warn "Weight optimization problem ill-conditioned. " *
#                  "Using uniform weights."
#             weights = ones(n)/n

#         else

#             # need regularization, `gamma`, between 0 and infinity:
#             if model.weight_regularization == 0
#                 gamma = 0
#             else
#                 gamma = exp(atanh(2*model.weight_regularization - 1))
#             end

#             # add regularization and augment linear system for constraint
#             # (weights sum to one)
#             AA = hcat(A + scale*gamma*Matrix(I, n, n), ones(n))
#             AA = vcat(AA, vcat(ones(n), [0.0])')
#             bb = b + scale*gamma*ones(n)/n
#             bb = vcat(bb, [1.0])

#             weights = (AA \ bb)[1:n] # drop Lagrange multiplier
#             verbosity < 1 || println("\r$n weights optimized.\n")

#         end

#     end

#     fitresult = WrappedEnsemble(model.atom, ensemble, weights)
#     report = Dict{Symbol, Any}()
#     report[falsermalized_weights] = weights*length(weights)

#     cache = (X, y, scheme_X, ensemble)

#     return fitresult, report, cache

# end

# predict(model::ProbabilisticEnsembleModel, fitresult, Xt, parallel, verbosity) =
#     predict(fitresult, Xt)

# function fit_weights!(mach::SupervisedMachine{WrappedEnsemble{R, Atom},
#                                               ProbabilisticEnsembleModel{R, Atom}};
#               verbosity=1, parallel=true) where {R, Atom <: Supervised{R}}

#     mach.n_iter != 0 || @error "Cannot fit weights to empty ensemble."

#     mach.fitresult, report, mach.cache =
#         fit(mach.model, mach.cache, false, parallel, verbosity;
#             optimize_weights_only=true)
#     merge!(mach.report, report)

#     return mach
# end

# function weight_regularization_curve(mach::SupervisedMachine{WrappedEnsemble{R, Atom},
#                                                            ProbabilisticEnsembleModel{R, Atom}},
#                                      test_rows;
#                                      verbosity=1, parallel=true,
#                                      range=range(0, stop=1, length=101),
#                                      raw=false) where {R, Atom <: Supervised{R}}

#     mach.n_iter > 0 || @error "No atoms in the ensemble. Run `fit!` first."
#     !raw || verbosity < 0 ||
#         @warn "Reporting errors for *transformed* target. Use `raw=false` "*
#              " to report true errors."

#     if parallel && nworkers() > 1
#         if verbosity >= 1
#             println("Optimizing weights in parallel on $(nworkers()) processors.")
#         end
#         errors = pmap(range) do w
#             verbosity < 2 || print("\rweight_regularization=$w       ")
#             mach.model.weight_regularization = w
#             fit_weights!(mach; parallel=false, verbosity=verbosity - 1)
#             err(mach, test_rows, raw=raw)
#         end
#     else
#         errors = Float64[]
#         for w in range
#             verbosity < 1 || print("\rweight_regularization=$w       ")
#             mach.model.weight_regularization = w
#             fit_weights!(mach; parallel= parallel, verbosity=verbosity - 1)
#             push!(errors, err(mach, test_rows, raw=raw))
#         end
#         verbosity < 1 || println()

#         mach.report[:weight_regularization_curve] = (range, errors)
#     end

#     return range, errors
# end


text:
using Revise
using Random
using SpecialFunctions
using PyCall
using PyPlot
using StatsPlots

function mySliceSampler(pdf_log, x_0, w, m, N, burnIn = 1000)
    # REQUIRES MODULES:
    # Random
    #
    # INPUT:
    # pdf_log --> log-pdf of the target distribution (a function)
    # x_init --> inital vector of x (dimension 1xD), where D dimension of pdf
    # w --> typical slice size, for a window (dimension 1xD)
    # m --> integer limiting the slice size to m*w
    # N --> number of sample points
    # burnIn --> number of burn-in samples (optional, results will be discarded)
    #
    # OUTPUT:
    # x_s --> sampled x values (NxD), only keep values after full permutations along dimensions
    # pdflog_x_s --> log-pdf of sampled values

    D = length(x_0)  # the dimension of the distribution (nr of coordinates)
    x_s = Array{Float64}(undef,N,D)   # Samples (NxD), preallocate array, with undefined values
    pdflog_x_s = Array{Float64}(undef,N)   # (unnormalised) log-prob of sampled points (Nx1)
                                           # for reducing evaluations
    pdflog_x_1 = 0.0;

    # to be able to input scalars:
    if typeof(x_0) == Float64
        x_0 = [x_0]
    end

    for ii in 1:(N+burnIn)
        # Update the new x progressively for each dimension:
        L   = 1*x_0  # R and L need to be of same dimension, because we apply the pdf upon those!
        R   = 1*x_0
        x_1 = 1*x_0  # "1*" needed, otherwise x_0 suddendly takes mysterious values
        for d in randperm(D)
            ### 1) Make the Slice
            # random "vertical" position
            if (ii == 1)
                #vertical = rand() * pdf_log(x_0) # in normal space
                z = pdf_log(x_0) + log(rand()) # in log-space
            else
                # saving evaluations, pdflog_x_1 was updated for current permutation in the end
                z = pdflog_x_1 + log(rand())
            end

            ### 2) Make the window, "Stepping-out" alogrithm
            # Randomly place a window, containing x_0:
            L[d] = x_0[d] - rand() * w[d]
            R[d] = L[d] + w[d]

            # Randomly share the max window size among left/right:
            J = floor(m*rand());
            K = (m-1) - J;

            # Extend window to the left, until outside slice or allowance seizes:
            while ((J > 0) && (z < pdf_log(L)))
                #println("Lefting")
                L[d] -= w[d]
                J -= 1
            end

            # Extend window to the right, until outside slice or allowance seizes:
            while ((K > 0) && (z < pdf_log(R)))
                #println("Righting")
                R[d] += w[d]
                K -= 1
            end

            ### 3) Sample from window
            # finding an allowable point:
            while true
                x_1[d] = L[d] + rand() * (R[d] - L[d])

                # this + breaking out of loop reduces the amount of log-pdf evaluations:
                pdflog_x_1 = pdf_log(x_1)

                if (pdflog_x_1 >= z)     # new value found
                    x_0[d] = x_1[d]     # update value for this dimension
                    break
                # Value was not within slice, shrink the interval:
                elseif (x_1[d] < x_0[d])
                    L[d] = x_1[d]
                elseif (x_1[d] > x_0[d])
                    R[d] = x_1[d]
                else
                    throw(ErrorException("Error during shrinking the interval"))
                end
            end

        end # end of permutating through dimensions

        ### 4) Update the chain:
        # Just overwrite the burnIn by using modular arithmetic:
        i = 1 + (ii-1)%N    # gives 1,2,3,...N,1,2,3,...,(N+burnIn)
        x_s[i,:] = x_1    # sample point, after one round of permutations
        pdflog_x_s[i] = pdflog_x_1 # now contains log_pdf of all new coordinates of new point, since has been updated at each permutation
        # note that x_0 has been updated during the permutations

    end

    return x_s, pdflog_x_s
end

#%%
################################################################

function log_beta(θ,a::Any,b::Any)
    #using SpecialFunctions
    θ = θ[1]
    if 0 < θ < 1
        # Numerator
        num = log(θ)*(a-1) + log(1-θ)*(b-1)
        # The denominator (Beta-function) can be integrated or made up of Gamma-functions
        # Since we only treat integer values here, Γ(a) = (a-1)!
        #den = log(factorial(a-1) * factorial(b-1) / factorial(a+b-1))
        den = log(SpecialFunctions.beta(a,b))
        return num - den
    else
        return -Inf
    end
end

# Jespers HDI:
function hdi(theta_samp,alpha=0.05)
    cred_mass = 1.0-alpha
    ci = zeros(2)
    if length(size(theta_samp))>1   # multidimensional theta
        K,N = size(theta_samp)
        cis = zeros(2,K)
        for k in 1:K    # for each dimension

            ts = theta_samp[k,:]
            sind = sortperm(ts)     # get series of indices that sorts theta values
            sts = ts[sind]  # sorted theta values

            # Shifting a 95% index block from far left to far right and use the
            # samllest interval that can be found during the shift
            N = length(sind)    # number of sample points
            length_ci = Inf     # start with inf HDI interval and shrink
            for i in 1:Int(floor(N*alpha))  # iterate over the 5% lowest values (low to high)
                i2 = Int(floor(N*cred_mass)+i)  # the 5% highest values (low to high), one 95% block away on the other side
                prop_ci = [sts[i],sts[i2]]      # proposed interval
                length_prop_ci = prop_ci[2]-prop_ci[1]
                if length_prop_ci < length_ci
                    ci = prop_ci
                    length_ci = ci[2]-ci[1]
                end
            end
            cis[:,k] = ci

        end
        return cis
    else
        N = length(theta_samp)

        ts = theta_samp
        sind = sortperm(ts)
        sts = ts[sind]

        N = length(sind)
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i],sts[i2]]
            length_prop_ci = prop_ci[2]-prop_ci[1]
            if length_prop_ci < length_ci
                ci = prop_ci
                length_ci = ci[2]-ci[1]
            end
        end
        return ci
    end
end

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

# Log-Bernoulli PDF:
function log_bernoulli_pdf(θ,y)
    if 0 < θ < 1 && (y==0 || y==1)    # Only admissible values
        return y * log(θ) + (1-y) * log(1-θ)
    else
        return -Inf
    end
end

# Log-Likelihood fct of a (series of) coin tosses with bias:
function log_likelihood_fct(θ, y)
    # Make array s.t. multiple theta values can be input:
    #likelihood = Array{Float64}(undef,length(θ))
    # Do that for every theta value requested:
    #for ii in 1:length(θ)
    #    # Elementwise application of the log_beronoulli_pdf over outcome vector y, then SUM together:
    #    likelihood[ii] = sum(log_bernoulli_pdf.(θ[ii],y))
    #end

    # Un array-ise the output:
    #if length(θ) == 1
    #    likelihood = likelihood[1]
    #end

    θ = θ[1]
    likelihood = sum(log_bernoulli_pdf.(θ,y))

    return likelihood
end

#%%
################################################################

### Part A
# Chain inputs:
N = 10^6
x_start = 0.5;
w = .5;  # typical window width
m = 100; # maximum number of window widths

#%%
#First sampling
a = 1   # Flat Prior
b = 1
y = [ones(Int64,11); zeros(Int64,3)];   # Given observaitons

log_prior(θ) = log_beta(θ,a,b)
log_lklhd(θ) = log_likelihood_fct(θ,y)
log_posterior(θ) = log_lklhd(θ) + log_prior(θ)

# Need to use "StatsPlots" or "Plots" histogram to add the density curve, PyPlot won't work
chain, pdf = mySliceSampler(log_posterior,x_start,w,m,N);
chain = reshape(chain,:)
#%%
histogram(chain, bins=100, normalize=:pdf, label="posterior")   # Comes from StatsPlots now
density!(chain,linewidth=3,label="posterior")

#%%
# Sample Mean
μ = sum(chain)/N

# Mode
ω = chain[argmax(pdf)]

# Equal Tail Interval:
# Find the 2.5% highest and lowest (equal tails)
# get their index and receive an interval
α = 0.05;
idcesForOrder = sortperm(chain) # Indices of the chain which make it sorted
chain_sorted = chain[idcesForOrder]

leftTailIdx = Int(floor(N*(α/2)))
rightTailIdx = Int(ceil(N*(1- α/2)))

ETI = [chain_sorted[leftTailIdx],chain_sorted[rightTailIdx]]

HDI = hdi(chain)

#%%
p_gthalf = count(i->(i>0.5), chain)/N

#%%

### Part B

# Chain inputs:
N = 10^6
x_start = 0.5;
w = .5;  # typical window width
m = 100; # maximum number of window widths

#%%
#Second dataset
a = 1   # Flat Prior
b = 1
y = [ones(Int64,3); zeros(Int64,7)];   # Given observtions

log_prior(θ) = log_beta(θ,a,b)
log_lklhd(θ) = log_likelihood_fct(θ,y)
log_posterior(θ) = log_lklhd(θ) + log_prior(θ)

# Need to use "StatsPlots" or "Plots" histogram to add the density curve, PyPlot won't work
chain2, pdf2 = mySliceSampler(log_posterior,x_start,w,m,N);
chain2 = reshape(chain2,:);
#%%
histogram(chain2, bins=100, normalize=:pdf, label="posterior")   # Comes from StatsPlots now
density!(chain2,linewidth=3,label="posterior")

#%%  Method 1
gt = (chain .> chain2)*1
sum(gt)/N

#%% Method 1.1, questionable...
#gt = Array{Float64,1}(undef,N)
#for i in 1:length(chain)
    # for each entry in chain 1, compare it to ALL entries in chain2
#    gt[i] = sum((chain[i] .> chain2) * 1)/length(chain2)
#end

#sum(gt)/length(chain)

#%% Method 1.1, questionable...
#gt = (μ .> chain2) * 1
#sum(gt)/N

#%% Method 1.2, questionable...
#μ2 = sum(chain2)/N
#ω2 = chain[argmax(pdf)]

#idcesForOrder = sortperm(chain2) # Indices of the chain which make it sorted
#chain2_sorted = chain[idcesForOrder]

#gt = (chain_sorted .> maximum(chain2_sorted))*1
#sum(gt)

#%% Method 2
dθ = (chain - chain2)
gtZero = (dθ .> 0)*1
sum(gtZero)/N
hdi(dθ)

#%% Method 3
dθ = (chain - chain2)
dθ = reshape(dθ,:)      # to get an Array{Float64,1}
histogram(dθ, bins=100, normalize=:pdf, label="d_theta")
density!(dθ,linewidth=3,label="d_theta")
#%%
diff_hdi = hdi(dθ)
#%%
using Plots
Plots.savefig("/home/johhub/Desktop/ABDA/A4/A4-B-hist.pdf")