text:
using .CApi

#-------------------------------------------------------------------------------
# FDBError
#-------------------------------------------------------------------------------
struct FDBError
    code::fdb_error_t
    desc::String
end

strerrdesc(cs::Cstring) = (cs == C_NULL) ? "unknown" : unsafe_string(cs)
FDBError(code::fdb_error_t) = FDBError(code, strerrdesc(fdb_get_error(code)))
function FDBError(future::fdb_future_ptr_t)
    desc = Ref{Cstring}(C_NULL)
    block_until(future)
    errcode = fdb_future_get_error(future, desc)
    if desc[] == C_NULL
        # try getting it with a different API
        desc[] = fdb_get_error(errcode)
    end
    FDBError(errcode, strerrdesc(desc[]))
end

function throw_on_error(result)
    (result == 0) || throw(FDBError(result))
    nothing
end

function throw_on_error(future::fdb_future_ptr_t)
    err = FDBError(future)
    if err.code != 0
        fdb_future_destroy(future)
        throw(err)
    end
    nothing
end

function show(io::IO, err::FDBError)
    code = Int(err.code)
    print(io, "FDB error ", code, " - ", err.desc)
end

#-------------------------------------------------------------------------------
# Working with futures and starting/stopping network
#-------------------------------------------------------------------------------
"""
An opaque type that represents a Future in the FoundationDB C API.

Most functions in the FoundationDB API are asynchronous, meaning that they may
return to the caller before actually delivering their result. These functions
always return FDBFuture*. An FDBFuture object represents a result value or error
to be delivered at some future time. You can wait for a Future to be “ready” –
to have a value or error delivered – by setting a callback function, or by
blocking a thread, or by polling. Once a Future is ready, you can extract either
an error code or a value of the appropriate type (the documentation for the
original function will tell you which fdb_future_get_*() function you should
call).
"""
const FDBFuture = Ref{fdb_future_ptr_t}

cancel(future::FDBFuture) = fdb_future_cancel(future[])

function block_until(future::fdb_future_ptr_t)
    Bool(fdb_future_is_ready(future)) && return
    wait_task = @async fdb_future_block_until_ready_in_thread(future)
    fetch(wait_task)
    future
end
block_until(errcode) = errcode

function with_err_check(on_success, result::Union{fdb_future_ptr_t,fdb_error_t}, on_error=throw_on_error)
    err = FDBError(result)
    (err.code == 0) ? on_success(result) : on_error(result)
end

function err_check(result::Union{fdb_future_ptr_t,fdb_error_t}, on_error=throw_on_error)
    with_err_check((x)->x, result, on_error)
end

"""
Encapsulates starting and stopping of the FDB Network
"""
struct FDBNetwork
    addr::String
    version::Cint
    task::Task

    function FDBNetwork(addr::String="127.0.0.1:4500", version::Cint=FDB_API_VERSION)
        throw_on_error(fdb_select_api_version(version))
        throw_on_error(fdb_setup_network(addr))
        network_task = @async fdb_run_network_in_thread()
        network = new(addr, version, network_task)
    end
end

const network = Ref{Union{FDBNetwork,Nothing}}(nothing)

is_client_running() = (network[] !== nothing) && !istaskdone((network[]).task)

function start_client(addr::String="127.0.0.1:4500", version::Cint=FDB_API_VERSION)
    if network[] === nothing
        network[] = FDBNetwork(addr, version)
    elseif istaskdone((network[]).task)
        error("Client stopped. Can only start one client in the lifetime of a process.")
    end
    nothing
end

function stop_client()
    if network[] !== nothing
        if !istaskdone((network[]).task)
            fdb_stop_network()
            fetch((network[]).task)
        end
    end
    nothing
end

#-------------------------------------------------------------------------------
# Opening cluster, database, transaction
#-------------------------------------------------------------------------------

"""
An opaque type that represents a Cluster in the FoundationDB C API.
"""
mutable struct FDBCluster
    cluster_file::String
    ptr::fdb_cluster_ptr_t

    function FDBCluster(cluster_file::String="/etc/foundationdb/fdb.cluster")
        cluster = new(cluster_file, C_NULL)
        finalizer((cluster)->close(cluster), cluster)
        cluster
    end
end

function show(io::IO, cluster::FDBCluster)
    print(io, "FDBCluster(", cluster.cluster_file, ") - ", (cluster.ptr == C_NULL) ? "closed" : "open")
end

open(fn::Function, cluster::FDBCluster) = try fn(open(cluster)) finally close(cluster) end
function open(cluster::FDBCluster)
    if !isopen(cluster)
        @assert is_client_running()
        handle = with_err_check(fdb_create_cluster(cluster.cluster_file)) do future
            h = Ref{fdb_cluster_ptr_t}(C_NULL)
            fdb_future_get_cluster(future, h)
            fdb_future_destroy(future)
            h[]
        end
        cluster.ptr = handle
    end
    cluster
end

function close(cluster::FDBCluster)
    if is_client_running() && isopen(cluster)
        fdb_cluster_destroy(cluster.ptr)
        cluster.ptr = C_NULL
    end
    nothing
end

"""
An opaque type that represents a database in the FoundationDB C API.

An FDBDatabase represents a FoundationDB database - a mutable, lexicographically
ordered mapping from binary keys to binary values. Modifications to a database
are performed via transactions.
"""
mutable struct FDBDatabase
    cluster::FDBCluster
    name::String
    ptr::fdb_database_ptr_t

    function FDBDatabase(cluster::FDBCluster, name::String="DB")
        db = new(cluster, name, C_NULL)
        finalizer((db)->close(db), db)
        db
    end
end

function show(io::IO, db::FDBDatabase)
    print(io, "FDBDatabase(", db.name, ") - ", (db.ptr == C_NULL) ? "closed" : "open")
end

open(fn::Function, db::FDBDatabase) = try fn(open(db)) finally close(db) end
function open(db::FDBDatabase)
    if !isopen(db)
        @assert is_client_running()
        cl = db.cluster.ptr
        name = convert(Vector{UInt8}, codeunits(db.name))
        lname = Cint(length(name))
        handle = with_err_check(fdb_cluster_create_database(cl, name, lname)) do future
            h = Ref{fdb_database_ptr_t}(C_NULL)
            fdb_future_get_database(future, h)
            fdb_future_destroy(future)
            h[]
        end
        db.ptr = handle
    end
    db
end

function close(db::FDBDatabase)
    if is_client_running() && isopen(db)
        fdb_database_destroy(db.ptr)
        db.ptr = C_NULL
    end
    nothing
end

"""
An opaque type that represents a transaction in the FoundationDB C API.

In FoundationDB, a transaction is a mutable snapshot of a database. All read and
write operations on a transaction see and modify an otherwise-unchanging version
of the database and only change the underlying database if and when the
transaction is committed. Read operations do see the effects of previous write
operations on the same transaction. Committing a transaction usually succeeds in
the absence of conflicts.

Applications must provide error handling and an appropriate retry loop around
the application code for a transaction. See the documentation for
fdb_transaction_on_error().

Transactions group operations into a unit with the properties of atomicity,
isolation, and durability. Transactions also provide the ability to maintain an
application’s invariants or integrity constraints, supporting the property of
consistency. Together these properties are known as ACID.

Transactions are also causally consistent: once a transaction has been
successfully committed, all subsequently created transactions will see the
modifications made by it.
"""
mutable struct FDBTransaction
    db::FDBDatabase
    ptr::fdb_transaction_ptr_t
    needscommit::Bool
    versionstamp::Union{Nothing,Vector{UInt8}}
    autocommit::Bool
    trackversionstamp::Bool
    versionstampfuture::Union{Nothing,fdb_future_ptr_t}

    function FDBTransaction(db::FDBDatabase; autocommit::Bool=true, trackversionstamp::Bool=false)
        tran = new(db, C_NULL, false, nothing, autocommit, trackversionstamp, nothing)
        finalizer((tran)->close(tran), tran)
        tran
    end
end

function show(io::IO, tran::FDBTransaction)
    print(io, "FDBTransaction - ", (tran.ptr == C_NULL) ? "closed" : "open")
end

function open(fn::Function, tran::FDBTransaction)
    try
        retry = true
        result = nothing
        while retry
            result = fn(open(tran))
            retry = tran.autocommit && tran.needscommit && !commit(tran)
        end
        result
    finally
        close(tran)
    end
end
function open(tran::FDBTransaction)
    if !isopen(tran)
        @assert is_client_running()
        db = tran.db.ptr
        h = Ref{fdb_database_ptr_t}(C_NULL)
        handle = with_err_check(fdb_database_create_transaction(db, h)) do result
            h[]
        end
        tran.ptr = handle
    end
    tran
end

function close(tran::FDBTransaction)
    if is_client_running() && isopen(tran)
        fdb_transaction_destroy(tran.ptr)
        tran.ptr = C_NULL
    end
    nothing
end

"""
Check if it is open.
"""
isopen(x::Union{FDBCluster,FDBDatabase,FDBTransaction}) = !(x.ptr === C_NULL)

#-------------------------------------------------------------------------------
# Transaction Ops
#-------------------------------------------------------------------------------

function reset(tran::FDBTransaction)
    fdb_transaction_reset(tran.ptr)
    nothing
end

function cancel(tran::FDBTransaction)
    fdb_transaction_cancel(tran.ptr)
    tran.needscommit = false
    nothing
end

"""
- returns true on success
- returns false on a retryable error
- throws error on non-retryable error
"""
function retry_on_error(tran::FDBTransaction, future::fdb_future_ptr_t)
    err = FDBError(future)
    if err.code != 0
        throw_on_error(fdb_transaction_on_error(tran.ptr, err.code))
        return false
    end
    true
end

function commit(tran::FDBTransaction, on_error=(result)->retry_on_error(tran,result))
    ret = err_check(fdb_transaction_commit(tran.ptr), on_error)

    if (!isa(ret, Bool) || ret) && tran.trackversionstamp && (tran.versionstampfuture !== nothing)
        keyptr = Ref{Ptr{UInt8}}(C_NULL)
        keylen = Ref{Cint}(0)
        err_check(fdb_future_get_key(tran.versionstampfuture, keyptr, keylen))
        tran.versionstamp = copy(unsafe_wrap(Array, keyptr[], (keylen[],), own=false))
        fdb_future_destroy(tran.versionstampfuture)
        tran.versionstampfuture = nothing
    end

    if isa(ret, Bool)
        if ret
            tran.needscommit = false
        end
        ret
    else
        tran.needscommit = false
        true
    end
end

function get_read_version(tran::FDBTransaction)
    ver = Ref{Int64}(0)
    with_err_check(fdb_transaction_get_read_version(tran.ptr)) do result
        fdb_future_get_version(result, ver)
        fdb_future_destroy(result)
    end
    ver[]
end

function set_read_version(tran::FDBTransaction, version)
    fdb_transaction_set_read_version(tran.ptr, Int64(version))
end

function get_committed_version(tran::FDBTransaction)
    ver = Ref{Int64}(0)
    err_check(fdb_transaction_get_committed_version(tran.ptr, ver))
    ver[]
end

#-------------------------------------------------------------------------------
# Get Set Ops
#-------------------------------------------------------------------------------
struct fdb_key_sel_t
    last_less_than::Symbol
    last_less_or_equal::Symbol
    first_greater_than::Symbol
    first_greater_or_equal::Symbol
end
const FDBKeySel = fdb_key_sel_t(:last_less_than, :last_less_or_equal, :first_greater_than, :first_greater_or_equal)

function keysel(mode::Symbol, key::Vector{UInt8})
    len = Cint(length(key))
    (mode === FDBKeySel.last_less_than)         ? (key, len, Cint(0), Cint(0)) :
    (mode === FDBKeySel.last_less_or_equal)     ? (key, len, Cint(1), Cint(0)) :
    (mode === FDBKeySel.first_greater_than)     ? (key, len, Cint(1), Cint(1)) :
    (mode === FDBKeySel.first_greater_or_equal) ? (key, len, Cint(0), Cint(1)) :
    error("unknown key selector $mode")
end

copyval(tran::FDBTransaction, present::Bool, val::Vector{UInt8}) = present ? copy(val) : nothing

function clearkey(tran::FDBTransaction, key::Vector{UInt8})
    ret = fdb_transaction_clear(tran.ptr, key, Cint(length(key)))
    tran.needscommit = true
    ret
end

function clearkeyrange(tran::FDBTransaction, begin_key::Vector{UInt8}, end_key::Vector{UInt8})
    ret = fdb_transaction_clear_range(tran.ptr, begin_key, Cint(length(begin_key)), end_key, Cint(length(end_key)))
    tran.needscommit = true
    ret
end

getkey(tran::FDBTransaction, keyselector) = getkey(copyval, tran, keyselector)
function getkey(fn::Function, tran::FDBTransaction, keyselector)
    keyptr = Ref{Ptr{UInt8}}(C_NULL)
    keylen = Ref{Cint}(0)
    with_err_check(fdb_transaction_get_key(tran.ptr, keyselector...)) do result
        err_check(fdb_future_get_key(result, keyptr, keylen))
        key = fn(tran, true, unsafe_wrap(Array, keyptr[], (keylen[],), own=false))
        fdb_future_destroy(result)
        key
    end::Vector{UInt8}
end

getval(tran::FDBTransaction, key::Vector{UInt8}) = getval(copyval, tran, key)
function getval(fn::Function, tran::FDBTransaction, key::Vector{UInt8})
    val = nothing
    present = Ref{fdb_bool_t}(false)
    valptr = Ref{Ptr{UInt8}}(C_NULL)
    vallen = Ref{Cint}(0)
    with_err_check(fdb_transaction_get(tran.ptr, key, Cint(length(key)))) do result
        err_check(fdb_future_get_value(result, present, valptr, vallen))
        val = fn(tran, Bool(present[]), unsafe_wrap(Array, valptr[], (vallen[],), own=false))
        fdb_future_destroy(result)
    end
    val
end

copykeyval(tran::FDBTransaction, key::Vector{UInt8}, val::Vector{UInt8}) = (copy(key),copy(val))
function copykeyval(tran::FDBTransaction, keyval::fdb_kv_t)
    key = unsafe_wrap(Array, convert(Ptr{UInt8}, keyval.key), (keyval.key_length,), own=false)
    val = unsafe_wrap(Array, convert(Ptr{UInt8}, keyval.value), (keyval.value_length,), own=false)
    copykeyval(tran, key, val)
end
collectkeyval(tran::FDBTransaction, keyval::fdb_kv_t, into) = push!(into, copykeyval(tran, keyval))
function collectkeyval(tran::FDBTransaction, keyvals::Vector{fdb_kv_t}, into)
    for keyval in keyvals
        collectkeyval(tran, keyval, into)
    end
    into
end
function getrange(tran::FDBTransaction, begin_keysel, end_keysel, into=Vector{Tuple{Vector{UInt8},Vector{UInt8}}}(); kwargs...)
    getrange(tran, begin_keysel, end_keysel; kwargs...) do tran, keyval
        collectkeyval(tran, keyval, into)
    end
end
function getrange(fn::Function, tran::FDBTransaction, begin_keysel, end_keysel; limit::Integer=0, target_bytes::Integer=0,
        streaming_mode::Cint=FDBStreamingMode.WANT_ALL, iteration::Integer=0, snapshot::Integer=0, reverse::Bool=false)
    more = Ref{fdb_bool_t}(false)
    count = Ref{Cint}(0)
    kvptr = Ref{Ptr{Nothing}}(C_NULL)

    fnoutput = with_err_check(fdb_transaction_get_range(tran.ptr,
        begin_keysel[1], begin_keysel[2], begin_keysel[3], begin_keysel[4],
        end_keysel[1], end_keysel[2], end_keysel[3], end_keysel[4],
        Cint(limit), Cint(target_bytes), Cint(streaming_mode), Cint(iteration), Cint(snapshot), fdb_bool_t(reverse ? 1 : 0))) do result

        fnoutput = nothing
        err_check(fdb_future_get_keyvalue_array(result, kvptr, count, more))
        if count[] > 0
            total_size = (sizeof(Ptr{Nothing}) + sizeof(Cint)) * 2 * count[]
            ptr = convert(Ptr{UInt8}, kvptr[])
            valarr = Vector{fdb_kv_t}()
            # the fields are byte packed
            for idx in 1:(count[])
                _key = unsafe_load(convert(Ptr{Ptr{Nothing}}, ptr))
                ptr += sizeof(Ptr{Nothing})
                _key_length = unsafe_load(convert(Ptr{Cint}, ptr))
                ptr += sizeof(Cint)
                _value = unsafe_load(convert(Ptr{Ptr{Nothing}}, ptr))
                ptr += sizeof(Ptr{Nothing})
                _value_length = unsafe_load(convert(Ptr{Cint}, ptr))
                ptr += sizeof(Cint)
                push!(valarr, fdb_kv_t(_key, _key_length, _value, _value_length))
            end
            fnoutput = fn(tran, valarr)
        end
        fdb_future_destroy(result)
        fnoutput
    end
    fnoutput, Bool(more[])
end

function setval(tran::FDBTransaction, key::Vector{UInt8}, val::Vector{UInt8})
    ret = fdb_transaction_set(tran.ptr, key, Cint(length(key)), val, Cint(length(val)))
    tran.needscommit = true
    ret
end

function watchkey_internal(fn::Function, tran::FDBTransaction, key::Vector{UInt8}, on_finish::Function, handle::Union{Nothing,FDBFuture})
    future = fdb_transaction_watch(tran.ptr, key, Cint(length(key)))
    (handle === nothing) || (handle[] = future)
    tran.needscommit = true
    fn(tran, key)
    err_check(future)
    nothing
end

function watchkey(tran::FDBTransaction, key::Vector{UInt8}; on_finish::Function=err_check, handle::Union{Nothing,FDBFuture}=nothing)
    watchstarted = Future()
    watch_task = @async watchkey_internal(tran, key, on_finish, handle) do tran, key
        put!(watchstarted, true)
    end
    wait(watchstarted)
    watch_task
end

function atomic(tran::FDBTransaction, key::Vector{UInt8}, param::Vector{UInt8}, op::Cint)
    tran.needscommit = true
    fdb_transaction_atomic_op(tran.ptr, key, Cint(length(key)), param, Cint(length(param)), op)
end

function atomic_add(tran::FDBTransaction, key::Vector{UInt8}, param::Integer)
    leparam = htol(param) # ensure param is little-endian, unsigned or signed in two's complement format
    param_bytes = convert(Vector{UInt8}, reinterpret(UInt8, [leparam]))
    atomic(tran, key, param_bytes, FDBMutationType.ADD)
end

function atomic_max(tran::FDBTransaction, key::Vector{UInt8}, param::Cuint)
    param_bytes = convert(Vector{UInt8}, reinterpret(UInt8, [param]))
    atomic(tran, key, param_bytes, FDBMutationType.MAX)
end

function atomic_min(tran::FDBTransaction, key::Vector{UInt8}, param::Cuint)
    param_bytes = convert(Vector{UInt8}, reinterpret(UInt8, [param]))
    atomic(tran, key, param_bytes, FDBMutationType.MIN)
end

"""
Loads the value read as an integer type.

Suitable for reading back values for keys used with `atomic_add`, `atomic_max`, `atomic_min` operations.
"""
atomic_integer(::Type{T}, val::Vector{UInt8}) where {T <: Integer} = unsafe_load(convert(Ptr{T}, pointer(val)))

atomic_and(tran::FDBTransaction, key::Vector{UInt8}, param::Vector{UInt8}) = atomic(tran, key, param, FDBMutationType.AND)
atomic_or(tran::FDBTransaction,  key::Vector{UInt8}, param::Vector{UInt8}) = atomic(tran, key, param, FDBMutationType.OR)
atomic_xor(tran::FDBTransaction, key::Vector{UInt8}, param::Vector{UInt8}) = atomic(tran, key, param, FDBMutationType.XOR)
atomic_max(tran::FDBTransaction, key::Vector{UInt8}, param::Vector{UInt8}) = atomic(tran, key, param, FDBMutationType.BYTE_MAX)
atomic_min(tran::FDBTransaction, key::Vector{UInt8}, param::Vector{UInt8}) = atomic(tran, key, param, FDBMutationType.BYTE_MIN)
function prep_atomic_key!(key::Vector{UInt8}, versionpos::Integer)
    pos = htol(UInt16(versionpos-1))  # make it 0 based index
    pos_bytes = reinterpret(UInt8, [pos])
    splice!(key, versionpos:(versionpos-1), zeros(UInt8, 10))
    splice!(key, (length(key)+1):length(key), pos_bytes)
    key
end
function atomic_setval(tran::FDBTransaction, key::Vector{UInt8}, val::Vector{UInt8}, op::Cint)
    @assert (op == FDBMutationType.SET_VERSIONSTAMPED_KEY) || (op == FDBMutationType.SET_VERSIONSTAMPED_VALUE)
    atomic(tran, key, val, op)
    if tran.trackversionstamp && (tran.versionstampfuture === nothing)
        tran.versionstampfuture = fdb_transaction_get_versionstamp(tran.ptr)
    end
    nothing
end

function conflict(tran::FDBTransaction, begin_key::Vector{UInt8}, end_key::Vector{UInt8}, conflict_type::fdb_conflict_range_type_t)
    err_check(fdb_transaction_add_conflict_range(tran.ptr, begin_key, Cint(length(begin_key)), end_key, Cint(length(end_key)), conflict_type))
end


text:
module ODMXMLTools
    using  EzXML, DataFrames, AbstractTrees
    #using Dates
    import AbstractTrees: children
    import Base: show, ht_keyindex

    export importxml,
    metadatalist,
    findclinicaldata,
    findstudy,
    findstudymetadata,
    findelement,
    findallelements,
    buildmetadata,
    eventlist,
    formlist,
    itemgrouplist,
    itemlist,
    itemgroupcontent,
    validateodm,
    clinicaldatatable,
    clinicaldatalist

    StrOrSym = Union{String, Symbol}

    const FILETYPE    = Set(["Snapshot", "Transactional"])
    const GRANULARITY = Set(["All", "Metadata", "AdminData", "ReferenceData", "AllClinicalData", "SingleSite", "SingleSubject"])
    const ARCHIVAL    = Set(["Yes", "No"])
    const ODMVERSION  = Set(["1.2", "1.2.1", "1.3", "1.3.1", "1.3.2"])

    const ODMNAMESPACE = Set([:Address
	:AdminData
	:Alias
	:Annotations
    :Annotation
	:ArchiveLayout
	:ArchiveLayoutRef
	:Association
	:AuditRecord
	:AuditRecords
	:BasicDefinitions
	:Certificate
	:CheckValue
	:City
	:ClinicalData
	:CodeList
	:CodeListItem
	:CodeListRef
	:Comment
	:ConditionDef
	:Country
	:CryptoBindingManifest
	:DateTimeStamp
	:Decode
	:Description
	:DisplayName
	:Email
	:EnumeratedItem
	:ErrorMessage
	:ExternalCodeList
	:ExternalQuestion
	:Fax
	:FirstName
	:Flag
	:FlagType
	:FlagValue
	:FormalExpression
	:FormData
	:FormDef
	:FormRef
	:FullName
	:GlobalVariables
	:ImputationMethod
	:Include
	:InvestigatorRef
	:ItemData
	:ItemDataAny
	:ItemDataBase64Binary
	:ItemDataBase64Float
	:ItemDataBoolean
	:ItemDataDate
	:ItemDataDatetime
	:ItemDataDouble
	:ItemDataDurationDatetime
	:ItemDataFloat
	:ItemDataHexBinary
	:ItemDataHexFloat
	:ItemDataIncompleteDate
	:ItemDataIncompleteDatetime
	:ItemDataIncompleteTime
	:ItemDataInteger
	:ItemDataIntervalDatetime
	:ItemDataPartialDate
	:ItemDataPartialDatetime
	:ItemDataPartialTime
	:ItemDataString
	:ItemDataTime
	:ItemDataURI
	:ItemDef
	:ItemGroupData
	:ItemGroupDef
	:ItemGroupRef
	:ItemRef
	:KeySet
	:LastName
	:LegalReason
	:Location
	:LocationRef
	:LoginName
	:Meaning
	:MeasurementUnit
	:MeasurementUnitRef
	:MetaDataVersion
	:MetaDataVersionRef
	:MethodDef
	:ODM
	:Organization
	:OtherText
	:Pager
	:Phone
	:Picture
	:PostalCode
	:Presentation
	:Protocol
	:ProtocolName
	:Question
	:RangeCheck
	:ReasonForChange
	:ReferenceData
	:Role
	:Signature
	:SignatureDef
	:SignatureRef
	:Signatures
	:SiteRef
	:SourceID
	:StateProv
	:StreetName
	:Study
	:StudyDescription
	:StudyEventData
	:StudyEventDef
	:StudyEventRef
	:StudyName
	:SubjectData
	:Symbol
	:TranslatedText
	:User
	:UserRef
	Symbol("ds:Signature")])

    const ODMATTRNAMESPACE = Set([:AddressAnnotationID
	:Archival
	:ArchiveLayoutOID
	:AsOfDateTime
	:AuditRecordID
	:Category
	:Code
	:CodedValue
	:CodeListOID
	:CollectionExceptionConditionOID
	:Comment
	:Comparator
	:Context
	:CreationDateTime
	:DataType
	:Description
	:Dictionary
	:Domain
	:EditPoint
	:EffectiveDate
	:FileOID
	:FileType
	:FormOID
	:FormRepeatKey
	:Granularity
	:ID
	:Id
	:ImageType
	:ImputationMethodOID
	:IsNull
	:IsReferenceData
	:ItemGroupOID
	:ItemGroupRepeatKey
	:ItemOID
	:KeySequence
	:Length
	:LocationOID
	:LocationType
	:Mandatory
	:MeasurementUnitOID
	:MetaDataVersionOID
	:MethodOID
	:Methodology
	:Name
	:ODMVersion
	:OID
	:OrderNumber
	:Origin
	:Originator
	:PdfFileName
	:PictureFileName
	:PresentationOID
	:PriorFileOID
	:Purpose
	:Rank
	:ref
	:Repeating
	:Role
	:RoleCodeListOID
	:SASDatasetName
	:SASFieldName
	:SASFormatName
	:SDSVarName
	:SeqNum
	:SignatureID
	:SignatureOID
	:SignificantDigits
	:SoftHard
	:SourceSystem
	:SourceSystemVersion
	:SponsorOrSite
	:StudyEventOID
	:StudyEventRepeatKey
	:StudyOID
	:SubjectKey
	:TransactionType
	:Type
	:UsedImputationMethod
	:UserOID
	:UserType
	:Value
	:Version
	Symbol("xml:lang")
	:xmlns])

    const ITEMDATATYPE = Set([:ItemDataAny
	:ItemDataBase64Binary
	:ItemDataBase64Float
	:ItemDataBoolean
	:ItemDataDate
	:ItemDataDatetime
	:ItemDataDouble
	:ItemDataDurationDatetime
	:ItemDataFloat
	:ItemDataHexBinary
	:ItemDataHexFloat
	:ItemDataIncompleteDate
	:ItemDataIncompleteDatetime
	:ItemDataIncompleteTime
	:ItemDataInteger
	:ItemDataIntervalDatetime
	:ItemDataPartialDate
	:ItemDataPartialDatetime
	:ItemDataPartialTime
	:ItemDataString
	:ItemDataTime
	:ItemDataURI])

    const CHNS = Dict(
    :Symbol => [:TranslatedText],
    )

    include("odmxml.jl")
    include("checknode.jl")
    include("spss.jl")
    include("ocl.jl")

end # module

ext: jl


text:
######################
# 1: The Julia type for ToricVarieties
######################
abstract type AbstractNormalToricVariety end

struct NormalToricVariety <: AbstractNormalToricVariety
           polymakeNTV::Polymake.BigObject
end
export NormalToricVariety

struct AffineNormalToricVariety <: AbstractNormalToricVariety
           polymakeNTV::Polymake.BigObject
end
export AffineNormalToricVariety

function pm_ntv(v::AbstractNormalToricVariety)
    return v.polymakeNTV
end

######################
# 2: Generic constructors
######################

@doc Markdown.doc"""
    AffineNormalToricVariety(C::Cone)

Construct the affine normal toric variety $U_{C}$ corresponding to a polyhedral
cone `C`.

# Examples
Set `C` to be the positive orthant in two dimensions.
```jldoctest
julia> C = Oscar.positive_hull([1 0; 0 1])
A polyhedral cone in ambient dimension 2

julia> antv = AffineNormalToricVariety(C)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function AffineNormalToricVariety(C::Cone)
    pmc = Oscar.pm_cone(C)
    fan = Polymake.fan.check_fan_objects(pmc)
    pmntv = Polymake.fulton.NormalToricVariety(fan)
    return AffineNormalToricVariety(pmntv)
end

@doc Markdown.doc"""
    NormalToricVariety(C::Cone)
Construct the (affine) normal toric variety $X_{\Sigma}$ corresponding to a
polyhedral fan $\Sigma = C$ consisting only of the cone `C`.
# Examples
Set `C` to be the positive orthant in two dimensions.
```jldoctest
julia> C = Oscar.positive_hull([1 0; 0 1])
A polyhedral cone in ambient dimension 2
julia> ntv = NormalToricVariety(C)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function NormalToricVariety(C::Cone)
    return AffineNormalToricVariety(C)
end

@doc Markdown.doc"""
    NormalToricVariety(PF::PolyhedralFan)

Construct the normal toric variety $X_{PF}$ corresponding to a polyhedral fan `PF`.

# Examples
Take `PF` to be the normal fan of the square.
```jldoctest
julia> square = Oscar.cube(2)
A polyhedron in ambient dimension 2

julia> nf = Oscar.normal_fan(square)
A polyhedral fan in ambient dimension 2

julia> ntv = NormalToricVariety(nf)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function NormalToricVariety(PF::PolyhedralFan)
    fan = Oscar.pm_fan(PF)
    pmntv = Polymake.fulton.NormalToricVariety(fan)
    if fan.N_MAXIMAL_CONES == 1
        return AffineNormalToricVariety( pmntv )
    end
    return NormalToricVariety(pmntv)
end

@doc Markdown.doc"""
    NormalToricVariety(P::Polyhedron)

Construct the normal toric variety $X_{\Sigma_P}$ corresponding to the normal
fan $\Sigma_P$ of the given polyhedron `P`.

Note that this only coincides with the projective variety associated to `P`, if
`P` is normal.

# Examples
Set `P` to be a square.
```jldoctest
julia> square = Oscar.cube(2)
A polyhedron in ambient dimension 2

julia> ntv = NormalToricVariety(square)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function NormalToricVariety(P::Polyhedron)
    fan = normal_fan(P)
    return NormalToricVariety(fan)
end

@doc Markdown.doc"""
    NormalToricVariety( r::Matrix{Int}, c::Vector{Vector{Int}} )

Construct the normal toric variety whose fan has ray generators `r` and maximal cones `c`.

# Examples
```jldoctest
julia> NormalToricVariety( [-1 5; 0 1; 1 0; 0 -1], [[1,2],[2,3],[3,4],[4,1]] )
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function NormalToricVariety( rays::Matrix{Int}, cones::Vector{Vector{Int}} )
    Incidence = Oscar.IncidenceMatrix(cones)
    arr = Polymake.@convert_to Array{Set{Int}} Polymake.common.rows(Incidence.pm_incidencematrix)
    pmntv = Polymake.fulton.NormalToricVariety(
        RAYS = Oscar.matrix_for_polymake(rays),
        MAXIMAL_CONES = arr,
    )
    if length( cones ) == 1
        return AffineNormalToricVariety( pmntv )
    end
    return NormalToricVariety( pmntv )
end

export NormalToricVariety

######################
# 3: Special constructors
######################

@doc Markdown.doc"""
    projective_space( d::Int )

Construct the projective space of dimension `d`.

# Examples
```jldoctest
julia> projective_space( 2 )
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function projective_space( d::Int )
    f = normal_fan(Oscar.simplex(d))
    pm_ntv = Polymake.fulton.NormalToricVariety(Oscar.pm_fan(f))
    return NormalToricVariety(pm_ntv)
end
export projective_space

@doc Markdown.doc"""
    hirzebruch_surface( r::Int )

Constructs the r-th Hirzebruch surface.

# Examples
```jldoctest
julia> hirzebruch_surface( 5 )
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function hirzebruch_surface( r::Int )
    Rays = [ 1 0; 0 1; -1 r; 0 -1]
    Cones = [[1,2],[2,3],[3,4],[4,1]]
    return NormalToricVariety( Rays, Cones )
end
export hirzebruch_surface

@doc Markdown.doc"""
    delPezzo( b::Int )

Constructs the delPezzo surface with b blowups for b at most 3.

# Examples
```jldoctest
julia> del_pezzo( 3 )
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2
```
"""
function del_pezzo( b::Int )
    if b < 0
        throw(ArgumentError("Number of blowups for construction of delPezzo surfaces must be non-negative."))
        return 0
    end
    if b == 0
        return projective_space( 2 )
    end
    if b == 1
        Rays = [ 1 0; 0 1; -1 0; -1 -1 ]
        Cones = [ [1,2],[2,3],[3,4],[4,1] ]
        return NormalToricVariety( Rays, Cones )
    end
    if b == 2
        Rays = [ 1 0; 0 1; -1 0; -1 -1; 0 -1 ]
        Cones = [ [1,2],[2,3],[3,4],[4,5],[5,1] ]
        return NormalToricVariety( Rays, Cones )
    end
    if b == 3
        Rays = [ 1 0; 1 1; 0 1; -1 0; -1 -1; 0 -1 ]
        Cones = [ [1,2],[2,3],[3,4],[4,5],[5,6],[6,1] ]
        return NormalToricVariety( Rays, Cones )
    end
    if b > 3
        throw(ArgumentError("delPezzo surfaces with more than 3 blowups are realized as subvarieties of toric ambient spaces. This is currently not supported."))
        return 0
    end
end
export del_pezzo

###############################################################################
###############################################################################
### Display
###############################################################################
###############################################################################
function Base.show(io::IO, ntv::AbstractNormalToricVariety)
    # fan = get_polyhedral_fan(ntv)
    pmntv = pm_ntv(ntv)
    ambdim = pmntv.FAN_AMBIENT_DIM
    print(io, "A normal toric variety corresponding to a polyhedral fan in ambient dimension $(ambdim)")
end


text:
## Implementation of the algorithm described in: "Invariant approximation of the minimal robust positively invariant set" S.V.Rakovic
# Invariant approximations of the minimal robust positively invariant set for and2D systems
#TODO 1D-systems

using LinearAlgebra
using LazySets
using Polyhedra
using LazySets.Approximations # overapproximate, SphericalDirections

""" Pontryagin Set Difference for HPolytope and HPolygon
    From Theory and Computation of Disturbance Invariant Sets for discrete time
    linear systems https://www.hindawi.com/journals/mpe/1998/934097/abs/"""
Base.:-(X::HPolytope, Y::HPolytope) = begin
    ax, bx, n = get_constraints(X)
    b_xminusy = bx .- ρ.(ax,Ref(Y))
    XminusY = HPolytope(vcat(ax'...), b_xminusy)
    return XminusY
end

Base.:-(X::AbstractPolytope, Y::AbstractPolytope) = begin
    ax, bx, n = get_constraints(X)
    b_xminusy = bx .- ρ.(ax,Ref(Y))
    XminusY = HPolytope(vcat(ax'...), b_xminusy)
    return XminusY
end

""" Pontryagin Set Difference for two Interval from IntervalArithmetic"""
Base.:-(a::Interval{Float64,LazySets.IntervalArithmetic.Interval{Float64}},
        b::Interval{Float64,LazySets.IntervalArithmetic.Interval{Float64}}) = begin
    low = a.dat.lo - b.dat.lo
    high = a.dat.hi - b.dat.hi
    return LazySets.Interval(low, high)
end

function get_support(P::AbstractPolytope)
    dims = LazySets.dim(P)
    x = Vector{Tuple{Float64,Float64}}(undef,dims)
    E = Matrix(I,dims,dims)
    for i=1:dims
        x[i] = (-ρ(-e(i,dims), P), ρ(e(i,dims), P))
    end
    return x
end

function get_constraints(𝒫::AbstractPolyhedron)
    # not all AbstractPolyhedron have 𝒫.constraints => constraints_list
    a = getfield.(constraints_list(𝒫),:a) # normal vector of halfspace
    b = getfield.(constraints_list(𝒫),:b) # support value of halfspace
    n = length(𝒫.constraints)
    return a,b, n
end

""" Equation (11) """
function αᴼ(s::Int64, A::AbstractMatrix, 𝒫::HPolytope)
    a, b, n_halfspaces = get_constraints(𝒫)
    hw_vec = zeros(n_halfspaces)
    for i=1:n_halfspaces
        hw_vec[i] = ρ((A^s)'*a[i], 𝒫) / b[i]
    end
    return maximum(hw_vec)
end

## Approximation to Invariant Set
"""  Equation (13) """
function M(s::Int64, A::AbstractMatrix, 𝒫::HPolytope)
    n = size(A, 2)
    E = Matrix{Float64}(I, n, n)
    hw_vec = zeros(n,2)
    for j=1:n
        eⱼ = E[:,j]
        sup_value_1 = 0
        sup_value_2 = 0
        for i=0:s-1
            sup_value_1 += ρ((A^i)'*eⱼ, 𝒫)
            sup_value_2 += ρ(-(A^i)'*eⱼ, 𝒫)
        end
        hw_vec[j,1] = sup_value_1
        hw_vec[j,2] = sup_value_2
    end
    return maximum(hw_vec)
end

# linear_map(A, P) is efficient, if A is invertible, since it does not have to convert
# to V-Representation!
function F(s::Int64, A::AbstractMatrix, 𝒫::HPolytope)
    # As defined in the paper, but different loop
    if s == 0
        return Ball2(zeros(size(A,1)), 0.0)
    else
        F = 𝒫
        for i=1:s-1
            F += linear_map(A^i,𝒫)
        end
        return F
    end
end

#TODO α and M for Interval type (1D)
# SinglePoint: A = Interval(0.5,0.5)
# or Apply just geometric series 1/(1-\rho)
function F(s::Int64, A::Real, 𝒫::Interval)
    Aint = Interval(A,A)
    # As defined in the paper, but different loop
    if s == 0
        return Interval(0.0,0.0)
    else
        F = 𝒫
        Ai = Aint
        for i=1:s-1
            F += 𝒫*Ai
            Ai = Ai*Aint
        end
        return F
    end
end

# The overapproximate case only works in the case of a HPolygon!
"""
    approx_mRPI(ϵ::Float64, A::AbstractMatrix, 𝒲::HPolyhedron;
                err_approx::Float64=1e-5, isLazy::Bool=false)

Calculate approximative mRPI set according to "Invariant Approximation of the mRPI set" by Rakovic et all."""
function approx_mRPI(ϵ::Float64, A::AbstractMatrix, 𝒲::HPolygon;
                    err_approx::Float64=1e-5, isLazy::Bool=false)
    s = 0; α = 1; m = 1
    while α > ϵ/(ϵ + m)
        s += 1
        α = αᴼ(s, A, 𝒲)
        m = M(s, A, 𝒲)
    end
    Fs = F(s, A, 𝒲)
    @show s, α
    lazyF_infty = (1 - α)*Fs # lazy (symoblic) approximative mRPI-set
    # Return the lazy (symbolic) set or an HPolygon
    (isLazy) && return lazyF_infty
    return overapproximate(lazyF_infty, err_approx) #overapproximate only works for 2D
end

# Only works for 3D 🙈
function approx_mRPI(ϵ::Float64, A::AbstractMatrix, 𝒲::HPolytope;
                    n_dirs::Int=10, isLazy::Bool=false)
    s = 0; α = 1; m = 1
    while α > ϵ/(ϵ + m)
        s += 1
        α = αᴼ(s, A, 𝒲)
        m = M(s, A, 𝒲)
    end
    Fs = F(s, A, 𝒲)
    @show s, α
    lazyF_infty = (1 - α)*Fs # lazy (symoblic) approximative mRPI-set
    # Return the lazy (symbolic) set or an HPolygon
    (isLazy) && return lazyF_infty
    dirs = SphericalDirections(n_dirs, n_dirs)
    return overapproximate(lazyF_infty, dirs)
end

function approx_mRPI(s::Int, A::Real, 𝒲::Interval)
    return F(s, A, 𝒲)
end

## use https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/master/examples/Minimal%20Robust%20Positively%20Invariant%20Set.ipynb for 3D and more!

# W = polyhedron(Wv, CDDLib.Library())

# using Polyhedra
# Wv = vrep([[x, y] for x in [-1.0, 1.0] for y in [-1.0, 1.0]])
# W = polyhedron(Wv)
#
# function vertices_to_polhedron(x_vec,y_vec)
#     Wv = vrep([[x, y] for x in x_vec for y in y_vec])
#     W = polyhedron(Wv)
# end
#
# function vertices_to_polhedron(x_vec,y_vec, z_vec)
#     Wv = vrep([[x, y,z] for x in x_vec for y in y_vec for z in z_vec])
#     W = polyhedron(Wv)
# end
#
# function Fs(s::Integer, verbose=1)
#     @assert s ≥ 1
#     F = W
#     A_W = W
#     for i in 1:(s-1)
#         A_W = A * A_W
#         F += A_W
#         if verbose ≥ 1
#             println("Number of points after adding A^$i * W: ", npoints(F))
#         end
#         removevredundancy!(F)
#         if verbose ≥ 1
#             println("Number of points after removing redundant ones: ", npoints(F))
#         end
#     end
#     return F
# end


text:
using JsonGrinder, JSON, Test, SparseArrays, Mill, HierarchicalUtils
import HierarchicalUtils: printtree
using Mill: nobs

@testset "ExtractMultipleRepresentation" begin
	@testset "with uniontypes" begin
		ex = MultipleRepresentation((
			ExtractCategorical(["Olda", "Tonda", "Milda"], true),
			ExtractString(true)))
		e = ex("Olda")

		@test nobs(e) == 1
		@test length(e.data) == 2
		@test size(e.data[1].data) == (4, 1)
		@test e.data[1].data == MaybeHotMatrix([2],4)
		@test size(e.data[2].data) == (2053, 1)
		@test findall(x->x > 0, SparseMatrixCSC(e.data[2].data)) .|> Tuple == [(206, 1), (272, 1), (624, 1), (738, 1), (1536, 1), (1676, 1)]

		e = ex(extractempty)
		@test nobs(e) == 0
		@test nobs(e[:e1]) == 0
		@test e[:e1].data isa MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}
		@test nobs(e[:e2]) == 0
		@test e[:e2].data isa NGramMatrix{Union{Missing, String},Vector{Union{Missing, String}},Union{Missing, Int64}}

		ex2 = MultipleRepresentation((
			ExtractCategorical(["Olda", "Tonda", "Milda"], true),
			JsonGrinder.ExtractString(true)))
		@test hash(ex) === hash(ex2)
		@test ex == ex2

		ex3 = MultipleRepresentation((
			ExtractCategorical(["Polda", "Tonda", "Milada"], true),
			ExtractString(true)))
		@test hash(ex) !== hash(ex3)
		@test ex != ex3

		ex = MultipleRepresentation((
			ExtractArray(JsonGrinder.extractscalar(Float32, true)),
			ExtractString(true)))
		e = ex("Olda")
		@test e[:e1].data.data == fill(1,1,0)
		@test e[:e2].data.s == ["Olda"]

		e = ex([1,2])
		@test e[:e1].data.data == [1 2]
		@test e[:e2].data.s ≃ [missing]

		ex = MultipleRepresentation([
        	ExtractCategorical(["packer", "plain"], true),
        	ExtractArray(ExtractString(true))
        ])

		@test ex(["sfx", "plain"]) ≃ ProductNode((;
			:e1=>ArrayNode(MaybeHotMatrix([missing], 3)),
			:e2=>BagNode(
				ArrayNode(NGramMatrix(["sfx", "plain"])),
				[1,1])
		))

		@test ex(["sfx", "plain"]) ≃ ex(Any["sfx", "plain"])
		@test ex(["sfx", "plain"]) ≃ ex(String["sfx", "plain"])
	end

	@testset "without uniontypes" begin
		ex = MultipleRepresentation((
			ExtractCategorical(["Olda", "Tonda", "Milda"], false),
			ExtractString(false)))
		e = ex("Olda")

		@test length(e.data) == 2
		@test size(e.data[1].data) == (4, 1)
		@test e.data[1].data == MaybeHotMatrix([2],4)
		@test size(e.data[2].data) == (2053, 1)
		@test findall(x->x > 0, SparseMatrixCSC(e.data[2].data)) .|> Tuple == [(206, 1), (272, 1), (624, 1), (738, 1), (1536, 1), (1676, 1)]

		e = ex(extractempty)
		@test nobs(e) == 0
		@test nobs(e[:e1]) == 0
		@test e[:e1].data isa OneHotMatrix{UInt32, UInt32(4), Vector{UInt32}}
		@test nobs(e[:e2]) == 0
		@test e[:e2].data isa NGramMatrix{String,Array{String,1},Int64}

		ex2 = MultipleRepresentation((
			ExtractCategorical(["Olda", "Tonda", "Milda"], false),
			ExtractString(false)))
		@test hash(ex) === hash(ex2)
		@test ex == ex2

		ex3 = MultipleRepresentation((
			ExtractCategorical(["Polda", "Tonda", "Milada"], false),
			ExtractString(false)))
		@test hash(ex) !== hash(ex3)
		@test ex != ex3

		ex = MultipleRepresentation([
        	ExtractCategorical(["packer", "plain"], false),
        	ExtractArray(ExtractString(false))
        ])
		@test_throws ErrorException ex(["sfx", "plain"])
	end
end

@testset "show" begin
	ex = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"], true),
		JsonGrinder.ExtractString(true)))
	e = ex("Olda")

	@test buf_printtree(ex, trav=true) ==
	"""
	MultiRepresentation [""]
	  ├── e1: Categorical d = 4 ["E"]
	  └── e2: String ["U"]
	"""

	@test buf_printtree(e, trav=true) ==
	"""
	ProductNode with 1 obs [""]
	  ├── e1: ArrayNode(4×1 MaybeHotMatrix with Union{Missing, Bool} elements) with 1 obs ["E"]
	  └── e2: ArrayNode(2053×1 NGramMatrix with Union{Missing, Int64} elements) with 1 obs ["U"]
	"""

	ex = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"], false),
  		JsonGrinder.ExtractString(false)))
  	e = ex("Olda")

  	@test buf_printtree(ex, trav=true) ==
	"""
	MultiRepresentation [""]
	  ├── e1: Categorical d = 4 ["E"]
	  └── e2: String ["U"]
	"""

	@test buf_printtree(e, trav=true) ==
	"""
	ProductNode with 1 obs [""]
	  ├── e1: ArrayNode(4×1 OneHotArray with Bool elements) with 1 obs ["E"]
	  └── e2: ArrayNode(2053×1 NGramMatrix with Int64 elements) with 1 obs ["U"]
	"""
end


text:
module SnakeAI
#=
TODO:
- Plotting in real time of training's progress
- Modifying the state's vector to include information about the snake's body
- Adding tests on memory types
=#

export play_snake,
    demo_agent,
    train_agent,
    Game,
    play_step!,
    reset!

include("constants.jl")
include("point.jl")
include("snake.jl")
include("game.jl")

# Model module
include("model/Model.jl")
using .Model
export save_model

# Agent module
include("agent/Agent.jl")
using .Agent
export SnakeAgent,
    CircularBufferMemory,
    train!

include("functions.jl")

end


text:
# resp
printstyled("  Type preservation after processing\n", color=:light_green)

for f in String["convert_seis!", "demean!", "detrend!", "filtfilt!", "merge!", "sync!", "taper!", "ungap!", "unscale!"]
  printstyled(string("    ", f, "\n"), color=:light_green)
  S = randSeisData(s=1.0)
  if f == "filtfilt!"
    while true
      i = findall(S.fs .< 30.0)
      deleteat!(S, i)
      if S.n > 0
        break
      end
      S = randSeisData(s=1.0)
    end
  end
  T = [eltype(S.x[i]) for i=1:S.n]
  id = S.id
  getfield(SeisIO, Symbol(f))(S)
  for i = 1:S.n
    j = findfirst(id.==S.id[i])
    @test T[j] == eltype(S.x[i])
  end
end

printstyled(string("    remove_resp!\n"), color=:light_green)
r = fctoresp(1.0, 1.0/sqrt(2.0))
r2 = fctoresp(0.0333, 1.0/sqrt(2.0))
S = randSeisData(3, s=1.0)
S.resp[1] = r
S.resp[2] = deepcopy(r)
S.resp[3] = r2
S.x[1] = randn(Float32, S.t[1][end,1])
T = [eltype(S.x[i]) for i=1:S.n]
remove_resp!(S)
for i=1:S.n
  @test T[i] == eltype(S.x[i])
end


text:
using Mustache

template = """
{{#jobs}}
julia --project experiments/execute.jl  {{{arg}}}
{{/jobs}}
"""

function generate_jobs(directory, data_dir, no_samples, build_no, exp_name)
    items = Dict()
    jobs = Vector{Dict}()
    for ix in 1:no_samples
        args = """ --build "$(data_dir)" "/tmp/scratch/" "/home/sdalvi/work/Surrogate/sample_configuration_$(exp_name)_$(build_no)/" $ix $build_no $exp_name """
        push!(jobs, Dict("arg" => args))
    end
    items["jobs"] = jobs
    filename = joinpath(directory, "jobs.txt")
    open(filename, "w") do io
        write(io, Mustache.render(template, items))
    end
    println("Wrote $filename")
    return
end


text:
"""
Defines a group which behaves as an elementwise container, and
where the group dimensions are the same as those of field X.

For instance, if:
```
struct Neurons <: Group{:v}
    v::Vector{Float32}
end
```
then `size(group::Neurons) = size(group.v)`.

Groups also support the broadcasting dot-syntax: given a function
`f(::Elem{Neurons})`, calling `f.(group)` is equivalent to calling
f on each Elem in the group.
"""
abstract type Group{X} end

@generated function Base.length(g::Group{X}) where {X}
    :(length(g.$X))
end

@generated function Base.size(g::Group{X}, args...) where {X}
    :(size(g.$X, args...))
end

@generated function Base.eachindex(g::Group{X}) where {X}
    :(eachindex(g.$X))
end

@generated function Base.getindex(g::Group{X},i...) where {X}
    :(Elem(g,i))
end

function Base.broadcast(f, g::Group)
    @inbounds @simd for i in eachindex(g)
        f(Elem(g,i))
    end
    nothing
end

function Base.broadcast(f, g::Group, g2::Group)
    @assert size(g) == size(g2)
    @inbounds @simd for i in eachindex(g)
        f(Elem(g,i), Elem(g2,i))
    end
    nothing
end


text:
using Test
using LinearAlgebra, SparseArrays
using ArnoldiMethod: partialschur, partialeigen
using Random

@testset "Schur to eigen $T take $i" for T in (Float64,ComplexF64), i in 1:10
    Random.seed!(i)
    A = spdiagm(0 => 1:100) + sprand(100, 100, 0.01)
    ε = 1e-7
    minim, maxim = 10, 20

    decomp, history = partialschur(A, nev=minim, tol=ε, restarts=200)
    @test history.converged

    vals, vecs = partialeigen(decomp)

    # This test seems a bit flaky sometimes -- have to dig into it.
    for i = 1 : minim
        @test norm(A * vecs[:,i] - vecs[:,i] * vals[i]) < ε * abs(vals[i])
    end
end


text:
# To run: `julia bfs_test.jl`
include("bfs.jl")

nodes = [
     1 => [2, 3, 4],
     2 => [1, 5, 6],
     3 => [1],
     4 => [1, 7, 8],
     5 => [2, 9, 10],
     6 => [2],
     7 => [4, 11, 12],
     8 => [4],
     9 => [5],
    10 => [5],
    11 => [7],
    12 => [7]
]

visited = Number[]
Bfs.visit(1, nodes, n -> push!(visited, n))

print_joined(STDOUT, visited, ", ")
println()

max_stars_repo_name: adlawson/search-algorithms

max_issues_repo_name: adlawson/search-algorithms

max_forks_repo_name: adlawson/search-algorithms


text:
# Methods for checking matrix sizes etc.
import LinearAlgebra: issymmetric, eigmin

issquare(x::Matrix) = size(x, 1) == size(x, 2)
ispossemidef(x::Matrix) = issymmetric(x) && (eigmin(x) >= 0)

function confirm_matrix_sizes(F, B, V, G, W)

    nx = size(F, 1)
    nu = size(B, 2)
    ny = size(G, 1)

    @assert size(F) == (nx, nx)
    @assert size(B) == (nx, nu)
    @assert size(V) == (nx, nx)

    @assert size(G) == (ny, nx)
    @assert size(W) == (ny, ny)

    return nx, ny, nu

end


text:
// multiple variables of the same type declared
// and possibly initialized in the same statement

int main() {
 int x, y = 7;
 x = -1234234;
 printInt(x);
 printInt(y);
 return 0 ;

}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
# # Network Matrices

# `PowerSystems.jl` is able to build classic power systems modeling network matrices such as
# [Ybus](https://en.wikipedia.org/wiki/Nodal_admittance_matrix), [PTDF](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Power_Transfer_Distribution_Factors.htm) and LODF.

# Check section [Network Matrices](@ref net_mat) for more details.

# ## Overview
#
# Network matrices are implemented in `PowerSystems.jl` as arrays that enable using Branch or Bus
# names as indexes to facilitate exploration and analysis. Ybus is stored as an SparseMatrix
# and the PTDF and LODF are stored as dense matrices. **Note*** Ybus is converted to a dense
# matrix for printing in the REPL,
# The network matrices code implements the Goderya algorithm to find islands.

using PowerSystems
const PSY = PowerSystems
DATA_DIR = "../../../data" #hide
system_data = System(joinpath(DATA_DIR, "matpower/case14.m"))

# ## Ybus
#
# The ybus can be calculated as follows:
ybus = Ybus(system_data)

# The matrix can be indexed using directly the bus numbers. In this example buses are numbered
# 1-14. However, in large systems buses don't usually follow sequential numbering. You can access
# the entries of the Ybus with direct indexing or using the buses.

ybus_entry = ybus[3,3]
#
bus3 = get_component(Bus, system_data, "Bus 3     HV")
ybus_entry = ybus[bus3, bus3]

# We recognize that many models require matrix operations. For those cases, you can access the
# sparse data as follows:

sparse_array = get_data(ybus)

# ## PTDF
#
# The PTDF matrix can be calculated as follows:
ptdf = PTDF(system_data)

# The entries of the matrix can also be indexed using the devices directly:
line = get_component(Line, system_data, "Bus 1     HV-Bus 2     HV-i_1")
bus = get_component(Bus, system_data, "Bus 3     HV")

ptdf_entry = ptdf[line, bus]

# Alternatively, the branch name and bus number can be used. For instance, for the same PTDF
# entry we can do:

ptdf_entry = ptdf["Bus 1     HV-Bus 2     HV-i_1", 3]

# PTDF also takes a vector of distributed slacks, for now this feature requires passing a
# vector of weights with the same number of elements are buses in the system. For more details
# check the API entry for [`PTDF`](@ref).

# In the same fashion as in the case of the Ybus, we recognize many models use matrix operations
# when using the PTDF matrix. The raw data is available using `get_data` and the lookup table
# between (names,numbers) and (i,j) matrix entries is available using [`get_lookup`](@ref)

branch_lookup, bus_lookup = get_lookup(ptdf);

# Branch lookup:

branch_lookup

# Bus Look up

bus_lookup

# ## LODF

# LODF and PTDF share the same characteristics in terms of indexing and calculation:

lodf = LODF(system_data)


text:
struct CyclicContainer <: AbstractVector{String}
    c::Vector{String}
    n::Int
end
CyclicContainer(c) = CyclicContainer(c, 0)

Base.length(c::CyclicContainer) = typemax(Int)
Base.size(c::CyclicContainer) = size(c.c)
Base.iterate(c::CyclicContainer, state=1) = Base.iterate(c.c, state)
Base.getindex(c::CyclicContainer, i) = c.c[(i-1)%length(c.c) + 1]
function Base.getindex(c::CyclicContainer)
    c.n += 1
    c[c.n]
end
Base.iterate(c::CyclicContainer, i = 1) = iterate(c.c, i)

COLORSCHEME = [
    "#6D44D0",
    "#2CB3BF",
    "#1B1B1B",
    "#DA5210",
    "#03502A",
    "#866373",
]

COLORS = CyclicContainer(COLORSCHEME)
LINESTYLES = CyclicContainer(["-", ":", "--", "-."])

using PyPlot
using3D()

PyPlot.rc("lines", lw = 2)
PyPlot.rc("axes", grid = true)
PyPlot.rc("grid", color = "0.75", alpha = 0.75)

PyPlot.rc("font", size = 22) # set default fontsize
PyPlot.rc("xtick", labelsize = 20)
PyPlot.rc("ytick", labelsize = 20)
PyPlot.rc("axes", labelsize = 24)
PyPlot.rc("legend", fontsize = 20)
# PyPlot.rc("font", family = "Times New Roman") # Serif main font
PyPlot.rc("font", family = "DejaVu Sans") # sans main font
# PyPlot.rc("mathtext", rm = "sanserif", fontset="dejavusans") # sans math font
PyPlot.rc("mathtext", rm = "serif", fontset="cm") # serif math font

for z in ("x", "y")
    PyPlot.rc("$(z)tick.major", size = 6, width = 1.2)
    PyPlot.rc("$(z)tick.minor", size = 3, visible = false)
end

figx = 9
figy = 7
PyPlot.rc("figure", figsize = (figx, figy))
PyPlot.rc("savefig", dpi = 300, transparent = false, format = "png")

# set default color cycle
PyPlot.rc("axes", prop_cycle = matplotlib.cycler(color=COLORS.c))


text:
struct Config
    action::String
    parameters::Dict
    storage::Dict
    image_parameters::Dict
    authorization::Dict
end

max_stars_repo_name: keboola/julia-component

max_issues_repo_name: keboola/julia-component

max_forks_repo_name: keboola/julia-component


text:
# This is a temporary module to validate `AbstractImageFilter` idea
# proposed in https://github.com/JuliaImages/ImagesAPI.jl/pull/3
module ThresholdAPI

using ImageBase
using MappedArrays

"""
    AbstractImageAlgorithm

The root of image algorithms type system
"""
abstract type AbstractImageAlgorithm end

"""
    AbstractImageFilter <: AbstractImageAlgorithm

Filters are image algorithms whose input and output are both images
"""
abstract type AbstractImageFilter <: AbstractImageAlgorithm end

include("find_threshold.jl")
include("build_histogram.jl")

export build_histogram

end  # module ThresholdAPI

max_stars_repo_name: zygmuntszpak/HistogramThresholding

max_forks_repo_name: zygmuntszpak/HistogramThresholding


text:
"""
    modularity(g, c)

Return a value representing Newman's modularity `Q` for the undirected graph
`g` given the partitioning vector `c`.
"""
function modularity end
@traitfn function modularity(g::::(!IsDirected), c::Vector)
    m = 2*ne(g)
    m == 0 && return 0.
    nc = maximum(c)
    a = zeros(Int,nc)
    Q = 0
    for u in vertices(g)
        for v in neighbors(g,u)
            if u <= v
                c1 = c[u]
                c2 = c[v]
                if c1 == c2
                    Q += 2
                end
                a[c1] += 1
                a[c2] += 1
            end
        end
    end
    Q = Q*m
    @inbounds for i=1:nc
        Q -= a[i]*a[i]
    end
    return Q/m/m
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
module ExampleWebSocket

using Bukdu
const ServerHost = "localhost"
const ServerPort = 8080

struct WSController <: ApplicationController
    conn::Conn
end

function index(c::WSController)
    render(HTML, """
<!DOCTYPE html>

  <meta charset="utf-8" />
  <title>WebSocket Test</title>

WebSocket Test

julia> ws = first(Bukdu.websockets())
HTTP.WebSockets.WebSocket{Sockets.TCPSocket}(Sockets.TCPSocket(RawFD(0x0000001a) active, 0 bytes waiting), 0x01, true, UInt8[0x57, 0x65, 0x62, 0x53, 0x6f, 0x63, 0x6b, 0x65, 0x74, 0x20, 0x72, 0x6f, 0x63, 0x6b, 0x73], UInt8[], false, false)

julia> write(ws, "hello")
5

""")
end

end # module ExampleWebSocket

import PushInterface: RTMIDI
import .RTMIDI: rtmidi_in_create_default, rtmidi_open_port, rtmidi_in_set_callback
import .RTMIDI: destroy, EventCB, rtmidi_callback_func

device = rtmidi_in_create_default()
rtmidi_open_port(device, 1, "Ableton Push 2 User Port")
if !(device[].ok)
    @error :rtmidi_open_port "could not found ableton push 2"
    exit(1)
end

cb_ptr = @cfunction rtmidi_callback_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{EventCB})
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
finalizer(destroy, device)

using Bukdu
import .ExampleWebSocket: WSController, ServerHost, ServerPort, index

plug(Plug.Logger, access_log=(path=normpath(@__DIR__, "access.log"),), formatter=Plug.LoggerFormatter.datetime_message)

routes() do
    get("/", WSController, index)
    plug(Plug.WebSocket) #
end

Bukdu.start(ServerPort, host=ServerHost)

function midi_callback(msg)
    if 0x90 == first(msg)
        websocks = Bukdu.websockets()
        if isempty(websocks)
            @info "visit the url"
        else
            ws = first(websocks)
            write(ws, repr(msg))
        end
    end
end

function _callback_async_loop(cond, r_ecb)
    while isopen(cond)
        wait(cond)
        msg = codeunits(unsafe_string(r_ecb[].message))
        midi_callback(msg)
    end
end

cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)

# Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule IGLWrap_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("IGLWrap")
JLLWrappers.@generate_main_file("IGLWrap", UUID("283677c1-8365-580c-84e5-ef4b5d190868"))
end  # module IGLWrap_jll


text:
# Copyright (c) 2021 California Institute of Technology (“Caltech”).
# U.S. Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# - Neither the name of Caltech nor its operating division, the Jet Propulsion
#   Laboratory, nor the names of its contributors may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# Here we collect the angular cryo-window functions.

function sphericalharmonicsy(l, m, θ, ϕ)
    #return scipy.special.sph_harm(m, l, ϕ, θ)  # scipy has ℓ,m and θ,ϕ reversed
    rY1 = sphevaluate(θ, ϕ, l, abs(m))
    rY2 = sphevaluate(θ, ϕ, l, -abs(m))
    if m < 0
        return (rY1 - im*rY2) / √2
    elseif m == 0
        return rY1 + 0*im
    else
        return (-1)^m * (rY1 + im*rY2) / √2
    end
end

function realsphericalharmonicsy(l, m, θ, ϕ)
    #Ylm = sphericalharmonicsy(l, m, θ, ϕ)
    #if m < 0
    #    return √2 * (-1)^m * imag(Ylm)
    #elseif m == 0
    #    return real(Ylm)
    #else
    #    return √2 * (-1)^m * real(Ylm)
    #end
    return sphevaluate.(θ, ϕ, l, m)
end

function get_wl_Y_j_lm(bmask, nside_hires=npix2nside(length(bmask)); lmax=Inf)
    npix = length(bmask)
    lmsize = npix
    pix = collect((1:npix)[bmask])
    nside = npix2nside(npix)
    pixinv = fill(0, npix)
    for i=1:length(pix)
        pixinv[pix[i]] = i
    end

    pw = Healpix.pixwin(nside)
    pixwin = Spline1D(0:length(pw)-1, pw)

    for i in pix
        θ, ϕ = Healpix.pix2angRing(nside, [i])
        n = Healpix.ang2pixRing(Resolution(nside), θ[1], ϕ[1])
        @assert n==i
    end

    bmask_hires = udgrade(HealpixMap{Bool,Healpix.RingOrder}(bmask), nside_hires)
    npix_hires = nside2npix(nside_hires)
    pix_hires = collect((1:npix_hires)[bmask_hires])
    θ, ϕ = Healpix.pix2angRing(nside_hires, pix_hires)
    pixmap = [Healpix.ang2pixRing(Resolution(nside), θ[i], ϕ[i]) for i=1:length(pix_hires)]

    wl_Yjlm = fill(0.0, length(pix), lmsize)
    l = 0
    m = 0
    for lm=1:lmsize
        (l > lmax) && continue
        Ylm = CryoFaBs.realsphericalharmonicsy(l, m, θ, ϕ)
        for i=1:length(pix_hires)
            n = pixinv[pixmap[i]]
            n==0 && @show i,pixmap[i],n
            wl_Yjlm[n,lm] += Ylm[i]
            if !isfinite(wl_Yjlm[n,lm])
                @error "Non-finite value in w_Yjlm! Setting to zero." l m n lm wl_Yjlm[n,lm] i Ylm[i] pixmap[i] θ[i] ϕ[i] CryoFaBs.sphericalharmonicsy(l, m, θ[i], ϕ[i])
                wl_Yjlm[n,lm] = 0
            end
        end
        wl_Yjlm[:,lm] .*= npix / npix_hires

        if m < 0
            m = -m
        else # m >= 0
            m = -m - 1
        end

        if m < -l
            l += 1
            m = 0
        end
    end
    return wl_Yjlm
end

function fill_l(n)
    ℓ = Int[]
    l = 0
    m = 0
    for i=1:n
        push!(ℓ, l)
        m += 1
        if m > l
            l += 1
            m = -l
        end
    end
    return ℓ
end

function angular_cryo_window_sum_m(w̃)
    ell = fill_l(size(w̃,2))
    ell_uniq = unique(ell)
    ww = fill(0.0, size(w̃,1), length(ell_uniq))
    for j=1:length(ell_uniq)
        m = (ell .== ell_uniq[j])
        for i=1:size(w̃,1)
            ww[i,j] = sum(abs2.(w̃[i,m]))
        end
    end
    return ell_uniq, ww
end

function angular_cryo_window(cfb::AngularCryoFaB; do_calc=true, lmax=Inf)
    nside = cfb.nside
    npix = nside2npix(nside)
    pwspl = Spline1D(0:4*nside, Healpix.pixwin(nside))  # Don't really need this to be a spline.
    pwlm = pwspl.(fill_l(npix))
    bmask = @. cfb.mask > 0
    if do_calc
        wjlm = pwlm' .* get_wl_Y_j_lm(bmask, nside; lmax=lmax)
        w̃ = cfb.TransInv * wjlm
    else
        w̃ = zeros(size(cfb,2), length(pwlm))
    end
    unique_ell, ww = angular_cryo_window_sum_m(w̃)
    return unique_ell, ww
end

function angular_cryo_window(cfb::AngRadCryoFaB; do_calc=true, lmax=Inf)
    return angular_cryo_window(cfb.angfab; do_calc=do_calc, lmax=lmax)
end

# vim: set sw=4 et sts=4 :


text:
using SparseArrays
using CUDA
using Random
using Adapt

function D2_cpu(idata,odata,Nx,Ny,h)
    for i in 1:Nx
        for j in 1:Ny
            if 2 <= i <= Nx-1 && 2 <= j <= Ny - 1
                odata[i,j] = (idata[i-1,j] + idata[i+1,j] + idata[i,j-1] + idata[i,j+1] - 4*idata[i,j])
            end
        end
    end
end

function D2_split_dev(idata,odata,Nx,Ny,h,::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = (i-1)*Ny+j

    if 0 <= i <= Nx && 1 <= j <= Ny
        odata[i,j] = 0
    end

    if 2 <= i <= Nx-1 && 2 <= j <= Ny - 1
        odata[i,j] = (idata[i-1,j] + idata[i+1,j] + idata[i,j-1] + idata[i,j+1] - 4*idata[i,j])
    end

    nothing
end

function matrix_free_cpu(idata,odata,Nx,Ny,h)
    # odata .= 0
    # for i in 1:Nx
    #     for j in 1:Ny
    #         if (i == 1) && (j == 1)
    #             odata[i,j] = (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+1]) / 4
    #         end
    #         if (i == 1) && (j == Ny)
    #             odata[i,j] = (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-1]) / 4
    #         end
    #         if (i == Nx) && (j == 1)
    #             odata[i,j] = (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+1]) / 4
    #         end
    #         # if (i == Nx) && (j == Ny)
    #         #     odata[i,j] = (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-1]) / 4
    #         # end
    #     end
    # end
    odata .= 0
    # alpha1 = alpha2 = alpha3 = alpha4 = beta = 1
    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1
    (i,j) = (1,1)

    odata[i,j] += (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+2]) / 4 # D2

    odata[i,j] += 2 * alpha3 * (( 1.5* idata[i,j] - 2*idata[i+1,j] + 0.5*idata[i+2,j])) / 4 # Neumann

    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha1 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j+1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j+2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (1,Ny)
    odata[i,j] += (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-2]) / 4 # D2

    odata[i,j] += 2 * alpha3 * (1.5 * idata[i,j] - 2*idata[i+1,j] + 0.5 * idata[i+2,j]) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha2 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j-1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j-2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,1)
    odata[i,j] += (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+2]) / 4 # D2

    odata[i,j] += 2 * alpha4 * (( 1.5* idata[i,j] - 2*idata[i-1,j] + 0.5*idata[i-2,j])) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha1 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j+1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j+2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,Ny)
    odata[i,j] += (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-2]) / 4 # D2

    odata[i,j] += 2 * alpha4 * (1.5 * idata[i,j] - 2*idata[i-1,j] + 0.5 * idata[i-2,j]) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha2 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j-1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j-2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (1,2:Ny-1)
    # odata[i,j] .+= ((idata[i,j] .- 2*idata[i+1,j] .+ idata[i+2,j] .+ idata[i,j .- 1] .- 2*idata[i,j] .+ idata[i,j .+ 1] ) ./ 2 # D2
    #                 # .+ 2 .* alpha3 .* (1.5 * idata[i,j] .- 2*idata[i+1,j] .+ 0.5 * idata[i+2,j]) ./ 2
    #                 )
    # odata[i,j] .+= 2 * alpha3 * (1.5 * idata[i,j] .- 2*idata[i+1,j] .+ 0.5 * idata[i+2,j]) ./ 2 # Neumann

    # test @view
    odata[i,j] .+= (view(idata,i,j) .- 2*view(idata,i+1,j) .+ view(idata,i+2,j) .+ view(idata,i,j .-1) .- 2*view(idata,i,j) .+ view(idata,i,j .+ 1)
                    .+ 2 * alpha3 .* (1.5 * view(idata,i,j) .- 2*view(idata,i+1,j) .+ 0.5 * view(idata,i+2,j))) ./ 2

    (i,j) = (Nx,2:Ny-1)
    # odata[i,j] .+= (idata[i,j] .- 2*idata[i-1,j] .+ idata[i-2,j] .+ idata[i,j .- 1] .- 2*idata[i,j] .+ idata[i,j .+ 1]) / 2 # D2
    # odata[i,j] .+= 2 * alpha4 * (1.5 * idata[i,j] .- 2*idata[i-1,j] .+ 0.5 * idata[i-2,j]) ./ 2 # Neumann
    odata[i,j] .+= ((view(idata,i,j) .- 2*view(idata,i-1,j) .+ view(idata,i-2,j) .+ view(idata,i,j.-1) .- 2*view(idata,i,j) .+ view(idata,i,j.+1))
                    .+ 2 * alpha4 * (1.5 * view(idata,i,j) .- 2*view(idata,i-1,j) .+ 0.5 * view(idata,i-2,j))) ./ 2

    (i,j) = (2:Nx-1,1)
    # odata[i,j] .+= (idata[i.-1,j] .- 2*idata[i,j] .+ idata[i.+1,j] .+ idata[i,j] .- 2*idata[i,j+1] .+ idata[i,j+2]) / 2 # D2
    odata[i,j] .+= (view(idata,i .- 1,j) .- 2*view(idata,i,j) .+ view(idata,i .+ 1,j) .+ view(idata,i,j) .- 2*view(idata,i,j+1) .+ view(idata,i,j+2)) ./ 2

    # odata[i,j] .+= (2 * beta * (1.5 * idata[i,j]) .+ 2 * alpha1 * (idata[i,j]) * h) ./ 2 # Dirichlet
    # odata[i,j+1] .+= (2 * beta * (-1 * idata[i,j])) # Dirichlet
    # odata[i,j+2] .+= (0.5 * beta * (idata[i,j]))  # Dirichlet
    odata[i,j] .+= (2 * beta * (1.5 * view(idata,i,j)) .+ 2 * alpha2 * view(idata,i,j) * h) ./ 2
    odata[i,j+1] .+= (2 * beta * (-1 * view(idata,i,j)))
    odata[i,j+2] .+= (0.5 * beta * (view(idata,i,j)))

    (i,j) = (2:Nx-1,Ny)
    # odata[i,j] .+= (idata[i.-1,j] .- 2*idata[i,j] .+ idata[i.+1,j] .+ idata[i,j] .- 2*idata[i,j-1] .+ idata[i,j - 2]) / 2 # D2
    odata[i,j] .+= (view(idata,i .- 1, j) .- 2*view(idata,i,j) .+ view(idata,i .+1,j) .+ view(idata,i,j) .- 2*view(idata,i,j-1) .+ view(idata,i,j-2)) ./ 2

    # odata[i,j] .+= (2 * beta * (1.5 * idata[i,j]) .+ 2 * alpha1 * (idata[i,j]) * h) ./ 2 # Dirichlet
    # # odata[i,j-1] .+= (2 * beta * (-1 * idata[i,j])) # Dirichlet
    # # odata[i,j-2] .+= (0.5 * beta * (idata[i,j])) # Dirichlet
    odata[i,j] .+= (2 * beta * (1.5 * view(idata,i,j)) .+ 2 * alpha1 * view(idata,i,j) * h) ./ 2
    odata[i,j-1] .+= (2 * beta * (-1 * view(idata,i,j)))
    odata[i,j-2] .+= (0.5 * beta * (view(idata,i,j)))
end

function matrix_free_cpu_v2(idata,odata,Nx,Ny,h)
    # odata .= 0
    # for i in 1:Nx
    #     for j in 1:Ny
    #         if (i == 1) && (j == 1)
    #             odata[i,j] = (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+1]) / 4
    #         end
    #         if (i == 1) && (j == Ny)
    #             odata[i,j] = (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-1]) / 4
    #         end
    #         if (i == Nx) && (j == 1)
    #             odata[i,j] = (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+1]) / 4
    #         end
    #         # if (i == Nx) && (j == Ny)
    #         #     odata[i,j] = (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-1]) / 4
    #         # end
    #     end
    # end
    odata .= 0
    # alpha1 = alpha2 = alpha3 = alpha4 = beta = 1
    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1
    (i,j) = (1,1)

    odata[i,j] += (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+2]) / 4 # D2

    odata[i,j] += 2 * alpha3 * (( 1.5* idata[i,j] - 2*idata[i+1,j] + 0.5*idata[i+2,j])) / 4 # Neumann

    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha1 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j+1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j+2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (1,Ny)
    odata[i,j] += (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-2]) / 4 # D2

    odata[i,j] += 2 * alpha3 * (1.5 * idata[i,j] - 2*idata[i+1,j] + 0.5 * idata[i+2,j]) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha2 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j-1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j-2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,1)
    odata[i,j] += (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+2]) / 4 # D2

    odata[i,j] += 2 * alpha4 * (( 1.5* idata[i,j] - 2*idata[i-1,j] + 0.5*idata[i-2,j])) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha1 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j+1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j+2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,Ny)
    odata[i,j] += (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-2]) / 4 # D2

    odata[i,j] += 2 * alpha4 * (1.5 * idata[i,j] - 2*idata[i-1,j] + 0.5 * idata[i-2,j]) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha2 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j-1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j-2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (1,2:Ny-1)
    odata[i,j] .+= (view(idata,i,j) .- 2*view(idata,i+1,j) .+ view(idata,i+2,j) .+ view(idata,i,j .-1) .- 2*view(idata,i,j) .+ view(idata,i,j .+ 1)
                    .+ 2 * alpha3 .* (1.5 * view(idata,i,j) .- 2*view(idata,i+1,j) .+ 0.5 * view(idata,i+2,j))) ./ 2

    (i,j) = (Nx,2:Ny-1)
    odata[i,j] .+= ((view(idata,i,j) .- 2*view(idata,i-1,j) .+ view(idata,i-2,j) .+ view(idata,i,j.-1) .- 2*view(idata,i,j) .+ view(idata,i,j.+1))
                    .+ 2 * alpha4 * (1.5 * view(idata,i,j) .- 2*view(idata,i-1,j) .+ 0.5 * view(idata,i-2,j))) ./ 2

    (i,j) = (2:Nx-1,1)
    odata[i,j] .+= (view(idata,i .- 1,j) .- 2*view(idata,i,j) .+ view(idata,i .+ 1,j) .+ view(idata,i,j) .- 2*view(idata,i,j+1) .+ view(idata,i,j+2)) ./ 2

    odata[i,j] .+= (2 * beta * (1.5 * view(idata,i,j)) .+ 2 * alpha2 * view(idata,i,j) * h) ./ 2
    odata[i,j+1] .+= (2 * beta * (-1 * view(idata,i,j)))
    odata[i,j+2] .+= (0.5 * beta * (view(idata,i,j)))

    (i,j) = (2:Nx-1,Ny)
    odata[i,j] .+= (view(idata,i .- 1, j) .- 2*view(idata,i,j) .+ view(idata,i .+1,j) .+ view(idata,i,j) .- 2*view(idata,i,j-1) .+ view(idata,i,j-2)) ./ 2

    odata[i,j] .+= (2 * beta * (1.5 * view(idata,i,j)) .+ 2 * alpha1 * view(idata,i,j) * h) ./ 2
    odata[i,j-1] .+= (2 * beta * (-1 * view(idata,i,j)))
    odata[i,j-2] .+= (0.5 * beta * (view(idata,i,j)))
end

function matrix_free_cpu_v3(idata,odata,Nx,Ny,h)
    # odata .= 0
    # for i in 1:Nx
    #     for j in 1:Ny
    #         if (i == 1) && (j == 1)
    #             odata[i,j] = (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+1]) / 4
    #         end
    #         if (i == 1) && (j == Ny)
    #             odata[i,j] = (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-1]) / 4
    #         end
    #         if (i == Nx) && (j == 1)
    #             odata[i,j] = (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+1]) / 4
    #         end
    #         # if (i == Nx) && (j == Ny)
    #         #     odata[i,j] = (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-1]) / 4
    #         # end
    #     end
    # end
    odata .= 0
    # alpha1 = alpha2 = alpha3 = alpha4 = beta = 1
    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1
    (i,j) = (1,1)

    odata[i,j] += (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+2]) / 4 # D2

    odata[i,j] += 2 * alpha3 * (( 1.5* idata[i,j] - 2*idata[i+1,j] + 0.5*idata[i+2,j])) / 4 # Neumann

    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha1 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j+1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j+2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (1,Ny)
    odata[i,j] += (idata[i,j] - 2*idata[i+1,j] + idata[i+2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-2]) / 4 # D2

    odata[i,j] += 2 * alpha3 * (1.5 * idata[i,j] - 2*idata[i+1,j] + 0.5 * idata[i+2,j]) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha2 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j-1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j-2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,1)
    odata[i,j] += (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j+1] + idata[i,j+2]) / 4 # D2

    odata[i,j] += 2 * alpha4 * (( 1.5* idata[i,j] - 2*idata[i-1,j] + 0.5*idata[i-2,j])) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha1 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j+1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j+2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,Ny)
    odata[i,j] += (idata[i,j] - 2*idata[i-1,j] + idata[i-2,j] + idata[i,j] - 2*idata[i,j-1] + idata[i,j-2]) / 4 # D2

    odata[i,j] += 2 * alpha4 * (1.5 * idata[i,j] - 2*idata[i-1,j] + 0.5 * idata[i-2,j]) / 4 # Neumann
    odata[i,j] += (2 * beta * (1.5 * idata[i,j]) + 2 * alpha2 * (idata[i,j]) * h) / 4 # Dirichlet
    odata[i,j-1] += (2 * beta * (-1 * idata[i,j])) / 2 # Dirichlet
    odata[i,j-2] += (0.5 * beta * (idata[i,j])) / 2 # Dirichlet

    # (i,j) = (1,2:Ny-1)
    i = 1
    idata_N = view(idata,1:3,1:Ny)
    # Threads.@threads for j in 2:Ny-1
    @inbounds for j in 2:Ny-1
        odata[1,j] += (idata_N[1,j] - 2*idata_N[2,j] + idata_N[3,j] + idata_N[1,j-1] - 2* idata_N[1,j] + idata_N[1,j+1] + 2 * alpha3 * (1.5 * idata_N[1,j] - 2*idata_N[2,j] + 0.5*idata_N[3,j])) / 2
    end
    # synchronize()

    i = Nx
    idata_S = view(idata,Nx-2:Nx,1:Ny)
    # Threads.@threads for j in 2:Ny-1
    @inbounds for j in 2:Ny-1
        odata[i,j] += (idata_S[3,j] - 2*idata_S[2,j] + idata_S[1,j] + idata_S[3,j-1] - 2* idata_S[3,j] + idata_S[3,j+1] + 2 * alpha4 * (1.5 * idata_S[3,j] - 2*idata_S[2,j] + 0.5*idata_S[1,j])) / 2
    end
    # synchronize()

    j = 1
    idata_W = view(idata,1:Nx,1:3)

    # @inbounds for i in 2:Nx-1
    #     odata[i,j] += (idata_W[i-1,1] - 2*idata_W[i,1] + idata_W[i+1,1] + idata_W[i,1] - 2*idata_W[i,2] + idata_W[i,3]) / 2
    #     # odata[i,j] += (2 * beta * (1.5 * idata_W[i,1]) + 2 * alpha2 * idata_W[i,1] * h) / 2
    #     # odata[i,j+1] += (2 * beta * (-1 * idata_W[i,1]))
    #     # odata[i,j+2] += (0.5 * beta * idata_W[i,1])
    # end

    CPU_W_T = copy(idata_W')
    odata_W_T= zeros(size(CPU_W_T))
    # @inbounds for i in 2:Nx-1
    @inbounds for i in 2:Nx-1
        odata_W_T[1,i] += (CPU_W_T[1,i-1] - 2*CPU_W_T[1,i] + CPU_W_T[1,i+1] + CPU_W_T[1,i] - 2*CPU_W_T[2,i] + CPU_W_T[3,i]) / 2
        odata_W_T[1,i] += (2 * beta * (1.5 * CPU_W_T[1,i]) + 2 * alpha2 * CPU_W_T[1,i] * h) / 2
        odata_W_T[2,i] += (2 * beta * (-1 * CPU_W_T[1,i]))
        odata_W_T[3,i] += (0.5 * beta * CPU_W_T[1,i])
    end

    # odata[:,1] .+= odata_W_T[1,:]

    j = Ny
    idata_E = view(idata,1:Nx,Ny-2:Ny)
    CPU_E_T = copy(idata_E')
    CPU_OUT_E_T = zeros(size(CPU_E_T))
    @inbounds for i in 2:Nx-1
        CPU_OUT_E_T[3,i] += (CPU_E_T[3,i-1] - 2*CPU_E_T[3,i] + CPU_E_T[3,i+1] + CPU_E_T[3,i] - 2*CPU_E_T[2,i] + CPU_E_T[1,i]) / 2
        CPU_OUT_E_T[3,i] += (2 * beta * (1.5 * CPU_E_T[3,i]) + 2 * alpha1 * CPU_E_T[3,i] * h) / 2
        CPU_OUT_E_T[2,i] += (2 * beta * (-1 * CPU_E_T[3,i]))
        CPU_OUT_E_T[1,i] += (0.5 * beta * CPU_E_T[3,i])
    end

    # @inbounds for i in 2:Nx-1
    #     odata[i,j] += (idata_E[i-1,3] - 2*idata_E[i,3] + idata_E[i+1,3] + idata_E[i,3] - 2*idata_E[i,2] + idata_E[i,1]) / 2
    #     odata[i,j] += (2 * beta * (1.5 * idata_E[i,3]) + 2 * alpha1 * idata_E[i,3] * h) / 2
    #     odata[i,j-1] += (2 * beta * (-1 * idata_E[i,3]))
    #     odata[i,j-2] += (0.5 * beta * idata_E[i,3])
    # end

    odata[:,1:3] .+= odata_W_T'
    odata[:,end-2:end] .+= CPU_OUT_E_T'
    nothing
end

function matrix_free_cpu_optimized(GPU_Array,odata,Nx,Ny,h)
    # input as idata_gpu
    # Calculation using CPU

    odata .= 0

    CPU_W = zeros(Nx,3)
    CPU_W_T = zeros(3,Nx)
    CPU_E = zeros(Nx,3)
    CPU_E_T = zeros(3,Nx)
    CPU_N = zeros(3,Ny)
    CPU_S = zeros(3,Ny)

    CPU_OUT_W = zeros(Nx,3)
    CPU_OUT_E = zeros(Nx,3)
    CPU_OUT_N = zeros(3,Ny)
    CPU_OUT_S = zeros(3,Ny)

    copyto!(CPU_W,GPU_Array[:,1:3])
    copyto!(CPU_E,GPU_Array[:,end-2:end])
    copyto!(CPU_N,GPU_Array[1:3,:])
    copyto!(CPU_S,GPU_Array[end-2:end,:])

    # copyto!(CPU_W,view(GPU_Array,1:Nx,1:3))
    # copyto!(CPU_E,view(GPU_Array,1:Nx,Ny-2:Ny))
    # copyto!(CPU_N,view(GPU_Array,1:3,:))
    # copyto!(CPU_S,view(GPU_Array,Nx-2:Nx,:))

    CPU_W_T .= CPU_W'
    CPU_E_T .= CPU_E'
    CPU_OUT_W_T = zeros(3,Ny)
    CPU_OUT_E_T = zeros(3,Ny)

    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1

    # (i,j) = (1,2:Ny-1)
    i = 1
    # Threads.@threads for j in 2:Ny-1
    # @inbounds for j in 2:Ny-1
    for j in 2:Ny-1
        CPU_OUT_N[1,j] += (CPU_N[1,j] - 2*CPU_N[2,j] + CPU_N[3,j] + CPU_N[1,j-1] - 2* CPU_N[1,j] + CPU_N[1,j+1] + 2 * alpha3 * (1.5 * CPU_N[1,j] - 2*CPU_N[2,j] + 0.5*CPU_N[3,j])) / 2
    end
    synchronize()

    i = Nx
    # # Threads.@threads for j in 2:Ny-1
    # @inbounds for j in 2:Ny-1
    for j in 2:Ny-1
        CPU_OUT_S[3,j] += (CPU_S[3,j] - 2*CPU_S[2,j] + CPU_S[1,j] + CPU_S[3,j-1] - 2* CPU_S[3,j] + CPU_S[3,j+1] + 2 * alpha4 * (1.5 * CPU_S[3,j] - 2*CPU_S[2,j] + 0.5*CPU_S[1,j])) / 2
    end
    # synchronize()

    j = 1

    # @inbounds for i in 2:Nx-1
    for i in 2:Nx-1
        CPU_OUT_W_T[1,i] += (CPU_W_T[1,i-1] - 2*CPU_W_T[1,i] + CPU_W_T[1,i+1] + CPU_W_T[1,i] - 2*CPU_W_T[2,i] + CPU_W_T[3,i]) / 2
        CPU_OUT_W_T[1,i] += (2 * beta * (1.5 * CPU_W_T[1,i]) + 2 * alpha2 * CPU_W_T[1,i] * h) / 2
        CPU_OUT_W_T[2,i] += (2 * beta * (-1 * CPU_W_T[1,i]))
        CPU_OUT_W_T[3,i] += (0.5 * beta * CPU_W_T[1,i])
    end

    # # odata[:,1] .+= odata_W_T[1,:]

    j = Ny
    # @inbounds for i in 2:Nx-1
    for i in 2:Nx-1
        CPU_OUT_E_T[3,i] += (CPU_E_T[3,i-1] - 2*CPU_E_T[3,i] + CPU_E_T[3,i+1] + CPU_E_T[3,i] - 2*CPU_E_T[2,i] + CPU_E_T[1,i]) / 2
        CPU_OUT_E_T[3,i] += (2 * beta * (1.5 * CPU_E_T[3,i]) + 2 * alpha1 * CPU_E_T[3,i] * h) / 2
        CPU_OUT_E_T[2,i] += (2 * beta * (-1 * CPU_E_T[3,i]))
        CPU_OUT_E_T[1,i] += (0.5 * beta * CPU_E_T[3,i])
    end

    (i,j) = (1,1)

    CPU_OUT_N[1,j] += (CPU_W[i,j] - 2*CPU_W[i+1,j] + CPU_W[i+2,j] + CPU_W[i,j] - 2*CPU_W[i,j+1] + CPU_W[i,j+2]) / 4 # D2

    CPU_OUT_N[1,j] += 2 * alpha3 * (( 1.5* CPU_W[i,j] - 2*CPU_W[i+1,j] + 0.5*CPU_W[i+2,j])) / 4 # Neumann

    CPU_OUT_N[1,j] += (2 * beta * (1.5 * CPU_W[i,j]) + 2 * alpha1 * (CPU_W[i,j]) * h) / 4 # Dirichlet
    CPU_OUT_N[1,j+1] += (2 * beta * (-1 * CPU_W[i,j])) / 2 # Dirichlet
    CPU_OUT_N[1,j+2] += (0.5 * beta * (CPU_W[i,j])) / 2 # Dirichlet

    (i,j) = (1,Ny)
    CPU_OUT_N[1,j] += (CPU_E[i,3] - 2*CPU_E[i+1,3] + CPU_E[i+2,3] + CPU_E[i,3] - 2*CPU_E[i,2] + CPU_E[i,1]) / 4 # D2

    CPU_OUT_N[1,j] += 2 * alpha3 * (1.5 * CPU_E[i,3] - 2*CPU_E[i+1,3] + 0.5 * CPU_E[i+2,3]) / 4 # Neumann
    CPU_OUT_N[1,j] += (2 * beta * (1.5 * CPU_E[i,3]) + 2 * alpha2 * (CPU_E[i,3]) * h) / 4 # Dirichlet
    CPU_OUT_N[1,j-1] += (2 * beta * (-1 * CPU_E[i,3])) / 2 # Dirichlet
    CPU_OUT_N[1,j-2] += (0.5 * beta * (CPU_E[i,3])) / 2 # Dirichlet

    (i,j) = (Nx,1)
    CPU_OUT_S[3,j] += (CPU_W[i,j] - 2*CPU_W[i-1,j] + CPU_W[i-2,j] + CPU_W[i,j] - 2*CPU_W[i,j+1] + CPU_W[i,j+2]) / 4 # D2

    CPU_OUT_S[3,j] += 2 * alpha4 * (( 1.5* CPU_W[i,j] - 2*CPU_W[i-1,j] + 0.5*CPU_W[i-2,j])) / 4 # Neumann
    CPU_OUT_S[3,j] += (2 * beta * (1.5 * CPU_W[i,j]) + 2 * alpha1 * (CPU_W[i,j]) * h) / 4 # Dirichlet
    CPU_OUT_S[3,j+1] += (2 * beta * (-1 * CPU_W[i,j])) / 2 # Dirichlet
    CPU_OUT_S[3,j+2] += (0.5 * beta * (CPU_W[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,Ny)
    CPU_OUT_S[3,j] += (CPU_E[Nx,3] - 2*CPU_E[Nx-1,3] + CPU_E[Nx-2,3] + CPU_E[Nx,3] - 2*CPU_E[Nx,2] + CPU_E[Nx,1]) / 4 # D2

    CPU_OUT_S[3,j] += 2 * alpha4 * (1.5 * CPU_E[Nx,3] - 2*CPU_E[Nx-1,3] + 0.5 * CPU_E[Nx-2,3]) / 4 # Neumann
    CPU_OUT_S[3,j] += (2 * beta * (1.5 * CPU_E[Nx,3]) + 2 * alpha2 * (CPU_E[Nx,3]) * h) / 4 # Dirichlet
    CPU_OUT_S[3,j-1] += (2 * beta * (-1 * CPU_E[Nx,3])) / 2 # Dirichlet
    CPU_OUT_S[3,j-2] += (0.5 * beta * (CPU_E[Nx,3])) / 2 # Dirichlet

    # CPU_OUT_cache = CuArray(CPU_OUT_W_T)'
    # CPU_OUT_cache = CuArray(CPU_OUT_E_T)'
    # CuArray(CPU_OUT_W_T')
    # CuArray(CPU_OUT_E_T')
    # @show CPU_OUT_W_T
    # @show CPU_OUT_E_T
    # CPU_OUT_W .= CPU_OUT_W_T'

    # Copy E & W boundary
    copyto!(view(odata,1:Nx,1:3),CuArray(CPU_OUT_W_T)')
    copyto!(view(odata,1:Nx,Ny-2:Ny),CuArray(CPU_OUT_E_T)')

    # Copy N & S boundary
    copyto!(view(odata,1:1,1:Ny),CuArray(CPU_OUT_N[1,:]))
    copyto!(view(odata,Nx:Nx,1:Ny),CuArray(CPU_OUT_S[end,:]))
    nothing
end

function test_matrix_free_boundary(level)
    Nx = Ny = 2^level+1
    h = 1/(Nx-1)
    Random.seed!(0)
    A = randn(Nx,Ny)
    A_sparse = spzeros(Nx,Ny)
    for i in 1:Nx
        for j in 1:Ny
            if !(( 4 <= i <= Nx - 4) && (4 <= j <= Ny-4))
                A_sparse[i,j] = A[i,j]
            end
        end
    end
    odata = spzeros(Nx,Ny)
    odata_v3 = spzeros(Nx,Ny)

    # # make sure matrix_free_cpu_v3 equals matrix_free_cpu
    # odata_test = spzeros(Nx,Ny)
    # matrix_free_cpu_v3(A,odata_test,Nx,Ny,h)
    # matrix_free_cpu(A,odata,Nx,Ny,h)
    # # @show Array(odata)
    # # @show Array(odata_test)
    # @assert odata ≈ odata_test
    # # end test

    # t_cpu = time()
    # iter_times = 1000
    # for i in 1:iter_times
    #     matrix_free_cpu_v2(A_sparse,odata,Nx,Ny,h)
    #     odata .= 0
    #     # matrix_free_cpu_v3(A,odata,Nx,Ny,h)
    # end
    # t_cpu = time() - t_cpu
    # @show t_cpu

    t_cpu_v3 = time()
    iter_times = 1000
    for i in 1:iter_times
        matrix_free_cpu_v3(A,odata_v3,Nx,Ny,h)
        odata_v3 .= 0
        # matrix_free_cpu_v3(A,odata,Nx,Ny,h)
    end
    synchronize()
    t_cpu_v3 = time() - t_cpu_v3
    @show t_cpu_v3

    t_convert = time()
    for i in 1:iter_times
        cu_out = CUDA.CUSPARSE.CuSparseMatrixCSC(odata)
    end
    synchronize()
    t_convert = time() - t_convert
    @show t_convert
    nothing
end

function test_D2_split(level)
    Nx = Ny = 2^level+1
    h = 1/(Nx-1)
    Random.seed!(0)
    A = randn(Nx,Ny)
    A = sparse(A)
    # cu_A = CUDA.CUSPARSE.CuSparseMatrixCSC(A)
    # cu_A = CUDA.CUSPARSE.CuSparseMatrixCSC(A)
    # cu_out = CUDA.CUSPARSE.CuSparseMatrixCSC(spzeros(Nx,Ny))
    cu_A = CuArray(A)
    cu_out = CuArray(zeros(Nx,Ny))
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)

    @cuda threads=blockdim blocks=griddim D2_split_dev(cu_A,cu_out,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))

    iter_times = 1000

    t_D2_start = time()
    for _ in 1:iter_times
        @cuda threads=blockdim blocks=griddim D2_split_dev(cu_A,cu_out,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
    end
    synchronize()
    t_D2 = time() - t_D2_start

    @show t_D2
    nothing
end

function matrix_free_A_v1(idata,odata)
    Nx,Ny = size(idata)
    h = 1/(Nx-1)
    odata_cpu = spzeros(Nx,Ny)
    # odata_gpu = CUDA.CUSPARSE.CuSparseMatrixCSC(spzeros(Nx,Ny))
    odata_gpu = CuArray(zeros(Nx,Ny))
    # idata_gpu = CUDA.CUSPARSE.CuSparseMatrixCSC(idata)
    idata_gpu = CUDA.CuArray(idata)
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)
    @cuda threads=blockdim blocks=griddim D2_split_dev(idata_gpu,odata_gpu,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
    matrix_free_cpu_v2(idata,odata_cpu,Nx,Ny,h)
    odata .= adapt(Array,odata_gpu) .+ odata_cpu
end

function matrix_free_A_v2(idata,odata)
    Nx,Ny = size(idata)
    h = 1/(Nx-1)
    # odata_gpu = CuArray(zeros(Nx,Ny))
    # idata_gpu = CUDA.CuArray(idata)
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)
    @cuda threads=blockdim blocks=griddim D2_split_dev(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
end

function matrix_free_A_v3(idata,odata,odata_D2_GPU,odata_boundary_GPU)
    # final version for CG
    # odata .= 0
    Nx,Ny = size(idata)
    # odata_D2 = CUDA.zeros(Nx,Ny)
    # odata_boundary = CUDA.zeros(Nx,Ny)
    # odata_boundary = CuArray(zeros(Nx,Ny))
    h = 1/(Nx-1)
    # odata_gpu = CuArray(zeros(Nx,Ny))
    # idata_gpu = CUDA.CuArray(idata)
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)
    @cuda threads=blockdim blocks=griddim D2_split_dev(idata,odata_D2_GPU,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
    matrix_free_cpu_optimized(idata,odata_boundary_GPU,Nx,Ny,h)
    synchronize()
    odata .= odata_D2_GPU .+ odata_boundary_GPU
    # copyto!(view(odata,1:3,:), view(odata_D2_GPU,1:3,:) + view(odata_boundary_GPU,1:3,:))
    # return odata
end

function matrix_free_A_v4(idata,odata)
    # odata .= 0
    Nx,Ny = size(idata)
    h = 1/(Nx-1)
    # odata = CUDA.zeros(Nx,Ny)
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)

    CPU_W = zeros(Nx,3)
    CPU_W_T = zeros(3,Nx)
    CPU_E = zeros(Nx,3)
    CPU_E_T = zeros(3,Nx)
    CPU_N = zeros(3,Ny)
    CPU_S = zeros(3,Ny)

    CPU_OUT_W = zeros(Nx,3)
    CPU_OUT_E = zeros(Nx,3)
    CPU_OUT_N = zeros(3,Ny)
    CPU_OUT_S = zeros(3,Ny)

    copyto!(CPU_W,idata[:,1:3])
    copyto!(CPU_E,idata[:,end-2:end])
    copyto!(CPU_N,idata[1:3,:])
    copyto!(CPU_S,idata[end-2:end,:])
    # copyto!(CPU_W,GPU_Array[:,1:3])
    # copyto!(CPU_E,GPU_Array[:,end-2:end])
    # copyto!(CPU_N,GPU_Array[1:3,:])
    # copyto!(CPU_S,GPU_Array[end-2:end,:])

    @cuda threads=blockdim blocks=griddim D2_split_dev(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))

    CPU_W_T .= CPU_W'
    CPU_E_T .= CPU_E'
    CPU_OUT_W_T = zeros(3,Ny)
    CPU_OUT_E_T = zeros(3,Ny)

    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1

    i = 1
    # Threads.@threads for j in 2:Ny-1
    # @inbounds for j in 2:Ny-1
    for j in 2:Ny-1
        CPU_OUT_N[1,j] += (CPU_N[1,j] - 2*CPU_N[2,j] + CPU_N[3,j] + CPU_N[1,j-1] - 2* CPU_N[1,j] + CPU_N[1,j+1] + 2 * alpha3 * (1.5 * CPU_N[1,j] - 2*CPU_N[2,j] + 0.5*CPU_N[3,j])) / 2
    end
    synchronize()

    i = Nx
    # # Threads.@threads for j in 2:Ny-1
    # @inbounds for j in 2:Ny-1
    for j in 2:Ny-1
        CPU_OUT_S[3,j] += (CPU_S[3,j] - 2*CPU_S[2,j] + CPU_S[1,j] + CPU_S[3,j-1] - 2* CPU_S[3,j] + CPU_S[3,j+1] + 2 * alpha4 * (1.5 * CPU_S[3,j] - 2*CPU_S[2,j] + 0.5*CPU_S[1,j])) / 2
    end
    # synchronize()

    j = 1

    # @inbounds for i in 2:Nx-1
    for i in 2:Nx-1
        CPU_OUT_W_T[1,i] += (CPU_W_T[1,i-1] - 2*CPU_W_T[1,i] + CPU_W_T[1,i+1] + CPU_W_T[1,i] - 2*CPU_W_T[2,i] + CPU_W_T[3,i]) / 2
        CPU_OUT_W_T[1,i] += (2 * beta * (1.5 * CPU_W_T[1,i]) + 2 * alpha2 * CPU_W_T[1,i] * h) / 2
        CPU_OUT_W_T[2,i] += (2 * beta * (-1 * CPU_W_T[1,i]))
        CPU_OUT_W_T[3,i] += (0.5 * beta * CPU_W_T[1,i])
    end

    # # odata[:,1] .+= odata_W_T[1,:]

    j = Ny
    # @inbounds for i in 2:Nx-1
    for i in 2:Nx-1
        CPU_OUT_E_T[3,i] += (CPU_E_T[3,i-1] - 2*CPU_E_T[3,i] + CPU_E_T[3,i+1] + CPU_E_T[3,i] - 2*CPU_E_T[2,i] + CPU_E_T[1,i]) / 2
        CPU_OUT_E_T[3,i] += (2 * beta * (1.5 * CPU_E_T[3,i]) + 2 * alpha1 * CPU_E_T[3,i] * h) / 2
        CPU_OUT_E_T[2,i] += (2 * beta * (-1 * CPU_E_T[3,i]))
        CPU_OUT_E_T[1,i] += (0.5 * beta * CPU_E_T[3,i])
    end

    (i,j) = (1,1)

    CPU_OUT_N[1,j] += (CPU_W[i,j] - 2*CPU_W[i+1,j] + CPU_W[i+2,j] + CPU_W[i,j] - 2*CPU_W[i,j+1] + CPU_W[i,j+2]) / 4 # D2

    CPU_OUT_N[1,j] += 2 * alpha3 * (( 1.5* CPU_W[i,j] - 2*CPU_W[i+1,j] + 0.5*CPU_W[i+2,j])) / 4 # Neumann

    CPU_OUT_N[1,j] += (2 * beta * (1.5 * CPU_W[i,j]) + 2 * alpha1 * (CPU_W[i,j]) * h) / 4 # Dirichlet
    CPU_OUT_N[1,j+1] += (2 * beta * (-1 * CPU_W[i,j])) / 2 # Dirichlet
    CPU_OUT_N[1,j+2] += (0.5 * beta * (CPU_W[i,j])) / 2 # Dirichlet

    (i,j) = (1,Ny)
    CPU_OUT_N[1,j] += (CPU_E[i,3] - 2*CPU_E[i+1,3] + CPU_E[i+2,3] + CPU_E[i,3] - 2*CPU_E[i,2] + CPU_E[i,1]) / 4 # D2

    CPU_OUT_N[1,j] += 2 * alpha3 * (1.5 * CPU_E[i,3] - 2*CPU_E[i+1,3] + 0.5 * CPU_E[i+2,3]) / 4 # Neumann
    CPU_OUT_N[1,j] += (2 * beta * (1.5 * CPU_E[i,3]) + 2 * alpha2 * (CPU_E[i,3]) * h) / 4 # Dirichlet
    CPU_OUT_N[1,j-1] += (2 * beta * (-1 * CPU_E[i,3])) / 2 # Dirichlet
    CPU_OUT_N[1,j-2] += (0.5 * beta * (CPU_E[i,3])) / 2 # Dirichlet

    (i,j) = (Nx,1)
    CPU_OUT_S[3,j] += (CPU_W[i,j] - 2*CPU_W[i-1,j] + CPU_W[i-2,j] + CPU_W[i,j] - 2*CPU_W[i,j+1] + CPU_W[i,j+2]) / 4 # D2

    CPU_OUT_S[3,j] += 2 * alpha4 * (( 1.5* CPU_W[i,j] - 2*CPU_W[i-1,j] + 0.5*CPU_W[i-2,j])) / 4 # Neumann
    CPU_OUT_S[3,j] += (2 * beta * (1.5 * CPU_W[i,j]) + 2 * alpha1 * (CPU_W[i,j]) * h) / 4 # Dirichlet
    CPU_OUT_S[3,j+1] += (2 * beta * (-1 * CPU_W[i,j])) / 2 # Dirichlet
    CPU_OUT_S[3,j+2] += (0.5 * beta * (CPU_W[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,Ny)
    CPU_OUT_S[3,j] += (CPU_E[Nx,3] - 2*CPU_E[Nx-1,3] + CPU_E[Nx-2,3] + CPU_E[Nx,3] - 2*CPU_E[Nx,2] + CPU_E[Nx,1]) / 4 # D2

    CPU_OUT_S[3,j] += 2 * alpha4 * (1.5 * CPU_E[Nx,3] - 2*CPU_E[Nx-1,3] + 0.5 * CPU_E[Nx-2,3]) / 4 # Neumann
    CPU_OUT_S[3,j] += (2 * beta * (1.5 * CPU_E[Nx,3]) + 2 * alpha2 * (CPU_E[Nx,3]) * h) / 4 # Dirichlet
    CPU_OUT_S[3,j-1] += (2 * beta * (-1 * CPU_E[Nx,3])) / 2 # Dirichlet
    CPU_OUT_S[3,j-2] += (0.5 * beta * (CPU_E[Nx,3])) / 2 # Dirichlet

    synchronize()

    # Copy W & E boundary
    copyto!(view(odata,1:Nx,1:3),view(odata,1:Nx,1:3) + CuArray(CPU_OUT_W_T)')
    copyto!(view(odata,1:Nx,Ny-2:Ny),view(odata,1:Nx,Ny-2:Ny) + CuArray(CPU_OUT_E_T)')

    # view(odata,1,1:Ny) + CuArray(CPU_OUT_N[1,:])
    # view(odata,Nx,1:Ny) + CuArray(CPU_OUT_S[end,:])
    # # @show size(view(odata,1:1,1:Ny))
    # # @show size(CuArray(CPU_OUT_N[1,:]))

    # Copy N & S boundary
    copyto!(view(odata,1,1:Ny),view(odata,1,1:Ny) + CuArray(CPU_OUT_N[1,:]))
    copyto!(view(odata,Nx,1:Ny),view(odata,Nx,1:Ny) + CuArray(CPU_OUT_S[end,:]))
    nothing
end

function CG_GPU_dev(b_reshaped_GPU,x_GPU)
    (Nx,Ny) = size(b_reshaped_GPU)
    odata = CUDA.zeros(Nx,Ny)
    # odata_D2_GPU = CUDA.zeros(Nx,Ny)
    # odata_boundary_GPU = CUDA.zeros(Nx,Ny)
    # matrix_free_A_v3(x_GPU,odata,odata_D2_GPU,odata_boundary_GPU)
    matrix_free_A_v4(x_GPU,odata)
    r_GPU = b_reshaped_GPU - odata
    p_GPU = copy(r_GPU)
    rsold_GPU = sum(r_GPU .* r_GPU)
    Ap_GPU = CUDA.zeros(Nx,Ny)
    # num_steps = 0
    # @show rsold_GPU
    for i in 1:Nx*Ny
    # for i in 1:20
        # @show i
        # @show rsold_GPU
        # matrix_free_A_v3(p_GPU,Ap_GPU,odata_D2_GPU,odata_boundary_GPU)
        # num_steps += 1
        matrix_free_A_v4(p_GPU,Ap_GPU)
        alpha_GPU = rsold_GPU / (sum(p_GPU .* Ap_GPU))
        x_GPU .= x_GPU .+ alpha_GPU * p_GPU
        r_GPU .= r_GPU .- alpha_GPU * Ap_GPU
        rsnew_GPU = sum(r_GPU .* r_GPU)
        if sqrt(rsnew_GPU) <  sqrt(eps(real(eltype(b_reshaped_GPU))))
            break
        end
        p_GPU .= r_GPU .+ (rsnew_GPU/rsold_GPU) * p_GPU
        rsold_GPU = rsnew_GPU
        # @show rsold_GPU
    end
    # @show num_steps
end

function CG_CPU_dev(A,b,x)
    r = b - A * x;
    p = r;
    rsold = r' * r
    # Ap = p
    for i = 1:length(b)
    # for i = 1:2
        # @show i
        # @show rsold
        Ap = A * p;
        alpha = rsold / (p' * Ap)
        x .= x .+ alpha * p;
        r .= r .- alpha * Ap;
        rsnew = r' * r
        if sqrt(rsnew) < sqrt(eps(real(eltype(b))))
              break
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew
        # @show rsold
    end
end

function test_CG()
    x = zeros(Nx*Ny)
    x_GPU = CUDA.zeros(Nx,Ny)
    CG_GPU(b_reshaped_GPU,x_GPU)
    CG_CPU(A,b,x)
end

function test_matrix_free_A(level)
    Nx = Ny = 2^level + 1
    h = 1/(Nx-1)
    println("")
    println("Starting Test")
    println("2D Domain Size: $Nx by $Ny")
    Random.seed!(0)
    # idata = sparse(randn(Nx,Ny))
    idata = CuArray(randn(Nx,Ny))
    # odata = spzeros(Nx,Ny)
    odata = CUDA.zeros(Nx,Ny)
    # odata_v4 = CUDA.zeros(Nx,Ny)
    # odata_D2_GPU = CUDA.zeros(Nx,Ny)
    odata_boundary_GPU = CUDA.zeros(Nx,Ny)
    println("Size of the solution matrix (GPU): ", sizeof(idata), " Bytes")

    idata_cpu = zeros(Nx,Ny)
    odata_cpu = spzeros(Nx,Ny)
    odata_D2 = zeros(Nx,Ny)
    copyto!(idata_cpu,idata)
    println("Size of the solution matrix (CPU): ", sizeof(idata_cpu), " Bytes")
    matrix_free_cpu(idata_cpu,odata_cpu,Nx,Ny,h)
    D2_cpu(idata_cpu,odata_D2,Nx,Ny,h)
    odata_cpu + odata_D2

    println("")

    println("Evaluting time to do one A * b")
    println("Timing results in ms")
    println("")

    iter_times = 1000

    #precompile functions

    # odata_boundary = CUDA.zeros(Nx,Ny)

    # matrix_free_A_v3(idata,odata,odata_D2_GPU,odata_boundary_GPU)
    matrix_free_A_v4(idata,odata)
    matrix_free_cpu_v3(idata_cpu,odata_cpu,Nx,Ny,h)
    matrix_free_cpu_optimized(idata,odata_boundary_GPU,Nx,Ny,h)

    @assert odata_cpu ≈ Array(odata_boundary_GPU)

    # Evaluating only D2
    t_start_D2 = time()
    for _ in 1:iter_times
        matrix_free_A_v2(idata,odata)
    end
    synchronize()
    t_D2 = (time() - t_start_D2) * 1000 / iter_times
    @show t_D2
    # End evaluating D2

    # Evaluating only boundary
    t_start_boundary = time()
    for _ in 1:iter_times
        # matrix_free_cpu_v3(idata_cpu,odata_cpu,Nx,Ny,h)
        matrix_free_cpu_v3(idata_cpu,odata_cpu,Nx,Ny,h)
    end
    t_boundary = (time() - t_start_boundary) * 1000 / iter_times
    @show t_boundary
    # End evaluating boundary

    # Evaluating Both in asynchronous way
    t_start_total = time()
    for _ in 1:iter_times
        matrix_free_A_v2(idata,odata)
    end

    for _ in 1:iter_times
        matrix_free_cpu_v3(idata_cpu,odata_cpu,Nx,Ny,h)
    end
    synchronize()
    t_total = (time() - t_start_total) * 1000 / iter_times
    @show t_total
    # End evaluating both in asynchrnous way

    # Test reduction
    t_reduction = time()
    for _ in 1:iter_times
        # reduce(+,idata)
        sum(idata)
        # sum(idata_cpu)
    end
    synchronize()
    t_reduction = (time() - t_reduction ) * 1000 / iter_times
    @show t_reduction
    # End reduction test

    # Evaluating time in Data IO
    t_copy_data = time()
    iter_times_copy_data = 100
    for _ in 1:iter_times_copy_data
        copyto!(idata_cpu,idata)
    end
    # End evaluating time in Data IO
    t_copy_data = ( time() - t_copy_data ) * 1000 / iter_times_copy_data
    @show t_copy_data

    t_boundary_data_part = time()
    iter_times_copy_data = 100
    for _ in 1:iter_times_copy_data
        # matrix_free_cpu_v4(idata,odata_cpu,Nx,Ny,h)
        matrix_free_cpu_optimized(idata,odata_boundary_GPU,Nx,Ny,h)
    end
    # End evaluating time in Data IO
    t_boundary_data_part = ( time() - t_boundary_data_part ) * 1000 / iter_times_copy_data
    @show t_boundary_data_part

     # Evaluating final matrix-free A

    #  t_start_A = time()
    #  for _ in 1:iter_times
    #      matrix_free_A_v3(idata,odata,odata_D2_GPU,odata_boundary_GPU)
    #  end
    #  synchronize()
    #  t_A = (time() - t_start_A) * 1000 / iter_times
    #  @show t_A
     # End evaluating matrix-free A in asynchrnous way

    #  iter_times = 1
    t_start_A_v4 = time()
    for _ in 1:iter_times
         matrix_free_A_v4(idata,odata)
    end
    synchronize()
    t_A_v4 = (time() - t_start_A_v4) * 1000 / iter_times
    @show t_A_v4

    CUDA.unsafe_free!(idata)
    CUDA.unsafe_free!(odata)
    CUDA.unsafe_free!(odata_boundary_GPU)
    nothing
end

# test_matrix_free_A(13)
# test_matrix_free_A(14)
# test_matrix_free_A(15)