text:
__precompile__()

module SortingAlgorithms

using Compat
using DataStructures
import Compat.view
using Base.Sort
using Base.Order

import Base.Sort: sort!
import DataStructures: heapify!, percolate_down!

export HeapSort, TimSort, RadixSort

immutable HeapSortAlg  <: Algorithm end
immutable TimSortAlg   <: Algorithm end
immutable RadixSortAlg <: Algorithm end

const HeapSort  = HeapSortAlg()
const TimSort   = TimSortAlg()
const RadixSort = RadixSortAlg()

## Heap sort

function sort!(v::AbstractVector, lo::Int, hi::Int, a::HeapSortAlg, o::Ordering)
    if lo > 1 || hi < length(v)
        return sort!(view(v, lo:hi), 1, length(v), a, o)
    end
    r = ReverseOrdering(o)
    heapify!(v, r)
    for i = length(v):-1:2
        # Swap the root with i, the last unsorted position
        x = v[i]
        v[i] = v[1]
        # The heap portion now ends at position i-1, but needs fixing up
        # starting with the root
        percolate_down!(v,1,x,r,i-1)
    end
    v
end

## Radix sort

# Map a bits-type to an unsigned int, maintaining sort order
uint_mapping(::ForwardOrdering, x::Unsigned) = x
for (signedty, unsignedty) in ((Int8, UInt8), (Int16, UInt16), (Int32, UInt32), (Int64, UInt64), (Int128, UInt128))
    # In Julia 0.4 we can just use unsigned() here
    @eval uint_mapping(::ForwardOrdering, x::$signedty) = reinterpret($unsignedty, xor(x, typemin(typeof(x))))
end
uint_mapping(::ForwardOrdering, x::Float32)  = (y = reinterpret(Int32, x); reinterpret(UInt32, ifelse(y < 0, ~y, xor(y, typemin(Int32)))))
uint_mapping(::ForwardOrdering, x::Float64)  = (y = reinterpret(Int64, x); reinterpret(UInt64, ifelse(y < 0, ~y, xor(y, typemin(Int64)))))

uint_mapping{Fwd}(rev::ReverseOrdering{Fwd}, x) = ~uint_mapping(rev.fwd, x)
uint_mapping{T<:Real}(::ReverseOrdering{ForwardOrdering}, x::T) = ~uint_mapping(Forward, x) # maybe unnecessary; needs benchmark

uint_mapping(o::By,   x     ) = uint_mapping(Forward, o.by(x))
uint_mapping(o::Perm, i::Int) = uint_mapping(o.order, o.data[i])
uint_mapping(o::Lt,   x     ) = error("uint_mapping does not work with general Lt Orderings")

const RADIX_SIZE = 11
const RADIX_MASK = 0x7FF

function sort!(vs::AbstractVector, lo::Int, hi::Int, ::RadixSortAlg, o::Ordering, ts=similar(vs))
    # Input checking
    if lo >= hi;  return vs;  end

    # Make sure we're sorting a bits type
    T = Base.Order.ordtype(o, vs)
    if !isbits(T)
        error("Radix sort only sorts bits types (got $T)")
    end

    # Init
    iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    bin = zeros(UInt32, 2^RADIX_SIZE, iters)
    if lo > 1;  bin[1,:] = lo-1;  end

    # Histogram for each element, radix
    for i = lo:hi
        v = uint_mapping(o, vs[i])
        for j = 1:iters
            idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    swaps = 0
    len = hi-lo+1
    for j = 1:iters
        # Unroll first data iteration, check for degenerate case
        v = uint_mapping(o, vs[hi])
        idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        cbin = cumsum(bin[:,j])
        ci = cbin[idx]
        ts[ci] = vs[hi]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in hi-1:-1:lo
            v = uint_mapping(o, vs[i])
            idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
            ci = cbin[idx]
            ts[ci] = vs[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        for i = lo:hi
            vs[i] = ts[i]
        end
    end
    vs
end

function Base.Sort.Float.fpsort!(v::AbstractVector, ::RadixSortAlg, o::Ordering)
    lo, hi = Base.Sort.Float.nans2end!(v,o)
    sort!(v, lo, hi, RadixSort, o)
end

# Implementation of TimSort based on the algorithm description at:
#
#   http://svn.python.org/projects/python/trunk/Objects/listsort.txt
#   http://en.wikipedia.org/wiki/Timsort
#
# Original author: @kmsquire

const Run = UnitRange{Int}

const MIN_GALLOP = 7

type MergeState
    runs::Vector{Run}
    min_gallop::Int
end
MergeState() = MergeState(Run[], MIN_GALLOP)

# Determine a good minimum run size for efficient merging
# For details, see "Computing minrun" in
# http://svn.python.org/projects/python/trunk/Objects/listsort.txt
function merge_compute_minrun(N::Int, bits::Int)
    r = 0
    max_val = 2^bits
    while N >= max_val
        r |= (N & 1)
        N >>= 1
    end
    N + r
end
merge_compute_minrun(N::Int) = merge_compute_minrun(N, 6)

# Galloping binary search starting at left
# Finds the last value in v <= x
function gallop_last(o::Ordering, v::AbstractVector, x, lo::Int, hi::Int)
    i = lo
    inc = 1
    lo = lo-1
    hi = hi+1
    while i < hi && !lt(o, x, v[i])
        lo = i
        i += inc
        inc <<= 1
    end
    hi = min(i+1, hi)
    # Binary search
    while lo < hi-1
        i = (lo+hi)>>>1
        if lt(o, x, v[i])
            hi = i
        else
            lo = i
        end
    end
    lo
end

# Galloping binary search starting at right
# Finds the last value in v <= x
function rgallop_last(o::Ordering, v::AbstractVector, x, lo::Int, hi::Int)
    i = hi
    dec = 1
    lo = lo-1
    hi = hi+1
    while i > lo && lt(o, x, v[i])
        hi = i
        i -= dec
        dec <<= 1
    end
    lo = max(lo, i-1)
    # Binary search
    while lo < hi-1
        i = (lo+hi)>>>1
        if lt(o, x, v[i])
            hi = i
        else
            lo = i
        end
    end
    lo
end

# Galloping binary search starting at left
# Finds the first value in v >= x
function gallop_first(o::Ordering, v::AbstractVector, x, lo::Int, hi::Int)
    i = lo
    inc = 1
    lo = lo-1
    hi = hi+1
    while i < hi && lt(o, v[i], x)
        lo = i
        i += inc
        inc <<= 1
    end
    hi = min(i+1, hi)
    # Binary search
    while lo < hi-1
        i = (lo+hi)>>>1
        if lt(o, v[i], x)
            lo = i
        else
            hi = i
        end
    end
    hi
end

# Galloping binary search starting at right
# Finds the first value in v >= x
function rgallop_first(o::Ordering, v::AbstractVector, x, lo::Int, hi::Int)
    i = hi
    dec = 1
    lo = lo-1
    hi = hi+1
    while i > lo && !lt(o, v[i], x)
        hi = i
        i -= dec
        dec <<= 1
    end
    lo = max(lo, i-1)
    # Binary search
    while lo < hi-1
        i = (lo+hi)>>>1
        if lt(o, v[i], x)
            lo = i
        else
            hi = i
        end
    end
    hi
end

# Get the next run
# Returns the v range a:b, or b:-1:a for a reversed sequence
function next_run(o::Ordering, v::AbstractVector, lo::Int, hi::Int)
    lo == hi && return lo:hi
    if !lt(o, v[lo+1], v[lo])
        for i = lo+2:hi
            if lt(o, v[i], v[i-1])
                return lo:i-1
            end
        end
        return lo:hi
    else
        for i = lo+2:hi
            if !lt(o, v[i], v[i-1])
                return i-1:-1:lo
            end
        end
        return hi:-1:lo
    end
end

function merge_at(o::Ordering, v::AbstractVector, state::MergeState, n::Integer)
    a = state.runs[n]
    b = state.runs[n+1]
    merge(o,v,a,b,state)
    state.runs[n] = first(a):last(b)
    deleteat!(state.runs, n+1)
    nothing
end

# Merge consecutive runs
# For A,B,C,D = last four lengths, merge_collapse!()
# maintains 3 invariants:
#
#  A > B + C
#  B > C + D
#  C > D
#
# If any of these are violated, a merge occurs to
# correct it
function merge_collapse(o::Ordering, v::AbstractVector, state::MergeState)
    while true
        n = length(state.runs)
        n <= 1 && break

        # Check invariants 1 and 2
        if (n >= 3 && length(state.runs[end-2]) <= length(state.runs[end-1]) + length(state.runs[end])) ||
            (n >= 4 && length(state.runs[end-3]) <= length(state.runs[end-2]) + length(state.runs[end-1]))

            if length(state.runs[end-2]) < length(state.runs[end])
                merge_at(o,v,state,n-2)
            else
                merge_at(o,v,state,n-1)
            end

        # Check invariant 3
        elseif length(state.runs[end-1]) <= length(state.runs[end])
            merge_at(o,v,state,n-1)

        else # Invariant is satisfied
            break
        end
    end
end

# Merge runs a and b in vector v
function merge(o::Ordering, v::AbstractVector, a::Run, b::Run, state::MergeState)

    # First elements in a <= b[1] are already in place
    a = gallop_last(o, v, v[first(b)], first(a), last(a))+1: last(a)

    if length(a) == 0  return  end

    # Last elements in b >= a[end] are already in place
    b = first(b) : rgallop_first(o, v, v[last(a)], first(b), last(b))-1

    # Choose merge_lo or merge_hi based on the amount
    # of temporary memory needed (smaller of a and b)
    if length(a) < length(b)
        merge_lo(o, v, a, b, state)
    else
        merge_hi(o, v, a, b, state)
    end
end

# Merge runs a and b in vector v (a is smaller)
function merge_lo(o::Ordering, v::AbstractVector, a::Run, b::Run, state::MergeState)

    # Copy a
    v_a = v[a]

    # Pointer into (sub)arrays
    i = first(a)
    from_a = 1
    from_b = first(b)

    mode = :normal
    while true
        if mode == :normal
            # Compare and copy element by element
            count_a = count_b = 0
            while from_a <= length(a) && from_b <= last(b)
                if lt(o, v[from_b], v_a[from_a])
                    v[i] = v[from_b]
                    from_b += 1
                    count_a = 0
                    count_b += 1
                else
                    v[i] = v_a[from_a]
                    from_a += 1
                    count_a += 1
                    count_b = 0
                end
                i += 1

                # Switch to galloping mode if a string of elements
                # has come from the same set
                if count_b >= state.min_gallop || count_a >= state.min_gallop
                    mode = :galloping
                    break
                end
            end
            # Finalize if we've exited the loop normally
            if mode == :normal
                mode = :finalize
            end
        end

        if mode == :galloping
            # Use binary search to find range to copy
            while from_a <= length(a) && from_b <= last(b)
                # Copy the next run from b
                b_run = from_b : gallop_first(o, v, v_a[from_a], from_b, last(b)) - 1
                i_end = i + length(b_run) - 1
                v[i:i_end] = v[b_run]
                i = i_end + 1
                from_b = last(b_run) + 1

                # ... then copy the first element from a
                v[i] = v_a[from_a]
                i += 1
                from_a += 1

                if from_a > length(a) || from_b > last(b) break end

                # Copy the next run from a
                a_run = from_a : gallop_last(o, v_a, v[from_b], from_a, length(a))
                i_end = i + length(a_run) - 1
                v[i:i_end] = v_a[a_run]
                i = i_end + 1
                from_a = last(a_run) + 1

                # ... then copy the first element from b
                v[i] = v[from_b]
                i += 1
                from_b += 1

                # Return to normal mode if we haven't galloped...
                if length(a_run) < MIN_GALLOP && length(b_run) < MIN_GALLOP
                    mode = :normal
                    break
                end
                # Reduce min_gallop if this gallop was successful
                state.min_gallop -= 1
            end
            if mode == :galloping
                mode = :finalize
            end
            state.min_gallop = max(state.min_gallop,0) + 2  # penalty for leaving gallop mode
        end

        if mode == :finalize
            # copy end of a
            i_end = i + (length(a) - from_a)
            v[i:i_end] = v_a[from_a:end]
            break
        end
    end
end

# Merge runs a and b in vector v (b is smaller)
function merge_hi(o::Ordering, v::AbstractVector, a::Run, b::Run, state::MergeState)

    # Copy b
    v_b = v[b]

    # Pointer into (sub)arrays
    i = last(b)
    from_a = last(a)
    from_b = length(b)

    mode = :normal
    while true
        if mode == :normal
            # Compare and copy element by element
            count_a = count_b = 0
            while from_a >= first(a) && from_b >= 1
                if !lt(o, v_b[from_b], v[from_a])
                    v[i] = v_b[from_b]
                    from_b -= 1
                    count_a = 0
                    count_b += 1
                else
                    v[i] = v[from_a]
                    from_a -= 1
                    count_a += 1
                    count_b = 0
                end
                i -= 1

                # Switch to galloping mode if a string of elements
                # has come from the same set
                if count_b >= state.min_gallop || count_a >= state.min_gallop
                   mode = :galloping
                   break
                end
            end
            # Finalize if we've exited the loop normally
            if mode == :normal
                mode = :finalize
            end
        end

        if mode == :galloping
            # Use binary search to find range to copy
            while from_a >= first(a) && from_b >= 1
                # Copy the next run from b
                b_run = rgallop_first(o, v_b, v[from_a], 1, from_b) : from_b
                i_start = i - length(b_run) + 1
                v[i_start:i] = v_b[b_run]
                i = i_start - 1
                from_b = first(b_run) - 1

                # ... then copy the first element from a
                v[i] = v[from_a]
                i -= 1
                from_a -= 1

                if from_a < first(a) || from_b < 1 break end

                # Copy the next run from a
                a_run = rgallop_last(o, v, v_b[from_b], first(a), from_a) + 1: from_a
                i_start = i - length(a_run) + 1
                v[i_start:i] = v[a_run]
                i = i_start - 1
                from_a = first(a_run) - 1

                # ... then copy the first element from b
                v[i] = v_b[from_b]
                i -= 1
                from_b -= 1

                # Return to normal mode if we haven't galloped...
                if length(a_run) < MIN_GALLOP && length(b_run) < MIN_GALLOP
                    mode = :normal
                    break
                end
                # Reduce min_gallop if this gallop was successful
                state.min_gallop -= 1
            end
            if mode == :galloping
                mode = :finalize
            end
            state.min_gallop = max(state.min_gallop, 0) + 2  # penalty for leaving gallop mode
        end

        if mode == :finalize
            # copy start of b
            i_start = i - from_b + 1
            v[i_start:i] = v_b[1:from_b]
            break
        end
    end
end

# TimSort main method
function sort!(v::AbstractVector, lo::Int, hi::Int, ::TimSortAlg, o::Ordering)
    minrun = merge_compute_minrun(hi-lo+1)
    state = MergeState()
    i = lo
    while i <= hi
        run_range = next_run(o, v, i, hi)
        count = length(run_range)
        if count < minrun
            # Make a run of length minrun
            count = min(minrun, hi-i+1)
            run_range = i:i+count-1
            sort!(v, i, i+count-1, DEFAULT_STABLE, o)
        else
            if !issorted(run_range)
                run_range = last(run_range):first(run_range)
                reverse!(view(v, run_range))
            end
        end
        # Push this run onto the queue and merge if needed
        push!(state.runs, run_range)
        i = i+count
        merge_collapse(o, v, state)
    end
    # Force merge at the end
    while true
        n = length(state.runs)
        n <= 1 && break
        merge_at(o, v, state, n-1)
    end
    return v
end

end # module


text:
@testset "Test default Grover search operation 5 qubits" begin
    num_qubits = 5
    use_aux_qubits = false
    bit_pattern = 10

    qq_cct = QuantExQASM.Algorithms.create_grover_circuit(num_qubits, use_aux_qubits, bit_pattern)

    @test begin
        sv = get_statevector_using_qiskit( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
    @test begin
        sv = get_statevector_using_picoquant( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
end

@testset "Test optimised Grover search operation 5+1+3 qubits (Ctrl, Tgt, Aux)" begin
    num_qubits = 9
    use_aux_qubits = true
    bit_pattern = 10

    # Convert circuit to Qiskit circuit
    qq_cct = QuantExQASM.Algorithms.create_grover_circuit(num_qubits, use_aux_qubits, bit_pattern)

    @test begin
        sv = get_statevector_using_qiskit( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
    @test begin
        sv = get_statevector_using_picoquant( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
end

@testset "Test default Grover search operation 6 qubits" begin
    num_qubits = 6
    use_aux_qubits = false
    bit_pattern = 37

    qq_cct = QuantExQASM.Algorithms.create_grover_circuit(num_qubits, use_aux_qubits, bit_pattern)

    @test begin
        sv = get_statevector_using_qiskit( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
    @test begin
        sv = get_statevector_using_picoquant( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
end

@testset "Test optimised Grover search operation 6+1+4 qubits (Ctrl, Tgt, Aux)" begin
    num_qubits = 11
    use_aux_qubits = true
    bit_pattern = 10

    qq_cct = QuantExQASM.Algorithms.create_grover_circuit(num_qubits, use_aux_qubits, bit_pattern)

    @test begin
        sv = get_statevector_using_qiskit( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
    @test begin
        sv = get_statevector_using_picoquant( qq_cct, big_endian=true )
        rho = abs.(sv).^2
        rho[bit_pattern+1] > 0.8
    end
end

ext: jl


text:
using OptiMo, Bazinga, ScSTO
using Printf
using PyPlot, PyCall

################################################################################
################################################################################
# unconstrained STO
################################################################################
################################################################################
# time interval
t0 = 0.0
tf = 12.0

# control sequence [nu x N]
uvec = Matrix([repeat([0.0; 1.0], 5, 1); 0.0]')
nu = size(uvec, 1)
N = size(uvec, 2)

# cost matrix
C = Matrix([1.0 0.0 -1.0 0.0; 0.0 1.0 0.0 -1.0])
Q = Matrix(C' * C)

# initial state
x0 = [0.5; 0.7; 1; 1] # nx x 1
nx = length(x0)

# system dynamics
# nx, nu  ->  nx
# nx, nu  -> nx x nx
function dynam(x::Vector{Float64}, u::Vector{Float64})
    n = length(x)
    f = zeros(n)
    f[1] = x[1] - x[1] * x[2] - 0.4 * x[1] * u[1]
    f[2] = -x[2] + x[1] * x[2] - 0.2 * x[2] * u[1]
    return f
end
function d_dynam(x::Vector{Float64}, u::Vector{Float64})
    n = length(x)
    df = zeros(n, n)
    df[1, 1] = 1.0 - x[2] - 0.4 * u[1]
    df[1, 2] = -x[1]
    df[2, 1] = x[2]
    df[2, 2] = -1.0 + x[1] - 0.2 * u[1]
    return df
end

# switching cost
swc = 0.1

# number of points on the fixed grid
ngrid_grid = [100; 200; 1000]
ng = length(ngrid_grid)

# solver
solver = Bazinga.ZEROFPR( max_iter = 50, verbose = true )

# allocation
nt = 2000;
swdelta = Array{Float64}(undef, N, ng);
swtau = Array{Float64}(undef, N - 1, ng);
xsim = Array{Float64}(undef, nx, nt, ng);
usim = Array{Float64}(undef, nu, nt, ng);

for k = 1:ng

    ngrid = ngrid_grid[k]

    prob = scstoproblem(
        x0,
        dynam,
        d_dynam,
        uvec,
        ngrid = ngrid,
        t0 = t0,
        tf = tf,
        Q = Q,
        swc = swc,
    )

    out = solver(prob)
    print(out)

    swdelta[:, k] .= out.x
    swtau[:, k], _ = gettau(prob, swdelta[:, k])

    xsim[:, :, k], _, _, t = simulate(prob, swtau[:, k], length = nt)
    usim[:, :, k], _ = simulateinput(prob, swtau[:, k], t)

    @printf "\n\n"

end

t = collect(range(t0, stop = tf, length = nt))

figure()
subplot(3, 1, 1)
for k = 1:ng
    plot(t, xsim[1, :, k])
end
ylim(0, 1.75)
xlim(0, 12)
yticks([0; 1])
ylabel(L"x_1")
subplot(3, 1, 2)
for k = 1:ng
    plot(t, xsim[2, :, k])
end
ylabel(L"x_2")
ylim(0, 1.75)
xlim(0, 12)
yticks([0; 1])
subplot(3, 1, 3)
for k = 1:ng
    plot(t, usim[1, :, k], label = "n = $(ngrid_grid[k])")
end
legend()
ylim(-0.2, 1.2)
yticks([0; 1])
xlim(0, 12)
ylabel(L"u")
xlabel(L"$\mathrm{Time}\; [s]$")
gcf()
#savefig(foldername * "ScSTO.jl/test/data/fishing_ngrid.pdf")

################################################################################
# make switching times vectors comparable and compare results

using LinearAlgebra

function cleandelta(delta::Vector{Float64})
    # !!! for binary controls only !!!
    # uvec = 0, 1, 0, 1, 0, 1, ...
    deltanew = copy( delta )
    n = length( delta )
    for k in 1:n-2
        if deltanew[k] > 0.0 && deltanew[k+1] == 0.0
            deltanew[k+2] += deltanew[k]
            deltanew[k] = 0.0
        end
    end
    return deltanew
end
function cleandelta(delta::Matrix{Float64})
    deltanew = copy( delta )
    m = size( delta,2 )
    for k in 1:m
        deltanew[:,k] .= cleandelta( delta[:,k] )
    end
    return deltanew
end

swdelta_clean = cleandelta( swdelta )
global swdelta_diff = Array{Float64}(undef, ng-1);
for k in 1:ng-1
    swdelta_diff[k] = norm( swdelta_clean[:,k] - swdelta_clean[:,ng] , Inf )
end


text:
using GaussianDistributions
using CovarianceFunctions: GradientKernel, ValueGradientKernel

k = CovarianceFunctions.EQ()
g = GradientKernel(k)


text:
module FieldsMocksTests

using Gridap

using ..FieldsMocks

f = MockField(2,Float64)

p1 = Point(1,2)
p2 = Point(0,3)
p3 = Point(7,4)
p4 = Point(8,1)
p = [p1,p2,p3,p4]
v = Float64[1,0,7,8]
g1 = VectorValue(1.0,0.0)
g = fill(g1,4)

test_field(f,p,v,g)

b = MockBasis(2,Float64)

v = zeros(Float64,3,4)
v1 = Float64[1,0,7,8]
v[1,:] = 1*v1
v[2,:] = 2*v1
v[3,:] = 3*v1

g = zeros(VectorValue{2,Float64},3,4)
g1 = fill(g1,4)
g[1,:] = 1*g1
g[2,:] = 2*g1
g[3,:] = 3*g1

test_basis(b,p,v,g)

f = MockGeomap(2,Int)
v= 3 .* p

g1 = 3*TensorValue(1,0,0,1)
g = fill(g1,4)

test_field(f,p,v,g)

end # module


text:
doc"""
`radpoisson(r,rho,l)`

solves the poisson equation for a radial density `rho` in
the `l`-angular momentum channel

computes the hartree potential

$v_{h,lm}(r) = \int_0^r dr_<\, r_<^2 g_l(r_<,r) \rho_{lm}(r_<) +
	\int_r^\infty dr_>\, r_>^2 g_l(r,r_>) \rho_{lm}(r_>)$,

where $g_l(r_<,r_>) = r_<^l/r_>^{l+1}$

uses a 3rd order adams-moulton multistep integrator to compute
each integral

`r` (LogarithmicGrid) - logarithmic grid (`loggrid`)

`rho` (AbstractVector) - radial density $\rho_{lm}(r)$ evaluated on `r`

`l` (Int) - l-quantum number

e.g.
```julia
l = 0
r = loggrid(-13,5,100)
rho = exp(-2r)
vh = radpoisson(r,rho,l)
q = multipolemoment(r,rho,l)
plot(r,vh)
plot!(r,q./r.^(l+1))  # asymptotic form of hartree potential
xlims!(0,10)
ylims!(0,0.3)
```

[`V. Blum et al., CPC 180 (2009)` p. 2185-2186]
"""
function radpoisson(r::LogarithmicGrid, rho::AbstractVector, l::Int)
    a = similar(r)
    b = similar(r)
    poisson_multistep_a(a,r,rho,r.dx,l)
    poisson_multistep_b(b,reverse(r),reverse(rho),-r.dx,l)
    a + reverse(b)
end

@adamsmoulton3_multistep 	poisson_multistep_a 	-(l+1)*y[j] + r[j]^2*rho[j]
@adamsmoulton3_multistep 	poisson_multistep_b 	     l*y[j] - r[j]^2*rho[j]

doc"""
`multipolemoment(r,rho,l)`

computes the multipole moment of radial density `rho` in
the `l`-angular momentum channel

$q = \int dr\, r^{l+2} \rho_{lm}(r)$.

`r` (LogarithmicGrid) - logarithmic grid (`loggrid`)

`rho` (AbstractVector) - radial density $\rho_{lm}(r)$ evaluated on `r`

`l` (Int) - l-quantum number

"""
multipolemoment(r,rho,l) = integrate(r,r.^(l+2).*rho)

# @fastmath function poisson_multistep_rule(f_jm1::Float64, f_j::Float64, f_jp1::Float64, y_jm1::Float64, y_j::Float64, h::Float64, c::Float64)
#     (h*h*(f_jp1 + 10f_j + f_jm1 + c*(10y_j + y_jm1)) + 24y_j - 12y_jm1)/(12 - c*h*h)
# end

# @fastmath @inbounds function poisson_multistep_integrator!(y::Vector{Float64}, f::Vector{Float64}, h::Float64,
#                                                      y1::Float64, y2::Float64, c::Float64)
#     y[1] = y1
#     y[2] = y2
#     for j in 2:length(f)-1
#         y[j+1] = poisson_multistep_rule(f[j-1], f[j], f[j+1], y[j-1], y[j], h, c)
#     end
# end

# @fastmath @inbounds function poisson_multistep_integrator(f::Vector{Float64}, h::Float64, y1::Float64, y2::Float64, c::Float64)
#     y = similar(f)
#     poisson_multistep_integrator!(y,f,h,y1,y2,c)
#     y
# end

# @multistep poisson_multistep_integrator y[j+1] = (h*h*(f[j+1] + 10f[j] + f[j-1] + c*(10y[j] + y[j-1])) + 24y[j] - 12y[j-1])/(12 - c*h*h)

# @multistep numerov_multistep y[j+1] = ( (2-5h*h*g[j]/6)*y[j] - (1+h*h*g[j-1]/12)*y[j-1] + h*h/12*(s[j+1]+10s[j]+s[j-1]) ) / (1+h*h*g[j+1]/12)

# function radpoisson(r::LogarithmicGrid, rho::AbstractVector)
#     const c = 0.25
#     N = integrate(r, 4pi*r.^2.*rho)
#     f = -4pi.*r.^2.5.*rho
#     U = similar(r)
#     poisson_multistep_integrator!(U, 0., 1e-16, f, c, r.dx)
#     U .*= sqrt(r)
#     R = r[end]
#     U += (N - U[end])/R.*r
#     vh = U./r
#     vh/4pi
# end

# @fastmath function multistep_rule1(f_jm1::Float64, f_j::Float64, f_jp1::Float64, y_jm1::Float64, y_j::Float64, h::Float64)
#     2y_j - y_jm1 + h*h*(f_jp1 + 10f_j + f_jm1)/12
# end

# @fastmath @inbounds function multistep_integrator1!(y::Vector{Float64}, f::Vector{Float64},
#                                                       h::Float64, y1::Float64, y2::Float64)
#     y[1] = y1
#     y[2] = y2
#     for j in 2:length(f)-1
#         y[j+1] = multistep_rule1(f[j-1], f[j], f[j+1], y[j-1], y[j], h)
#     end
# end

# @fastmath @inbounds function multistep_integrator1(f::Vector{Float64}, h::Float64, y1::Float64, y2::Float64)
#     y = similar(f)
#     multistep_integrator1!(y,f,h,y1,y2)
#     y
# end

# uniform radial grid
# function poissonrad(rho,r)
#     N = integrate(r, 4pi*r.^2.*rho)
#     dr = r[1] - r[0]
#     f = -4pi.*r.*rho
#     U = similar(r)
#     multistep_integrator1!(U, f, dr, 0., 1e-16)
#     R = r[end]
#     U += (N - U[end])/R.*r
#     vh = U./r
#     vh
# end

# function hartreeatomicrad(r::Grid, rho::AbstractArray, l::Int)
#     cumsum(r.w.*r.^(2+l).*rho)./r.^(l+1) + r.^l.*flipdim(cumsum(flipdim((r.w.*r.^(2-l-1).*rho),1)),1)
# end

max_stars_repo_name: johncsnyder/Element


text:
include("src/A1.jl")
using .A1

using GRUtils

a = 0
b = 1

p = Problem(
    (x) -> 2 / √(x + 1),
    (x) -> -2 * √(x^2 + 1),
    (x) -> 10 / (1 + x * (1 - x)),
    a,
    0,
    0.5,
    b,
    -1,
    0,
)

open("result", "w") do io
    println(io, solve(p))
end

for n in (10, 100)
    h = (b - a) / n
    x = [a + i * h for i in 0:n]
    y = solve(p, Options(n))
    if n == 10
        plot(x, y, "+", x, y, "b--", hold=true)
    else
        plot(x, y, hold=true)
    end
end

legend("n = 10", "n = 10", "n = 100", location = "upper left")
savefig("figures/comparison.pdf")


text: require("core.jl")


text:
export CuSubArray, CuSubVector, CuSubMatrix, CuSubVecOrMat

struct CuSubArray{T,N}
    parent::CuArray{T}
    indices::NTuple{N,UnitRange{Int}}
    offset::Int
    linear::Bool
end

const CuSubVector{T} = CuSubArray{T,1}
const CuSubMatrix{T} = CuSubArray{T,2}
const CuSubVecOrMat{T} = Union{CuSubVector{T},CuSubMatrix{T}}

CuArray(x::CuSubArray) = copyto!(similar(x), x)
Base.Array(x::CuSubArray) = Array(CuArray(x))
Base.similar(x::CuSubArray{T}) where T = CuArray{T}(size(x))
Base.size(x::CuSubArray) = map(length, x.indices)
Base.size(x::CuSubArray, i::Int) = length(x.indices[i])
Base.strides(x::CuSubArray) = strides(x.parent)
Base.stride(x::CuSubArray, i::Int) = stride(x.parent, i)
Base.length(x::CuSubArray) = prod(size(x))
Base.pointer(x::CuSubArray, index::Int=1) = pointer(x.parent, x.offset+index)
Base.cconvert(::Type{Ptr{T}}, x::CuSubArray) where T = Base.cconvert(Ptr{T}, pointer(x))

function Base.view(x::CuArray{T,N}, I...) where {T,N}
    @assert length(I) == N
    I = UnitRange.(Base.to_indices(x,I))
    offset = 0
    strides = Base.strides(x)
    for i = 1:length(strides)
        offset += strides[i] * (first(I[i])-1)
    end
    linear = true
    for i = 1:length(I)-1
        length(I[i]) == size(x,i) && continue
        linear = false
        break
    end
    CuSubArray(x, I, offset, linear)
end

##### IO #####
Base.show(io::IO, x::CuSubArray) = print(io, Array(x))

@generated function Base.copyto!(dest::CuArray{T}, src::CuSubArray{T,N}) where {T,N}
    Ct = cstring(T)
    k = Kernel("""
    __global__ void copyto($Ct *dest, Array<$Ct,$N> src) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= src.length()) return;
        dest[idx] = src(idx);
    }""")
    quote
        @assert length(dest) == length(src)
        src.linear && return copyto!(dest, 1, src.parent, src.offset+1, length(dest))
        gdims, bdims = cudims(length(dest))
        $k(gdims, bdims, pointer(dest), src)
        dest
    end
end

@generated function Base.copyto!(dest::CuSubArray{T,N}, src::CuArray{T}) where {T,N}
    Ct = cstring(T)
    k = Kernel("""
    __global__ void copyto(Array<$Ct,$N> dest, $Ct *src) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= dest.length()) return;
        dest(idx) = src[idx];
    }""")
    quote
        @assert length(dest) == length(src)
        dest.linear && return copyto!(dest.parent, dest.offset+1, src, 1, length(src))
        gdims, bdims = cudims(length(dest))
        $k(gdims, bdims, dest, pointer(src))
        dest
    end
end


text:
__precompile__(true)

module MLSuite

using Statistics, Random, Distributed, Printf, LinearAlgebra, DelimitedFiles, Mmap
using DataStructures, Parameters, ProgressMeter, PyCall, BSON, Requires
using PyCallUtils, PandasLite, HDF5Utils, SVMLightWriter, Iconv
using ScikitLearnBase: BaseEstimator, BaseClassifier, BaseRegressor
import ScikitLearnBase: fit!, predict, predict_proba, is_classifier
import MLSuiteBase: gridparams, reset!, istree, isrnn, is_ranker, support_multiclass, modelhash, signone

const JULIA = joinpath(Sys.BINDIR, Base.julia_exename())
const VW = joinpath(DEPOT_PATH[1], "mlsuite/bin/vw")
const VW_VARINFO = joinpath(DEPOT_PATH[1], "mlsuite/bin/vw-varinfo")
const SVM_RANK_LEARN = joinpath(DEPOT_PATH[1], "mlsuite/bin/svm_rank_learn")
const SVM_RANK_CLASSIFY = joinpath(DEPOT_PATH[1], "mlsuite/bin/svm_rank_classify")
const RANKLIB = joinpath(DEPOT_PATH[1], "mlsuite/bin/ranklib")
const LIGHTGBM = joinpath(DEPOT_PATH[1], "mlsuite/bin/lightgbm")

export accuracy_score, r2_score
accuracy_score(a...; ka...) = pyimport("sklearn.metrics").accuracy_score(a...; ka...)
r2_score(a...; ka...) = pyimport("sklearn.metrics").r2_score(a...; ka...)

include("util.jl")

include("base/vw.jl")
include("base/h2o.jl")
include("base/dai.jl")
include("base/gbm.jl")
include("base/lgbm.jl")

include("class/scikit.jl")

include("regress/grid.jl")
include("regress/scikit.jl")

include("rank/svm.jl")
include("rank/ranklib.jl")
include("rank/tf.jl")

function __init__()
    boostlib = joinpath(DEPOT_PATH[1], "mlsuite/boost/lib")
    ENV["LD_LIBRARY_PATH"] = boostlib * ":" * get(ENV, "LD_LIBRARY_PATH", "")
    @require Elemental="902c3f28-d1ec-5e7e-8399-a24c3845ee38" include("regress/el.jl")
end

end


text:
"""
$(TYPEDEF)
Abstract parent type for all distributions.
"""
abstract type Distribution end

"""
$(TYPEDEF)
One dimensional normal distribution.
"""
Base.@kwdef mutable struct Normal <: Distribution
  mu::Float64
  sigma::Float64
  _type::String = "IncrementalInference.PackedNormal"
end
Normal(mu::Number, sigma::Number) = Normal(;mu, sigma)

"""
$(TYPEDEF)
One dimensional Rayleigh distribution.
"""
Base.@kwdef mutable struct Rayleigh <: Distribution
  sigma::Float64
  _type::String = "IncrementalInference.PackedRayleigh"
end
Rayleigh(sigma::Number) = Rayleigh(sigma=sigma)

"""
$(TYPEDEF)
Multidimensional normal distribution specified by means and a covariance matrix.
"""
Base.@kwdef mutable struct FullNormal <: Distribution
  mu::Vector{Float64}
  cov::Vector{Float64}
  _type::String = "IncrementalInference.PackedFullNormal"
end
# TODO: Generalize this to any number type.
FullNormal(mu::Vector{Float64}, cov::Matrix{Float64}) = FullNormal(;mu, cov=vec(cov))

"""
$(TYPEDEF)
One dimensional uniform distribution.
"""
Base.@kwdef mutable struct Uniform <: Distribution
  a::Float64
  b::Float64
  _type::String = "IncrementalInference.PackedUniform"
end
Uniform(a::Number, b::Number) = Uniform(;a, b)

"""
$(TYPEDEF)
Categorical distribution specified by a set of probabilities summing up to 1.
"""
Base.@kwdef mutable struct Categorical <: Distribution
  p::Vector{Float64}
  _type::String = "IncrementalInference.PackedCategorical"
end
# TODO: Generalize this to any number type.
Categorical(p::Vector{Float64}) = Categorical(;p)

Base.@kwdef mutable struct ManifoldKernelDensity <: Distribution
  varType::String
  pts::Vector{Vector{Float64}}
  bw::Vector{Float64}   = Float64[]
  partial::Vector{Int}  = Int[]
  infoPerCoord::Vector{Float64} = zeros(length(pts[1]))
  _type::String         = "IncrementalInference.PackedManifoldKernelDensity"
end
ManifoldKernelDensity(varType::String, pts::AbstractVector{<:AbstractVector{<:Real}}; kw...) = ManifoldKernelDensity(;varType,pts,kw...)


text:
using Plots

dt = 0.001
T = 20.0
vx_0 = -2
vy_0 = 1
x_0 = 1
y_0 = 1

q1_x = 5
q1_y = 4
q1 = 4
q2_x = -5
q2_y = -5
q2 = 10

function efield()
    t = 0
    x = x_0
    y = y_0
    vx = 0
    vy = 0

    xs = []
    ys = []

    while (t <= T)
        push!(xs, x)
        push!(ys, y)

        t += dt
        rx_1 = q1_x - x
        ry_1 = q1_y - y
        r_1 = sqrt(rx_1^2+ry_1^2)
        vx += ((rx_1*q1)/r_1^3)*dt
        vy += ((ry_1*q1)/r_1^3)*dt

        rx_2 = q2_x - x
        ry_2 = q2_y - y
        r_2 = sqrt(rx_2^2+ry_2^2)
        vx += ((rx_2*q2)/r_2^3)*dt
        vy += ((ry_2*q2)/r_2^3)*dt
        print(vx)

        x += vx * dt
        y += vy * dt

    end

    (xs, ys)
end

xs, ys = efield()

plot(xs, ys, label="", xlabel="x",ylabel="y", color=:cornflowerblue)
savefig("efield.png")


text:
# multi selection widget

SELECTEDORDERABLE    = 1 # whether selected items are orderable, selected always on top
SELECTSUBSTR         = 2 # search by substring (default by prefix)

defaultMultiSelectHelpText = utf8("""
arrows : move cursor
home   : move to top
end    : move to bottom
shft-up: move the item up
shft-dn: move the item down
space  : toggle selection
enter  : finalize selection

Search box:
alt-space: enter space inside search box
ctrl-a : move search cursor to start
ctrl-e : move search cursor to end
ctrl-k : empty search entry
ctrl-r : Toggle insertion/overwrite mode

ctrl-n : move to the next matched item
ctrl-p : move to the previous matched item
""")

type TwMultiSelectData
    choices::Array{UTF8String,1}
    selected::Array{UTF8String,1}
    datalist::Array{Any, 1}
    maxchoicelength::Int
    searchbox::Any
    currentLine::Int
    currentLeft::Int
    currentTop::Int
    selectmode::Int
    helpText::UTF8String
    TwMultiSelectData( arr::Array{UTF8String,1}, selected::Array{UTF8String,1} ) = new( arr, selected, Any[], 0, nothing, 1, 1, 1, 0, utf8("") )
end
TwMultiSelectData{T<:AbstractString,T2<:AbstractString}( arr::Array{T,1}, selected::Array{T2,1} ) = TwMultiSelectData( map( x->utf8( x ), arr ), map( x->utf8(x), selected ) )

# the ways to use it:
# standalone panel
# as a subwin as part of another widget (see next function)
# w include title width, if it's shown on the left
function newTwMultiSelect{T<:AbstractString}( scr::TwObj, arr::Array{T,1};
        posy::Any = :center,posx::Any = :center,
        selected = UTF8String[],
        title = utf8(""), maxwidth = 50, maxheight = 20, minwidth = 25,
        orderable = false, substrsearch=false )
    obj = TwObj( TwMultiSelectData( arr, UTF8String[ utf8(string(_)) for _ in selected ] ), Val{ :MultiSelect } )
    obj.box = true
    obj.title = title
    obj.borderSizeV= 1
    obj.borderSizeH= 1
    if  orderable
        obj.data.selectmode |= SELECTEDORDERABLE
    end
    if substrsearch
        obj.data.selectmode |= SELECTSUBSTR
    end
    rebuild_select_datalist( obj )
    obj.data.helpText = defaultMultiSelectHelpText
    obj.data.maxchoicelength = 0
    if !isempty(arr)
        obj.data.maxchoicelength = maximum( map(_->length(_), arr ) )
    end

    h = 2 + min( length( arr ), maxheight )
    w = 4 + max( min( max( length( title ), obj.data.maxchoicelength ), maxwidth ), minwidth )

    link_parent_child( scr, obj, h,w, posy, posx )

    obj.data.searchbox = newTwEntry( obj, UTF8String, width=minwidth, posy=:bottom, posx=1, box=false )
    obj.data.searchbox.title = utf8("?")
    obj.data.searchbox.hasFocus = false
    obj
end

function rebuild_select_datalist( o::TwObj{TwMultiSelectData} )
    o.data.datalist = Any[]
    if o.data.selectmode & SELECTEDORDERABLE != 0
        for s in o.data.selected
            push!( o.data.datalist, [ s, true ] )
        end
        for s in o.data.choices
            if !in( s, o.data.selected )
                push!( o.data.datalist, [s, false ] )
            end
        end
    else
        for s in o.data.choices
            push!( o.data.datalist, [s, in( s, o.data.selected ) ] )
        end
    end
end

function draw( o::TwObj{TwMultiSelectData} )
    werase( o.window )
    if o.box
        box( o.window, 0,0 )
    end
    if !isempty( o.title ) && o.box
        mvwprintw( o.window, 0, (@compat round(Int, ( o.width - length(o.title) )/2 )), "%s", o.title )
    end
    starty = o.borderSizeV
    viewContentHeight = o.height - o.borderSizeV * 2
    viewContentWidth  = o.width - o.borderSizeH * 2
    n = length( o.data.datalist )
    for r in o.data.currentTop:min( o.data.currentTop + viewContentHeight-1, n )
        flag = 0
        if r == o.data.currentLine
            flag = A_BOLD | COLOR_PAIR(o.hasFocus ? 15 : 30 )
        end
        s = o.data.datalist[r][1]
        if o.data.datalist[r][2]
            s = string( '\U2612' ) * " " * s
        else
            s = string( '\U2610' ) * " " * s
        end
        s = substr_by_width( s, o.data.currentLeft-1, viewContentWidth )

        wattron( o.window, flag )
        mvwprintw( o.window, r - o.data.currentTop + starty, o.borderSizeH, "%s", s )
        wattroff( o.window, flag )
    end
    draw( o.data.searchbox )
end

function select_search_next( o::TwObj{TwMultiSelectData}, step::Int, trivialstop::Bool )
    st = o.data.currentLine
    tmpstr = lowercase(o.data.searchbox.data.inputText)
    if length(tmpstr) == 0
        TermWin.beep()
        return 0
    end

    n = length( o.data.datalist )

    local i::Int = trivialstop ? st : mod1( st+step,n )
    local usesubstr::Bool = o.data.selectmode * SELECTSUBSTR != 0
    while true
        if usesubstr
            if contains( lowercase( o.data.datalist[i][1] ), tmpstr )
                o.data.currentLine = i
                return i
            end
        else
            if startswith( lowercase( o.data.datalist[i][1] ), tmpstr )
                o.data.currentLine = i
                return i
            end
        end
        i = mod1(i+step,n )
        if i == st
            TermWin.beep()
            return 0
        end
    end
end

function inject( o::TwObj{TwMultiSelectData}, token )
    @lintpragma( "Ignore incompatible type comparison")
    dorefresh = false
    retcode = :got_it # default behavior is that we know what to do with it

    viewContentWidth = o.width - o.borderSizeH * 2
    viewContentHeight = o.height - 2 * o.borderSizeV

    checkTop = () -> begin
        if o.data.currentTop > o.data.currentLine
            o.data.currentTop = o.data.currentLine
        elseif o.data.currentLine - o.data.currentTop > viewContentHeight - 1
            o.data.currentTop = o.data.currentLine - viewContentHeight + 1
        end
    end
    moveby = n -> begin
        sz = length( o.data.datalist )

        o.data.currentLine = max(1, min( sz, o.data.currentLine + n) )
        checkTop()
        return true
    end

    if token != :F1 && token != " "
        inputText = o.data.searchbox.data.inputText
        if token == " " # no break space
            token = " "
        end
        result = inject( o.data.searchbox, token )

        if result == :got_it
            if inputText != o.data.searchbox.data.inputText
                select_search_next( o, 1, true )
                checkTop()
            end
            refresh( o )
            return result
        end
    end

    if token == :esc
        retcode = :exit_nothing
    elseif token == :up
        dorefresh = moveby(-1)
    elseif token == :down
        dorefresh = moveby(1)
    elseif token == :left
        if o.data.currentLeft > 1
            o.data.currentLeft -= 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :right
        if o.data.currentLeft + viewContentWidth < o.data.maxchoicelength
            o.data.currentLeft += 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :shift_left # TODO ctrl-left
        if o.data.currentLeft > 1
            o.data.currentLeft -= 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :ctrlshift_left # TODO ctrl-left
        if o.data.currentLeft > 1
            o.data.currentLeft = 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :shift_right
        if o.data.currentLeft + viewContentWidth < o.data.maxchoicelength
            o.data.currentLeft += 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :ctrlshift_right
        if o.data.currentLeft + viewContentWidth < o.data.maxchoicelength
            o.data.currentLeft = o.data.maxchoicelength - viewContentWidth
            dorefresh = true
        else
            beep()
        end
    elseif token == :pageup
        dorefresh = moveby( -viewContentHeight )
    elseif token == :pagedown
        dorefresh = moveby( viewContentHeight )
    elseif token == :ctrl_n
        select_search_next( o, 1, false )
        checkTop()
        dorefresh = true
    elseif token == :ctrl_p
        select_search_next( o, -1, false )
        checkTop()
        dorefresh = true
    elseif token == :KEY_MOUSE
        (mstate,x,y, bs ) = getmouse()
        if mstate == :scroll_up
            dorefresh = moveby( -(@compat round(Int, viewContentHeight/10 )) )
        elseif mstate == :scroll_down
            dorefresh = moveby( (@compat round(Int, viewContentHeight/10 )) )
        elseif mstate == :button1_pressed && o.data.trackLine
            (rely,relx) = screen_to_relative( o.window, y, x )
            if 0<=relx<o.width && 0<=rely<o.height
                o.data.currentLine = o.data.currentTop + rely - o.borderSizeH + 1
                dorefresh = true
            else
                retcode = :pass
            end
        end
    elseif  token == :home
        if o.data.currentTop != 1 || o.data.currentLeft != 1 || o.data.currentLine != 1
            o.data.currentTop = 1
            o.data.currentLeft = 1
            o.data.currentLine = 1
            dorefresh = true
        else
            beep()
        end
    elseif in( token, Any[ Symbol("end") ] )
        n = length( o.data.datalist )
        if o.data.currentLine != n
            o.data.currentLine = n
            checkTop()
            dorefresh = true
        else
            beep()
        end
    elseif o.data.selectmode & SELECTEDORDERABLE != 0 && token == :shift_up
        currstatus = o.data.datalist[ o.data.currentLine][2]
        if currstatus && o.data.currentLine > 1
            currstr = o.data.datalist[o.data.currentLine][1]
            idx = findfirst( o.data.selected, currstr )
            o.data.selected[ idx-1], o.data.selected[ idx ] = (o.data.selected[ idx ], o.data.selected[ idx - 1] )
            o.data.currentLine -=1
            rebuild_select_datalist( o )
            checkTop()
            dorefresh = true
        else
            beep()
        end
    elseif o.data.selectmode & SELECTEDORDERABLE != 0 && token == :shift_down
        currstatus = o.data.datalist[ o.data.currentLine][2]
        if currstatus && o.data.currentLine < length( o.data.selected )
            currstr = o.data.datalist[o.data.currentLine][1]
            idx = findfirst( o.data.selected, currstr )
            o.data.selected[ idx+1], o.data.selected[ idx ] = (o.data.selected[ idx ], o.data.selected[ idx + 1] )
            o.data.currentLine +=1
            rebuild_select_datalist( o )
            checkTop()
            dorefresh = true
        else
            beep()
        end
    elseif token == " "
        currstr = o.data.datalist[o.data.currentLine][1]
        currstatus = o.data.datalist[ o.data.currentLine][2]
        if !currstatus # we are selecting it
            push!( o.data.selected, currstr )
            if o.data.selectmode & SELECTEDORDERABLE != 0
                rebuild_select_datalist( o )
            else
                o.data.datalist[o.data.currentLine][2] = true
            end
        else # we are de-selecting it
            idx = findfirst( o.data.selected, currstr )
            deleteat!( o.data.selected, idx )
            if o.data.selectmode & SELECTEDORDERABLE != 0
                rebuild_select_datalist( o )
            else
                o.data.datalist[o.data.currentLine][2] = false
            end
        end
        dorefresh = true
    elseif token == :enter || token == Symbol( "return" )
        o.value = o.data.selected
        retcode = :exit_ok
    else
        retcode = :pass # I don't know what to do with it
    end

    if dorefresh
        refresh(o)
    end

    return retcode
end

function helptext( o::TwObj{TwMultiSelectData} )
    o.data.helpText
end


text:
using NewtonsMethod
using Test

#Test method 1
f(x)=x^2-4.0
f′(x)=2x
@test newtonroot(f,f′,x₀=0.2)[1]≈2.0 atol=0.000001

f(x)=log(x)-20
f′(x)=1/x
@test newtonroot(f,f′,x₀=0.2)[1]≈4.851651954097909e8 atol=0.000001

f(x)=3x^2-5x+1
f′(x)=6x-5
@test newtonroot(f,f′,x₀=0.2)[1]≈0.23240824122077 atol=0.000001

f(x)=-10x^3-5x^2+20
f′(x)=-20x^2-10x
@test newtonroot(f,f′,x₀=0.2)[1]≈1.11338619006481 atol=0.000001

#Test method 2
f(x)=x^2-4.0
@test newtonroot(f,x₀=0.2)[1]≈2.0 atol=0.000001

f(x)=log(x)-20
@test newtonroot(f,x₀=0.2)[1]≈4.851651954097909e8 atol=0.000001

f(x)=3x^2-5x+1
@test newtonroot(f,x₀=0.2)[1]≈0.23240824122077 atol=0.000001

f(x)=-10x^3-5x^2+20
@test newtonroot(f,x₀=0.2)[1]≈1.11338619006481 atol=0.000001

#Test Big Float (Not really sure what to do here. What I did is compare the root with a BigFloat)
f(x)=3x^2-5x+1
a=BigFloat(0.23240824122077)
@testset "BigFloat" begin
  @test newtonroot(f,x₀=0.2)[1]≈0.23240824122077 atol=0.000001
  @test newtonroot(f,x₀=0.2)[1]≈a atol=0.000001
 end;

#(I also set all parameters to be BigFloat)
f(x)=3x^2-5x+1
@test newtonroot(f,x₀= BigFloat(0.2),tolerance = BigFloat(1E-7), maxiter = BigFloat(1000))[1]≈0.23240824122077 atol=0.000001

#Test non-convergence (return nothing)
f(x)=2+x^2
@test newtonroot(f,x₀=0.2)==nothing

#Test maxiter (first should return the value, second should return nothing)
f(x)=log(x)-20
a=newtonroot(f,x₀=0.2)[1] #Algorithm needs 17 iterations in this case
b=newtonroot(f,x₀=0.2,maxiter=5)
@testset "maxiter" begin
  @test a≈4.851651954097909e8 atol=0.000001
  @test b==nothing
 end;

#Test tolerance (As tolerance parameter increases, the code is less accurate)
f(x)=3x^2-5x+1
a=newtonroot(f,x₀=1)[1]
b=newtonroot(f,x₀=1, tolerance=0.0005)[1]
c=newtonroot(f,x₀=1, tolerance=0.05)[1]
@test f(a)<f(b)<f(c)


text:
# Counts the number of local maxima.
function count_maxima(histogram::AbstractArray)
    maxima_count = 0
    valid_range = firstindex(histogram)+1:lastindex(histogram)-1
    for i in valid_range
        if  histogram[i-1] < histogram[i] > histogram[i+1]
            maxima_count += 1
        end
    end
    return maxima_count
end

function find_maxima_indices(histogram::AbstractArray)
    indices = Array{Int, 1}()
    t = 0
    valid_range = firstindex(histogram)+1:lastindex(histogram)-1
    for i in valid_range
        if histogram[i-1] < histogram[i] > histogram[i+1]
            push!(indices, i)
        end
    end
    return indices
end

function smooth_histogram(histogram::AbstractArray, max_iterations::Int)
    histogram_local = convert.(Float32, histogram)
    maxima_count = count_maxima(histogram)
    smooth_histogram = similar(histogram_local)
    lb = firstindex(histogram_local)
    ub = lastindex(histogram_local)
    iterations = 0
    # Smooth histogram until at most two peaks remain or max_iterations is reached.
    while maxima_count > 2 && iterations < max_iterations
        smooth_histogram[lb] = (2*histogram_local[lb] + histogram_local[lb+1]) / 3
        smooth_histogram[ub] = (2*histogram_local[ub] + histogram_local[ub-1]) / 3
        for i in lb+1:ub-1
            smooth_histogram[i] = (histogram_local[i-1] + histogram_local[i] + histogram_local[i+1]) / 3
        end
        maxima_count = count_maxima(smooth_histogram)
        copyto!(histogram_local, smooth_histogram)
        iterations += 1
    end
    return histogram_local
end


text:
"""
    DockerRunner

Use `docker` as an execution engine; a reasonable backup for platforms that do
not have user namespaces (e.g. MacOS, Windows).
"""
mutable struct DockerRunner <: Runner
    docker_cmd::Cmd
    env::Dict{String, String}
    platform::AbstractPlatform

    shards::Vector{CompilerShard}
    workspace_root::String
end

docker_image(version::VersionNumber) = "julia_binarybuilder_rootfs:v$(version)"
function docker_image(rootfs::CompilerShard)
    name = artifact_name(rootfs)
    artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
    hash = artifact_hash(name, artifacts_toml; platform=rootfs.host)
    return string(
        "julia_binarybuilder_rootfs:",
        "v",
        string(rootfs.version),
        "-",
        string(bytes2hex(hash.bytes))[end-5:end],
    )
end

"""
    import_docker_image(rootfs::CompilerShard; verbose::Bool = false)

Checks to see if the given rootfs has been imported into docker yet; if it
hasn't, then do so so that we can run things like:

    docker run -ti binarybuilder_rootfs:v2018.08.27 /bin/bash

Which, after all, is the foundation upon which this whole doodad is built.
"""
function import_docker_image(rootfs::CompilerShard, workspace_root::String; verbose::Bool = false)
    # Does this image already exist?  If so, we're done!
    if success(`docker inspect --type=image $(docker_image(rootfs))`)
        if verbose
            @info("Docker base image already exists, skipping import...")
        end
        return
    end

    # Otherwise, import it!
    dockerfile_cmds = "ENTRYPOINT [\"/docker_entrypoint.sh\"]"
    rootfs_path = mount(rootfs, workspace_root; verbose=verbose)
    if verbose
        @info("Importing docker base image from $(rootfs_path) to $(docker_image(rootfs))")
    end
    try
        run(pipeline(pipeline(
            `tar -c -C $(rootfs_path) .`,
            `docker import - -c $(dockerfile_cmds) $(docker_image(rootfs))`;
        ); stdout=devnull))
    finally
        unmount(rootfs, workspace_root)
    end
    return
end

# Helper function to delete a previously-imported docker image
delete_docker_image() = delete_docker_image(first(choose_shards(HostPlatform())))
delete_docker_image(rootfs::CompilerShard) = success(`docker rmi -f $(docker_image(rootfs))`)

function DockerRunner(workspace_root::String;
                      cwd = nothing,
                      platform::AbstractPlatform = HostPlatform(),
                      workspaces::Vector = [],
                      extra_env=Dict{String, String}(),
                      verbose::Bool = false,
                      compiler_wrapper_path::String = mktempdir(),
                      toolchains_path::String = mktempdir(),
                      src_name::AbstractString = "",
                      shards = nothing,
                      kwargs...)
    global use_ccache

    # Check to make sure we're not going to try and bindmount within an
    # encrypted directory, as that can trigger kernel bugs
    check_encryption(workspace_root; verbose=verbose)

    # Construct environment variables we'll use from here on out
    platform = get_concrete_platform(platform; extract_kwargs(kwargs, (:preferred_gcc_version,:preferred_llvm_version,:compilers))...)
    envs = merge(platform_envs(platform, src_name; verbose=verbose), extra_env)

    # JIT out some compiler wrappers, add it to our mounts
    generate_compiler_wrappers!(platform; bin_path=compiler_wrapper_path, extract_kwargs(kwargs, (:compilers,:allow_unsafe_flags,:lock_microarchitecture))...)
    push!(workspaces, compiler_wrapper_path => "/opt/bin")

    # Generate CMake and Meson files
    generate_toolchain_files!(platform; toolchains_path=toolchains_path)
    push!(workspaces, toolchains_path => "/opt/toolchains")

    # the workspace_root is always a workspace, and we always mount it first
    insert!(workspaces, 1, workspace_root => "/workspace")

    # If we're enabling ccache, then map in a named docker volume for it
    if use_ccache
        if !isdir(ccache_dir())
            mkpath(ccache_dir())
        end
        push!(workspaces, "binarybuilder_ccache" => "/root/.ccache")
    end

    if isnothing(shards)
        # Choose the shards we're going to mount
        shards = choose_shards(platform; extract_kwargs(kwargs, (:preferred_gcc_version,:preferred_llvm_version,:bootstrap_list,:compilers))...)
    end

    # Import docker image
    import_docker_image(shards[1], workspace_root; verbose=verbose)

    # Construct docker command
    docker_cmd = `docker run --rm --privileged --user $(getuid):$(getgid)`#--cap-add SYS_ADMIN`

    if cwd != nothing
        docker_cmd = `$docker_cmd -w /$(abspath(cwd))`
    end

    # Add in read-only mappings and read-write workspaces
    for shard in shards[2:end]
        outside = mount_path(shard, workspace_root)
        inside = map_target(shard)
        docker_cmd = `$docker_cmd -v $(outside):$(inside):ro`
    end
    for (outside, inside) in workspaces
        if isdir(outside) || isfile(outside)
            outside = realpath(outside)
        end
        docker_cmd = `$docker_cmd -v $(outside):$inside`
    end

    # Build up environment mappings
    for (k, v) in envs
        docker_cmd = `$docker_cmd -e $k=$v`
    end

    # Finally, return the DockerRunner in all its glory
    return DockerRunner(docker_cmd, envs, platform, shards, workspace_root)
end

"""
    chown_cleanup(dr::DockerRunner)

On Linux, the user id inside of the docker container doesn't correspond to ours
on the outside, so permissions get all kinds of screwed up.  To fix this, we
have to `chown -R \$(id -u):\$(id -g) \$prefix`, which really sucks, but is
still better than nothing.  This is why we prefer the UserNSRunner on Linux.
"""
function chown_cleanup(dr::DockerRunner; verbose::Bool = false)
    # if !Sys.islinux()
    #     return
    # end

    # if verbose
    #     @info("chown'ing prefix back to us...")
    # end
    # run(`$(sudo_cmd()) chown $(getuid()):$(getgid()) -R $(dr.workspace_root)`)
end

function Base.run(dr::DockerRunner, cmd, logger::IO=stdout; verbose::Bool=false, tee_stream=stdout)
    did_succeed = true
    docker_cmd = `$(dr.docker_cmd) $(docker_image(dr.shards[1])) $(cmd)`
    @debug("About to run: $(docker_cmd)")

    try
        mount_shards(dr; verbose=verbose)
        oc = OutputCollector(docker_cmd; verbose=verbose, tee_stream=tee_stream)
        did_succeed = wait(oc)

        # First write out the actual command, then the command output
        println(logger, cmd)
        print(logger, merge(oc))
    finally
        unmount_shards(dr; verbose=verbose)

        # Cleanup permissions, if we need to.
        chown_cleanup(dr; verbose=verbose)
    end

    # Return whether we succeeded or not
    return did_succeed
end

function Base.read(dr::DockerRunner, cmd; verbose=false)
    did_succeed = true
    docker_cmd = `$(dr.docker_cmd) $(docker_image(dr.shards[1])) $(cmd)`
    @debug("About to run: $(docker_cmd)")

    local oc
    did_succeed = false
    try
        mount_shards(dr; verbose=verbose)
        oc = OutputCollector(docker_cmd; verbose=verbose)
        did_succeed = wait(oc)
    finally
        unmount_shards(dr; verbose=verbose)

        # Cleanup permissions, if we need to.
        chown_cleanup(dr; verbose=verbose)
    end

    if !did_succeed
        print(stderr, collect_stderr(oc))
        return nothing
    end

    return collect_stdout(oc)
end

function run_interactive(dr::DockerRunner, cmd::Cmd; stdin = nothing, stdout = nothing, stderr = nothing, verbose::Bool = false)
    function tty_or_nothing((s, bs))
        # If the `s` stream is `nothing`, we need to check whether the
        # corresponding Base stream `bs` is a TTY.
        s === nothing && (s = bs)
        return typeof(s) <: Base.TTY
    end
    run_flags = all(tty_or_nothing.(((stdin, Base.stdin), (stdout, Base.stdout), (stderr, Base.stderr)))) ? "-ti" : "-i"
    docker_cmd = `$(dr.docker_cmd) $(run_flags) -i $(docker_image(dr.shards[1])) $(cmd.exec)`
    if cmd.ignorestatus
        docker_cmd = ignorestatus(docker_cmd)
    end
    @debug("About to run: $(docker_cmd)")

    if stdin isa AnyRedirectable
        docker_cmd = pipeline(docker_cmd, stdin=stdin)
    end
    if stdout isa AnyRedirectable
        docker_cmd = pipeline(docker_cmd, stdout=stdout)
    end
    if stderr isa AnyRedirectable
        docker_cmd = pipeline(docker_cmd, stderr=stderr)
    end

    try
        mount_shards(dr; verbose=verbose)
        if stdout isa IOBuffer
            if !(stdin isa IOBuffer)
                stdin = devnull
            end
            process = open(docker_cmd, "r", stdin)
            @async begin
                while !eof(process)
                    write(stdout, read(process))
                end
            end
            wait(process)
            return success(process)
        else
            return success(run(docker_cmd))
        end
    finally
        unmount_shards(dr; verbose=verbose)
        # Cleanup permissions, if we need to.
        chown_cleanup(dr; verbose=verbose)
    end
end


text:
using MAT
using Random
include("../src/hypergraph-affinity-functions.jl")
D = matread("DBLP-Author-Gender_H.mat")
H = D["H"]
classes = D["class"]
order = vec(sum(H,dims = 2))
r = 21
alpha, B1, B2, H1, H2, R1, R2, N = Hypergraph_to_Scores(H,classes,r)
class1 = "female"
class2 = "male"

# homophily indices for clique expanded hypergraph

hw, hm = clique_expansion_homophily(N,10)

## Some plots
using Plots

lw = 2.5
ms = 9
gfs = 15
tfs = 15
titlesize = 14
s1 = 350
s2 = 300
k = 4
color1 = :green
color2 = :blue
title = "$k-author papers"
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity / Baseline "
p = Plots.plot([1; k],[1; 1],fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = 1:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

plot!(p,1:k,R2[k,1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
plot!(p,1:k,R1[k,1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)

savefig("Figures/Affinities_$(k)_DBLP.pdf")

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
"""
Reproduction of the empirical results concerning an empirical application
in the prediction of transportation mode choice using the best subset
maximum score approach of Chen and Lee (2017).
The work-trip mode choice dataset of Horowitz (1993) is included in the package.
"""
module transportationMode

using DataFrames, Compat
using MathProgBase, Gurobi
using ..BSBP

""" Best subset maximum score approach on the work-trip mode choice dataset of Horowitz.
Args:
    warm_start (int,optional):  Set to 1 for warm start strategy. Default = 1.
    tau (float,optional):       Tuning paramater to construct the refined bound used in the warm start approach. Default = 1.5.
    mio (int,optional):         Set to 1 for using Method 1 for the MIO formulation. Set to 2 for Method 2. Default to 1.
    q (int, optional):          Value of the variable selection bound. Default = 1.
    series_exp (int, optional): Set to 1 to use quadratic expansion terms as covariates. Default = 0.
    b (int, optional):          Bound value. Default = 1.
"""
function TransportationMode(;warm_start=0,tau=1.5,mio=2,series_exp=0,beta0=1,b=10,q=1,time_limit=86400)
  # load data
  data = readdlm("./data_horowitz.csv",',')

  bhat = []
  score = 0
  gap = 0
  rtime = 0
  ncount = 0

  println("estimation based on full sample")
  Y_tr=data[:,1]
  temp=data[:,2:end]

  n_tr=length(Y_tr)
  # [DCOST CARS DOVTT DIVTT]
  x_std=(temp.-repmat(mean(temp,1),n_tr,1))./repmat(std(temp,1),n_tr,1)
  x_foc = [x_std[:,1] ones(n_tr,1)] # [DCOST Intercept]

  if series_exp == 1
      z2 = x_std[:,2]
      z3 = x_std[:,3]
      z4 = x_std[:,4]
      x_aux1 = [z2, z3, z4] # linear terms
      x_aux2 = [z2.*z3, z3.*z4, z2.*z4]
      x_aux3 = [z2.*z2,z3.*z3, z4.*z4]
      x_aux = [x_aux1, x_aux2, x_aux3]
  else
      x_aux = x_std[:,2:4] # [CARS DOVTT DIVTT]
  end

  k=size(x_foc)[2]
  d=size(x_aux)[2]

  bnd= hcat(-b*ones(k-1+d,1),b*ones(k-1+d,1)) # set the initial parameter bounds

  tol = floor(sqrt(log(n_tr)*n_tr)/2) # set the tolerance level value
  println("tolerance level: ", tol/n_tr)

  if warm_start == 1 # warm start MIO
      bhat,score,gap,rtime,ncount  = WarmStartMaxScore(Y_tr,x_foc,x_aux,beta0,q,time_limit,tol,bnd,mio,tau)
  else # cold start MIO
      bhat,score,gap,rtime,ncount  = MaxScore(Y_tr,x_foc,x_aux,beta0,q,time_limit,tol,bnd,mio)
  end

  println("parameter estimates: ", bhat)
  println("avg_score gap time node_count:", score, " , ", gap, " , ", rtime, " , ", ncount)
  return score,gap,rtime,ncount
end

export TransportationMode
end


text:
#=
operate_df:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-14
=#

module DfOperator
    # packages
    using DataFrames
    using Statistics

    # dataframe
    global df = DataFrame(
            A = 1:5,
            B = rand(5),
            C = ["A", "B", "C", "A", "B"])

    # method
    function get_first_last(n)
        println("Dataframe first $(n) columns = ", first(df, n))
        println("")
        println("Dataframe last $(n) columns = ", last(df, n))
        println("")
    end

    function get_info()
        println("Data size = ", size(df))
        println("Column Names = ", names(df))
        println("Data type = ", eltype.(eachcol(df)))
        println("Describe = ", describe(df))
        println("")
    end

    function extract()
        println("df.A = ", df.A)
        println("")
        println("df[:, [:A, :B]] = ", df[:, [:A, :B]])
        println("")
        df.D = df.A .* 40
        println("df.D = df.A .* 40: ", df)
        println("")
        println("df[2, :] ", df[2, :])
        println("")
        println("df[2:5, :] ", df[2:5, :])
        println("")
        println("df[df.A.>3, :] ", df[df.A.>3, :])
        println("")
        println("df[(df.A.>3).&(df.C.=='B'), :] ", df[(df.A.>3).&(df.C.=="B"), :])
        println("")
        df[(df.A.>3).&(df.C.=="B"), :C] .= "X"
        println("df[(df.A.>3).&(df.C.=='B'), :C] .= 'X' ", df)
        println("")
    end

    function union_join()
        df_2 = DataFrame(
            A = [10],
            B = [0.1],
            C = ["C"],
            D = [20])
        println("append!(df, df_2) = ", append!(df, df_2))
        println("")
        df_3 = DataFrame(C = ["A", "B"], E = ["abs", "xyz"])
        println("Inner join(df, df_3, on=:C) = ", join(df, df_3, on=:C))
        println("")
        println("Left join(df, df_3, on=:C, kind=:left) = ", join(df, df_3, on=:C, kind=:left))
        println("")
    end

    function aggregate_by()
        println("aggregate(df, :C, sum) = ", aggregate(df, :C, sum))
        println("")
        println("by(df, :C, x->mean(x[:A])) = ", by(df, :C, x->mean(x[:A])))
        println("")
        println("mean(df[:B]) = ", mean(df[:B]))
        println("")
    end

    function sort_up_down()
        println("sort!(df, :B) = ", sort!(df, :B))
        println("")
        println("sort!(df, :B, rev=true) = ", sort!(df, :B, rev=true))
        println("")
    end

    function change_colname()
        new_col_list = [:W, :X, :Y, :Z]
        println("Change Col Name by List = ", rename!(df, new_col_list))
        println("")
        new_col_dict = Dict(:W=>:A, :X=>:B, :Y=>:C, :Z=>:D)
        println("Change Col Name by Dict = ", rename!(df, new_col_dict))
        println("")
    end

    function v_h_cat()
        df_4 = DataFrame(A = [1, 2, 3], B = [4, 5, 6])
        df_5 = DataFrame(A = [11, 12, 13], B = [14, 15, 16])
        println("vcat(df_4, df_5) = ", vcat(df_4, df_5))
        println("")
        df_6 = DataFrame(A = [1, 2, 3], B = [4, 5, 6])
        df_7 = DataFrame(C = [11, 12, 13], D = [14, 15, 16])
        println("hcat(df_6, df_7) = ", hcat(df_6, df_7))
        println("")
    end

    function main()
        get_first_last(3)
        get_info()
        extract()
        union_join()
        aggregate_by()
        sort_up_down()
        change_colname()
        v_h_cat()
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .DfOperator
    DfOperator.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
@testset "203.remove-linked-list-elements.jl" begin
    @test remove_elements!(ListNode{Int}([1, 2, 6, 3, 4, 5, 6]), 6) ==
        ListNode{Int}([1, 2, 3, 4, 5])
    @test remove_elements!(nothing, 1) === nothing
    @test remove_elements!(ListNode{Int}([7, 7, 7, 7]), 7) === nothing
end


text:
"""
    function video_trajdata_replay

Makes a video of the trajdata taking frame range as input

# Example
```julia
f = FilteringEnvironment()
video_trajdata_replay(f,id_list=[2,7,10,18,25],range=1:30,filename="media/lower_1.gif")
```
"""
function video_trajdata_replay(f::FilteringEnvironment;id_list = [],range=nothing,filename)
    frames = Frames(MIME("image/png"), fps=10)
    mr = MergingRoadway(f.roadway) # Wrapper for specialized render for merging lanes
    mp = VecE2(1064.5227, 959.1559)
    for i in range
        temp_scene = deepcopy(f.traj[i])
        if !isempty(id_list) keep_vehicle_subset!(temp_scene,id_list) end
        renderables = [
            mr,
            (FancyCar(car=temp_scene[j]) for j in 1:length(temp_scene))...,
            IDOverlay(scene=temp_scene),
            TextOverlay(text=["frame=$(i)"],font_size=12)
        ]
        scene_visual = render(renderables,camera=StaticCamera(position=mp,zoom=5.))

        push!(frames,scene_visual)
    end
    write(filename,frames)
    print("Making interaction replay video called $filename")

    return nothing
end

# function: make a video from a list of scenes
"""
    function scenelist2video(scene_list;filename = "media/scenelist_to_video.mp4")
- Make video from a list of scenes (generally generated using `get_hallucination_scenes`

# Examples
```julia
scenelist2video(f,scenelist,filename="media/test_scenelist2video.mp4")
```
"""
function scenelist2video(f::FilteringEnvironment,scene_list;filename)
    frames = Frames(MIME("image/png"),fps = 10)
    mr = MergingRoadway(f.roadway) # Wrapper for specialized render for merging lanes
    mp = VecE2(1064.5227, 959.1559)

    # Loop over list of scenes and convert to video
    for i in 1:length(scene_list)
        temp_scene = deepcopy(scene_list[i])
        renderables = [
            mr,
            (FancyCar(car=temp_scene[j]) for j in 1:length(temp_scene))...,
            IDOverlay(scene=temp_scene),
            TextOverlay(text=["frame=$(i)"],font_size=12)
        ]
        scene_visual = render(renderables,camera=StaticCamera(position=mp,zoom=5.))
        push!(frames,scene_visual)
    end
    print("scenelist2video says: Making video filename: $(filename)\n")
    write(filename,frames)
    return nothing
end

"""
function video_overlay_scenelists(scene_list_1,scene_list_2,roadway,filename)

# Example
```julia
scene_list_1 = run_vehicles(id_list=[6,8,19,28,29],roadway=road_ext,traj=traj_ext,
filename="model_driven.mp4")
scene_list_2 = f.traj[1:length(scene_list_1)]
video_overlay_scenes(scene_list_1,scene_list_2,id_list=[6,8,13,28,29],
roadway=road_ext,filename="imit_real.mp4")
```
"""
function video_overlay_scenelists(scene_list_1,scene_list_2;
    id_list=[],roadway,filename)
    @assert length(scene_list_1) == length(scene_list_2)
    frames = Frames(MIME("image/png"),fps = 10)
    mr = MergingRoadway(roadway) # Wrapper for specialized render for merging lanes
    mp = VecE2(1064.5227, 959.1559)

    # Loop over list of scenes and convert to video
    for i in 1:length(scene_list_1)
        if !isempty(id_list) keep_vehicle_subset!(scene_list_1[i],id_list) end
        if !isempty(id_list) keep_vehicle_subset!(scene_list_2[i],id_list) end
        scene1 = scene_list_1[i]
        scene2 = scene_list_2[i]

        # By default, all the cars are assigned random colors
        # Explicitly color code the cars to have only 2 colors
        renderables = [
            mr,
            (FancyCar(car=scene1[j],color=colorant"blue") for j in 1:length(scene1))...,
            IDOverlay(scene=scene1),
            (FancyCar(car=scene2[j],color=colorant"red") for j in 1:length(scene2))...,
            IDOverlay(scene=scene2),
            TextOverlay(text=["frame=$(i)"],font_size=12)
        ]
        scene_visual = render(renderables,camera=StaticCamera(position=mp,zoom=5.))
        push!(frames,scene_visual)
    end
    print("video_overlay_scenelists: List 1 in blue, list 2 in red. Filename: $(filename)\n")
    write(filename,frames)
    return nothing
end

# function: make a video from a list of scenes with curvepts overlayed
"""
    function scenelist2video(scene_list;filename = "media/scenelist_to_video.mp4")
- Make video from a list of scenes (generally generated using `get_hallucination_scenes`

# Examples
```julia
scene = Scene(500)
get!(scene,traj_interaction,1);
scene_list = get_hallucination_scenes(scene,models=models);
scenelist2video(scene_list,filename="media/scenelist_to_video.mp4")
```
"""
function scenelist2video_curvepts(scene_list;id_list=[],roadway,filename)
    frames = Frames(MIME("image/png"),fps = 10)

    # Loop over list of scenes and convert to video
    for i in 1:length(scene_list)
        if !isempty(id_list) keep_vehicle_subset!(scene_list[i],id_list) end
        scene_visual = render(scene_list[i],roadway,
        [IDOverlay(),TextOverlay(text=["frame=$(i)"],font_size=12),curvepts_overlay(roadway,colorant"yellow",0.05)],
        cam=FitToContentCamera(0.),
        #cam = SceneFollowCamera(10.)
        )
        push!(frames,scene_visual)
    end
    print("Making video filename: $(filename)\n")
    write(filename,frames)
    return nothing
end

function scenelist2video_mergeoverlay(scene_list;id_list=[],roadway,filename)
    frames = Frames(MIME("image/png"),fps = 10)
    env = MergingEnvironment()
    # Loop over list of scenes and convert to video
    for i in 1:length(scene_list)
        if !isempty(id_list) keep_vehicle_subset!(scene_list[i],id_list) end
        scene_visual = render(scene_list[i],roadway,
        [IDOverlay(),TextOverlay(text=["frame=$(i)"],font_size=12),MergeOverlay(env)],
        cam=FitToContentCamera(0.),
        #cam = SceneFollowCamera(10.)
        )
        push!(frames,scene_visual)
    end
    print("Making video filename: $(filename)\n")
    write(filename,frames)
    return nothing
end

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
module MultiplicationPartialBase

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export multiplication

function multiplication(model)

    #tot_cover_of_i += x_axis_covered*y_axis_covered

    @expression(model, static_new_variable_1[i=1:length(boxes), j=1:length(boxes); i!=j], 0.5*(model[:x_axis_covered][i, j] + model[:y_axis_covered][i, j]))
    @expression(model, static_new_variable_2[i=1:length(boxes), j=1:length(boxes); i!=j], 0.5*(model[:x_axis_covered][i, j] - model[:y_axis_covered][i, j]))

    # Piecewise linearisation start, need two, one for new each variable
    @variable(model, piecewise1[i=1:length(boxes), j=1:length(boxes); i!=j])
    @variable(model, static_lambda_box_cover1[i=1:length(boxes), j=1:length(boxes), k=1:10; i!=j] >= 0)
    @expression(model, static_lambda_box_cover1_row[i=1:length(boxes), j=1:length(boxes); i!=j],
        [
        model[:static_lambda_box_cover1][i, j, 1],
        model[:static_lambda_box_cover1][i, j, 2],
        model[:static_lambda_box_cover1][i, j, 3],
        model[:static_lambda_box_cover1][i, j, 4],
        model[:static_lambda_box_cover1][i, j, 5],
        model[:static_lambda_box_cover1][i, j, 6],
        model[:static_lambda_box_cover1][i, j, 7],
        model[:static_lambda_box_cover1][i, j, 8],
        model[:static_lambda_box_cover1][i, j, 9],
        model[:static_lambda_box_cover1][i, j, 10]
        ]
    )
    @constraint(model, static_lambda_box_cover_constraint1[i=1:length(boxes), j=1:length(boxes); i!=j],
        sum(model[:static_lambda_box_cover1_row][i, j]) == 1
    )
    @constraint(model, static_lambda_box_cover1_sos[i=1:length(boxes), j=1:length(boxes); i!=j],
        model[:static_lambda_box_cover1_row][i, j] in SOS2()
    )

    @constraint(model, static_product_lambda_constraint1[i=1:length(boxes), j=1:length(boxes); i!=j], (
        -model[:static_new_variable_1][i, j] +
        lin_plus[1]*model[:static_lambda_box_cover1][i, j, 2] +
        lin_plus[2]*model[:static_lambda_box_cover1][i, j, 3] +
        lin_plus[3]*model[:static_lambda_box_cover1][i, j, 4] +
        lin_plus[4]*model[:static_lambda_box_cover1][i, j, 5] +
        lin_plus[5]*model[:static_lambda_box_cover1][i, j, 6] +
        lin_plus[6]*model[:static_lambda_box_cover1][i, j, 7] +
        lin_plus[7]*model[:static_lambda_box_cover1][i, j, 8] +
        lin_plus[8]*model[:static_lambda_box_cover1][i, j, 9] +
        lin_plus[9]*model[:static_lambda_box_cover1][i, j, 10] == 0)
    )
    @constraint(model, static_piecewise_lambda_constraint1[i=1:length(boxes), j=1:length(boxes); i!=j], (
        -model[:piecewise1][i, j] +
        lin_plus_2[1]*model[:static_lambda_box_cover1][i, j, 2] +
        lin_plus_2[2]*model[:static_lambda_box_cover1][i, j, 3] +
        lin_plus_2[3]*model[:static_lambda_box_cover1][i, j, 4] +
        lin_plus_2[4]*model[:static_lambda_box_cover1][i, j, 5] +
        lin_plus_2[5]*model[:static_lambda_box_cover1][i, j, 6] +
        lin_plus_2[6]*model[:static_lambda_box_cover1][i, j, 7] +
        lin_plus_2[7]*model[:static_lambda_box_cover1][i, j, 8] +
        lin_plus_2[8]*model[:static_lambda_box_cover1][i, j, 9] +
        lin_plus_2[9]*model[:static_lambda_box_cover1][i, j, 10] == 0)
    )

    @variable(model, piecewise2[i=1:length(boxes), j=1:length(boxes); i!=j], )
    @variable(model, static_lambda_box_cover2[i=1:length(boxes), j=1:length(boxes), k=1:10; i!=j] >= 0)
    @expression(model, static_lambda_box_cover2_row[i=1:length(boxes), j=1:length(boxes); i!=j],
        [
        model[:static_lambda_box_cover2][i, j, 1],
        model[:static_lambda_box_cover2][i, j, 2],
        model[:static_lambda_box_cover2][i, j, 3],
        model[:static_lambda_box_cover2][i, j, 4],
        model[:static_lambda_box_cover2][i, j, 5],
        model[:static_lambda_box_cover2][i, j, 6],
        model[:static_lambda_box_cover2][i, j, 7],
        model[:static_lambda_box_cover2][i, j, 8],
        model[:static_lambda_box_cover2][i, j, 9],
        model[:static_lambda_box_cover2][i, j, 10]
        ]
    )
    @constraint(model, static_lambda_box_cover_constraint2[i=1:length(boxes), j=1:length(boxes); i!=j],
        sum(model[:static_lambda_box_cover2_row][i, j]) == 1
    )
    @constraint(model, static_lambda_box_cover2_sos[i=1:length(boxes), j=1:length(boxes); i!=j],
        model[:static_lambda_box_cover2_row][i, j] in SOS2()
    )

    @constraint(model, static_product_lambda_constraint2[i=1:length(boxes), j=1:length(boxes); i!=j],
        (
        -model[:static_new_variable_2][i, j] +
        lin_min[1]*model[:static_lambda_box_cover2][i, j, 2] +
        lin_min[2]*model[:static_lambda_box_cover2][i, j, 3] +
        lin_min[3]*model[:static_lambda_box_cover2][i, j, 4] +
        lin_min[4]*model[:static_lambda_box_cover2][i, j, 5] +
        lin_min[5]*model[:static_lambda_box_cover2][i, j, 6] +
        lin_min[6]*model[:static_lambda_box_cover2][i, j, 7] +
        lin_min[7]*model[:static_lambda_box_cover2][i, j, 8] +
        lin_min[8]*model[:static_lambda_box_cover2][i, j, 9] +
        lin_min[9]*model[:static_lambda_box_cover2][i, j, 10] == 0
        )
    )
    @constraint(model, static_piecewise_lambda_constraint2[i=1:length(boxes), j=1:length(boxes); i!=j],
        (
        -model[:piecewise2][i, j] +
        lin_min_2[1]*model[:static_lambda_box_cover2][i, j, 2] +
        lin_min_2[2]*model[:static_lambda_box_cover2][i, j, 3] +
        lin_min_2[3]*model[:static_lambda_box_cover2][i, j, 4] +
        lin_min_2[4]*model[:static_lambda_box_cover2][i, j, 5] +
        lin_min_2[5]*model[:static_lambda_box_cover2][i, j, 6] +
        lin_min_2[6]*model[:static_lambda_box_cover2][i, j, 7] +
        lin_min_2[7]*model[:static_lambda_box_cover2][i, j, 8] +
        lin_min_2[8]*model[:static_lambda_box_cover2][i, j, 9] +
        lin_min_2[9]*model[:static_lambda_box_cover2][i, j, 10] == 0
        )
    )

    # Finding total partial stability of one box start. (length(boxes)-1 as we do not want to get support from sentinel).
    # Multiply by if, as if tells whether i,j is to be considered.
    @expression(model, pbs_box_ij[i=1:length(boxes), j=1:length(boxes); i!=j],
        model[:if1_pbs][i, j] * (model[:piecewise1][i, j] - model[:piecewise2][i, j])
    )
    @expression(model, pbs_box_i[i=1:length(boxes)],
        sum([model[:pbs_box_ij][i, j] for j in 1:length(boxes) if i!=j])
    )

    return model
end

end


text:
export TimeScaleMODWT
import Wavelets
import Wavelets: wavelet, maxdyadiclevel, modwt

"""
    TimeScaleMODWT <: WaveletProbabilitiesEstimator
    TimeScaleMODWT(wl::Wavelets.WT.OrthoWaveletClass = Wavelets.WT.Daubechies{12}())

Apply the maximal overlap discrete wavelet transform (MODWT) to a
signal, then compute probabilities/entropy from the energies at different
wavelet scales. This implementation is based on Rosso et
al. (2001)[^Rosso2001].
Optionally specify a wavelet to be used.

The probability `p[i]` is the relative/total energy for the i-th wavelet scale.

## Example

Manually picking a wavelet is done as follows.

```julia
using Entropies, Wavelets
N = 200
a = 10
t = LinRange(0, 2*a*π, N)
x = sin.(t .+  cos.(t/0.1)) .- 0.1;

# Pick a wavelet (if no wavelet provided, defaults to Wavelets.WL.Daubechies{12}())
wl = Wavelets.WT.Daubechies{12}()

# Compute the probabilities (relative energies) at the different wavelet scales
probabilities(x, TimeScaleMODWT(wl))
```

[^Rosso2001]: Rosso, O. A., Blanco, S., Yordanova, J., Kolev, V., Figliola, A., Schürmann, M., & Başar, E. (2001). Wavelet entropy: a new tool for analysis of short duration brain electrical signals. Journal of neuroscience methods, 105(1), 65-75.
"""
struct TimeScaleMODWT <: WaveletProbabilitiesEstimator
    wl::Wavelets.WT.OrthoWaveletClass
    function TimeScaleMODWT(wl::Wavelets.WT.OrthoWaveletClass = Wavelets.WT.Daubechies{12}())
        new(wl)
    end
end

function get_modwt(x::AbstractVector{T}, wl::Wavelets.WT.OrthoWaveletClass = Wavelets.WT.Daubechies{12}()) where T<:Real
    orthofilter = wavelet(wl)
    nscales = maxdyadiclevel(x)
    W = modwt(x, orthofilter, nscales)
end

function energy_at_scale(W::AbstractArray{T, 2}, j::Int) where T<:Real
    1 <= j <= size(W, 2) || error("Scale j does not exist in wave coefficient matrix W. Available scales are j=1:$(size(W, 2))")
    Eⱼ = sum(W[:, j] .^ 2)
end

function energy_at_time(W::AbstractArray{T, 2}, t::Int) where T<:Real
    1 <= t <= size(W, 1) || error("Time t does not exist in wave coefficient matrix W. Available times are t=1:$(size(W, 1))")
    Eⱼ = sum(W[t, :] .^ 2)
end

function energy_total(W::AbstractArray{T, 2}) where T<:Real
    Etot = sum(W .^ 2)
end

function relative_wavelet_energy(W::AbstractArray{T, 2}, j::Int) where T<:Real
    energy_at_scale(W, j) / energy_total(W)
end

function relative_wavelet_energies(W::AbstractArray{T, 2}, js = 1:size(W, 2)) where T<:Real
    all(1 .<= js .<= size(W, 2)) || error(ArgumentError("scales $(js) contains scales not present in wavelet coefficient matrix with scales j=1:$(size(W, 2))"))
    [energy_at_scale(W, j) / energy_total(W) for j in js]
end

function time_scale_density(x::AbstractVector{T}, wl::Wavelets.WT.OrthoWaveletClass = WT.Daubechies{12}()) where T
    W = get_modwt(x, wl)
    Pⱼs = relative_wavelet_energies(W)
end

function probabilities(x::AbstractVector{T}, est::TimeScaleMODWT) where T<:Real
    Probabilities(time_scale_density(x, est.wl))
end


text:
export Generator, sampleUQ

# function lrelu(x, th=0.2)
#   return tf.maximum(th * x, x)
# end

# function Generator(x, isTrain=true; vmin=-1, vmax=1)
#   local o
#   variable_scope("generator") do

#       # 1st hidden layer
#       # conv1 = tf.layers.conv2d_transpose(x, 128, [4, 4], strides=(1, 1), padding="valid")
#       lrelu1 = lrelu(tf.layers.batch_normalization(conv1, training=isTrain), 0.2)
#       # lrelu1 = lrelu(conv1, 0.2)

#       # 2nd hidden layer
#       conv2 = tf.layers.conv2d_transpose(lrelu1, 128, [8, 8], strides=(2, 2), padding="same")
#       lrelu2 = lrelu(tf.layers.batch_normalization(conv2, training=isTrain), 0.2)
#       # lrelu2 = lrelu(conv2, 0.2)

#       # 3rd hidden layer
#       conv3 = tf.layers.conv2d_transpose(lrelu2, 128, [8, 8], strides=(2, 2), padding="same")
#       lrelu3 = lrelu(tf.layers.batch_normalization(conv3, training=isTrain), 0.2)
#       # lrelu3 = lrelu(conv3, 0.2)

#       # 4th hidden layer
#       conv4 = tf.layers.conv2d_transpose(lrelu3, 128, [8, 8], strides=(2, 2), padding="same")
#       lrelu4 = lrelu(tf.layers.batch_normalization(conv4, training=isTrain), 0.2)
#       # lrelu4 = lrelu(conv4, 0.2)

#       # output layer
#       conv5 = tf.layers.conv2d_transpose(lrelu4, 1, [8, 8], strides=(2, 2), padding="same")
#       # o = (vmax - vmin)/2. * (tf.nn.tanh(conv5) - (-1)) + vmin
#       o = conv5
#   end
#   return o
# end

function Generator(z, isTrain=true; num_layer=5, h0=4, w0=8, vmin=nothing, vmax=nothing)
  local o
  # activation = tf.keras.activations.relu
  # activation = tf.keras.activations.tanh
  activation = tf.keras.layers.LeakyReLU(alpha=0.1)
  variable_scope("generator") do
    x = tf.keras.layers.Dense(units = h0 * w0 * 8, use_bias=false)(z)
    x = tf.reshape(x, shape=[-1, w0, h0, 8])
    x = tf.keras.activations.tanh(x)
    # x = activation(x)

    # x = tf.keras.layers.Conv2DTranspose(64, (5,5), strides=(2, 2), padding="same", use_bias=false)(x)
    # x = tf.pad(x, ((0,0), (1,2), (1,2), (0,0)), "REFLECT")
    # x = tf.keras.layers.Conv2D(16, [4, 4], strides=(1, 1), padding="valid")(x)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(16, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
    # x = activation(x)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(32, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
    # x = activation(x)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(64, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
    # x = activation(x)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(32, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
    # x = activation(x)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(16, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
    # x = activation(x)

    for i = 1:num_layer
      x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
      # x = tf.keras.layers.Conv2D(min(2^(i+3), 2^(num_layer-i+4)), [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
      x = tf.keras.layers.Conv2D(2^(num_layer-i+4), [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
      x = activation(x)
    end

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(1, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)

    # x = tf.keras.layers.Conv2D(1, [1, 1], strides=(1, 1), padding="same", use_bias=false)(x)
    x = tf.keras.layers.Conv2D(1, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)

    if !isnothing(vmax) && !isnothing(vmin)
        o = ((vmax - vmin) * tf.keras.activations.tanh(x) + (vmax + vmin))/2.0
        @info "vmin=$vmin, vmax=$vmax"
    else
        o = x
    end
    o = tf.squeeze(o)
    # o = cast(o, Float64)
  end
  return o
end

function Generator(z, isTrain, dropout_rate; num_layer=5, h0=4, w0=8, vmin=nothing, vmax=nothing)
  local o
  # activation = tf.keras.activations.relu
  # activation = tf.keras.activations.tanh
  activation = tf.keras.layers.LeakyReLU(alpha=0.1)
  variable_scope("generator") do
    x = tf.keras.layers.Dense(units = h0 * w0 * 8, use_bias=false)(z)
    x = tf.reshape(x, shape=[-1, w0, h0, 8])
    x = tf.keras.activations.tanh(x)
    # x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)

    # x = activation(x)
    # x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)
    # x = tf.keras.layers.SpatialDropout2D(dropout_rate)(x, isTrain)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(32, [4, 4], strides=(1, 1), padding="same")(x)
    # x = activation(x)
    # x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)
    # x = tf.keras.layers.SpatialDropout2D(dropout_rate)(x, isTrain)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(64, [4, 4], strides=(1, 1), padding="same")(x)
    # x = activation(x)
    # x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)
    # x = tf.keras.layers.SpatialDropout2D(dropout_rate)(x, isTrain)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(32, [4, 4], strides=(1, 1), padding="same")(x)
    # x = activation(x)
    # x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)
    # x = tf.keras.layers.SpatialDropout2D(dropout_rate)(x, isTrain)

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(16, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
    # x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)
    # x = activation(x)

    for i = 1:num_layer
      x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
      # x = tf.keras.layers.Conv2D(min(2^(i+3), 2^(num_layer-i+4)), [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
      x = tf.keras.layers.Conv2D(2^(num_layer-i+4), [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)
      x = activation(x)
      x = tf.keras.layers.Dropout(dropout_rate)(x, isTrain)
    end

    # x = tf.keras.layers.UpSampling2D((2, 2), interpolation="bilinear")(x)
    # x = tf.keras.layers.Conv2D(1, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)

    # x = tf.keras.layers.Conv2D(1, [1, 1], strides=(1, 1), padding="same", use_bias=false)(x)
    x = tf.keras.layers.Conv2D(1, [4, 4], strides=(1, 1), padding="same", use_bias=false)(x)

    if !isnothing(vmax) && !isnothing(vmin)
        o = ((vmax - vmin) * tf.keras.activations.tanh(x) + (vmax + vmin))/2.0
        @info "vmin=$vmin, vmax=$vmax"
    else
        o = x
    end
    o = tf.squeeze(o)
    # o = cast(o, Float64)
  end
  return o
end

function sampleUQ(batch_size, sample_size, rcv_size; z_size=8, h0=4, w0=8, vmin=nothing, vmax=nothing)
  isTrain = placeholder(Bool)
  y = placeholder(Float64, shape=(sample_size, rcv_size...))
  z = placeholder(Float32, shape=(batch_size, z_size))
  G_z = cast(Float64, Generator(z, isTrain, h0=base, w0=ratio, vmin=vmin, vmax=vmax))
  return G_z, isTrain, y, z
end