text:
"""
    module Mutts
Explorations of a mutable-until-shared discipline in Julia. (MUTable 'Til Shared)
"""
module Mutts

import MacroTools
using MacroTools: postwalk, @capture

#=
Notes

- I've removed the `Base.isimmutable` overload, since I _think_ that's referring to
  specifically the julia "immutable struct" meaning. It may break things to "lie" about
  this. Also, if we do want to support this, we'd have to add an overload for each type
  (exactly like what was done for setproperty!).

TODO

- Need a LICENSE file for this repo
- Need some mechanism to ensure that Mutt types can't be shared
  outside the current Task without being marked immutable.
    - Perhaps by injecting a check into `put!(::Channel, ::Mutt)`,
      or by dasserting in setters if the current task has changed.
- Add special casing for empty structs? They are always immutable so don't need the bool..
- Add check (warning or error?) that all fields of Mutts type must also be either Mutts or
  immutable.
- Should we autogenerate a `copy()` function for @mutt types?
- Try out using different runtime types for mutable & immutable versions instead of
  injecting a boolean: measure performance difference.
=#

export @mutt, ismuttstype, branch!, ismutable, markimmutable!, getmutableversion!

"""
    ismuttstype(t::Type) -> Bool

Returns true if type `t` was created via the `@mutt` macro, meaning it is a Mutable Til
Shared type, and implements the _mutable-until-shared_ discipline. These types start out
mutable, and can be saved in a frozen version via `markimmutable!` and `branch!`.
"""
ismuttstype(t::Type) = mutts_trait(t) == MuttsType()

# --- Trait dispatch for Mutt types -----------------------
# Since Mutts types don't inherit from a common abstract type (to allow them to inherit
# from their own base types, e.g. AbstractDict), we use trait-inheritance to write functions
# that want different methods for Mutts and non-Mutts types.
# (For more on Traits in Julia ("Holy Traits"), see this blog post):
# https://invenia.github.io/blog/2019/11/06/julialang-features-part-2/
struct MuttsType end
struct NonMuttsType end
# Types are consider non-Mutts by default, and the @mutt macro will create an overload for
# `mutts_trait(T)` marking them as MuttsType.
mutts_trait(::Type) = NonMuttsType()
mutts_trait(v::T) where T = mutts_trait(T)

ismutable(obj::T) where T = ismutable(mutts_trait(T), obj)
ismutable(::MuttsType, obj) = obj.__mutt_mutable

getmutableversion!(obj::T) where T = getmutableversion!(mutts_trait(T), obj)
function getmutableversion!(::MuttsType, obj)
    ismutable(obj) ? obj : branch!(obj)
end

"""
    branchactions(obj) = nothing

This callback function is called immediately before a Mutt object is branched from (even if
it was already immutable). Users can add methods to this callback to perform arbitrary
actions right before an object is branched.
"""
branchactions(obj) = nothing

"""
    markimmutable!(obj)

Freeeze `obj` from further mutations, making it eligible to pass to
other Tasks, branch! from it, or otherwise share it.
"""
function markimmutable! end

markimmutable!(o::T) where T = markimmutable!(mutts_trait(T), o)

markimmutable!(::NonMuttsType, a) = a

@generated function markimmutable!(::MuttsType, obj :: T) where T
    as = map(fieldnames(T)) do sym
        :( markimmutable!(getfield(obj, $(QuoteNode(sym)))) )
    end

    return quote
        if ismutable(obj)
            # Mark all Mutt fields immutable
            $(as...)

            # Then mark this object immutable
            obj.__mutt_mutable = false
        end
        obj
    end
end

"""
    branch!(obj)

Return a mutable shallow copy of Mutts object `obj`, whose children are all still immutable.
"""
branch!(obj::T) where T = branch!(mutts_trait(T), obj)
function branch!(::MuttsType, obj)
    branchactions(obj)
    markimmutable!(obj)

    obj = copy(obj)
    obj.__mutt_mutable = true
    obj
end

# Overload setproperty! for Mutts types to throw exception if attempting to modify a Mutt
# once it's been marked immutable.
# NOTE: The @mutt macro will overload setproperty!(obj::T, name, x) to call this method.
function Base.setproperty!(::MuttsType, obj, name::Symbol, x)
    @assert ismutable(obj)
    setfield!(obj, name, x)
end

"""
    @mutt struct MyType
        fields...
    end

Macro to define a _mutable-until-shared_ data type. Mutable-until-shared types
are essentially immutable data structures, that give the flexibility of mutating
them until they are "finished", at which point they are free to be shared with
other Tasks, or other parts of the code.

`Mutt` types act like mutable structs, until the user calls `markimmutable!(obj)`,
after which they act like purely immutable types.

The complete API includes:
 - [`markimmutable!(obj)`](@ref): Freeze `obj`, preventing any future mutations.
 - [`branch!(obj)`](@ref): Make a _mutable_ shallow copy of `obj`.
 - [`getmutableversion!(obj)`](@ref): Return a mutable version of `obj`, either
   `obj` itself if already mutable, or a [`branch!`ed](@ref branch!) copy.
 - [`branchactions(obj::Mutt)`](@ref): Users can override this callback for their
   type with any actions that need to occur when it is branched.

This macro modifies the definition of `S` to include an extra first parameter:
`__mutt_mutable :: Bool`, which tracks at runtime when the value becomes immutable. Inner
constructors are handled automatically, so you not need to construct this generated field.

Example:
```julia
@mutt struct S
    x :: Int
    y
    S(x, y=x+1) = new(x,y)
end
```
"""
macro mutt(expr)
    return _mutt_macro(expr)
end

# Turns `@mutt struct S x end` into:
# ```
# mutable struct S <: Mutt
#     __mutt_mutable :: Bool
#     x
#     ... constructors ...
# end
# ```
function _mutt_macro(expr)
    if MacroTools.isstructdef(expr)
        def = MacroTools.splitstructdef(expr)
        typename = def[:name]
        # Mutts structs are julia-mutable
        def[:mutable] = true
        # Add `__mutt_mutable` field to the struct.
        # (Put our inserted variable first so the user's constructor can leave undefined
        # fields if that's a thing they're into.)
        pushfirst!(def[:fields], (:__mutt_mutable, Bool))
        # Initialize the new `__mutt_mutable` boolean in the construtors.
        if isempty(def[:constructors])
            # Add the default constructor(s), if none exist, to initialize __mutt_mutable.
            append!(def[:constructors], default_constructors(typename, def[:params], def[:fields][2:end]))
        else
            # Inject `true` to initialize __mutt_mutable in `new()` expressions
            def[:constructors] = map(inject_bool_into_constructor!, def[:constructors])
        end
        # Mark this type as a Mutt (Register with the trait dispatch).
        push!(def[:constructors],
              # (Note the <:T which covers the case where T is paramaterized)
              :($(@__MODULE__).mutts_trait(::Type{<:$typename}) = $MuttsType()))
        # Override setproperty! to prevent mutating a Mutt once it's been marked immutable.
        push!(def[:constructors],
        :(function $Base.setproperty!(obj::$typename, name::Symbol, x)
            setproperty!($MuttsType(), obj, name, x)
        end))

        return esc(MacroTools.combinestructdef(def))
   else
       throw(ArgumentError("@mutt macro must be called with a struct definition: @mutt struct S ... end"))
   end
end

function default_constructors(typename, typeparams, fields)
    # TODO: if empty, no bool...?
    typed_args = Tuple(if f[2] == Any f[1] else MacroTools.combinearg(f..., false, nothing) end
                       for f in fields)
    untyped_args = Tuple(f[1] for f in fields)
    function make_constructor(args, argnames)
        ps = typeparams
        if !isempty(ps)
            full_name = :($typename{$(ps...)})
            new_expr = :(new{$(ps...)})
            :(function $full_name($(args...)) where {$(ps...)} ; $new_expr(true, $(argnames...)) end)
        else
            :($typename($(args...)) = new(true, $(argnames...)))
        end
    end
    new_expr = :new
    # If none of the args have types, these will be the same.
    if typed_args == untyped_args
        [make_constructor(untyped_args, untyped_args)]
    else
        [make_constructor(typed_args,   untyped_args),
         make_constructor(untyped_args, untyped_args)]
    end
end
function inject_bool_into_constructor!(constructor)
    function inject_bool_into_new!(expr)
        if @capture(expr, new(args__))
            expr.args = [:new, true, args...]
            expr
        elseif @capture(expr, new{T__}(args__))
            insert!(expr.args, 2, true)
            expr
        else
            expr
        end
    end
    postwalk(inject_bool_into_new!, constructor)
end

end # module


text:
export ValidationPerformance
export register

type ValidationPerformance <: Coffee
  validation_net :: Net

  ValidationPerformance(net::Net) = new(net, Function[])

  # listeners will be notified each time we compute
  # performance on the validation set
  listeners :: Vector{Function}
end

# Listener should be a function that takes the coffee lounge, the
# net and solver state as parameters. The listener can query the coffee lounge
# for any updated statistics.
function register(coffee::ValidationPerformance, listener::Function)
  push!(coffee.listeners, listener)
end

function init(coffee::ValidationPerformance, ::Net)
  init(coffee.validation_net)
end
function enjoy(lounge::CoffeeLounge, coffee::ValidationPerformance, net::Net, state::SolverState)
  reset_statistics(coffee.validation_net)
  forward_epoch(coffee.validation_net)

  m_info("")
  m_info("## Performance on Validation Set after $(state.iter) iterations")
  m_info("---------------------------------------------------------")
  dump_statistics(lounge, coffee.validation_net, show=true)
  m_info("---------------------------------------------------------")
  m_info("")

  # notify listeners
  for listener in coffee.listeners
    listener(lounge, net, state)
  end
end
function destroy(coffee::ValidationPerformance, ::Net)
  # We don't destroy here as we didn't construct the network
  # Let the user destroy the network when they are done.
  # destroy(coffee.validation_net)
end

ext: jl


text:
# A simple example to demonstrate the use of CUDA.jl

using CUDA

try

	println("Devices")
	list_devices()
	println()

	dev = CuDevice(0)
	println("create context")
	ctx = create_context(dev)

	println("load module from vadd.ptx")
	md = CuModule("vadd.ptx")

	println("get function vadd")
	f = CuFunction(md, "vadd")

	siz = (3, 4)
	len = prod(siz)

	println("load array a to GPU")
	a = round(rand(Float32, siz) * 100)
	ga = CuArray(a)

	println("load array b to GPU")
	b = round(rand(Float32, siz) * 100)
	gb = CuArray(b)

	println("create array c on GPU")
	gc = CuArray(Float32, siz)

	println("launch kernel")
	launch(f, len, 1, (ga, gb, gc))

	println("fetch results from GPU")
	c = to_host(gc)

	println("free GPU memory")
	free(ga)
	free(gb)
	free(gc)

	println("Results:")
	println("a = \n$a")
	println("b = \n$b")
	println("c = \n$c")

	println("c == (a + b) ? $(c == (a + b))")

	println("unload module")
	unload(md)

	println("destroy context")
	destroy(ctx)

catch err

	if isa(err, CuDriverError)
		println("$err: $(description(err))")
	else
		throw(err)
	end
end


text:
struct Stop
    stop_lat::Float64
    stop_lon::Float64
end


text:
include("Graph.jl")

using SparseArrays
using LinearAlgebra

function OutDegree(G)
    C = zeros(G.n)
    for (u, v, w) in G.E
        C[u] += w
    end
    return C
end

function spA(G)
    Is = zeros(Int, G.m)
    Js = zeros(Int, G.m)
    Vs = zeros(G.m)
    ID = 0
    for (u, v, w) in G.E
        ID += 1
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = w
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function alpP(G, alp)
    d = OutDegree(G)
    Is = zeros(Int, G.m)
    Js = zeros(Int, G.m)
    Vs = zeros(G.m)
    ID = 0
    for (u, v, w) in G.E
        ID += 1
        Is[ID] = v
        Js[ID] = u
        Vs[ID] = alp * w / d[u]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function PageRank(G; alpha = 0.85)
    aP = alpP(G, alpha)
    C = zeros(G.n)
    C[1] = 1.0
    adC = ((1.0 - alpha) / G.n) * ones(G.n)

    while true
        pC = copy(C)
        C = aP * C + adC
        if norm(C - pC) < 1e-12
            break
        end
    end

    return C
end

function EigenvectorCentrality(G)
    C = zeros(G.n)
    C[1] = 1.0
    A = spA(G)

    while true
        pC = copy(C)
        C = A * C
        C ./= C[argmax(C)]
        if norm(C - pC) < 1e-9
            break
        end
    end

    return C
end

function getPc(C; eps = 1e-9)
    A = copy(C)
    n = size(A, 1)
    sort!(A)
    sum = 0.0
    for i = 1 : n
        ts = 0.0
        for j = i-1 : -1 : 1
            if abs(A[i] - A[j]) < eps
                ts += 1.0
            else
                break
            end
        end
        for j = i+1 : n
            if abs(A[i] - A[j]) < eps
                ts += 1.0
            else
                break
            end
        end
        sum += ts
    end
    sum /= (n * (n - 1))
    return (1.0 - sum) * 100
end


text:
using Test
include("tests.jl")


text:
using Mill: NGramMatrix

"""
	struct ExtractString{T} <: AbstractExtractor
		n::Int
		b::Int
		m::Int
		uniontypes::Bool
	end

Represents `String` as `n-`grams (`NGramMatrix` from `Mill.jl`) with base `b` and modulo `m`.

The `uniontypes` field determines whether extractor may or may not accept `missing`.
If `uniontypes` is false, it does not accept missing values. If `uniontypes` is true, it accepts missing values,
and always returns Mill structure of type Union{Missing, T} due to type stability reasons.

# Example
```jldoctest
julia> using Mill: catobs

julia> ExtractString(true)("hello")
2053×1 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 "hello"

julia> mapreduce(ExtractString(true), catobs, (["hello", "world"]))
2053×2 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 "hello"
 "world"

julia> mapreduce(ExtractString(true), catobs, ["hello", missing])
2053×2 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 "hello"
 missing

julia> ExtractString(true)(missing)
2053×1 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 missing

julia> ExtractString(false)("hello")
2053×1 Mill.ArrayNode{Mill.NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "hello"

julia> mapreduce(ExtractString(false), catobs, (["hello", "world"]))
2053×2 Mill.ArrayNode{Mill.NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "hello"
 "world"

julia> ExtractString(false)(["hello", "world"])
ERROR: This extractor does not support missing values

```
"""
struct ExtractString <: AbstractExtractor
	n::Int
	b::Int
	m::Int
	uniontypes::Bool
end

const MaybeString = Union{Missing, String}

ExtractString(n::Int, b::Int, m::Int) = ExtractString(n, b, m, true)
ExtractString(uniontypes::Bool) = ExtractString(3, 256, 2053, uniontypes)
ExtractString() = ExtractString(3, 256, 2053, true)
make_missing_string(s::ExtractString, v, store_input) =
	s.uniontypes ?
	_make_array_node(NGramMatrix(MaybeString[missing], s.n, s.b, s.m), [v], store_input) :
	error("This extractor does not support missing values")
stabilize_types_string(s::ExtractString, x) = s.uniontypes ? Vector{MaybeString}(x) : x
make_empty_string(s::ExtractString, store_input) where {T} =
	ArrayNode(NGramMatrix(s.uniontypes ? Vector{MaybeString}() : Vector{String}(), s.n, s.b, s.m))

(s::ExtractString)(v::String; store_input=false) =
	_make_array_node(NGramMatrix(stabilize_types_string(s, [v]), s.n, s.b, s.m), [v], store_input)
(s::ExtractString)(v::AbstractString; store_input=false) = s(String(v); store_input)
(s::ExtractString)(v::MissingOrNothing; store_input=false) = make_missing_string(s, v, store_input)
(s::ExtractString)(v::ExtractEmpty; store_input=false) = make_empty_string(s, store_input)
(s::ExtractString)(v; store_input=false) = make_missing_string(s, v, store_input)
(s::ExtractString)(v::Symbol; store_input=false) = s(String(v); store_input)

"""
	extractscalar(Type{String}, n = 3, b = 256, m = 2053)

represents strings as ngrams with
- `n` (the degree of ngram),
- `b` base of string,
- `m` modulo on index of the token to reduce dimension

	extractscalar(Type{Number}, m = 0, s = 1)

extracts number subtracting `m` and multiplying by `s`

# Example
```jldoctest
julia> JsonGrinder.extractscalar(String, 3, 256, 2053, true)("5")
2053×1 Mill.ArrayNode{Mill.NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}, Nothing}:
 "5"

julia> JsonGrinder.extractscalar(Int32, 3, 256, true)("5")
1×1 Mill.ArrayNode{Matrix{Union{Missing, Int32}}, Nothing}:
 512

julia> JsonGrinder.extractscalar(String, 3, 256, 2053, false)("5")
2053×1 Mill.ArrayNode{Mill.NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "5"

julia> JsonGrinder.extractscalar(Int32, 3, 256, false)("5")
1×1 Mill.ArrayNode{Matrix{Int32}, Nothing}:
 512
```
"""
extractscalar(::Type{<:AbstractString}, n = 3, b = 256, m = 2053, uniontypes = true) = ExtractString(n, b, m, uniontypes)

Base.hash(e::ExtractString, h::UInt) = hash((e.n, e.b, e.m, e.uniontypes), h)
Base.:(==)(e1::ExtractString, e2::ExtractString) =
	e1.n === e2.n && e1.b === e2.b && e1.m === e2.m && e1.uniontypes === e2.uniontypes

# to be compatible with the number version
extractscalar(::Type{<:AbstractString}, e::Entry, uniontypes::Bool = true) = ExtractString(uniontypes)


text:
using st16SseqJuliaTools

if abspath(PROGRAM_FILE) == @__FILE__
    uc2jc(ARGS[1])
end


text:
export min

import Base: min
import Base: show

"""
    min(; from = nothing)

Creates a min operator, which emits a single item: the item with the smallest value.

# Arguments
- `from`: optional initial minimal value, if `nothing` first item from the source will be used as initial instead

# Producing

Stream of type `<: Subscribable{Union{L, Nothing}}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ i for i in 1:42 ])
subscribe!(source |> min(), logger())
;

# output

[LogActor] Data: 1
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
min(; from::D = nothing) where D = MinOperator{D}(from)

struct MinOperator{D} <: InferableOperator
    from :: D
end

function on_call!(::Type{L}, ::Type{Union{L, Nothing}}, operator::MinOperator, source) where L
    return proxy(Union{L, Nothing}, source, MinProxy(operator.from !== nothing ? convert(L, operator.from) : nothing))
end

operator_right(operator::MinOperator, ::Type{L}) where L = Union{L, Nothing}

struct MinProxy{D} <: ActorProxy
    from :: D
end

actor_proxy!(::Type{Union{L, Nothing}}, proxy::MinProxy, actor::A) where { L, A } = MinActor{L, A}(actor, MinActorProps{L}(proxy.from))

mutable struct MinActorProps{L}
    current :: Union{L, Nothing}
end

struct MinActor{L, A} <: Actor{L}
    actor   :: A
    props   :: MinActorProps{L}
end

getcurrent(actor::MinActor)         = actor.props.current
setcurrent!(actor::MinActor, value) = actor.props.current = value

function on_next!(actor::MinActor{L}, data::L) where L
    current = getcurrent(actor)
    if current === nothing || data < current
        setcurrent!(actor, data)
    end
end

function on_error!(actor::MinActor, err)
    error!(actor.actor, err)
end

function on_complete!(actor::MinActor)
    next!(actor.actor, getcurrent(actor))
    complete!(actor.actor)
end

Base.show(io::IO, ::MinOperator)         = print(io, "MinOperator()")
Base.show(io::IO, ::MinProxy)            = print(io, "MinProxy()")
Base.show(io::IO, ::MinActor{L}) where L = print(io, "MinActor($L)")


text:
using Plots
using DelimitedFiles
using PyCall
include("../simulation.jl")
np = pyimport("numpy")

const SRC_CONST = 86400.0 # turns seconds into days
const GRAV_CONST = 0.0    # gravity constant
const CONV_CONST = 2000   # linear relation between permeability and porosity

################# --- Run Simulation --- #################

Nz = 1              # 1 cell in z direction
Nx = 101            # 100 cell in x direction   {  assert: typeof( (Nx-1) / (Nx_coloc-1) ) == Int  }
Nt = 2300           # number of timesteps       {  assert: typeof( (Nt) / (Nt_coloc-1) ) == Int  }
Δx = 30.0           # = Δz; cell size in x and z direction [m]
Δt = 0.1            # timespet size in [d]

z = (1:Nz)*Δx|>collect        # array of z coordinates
x = (1:Nx)*Δx|>collect        # array of x coordinates
X, Z = np.meshgrid(x, z)      # mesh grid

ρw = 1000.0     # density if wetting fluid [m^3/d]
ρo = 800.0      # density if produciton fluid [m^3/d]
μw = 1.0        # viscosity if wetting fluid [mPa*s]
μo = 10.7       # viscosity if production fluid [mPa*s]

g = GRAV_CONST  # gravitational acceleration

flow = 0.005 * (1/Δx^2)/10.0 * SRC_CONST # injection rate [m^3/d]

qw = zeros(Nt, Nz, Nx)      # injection rate
qw[:, 1, 1] .= flow
qo = zeros(Nt, Nz, Nx)      # produciton rate
qo[:,1,end] .= -flow

sw0 = zeros(Nz, Nx)         # initial saturation [%]

ϕ_real = ones(Nz, Nx)                # "real" porosity field [%]
ϕ_real[1, :] = ([610.0, 615.0, 620.0, 623.0, 628.0, 631.0, 635.0, 637.0, 638.0,
                640.0, 108.0, 111.0, 113.0, 120.0, 123.0, 125.0, 130.0, 139.0,
                145.0, 151.0, 155.0, 162.0, 515.0, 511.0, 506.0, 503.0, 498.0,
                495.0, 490.0, 910.0, 901.0, 892.0, 883.0, 875.0, 867.0, 859.0,
                852.0, 845.0, 839.0, 710.0, 709.0, 707.0, 704.0, 700.0, 695.0,
                689.0, 682.0, 674.0, 665.0, 655.0, 644.0, 633.0, 622.0, 610.0,
                689.0, 682.0, 674.0, 665.0, 655.0, 644.0, 633.0, 622.0, 610.0,
                610.0, 610.0, 610.0, 610.0, 630.0, 640.0, 650.0, 660.0, 310.0,
                310.0, 310.0, 330.0, 320.0, 330.0, 330.0, 340.0, 345.0, 350.0,
                355.0, 360.0, 365.0, 369.0, 373.0, 377.0, 381.0, 384.0, 386.0,
                388.0, 689.0, 388.0, 386.0, 384.0, 381.0, 377.0, 373.0, 369.0,
                365.0, 360.0])./ CONV_CONST
K_real = ones(Nz, Nx)
K_real[1, :] = ϕ_real[1, :] .* CONV_CONST       # "real" permeability field [milidarcy]

ϕ_init = minimum(ϕ_real) .* ones(Nz, Nx)    # initial porosity guess

# to compute error in longer simulation than what NN predicted change Nt_dob
Nt_dob = Nt
qw_dob = zeros(Nt_dob, Nz, Nx)      # injection rate
qw_dob[:, 1, 1] .= flow
qo_dob = zeros(Nt_dob, Nz, Nx)      # produciton rate
qo_dob[:,1,end] .= -flow

tf_param_real = tfCtxGen(Nz, Nx, Δx, Nt_dob, Δt, Z, X, ρw, ρo, μw, μo, K_real, g, ϕ_real, qw_dob, qo_dob, sw0, true)
tf_sw_real, tf_p_real = imseq(tf_param_real)

println("\nInverting field parameters ...")

# turn all simulation parameters into tf.Tensors
tf_param_init = tfCtxGen(Nz, Nx, Δx, Nt, Δt, Z, X, ρw, ρo, μw, μo, 0, g, ϕ_init, qw, qo, sw0, false)
# create computational graph for simulation
tf_sw_init, tf_p_init = imseq(tf_param_init)

# reading neural network predicted saturation data
sw_pred = readdlm("./results/sw_pred.txt")
sw_pred_tf = constant(reshape(sw_pred, Nt+1, Nz, Nx))

loss = sum((tf_sw_init - sw_pred_tf)^2) / (Nt+1 * Nx)

# Adam optimzer
optimizer_Adam = tf.compat.v1.train.AdamOptimizer(learning_rate=0.01)
train_op = optimizer_Adam.minimize(loss)

sess = Session(); init(sess)

sw_real = run(sess, tf_sw_real)

sw_error = 1.0
ϕ_error = 1.0
K_error = 1.0
ϕ_inv = ϕ_init
K_inv = ϕ_init * CONV_CONST
opt_iter = 0
sw_errs =[]

epochs = 100
losses = []

for i in 1:epochs
    run(sess, train_op)
    if i % 1 == 0
        print("Epoch: ", i, "\t Loss: ", run(sess, loss))

        ϕ_tmp = run(sess, tf_param_init.ϕ)
        K_tmp = run(sess, tf_param_init.K)

        tf_param_inv = tfCtxGen(Nz, Nx, Δx, Nt_dob, Δt, Z, X, ρw, ρo, μw, μo, K_tmp, g, ϕ_tmp, qw_dob, qo_dob, sw0, true)
        tf_sw_inv, tf_p_inv = imseq(tf_param_inv)

        sess_tmp = Session(); init(sess_tmp)
        sw_inv = run(sess_tmp, tf_sw_inv)
        sess_tmp.close()

        sw_error_tmp = sum(@. sqrt((sw_real - sw_inv)^2)) / sum(@. sqrt(sw_real^2))
        println("\tsw error: $sw_error_tmp")

        global sw_error
        global ϕ_error
        global K_error
        global ϕ_inv
        global K_inv
        if sw_error_tmp < sw_error
            ϕ_error_tmp = sum(@. sqrt((ϕ_real - ϕ_tmp)^2)) / sum(@. sqrt(ϕ_real^2))
            K_error_tmp = sum(@. sqrt((K_real - K_tmp)^2)) / sum(@. sqrt(K_real^2))

            sw_error = sw_error_tmp
            ϕ_error = ϕ_error_tmp
            K_error = K_error_tmp

            ϕ_inv = ϕ_tmp
            K_inv = K_tmp
            opt_iter = i
        end

        push!(sw_errs, sw_error_tmp)
    end
    push!(losses, run(sess, loss))
end
# Adam optimzer

sess.close()                                # close simulation session
println("\nInversion finished!")

println("\n\nBest results at epoch $opt_iter")
println("sw error:\t $sw_error")
println("ϕ error:\t $ϕ_error")
println("K error:\t $K_error")

# write results to file

open("./results/inv_results/inv_error.txt", "w") do file
    write(file, "sw error:\t $sw_error")
    write(file, "\nϕ error:\t $ϕ_error")
    write(file, "\nK error:\t $K_error")
    write(file, "\nAdam epochs:\t $(epochs)")
    write(file, "\nBest results at epoch $opt_iter")
end

open("./results/inv_results/inv_poro.txt", "w") do file
    writedlm(file, ϕ_inv)
end

open("./results/inv_results/inv_perm.txt", "w") do file
    writedlm(file, K_inv)
end

# get simulation with best field values
tf_param_inv = tfCtxGen(Nz, Nx, Δx, Nt_dob, Δt, Z, X, ρw, ρo, μw, μo, K_inv, g, ϕ_inv, qw_dob, qo_dob, sw0, true)
tf_sw_inv, tf_p_inv = imseq(tf_param_inv)

sess = Session(); init(sess)
sw_inv = run(sess, tf_sw_inv)
sess.close()

# plot predicted and "real" saturation
for i in 1:Int((Nt_dob)/10):Nt_dob+1
    close("all")
    day = (i-1) * Δt
    println(day)
    Plots.plot(x, sw_real[i, 1, :], label="simulation")
    if i < Nt
        Plots.plot!(x, sw_pred[i, :], label="NN", line=:dash, color=:orange)
    end
    Plots.plot!(x, sw_inv[i, 1, :], label="inverted", line=:dash, color=:green)
    Plots.ylims!((0.0, 1.0))
    name = "./results/inv_results/inv_$(day).png"
    Plots.savefig(name)
end

# plot inverted porosity
close("all")
Plots.plot(x, ϕ_real[1, :], label="Porosity real")
Plots.plot!(x, ϕ_inv[1, :], label="Porosity inverted")
Plots.ylims!((0.0, 1.0))
name = "./results/inv_results/inv_poro_field.png"
Plots.savefig(name)

# plot inverted permeability
close("all")
Plots.plot(x, K_real[1, :], label="Permeability real")
Plots.plot!(x, K_inv[1, :], label="Permeability inverted")
Plots.ylims!((0.0, 1000.0))
name = "./results/inv_results/inv_perm_field.png"
Plots.savefig(name)

# plot inversion losses
close("all")
Plots.plot(1:epochs, losses, label="inversion loss", yaxis=:log)
name = "./results/inv_results/inv_loss.png"
Plots.savefig(name)

# plot inversion error
close("all")
Plots.plot(1:epochs, sw_errs, label="inversion errors", yaxis=:log)
name = "./results/inv_results/inv_errors.png"
Plots.savefig(name)

max_stars_repo_path: 5. calibration heterogeneou single lambda/inversion.jl

max_issues_repo_path: 6. calibration heterogeneous lambda vector/inversion.jl


text:
using MatrixChainMultiply
using Base.Test

# case on readme
function readme_example(use_arrayfire=false)
  if use_arrayfire
    T = AFArray{Float64}
  else
    T = Float64
  end
  a = rand(T, 1000, 1000)
  b = rand(T, 1000, 100)
  c = rand(T, 100, 500)
  d = rand(T, 500)

  # get results
  result1 = *(a,b,c,d)
  result2 = matrixchainmultiply(a,b,c,d)
  difference = result1 - result2

  # do again to get timings
  tic(); *(a,b,c,d); stdtime = toc()
  tic(); matrixchainmultiply(a,b,c,d); chaintime = toc()
  speedup = stdtime / chaintime
  print("Speedup = $speedup")

  difference, speedup
end

difference, speedup = readme_example()

@test norm(difference) < 1e-5
@test speedup > 1.0


text:
@testset "Continuous identity system" begin
    for sd in 1:3
        s = ContinuousIdentitySystem(sd)
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Continuous constrained identity system" begin
    for sd in 1:3
        X = Singleton(ones(sd))
        s = ConstrainedContinuousIdentitySystem(sd, X)
        @test statedim(s) == sd
        @test stateset(s) == X
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Continuous linear system" begin
    for sd in 1:3
        s = LinearContinuousSystem(zeros(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test islinear(s) && isaffine(s)
    end
end

@testset "Continuous affine system" begin
    for sd in 1:3
        s = AffineContinuousSystem(zeros(sd, sd), zeros(sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test !islinear(s) && isaffine(s)
    end
end

@testset "Continuous affine control system with state constraints" begin
    A = zeros(2, 2)
    B = Matrix([0.5 1.5]')
    c = [0.0, 1.0]
    X = Line([1., -1], 0.) # line x = y
    U = Interval(-1.0, 1.0)  # -1 <= u <= 1
    s = ConstrainedAffineControlContinuousSystem(A, B, c, X, U)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 1
    @test inputset(s) == U
    @test !islinear(s) && isaffine(s)
end

@testset "Continuous linear control system" begin
    for sd in 1:3
        s = LinearControlContinuousSystem(zeros(sd, sd), ones(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == sd
        @test islinear(s) && isaffine(s)
    end
end

@testset "Continuous constrained linear system" begin
    A = [1. 1; 1 -1]
    X = Line([1., -1], 0.) # line x = y
    s = ConstrainedLinearContinuousSystem(A, X)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 0
    @test islinear(s) && isaffine(s)
end

@testset "Continuous constrained affine system" begin
    X = Line([1., -1], 0.) # line x = y
    s = ConstrainedAffineContinuousSystem(zeros(2, 2), zeros(2), X)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 0
    @test !islinear(s) && isaffine(s)
end

@testset "Continuous constrained linear control system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    X = Line([1., -1], 0.)
    U = Interval(0.9, 1.1)
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 1
    @test inputset(s) == U
    @test islinear(s) && isaffine(s)

    # initial value problem composite type
    x0 = Singleton([1.5, 2.0])
    p = InitialValueProblem(s, x0)
    @test inputset(p) == U
    @test statedim(p) == 2
    @test inputdim(p) == 1

    # check that the matrices A and B need not be of the same type
    # with A dense and B a lazy adjoint
    B = [0.5 1.5]'
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)

    # check that the matrices A and B need not be of the same type
    # with A sparse and B a lazy adjoint
    A = sparse([1], [2], [1.0], 2, 2) # sparse matrix
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)
end

@testset "Continuous linear algebraic system" begin
    for sd in 1:3
        s = LinearAlgebraicContinuousSystem(zeros(sd, sd), zeros(sd, sd))
        @test islinear(s) && isaffine(s)
        @test statedim(s) == sd
        @test inputdim(s) == 0
    end
end

@testset "Continuous constrained linear algebraic system" begin
    A = [1. 1; 1 -1]
    E = [0. 1; 1 0]
    X = LinearConstraint([0, -1.], 0.) # the set y ≥ 0
    s = ConstrainedLinearAlgebraicContinuousSystem(A, E, X)
    @test islinear(s) && isaffine(s)
    @test statedim(s) == 2
    @test stateset(s) == X
    @test inputdim(s) == 0
end

@testset "Initial value problem for a continuous constrained linear algebraic system" begin
    A = [1. 1; 1 -1]
    E = [0. 1; 1 0]
    X = LinearConstraint([0, -1.], 0.) # the set y ≥ 0
    s = ConstrainedLinearAlgebraicContinuousSystem(A, E, X)
    @test islinear(s) && isaffine(s)
    x0 = Singleton([1.5, 2.0])
    p = IVP(s, x0)
    @test statedim(p) == 2
    @test inputdim(p) == 0
end

@testset "Polynomial system in continuous time" begin
    @polyvar x y
    p = 2x^2 - 3x + y

    # default constructor for scalar p and
    s = PolynomialContinuousSystem(p)
    @test !islinear(s) && !isaffine(s)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test TypedPolynomials.nvariables(s) == 2
    @test TypedPolynomials.variables(s) == (x, y)

    # constructor for scalar p and given number of variables
    s = PolynomialContinuousSystem([p], 2)
end

@testset "Polynomial system in continuous time with state constraints" begin
    @polyvar x y
    p = 2x^2 - 3x + y
    X = BallInf(zeros(2), 0.1)

    # default constructor for scalar p and
    s = ConstrainedPolynomialContinuousSystem(p, X)
    @test !islinear(s) && !isaffine(s)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test dim(stateset(s)) == dim(X)
    @test TypedPolynomials.nvariables(s) == 2
    @test TypedPolynomials.variables(s) == (x, y)
end

function vanderpol!(x, dx)
    dx[1] = x[2]
    dx[2] = x[2] * (1-x[1]^2) - x[1]
    return dx
end

@testset "Continuous system defined by a function" begin
    s = BlackBoxContinuousSystem(vanderpol!, 2)
    f! = (x, dx) -> s.f(x, dx)
    x = [1.0, 0.0]
    dx = similar(x)
    f!(x, dx)
    @test dx ≈ [0.0, -1.0]
end

@testset "Continuous system defined by a function with state constraints" begin
    H = HalfSpace([1.0, 0.0], 0.0) # x <= 0
    s = ConstrainedBlackBoxContinuousSystem(vanderpol!, 2, H)
    f! = (x, dx) -> s.f(x, dx)
    x = [1.0, 0.0]
    dx = similar(x)
    f!(x, dx)
    @test dx ≈ [0.0, -1.0]
    @test stateset(s) == H
end

function vanderpol_controlled!(x, u, dx)
    dx[1] = x[2]
    dx[2] = x[2] * (1-x[1]^2) - x[1] + u[1]
    return dx
end

@testset "Continuous control system defined by a function with state constraints" begin
    H = HalfSpace([1.0, 0.0], 0.0) # x <= 0
    U = Interval(-0.1, 0.1)
    s = ConstrainedBlackBoxControlContinuousSystem(vanderpol_controlled!, 2, 1, H, U)
    f! = (x, u, dx) -> s.f(x, u, dx)
    x = [1.0, 0.0]
    u = an_element(U)
    dx = similar(x)
    f!(x, u, dx)
    @test dx ≈ [0.0, -1.0 + u[1]]
    @test stateset(s) == H
    @test inputdim(s) == 1
    @test inputset(s) == U
end


text:
using Distributions

"""
Example of a computation of beta credible interval with MC sampling.
"""
function mcBetaQuantileDemo()
    S = 47
    N = 100
    a = S+1
    b = (N-S)+1
    α = .05

    SS = 1000
    X = rand(Beta(a, b), SS)
    sort!(X)

    Xl = X[round(Int, SS*α/2)]
    Xu = X[round(Int, SS*(1-α/2))]
    CIhat = [Xl, Xu]

    CIhat2 = quantile(EmpiricalUnivariateDistribution(X), [α/2, 1-α/2])
    (CIhat, CIhat2)
end

"""
Example of a computation of normal credible interval with MC sampling.
"""
function mcNormalQuantileDemo()
    μ = 0
    Σ = 1
    S = 1000
    rn = Normal(μ, Σ)
    X = rand(rn, S)
    qs = [0.025 0.5 0.975]

    qexact = quantile(rn, qs)

    qapprox = quantile(EmpiricalUnivariateDistribution(X), qs)
end


text:
using GLFW
using ModernGL
using CSyntax
using CSyntax.CStatic

using CImGui
using CImGui.LibCImGui
using CImGui.GLFWBackend
using CImGui.OpenGLBackend

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    tmp_glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    tmp_glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    GLFW.WindowHint(GLFW.FOCUSED,GL_TRUE)
end
const glsl_version=tmp_glsl_version

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(1100, 720, "PlayGround")

@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsClassic()

#
# start to build the font:
#

#fonts_dir = joinpath(pathof(CImGui), "..","..","fonts")
fonts_dir = raw"c:\Windows\Fonts";
fonts = CImGui.GetIO().Fonts

builder=ImFontGlyphRangesBuilder_ImFontGlyphRangesBuilder()

#ImFontGlyphRangesBuilder_AddText(builder,"我是中文","")
ImFontGlyphRangesBuilder_AddChar(builder,'π')
ImFontGlyphRangesBuilder_AddChar(builder,'ϕ')
ImFontGlyphRangesBuilder_AddChar(builder,'θ')
ImFontGlyphRangesBuilder_AddChar(builder,'♐')
ImFontGlyphRangesBuilder_AddChar(builder,'☽')
ImFontGlyphRangesBuilder_AddRanges(builder, ImFontAtlas_GetGlyphRangesDefault(fonts))

ranges=ImVector_ImWchar_create()
ImVector_ImWchar_Init(ranges)
ImFontGlyphRangesBuilder_BuildRanges(builder, ranges)

#set MergeMode to true
config=ImFontConfig_ImFontConfig()
let offset=0
    for i in eachindex(Base.datatype_fieldtypes(ImFontConfig))
        fieldnames(ImFontConfig)[i] == :MergeMode ? break : nothing
        offset+=sizeof(Base.datatype_fieldtypes(ImFontConfig)[i])
    end
    unsafe_store!(Ptr{Bool}(config+offset),true,1)
end

r = unsafe_wrap(Vector{ImVector_ImWchar}, ranges, 1)

CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "JuliaMono-Regular.ttf"), 16, config, r[1].Data )
#CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "arial.ttf"), 16, C_NULL, r[1].Data )
#CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "arial.ttf"), 16, C_NULL, CImGui.GetGlyphRangesChineseSimplifiedCommon(fonts) )

#a=[0xf000, 0xf3ff, 0x0]
#a_range=Core.Ptr{ImVector_ImWchar}(Base.pointer_from_objref(a))
#CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "arial.ttf"), 16, config, a_range )

CImGui.Build(fonts)
default_font = CImGui.AddFontDefault(fonts)
#@assert default_font != C_NULL

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    while  !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()

        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()
        CImGui.SetNextWindowSize(CImGui.ImVec2(350.0,200.0),CImGui.ImGuiCond_Once);

        CImGui.Begin("Special strings")
        CImGui.Text("我是中文")
        CImGui.Text("π ϕ θ ♐ ☽")
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
    display_close(display)
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    ImVector_ImWchar_destroy(ranges)
    ImFontGlyphRangesBuilder_destroy(builder)
    ImFontConfig_destroy(config)
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end


text:
include("utils.jl")
include("task.jl")
include("fairness_problems.jl")
include("evaluate.jl")