text:
macro c(ret_type, func, arg_types, lib)
  local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
  quote
    $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, lib)) ), $ret_type, $arg_types, $(args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  quote
    typealias $fake_t $real_t
  end
end

@ctypedef CURL None
@ctypedef curl_socket_t Int32
type curl_httppost
  next::Ptr{Void}
  name::Ptr{Uint8}
  namelength::Int32
  contents::Ptr{Uint8}
  contentslength::Int32
  buffer::Ptr{Uint8}
  bufferlength::Int32
  contenttype::Ptr{Uint8}
  contentheader::Ptr{Void}
  more::Ptr{Void}
  flags::Int32
  showfilename::Ptr{Uint8}
  userp::Ptr{None}
end
@ctypedef curl_progress_callback Ptr{Void}
@ctypedef curl_write_callback Ptr{Void}
# enum curlfiletype
const CURLFILETYPE_FILE = 0
const CURLFILETYPE_DIRECTORY = 1
const CURLFILETYPE_SYMLINK = 2
const CURLFILETYPE_DEVICE_BLOCK = 3
const CURLFILETYPE_DEVICE_CHAR = 4
const CURLFILETYPE_NAMEDPIPE = 5
const CURLFILETYPE_SOCKET = 6
const CURLFILETYPE_DOOR = 7
const CURLFILETYPE_UNKNOWN = 8
# end
@ctypedef curlfiletype Int32
type curl_fileinfo
  filename::Ptr{Uint8}
  filetype::curlfiletype
  time::time_t
  perm::Uint32
  uid::Int32
  gid::Int32
  size::curl_off_t
  hardlinks::Int32
  strings::Void
  flags::Uint32
  b_data::Ptr{Uint8}
  b_size::size_t
  b_used::size_t
end
@ctypedef curl_chunk_bgn_callback Ptr{Void}
@ctypedef curl_chunk_end_callback Ptr{Void}
@ctypedef curl_fnmatch_callback Ptr{Void}
@ctypedef curl_seek_callback Ptr{Void}
@ctypedef curl_read_callback Ptr{Void}
# enum curlsocktype
const CURLSOCKTYPE_IPCXN = 0
const CURLSOCKTYPE_LAST = 1
# end
@ctypedef curlsocktype Int32
@ctypedef curl_sockopt_callback Ptr{Void}
type curl_sockaddr
  family::Int32
  socktype::Int32
  protocol::Int32
  addrlen::Uint32
  addr::Void
end
@ctypedef curl_opensocket_callback Ptr{Void}
@ctypedef curl_closesocket_callback Ptr{Void}
# enum curlioerr
const CURLIOE_OK = 0
const CURLIOE_UNKNOWNCMD = 1
const CURLIOE_FAILRESTART = 2
const CURLIOE_LAST = 3
# end
@ctypedef curlioerr Int32
# enum curliocmd
const CURLIOCMD_NOP = 0
const CURLIOCMD_RESTARTREAD = 1
const CURLIOCMD_LAST = 2
# end
@ctypedef curliocmd Int32
@ctypedef curl_ioctl_callback Ptr{Void}
@ctypedef curl_malloc_callback Ptr{Void}
@ctypedef curl_free_callback Ptr{Void}
@ctypedef curl_realloc_callback Ptr{Void}
@ctypedef curl_strdup_callback Ptr{Void}
@ctypedef curl_calloc_callback Ptr{Void}
# enum curl_infotype
const CURLINFO_TEXT = 0
const CURLINFO_HEADER_IN = 1
const CURLINFO_HEADER_OUT = 2
const CURLINFO_DATA_IN = 3
const CURLINFO_DATA_OUT = 4
const CURLINFO_SSL_DATA_IN = 5
const CURLINFO_SSL_DATA_OUT = 6
const CURLINFO_END = 7
# end
@ctypedef curl_infotype Int32
@ctypedef curl_debug_callback Ptr{Void}
# enum CURLcode
const CURLE_OK = 0
const CURLE_UNSUPPORTED_PROTOCOL = 1
const CURLE_FAILED_INIT = 2
const CURLE_URL_MALFORMAT = 3
const CURLE_NOT_BUILT_IN = 4
const CURLE_COULDNT_RESOLVE_PROXY = 5
const CURLE_COULDNT_RESOLVE_HOST = 6
const CURLE_COULDNT_CONNECT = 7
const CURLE_FTP_WEIRD_SERVER_REPLY = 8
const CURLE_REMOTE_ACCESS_DENIED = 9
const CURLE_FTP_ACCEPT_FAILED = 10
const CURLE_FTP_WEIRD_PASS_REPLY = 11
const CURLE_FTP_ACCEPT_TIMEOUT = 12
const CURLE_FTP_WEIRD_PASV_REPLY = 13
const CURLE_FTP_WEIRD_227_FORMAT = 14
const CURLE_FTP_CANT_GET_HOST = 15
const CURLE_OBSOLETE16 = 16
const CURLE_FTP_COULDNT_SET_TYPE = 17
const CURLE_PARTIAL_FILE = 18
const CURLE_FTP_COULDNT_RETR_FILE = 19
const CURLE_OBSOLETE20 = 20
const CURLE_QUOTE_ERROR = 21
const CURLE_HTTP_RETURNED_ERROR = 22
const CURLE_WRITE_ERROR = 23
const CURLE_OBSOLETE24 = 24
const CURLE_UPLOAD_FAILED = 25
const CURLE_READ_ERROR = 26
const CURLE_OUT_OF_MEMORY = 27
const CURLE_OPERATION_TIMEDOUT = 28
const CURLE_OBSOLETE29 = 29
const CURLE_FTP_PORT_FAILED = 30
const CURLE_FTP_COULDNT_USE_REST = 31
const CURLE_OBSOLETE32 = 32
const CURLE_RANGE_ERROR = 33
const CURLE_HTTP_POST_ERROR = 34
const CURLE_SSL_CONNECT_ERROR = 35
const CURLE_BAD_DOWNLOAD_RESUME = 36
const CURLE_FILE_COULDNT_READ_FILE = 37
const CURLE_LDAP_CANNOT_BIND = 38
const CURLE_LDAP_SEARCH_FAILED = 39
const CURLE_OBSOLETE40 = 40
const CURLE_FUNCTION_NOT_FOUND = 41
const CURLE_ABORTED_BY_CALLBACK = 42
const CURLE_BAD_FUNCTION_ARGUMENT = 43
const CURLE_OBSOLETE44 = 44
const CURLE_INTERFACE_FAILED = 45
const CURLE_OBSOLETE46 = 46
const CURLE_TOO_MANY_REDIRECTS = 47
const CURLE_UNKNOWN_OPTION = 48
const CURLE_TELNET_OPTION_SYNTAX = 49
const CURLE_OBSOLETE50 = 50
const CURLE_PEER_FAILED_VERIFICATION = 51
const CURLE_GOT_NOTHING = 52
const CURLE_SSL_ENGINE_NOTFOUND = 53
const CURLE_SSL_ENGINE_SETFAILED = 54
const CURLE_SEND_ERROR = 55
const CURLE_RECV_ERROR = 56
const CURLE_OBSOLETE57 = 57
const CURLE_SSL_CERTPROBLEM = 58
const CURLE_SSL_CIPHER = 59
const CURLE_SSL_CACERT = 60
const CURLE_BAD_CONTENT_ENCODING = 61
const CURLE_LDAP_INVALID_URL = 62
const CURLE_FILESIZE_EXCEEDED = 63
const CURLE_USE_SSL_FAILED = 64
const CURLE_SEND_FAIL_REWIND = 65
const CURLE_SSL_ENGINE_INITFAILED = 66
const CURLE_LOGIN_DENIED = 67
const CURLE_TFTP_NOTFOUND = 68
const CURLE_TFTP_PERM = 69
const CURLE_REMOTE_DISK_FULL = 70
const CURLE_TFTP_ILLEGAL = 71
const CURLE_TFTP_UNKNOWNID = 72
const CURLE_REMOTE_FILE_EXISTS = 73
const CURLE_TFTP_NOSUCHUSER = 74
const CURLE_CONV_FAILED = 75
const CURLE_CONV_REQD = 76
const CURLE_SSL_CACERT_BADFILE = 77
const CURLE_REMOTE_FILE_NOT_FOUND = 78
const CURLE_SSH = 79
const CURLE_SSL_SHUTDOWN_FAILED = 80
const CURLE_AGAIN = 81
const CURLE_SSL_CRL_BADFILE = 82
const CURLE_SSL_ISSUER_ERROR = 83
const CURLE_FTP_PRET_FAILED = 84
const CURLE_RTSP_CSEQ_ERROR = 85
const CURLE_RTSP_SESSION_ERROR = 86
const CURLE_FTP_BAD_FILE_LIST = 87
const CURLE_CHUNK_FAILED = 88
const CURL_LAST = 89
# end
@ctypedef CURLcode Int32
@ctypedef curl_conv_callback Ptr{Void}
@ctypedef curl_ssl_ctx_callback Ptr{Void}
# enum curl_proxytype
const CURLPROXY_HTTP = 0
const CURLPROXY_HTTP_1_0 = 1
const CURLPROXY_SOCKS4 = 4
const CURLPROXY_SOCKS5 = 5
const CURLPROXY_SOCKS4A = 6
const CURLPROXY_SOCKS5_HOSTNAME = 7
# end
@ctypedef curl_proxytype Int32
type curl_khkey
  key::Ptr{Uint8}
  len::size_t
  keytype::Int32
end
# enum curl_khstat
const CURLKHSTAT_FINE_ADD_TO_FILE = 0
const CURLKHSTAT_FINE = 1
const CURLKHSTAT_REJECT = 2
const CURLKHSTAT_DEFER = 3
const CURLKHSTAT_LAST = 4
# end
# enum curl_khmatch
const CURLKHMATCH_OK = 0
const CURLKHMATCH_MISMATCH = 1
const CURLKHMATCH_MISSING = 2
const CURLKHMATCH_LAST = 3
# end
@ctypedef curl_sshkeycallback Ptr{Void}
# enum curl_usessl
const CURLUSESSL_NONE = 0
const CURLUSESSL_TRY = 1
const CURLUSESSL_CONTROL = 2
const CURLUSESSL_ALL = 3
const CURLUSESSL_LAST = 4
# end
@ctypedef curl_usessl Int32
# enum curl_ftpccc
const CURLFTPSSL_CCC_NONE = 0
const CURLFTPSSL_CCC_PASSIVE = 1
const CURLFTPSSL_CCC_ACTIVE = 2
const CURLFTPSSL_CCC_LAST = 3
# end
@ctypedef curl_ftpccc Int32
# enum curl_ftpauth
const CURLFTPAUTH_DEFAULT = 0
const CURLFTPAUTH_SSL = 1
const CURLFTPAUTH_TLS = 2
const CURLFTPAUTH_LAST = 3
# end
@ctypedef curl_ftpauth Int32
# enum curl_ftpcreatedir
const CURLFTP_CREATE_DIR_NONE = 0
const CURLFTP_CREATE_DIR = 1
const CURLFTP_CREATE_DIR_RETRY = 2
const CURLFTP_CREATE_DIR_LAST = 3
# end
@ctypedef curl_ftpcreatedir Int32
# enum curl_ftpmethod
const CURLFTPMETHOD_DEFAULT = 0
const CURLFTPMETHOD_MULTICWD = 1
const CURLFTPMETHOD_NOCWD = 2
const CURLFTPMETHOD_SINGLECWD = 3
const CURLFTPMETHOD_LAST = 4
# end
@ctypedef curl_ftpmethod Int32
# enum CURLoption
const CURLOPT_FILE = 10001
const CURLOPT_URL = 10002
const CURLOPT_PORT = 3
const CURLOPT_PROXY = 10004
const CURLOPT_USERPWD = 10005
const CURLOPT_PROXYUSERPWD = 10006
const CURLOPT_RANGE = 10007
const CURLOPT_INFILE = 10009
const CURLOPT_ERRORBUFFER = 10010
const CURLOPT_WRITEFUNCTION = 20011
const CURLOPT_READFUNCTION = 20012
const CURLOPT_TIMEOUT = 13
const CURLOPT_INFILESIZE = 14
const CURLOPT_POSTFIELDS = 10015
const CURLOPT_REFERER = 10016
const CURLOPT_FTPPORT = 10017
const CURLOPT_USERAGENT = 10018
const CURLOPT_LOW_SPEED_LIMIT = 19
const CURLOPT_LOW_SPEED_TIME = 20
const CURLOPT_RESUME_FROM = 21
const CURLOPT_COOKIE = 10022
const CURLOPT_HTTPHEADER = 10023
const CURLOPT_HTTPPOST = 10024
const CURLOPT_SSLCERT = 10025
const CURLOPT_KEYPASSWD = 10026
const CURLOPT_CRLF = 27
const CURLOPT_QUOTE = 10028
const CURLOPT_WRITEHEADER = 10029
const CURLOPT_COOKIEFILE = 10031
const CURLOPT_SSLVERSION = 32
const CURLOPT_TIMECONDITION = 33
const CURLOPT_TIMEVALUE = 34
const CURLOPT_CUSTOMREQUEST = 10036
const CURLOPT_STDERR = 10037
const CURLOPT_POSTQUOTE = 10039
const CURLOPT_WRITEINFO = 10040
const CURLOPT_VERBOSE = 41
const CURLOPT_HEADER = 42
const CURLOPT_NOPROGRESS = 43
const CURLOPT_NOBODY = 44
const CURLOPT_FAILONERROR = 45
const CURLOPT_UPLOAD = 46
const CURLOPT_POST = 47
const CURLOPT_DIRLISTONLY = 48
const CURLOPT_APPEND = 50
const CURLOPT_NETRC = 51
const CURLOPT_FOLLOWLOCATION = 52
const CURLOPT_TRANSFERTEXT = 53
const CURLOPT_PUT = 54
const CURLOPT_PROGRESSFUNCTION = 20056
const CURLOPT_PROGRESSDATA = 10057
const CURLOPT_AUTOREFERER = 58
const CURLOPT_PROXYPORT = 59
const CURLOPT_POSTFIELDSIZE = 60
const CURLOPT_HTTPPROXYTUNNEL = 61
const CURLOPT_INTERFACE = 10062
const CURLOPT_KRBLEVEL = 10063
const CURLOPT_SSL_VERIFYPEER = 64
const CURLOPT_CAINFO = 10065
const CURLOPT_MAXREDIRS = 68
const CURLOPT_FILETIME = 69
const CURLOPT_TELNETOPTIONS = 10070
const CURLOPT_MAXCONNECTS = 71
const CURLOPT_CLOSEPOLICY = 72
const CURLOPT_FRESH_CONNECT = 74
const CURLOPT_FORBID_REUSE = 75
const CURLOPT_RANDOM_FILE = 10076
const CURLOPT_EGDSOCKET = 10077
const CURLOPT_CONNECTTIMEOUT = 78
const CURLOPT_HEADERFUNCTION = 20079
const CURLOPT_HTTPGET = 80
const CURLOPT_SSL_VERIFYHOST = 81
const CURLOPT_COOKIEJAR = 10082
const CURLOPT_SSL_CIPHER_LIST = 10083
const CURLOPT_HTTP_VERSION = 84
const CURLOPT_FTP_USE_EPSV = 85
const CURLOPT_SSLCERTTYPE = 10086
const CURLOPT_SSLKEY = 10087
const CURLOPT_SSLKEYTYPE = 10088
const CURLOPT_SSLENGINE = 10089
const CURLOPT_SSLENGINE_DEFAULT = 90
const CURLOPT_DNS_USE_GLOBAL_CACHE = 91
const CURLOPT_DNS_CACHE_TIMEOUT = 92
const CURLOPT_PREQUOTE = 10093
const CURLOPT_DEBUGFUNCTION = 20094
const CURLOPT_DEBUGDATA = 10095
const CURLOPT_COOKIESESSION = 96
const CURLOPT_CAPATH = 10097
const CURLOPT_BUFFERSIZE = 98
const CURLOPT_NOSIGNAL = 99
const CURLOPT_SHARE = 10100
const CURLOPT_PROXYTYPE = 101
const CURLOPT_ACCEPT_ENCODING = 10102
const CURLOPT_PRIVATE = 10103
const CURLOPT_HTTP200ALIASES = 10104
const CURLOPT_UNRESTRICTED_AUTH = 105
const CURLOPT_FTP_USE_EPRT = 106
const CURLOPT_HTTPAUTH = 107
const CURLOPT_SSL_CTX_FUNCTION = 20108
const CURLOPT_SSL_CTX_DATA = 10109
const CURLOPT_FTP_CREATE_MISSING_DIRS = 110
const CURLOPT_PROXYAUTH = 111
const CURLOPT_FTP_RESPONSE_TIMEOUT = 112
const CURLOPT_IPRESOLVE = 113
const CURLOPT_MAXFILESIZE = 114
const CURLOPT_INFILESIZE_LARGE = 30115
const CURLOPT_RESUME_FROM_LARGE = 30116
const CURLOPT_MAXFILESIZE_LARGE = 30117
const CURLOPT_NETRC_FILE = 10118
const CURLOPT_USE_SSL = 119
const CURLOPT_POSTFIELDSIZE_LARGE = 30120
const CURLOPT_TCP_NODELAY = 121
const CURLOPT_FTPSSLAUTH = 129
const CURLOPT_IOCTLFUNCTION = 20130
const CURLOPT_IOCTLDATA = 10131
const CURLOPT_FTP_ACCOUNT = 10134
const CURLOPT_COOKIELIST = 10135
const CURLOPT_IGNORE_CONTENT_LENGTH = 136
const CURLOPT_FTP_SKIP_PASV_IP = 137
const CURLOPT_FTP_FILEMETHOD = 138
const CURLOPT_LOCALPORT = 139
const CURLOPT_LOCALPORTRANGE = 140
const CURLOPT_CONNECT_ONLY = 141
const CURLOPT_CONV_FROM_NETWORK_FUNCTION = 20142
const CURLOPT_CONV_TO_NETWORK_FUNCTION = 20143
const CURLOPT_CONV_FROM_UTF8_FUNCTION = 20144
const CURLOPT_MAX_SEND_SPEED_LARGE = 30145
const CURLOPT_MAX_RECV_SPEED_LARGE = 30146
const CURLOPT_FTP_ALTERNATIVE_TO_USER = 10147
const CURLOPT_SOCKOPTFUNCTION = 20148
const CURLOPT_SOCKOPTDATA = 10149
const CURLOPT_SSL_SESSIONID_CACHE = 150
const CURLOPT_SSH_AUTH_TYPES = 151
const CURLOPT_SSH_PUBLIC_KEYFILE = 10152
const CURLOPT_SSH_PRIVATE_KEYFILE = 10153
const CURLOPT_FTP_SSL_CCC = 154
const CURLOPT_TIMEOUT_MS = 155
const CURLOPT_CONNECTTIMEOUT_MS = 156
const CURLOPT_HTTP_TRANSFER_DECODING = 157
const CURLOPT_HTTP_CONTENT_DECODING = 158
const CURLOPT_NEW_FILE_PERMS = 159
const CURLOPT_NEW_DIRECTORY_PERMS = 160
const CURLOPT_POSTREDIR = 161
const CURLOPT_SSH_HOST_PUBLIC_KEY_MD5 = 10162
const CURLOPT_OPENSOCKETFUNCTION = 20163
const CURLOPT_OPENSOCKETDATA = 10164
const CURLOPT_COPYPOSTFIELDS = 10165
const CURLOPT_PROXY_TRANSFER_MODE = 166
const CURLOPT_SEEKFUNCTION = 20167
const CURLOPT_SEEKDATA = 10168
const CURLOPT_CRLFILE = 10169
const CURLOPT_ISSUERCERT = 10170
const CURLOPT_ADDRESS_SCOPE = 171
const CURLOPT_CERTINFO = 172
const CURLOPT_USERNAME = 10173
const CURLOPT_PASSWORD = 10174
const CURLOPT_PROXYUSERNAME = 10175
const CURLOPT_PROXYPASSWORD = 10176
const CURLOPT_NOPROXY = 10177
const CURLOPT_TFTP_BLKSIZE = 178
const CURLOPT_SOCKS5_GSSAPI_SERVICE = 10179
const CURLOPT_SOCKS5_GSSAPI_NEC = 180
const CURLOPT_PROTOCOLS = 181
const CURLOPT_REDIR_PROTOCOLS = 182
const CURLOPT_SSH_KNOWNHOSTS = 10183
const CURLOPT_SSH_KEYFUNCTION = 20184
const CURLOPT_SSH_KEYDATA = 10185
const CURLOPT_MAIL_FROM = 10186
const CURLOPT_MAIL_RCPT = 10187
const CURLOPT_FTP_USE_PRET = 188
const CURLOPT_RTSP_REQUEST = 189
const CURLOPT_RTSP_SESSION_ID = 10190
const CURLOPT_RTSP_STREAM_URI = 10191
const CURLOPT_RTSP_TRANSPORT = 10192
const CURLOPT_RTSP_CLIENT_CSEQ = 193
const CURLOPT_RTSP_SERVER_CSEQ = 194
const CURLOPT_INTERLEAVEDATA = 10195
const CURLOPT_INTERLEAVEFUNCTION = 20196
const CURLOPT_WILDCARDMATCH = 197
const CURLOPT_CHUNK_BGN_FUNCTION = 20198
const CURLOPT_CHUNK_END_FUNCTION = 20199
const CURLOPT_FNMATCH_FUNCTION = 20200
const CURLOPT_CHUNK_DATA = 10201
const CURLOPT_FNMATCH_DATA = 10202
const CURLOPT_RESOLVE = 10203
const CURLOPT_TLSAUTH_USERNAME = 10204
const CURLOPT_TLSAUTH_PASSWORD = 10205
const CURLOPT_TLSAUTH_TYPE = 10206
const CURLOPT_TRANSFER_ENCODING = 207
const CURLOPT_CLOSESOCKETFUNCTION = 20208
const CURLOPT_CLOSESOCKETDATA = 10209
const CURLOPT_GSSAPI_DELEGATION = 210
const CURLOPT_DNS_SERVERS = 10211
const CURLOPT_ACCEPTTIMEOUT_MS = 212
const CURLOPT_TCP_KEEPALIVE = 213
const CURLOPT_TCP_KEEPIDLE = 214
const CURLOPT_TCP_KEEPINTVL = 215
const CURLOPT_SSL_OPTIONS = 216
const CURLOPT_MAIL_AUTH = 10217
const CURLOPT_LASTENTRY = 10218
# end
@ctypedef CURLoption Int32
# enum ANONYMOUS
const CURL_HTTP_VERSION_NONE = 0
const CURL_HTTP_VERSION_1_0 = 1
const CURL_HTTP_VERSION_1_1 = 2
const CURL_HTTP_VERSION_LAST = 3
# end
# enum CURL_NETRC_OPTION
const CURL_NETRC_IGNORED = 0
const CURL_NETRC_OPTIONAL = 1
const CURL_NETRC_REQUIRED = 2
const CURL_NETRC_LAST = 3
# end
# enum CURL_TLSAUTH
const CURL_TLSAUTH_NONE = 0
const CURL_TLSAUTH_SRP = 1
const CURL_TLSAUTH_LAST = 2
# end
# enum curl_TimeCond
const CURL_TIMECOND_NONE = 0
const CURL_TIMECOND_IFMODSINCE = 1
const CURL_TIMECOND_IFUNMODSINCE = 2
const CURL_TIMECOND_LASTMOD = 3
const CURL_TIMECOND_LAST = 4
# end
@ctypedef curl_TimeCond Int32
# enum CURLformoption
const CURLFORM_NOTHING = 0
const CURLFORM_COPYNAME = 1
const CURLFORM_PTRNAME = 2
const CURLFORM_NAMELENGTH = 3
const CURLFORM_COPYCONTENTS = 4
const CURLFORM_PTRCONTENTS = 5
const CURLFORM_CONTENTSLENGTH = 6
const CURLFORM_FILECONTENT = 7
const CURLFORM_ARRAY = 8
const CURLFORM_OBSOLETE = 9
const CURLFORM_FILE = 10
const CURLFORM_BUFFER = 11
const CURLFORM_BUFFERPTR = 12
const CURLFORM_BUFFERLENGTH = 13
const CURLFORM_CONTENTTYPE = 14
const CURLFORM_CONTENTHEADER = 15
const CURLFORM_FILENAME = 16
const CURLFORM_END = 17
const CURLFORM_OBSOLETE2 = 18
const CURLFORM_STREAM = 19
const CURLFORM_LASTENTRY = 20
# end
@ctypedef CURLformoption Int32
type curl_forms
  option::CURLformoption
  value::Ptr{Uint8}
end
# enum CURLFORMcode
const CURL_FORMADD_OK = 0
const CURL_FORMADD_MEMORY = 1
const CURL_FORMADD_OPTION_TWICE = 2
const CURL_FORMADD_NULL = 3
const CURL_FORMADD_UNKNOWN_OPTION = 4
const CURL_FORMADD_INCOMPLETE = 5
const CURL_FORMADD_ILLEGAL_ARRAY = 6
const CURL_FORMADD_DISABLED = 7
const CURL_FORMADD_LAST = 8
# end
@ctypedef CURLFORMcode Int32
@ctypedef curl_formget_callback Ptr{Void}
type curl_slist
  data::Ptr{Uint8}
  next::Ptr{Void}
end
type curl_certinfo
  num_of_certs::Int32
  certinfo::Ptr{Ptr{Void}}
end
# enum CURLINFO
const CURLINFO_NONE = 0
const CURLINFO_EFFECTIVE_URL = 1048577
const CURLINFO_RESPONSE_CODE = 2097154
const CURLINFO_TOTAL_TIME = 3145731
const CURLINFO_NAMELOOKUP_TIME = 3145732
const CURLINFO_CONNECT_TIME = 3145733
const CURLINFO_PRETRANSFER_TIME = 3145734
const CURLINFO_SIZE_UPLOAD = 3145735
const CURLINFO_SIZE_DOWNLOAD = 3145736
const CURLINFO_SPEED_DOWNLOAD = 3145737
const CURLINFO_SPEED_UPLOAD = 3145738
const CURLINFO_HEADER_SIZE = 2097163
const CURLINFO_REQUEST_SIZE = 2097164
const CURLINFO_SSL_VERIFYRESULT = 2097165
const CURLINFO_FILETIME = 2097166
const CURLINFO_CONTENT_LENGTH_DOWNLOAD = 3145743
const CURLINFO_CONTENT_LENGTH_UPLOAD = 3145744
const CURLINFO_STARTTRANSFER_TIME = 3145745
const CURLINFO_CONTENT_TYPE = 1048594
const CURLINFO_REDIRECT_TIME = 3145747
const CURLINFO_REDIRECT_COUNT = 2097172
const CURLINFO_PRIVATE = 1048597
const CURLINFO_HTTP_CONNECTCODE = 2097174
const CURLINFO_HTTPAUTH_AVAIL = 2097175
const CURLINFO_PROXYAUTH_AVAIL = 2097176
const CURLINFO_OS_ERRNO = 2097177
const CURLINFO_NUM_CONNECTS = 2097178
const CURLINFO_SSL_ENGINES = 4194331
const CURLINFO_COOKIELIST = 4194332
const CURLINFO_LASTSOCKET = 2097181
const CURLINFO_FTP_ENTRY_PATH = 1048606
const CURLINFO_REDIRECT_URL = 1048607
const CURLINFO_PRIMARY_IP = 1048608
const CURLINFO_APPCONNECT_TIME = 3145761
const CURLINFO_CERTINFO = 4194338
const CURLINFO_CONDITION_UNMET = 2097187
const CURLINFO_RTSP_SESSION_ID = 1048612
const CURLINFO_RTSP_CLIENT_CSEQ = 2097189
const CURLINFO_RTSP_SERVER_CSEQ = 2097190
const CURLINFO_RTSP_CSEQ_RECV = 2097191
const CURLINFO_PRIMARY_PORT = 2097192
const CURLINFO_LOCAL_IP = 1048617
const CURLINFO_LOCAL_PORT = 2097194
const CURLINFO_LASTONE = 42
# end
@ctypedef CURLINFO Int32
# enum curl_closepolicy
const CURLCLOSEPOLICY_NONE = 0
const CURLCLOSEPOLICY_OLDEST = 1
const CURLCLOSEPOLICY_LEAST_RECENTLY_USED = 2
const CURLCLOSEPOLICY_LEAST_TRAFFIC = 3
const CURLCLOSEPOLICY_SLOWEST = 4
const CURLCLOSEPOLICY_CALLBACK = 5
const CURLCLOSEPOLICY_LAST = 6
# end
@ctypedef curl_closepolicy Int32
# enum curl_lock_data
const CURL_LOCK_DATA_NONE = 0
const CURL_LOCK_DATA_SHARE = 1
const CURL_LOCK_DATA_COOKIE = 2
const CURL_LOCK_DATA_DNS = 3
const CURL_LOCK_DATA_SSL_SESSION = 4
const CURL_LOCK_DATA_CONNECT = 5
const CURL_LOCK_DATA_LAST = 6
# end
@ctypedef curl_lock_data Int32
# enum curl_lock_access
const CURL_LOCK_ACCESS_NONE = 0
const CURL_LOCK_ACCESS_SHARED = 1
const CURL_LOCK_ACCESS_SINGLE = 2
const CURL_LOCK_ACCESS_LAST = 3
# end
@ctypedef curl_lock_access Int32
@ctypedef curl_lock_function Ptr{Void}
@ctypedef curl_unlock_function Ptr{Void}
@ctypedef CURLSH None
# enum CURLSHcode
const CURLSHE_OK = 0
const CURLSHE_BAD_OPTION = 1
const CURLSHE_IN_USE = 2
const CURLSHE_INVALID = 3
const CURLSHE_NOMEM = 4
const CURLSHE_NOT_BUILT_IN = 5
const CURLSHE_LAST = 6
# end
@ctypedef CURLSHcode Int32
# enum CURLSHoption
const CURLSHOPT_NONE = 0
const CURLSHOPT_SHARE = 1
const CURLSHOPT_UNSHARE = 2
const CURLSHOPT_LOCKFUNC = 3
const CURLSHOPT_UNLOCKFUNC = 4
const CURLSHOPT_USERDATA = 5
const CURLSHOPT_LAST = 6
# end
@ctypedef CURLSHoption Int32
# enum CURLversion
const CURLVERSION_FIRST = 0
const CURLVERSION_SECOND = 1
const CURLVERSION_THIRD = 2
const CURLVERSION_FOURTH = 3
const CURLVERSION_LAST = 4
# end
@ctypedef CURLversion Int32
type curl_version_info_data
  age::CURLversion
  version::Ptr{Uint8}
  version_num::Uint32
  host::Ptr{Uint8}
  features::Int32
  ssl_version::Ptr{Uint8}
  ssl_version_num::Int32
  libz_version::Ptr{Uint8}
  protocols::Ptr{Ptr{Uint8}}
  ares::Ptr{Uint8}
  ares_num::Int32
  libidn::Ptr{Uint8}
  iconv_ver_num::Int32
  libssh_version::Ptr{Uint8}
end
@ctypedef CURLM None
# enum CURLMcode
const CURLM_CALL_MULTI_PERFORM = -1
const CURLM_OK = 0
const CURLM_BAD_HANDLE = 1
const CURLM_BAD_EASY_HANDLE = 2
const CURLM_OUT_OF_MEMORY = 3
const CURLM_INTERNAL_ERROR = 4
const CURLM_BAD_SOCKET = 5
const CURLM_UNKNOWN_OPTION = 6
const CURLM_LAST = 7
# end
@ctypedef CURLMcode Int32
# enum CURLMSG
const CURLMSG_NONE = 0
const CURLMSG_DONE = 1
const CURLMSG_LAST = 2
# end
@ctypedef CURLMSG Int32
type CURLMsg
  msg::CURLMSG
  easy_handle::Ptr{CURL}
  data::Void
end
@ctypedef curl_socket_callback Ptr{Void}
@ctypedef curl_multi_timer_callback Ptr{Void}
# enum CURLMoption
const CURLMOPT_SOCKETFUNCTION = 20001
const CURLMOPT_SOCKETDATA = 10002
const CURLMOPT_PIPELINING = 3
const CURLMOPT_TIMERFUNCTION = 20004
const CURLMOPT_TIMERDATA = 10005
const CURLMOPT_MAXCONNECTS = 6
const CURLMOPT_LASTENTRY = 7
# end
@ctypedef CURLMoption Int32


text:
"Omega Expression"
abstract type OExpr end

struct LambdaExpr <: OExpr
  var::Symbol
  body::OExpr
end

"Function Application Expr"
struct AppExpr <: OExpr
  f::Expr
  args::Vector{Expr}
end

struct LetExpr <: OExpr
  var::Symbol
  val::OExpr
  body::OExpr
end

ext: jl


text:
module GridArrays

using DomainSets, StaticArrays, RecipesBase, Test, FastGaussQuadrature,
        GaussQuadrature, FillArrays

using CompositeTypes, CompositeTypes.Display

using DomainSets: endpoints

@deprecate AbstractSubGrid SubGrid
@deprecate MaskedSubGrid MaskedGrid
@deprecate mapping forward_map
@deprecate mapped_grid map_grid

## List of imports

import Base: size, length, @propagate_inbounds, step, ndims,
    checkbounds, IndexStyle, ==, ≈, getindex, eachindex,
    convert, in, ^, string, axes,
    show

import CompositeTypes:
        iscomposite, component, components, ncomponents

import DomainSets:
        factors, factor, nfactors,
        ×, cross,
        dimension, prectype, numtype,
        minimum, maximum,
        boundary,
        mapto_canonical, mapfrom_canonical,
        forward_map, inverse_map

## List of exports

# from generic/grid.jl
export AbstractGrid, AbstractGrid1d,
        numtype, prectype,
        dimension,
        resize,
        covering,
        isperiodic,
        boundary,
        canonicalgrid

# from generic/product.jl
export ProductGrid, productgrid,
        VcatGrid, FlatProductGrid

# from subgrid/subgrid.jl
export subindices, supergrid, issubindex, similar_subgrid

# from generic/mapped.jl
export MappedGrid, map_grid, apply_map

# from subgrid/subgrid.jl
export subgrid, mask, ProductSubGrid

# from domains/interval.jl
export AbstractIntervalGrid,
        AbstractEquispacedGrid,
        EquispacedGrid, PeriodicEquispacedGrid, MidpointEquispacedGrid,
        UnitEquispacedGrid, UnitPeriodicEquispacedGrid, UnitMidpointEquispacedGrid,
        FourierGrid,
        step

# from domains/scattered.jl
export ScatteredGrid, randomgrid

# from domains/discretize.jl
export discretize_togrid

# from applications/gauss.jl
export ChebyshevTNodes, ChebyshevUNodes, ChebyshevNodes,
        ChebyshevGrid, ChebyshevPoints, ChebyshevExtremae,
        LaguerreNodes, HermiteNodes, LegendreNodes, JacobiNodes

include("util/PeriodicCartesianIndices.jl")
using ..PeriodicIndexing

include("util/common.jl")

include("generic/grid.jl")
include("generic/composite.jl")
include("generic/product.jl")
include("generic/mapped.jl")

include("domains/interval.jl")
include("domains/scattered.jl")
include("domains/discretize.jl")

include("subgrid/subgrid.jl")

include("applications/gauss.jl")

include("util/recipes.jl")

# We define some testing routines inside the package, so that they
# can also be used in the tests of other packages that extend grids
include("test/Test.jl")

export coverdomain
coverdomain(g) = covering(g)

end # module


text:
using Distributed, JSON

@everywhere begin
    using Pkg
    Pkg.activate(".")
end

@everywhere using CellularInfo

@everywhere function transfer(traj, src, dst, k)
    #  back = filter(f -> f != src && f != dst, 1:size(traj,2))
    #  cond = traj[:,back]
    #  transferentropy(traj[:,src], traj[:,dst], cond', k, 1000)
    transferentropy(traj[:,src], traj[:,dst], k, 1000)
end

function main()
    meante = zeros(Float64, 256)
    for rule in 0:255
        @time begin
            @info "Building graph for $(rule)..."
            k = 5
            width = 101
            duration = 200 + (width ÷ 2)

            a = zeros(Int, width)
            b = ones(Int, width)
            a[width ÷ 2 + 1] = 1
            b[width ÷ 2 + 1] = 0

            traja = simulate(rule, a, duration)[width ÷ 2 + 1:end,:]
            trajb = simulate(rule, b, duration)[width ÷ 2 + 1:end,:]

            N = 2*width*(width - 1)
            μ = 0.0
            for src in 1:width
                μ += @distributed (+) for dst in 1:width
                    if src == dst
                        0.0
                    else
                        tea = transfer(traja, src, dst, k)
                        teb = transfer(trajb, src, dst, k)

                        te = (tea.p < 0.05) ? tea.value : 0.0
                        te += (teb.p < 0.05) ? teb.value : 0.0
                    end
                end
            end
            meante[rule + 1] = μ / N
        end
    end

    open("onepoint.json", "w") do handle
        JSON.print(handle, meante)
    end
end

main()

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text: include("rescorrs.jl")


text:
# ---
# title: 198. House Robber
# id: problem198
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/house-robber/description/>
# hidden: true
# ---
#
# You are a professional robber planning to rob houses along a street. Each
# house has a certain amount of money stashed, the only constraint stopping you
# from robbing each of them is that adjacent houses have security system
# connected and **it will automatically contact the police if two adjacent
# houses were broken into on the same night**.
#
# Given a list of non-negative integers representing the amount of money of each
# house, determine the maximum amount of money you can rob tonight **without
# alerting the police**.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,1]
#     Output: 4
#     Explanation: Rob house 1 (money = 1) and then rob house 3 (money = 3).
#                  Total amount you can rob = 1 + 3 = 4.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,7,9,3,1]
#     Output: 12
#     Explanation: Rob house 1 (money = 2), rob house 3 (money = 9) and rob house 5 (money = 1).
#                  Total amount you can rob = 2 + 9 + 1 = 12.
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 100`
#   * `0 <= nums[i] <= 400`
#
#
## @lc code=start
using LeetCode

function rob_198(nums::Vector{Int})::Int
    n = length(nums)
    n == 1 && return nums[1]
    first, second = nums[1], max(nums[1], nums[2])
    for num in @view nums[3:end]
        first, second = second, max(second, first + num)
    end
    return second
end

## @lc code=end


text:
using Test
using ComparativeAutograder

@testset "Test runtestcase on built-in functions" begin
    @testset "testcase for println" begin
        case = TestCase("foo")
        result = runtestcase(println, case)
        @test result.result === nothing
        @test result.exception === nothing
        @test result.stdout == "foo\n"
        @test result.stderr == ""
    end

    @testset "Test runtestcase on built-in +" begin
        case = TestCase(1, 2)
        result = runtestcase(+, case)
        @test result.result === 3
        @test result.exception === nothing
        @test result.stdout == ""
        @test result.stderr == ""

        case = TestCase(1, "foo")
        result = runtestcase(+, case)
        @test result.result === nothing
        @test occursin("MethodError", result.exception)
        @test occursin("no method matching", result.backtrace)
        @test result.stdout == ""
        @test result.stderr == ""
    end
end

@testset "Test runtestcase on custom functions" begin
    @testset "testcase on custom add" begin
        add = (x, y) -> x + y
        case = TestCase(3, 4)
        result = runtestcase(add, case)
        @test result.result == 7
        @test result.exception === nothing
    end

    @testset "Test runtestcase on error throwing function" begin
        func = () -> error("I'm a error!")
        case = TestCase()
        result = runtestcase(func, case)
        @test result.result === nothing
        @test occursin("ErrorException", result.exception)
    end
end


text:
struct EOF <: Exception end

const _one = one(UInt32)

"""
    Simple datastruct that uses a dictionary as tape with a potentially multiple cell pointers.
    The time travel aspect is implemented by storing the values & pointers of each change in the tape
    in a stack-like vector and computing the prior state when needed.
    Each Timeline also keeps track of its instruction pointer
"""
mutable struct Timeline
    tape::Dict{Int, UInt8}
    memory_pointers::Vector{UInt32}
    instruction_pointer::Int
    history::Vector{Dict{Int, UInt8}}
    Timeline(tape::Dict{Int, UInt8} = Dict(1 => 0x0), mems = UInt32[_one], instr = 1, history = []) = new(tape, mems, instr, history)
end

function error(msg::AbstractString)
    printstyled("error: ", msg, '\n', color=:red)
    exit(-1)
end

"""
    Throws an error if the input is invalid. Otherwise it returns a char array of the input
"""
function assert()
    length(ARGS) < 1 && error("missing source file.\nusage: julia interpreter.jl src.5dbfwmtt")
    !isfile(ARGS[1]) && error("given source file is not a valid file")
    if !endswith(ARGS[1], ".5dbfwmtt") && !endswith(ARGS[1], ".bf")
        printstyled("warning: source file is not a .5dbfwmtt file", color=:red)
    end
    open(ARGS[1], "r") do f
        res = Char[]
        while !eof(f)
            c = read(f, Char)::Char
            c ∈ "><+-.,[]()~@^v" && push!(res, c)
        end
        res
    end  # O( l )
end

"""
    Create a minimal snapshot of all the cells currently pointed to
"""
snapshot(tl::Timeline) = push!(tl.history, Dict([cell => tl.tape[cell] for cell in tl.memory_pointers]))

"""
    Move all memory pointers in this timeline 1 cell to the right.
"""
function move(tl::Timeline, ::Int, ::Vector{Timeline})
    tl.memory_pointers .+= _one
    # initializing cells to zero (if they have no value)
    get!.(Ref(tl.tape), tl.memory_pointers, 0x0)
end

"""
    Move all memory pointers in this timeline 1 cell to the left.
    (note that bf interpreters usually start pointing at the leftmost cell)
"""
function back(tl::Timeline, ::Int, ::Vector{Timeline})
    any([m === _one for m in tl.memory_pointers]) && error("out of bounds")
    tl.memory_pointers .-= _one
    get!.(Ref(tl.tape), tl.memory_pointers, 0x0)
end

"""
    Increment all cells pointed to in this timeline.
"""
function incr(tl::Timeline, ::Int, ::Vector{Timeline})
    snapshot(tl)
    (cell -> tl.tape[cell] += 0x1).(tl.memory_pointers)
end

"""
    Decrement all cells pointed to in this timeline.
"""
function decr(tl::Timeline, ::Int, ::Vector{Timeline})
    snapshot(tl)
    (cell -> tl.tape[cell] -= 0x1).(tl.memory_pointers)
end

"""
    Output a character for all cells pointed to in this timeline.
    (using ascii encoding)
"""
function putc(tl::Timeline, ::Int, ::Vector{Timeline})
    print.(Char.(getindex.(Ref(tl.tape), tl.memory_pointers)))
end

"""
    Input a character and store it in all cells pointed to in this timeline.
"""
function getc(tl::Timeline, ::Int, ::Vector{Timeline})
    snapshot(tl)
    setindex!.(Ref(tl.tape), read(stdin, UInt8), tl.memory_pointers)
end

"""
    Move this instruction pointer to the matching ] if all cells pointed to in this timeline are 0.
"""
function loop(tl::Timeline, val::UInt32)
    if all([tl.tape[cell] == 0 for cell in tl.memory_pointers])
        tl.instruction_pointer = val # massive speed gain right here
    end
end

"""
    Move this instruction pointer back to the matching [ if any cells pointed to in this timeline are nonzero.
"""
function goto(tl::Timeline, val::UInt32)
    if any([tl.tape[cell] ≠ 0 for cell ∈ tl.memory_pointers])
        tl.instruction_pointer = val
    end
end

"""
    Rewind the current tape back in time by 1 step.
"""
function time(tl::Timeline, ::Int, ::Vector{Timeline})
    length(tl.history) > 0 && merge!(tl.tape, pop!(tl.history))
end

"""
    Spawn a parallel timeline below the current timeline, with a copy of the tape and all pointers in it.
    This instruction pointer jumps to the matching ).
    Spawn a new instruction pointer within the newly spawned timeline,
    beginning execution immediately after this instruction.
"""
function copy(tl::Timeline, tl_pointer::Int, timelines::Vector{Timeline}, val::UInt32)
    ntl = deepcopy(tl)
    ntl.instruction_pointer += 1
    tl.instruction_pointer = val # jump to the matching )
    insert!(timelines, tl_pointer + 1, ntl)
end

"""
    If this is executed outside of the main timeline,
    kill this timeline and all the memory/instruction pointers currently in it.
    Otherwise, do nothing.
"""
function kill(::Timeline, tl_pointer::Int, timelines::Vector{Timeline})
    if tl_pointer ≠ 1
        popat!(timelines, tl_pointer)
    end
end

"""
    Move all memory pointers in this timeline to the same location in the next ("lower") parallel universe.
"""
function next(tl::Timeline, tl_pointer::Int, timelines::Vector{Timeline})
    if tl_pointer ≠ length(timelines)
        append!(timelines[tl_pointer + 1].memory_pointers, tl.memory_pointers)
        get!.(Ref(timelines[tl_pointer + 1].tape), tl.memory_pointers, 0x0)
    end
    tl.memory_pointers = UInt32[]
end

"""
    Move all memory pointers in this timeline to the same location in the previous ("higher") parallel universe.
"""
function prev(tl::Timeline, tl_pointer::Int, timelines::Vector{Timeline})
    if tl_pointer ≠ 1
        append!(timelines[tl_pointer - 1].memory_pointers, tl.memory_pointers)
        # defaulting new cells to 0
        get!.(Ref(timelines[tl_pointer - 1].tape), tl.memory_pointers, 0x0)
    end
    tl.memory_pointers = UInt32[]
end

"""
    If the lower timeline has a nonzero amount of memory pointers,
    freeze this program pointer for the duration of this step.
    Otherwise, continue forward.
"""
function wait(tl::Timeline, tl_pointer::Int, timelines::Vector{Timeline})
    # has a lower tl & that tl has mem pointers
    if length(timelines) > tl_pointer && length(timelines[tl_pointer + 1].memory_pointers) > 0
        tl.instruction_pointer -= 1
    end
end

"""
    Ends this program if the current timeline is
    the main timeline. Otherwise it terminates the timeline
"""
done(::Timeline, tl_pointer::Int, timelines::Vector{Timeline}) = tl_pointer == 1 ? exit(0) : popat!(timelines, tl_pointer)

"""
    Converts all chars into their respective token.
    This is kinda useless for a language that consists of only 1-char instructions
    But I have it here for proper comments between # and eol's
"""
function lexer(src::Vector{Char})
    # vanilla bf
    d = Dict('>'=>move,'<'=>back,'+'=>incr,'-'=>decr,'.'=>putc,','=>getc,'['=>loop,']'=>goto,
    # 5D + time travel extension
    '~'=>time,'('=>copy,')'=>kill,'v'=>next,'^'=>prev,'@'=>wait)
    res = Function[]
    k = _one
    loops = UInt32[]
    while k ≤ length(src)
        p = src[k]
        if p ∈ "><+-.,~v^@)"
            push!(res, d[p])
        elseif p === '['
            n = 1
            f = _one
            while n > 0
                f + k > length(src) && error("expected closing bracket")
                n += (src[f + k] === '[') - (src[f + k] === ']')
                f += _one
            end
            push!(loops, k)
            v = f + k - _one
            push!(res, (tl::Timeline, ::Int, ::Vector{Timeline}) -> loop(tl, v))
        elseif p === ']'
            isempty(loops) && error("unexpected closing bracket")
            v = pop!(loops)
            push!(res, (tl::Timeline, ::Int, ::Vector{Timeline}) -> goto(tl, v))
        elseif p === '('
            n = 1
            f = _one
            while n > 0
                f + k > length(src) && error("expected closing parenthesis")
                n += (src[f + k] === '(') - (src[f + k] === ')')
                f += _one
            end
            v = f + k - _one
            push!(res, (tl::Timeline, tl_pointer::Int, timelines::Vector{Timeline}) -> copy(tl, tl_pointer, timelines, v))
        end
        k += _one
    end
    return push!(res, done)
end

function interpret(src::Vector{Function}) # worst: O( i d m ) average: O( i d ) best: O( i )
    timelines = Timeline[Timeline()]
    # O( i ) ; i = total instructions
    while true # program ends when the main timeline calls done()
        tl_pointer = 1
        # O( d ) ; d = depths of timelines
        while tl_pointer ≤ length(timelines)
            tl = timelines[tl_pointer]
            # all of these functions have either constant time O( 1 )
            # or vectorized & efficient O( m ) ; m = amount of memory pointers
            src[tl.instruction_pointer](tl, tl_pointer, timelines)
            tl_pointer += 1
            tl.instruction_pointer += 1
        end
    end
end

(interpret ∘ lexer ∘ assert)()


text:
module Metrics

"""
    creategrid!(x1, elemtocoord, r)

Create a 1-D grid using `elemtocoord` (see [`brickmesh`](@ref)) using the 1-D
`(-1, 1)` reference coordinates `r`. The element grids are filled using linear
interpolation of the element coordinates.

If `Nq = length(r)` and `nelem = size(elemtocoord, 3)` then the preallocated
array `x1` should be `Nq * nelem == length(x1)`.
"""
function creategrid!(x1, e2c, r)
    (d, nvert, nelem) = size(e2c)
    @assert d == 1
    Nq = length(r)
    x1 = reshape(x1, (Nq, nelem))

    # linear blend
    @inbounds for e in 1:nelem
        for i in 1:Nq
            x1[i, e] =
                ((1 - r[i]) * e2c[1, 1, e] + (1 + r[i]) * e2c[1, 2, e]) / 2
        end
    end
    nothing
end

"""
    creategrid!(x1, x2, elemtocoord, r)

Create a 2-D tensor product grid using `elemtocoord` (see [`brickmesh`](@ref))
using the 1-D `(-1, 1)` reference coordinates `r`. The element grids are filled
using bilinear interpolation of the element coordinates.

If `Nq = length(r)` and `nelem = size(elemtocoord, 3)` then the preallocated
arrays `x1` and `x2` should be `Nq^2 * nelem == size(x1) == size(x2)`.
"""
function creategrid!(x1, x2, e2c, r)
    (d, nvert, nelem) = size(e2c)
    @assert d == 2
    Nq = length(r)
    x1 = reshape(x1, (Nq, Nq, nelem))
    x2 = reshape(x2, (Nq, Nq, nelem))

    # bilinear blend of corners
    @inbounds for (f, n) in zip((x1, x2), 1:d)
        for e in 1:nelem
            for j in 1:Nq
                for i in 1:Nq
                    f[i, j, e] =
                        (
                            (1 - r[i]) * (1 - r[j]) * e2c[n, 1, e] +
                            (1 + r[i]) * (1 - r[j]) * e2c[n, 2, e] +
                            (1 - r[i]) * (1 + r[j]) * e2c[n, 3, e] +
                            (1 + r[i]) * (1 + r[j]) * e2c[n, 4, e]
                        ) / 4
                end
            end
        end
    end
    nothing
end

"""
    creategrid!(x1, x2, x3, elemtocoord, r)

Create a 3-D tensor product grid using `elemtocoord` (see [`brickmesh`](@ref))
using the 1-D `(-1, 1)` reference coordinates `r`. The element grids are filled
using trilinear interpolation of the element coordinates.

If `Nq = length(r)` and `nelem = size(elemtocoord, 3)` then the preallocated
arrays `x1`, `x2`, and `x3` should be `Nq^3 * nelem == size(x1) == size(x2) ==
size(x3)`.
"""
function creategrid!(x1, x2, x3, e2c, r)
    (d, nvert, nelem) = size(e2c)
    @assert d == 3
    # TODO: Add asserts?
    Nq = length(r)
    x1 = reshape(x1, (Nq, Nq, Nq, nelem))
    x2 = reshape(x2, (Nq, Nq, Nq, nelem))
    x3 = reshape(x3, (Nq, Nq, Nq, nelem))

    # trilinear blend of corners
    @inbounds for (f, n) in zip((x1, x2, x3), 1:d)
        for e in 1:nelem
            for k in 1:Nq
                for j in 1:Nq
                    for i in 1:Nq
                        f[i, j, k, e] =
                            (
                                (1 - r[i]) *
                                (1 - r[j]) *
                                (1 - r[k]) *
                                e2c[n, 1, e] +
                                (1 + r[i]) *
                                (1 - r[j]) *
                                (1 - r[k]) *
                                e2c[n, 2, e] +
                                (1 - r[i]) *
                                (1 + r[j]) *
                                (1 - r[k]) *
                                e2c[n, 3, e] +
                                (1 + r[i]) *
                                (1 + r[j]) *
                                (1 - r[k]) *
                                e2c[n, 4, e] +
                                (1 - r[i]) *
                                (1 - r[j]) *
                                (1 + r[k]) *
                                e2c[n, 5, e] +
                                (1 + r[i]) *
                                (1 - r[j]) *
                                (1 + r[k]) *
                                e2c[n, 6, e] +
                                (1 - r[i]) *
                                (1 + r[j]) *
                                (1 + r[k]) *
                                e2c[n, 7, e] +
                                (1 + r[i]) *
                                (1 + r[j]) *
                                (1 + r[k]) *
                                e2c[n, 8, e]
                            ) / 8
                    end
                end
            end
        end
    end
    nothing
end

"""
    computemetric!(x1, J, ξ1x1, sJ, n1, D)

Compute the 1-D metric terms from the element grid arrays `x1`. All the arrays
are preallocated by the user and the (square) derivative matrix `D` should be
consistent with the reference grid `r` used in [`creategrid!`](@ref).

If `Nq = size(D, 1)` and `nelem = div(length(x1), Nq)` then the volume arrays
`x1`, `J`, and `ξ1x1` should all have length `Nq * nelem`.  Similarly, the face
arrays `sJ` and `n1` should be of length `nface * nelem` with `nface = 2`.
"""
function computemetric!(x1, J, ξ1x1, sJ, n1, D)
    Nq = size(D, 1)
    nelem = div(length(J), Nq)
    x1 = reshape(x1, (Nq, nelem))
    J = reshape(J, (Nq, nelem))
    ξ1x1 = reshape(ξ1x1, (Nq, nelem))
    nface = 2
    n1 = reshape(n1, (1, nface, nelem))
    sJ = reshape(sJ, (1, nface, nelem))

    d = 1

    @inbounds for e in 1:nelem
        J[:, e] = D * x1[:, e]
    end
    ξ1x1 .= 1 ./ J

    n1[1, 1, :] .= -sign.(J[1, :])
    n1[1, 2, :] .= sign.(J[Nq, :])
    sJ .= 1
    nothing
end

"""
    computemetric!(x1, x2, J, ξ1x1, ξ2x1, ξ1x2, ξ2x2, sJ, n1, n2, D)

Compute the 2-D metric terms from the element grid arrays `x1` and `x2`. All the
arrays are preallocated by the user and the (square) derivative matrix `D`
should be consistent with the reference grid `r` used in [`creategrid!`](@ref).

If `Nq = size(D, 1)` and `nelem = div(length(x1), Nq^2)` then the volume arrays
`x1`, `x2`, `J`, `ξ1x1`, `ξ2x1`, `ξ1x2`, and `ξ2x2` should all be of size `(Nq, Nq,
nelem)`.  Similarly, the face arrays `sJ`, `n1`, and `n2` should be of size
`(Nq, nface, nelem)` with `nface = 4`.
"""
function computemetric!(x1, x2, J, ξ1x1, ξ2x1, ξ1x2, ξ2x2, sJ, n1, n2, D)
    T = eltype(x1)
    Nq = size(D, 1)
    nelem = div(length(J), Nq^2)
    d = 2
    x1 = reshape(x1, (Nq, Nq, nelem))
    x2 = reshape(x2, (Nq, Nq, nelem))
    J = reshape(J, (Nq, Nq, nelem))
    ξ1x1 = reshape(ξ1x1, (Nq, Nq, nelem))
    ξ2x1 = reshape(ξ2x1, (Nq, Nq, nelem))
    ξ1x2 = reshape(ξ1x2, (Nq, Nq, nelem))
    ξ2x2 = reshape(ξ2x2, (Nq, Nq, nelem))
    nface = 4
    n1 = reshape(n1, (Nq, nface, nelem))
    n2 = reshape(n2, (Nq, nface, nelem))
    sJ = reshape(sJ, (Nq, nface, nelem))

    @inbounds for e in 1:nelem
        for j in 1:Nq, i in 1:Nq
            x1ξ1 = x1ξ2 = zero(T)
            x2ξ1 = x2ξ2 = zero(T)
            for n in 1:Nq
                x1ξ1 += D[i, n] * x1[n, j, e]
                x1ξ2 += D[j, n] * x1[i, n, e]
                x2ξ1 += D[i, n] * x2[n, j, e]
                x2ξ2 += D[j, n] * x2[i, n, e]
            end
            J[i, j, e] = x1ξ1 * x2ξ2 - x2ξ1 * x1ξ2
            ξ1x1[i, j, e] = x2ξ2 / J[i, j, e]
            ξ2x1[i, j, e] = -x2ξ1 / J[i, j, e]
            ξ1x2[i, j, e] = -x1ξ2 / J[i, j, e]
            ξ2x2[i, j, e] = x1ξ1 / J[i, j, e]
        end

        for i in 1:Nq
            n1[i, 1, e] = -J[1, i, e] * ξ1x1[1, i, e]
            n2[i, 1, e] = -J[1, i, e] * ξ1x2[1, i, e]
            n1[i, 2, e] = J[Nq, i, e] * ξ1x1[Nq, i, e]
            n2[i, 2, e] = J[Nq, i, e] * ξ1x2[Nq, i, e]
            n1[i, 3, e] = -J[i, 1, e] * ξ2x1[i, 1, e]
            n2[i, 3, e] = -J[i, 1, e] * ξ2x2[i, 1, e]
            n1[i, 4, e] = J[i, Nq, e] * ξ2x1[i, Nq, e]
            n2[i, 4, e] = J[i, Nq, e] * ξ2x2[i, Nq, e]

            for n in 1:4
                sJ[i, n, e] = hypot(n1[i, n, e], n2[i, n, e])
                n1[i, n, e] /= sJ[i, n, e]
                n2[i, n, e] /= sJ[i, n, e]
            end
        end
    end

    nothing
end

"""
    computemetric!(x1, x2, x3, J, ξ1x1, ξ2x1, ξ3x1, ξ1x2, ξ2x2, ξ3x2, ξ1x3,
                   ξ2x3, ξ3x3, sJ, n1, n2, n3, D)

Compute the 3-D metric terms from the element grid arrays `x1`, `x2`, and `x3`.
All the arrays are preallocated by the user and the (square) derivative matrix
`D` should be consistent with the reference grid `r` used in
[`creategrid!`](@ref).

If `Nq = size(D, 1)` and `nelem = div(length(x1), Nq^3)` then the volume arrays
`x1`, `x2`, `x3`, `J`, `ξ1x1`, `ξ2x1`, `ξ3x1`, `ξ1x2`, `ξ2x2`, `ξ3x2`, `ξ1x3`,
`ξ2x3`, and `ξ3x3` should all be of length `Nq^3 * nelem`.  Similarly, the face
arrays `sJ`, `n1`, `n2`, and `n3` should be of size `Nq^2 * nface * nelem` with
`nface = 6`.

The curl invariant formulation of Kopriva (2006), equation 37, is used.

Reference:
  Kopriva, David A. "Metric identities and the discontinuous spectral element
  method on curvilinear meshes." Journal of Scientific Computing 26.3 (2006):
  301-327. <https://doi.org/10.1007/s10915-005-9070-8>
"""
function computemetric!(
    x1,
    x2,
    x3,
    J,
    ξ1x1,
    ξ2x1,
    ξ3x1,
    ξ1x2,
    ξ2x2,
    ξ3x2,
    ξ1x3,
    ξ2x3,
    ξ3x3,
    sJ,
    n1,
    n2,
    n3,
    D,
)
    T = eltype(x1)

    Nq = size(D, 1)
    nelem = div(length(J), Nq^3)

    x1 = reshape(x1, (Nq, Nq, Nq, nelem))
    x2 = reshape(x2, (Nq, Nq, Nq, nelem))
    x3 = reshape(x3, (Nq, Nq, Nq, nelem))
    J = reshape(J, (Nq, Nq, Nq, nelem))
    ξ1x1 = reshape(ξ1x1, (Nq, Nq, Nq, nelem))
    ξ2x1 = reshape(ξ2x1, (Nq, Nq, Nq, nelem))
    ξ3x1 = reshape(ξ3x1, (Nq, Nq, Nq, nelem))
    ξ1x2 = reshape(ξ1x2, (Nq, Nq, Nq, nelem))
    ξ2x2 = reshape(ξ2x2, (Nq, Nq, Nq, nelem))
    ξ3x2 = reshape(ξ3x2, (Nq, Nq, Nq, nelem))
    ξ1x3 = reshape(ξ1x3, (Nq, Nq, Nq, nelem))
    ξ2x3 = reshape(ξ2x3, (Nq, Nq, Nq, nelem))
    ξ3x3 = reshape(ξ3x3, (Nq, Nq, Nq, nelem))

    nface = 6
    #= This code is broken when views are used
    n1 = reshape(n1, Nq, Nq, nface, nelem)
    n2 = reshape(n2, Nq, Nq, nface, nelem)
    n3 = reshape(n3, Nq, Nq, nface, nelem)
    sJ = reshape(sJ, Nq, Nq, nface, nelem)
    =#

    JI2 = similar(J, (Nq, Nq, Nq))
    (yzr, yzs, yzt) = (similar(JI2), similar(JI2), similar(JI2))
    (zxr, zxs, zxt) = (similar(JI2), similar(JI2), similar(JI2))
    (xyr, xys, xyt) = (similar(JI2), similar(JI2), similar(JI2))

    ξ1x1 .= zero(T)
    ξ2x1 .= zero(T)
    ξ3x1 .= zero(T)
    ξ1x2 .= zero(T)
    ξ2x2 .= zero(T)
    ξ3x2 .= zero(T)
    ξ1x3 .= zero(T)
    ξ2x3 .= zero(T)
    ξ3x3 .= zero(T)

    @inbounds for e in 1:nelem
        for k in 1:Nq, j in 1:Nq, i in 1:Nq
            x1ξ1 = x1ξ2 = x1ξ3 = zero(T)
            x2ξ1 = x2ξ2 = x2ξ3 = zero(T)
            x3ξ1 = x3ξ2 = x3ξ3 = zero(T)
            for n in 1:Nq
                x1ξ1 += D[i, n] * x1[n, j, k, e]
                x1ξ2 += D[j, n] * x1[i, n, k, e]
                x1ξ3 += D[k, n] * x1[i, j, n, e]
                x2ξ1 += D[i, n] * x2[n, j, k, e]
                x2ξ2 += D[j, n] * x2[i, n, k, e]
                x2ξ3 += D[k, n] * x2[i, j, n, e]
                x3ξ1 += D[i, n] * x3[n, j, k, e]
                x3ξ2 += D[j, n] * x3[i, n, k, e]
                x3ξ3 += D[k, n] * x3[i, j, n, e]
            end
            J[i, j, k, e] = (
                x1ξ1 * (x2ξ2 * x3ξ3 - x3ξ2 * x2ξ3) +
                x2ξ1 * (x3ξ2 * x1ξ3 - x1ξ2 * x3ξ3) +
                x3ξ1 * (x1ξ2 * x2ξ3 - x2ξ2 * x1ξ3)
            )

            JI2[i, j, k] = 1 / (2 * J[i, j, k, e])

            yzr[i, j, k] = x2[i, j, k, e] * x3ξ1 - x3[i, j, k, e] * x2ξ1
            yzs[i, j, k] = x2[i, j, k, e] * x3ξ2 - x3[i, j, k, e] * x2ξ2
            yzt[i, j, k] = x2[i, j, k, e] * x3ξ3 - x3[i, j, k, e] * x2ξ3
            zxr[i, j, k] = x3[i, j, k, e] * x1ξ1 - x1[i, j, k, e] * x3ξ1
            zxs[i, j, k] = x3[i, j, k, e] * x1ξ2 - x1[i, j, k, e] * x3ξ2
            zxt[i, j, k] = x3[i, j, k, e] * x1ξ3 - x1[i, j, k, e] * x3ξ3
            xyr[i, j, k] = x1[i, j, k, e] * x2ξ1 - x2[i, j, k, e] * x1ξ1
            xys[i, j, k] = x1[i, j, k, e] * x2ξ2 - x2[i, j, k, e] * x1ξ2
            xyt[i, j, k] = x1[i, j, k, e] * x2ξ3 - x2[i, j, k, e] * x1ξ3
        end

        for k in 1:Nq, j in 1:Nq, i in 1:Nq
            for n in 1:Nq
                ξ1x1[i, j, k, e] += D[j, n] * yzt[i, n, k]
                ξ1x1[i, j, k, e] -= D[k, n] * yzs[i, j, n]
                ξ2x1[i, j, k, e] += D[k, n] * yzr[i, j, n]
                ξ2x1[i, j, k, e] -= D[i, n] * yzt[n, j, k]
                ξ3x1[i, j, k, e] += D[i, n] * yzs[n, j, k]
                ξ3x1[i, j, k, e] -= D[j, n] * yzr[i, n, k]
                ξ1x2[i, j, k, e] += D[j, n] * zxt[i, n, k]
                ξ1x2[i, j, k, e] -= D[k, n] * zxs[i, j, n]
                ξ2x2[i, j, k, e] += D[k, n] * zxr[i, j, n]
                ξ2x2[i, j, k, e] -= D[i, n] * zxt[n, j, k]
                ξ3x2[i, j, k, e] += D[i, n] * zxs[n, j, k]
                ξ3x2[i, j, k, e] -= D[j, n] * zxr[i, n, k]
                ξ1x3[i, j, k, e] += D[j, n] * xyt[i, n, k]
                ξ1x3[i, j, k, e] -= D[k, n] * xys[i, j, n]
                ξ2x3[i, j, k, e] += D[k, n] * xyr[i, j, n]
                ξ2x3[i, j, k, e] -= D[i, n] * xyt[n, j, k]
                ξ3x3[i, j, k, e] += D[i, n] * xys[n, j, k]
                ξ3x3[i, j, k, e] -= D[j, n] * xyr[i, n, k]
            end
            ξ1x1[i, j, k, e] *= JI2[i, j, k]
            ξ2x1[i, j, k, e] *= JI2[i, j, k]
            ξ3x1[i, j, k, e] *= JI2[i, j, k]
            ξ1x2[i, j, k, e] *= JI2[i, j, k]
            ξ2x2[i, j, k, e] *= JI2[i, j, k]
            ξ3x2[i, j, k, e] *= JI2[i, j, k]
            ξ1x3[i, j, k, e] *= JI2[i, j, k]
            ξ2x3[i, j, k, e] *= JI2[i, j, k]
            ξ3x3[i, j, k, e] *= JI2[i, j, k]
        end

        for j in 1:Nq, i in 1:Nq
            #= This code is broken when views are used
            n1[i, j, 1, e] = -J[ 1, i, j, e] * ξ1x1[ 1, i, j, e]
            n1[i, j, 2, e] =  J[Nq, i, j, e] * ξ1x1[Nq, i, j, e]
            n1[i, j, 3, e] = -J[ i, 1, j, e] * ξ2x1[ i, 1, j, e]
            n1[i, j, 4, e] =  J[ i,Nq, j, e] * ξ2x1[ i,Nq, j, e]
            n1[i, j, 5, e] = -J[ i, j, 1, e] * ξ3x1[ i, j, 1, e]
            n1[i, j, 6, e] =  J[ i, j,Nq, e] * ξ3x1[ i, j,Nq, e]
            n2[i, j, 1, e] = -J[ 1, i, j, e] * ξ1x2[ 1, i, j, e]
            n2[i, j, 2, e] =  J[Nq, i, j, e] * ξ1x2[Nq, i, j, e]
            n2[i, j, 3, e] = -J[ i, 1, j, e] * ξ2x2[ i, 1, j, e]
            n2[i, j, 4, e] =  J[ i,Nq, j, e] * ξ2x2[ i,Nq, j, e]
            n2[i, j, 5, e] = -J[ i, j, 1, e] * ξ3x2[ i, j, 1, e]
            n2[i, j, 6, e] =  J[ i, j,Nq, e] * ξ3x2[ i, j,Nq, e]
            n3[i, j, 1, e] = -J[ 1, i, j, e] * ξ1x3[ 1, i, j, e]
            n3[i, j, 2, e] =  J[Nq, i, j, e] * ξ1x3[Nq, i, j, e]
            n3[i, j, 3, e] = -J[ i, 1, j, e] * ξ2x3[ i, 1, j, e]
            n3[i, j, 4, e] =  J[ i,Nq, j, e] * ξ2x3[ i,Nq, j, e]
            n3[i, j, 5, e] = -J[ i, j, 1, e] * ξ3x3[ i, j, 1, e]
            n3[i, j, 6, e] =  J[ i, j,Nq, e] * ξ3x3[ i, j,Nq, e]

            for n = 1:6
              sJ[i, j, n, e] = hypot(n1[i, j, n, e], n2[i, j, n, e], n3[i, j, n, e])
              n1[i, j, n, e] /= sJ[i, j, n, e]
              n2[i, j, n, e] /= sJ[i, j, n, e]
              n3[i, j, n, e] /= sJ[i, j, n, e]
            end
            =#

            ije = i + (j - 1) * Nq + (e - 1) * nface * Nq^2
            n1[ije + (1 - 1) * Nq^2] = -J[1, i, j, e] * ξ1x1[1, i, j, e]
            n1[ije + (2 - 1) * Nq^2] = J[Nq, i, j, e] * ξ1x1[Nq, i, j, e]
            n1[ije + (3 - 1) * Nq^2] = -J[i, 1, j, e] * ξ2x1[i, 1, j, e]
            n1[ije + (4 - 1) * Nq^2] = J[i, Nq, j, e] * ξ2x1[i, Nq, j, e]
            n1[ije + (5 - 1) * Nq^2] = -J[i, j, 1, e] * ξ3x1[i, j, 1, e]
            n1[ije + (6 - 1) * Nq^2] = J[i, j, Nq, e] * ξ3x1[i, j, Nq, e]
            n2[ije + (1 - 1) * Nq^2] = -J[1, i, j, e] * ξ1x2[1, i, j, e]
            n2[ije + (2 - 1) * Nq^2] = J[Nq, i, j, e] * ξ1x2[Nq, i, j, e]
            n2[ije + (3 - 1) * Nq^2] = -J[i, 1, j, e] * ξ2x2[i, 1, j, e]
            n2[ije + (4 - 1) * Nq^2] = J[i, Nq, j, e] * ξ2x2[i, Nq, j, e]
            n2[ije + (5 - 1) * Nq^2] = -J[i, j, 1, e] * ξ3x2[i, j, 1, e]
            n2[ije + (6 - 1) * Nq^2] = J[i, j, Nq, e] * ξ3x2[i, j, Nq, e]
            n3[ije + (1 - 1) * Nq^2] = -J[1, i, j, e] * ξ1x3[1, i, j, e]
            n3[ije + (2 - 1) * Nq^2] = J[Nq, i, j, e] * ξ1x3[Nq, i, j, e]
            n3[ije + (3 - 1) * Nq^2] = -J[i, 1, j, e] * ξ2x3[i, 1, j, e]
            n3[ije + (4 - 1) * Nq^2] = J[i, Nq, j, e] * ξ2x3[i, Nq, j, e]
            n3[ije + (5 - 1) * Nq^2] = -J[i, j, 1, e] * ξ3x3[i, j, 1, e]
            n3[ije + (6 - 1) * Nq^2] = J[i, j, Nq, e] * ξ3x3[i, j, Nq, e]

            for n in 1:6
                sJ[ije + (n - 1) * Nq^2] = hypot(
                    n1[ije + (n - 1) * Nq^2],
                    n2[ije + (n - 1) * Nq^2],
                    n3[ije + (n - 1) * Nq^2],
                )
                n1[ije + (n - 1) * Nq^2] /= sJ[ije + (n - 1) * Nq^2]
                n2[ije + (n - 1) * Nq^2] /= sJ[ije + (n - 1) * Nq^2]
                n3[ije + (n - 1) * Nq^2] /= sJ[ije + (n - 1) * Nq^2]
            end
        end
    end

    nothing
end

creategrid1d(elemtocoord, r) = creategrid(Val(1), elemtocoord, r)
creategrid2d(elemtocoord, r) = creategrid(Val(2), elemtocoord, r)
creategrid3d(elemtocoord, r) = creategrid(Val(3), elemtocoord, r)

"""
    creategrid(::Val{1}, elemtocoord::AbstractArray{S, 3},
               r::AbstractVector{T}) where {S, T}

Create a grid using `elemtocoord` (see [`brickmesh`](@ref)) using the 1-D `(-1,
1)` reference coordinates `r`. The element grids are filled using bilinear
interpolation of the element coordinates.

The grid is returned as a tuple of with `x1` array
"""
function creategrid(
    ::Val{1},
    e2c::AbstractArray{S, 3},
    r::AbstractVector{T},
) where {S, T}
    (d, nvert, nelem) = size(e2c)
    @assert d == 1
    Nq = length(r)
    x1 = Array{T, 2}(undef, Nq, nelem)
    creategrid!(x1, e2c, r)
    (x1 = x1,)
end

"""
  creategrid(::Val{2}, elemtocoord::AbstractArray{S, 3},
             r::AbstractVector{T}) where {S, T}

Create a 2-D tensor product grid using `elemtocoord` (see [`brickmesh`](@ref))
using the 1-D `(-1, 1)` reference coordinates `r`. The element grids are filled
using bilinear interpolation of the element coordinates.

The grid is returned as a tuple of the `x1` and `x2` arrays
"""
function creategrid(
    ::Val{2},
    e2c::AbstractArray{S, 3},
    r::AbstractVector{T},
) where {S, T}
    (d, nvert, nelem) = size(e2c)
    @assert d == 2
    Nq = length(r)
    x1 = Array{T, 3}(undef, Nq, Nq, nelem)
    x2 = Array{T, 3}(undef, Nq, Nq, nelem)
    creategrid!(x1, x2, e2c, r)
    (x1 = x1, x2 = x2)
end

"""
  creategrid(::Val{3}, elemtocoord::AbstractArray{S, 3},
             r::AbstractVector{T}) where {S, T}

Create a 3-D tensor product grid using `elemtocoord` (see [`brickmesh`](@ref))
using the 1-D `(-1, 1)` reference coordinates `r`. The element grids are filled
using bilinear interpolation of the element coordinates.

The grid is returned as a tuple of the `x1`, `x2`, `x3` arrays
"""
function creategrid(
    ::Val{3},
    e2c::AbstractArray{S, 3},
    r::AbstractVector{T},
) where {S, T}
    (d, nvert, nelem) = size(e2c)
    @assert d == 3
    Nq = length(r)
    x1 = Array{T, 4}(undef, Nq, Nq, Nq, nelem)
    x2 = Array{T, 4}(undef, Nq, Nq, Nq, nelem)
    x3 = Array{T, 4}(undef, Nq, Nq, Nq, nelem)
    creategrid!(x1, x2, x3, e2c, r)
    (x1 = x1, x2 = x2, x3 = x3)
end

"""
  computemetric(x1::AbstractArray{T, 2}, D::AbstractMatrix{T}) where T

Compute the 1-D metric terms from the element grid array `x1` using the
derivative matrix `D`. The derivative matrix `D` should be consistent with the
reference grid `r` used in [`creategrid!`](@ref).

The metric terms are returned as a 'NamedTuple` of the following arrays:

- `J` the Jacobian determinant
- `ξ1x1` derivative ∂r / ∂x1'
- `sJ` the surface Jacobian
- 'n1` outward pointing unit normal in \$x1\$-direction
"""
function computemetric(x1::AbstractArray{T, 2}, D::AbstractMatrix{T}) where {T}

    Nq = size(D, 1)
    nelem = size(x1, 2)
    nface = 2

    J = similar(x1)
    ξ1x1 = similar(x1)

    sJ = Array{T, 3}(undef, 1, nface, nelem)
    n1 = Array{T, 3}(undef, 1, nface, nelem)

    computemetric!(x1, J, ξ1x1, sJ, n1, D)

    (J = J, ξ1x1 = ξ1x1, sJ = sJ, n1 = n1)
end

"""
  computemetric(x1::AbstractArray{T, 3}, x2::AbstractArray{T, 3},
                D::AbstractMatrix{T}) where T

Compute the 2-D metric terms from the element grid arrays `x1` and `x2` using
the derivative matrix `D`. The derivative matrix `D` should be consistent with
the reference grid `r` used in [`creategrid!`](@ref).

The metric terms are returned as a 'NamedTuple` of the following arrays:

- `J` the Jacobian determinant
- `ξ1x1` derivative ∂ξ1 / ∂x1'
- `ξ2x1` derivative ∂ξ2 / ∂x1'
 - `ξ1x2` derivative ∂ξ1 / ∂x2'
 - `ξ2x2` derivative ∂ξ2 / ∂x2'
 - `sJ` the surface Jacobian
 - 'n1` outward pointing unit normal in \$x1\$-direction
 - 'n2` outward pointing unit normal in \$x2\$-direction
"""
function computemetric(
    x1::AbstractArray{T, 3},
    x2::AbstractArray{T, 3},
    D::AbstractMatrix{T},
) where {T}
    @assert size(x1) == size(x2)
    Nq = size(D, 1)
    nelem = size(x1, 3)
    nface = 4

    J = similar(x1)
    ξ1x1 = similar(x1)
    ξ2x1 = similar(x1)
    ξ1x2 = similar(x1)
    ξ2x2 = similar(x1)

    sJ = Array{T, 3}(undef, Nq, nface, nelem)
    n1 = Array{T, 3}(undef, Nq, nface, nelem)
    n2 = Array{T, 3}(undef, Nq, nface, nelem)

    computemetric!(x1, x2, J, ξ1x1, ξ2x1, ξ1x2, ξ2x2, sJ, n1, n2, D)

    (
        J = J,
        ξ1x1 = ξ1x1,
        ξ2x1 = ξ2x1,
        ξ1x2 = ξ1x2,
        ξ2x2 = ξ2x2,
        sJ = sJ,
        n1 = n1,
        n2 = n2,
    )
end

"""
    computemetric(x1::AbstractArray{T, 3}, x2::AbstractArray{T, 3},
                  D::AbstractMatrix{T}) where T

Compute the 3-D metric terms from the element grid arrays `x1`, `x2`, and `x3`
using the derivative matrix `D`. The derivative matrix `D` should be consistent
with the reference grid `r` used in [`creategrid!`](@ref).

The metric terms are returned as a 'NamedTuple` of the following arrays:

 - `J` the Jacobian determinant
 - `ξ1x1` derivative ∂ξ1 / ∂x1'
 - `ξ2x1` derivative ∂ξ2 / ∂x1'
 - `ξ3x1` derivative ∂ξ3 / ∂x1'
 - `ξ1x2` derivative ∂ξ1 / ∂x2'
 - `ξ2x2` derivative ∂ξ2 / ∂x2'
 - `ξ3x2` derivative ∂ξ3 / ∂x2'
 - `ξ1x3` derivative ∂ξ1 / ∂x3'
 - `ξ2x3` derivative ∂ξ2 / ∂x3'
 - `ξ3x3` derivative ∂ξ3 / ∂x3'
 - `sJ` the surface Jacobian
 - 'n1` outward pointing unit normal in \$x1\$-direction
 - 'n2` outward pointing unit normal in \$x2\$-direction
 - 'n3` outward pointing unit normal in \$x3\$-direction

!!! note

   The storage of the volume terms and surface terms from this function are
   slightly different. The volume terms used Cartesian indexing whereas the
   surface terms use linear indexing.
"""
function computemetric(
    x1::AbstractArray{T, 4},
    x2::AbstractArray{T, 4},
    x3::AbstractArray{T, 4},
    D::AbstractMatrix{T},
) where {T}

    @assert size(x1) == size(x2) == size(x3)
    Nq = size(D, 1)
    nelem = size(x1, 4)
    nface = 6

    J = similar(x1)
    ξ1x1 = similar(x1)
    ξ2x1 = similar(x1)
    ξ3x1 = similar(x1)
    ξ1x2 = similar(x1)
    ξ2x2 = similar(x1)
    ξ3x2 = similar(x1)
    ξ1x3 = similar(x1)
    ξ2x3 = similar(x1)
    ξ3x3 = similar(x1)

    sJ = Array{T, 3}(undef, Nq^2, nface, nelem)
    n1 = Array{T, 3}(undef, Nq^2, nface, nelem)
    n2 = Array{T, 3}(undef, Nq^2, nface, nelem)
    n3 = Array{T, 3}(undef, Nq^2, nface, nelem)

    computemetric!(
        x1,
        x2,
        x3,
        J,
        ξ1x1,
        ξ2x1,
        ξ3x1,
        ξ1x2,
        ξ2x2,
        ξ3x2,
        ξ1x3,
        ξ2x3,
        ξ3x3,
        sJ,
        n1,
        n2,
        n3,
        D,
    )

    (
        J = J,
        ξ1x1 = ξ1x1,
        ξ2x1 = ξ2x1,
        ξ3x1 = ξ3x1,
        ξ1x2 = ξ1x2,
        ξ2x2 = ξ2x2,
        ξ3x2 = ξ3x2,
        ξ1x3 = ξ1x3,
        ξ2x3 = ξ2x3,
        ξ3x3 = ξ3x3,
        sJ = sJ,
        n1 = n1,
        n2 = n2,
        n3 = n3,
    )
end

end # module


text:
tests = [
	"common",
	"linearreg",
	"regproblems",
	"solve",
	"lasso"
]

for t in tests
	f = "$t.jl"
	println("* running $f ...")
	include(f)
end


text:
"""
    MapFun <: Augmentor.Operation

Description
--------------

Maps the given function over all individual array elements.

This means that the given function is called with an individual
elements and is expected to return a transformed element that
should take the original's place. This further implies that the
function is expected to be unary. It is encouraged that the
function should be consistent with its return type and
type-stable.

Usage
--------------

    MapFun(fun)

Arguments
--------------

- **`fun`** : The unary function that should be mapped over all
    individual array elements.

See also
--------------

[`AggregateThenMapFun`](@ref), [`ConvertEltype`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor, ColorTypes
img = testpattern()

# subtract the constant RGBA value from each pixel
augment(img, MapFun(px -> px - RGBA(0.5, 0.3, 0.7, 0.0)))

# separate channels to scale each numeric element by a constant value
pl = SplitChannels() |> MapFun(el -> el * 0.5) |> CombineChannels(RGBA)
augment(img, pl)
```
"""
struct MapFun{T} <: Operation
    fun::T
end

@inline supports_lazy(::Type{<:MapFun}) = true

function applyeager(op::MapFun, img::AbstractArray, param)
    plain_array(map(op.fun, img))
end

function applylazy(op::MapFun, img::AbstractArray, param)
    mappedarray(op.fun, img)
end

function showconstruction(io::IO, op::MapFun)
    print(io, nameof(typeof(op)), '(', op.fun, ')')
end

function Base.show(io::IO, op::MapFun)
    if get(io, :compact, false)
        print(io, "Map function \"", op.fun, "\" over image")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end

# --------------------------------------------------------------------

"""
    AggregateThenMapFun <: Augmentor.Operation

Description
--------------

Compute some aggregated value of the current image using the
given function `aggfun`, and map that value over the current
image using the given function `mapfun`.

This is particularly useful for achieving effects such as
per-image normalization.

Usage
--------------

    AggregateThenMapFun(aggfun, mapfun)

Arguments
--------------

- **`aggfun`** : A function that takes the whole current image as
    input and which result will also be passed to `mapfun`. It
    should have a signature of `img -> agg`, where `img` will the
    the current image. What type and value `agg` should be is up
    to the user.

- **`mapfun`** : The binary function that should be mapped over
    all individual array elements. It should have a signature of
    `(px, agg) -> new_px` where `px` is a single element of the
    current image, and `agg` is the output of `aggfun`.

See also
--------------

[`MapFun`](@ref), [`ConvertEltype`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor
img = testpattern()

# subtract the average RGB value of the current image
augment(img, AggregateThenMapFun(img -> mean(img), (px, agg) -> px - agg))
```
"""
struct AggregateThenMapFun{A,M} <: Operation
    aggfun::A
    mapfun::M
end

@inline supports_lazy(::Type{<:AggregateThenMapFun}) = true

function applyeager(op::AggregateThenMapFun, img::AbstractArray, param)
    # promote the result to float type because `aggfun` might not widely support `Gray{N0f8}`/`N0f8` types well
    # e.g, type instability in https://github.com/JuliaGraphics/ColorVectorSpace.jl/issues/134
    agg = convert(floattype(eltype(img)), op.aggfun(img))
    plain_array(map(x -> op.mapfun(x, agg), img))
end

function applylazy(op::AggregateThenMapFun, img::AbstractArray, param)
    agg = convert(floattype(eltype(img)), op.aggfun(img))
    mappedarray(x -> op.mapfun(x, agg), img)
end

function showconstruction(io::IO, op::AggregateThenMapFun)
    print(io, nameof(typeof(op)), '(', op.aggfun, ", ", op.mapfun, ')')
end

function Base.show(io::IO, op::AggregateThenMapFun)
    if get(io, :compact, false)
        print(io, "Map result of \"", op.aggfun, "\" using \"", op.mapfun, "\" over image")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end


text:
# Multithreaded matrix-vector multiplication for 'Transpose{<:Any,<:SparseMatrixCSC}'
# using Threads.@threads
# Why the unused α and β? See: https://github.com/JuliaLang/julia/blob/381693d3dfc9b7072707f6d544f82f6637fc5e7c/stdlib/LinearAlgebra/src/matmul.jl#L202
function mul!(y::StridedVector, transA::Transpose{<:Any, <:SparseMatrixCSC}, x::StridedVector, α::Number, β::Number)
    size(transA)[2] == length(x) || throw(DimensionMismatch("transA.n != length(x)"))
    A = transA.parent
    nzval = A.nzval
    colval = A.rowval
    rowptr = A.colptr
    m = A.n
    @inbounds Threads.@threads for i = 1:m
        y[i] = 0
        @inbounds for j = rowptr[i]:(rowptr[i + 1] - 1)
            y[i] += nzval[j] * x[colval[j]]
        end
    end
    y
end
*(transA::Transpose{<:Any, <:SparseMatrixCSC}, x::StridedVector{Tx}) where {Tx} =
    (T = promote_op(matprod, eltype(transA), Tx); mul!(similar(x, T, size(transA, 1)), transA, x, true, false))

# Multithreaded matrix-vector multiplication for 'Adjoint{<:Any,<:SparseMatrixCSC}'
# using Threads.@threads
# Why the unused α and β? See: https://github.com/JuliaLang/julia/blob/381693d3dfc9b7072707f6d544f82f6637fc5e7c/stdlib/LinearAlgebra/src/matmul.jl#L202
function mul!(y::StridedVector, adjA::Adjoint{<:Any, <:SparseMatrixCSC}, x::StridedVector, α::Number, β::Number)
    size(adjA)[2] == length(x) || throw(DimensionMismatch("adjA.n != length(x)"))
    A = adjA.parent
    nzval = A.nzval
    colval = A.rowval
    rowptr = A.colptr
    m = A.n
    @inbounds Threads.@threads for i = 1:m
        y[i] = 0
        @inbounds for j = rowptr[i]:(rowptr[i + 1] - 1)
            y[i] += nzval[j] * x[colval[j]]
        end
    end
    y
end
*(transA::Adjoint{<:Any, <:SparseMatrixCSC}, x::StridedVector{Tx}) where {Tx} =
    (T = promote_op(matprod, eltype(transA), Tx); mul!(similar(x, T, size(transA, 1)), adjA, x, true, false))


text:
# I got 99 ...
problems = Dict()
problems["unconstrained"] = Dict()

include("unconstrained/from_optim.jl")
include("unconstrained/twicediffed.jl")


text:
using Documenter,NeuroAnalysis

makedocs(sitename="NeuroAnalysis.jl",modules=[NeuroAnalysis],
    pages=[])

deploydocs(repo = "github.com/Experica/NeuroAnalysis.jl.git")


text:
@testset "OneHotOneCold" begin
    ohoc_cat = OneHotOneCold(true, ['a', 'b', 'c'])
    ohoc_cont = OneHotOneCold(false, 0:3)

    @testset "Encode" begin
        @test encode(ohoc_cat, 'a') == [1, 0, 0]
        @test encode(ohoc_cat, ['a', 'c']) == [[1, 0, 0], [0, 0, 1]]
        @test_throws ArgumentError encode(ohoc_cat, 'd')

        @test encode(ohoc_cont, 0) == [1, 0, 0]
        @test encode(ohoc_cont, 1.5) == [0, 1, 0]
        @test encode(ohoc_cont, 3) == [0, 0, 1]
        @test_throws AssertionError encode(ohoc_cont, 4)
    end

    @testset "Decode" begin
        @test decode(ohoc_cat, [1, 0, 0]) == 'a'
        @test_throws AssertionError decode(ohoc_cat, [1, 0])

        @test decode(ohoc_cont, [0, 1, 0]) == (1, 2)
        @test decode(ohoc_cat, [1 1 1; 0 0 0; 0 0 0]) == ['a', 'a', 'a']
        @test decode(ohoc_cont, [1 1 1; 0 0 0; 0 0 0]) == [(0, 1), (0, 1), (0, 1)]
        arr = zeros(2, 2, 3)
        arr[:, :, 2] .= 1
        @test decode(ohoc_cat, arr) == ['b' 'b'; 'b' 'b']
        @test decode(ohoc_cont, arr) == [(1, 2) (1, 2); (1, 2) (1, 2)]
    end

end


text:
using Documenter
using BarycentricLagrange

makedocs(sitename="BarycentricLagrange.jl")

deploydocs(
    repo = "github.com/ludvigak/BarycentricLagrange.jl.git",
)


text:
using Replay

instructions = [
    """
    function greet(msg::String="Hello!")
        println(msg)
    end
    """,
    "greet()",
    """greet("Hello, World!")""",
]

replay(instructions, use_ghostwriter=true)


text:
module TestTools

using Test
using Mimi

import Mimi:
    getproperty, reset_compdefs

reset_compdefs()

#utils: prettify
@test Mimi.prettify("camelCaseBasic") == Mimi.prettify(:camelCaseBasic) == "Camel Case Basic"
@test Mimi.prettify("camelWithAOneLetterWord") == Mimi.prettify(:camelWithAOneLetterWord) == "Camel With A One Letter Word"
@test Mimi.prettify("snake_case_basic") == Mimi.prettify(:snake_case_basic) == "Snake Case Basic"
@test Mimi.prettify("_snake__case__weird_") == Mimi.prettify(:_snake__case__weird_) == "Snake Case Weird"

#utils: interpolate
stepsize = 2       # N.B. ERROR: cannot assign variable Base.step from module Main
final = 10         # N.B. ERROR: cannot assign variable Base.last from module Main
ts = 10
@test Mimi.interpolate(collect(0:stepsize:final), ts) == collect(0:stepsize/ts:final)

@defcomp Foo begin
    input = Parameter()
    intermed = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.intermed[t] = p.input
    end
end

@defcomp Bar begin
    intermed = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.output[t] = p.intermed[t]
    end
end

m = Model()
set_dimension!(m, :time, 2)
foo = add_comp!(m, Foo)
bar = add_comp!(m, Bar)

foo[:input] = 3.14
bar[:intermed] = foo[:intermed]

run(m)

end #module


text:
function petk07()
#
# M. Petkovic testing polynomials, page 147
#
    y = [1.0*[1;1;1];-2+im; -2-im; 5*im;5*im; -5*im; -5*im];
    p = poly(y);
    z = [1.0 3; -2+im 1;  -2-im 1;  5*im 2; -5*im 2];

    p, PolyZeros(z)
end


text:
module DenseHashDictionaries

import Base: @propagate_inbounds
using Dictionaries
export DenseHashIndices, DenseHashDictionary

perfect_hash(::Any) = false
perfect_hash(::Union{Bool, UInt8, UInt16, UInt32, UInt, Int8, Int16, Int32, Int, Char}) = true

struct DenseHashIndices{I} <: AbstractIndices{I}
    # The hash table
    slots::Vector{Int}

    # Densely ordered hashes and values
    hashes::Vector{UInt}
    values::Vector{I}
end

function DenseHashIndices{I}(; sizehint = 8) where {I}
    @assert sizehint > 0
    DenseHashIndices{I}(fill(0, sizehint), Vector{UInt}(), Vector{I}())
end

function rehash!(indices::DenseHashIndices{I}, newsize::Integer) where {I}
    slots = resize!(indices.slots, newsize)
    fill!(slots, 0)
    bit_mask = newsize - one(typeof(newsize)) # newsize is a power of two

    for (index, full_hash) in enumerate(indices.hashes)
        trial_slot = reinterpret(Int, full_hash) & bit_mask
        @inbounds while true
            trial_slot = (trial_slot + 1)
            if slots[trial_slot] == 0
                slots[trial_slot] = index
                break
            else
                trial_slot = trial_slot & bit_mask
            end
            # This is potentially an infinte loop and care must be taken by the callee not
            # to overfill the container
        end
    end
end

Base.length(indices::DenseHashIndices) = length(indices.values)

# Token interface
Dictionaries.istokenizable(::DenseHashIndices) = true

Dictionaries.tokentype(::DenseHashIndices) = Int

# Duration iteration the token cannot be used for deletion - we do not worry about the slots
@propagate_inbounds function Dictionaries.iteratetoken(indices::DenseHashIndices)
    if isempty(indices.values)
        return nothing
    end
    return ((0, 1), 1)
end

@propagate_inbounds function Dictionaries.iteratetoken(indices::DenseHashIndices, index::Int)
    if index == length(indices.values)
        return nothing
    end
    index = index + 1
    return ((0, index), index)
end

function Dictionaries.gettoken(indices::DenseHashIndices{I}, i::I) where {I}
    full_hash = hash(i)
    n_slots = length(indices.slots)
    bit_mask = n_slots - 1 # n_slots is always a power of two

    trial_slot = reinterpret(Int, full_hash) & bit_mask
    @inbounds while true
        trial_slot = (trial_slot + 1)
        trial_index = indices.slots[trial_slot]
        if trial_index == 0
            return (false, (0, 0))
        end

        if full_hash == indices.hashes[trial_index] && (perfect_hash(i) || isequal(i, indices.values[trial_index]))
            return (true, (trial_slot, trial_index))
        end

        trial_slot = trial_slot & bit_mask
        # This is potentially an infinte loop and care must be taken upon insertion not
        # to completely fill the container
    end
end

@propagate_inbounds function Dictionaries.gettokenvalue(indices::DenseHashIndices, (_slot, index))
    return indices.values[index]
end

# Insertion interface
Dictionaries.isinsertable(::DenseHashIndices) = true

function Dictionaries.gettoken!(indices::DenseHashIndices{I}, i::I) where {I}
    full_hash = hash(i)
    n_slots = length(indices.slots)
    bit_mask = n_slots - 1 # n_slots is always a power of two
    n_values = length(indices.values)

    trial_slot = reinterpret(Int, full_hash) & bit_mask
    trial_index = 0
    @inbounds while true
        trial_slot = (trial_slot + 1)
        trial_index = indices.slots[trial_slot]
        if trial_index <= 0
            indices.slots[trial_slot] = n_values + 1
            break
        end

        trial_hash = indices.hashes[trial_index]

        if trial_hash == full_hash && (perfect_hash(i) || isequal(i, indices.values[trial_index]))
            return (true, (trial_slot, trial_index))
        end

        trial_slot = trial_slot & bit_mask
        # This is potentially an infinte loop and care must be taken upon insertion not
        # to completely fill the container
    end

    push!(indices.hashes, full_hash)
    push!(indices.values, i)

    # Expand the hash map when it reaches 2/3rd full
    if 3 * (n_values + 1) > 2 * n_slots
        # Grow faster for small hash maps than for large ones
        rehash!(indices, n_slots > 16000 ? 2 * n_slots : 4 * n_slots)
    end

    return (false, (trial_slot, n_values + 1))
end

@propagate_inbounds function Dictionaries.deletetoken!(indices::DenseHashIndices, (slot, index))
    indices.slots[slot] = -1
    splice!(indices.hashes, index)
    splice!(indices.values, index)

    # Shrink the hash map when it is less than 1/4th full
    n_slots = length(indices.slots)
    n_values = length(indices.values)
    if 4 * n_values < n_slots && n_slots > 8
        # Halve each time
        rehash!(indices, n_slots >> 0x01)
    end

    return indices
end

# Factories

Base.empty(::DenseHashIndices, ::Type{I}) where {I} = DenseHashIndices{I}()

########

struct DenseHashDictionary{I, T} <: AbstractDictionary{I, T}
    indices::DenseHashIndices{I}
    values::Vector{T}
end

function DenseHashDictionary{I, T}(; sizehint = 8) where {I, T}
    DenseHashDictionary{I, T}(DenseHashIndices{I}(; sizehint = sizehint), Vector{T}())
end

# indices

Base.keys(dict::DenseHashDictionary) = dict.indices

# tokens

Dictionaries.tokenized(dict::DenseHashDictionary) = dict.values

# values

function Dictionaries.istokenassigned(dict::DenseHashDictionary, (_slot, index))
    return isassigned(dict.values, index)
end

@propagate_inbounds function Dictionaries.gettokenvalue(dict::DenseHashDictionary, (_slot, index))
    return dict.values[index]
end

Dictionaries.issettable(::DenseHashDictionary) = true

@propagate_inbounds function Dictionaries.settokenvalue!(dict::DenseHashDictionary{<:Any, T}, (_slot, index), value::T) where {T}
    dict.values[index] = value
    return dict
end

# insertion

function Dictionaries.gettoken!(dict::DenseHashDictionary{I}, i::I) where {I}
    (hadtoken, (slot, index)) = gettoken!(keys(dict), i)
    if !hadtoken
        resize!(dict.values, length(dict.values) + 1)
    end
    return (hadtoken, (slot, index))
end

function Dictionaries.deletetoken!(dict::DenseHashDictionary, (slot, index))
    deletetoken!(dict.indices, (slot, index))
    splice!(dict.values, index)
    return dict
end

# Factories

Base.empty(::DenseHashIndices, ::Type{I}, ::Type{T}) where {I, T} = DenseHashDictionary{I, T}()

function Base.similar(indices::DenseHashIndices{I}, ::Type{T}) where {I, T}
    return DenseHashDictionary(indices, Vector{T}(undef, length(indices)))
end

end # module


text:
@inline function writeshortest(buf::Vector{UInt8}, pos, x::T,
    plus=false, space=false, hash=true,
    precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.'), typed=false) where {T}
    @assert 0 < pos <= length(buf)
    neg = signbit(x)
    # special cases
    if x == 0
        if neg
            buf[pos] = UInt8('-')
            pos += 1
        elseif plus
            buf[pos] = UInt8('+')
            pos += 1
        elseif space
            buf[pos] = UInt8(' ')
            pos += 1
        end
        buf[pos] = UInt8('0')
        pos += 1
        if hash
            buf[pos] = decchar
            pos += 1
        end
        if precision == -1
            buf[pos] = UInt8('0')
            pos += 1
            if typed && x isa Float32
                buf[pos] = UInt8('f')
                buf[pos + 1] = UInt8('0')
                pos += 2
            end
            return pos
        end
        while precision > 1
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
        return pos
    elseif isnan(x)
        buf[pos] = UInt8('N')
        buf[pos + 1] = UInt8('a')
        buf[pos + 2] = UInt8('N')
        if typed
            if x isa Float32
                buf[pos + 3] = UInt8('3')
                buf[pos + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + 3] = UInt8('1')
                buf[pos + 4] = UInt8('6')
            end
        end
        return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    elseif !isfinite(x)
        if neg
            buf[pos] = UInt8('-')
        end
        buf[pos + neg] = UInt8('I')
        buf[pos + neg + 1] = UInt8('n')
        buf[pos + neg + 2] = UInt8('f')
        if typed
            if x isa Float32
                buf[pos + neg + 3] = UInt8('3')
                buf[pos + neg + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + neg + 3] = UInt8('1')
                buf[pos + neg + 4] = UInt8('6')
            end
        end
        return pos + neg + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    end

    bits = uint(x)
    mant = bits & (oftype(bits, 1) << mantissabits(T) - oftype(bits, 1))
    exp = Int((bits >> mantissabits(T)) & ((Int64(1) << exponentbits(T)) - 1))
    m2 = oftype(bits, Int64(1) << mantissabits(T)) | mant
    e2 = exp - bias(T) - mantissabits(T)
    fraction = m2 & ((oftype(bits, 1) << -e2) - 1)
    if e2 > 0 || e2 < -52 || fraction != 0
        if exp == 0
            e2 = 1 - bias(T) - mantissabits(T) - 2
            m2 = mant
        else
            e2 -= 2
        end
        even = (m2 & 1) == 0
        mv = oftype(m2, 4 * m2)
        mp = oftype(m2, mv + 2)
        mmShift = mant != 0 || exp <= 1
        mm = oftype(m2, mv - 1 - mmShift)
        vmIsTrailingZeros = false
        vrIsTrailingZeros = false
        lastRemovedDigit = 0x00
        if e2 >= 0
            q = log10pow2(e2) - (T == Float64 ? (e2 > 3) : 0)
            e10 = q
            k = pow5_inv_bitcount(T) + pow5bits(q) - 1
            i = -e2 + q + k
            vr, vp, vm = mulshiftinvsplit(T, mv, mp, mm, q, i)
            if T == Float32 || T == Float16
                if q != 0 && div(vp - 1, 10) <= div(vm, 10)
                    l = pow5_inv_bitcount(T) + pow5bits(q - 1) - 1
                    mul = T == Float32 ? FLOAT_POW5_INV_SPLIT[q] : HALF_POW5_INV_SPLIT[q]
                    lastRemovedDigit = (mulshift(mv, mul, -e2 + q - 1 + l) % 10) % UInt8
                end
            end
            if q <= qinvbound(T)
                if ((mv % UInt32) - 5 * div(mv, 5)) == 0
                    vrIsTrailingZeros = pow5(mv, q)
                elseif even
                    vmIsTrailingZeros = pow5(mm, q)
                else
                    vp -= pow5(mp, q)
                end
            end
        else
            q = log10pow5(-e2) - (T == Float64 ? (-e2 > 1) : 0)
            e10 = q + e2
            i = -e2 - q
            k = pow5bits(i) - pow5_bitcount(T)
            j = q - k
            vr, vp, vm = mulshiftsplit(T, mv, mp, mm, i, j)
            if T == Float32 || T == Float16
                if q != 0 && div(vp - 1, 10) <= div(vm, 10)
                    j = q - 1 - (pow5bits(i + 1) - pow5_bitcount(T))
                    mul = T == Float32 ? FLOAT_POW5_SPLIT[i + 2] : HALF_POW5_SPLIT[i + 2]
                    lastRemovedDigit = (mulshift(mv, mul, j) % 10) % UInt8
                end
            end
            if q <= 1
                vrIsTrailingZeros = true
                if even
                    vmIsTrailingZeros = mmShift
                else
                    vp -= 1
                end
            elseif q < qbound(T)
                vrIsTrailingZeros = pow2(mv, q - (T != Float64))
            end
        end
        removed = 0
        if vmIsTrailingZeros || vrIsTrailingZeros
            while true
                vpDiv10 = div(vp, 10)
                vmDiv10 = div(vm, 10)
                vpDiv10 <= vmDiv10 && break
                vmMod10 = (vm % UInt32) - UInt32(10) * (vmDiv10 % UInt32)
                vrDiv10 = div(vr, 10)
                vrMod10 = (vr % UInt32) - UInt32(10) * (vrDiv10 % UInt32)
                vmIsTrailingZeros &= vmMod10 == 0
                vrIsTrailingZeros &= lastRemovedDigit == 0
                lastRemovedDigit = vrMod10 % UInt8
                vr = vrDiv10
                vp = vpDiv10
                vm = vmDiv10
                removed += 1
            end
            if vmIsTrailingZeros
                while true
                    vmDiv10 = div(vm, 10)
                    vmMod10 = (vm % UInt32) - UInt32(10) * (vmDiv10 % UInt32)
                    vmMod10 != 0 && break
                    vpDiv10 = div(vp, 10)
                    vrDiv10 = div(vr, 10)
                    vrMod10 = (vr % UInt32) - UInt32(10) * (vrDiv10 % UInt32)
                    vrIsTrailingZeros &= lastRemovedDigit == 0
                    lastRemovedDigit = vrMod10 % UInt8
                    vr = vrDiv10
                    vp = vpDiv10
                    vm = vmDiv10
                    removed += 1
                end
            end
            if vrIsTrailingZeros && lastRemovedDigit == 5 && vr % 2 == 0
                lastRemovedDigit = UInt8(4)
            end
            output = vr + ((vr == vm && (!even || !vmIsTrailingZeros)) || lastRemovedDigit >= 5)
        else
            roundUp = false
            vpDiv100 = div(vp, 100)
            vmDiv100 = div(vm, 100)
            if vpDiv100 > vmDiv100
                vrDiv100 = div(vr, 100)
                vrMod100 = (vr % UInt32) - UInt32(100) * (vrDiv100 % UInt32)
                roundUp = vrMod100 >= 50
                vr = vrDiv100
                vp = vpDiv100
                vm = vmDiv100
                removed += 2
            end
            while true
                vpDiv10 = div(vp, 10)
                vmDiv10 = div(vm, 10)
                vpDiv10 <= vmDiv10 && break
                vrDiv10 = div(vr, 10)
                vrMod10 = (vr % UInt32) - UInt32(10) * (vrDiv10 % UInt32)
                roundUp = vrMod10 >= 5
                vr = vrDiv10
                vp = vpDiv10
                vm = vmDiv10
                removed += 1
            end
            output = vr + (vr == vm || roundUp || lastRemovedDigit >= 5)
        end
        nexp = e10 + removed
    else
        output = m2 >> -e2
        nexp = 0
        while true
            q = div(output, 10)
            r = (output % UInt32) - UInt32(10) * (q % UInt32)
            r != 0 && break
            output = q
            nexp += 1
        end
    end

    if typed && x isa Float16
        buf[pos] = UInt8('F')
        buf[pos + 1] = UInt8('l')
        buf[pos + 2] = UInt8('o')
        buf[pos + 3] = UInt8('a')
        buf[pos + 4] = UInt8('t')
        buf[pos + 5] = UInt8('1')
        buf[pos + 6] = UInt8('6')
        buf[pos + 7] = UInt8('(')
        pos += 8
    end
    if neg
        buf[pos] = UInt8('-')
        pos += 1
    elseif plus
        buf[pos] = UInt8('+')
        pos += 1
    elseif space
        buf[pos] = UInt8(' ')
        pos += 1
    end

    olength = decimallength(output)
    exp_form = true
    pt = nexp + olength
    if -4 < pt <= (precision == -1 ? (T == Float16 ? 5 : 6) : precision)
        exp_form = false
        if pt <= 0
            buf[pos] = UInt8('0')
            pos += 1
            buf[pos] = decchar
            pos += 1
            for _ = 1:abs(pt)
                buf[pos] = UInt8('0')
                pos += 1
            end
        # elseif pt >= olength
            # nothing to do at this point
        # else
            # nothing to do at this point
        end
    else
        pos += 1
    end
    i = 0
    ptr = pointer(buf)
    ptr2 = pointer(DIGIT_TABLE)
    if (output >> 32) != 0
        q = output ÷ 100000000
        output2 = (output % UInt32) - UInt32(100000000) * (q % UInt32)
        output = q

        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        d = output2 % UInt32(10000)
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        d0 = (d % 100) << 1
        d1 = (d ÷ 100) << 1
        memcpy(ptr, pos + olength - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - 4, ptr2, c1 + 1, 2)
        memcpy(ptr, pos + olength - 6, ptr2, d0 + 1, 2)
        memcpy(ptr, pos + olength - 8, ptr2, d1 + 1, 2)
        i += 8
    end
    output2 = output % UInt32
    while output2 >= 10000
        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        memcpy(ptr, pos + olength - i - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - i - 4, ptr2, c1 + 1, 2)
        i += 4
    end
    if output2 >= 100
        c = (output2 % UInt32(100)) << 1
        output2 = div(output2, UInt32(100))
        memcpy(ptr, pos + olength - i - 2, ptr2, c + 1, 2)
        i += 2
    end
    if output2 >= 10
        c = output2 << 1
        buf[pos + 1] = DIGIT_TABLE[c + 2]
        buf[pos - exp_form] = DIGIT_TABLE[c + 1]
    else
        buf[pos - exp_form] = UInt8('0') + (output2 % UInt8)
    end

    if !exp_form
        if pt <= 0
            pos += olength
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        elseif pt >= olength
            pos += olength
            precision -= olength
            for _ = 1:nexp
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
            if hash
                buf[pos] = decchar
                pos += 1
                if precision < 0
                    buf[pos] = UInt8('0')
                    pos += 1
                end
                while precision > 0
                    buf[pos] = UInt8('0')
                    pos += 1
                    precision -= 1
                end
            end
        else
            pointoff = olength - abs(nexp)
            memmove(ptr, pos + pointoff + 1, ptr, pos + pointoff, olength - pointoff + 1)
            buf[pos + pointoff] = decchar
            pos += olength + 1
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
    else
        if olength > 1 || hash
            buf[pos] = decchar
            pos += olength
            precision -= olength
        end
        if hash && olength == 1
            buf[pos] = UInt8('0')
            pos += 1
        end
        while hash && precision > 0
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end

        buf[pos] = expchar
        pos += 1
        exp2 = nexp + olength - 1
        if exp2 < 0
            buf[pos] = UInt8('-')
            pos += 1
            exp2 = -exp2
        elseif padexp
            buf[pos] = UInt8('+')
            pos += 1
        end

        if exp2 >= 100
            c = exp2 % 10
            memcpy(ptr, pos, ptr2, 2 * div(exp2, 10) + 1, 2)
            buf[pos + 2] = UInt8('0') + (c % UInt8)
            pos += 3
        elseif exp2 >= 10
            memcpy(ptr, pos, ptr2, 2 * exp2 + 1, 2)
            pos += 2
        else
            if padexp
                buf[pos] = UInt8('0')
                pos += 1
            end
            buf[pos] = UInt8('0') + (exp2 % UInt8)
            pos += 1
        end
    end
    if typed && x isa Float16
        buf[pos] = UInt8(')')
        pos += 1
    end

    return pos
end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: milesfrain/julia

max_forks_repo_name: milesfrain/julia


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

## original example, following here

# custom shape wrapped in a function
function shape_curve_rectangle(cr::CairoContext,x0::Real,y0::Real,
    rect_width::Real, rect_height::Real, radius::Real)

    save(cr);

    x1=x0+rect_width;
    y1=y0+rect_height;

    if (rect_width == 0 || rect_height == 0)
        return;
        end
    if (rect_width/2 < radius)
        if (rect_height/2 < radius)
            move_to(cr, x0, (y0 + y1)/2);
            curve_to(cr, x0 ,y0, x0, y0, (x0 + x1)/2, y0);
            curve_to(cr, x1, y0, x1, y0, x1, (y0 + y1)/2);
            curve_to(cr, x1, y1, x1, y1, (x1 + x0)/2, y1);
            curve_to(cr, x0, y1, x0, y1, x0, (y0 + y1)/2);
        else
            move_to(cr, x0, y0 + radius);
            curve_to(cr, x0 ,y0, x0, y0, (x0 + x1)/2, y0);
            curve_to(cr, x1, y0, x1, y0, x1, y0 + radius);
            line_to(cr, x1 , y1 - radius);
            curve_to(cr, x1, y1, x1, y1, (x1 + x0)/2, y1);
            curve_to(cr, x0, y1, x0, y1, x0, y1- radius);
        end

    else
        if (rect_height/2 < radius)
            move_to(cr, x0, (y0 + y1)/2);
            curve_to(cr, x0 , y0, x0 , y0, x0 + radius, y0);
            line_to(cr, x1 - radius, y0);
            curve_to(cr, x1, y0, x1, y0, x1, (y0 + y1)/2);
            curve_to(cr, x1, y1, x1, y1, x1 - radius, y1);
            line_to(cr, x0 + radius, y1);
            curve_to(cr, x0, y1, x0, y1, x0, (y0 + y1)/2);
        else
            move_to(cr, x0, y0 + radius);
            curve_to(cr, x0 , y0, x0 , y0, x0 + radius, y0);
            line_to(cr, x1 - radius, y0);
            curve_to(cr, x1, y0, x1, y0, x1, y0 + radius);
            line_to(cr, x1 , y1 - radius);
            curve_to(cr, x1, y1, x1, y1, x1 - radius, y1);
            line_to(cr, x0 + radius, y1);
            curve_to(cr, x0, y1, x0, y1, x0, y1- radius);
        end
    end
    close_path(cr);

    set_source_rgb(cr, 0.5, 0.5, 1);
    fill_preserve(cr);
    set_source_rgba(cr, 0.5, 0, 0, 0.5);
    set_line_width(cr, 10.0);
    stroke(cr);
    restore(cr);
    return;
    end

shape_curve_rectangle(cr,25.6,25.6,204.8,204.8,102.4);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_curve_rectangle.png");


text:
"""
Return n indices in order
Cannot find a julia method which does, this.
If a faster method is found, replace it later.
"""
function indminn!(ivec, x::AbstractVector{Float64}, val::Float64)
	# using enumerate to avoid indexing
	n=length(ivec)
	if(length(x)<n)
		# if the desired indices is greater than length of x, what to do..
		fill!(ivec,one(eltype(ivec)))
	else
		fill!(ivec,zero(eltype(ivec)))
		for inn in 1:n
			ivec[inn] = indminimum(x, val, ivec)
		end
		sort!(ivec)
	end
	return ivec
end

"""
minimum index using imask
"""
function indminimum(x, val, imask=[])
	min_i=0
	min_x=typemax(Float64)
	for (i, xi) in enumerate(x)
		dist = abs(xi - val)
		if ((dist < min_x) && (i ∉ imask))
			min_x = dist
			min_i = i
		end
	end
	return min_i
end

function indminn(x::AbstractVector{Float64}, val::Float64, n::Int64=1)
	ivec=fill(0,n)
	indminn!(ivec,x,val)
	return ivec
end

"return index such that "
function indminn_inside(x, valbeg, valend)
	iibeg1 = indminimum(x, valbeg)
	iibeg2 = indminimum(x, valbeg,iibeg1)
	iiend1 = indminimum(x, valend)
	iiend2 = indminimum(x, valend,iiend1)

	if((valbeg-x[iibeg1])*(valend-x[iibeg1]) <= 0.0)
		iibeg=iibeg1
	else
		iibeg=iibeg2
	end
	if((valbeg-x[iiend1])*(valend-x[iiend1])<=0.0)
		iiend=iiend1
	else
		iiend=iiend2
	end

	return iibeg, iiend
end


text:
import NaiveNASlib
import InteractiveUtils:subtypes
using Test

@testset "Basic vertex tests" begin

    @testset "Method contracts $subtype" for subtype in implementations(AbstractVertex)
        @test hasmethod(inputs, (subtype,))
        @test hasmethod(clone, (subtype, Vararg{AbstractVertex}))
    end

    @testset "InputVertex tests" begin
        iv1 = InputVertex(1)
        iv2 = InputVertex("name")

        @test iv1.name == 1
        @test iv2.name == "name"

        @test iv1 == iv1
        @test iv2 == iv2
        @test iv1 != iv2
        @test iv1 == clone(iv1)
        @test iv2 == clone(iv2)
        @test_throws ErrorException clone(iv1, iv1)
    end

    @testset "CompVertex tests" begin

        # Helper function to create vertices with one input
        cv(f, i=InputVertex(1)) = CompVertex(f, i)
        @testset "Scalar operations" begin
            @test cv(x -> 5)(2) == 5
            @test cv(x -> 5x)(2) == 10
            @test cv(sum)([2,3]) == 5
            @test cv(+)(3,4) == 7
            @test cv(vcat)(5,6) == [5,6]
        end

        @testset "Array operations" begin
            @test cv(x -> x .* 5)(ones(2,3)) == [5 5 5; 5 5 5]
            @test cv(x -> [3 4] * x)(ones(2,3)) == [7.0 7.0 7.0]
        end

        @testset "Copy" begin
            iv = InputVertex(1)
            cv_orig = cv(x -> 3x, iv)
            cv_copy = clone(cv_orig, clone(iv, ))

            @test issame(cv_orig, cv_copy)
            @test cv_orig(3) == cv_copy(3)
        end
    end

    @testset "Pretty printing" begin

        @testset "Show CompVertex" begin
            cv = CompVertex(+, InputVertex.(1:2))
            @test showstr(show, cv) == "CompVertex(+), inputs=[InputVertex(1), InputVertex(2)]"
        end

        @testset "Info string CompVertex" begin
            cv = CompVertex(+, InputVertex.(["input1", "input2"]))

            @test infostr(RawInfoStr(), cv) == "CompVertex(+), inputs=[input1, input2]"
            @test infostr(NameInfoStr(), cv) == "CompVertex"
            @test infostr(MutationTraitInfoStr(), cv) == "Immutable()"
            @test infostr(InputsInfoStr(), cv) == "[input1, input2]"
            @test infostr(NameAndInputsInfoStr(), cv) == "CompVertex, inputs=[input1, input2]"
        end

    end
end


text:
function make_text(
    args...;
    body = "stdtxt",
    size = 200,
    pos = O,
    text_frames = 1:100,
    action_frames = 1:100,
    valign = :center,
    halign = :baseline,
    kwargs...,
)
    ob = Object(
        text_frames,
        (args...) -> draw_text(
            body = body,
            size = size,
            pos = pos,
            valign = valign,
            halign = halign,
        ),
    )
    act!(ob, Action(action_frames, appear(:draw_text)))
    ob
end

function make_lamp(args...; pos, r, frames = 1:100, show_drawing = false)
    ob = Object(
        frames,
        (args...) -> draw_lamp(pos - r, r, max_frames, show_drawing = show_drawing),
        pos,
    )
    return ob
end

function make_lens(args...; pos, size, frames = 1:100)
    ob = Object(frames, (args...) -> draw_lens(pos = pos, size = size))
    return ob
end

function make_doc(args...; pos, sizey, sizex = sizey / φ, frames = 1:100, start_pos = O)
    ob = Object(
        frames,
        (args...) -> draw_doc(pos = pos, sizex = sizex, sizey = sizey),
        start_pos,
    )
    return ob
end

function make_banner(
    args...;
    pos,
    sizex,
    sizey = sizex / φ,
    frames = 1:100,
    words = ["NEW", "IDEAS"],
    rng = 42,
)
    Random.seed!(rng)
    ob = Object(
        frames,
        (args...) -> draw_banner(
            pos = pos,
            sizex = sizex,
            sizey = sizey,
            words = words,
            rng = rng,
        ),
    )
    return ob
end


text:
# Run package tests
# println("Testing Silo.jl in Julia version ", VERSION)

@static if VERSION >= v"0.7-"
  using Test
else
  using Base.Test
  range{T}(start::T; stop::T=one(T), length::Int = 100) = linspace(start, stop, length)
end
using Silo
using Compat

@testset begin
    @test string(Silo._VERSION()) == Silo._VERSION_STRING()
    include("test1dwriteInt.jl")
    # include("test1dreadwrite.jl")
end