text:
"""
    OffsetAxis(keys::AbstractUnitRange{<:Integer}, parent::AbstractUnitRange{<:Integer})
    OffsetAxis(offset::Integer, parent::AbstractUnitRange{<:Integer})

An axis that has the indexing behavior of an [`AbstractOffsetAxis`](@ref) and retains an
offset from its underlying indices in its keys. Note that `offset` is only the offset from
the parent indices. If `OffsetAxis` is part of an `AxisArray`, the number returned by
`ArrayInterface.offsets` refers to the offset from zero, not the offset found in this axis.

## Examples

Users may construct an `OffsetAxis` by providing an from a set of indices.
```jldoctest offset_axis_examples
julia> using AxisIndices

julia> axis = AxisIndices.OffsetAxis(-2, 1:3)
offset(-2)(SimpleAxis(1:3))

```

In this instance the first index of the wrapped indices is 1 (`firstindex(indices(axis))`)
but adding the offset (`-2`) moves it to `-1`.
```jldoctest offset_axis_examples
julia> firstindex(axis)
-1

julia> axis[-1]
-1
```

Similarly, the last index is move by `-2`.
```jldoctest offset_axis_examples
julia> lastindex(axis)
1

julia> axis[1]
1

```

This means that traditional one based indexing no longer applies and may result in
errors.
```jldoctest offset_axis_examples
julia> axis[3]
ERROR: BoundsError: attempt to access offset(-2)(SimpleAxis(1:3)) at index [3]
[...]
```

When an `OffsetAxis` is reconstructed the offset from indices are presserved.
```jldoctest offset_axis_examples
julia> axis[0:1]  # offset of -2 still applies
offset(-2)(SimpleAxis(2:3))

```
"""
struct OffsetAxis{I,Inds<:AbstractAxis,F} <: AbstractOffsetAxis{I,Inds,F}
    offset::F
    parent::Inds

    function OffsetAxis{I,Inds,F}(f::Integer, inds::AbstractUnitRange) where {I,Inds,F}
        if inds isa Inds && f isa F
            return new{I,Inds,F}(f, inds)
        else
            return OffsetAxis(f, Inds(inds))
        end
    end

    function OffsetAxis{I,Inds,F}(ks::AbstractUnitRange, inds::AbstractUnitRange) where {I,Inds,F}
        check_axis_length(ks, inds)
        return OffsetAxis{I,Inds,F}(static_first(ks) - static_first(inds), inds)
    end

    # OffsetAxis{I,Inds}
    function OffsetAxis{I,Inds}(f::Integer, inds::AbstractUnitRange) where {I,Inds}
        if inds isa Inds
            return OffsetAxis{I,Inds,typeof(f)}(f, inds)
        else
            return OffsetAxis{I,Inds}(f, Inds(inds))
        end
    end
    @inline function OffsetAxis{I,Inds}(ks::AbstractUnitRange, inds::AbstractUnitRange) where {I,Inds}
        check_axis_length(ks, inds)
        return OffsetAxis{I,Inds}(static_first(ks) - static_first(inds), inds)
    end
    @inline function OffsetAxis{I,Inds}(ks::AbstractUnitRange) where {I,Inds}
        f = static_first(ks)
        return OffsetAxis{I}(f - one(f), Inds(OneTo(static_length(ks))))
    end

    # OffsetAxis{I}
    function OffsetAxis{I}(f::Integer, inds::AbstractAxis) where {I}
        return OffsetAxis{I,typeof(inds)}(f, inds)
    end
    function OffsetAxis{I}(f::Integer, inds::AbstractArray) where {I}
        return OffsetAxis{I}(f, compose_axis(inds))
    end
    function OffsetAxis{I}(f::AbstractUnitRange, inds::AbstractArray) where {I}
        return OffsetAxis{I}(f, compose_axis(inds))
    end
    function OffsetAxis{I}(ks::AbstractUnitRange, inds::AbstractAxis) where {I}
        check_axis_length(ks, inds)
        return OffsetAxis{I}(static_first(ks) - static_first(inds), inds)
    end
    function OffsetAxis{I}(ks::AbstractUnitRange) where {I}
        f = static_first(ks)
        return OffsetAxis{I}(f - one(f), SimpleAxis(One():static_length(ks)))
    end
    function OffsetAxis{I}(ks::AbstractUnitRange, inds::AbstractOffsetAxis) where {I}
        check_axis_length(ks, inds)
        p = parent(inds)
        return OffsetAxis{I}(static_first(ks) + static_first(inds) - static_first(p), p)
    end
    function OffsetAxis{I}(f::Integer, inds::AbstractOffsetAxis) where {I}
        p = parent(inds)
        return OffsetAxis{I}(f + static_first(inds) - static_first(p), parent(inds))
    end

    # OffsetAxis
    OffsetAxis(f::Integer, inds::AbstractAxis) = OffsetAxis{eltype(inds)}(f, inds)
    OffsetAxis(f::AbstractUnitRange, inds::AbstractAxis) = OffsetAxis{eltype(inds)}(f, inds)
    OffsetAxis(f::Integer, inds::AbstractArray) = OffsetAxis(f, compose_axis(inds))
    OffsetAxis(ks::AbstractUnitRange, inds::AbstractArray) = OffsetAxis(ks, compose_axis(inds))
    function OffsetAxis(ks::Ks) where {Ks}
        fst = static_first(ks)
        if can_change_size(ks)
            return OffsetAxis(fst - one(fst), SimpleAxis(OneToMRange(length(ks))))
        else
            return OffsetAxis(fst - one(fst), SimpleAxis(One():static_length(ks)))
        end
    end

    OffsetAxis(axis::OffsetAxis) = axis
end

@inline Base.getproperty(axis::OffsetAxis, k::Symbol) = getproperty(parent(axis), k)

ArrayInterface.known_first(::Type{T}) where {T<:OffsetAxis{<:Any,<:Any,<:Any}} = nothing
function ArrayInterface.known_first(::Type{T}) where {Inds,F,T<:OffsetAxis{<:Any,Inds,StaticInt{F}}}
    if known_first(Inds) === nothing
        return nothing
    else
        return known_first(Inds) + F
    end
end
Base.first(axis::OffsetAxis) = first(parent(axis)) + getfield(axis, :offset)

ArrayInterface.known_last(::Type{T}) where {T<:OffsetAxis{<:Any,<:Any,<:Any}} = nothing
function ArrayInterface.known_last(::Type{T}) where {Inds,F,T<:OffsetAxis{<:Any,Inds,StaticInt{F}}}
    if known_last(Inds) === nothing
        return nothing
    else
        return known_last(Inds) + F
    end
end
Base.last(axis::OffsetAxis) = last(parent(axis)) + getfield(axis, :offset)

function ArrayInterface.unsafe_reconstruct(axis::OffsetAxis, inds)
    if inds isa AbstractOffsetAxis
        f_axis = offsets(axis, 1)
        f_inds = offsets(inds, 1)
        if f_axis === f_inds
            return OffsetAxis(offsets(axis, 1), unsafe_reconstruct(parent(axis), parent(inds)))
        else
            return OffsetAxis(f_axis + f_inds, unsafe_reconstruct(parent(axis), parent(inds)))
        end
    else
        return OffsetAxis(getfield(axis, :offset), unsafe_reconstruct(parent(axis), inds))
    end
end

"""
    offset([collection,] x)

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisArray(ones(3), offset(2))
3-element AxisArray(::Vector{Float64}
  • axes:
     1 = 3:5
)
     1
  3  1.0
  4  1.0
  5  1.0

```
"""
struct Offset <: AxisInitializer end
axis_method(::Offset, x, inds) = OffsetAxis(x, inds)

const offset = Offset()

"""
    OffsetArray

An array whose axes are all `OffsetAxis`
"""
const OffsetArray{T,N,P,A<:Tuple{Vararg{<:OffsetAxis}}} = AxisArray{T,N,P,A}

"""
    OffsetVector

A vector whose axis is `OffsetAxis`
"""
const OffsetVector{T,P<:AbstractVector{T},Ax1<:OffsetAxis} = OffsetArray{T,1,P,Tuple{Ax1}}

OffsetArray(A::AbstractArray{T,N}, inds::Vararg) where {T,N} = OffsetArray(A, inds)

OffsetArray(A::AbstractArray{T,N}, inds::Tuple) where {T,N} = OffsetArray{T,N}(A, inds)

OffsetArray(A::AbstractArray{T,0}, ::Tuple{}) where {T} = AxisArray(A)
OffsetArray(A::AbstractArray) = OffsetArray(A, axes(A))

# OffsetVector constructors
OffsetVector(A::AbstractVector, arg) = OffsetArray{eltype(A)}(A, arg)

function OffsetVector{T}(init::ArrayInitializer, arg) where {T}
    return OffsetVector(Vector{T}(init, length(arg)), arg)
end
OffsetArray{T}(A, inds::Tuple) where {T} = OffsetArray{T,length(inds)}(A, inds)
OffsetArray{T}(A, inds::Vararg) where {T} = OffsetArray{T,length(inds)}(A, inds)

function OffsetArray{T,N}(init::ArrayInitializer, inds::Tuple=()) where {T,N}
    return OffsetArray{T,N}(Array{T,N}(init, map(length, inds)), inds)
end
OffsetArray{T,N}(A, inds::Vararg) where {T,N} = OffsetArray{T,N}(A, inds)
function OffsetArray{T,N}(A::AbstractArray{T,N}, inds::Tuple) where {T,N}
    return OffsetArray{T,N,typeof(A)}(A, inds)
end
function OffsetArray{T,N}(A::AbstractArray{T2,N}, inds::Tuple) where {T,T2,N}
    return OffsetArray{T,N}(copyto!(Array{T}(undef, size(A)), A), inds)
end
function OffsetArray{T,N}(A::AxisArray, inds::Tuple) where {T,N}
    return OffsetArray{T,N}(parent(A), inds)
end

function OffsetArray{T,N,P}(A::AbstractArray, inds::NTuple{M,Any}) where {T,N,P<:AbstractArray{T,N},M}
    return OffsetArray{T,N,P}(convert(P, A))
end

OffsetArray{T,N,P}(A::OffsetArray{T,N,P}) where {T,N,P} = A

function OffsetArray{T,N,P}(A::OffsetArray) where {T,N,P}
    return initialize_axis_array(convert(P, parent(A)), axes(A))
end

function OffsetArray{T,N,P}(A::P, inds::Tuple{Vararg{<:Any,N}}) where {T,N,P<:AbstractArray{T,N}}
    if N === 1
        if can_change_size(P)
            axs = (OffsetAxis(first(inds), SimpleAxis(OneToMRange(axes(A, 1)))),)
        else
            axs = (OffsetAxis(first(inds), axes(A, 1)),)
        end
    else
        axs = map((f, axis) -> OffsetAxis(f, axis), inds, axes(A))
    end
    return initialize_axis_array(A, axs)
end

function print_axis(io::IO, axis::OffsetAxis)
    print(io, "offset($(Int(getfield(axis, :offset))))($(parent(axis)))")
end


text:
# minimal keyboard event support
"""
    read_key() -> control_value

read control key from keyboard input.

# Reference table

| value               | control_value     | effect                 |
| ------------------- | ----------------- | -------------------    |
| UP, LEFT, b         | :CONTROL_BACKWARD | show previous frame    |
| DOWN, RIGHT, f      | :CONTROL_FORWARD  | show next frame        |
| SPACE, p            | :CONTROL_PAUSE    | pause/resume play      |
| CTRL-c, q           | :CONTROL_EXIT     | exit current play      |
| others...           | :CONTROL_VOID     | no effect              |
"""
function read_key(io=stdin)
    control_value = :CONTROL_VOID
    try
        _setraw!(io, true)
        keyin = read(io, Char)
        if keyin == '\e'
            # some special keys are more than one byte, e.g., left key is `\e[D`
            # reference: https://en.wikipedia.org/wiki/ANSI_escape_code
            keyin = read(io, Char)
            if keyin == '['
                keyin = read(io, Char)
                if keyin in ['A', 'D'] # up, left
                    control_value = :CONTROL_BACKWARD
                elseif keyin in ['B', 'C'] # down, right
                    control_value = :CONTROL_FORWARD
                end
            end
        elseif 'A' <= keyin <= 'Z' || 'a' <= keyin <= 'z'
            keyin = lowercase(keyin)
            if keyin == 'p'
                control_value = :CONTROL_PAUSE
            elseif keyin == 'q'
                control_value = :CONTROL_EXIT
            elseif keyin == 'f'
                control_value = :CONTROL_FORWARD
            elseif keyin == 'b'
                control_value = :CONTROL_BACKWARD
            end
        elseif keyin == ' '
            control_value = :CONTROL_PAUSE
        end
    catch e
        if e isa InterruptException # Ctrl-C
            control_value = :CONTROL_EXIT
        else
            rethrow(e)
        end
    finally
        _setraw!(io, false)
    end
    return control_value
end

_setraw!(io::Base.TTY, raw) = ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid},Int32), io.handle, raw)
_setraw!(::IO, raw) = nothing

# keytask handling
# This works, although I have no idea how it is working internally -- Johnny Chen
# https://discourse.julialang.org/t/how-to-kill-thread/34236/8
if VERSION >= v"1.3"
    struct Stop end
    struct Continue end
    function safe_hanging(ch::Channel)
        signal = Continue()
        while true
            isready(ch) && (signal = take!(ch))
            signal == Stop() && break
            sleep(1e-4)
        end
    end

    make_channel(task) = Channel(safe_hanging, taskref=task, spawn=true)
    function safe_kill(task, ch::Channel)
        put!(ch, Stop())
        wait(task)
    end
else
    # non-safe quit
    make_channel(task) = task
    function safe_kill(task, ch)
        # this is not safe at all
        @async Base.throwto(keytask, InterruptException())
        wait(task)
    end
end

ext: jl


text:
"""
	2D
"""

function integralArray!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,2}, fun=(x)->(x) )

	@inbounds for c in 2:size(img,2)+1, r in 2:size(img,1)+1
		intArr[r,c] = fun(img[r-1,c-1]) + intArr[r-1,c] + intArr[r,c-1] - intArr[r-1,c-1]
	end
	return intArr
end

function integralArray( img::AbstractArray{<:Real,2}; type=Float32, fun=(x)->(x) )

	return integralArray!( img, zeros( type, size(img) .+ 1 ), fun )
end

function integralArea( intArr::AbstractArray{<:Real,2}, TL, BR )
	TL   = TL .+ 1;
	BR   = BR .+ 1;
	area = intArr[BR[1],BR[2]] - intArr[BR[1],TL[2]] - intArr[TL[1],BR[2]] + intArr[TL[1],TL[2]]
end

# Reversed variant, which is used during normalized cross-correlation
function integralArrayRev!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,2}, fun=(x)->(x) )
	h, w = size(img)
	@inbounds for c in 2:w+1, r in 2:h+1
		intArr[r,c] = fun(img[h-(r-2),w-(c-2)]) + intArr[r-1,c] + intArr[r,c-1] - intArr[r-1,c-1]
	end
	return intArr
end

function integralArrayRev( img::AbstractArray{<:Real,2}; type=Float32, fun=(x)->(x) )

	return integralArrayRev!( img, zeros( type, size(img) .+ 1 ), fun )
end

"""
	3D: optimized integral array by reusing previous operations
"""
# Below is the unoptimized code, in case I want to check again:
#    for z in 2:d+1, c in 2:w+1, r in 2:h+1
#        arr[r,c,z] = fun(vol[r-1,c-1,z-1]) + arr[r-1,c,z] + arr[r,c-1,z] + arr[r,c,z-1] -
#                     arr[r,c-1,z-1] - arr[r-1,c,z-1] - arr[r-1,c-1,z] + arr[r-1,c-1,z-1]
#    end

function integralArray!( vol::AbstractArray{<:Real,3}, intArr::AbstractArray{<:AbstractFloat,3}, fun=(x)->(x) )

	@inbounds for z in 2:size(vol,3)+1, c in 2:size(vol,2)+1
		reuse = 0.0
		for r in 2:size(vol,1)+1
			tmp = fun(vol[r-1,c-1,z-1])
			intArr[r,c,z] = tmp + intArr[r,c-1,z] + intArr[r,c,z-1] - intArr[r,c-1,z-1] + reuse;
			reuse += tmp
		end
	end
	return intArr
end

function integralArray( vol::AbstractArray{<:Real,3}; type=Float32, fun=(x)->(x) )

	return integralArray!( vol, zeros( type, size(vol) .+ 1 ), fun )
end

function integralArea( intArr::AbstractArray{<:Real,3}, TL, BR )
	TL = TL .+ 1;
	BR = BR .+ 1;
	area  = intArr[BR[1],BR[2],BR[3]] - intArr[TL[1],TL[2],TL[3]]
	area -= intArr[TL[1],BR[2],BR[3]] + intArr[BR[1],TL[2],BR[3]] + intArr[BR[1],BR[2],TL[3]]
	area += intArr[BR[1],TL[2],TL[3]] + intArr[TL[1],BR[2],TL[3]] + intArr[TL[1],TL[2],BR[3]]
    return area
end

# Reversed variant, which is used during normalized cross-correlation
function integralArrayRev!( vol::AbstractArray{<:Real,3}, intArr::AbstractArray{<:AbstractFloat,3}, fun=(x)->(x) )

	h, w, d = size(vol)
	@inbounds for z in 2:d+1, c in 2:w+1
		reuse = 0.0
		for r in 2:h+1
			tmp = fun(vol[h-(r-2),w-(c-2),d-(z-2)])
			intArr[r,c,z] = tmp + intArr[r,c-1,z] + intArr[r,c,z-1] - intArr[r,c-1,z-1] + reuse;
			reuse += tmp
		end
	end
	return intArr
end

function integralArrayRev( vol::AbstractArray{<:Real,3}; type=Float32, fun=(x)->(x) )

	return integralArrayRev!( vol, zeros( type, size(vol) .+ 1 ), fun )
end

"""
	Derived functions
"""
function scan4Simmetry( intArr, rad )

	score = zeros( length( rad+1:size(intArr,1)-rad), length( rad+1:size(intArr,2)-rad ) )

	cont = 1;
	for c in rad+1:size(intArr,2)-rad;     c1, c2 = c - rad, c + rad
		for r in rad+1:size(intArr,1)-rad; r1, r2 = r - rad, r + rad

			Q1  = integralArea( intArr, (r1,c ), ( r-1,c2-1) )
			Q2  = integralArea( intArr, (r1,c1), ( r-1,c -1) )
			Q3  = integralArea( intArr, (r ,c1), (r2-1,c -1) )
			Q4  = integralArea( intArr, (r ,c ), (r2-1,c2-1) )

			acc = Q1 + Q2 + Q3 + Q4;
			ps  = [ Q1, Q2, Q3, Q4 ] ./ acc

			score[ cont ] = 1/( 1 + sum( abs.( ps .- 0.25 ) ) )
			score[ cont ] = isnan( score[cont] ) ? 0.0 : score[cont];
			cont += 1
	end	end

	return score
end

function followCentroids( data, intArr, rad; iters=5 )

	# Goal: At any position, compute centroid, move to centroid
	# recompute centroid and so on until converging

	score = zeros( Int64, size(data) )

	xx = [ x for y in -rad:rad, x in -rad:rad ]
	yy = [ y for y in -rad:rad, x in -rad:rad ]

	h, w = size( data );
	for c in 1:w;
		for r in 1:h;

			pos = r, c;
			#println( pos )
			# Contemplating the case where the position is closer than "rad"
			# to the edges of the images. As we iterate and update the position
			# of the centroid, it becomes necessary. Here is not needed though.
			left   = rad + min(0, pos[1]-1-rad)
			up     = rad + min(0, pos[2]-1-rad)
			right  = min(  h - pos[ 1 ], rad  )
			down   = min(  w - pos[ 2 ], rad  )

			TL     = pos .- ( left , up   ) .- 1
			BR     = pos .+ ( right, down )

			acc    = integralArea( intArr, TL, BR )
			( acc == 0 ) && ( continue; )

			subimg = pos[1]-left:pos[1]+right, pos[2]-up:pos[2]+down
			subxy  =  rad+1-left:rad+1+right ,  rad+1-up:rad+1+down

			xmean  = sum( (data[ subimg... ] ./ acc) .* xx[ subxy... ] )
			ymean  = sum( (data[ subimg... ] ./ acc) .* yy[ subxy... ] )

			for i in 1:iters

				offset = round.( Int64, (ymean,xmean) )

				if offset[1] == 0 && offset[2] == 0
					break
				end

				pos    = pos .+ offset
				#println( "\t", pos )
				left   = rad + min(0, pos[1]-1-rad)
				up     = rad + min(0, pos[2]-1-rad)
				right  = min(  h - pos[ 1 ], rad  )
				down   = min(  w - pos[ 2 ], rad  )

				TL     = pos .- ( left , up   ) .- 1
				BR     = pos .+ ( right, down )

				acc    = integralArea( intArr, TL, BR )
				( acc == 0 ) && ( break; )

				subimg = pos[1]-left:pos[1]+right, pos[2]-up:pos[2]+down
				subxy  =  rad+1-left:rad+1+right ,  rad+1-up:rad+1+down

				xmean  = sum( (data[ subimg... ] ./ acc) .* xx[ subxy... ] )
				ymean  = sum( (data[ subimg... ] ./ acc) .* yy[ subxy... ] )
			end

			pos = pos .+ round.( Int64, (ymean,xmean) )

			score[ pos... ] += 1;
	end	end

	return score
end

function integralAreaZNCC( sumS::Array{T,2}, sumS2::Array{T,2},
                           row, col, ss, si ) where {T<:Real}
    pad = 1;
    row = row + pad;
    col = col + pad;
    r0, r1 = max( 1, row - si[1] ), min( ss[1] + pad, row );
    c0, c1 = max( 1, col - si[2] ), min( ss[2] + pad, col );

    opS  = 0.0
    opS += sumS[ r1, c1 ];
    opS -= sumS[ r1, c0 ];
    opS -= sumS[ r0, c1 ];
    opS += sumS[ r0, c0 ];

    opS2  = 0.0
    opS2 += sumS2[ r1, c1 ];
    opS2 -= sumS2[ r1, c0 ];
    opS2 -= sumS2[ r0, c1 ];
    opS2 += sumS2[ r0, c0 ];

    return opS, opS2
end


text:
using MLStyle

"""
    @yao_str
    yao"..."

The mark up language for quantum circuit.
"""
macro yao_str(str::String)
    yaofromstring(str)
end

function yaofromstring(x::String)
    ex = Meta.parse(x)
    parse_yaoscript(ex)
end

macro yaoscript(ex)
    parse_yaoscript(ex)
end

function parse_yaoscript(ex)
    @match ex begin
        Expr(:let, header, body) => begin
            info = ParseInfo(-1, "")
            if header.head != :block
                header = [header]
            else
                header = header.args
            end
            for hd in header
                @match hd begin
                    :(nqubits = $n) => (info.nbit = Int(n))
                    :(version = $v) => (info.version = String(v))
                    _ => error("unknown configuration $header")
                end
            end
            parse_ex(body |> rmlines, info)
        end
        _ => error("wrong format, expect expression like `let nqubits=5 GATEDEF end`, got $ex")
    end
end

yaofromfile(x::String) = yaofromstring(read(x, String))
parse_ex(ex, nbit::Int) = parse_ex(ex, ParseInfo(nbit, ""))

mutable struct ParseInfo
    nbit::Int
    version::String
end

function parse_param(x)
    @match x begin
        ::Number => x
        :(:π) => x
        _ => error("can not parse parameter $x.")
    end
end

function parse_ex(ex, info::ParseInfo)
    @match ex begin
        :(version = $vnumber) => (info.version = String(vnumber); nothing)
        :(nqubits = $x) => (info.nbit = Int(x); nothing)
        ::Nothing => nothing
        :($g') => :($(parse_ex(g, info))')
        :($a * $b) => :($(Number(a)) * $(parse_ex(b, info)))
        :(kron($(args...))) => :(kron($(parse_ex.(args, Ref(ParseInfo(1, info.version)))...)))
        :(repeat($(exloc...)) => $g) => begin
            loc = render_loc((exloc...,), info.nbit)
            :(repeat($(info.nbit), $(parse_ex(g, ParseInfo(1, info.version))), $loc))
        end
        :(cache($g)) => :(cache($(parse_ex(g, info))))
        :(rot($g, $theta)) => :(rot($(parse_ex(g, info)), $(parse_param(theta))))
        :(time($dt) => $h) => :(time_evolve($(parse_ex(h, info)), $(parse_param(dt))))
        :($exloc => Measure) => parse_ex(:($exloc => Measure(nothing) => nothing), info)
        :($exloc => Measure($op)) => parse_ex(:($exloc => Measure($op) => nothing), info)
        :($exloc => Measure => $post) => parse_ex(:($exloc => Measure(nothing) => $post), info)
        :($exloc => Measure($op) => $post) => begin
            locs = exloc == :ALL ? :(AllLocs()) : render_loc(exloc, info.nbit)
            op =
                op isa Nothing || op == :nothing ? :(ComputationalBasis()) :
                parse_ex(op, exloc == :ALL ? info : ParseInfo(length(locs), info.version))
            @match post begin
                ::Nothing || :nothing => :(Measure($(info.nbit); locs = $locs, operator = $(op)))
                :(resetto($(rbits...))) => begin
                    cb = bit_literal(render_bitstring.(rbits)...)
                    :(Measure($(info.nbit); locs = $locs, operator = $(op), resetto = $cb))
                end
                :remove => :(Measure($(info.nbit); locs = $locs, operator = $(op), remove = true))
            end
        end
        :(+($(args...))) => :(+($(parse_ex.(args, Ref(info))...)))
        :(focus($(exloc...)) => $g) => begin
            loc = render_loc((exloc...,), info.nbit)
            :(subroutine($(info.nbit), $(parse_ex(g, ParseInfo(length(loc), info.version))), $loc))
        end
        :(
            begin
                $(cargs...)
            end
        ) => begin
            args = filter(x -> x !== nothing, [parse_ex(arg, info) for arg in cargs])
            :(chain($(info.nbit), [$(args...)]))
        end
        :($exloc => $gate) => begin
            loc = render_loc(exloc, info.nbit)
            :(put($(info.nbit), $loc => $(parse_ex(gate, ParseInfo(length(loc), info.version)))))
        end
        :($(cargs...), $exloc => $gate) => begin
            loc = render_loc(exloc, info.nbit)
            cbits = render_cloc.(cargs, Ref(info))
            if cbits[1] isa Integer
                :(control(
                    $(info.nbit),
                    $cbits,
                    $loc => $(parse_ex(gate, ParseInfo(length(loc), info.version))),
                ))
            else
                :(kron($(info.nbit), $(cbits...), $loc => $(parse_ex(gate, ParseInfo(1, info.version)))))
            end
        end
        :($f($(args...))) => gate_expr(Val(Symbol(f)), args, info)

        ::LineNumberNode => nothing
        ::Symbol => ex  # const gate
        _ => error("scipt format error! got $ex of type $(typeof(ex))")
    end
end

function check_dumpload(gate::AbstractBlock{N}) where {N}
    gate2 = eval(parse_ex(dump_gate(gate), N))
    gate2 == gate || mat(gate2) ≈ mat(gate)
end

render_bitstring(ex) = @match ex begin
    ::Number => begin
        if ex == 1 || ex == 0
            ex
        else
            error("expect a bitstring like `1` or `(1,0)`, got $ex")
        end
    end
    ::Tuple => render_bitstring.(ex)
    :($(args...),) => (render_bitstring.(args)...,)
    _ => error("expect a bitstring like `1` or `(1,0)`, got $ex")
end

render_loc(ex, nbit::Int) = @match ex begin
    :($(args...),) => (render_loc.(args, nbit)...,)
    ::Number => Int(ex)
    :($a:$b) => Int(a):Int(b)
    :ALL => (1:nbit...,)
    :($a:$step:$b) => Int(a):Int(step):Int(b)
    ::Tuple => Int.(ex)
    ::UnitRange{Int} => ex
    _ => error("expect a location specification like `2`, `2:5` or `(2,3)`, got $ex")
end

render_cloc(ex, info) = @match ex begin
    :($a => C($b)) => begin
        if !(b == 1 || b == 0)
            error("expect a control values `0` or `1`, got $ex")
        end
        Int(a) * (2 * Int(b) - 1)
    end
    :($a => C) => render_cloc(:($a => C(1)), info)
    :($a => $g) => :($(render_loc(a, info.nbit)) => $(parse_ex(g, ParseInfo(1, info.version))))
    _ => error("expect a control location specification like `2=>0` or `3=>1`, got $ex")
end

"""
    gate_expr(::Val{G}, args, info)

Obtain the gate constructior from its YaoScript expression.
`G` is a symbol for the gate type,
the default constructor is `G(args...)`.
`info` contains the informations about the number of qubit and Yao version.
"""
function gate_expr(::Val{G}, args, info) where {G}
    throw(NotImplementedError(:gate_expr, (Val(G), args, info)))
end

render_arg(ex, info) = @match ex begin
    # locs
    :($(args...),) => (render_loc.(args, info.nbit)...,)
    ::Integer => Int(ex)
    :($a:$b) => Int(a):Int(b)
    :ALL => (1:nbit...,)
    :($a:$step:$b) => Int(a):Int(step):Int(b)
    ::Tuple => Int.(ex)

    # floating point numbers
    ::Number => Number(ex)

    # Pair of (loc => gate)
    :($a => $g) => begin
        locs = render_loc(a, info.nbit)
        :($(locs) => $(parse_ex(g, ParseInfo(length(locs), info.version))))
    end

    # try parse a gate
    _ => parse_ex(ex, info)
end

for (G, F) in [(:ShiftGate, :shift), (:PhaseGate, :phase)]
    fname = QuoteNode(F)
    @eval function dump_gate(blk::$G)
        Expr(:call, $fname, tokenize_param(blk.theta))
    end
    @eval function gate_expr(::Val{$(QuoteNode(F))}, args, info)
        Expr(:call, $fname, parse_param(args[1]))
    end
end

max_stars_repo_name: QuantumBFS/YaoStdCircuit

max_issues_repo_name: QuantumBFS/YaoStdCircuit


text:
"""
   update_feed!(state, agent_idx, config)

Update an agent's feed by applying deletions, weight decays and sorting after post weight

# Arguments
- `state`: a tuple of the current graph and agent_list
- `agent_idx`: agent index
- `config`: Config object as provided by Config()

See also: [`Config`](@ref), [`update__perceiv_publ_opinion!`](@ref)
"""
function update_feed!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config
)

    graph, agent_list = state
    this_agent = agent_list[agent_idx]

    unique!(this_agent.feed)
    deleted_posts = Integer[]

    for (index, post) in enumerate(this_agent.feed)
        if (
            post.weight == -1
            || (
                !(post.source_agent in inneighbors(graph, agent_idx))
                && post.share_count == 0)
            )
            push!(deleted_posts, index)
        else
            post.weight = config.feed_props.post_decay * post.weight
        end
    end

    deleteat!(this_agent.feed, deleted_posts)
    sort!(this_agent.feed, lt=isless, rev=true)

    if length(this_agent.feed) > config.feed_props.feed_size
        this_agent.feed = this_agent.feed[1:config.feed_props.feed_size]
    end

    for post in this_agent.feed
        if !(agent_idx in post.seen_by)
            push!(post.seen_by, agent_idx)
        end
    end

    return state
end

# suppress output of include()
;


text:
module sparsey

include("hamiltonian.jl")

import Plots

export diag
export plotEig

struct HOscParams
    omegaXY::Array{Float64}
    mass::Float64
    gridDivs::Array{Integer}
    gridScale::Array{Float64}
end

function diag(osc::HOscParams, numEig::Integer)
    # Set number of eigenvalues, and spatial dim of grid

    Ngx = [osc.gridDivs[1], osc.gridDivs[2]]
    dim = length(Ngx)
    gMax = [osc.gridScale[1], osc.gridScale[2]]
    gMin = [-osc.gridScale[1], -osc.gridScale[2]]

    g0 = SpatialGrid(Ngx, gMax, gMin)
    V = VHarm(g0, Ngx, [osc.omegaXY[1], osc.omegaXY[2]], osc.mass)

    # Returns sparse diagonal matrix with linear indexed (lexicographical)
    # values of the vectors V, with dimensions Ngx
    VD=kronSum(Ngx,V)

    dx = abs.(g0[1][1] - g0[1][2])
    KD = -(1.0/(2.0*osc.mass))*l2d(Ngx)./dx^2

    H = VD + KD
    val,vec = eigs(H, nev=numEig, which=:SM); #Calculate 8 of the smallest eigenvectors
    return val, vec
end

function plotEig(val::Array{Float64}, vec::Array{Float64}, osc::HOscParams, numEig::Integer)
    # Reshape and plot results
    p1 = Plots.plot(val, lw=2, legend = :none)
    plot_arr = Array{typeof(p1)}(undef, numEig+1);
    for i in 1:numEig
        plot_arr[i] = Plots.heatmap(  reshape( abs.(vec[:,i]).^2, osc.gridDivs[1], osc.gridDivs[2]),
                                legend = :none,
                                aspect_ratio=:equal,
                                grid=false,
                                border=nothing,
                                axiscolor=nothing )
    end
    plot_arr[1 + numEig] = p1;
    Plots.plot(plot_arr..., size=(800,800))
end

end


text:
# Statistical Models

abstract type StatisticalModel end

"""
    coef(obj::StatisticalModel)

Return the coefficients of the model.
"""
coef(obj::StatisticalModel) = error("coef is not defined for $(typeof(obj)).")

"""
    coefnames(obj::StatisticalModel)

Return the names of the coefficients.
"""
coefnames(obj::StatisticalModel) = error("coefnames is not defined for $(typeof(obj)).")

"""
    coeftable(obj::StatisticalModel; level::Real=0.95)

Return a table of class `CoefTable` with coefficients and related statistics.
`level` determines the level for confidence intervals (by default, 95%).
"""
coeftable(obj::StatisticalModel) = error("coeftable is not defined for $(typeof(obj)).")

"""
    confint(obj::StatisticalModel; level::Real=0.95)

Compute confidence intervals for coefficients, with confidence level `level` (by default 95%).
"""
confint(obj::StatisticalModel) = error("confint is not defined for $(typeof(obj)).")

"""
    deviance(obj::StatisticalModel)

Return the deviance of the model relative to a reference, which is usually when applicable
the saturated model. It is equal, *up to a constant*, to ``-2 \\log L``, with ``L``
the likelihood of the model.
"""
deviance(obj::StatisticalModel) = error("deviance is not defined for $(typeof(obj)).")

"""
    islinear(obj::StatisticalModel)

Indicate whether the model is linear.
"""
islinear(obj::StatisticalModel) = error("islinear is not defined for $(typeof(obj)).")

"""
    nulldeviance(obj::StatisticalModel)

Return the deviance of the null model, that is the one including only the intercept.
"""
nulldeviance(obj::StatisticalModel) = error("nulldeviance is not defined for $(typeof(obj)).")

"""
    loglikelihood(obj::StatisticalModel)

Return the log-likelihood of the model.
"""
loglikelihood(obj::StatisticalModel) = error("loglikelihood is not defined for $(typeof(obj)).")

"""
    loglikelihood(obj::StatisticalModel)

Return the log-likelihood of the null model corresponding to model `obj`.
This is usually the model containing only the intercept.
"""
nullloglikelihood(obj::StatisticalModel) = error("nullloglikelihood is not defined for $(typeof(obj)).")

"""
    score(obj::StatisticalModel)

Return the score of the statistical model. The score is the gradient of the
log-likelihood with respect to the coefficients.
"""
score(obj::StatisticalModel) = error("score is not defined for $(typeof(obj)).")

"""
    nobs(obj::StatisticalModel)

Return the number of independent observations on which the model was fitted. Be careful
when using this information, as the definition of an independent observation may vary
depending on the model, on the format used to pass the data, on the sampling plan
(if specified), etc.
"""
nobs(obj::StatisticalModel) = error("nobs is not defined for $(typeof(obj)).")

"""
    dof(obj::StatisticalModel)

Return the number of degrees of freedom consumed in the model, including
when applicable the intercept and the distribution's dispersion parameter.
"""
dof(obj::StatisticalModel) = error("dof is not defined for $(typeof(obj)).")

"""
    mss(obj::StatisticalModel)

Return the model sum of squares.
"""
mss(obj::StatisticalModel) = error("mss is not defined for $(typeof(obj)).")

"""
    rss(obj::StatisticalModel)

Return the residual sum of squares.
"""
rss(obj::StatisticalModel) = error("rss is not defined for $(typeof(obj)).")

"""
    informationmatrix(model::StatisticalModel; expected::Bool = true)

Return the information matrix. By default the Fisher information matrix is returned,
while the observed information matrix can be requested with `expected = false`.
"""
informationmatrix(model::StatisticalModel; expected::Bool = true) =
    error("informationmatrix is not defined for $(typeof(obj)).")

"""
    stderror(obj::StatisticalModel)

Return the standard errors for the coefficients of the model.
"""
stderror(obj::StatisticalModel) = sqrt.(diag(vcov(obj)))

"""
    vcov(obj::StatisticalModel)

Return the variance-covariance matrix for the coefficients of the model.
"""
vcov(obj::StatisticalModel) = error("vcov is not defined for $(typeof(obj)).")

"""
    weights(obj::StatisticalModel)

Return the weights used in the model.
"""
weights(obj::StatisticalModel) = error("weights is not defined for $(typeof(obj)).")

"""
    isfitted(obj::StatisticalModel)

Indicate whether the model has been fitted.
"""
isfitted(obj::StatisticalModel) = error("isfitted is not defined for $(typeof(obj)).")

"""
Fit a statistical model.
"""
fit(obj::StatisticalModel, args...) = error("fit is not defined for $(typeof(obj)).")

"""
Fit a statistical model in-place.
"""
fit!(obj::StatisticalModel, args...) = error("fit! is not defined for $(typeof(obj)).")

"""
    aic(obj::StatisticalModel)

Akaike's Information Criterion, defined as ``-2 \\log L + 2k``, with ``L`` the likelihood
of the model, and `k` its number of consumed degrees of freedom
(as returned by [`dof`](@ref)).
"""
aic(obj::StatisticalModel) = -2loglikelihood(obj) + 2dof(obj)

"""
    aicc(obj::StatisticalModel)

Corrected Akaike's Information Criterion for small sample sizes (Hurvich and Tsai 1989),
defined as ``-2 \\log L + 2k + 2k(k-1)/(n-k-1)``, with ``L`` the likelihood of the model,
``k`` its number of consumed degrees of freedom (as returned by [`dof`](@ref)),
and ``n`` the number of observations (as returned by [`nobs`](@ref)).
"""
function aicc(obj::StatisticalModel)
    k = dof(obj)
    n = nobs(obj)
    -2loglikelihood(obj) + 2k + 2k*(k+1)/(n-k-1)
end

"""
    bic(obj::StatisticalModel)

Bayesian Information Criterion, defined as ``-2 \\log L + k \\log n``, with ``L``
the likelihood of the model,  ``k`` its number of consumed degrees of freedom
(as returned by [`dof`](@ref)), and ``n`` the number of observations
(as returned by [`nobs`](@ref)).
"""
bic(obj::StatisticalModel) = -2loglikelihood(obj) + dof(obj)*log(nobs(obj))

"""
    r2(obj::StatisticalModel)
    r²(obj::StatisticalModel)

Coefficient of determination (R-squared).

For a linear model, the R² is defined as ``ESS/TSS``, with ``ESS`` the explained sum of squares
and ``TSS`` the total sum of squares.
"""
function r2(obj::StatisticalModel)
    Base.depwarn("The default r² method for linear models is deprecated. " *
                 "Packages should define their own methods.", :r2)

    mss(obj) / deviance(obj)
end

"""
    r2(obj::StatisticalModel, variant::Symbol)
    r²(obj::StatisticalModel, variant::Symbol)

Pseudo-coefficient of determination (pseudo R-squared).

For nonlinear models, one of several pseudo R² definitions must be chosen via `variant`.
Supported variants are:
- `:MacFadden` (a.k.a. likelihood ratio index), defined as ``1 - \\log (L)/\\log (L_0)``;
- `:CoxSnell`, defined as ``1 - (L_0/L)^{2/n}``;
- `:Nagelkerke`, defined as ``(1 - (L_0/L)^{2/n})/(1 - L_0^{2/n})``.

In the above formulas, ``L`` is the likelihood of the model,
``L_0`` is the likelihood of the null model (the model with only an intercept),
``n`` is the number of observations, ``y_i`` are the responses,
``\\hat{y}_i`` are fitted values and ``\\bar{y}`` is the average response.

Cox and Snell's R² should match the classical R² for linear models.
"""
function r2(obj::StatisticalModel, variant::Symbol)
    ll = loglikelihood(obj)
    ll0 = nullloglikelihood(obj)
    if variant == :McFadden
        1 - ll/ll0
    elseif variant == :CoxSnell
        1 - exp(2 * (ll0 - ll) / nobs(obj))
    elseif variant == :Nagelkerke
        (1 - exp(2 * (ll0 - ll) / nobs(obj))) / (1 - exp(2 * ll0 / nobs(obj)))
    else
        error("variant must be one of :McFadden, :CoxSnell or :Nagelkerke")
    end
end

const r² = r2

"""
    adjr2(obj::StatisticalModel)
    adjr²(obj::StatisticalModel)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as ``1 - (1 - (1-R^2)(n-1)/(n-p))``, with ``R^2``
the coefficient of determination, ``n`` the number of observations, and ``p`` the number of
coefficients (including the intercept). This definition is generally known as the Wherry Formula I.
"""
adjr2(obj::StatisticalModel) = error("adjr2 is not defined for $(typeof(obj)).")

"""
    adjr2(obj::StatisticalModel, variant::Symbol)
    adjr²(obj::StatisticalModel, variant::Symbol)

Adjusted pseudo-coefficient of determination (adjusted pseudo R-squared).

For nonlinear models, one of the several pseudo R² definitions must be chosen via `variant`.
The only currently supported variant is `:MacFadden`, defined as ``1 - (\\log (L) - k)/\\log (L0)``.
In this formula, ``L`` is the likelihood of the model, ``L0`` that of the null model
(the model including only the intercept), and ``k`` is the number of consumed degrees of freedom
of the model (as returned by [`dof`](@ref)).
"""
function adjr2(obj::StatisticalModel, variant::Symbol)
    ll = loglikelihood(obj)
    ll0 = nullloglikelihood(obj)
    k = dof(obj)
    if variant == :McFadden
        1 - (ll - k)/ll0
    else
        error(":McFadden is the only currently supported variant")
    end
end

const adjr² = adjr2

abstract type RegressionModel <: StatisticalModel end

"""
    fitted(obj::RegressionModel)

Return the fitted values of the model.
"""
fitted(obj::RegressionModel) = error("fitted is not defined for $(typeof(obj)).")

"""
    response(obj::RegressionModel)

Return the model response (a.k.a. the dependent variable).
"""
response(obj::RegressionModel) = error("response is not defined for $(typeof(obj)).")

"""
    responsename(obj::RegressionModel)

Return the name of the model response (a.k.a. the dependent variable).
"""
responsename(obj::RegressionModel) = error("responsename is not defined for $(typeof(obj)).")

"""
    meanresponse(obj::RegressionModel)

Return the mean of the response.
"""
meanresponse(obj::RegressionModel) = error("meanresponse is not defined for $(typeof(obj)).")

"""
    modelmatrix(obj::RegressionModel)

Return the model matrix (a.k.a. the design matrix).
"""
modelmatrix(obj::RegressionModel) = error("modelmatrix is not defined for $(typeof(obj)).")

"""
    crossmodelmatrix(obj::RegressionModel)

Return `X'X` where `X` is the model matrix of `obj`.
This function will return a pre-computed matrix stored in `obj` if possible.
"""
crossmodelmatrix(obj::RegressionModel) = (x = modelmatrix(obj); Symmetric(x' * x))

"""
    leverage(obj::RegressionModel)

Return the diagonal of the projection matrix.
"""
leverage(obj::RegressionModel) = error("leverage is not defined for $(typeof(obj)).")

"""
    residuals(obj::RegressionModel)

Return the residuals of the model.
"""
residuals(obj::RegressionModel) = error("residuals is not defined for $(typeof(obj)).")

"""
    predict(obj::RegressionModel, [newX])

Form the predicted response of model `obj`. An object with new covariate values `newX` can be supplied,
which should have the same type and structure as that used to fit `obj`; e.g. for a GLM
it would generally be a `DataFrame` with the same variable names as the original predictors.
"""
function predict end

predict(obj::RegressionModel) = error("predict is not defined for $(typeof(obj)).")

"""
    predict!

In-place version of [`predict`](@ref).
"""
function predict! end

predict!(obj::RegressionModel) = error("predict! is not defined for $(typeof(obj)).")

"""
    dof_residual(obj::RegressionModel)

Return the residual degrees of freedom of the model.
"""
dof_residual(obj::RegressionModel) = error("dof_residual is not defined for $(typeof(obj)).")

"""
    params(obj)

Return all parameters of a model.
"""
params(obj) = error("params is not defined for $(typeof(obj))")
function params! end

## coefficient tables with specialized show method

mutable struct CoefTable
    cols::Vector
    colnms::Vector
    rownms::Vector
    pvalcol::Int
    teststatcol::Int
    function CoefTable(cols::Vector,colnms::Vector,rownms::Vector,
                       pvalcol::Int=0,teststatcol::Int=0)
        nc = length(cols)
        nrs = map(length,cols)
        nr = nrs[1]
        length(colnms) in [0,nc] || throw(ArgumentError("colnms should have length 0 or $nc"))
        length(rownms) in [0,nr] || throw(ArgumentError("rownms should have length 0 or $nr"))
        all(nrs .== nr) || throw(ArgumentError("Elements of cols should have equal lengths, but got $nrs"))
        pvalcol in 0:nc || throw(ArgumentError("pvalcol should be between 0 and $nc"))
        teststatcol in 0:nc || throw(ArgumentError("teststatcol should be between 0 and $nc"))
        new(cols,colnms,rownms,pvalcol,teststatcol)
    end

    function CoefTable(mat::Matrix,colnms::Vector,rownms::Vector,
                       pvalcol::Int=0,teststatcol::Int=0)
        nc = size(mat,2)
        cols = Any[mat[:, i] for i in 1:nc]
        CoefTable(cols,colnms,rownms,pvalcol,teststatcol)
    end
end

"""
Show a p-value using 6 characters, either using the standard 0.XXXX
representation or as <Xe-YY.
"""
struct PValue
    v::Real
    function PValue(v::Real)
        0 <= v <= 1 || isnan(v) || error("p-values must be in [0; 1]")
        new(v)
    end
end

function show(io::IO, pv::PValue)
    v = pv.v
    if isnan(v)
        @printf(io,"%d", v)
    elseif v >= 1e-4
        @printf(io,"%.4f", v)
    else
        @printf(io,"<1e%2.2d", ceil(Integer, max(nextfloat(log10(v)), -99)))
    end
end

"""Show a test statistic using 2 decimal digits"""
struct TestStat <: Real
    v::Real
end

show(io::IO, x::TestStat) = @printf(io, "%.2f", x.v)

"""Wrap a string so that show omits quotes"""
struct NoQuote
    s::String
end

show(io::IO, n::NoQuote) = print(io, n.s)

function show(io::IO, ct::CoefTable)
    cols = ct.cols; rownms = ct.rownms; colnms = ct.colnms;
    nc = length(cols)
    nr = length(cols[1])
    if length(rownms) == 0
        rownms = [lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
    end
    mat = [j == 1 ? NoQuote(rownms[i]) :
           j-1 == ct.pvalcol ? PValue(cols[j-1][i]) :
           j-1 in ct.teststatcol ? TestStat(cols[j-1][i]) :
           cols[j-1][i] isa AbstractString ? NoQuote(cols[j-1][i]) : cols[j-1][i]
           for i in 1:nr, j in 1:nc+1]
    # Code inspired by print_matrix in Base
    io = IOContext(io, :compact=>true, :limit=>false)
    A = Base.alignment(io, mat, 1:size(mat, 1), 1:size(mat, 2),
                       typemax(Int), typemax(Int), 3)
    nmswidths = pushfirst!(length.(colnms), 0)
    A = [nmswidths[i] > sum(A[i]) ? (A[i][1]+nmswidths[i]-sum(A[i]), A[i][2]) : A[i]
         for i in 1:length(A)]
    totwidth = sum(sum.(A)) + 2 * (length(A) - 1)
    println(io, repeat('─', totwidth))
    print(io, repeat(' ', sum(A[1])))
    for j in 1:length(colnms)
        print(io, "  ", lpad(colnms[j], sum(A[j+1])))
    end
    println(io, '\n', repeat('─', totwidth))
    for i in 1:size(mat, 1)
        Base.print_matrix_row(io, mat, A, i, 1:size(mat, 2), "  ")
        i != size(mat, 1) && println(io)
    end
    print(io, '\n', repeat('─', totwidth))
    nothing
end

"""
    ConvergenceException(iters::Int, lastchange::Real=NaN, tol::Real=NaN)

The fitting procedure failed to converge in `iters` number of iterations,
i.e. the `lastchange` between the cost of the final and penultimate iteration was greater than
specified tolerance `tol`.
"""
struct ConvergenceException{T<:Real} <: Exception
    iters::Int
    lastchange::T
    tol::T
    msg::String
    function ConvergenceException{T}(iters, lastchange::T, tol::T, msg::String) where T<:Real
        if tol > lastchange
            throw(ArgumentError("Change must be greater than tol."))
        else
            new(iters, lastchange, tol, msg)
        end
    end
end

ConvergenceException(iters, lastchange::T=NaN, tol::T=NaN,
                     msg::AbstractString="") where {T<:Real} =
    ConvergenceException{T}(iters, lastchange, tol, String(msg))

function Base.showerror(io::IO, ce::ConvergenceException)
    print(io, "failure to converge after $(ce.iters) iterations.")
    if !isnan(ce.lastchange)
        print(io, " Last change ($(ce.lastchange)) was greater than tolerance ($(ce.tol)).")
    end
    if !isempty(ce.msg)
        print(io, ' ', ce.msg)
    end
end


text:
module Genesys
# TODO: break into submodules
# Optimisation
using JuMP, Cbc, Metaheuristics, SDDP
# Math
using Statistics, StatsBase, MultivariateStats, Clustering, Distributions, Distances, LinearAlgebra
# Others
using Seaborn, ProgressMeter, Dates, Distributed, SharedArrays, CSV, DataFrames, JLD
# Assets
include(joinpath("assets","microgrid.jl"))
include(joinpath("assets","carriers.jl"))
include(joinpath("assets","liion.jl"))
include(joinpath("assets","tes.jl"))
include(joinpath("assets","h2tank.jl"))
include(joinpath("assets","electrolyzer.jl"))
include(joinpath("assets","fuelcell.jl"))
include(joinpath("assets","heater.jl"))
include(joinpath("assets","grid.jl"))
include(joinpath("assets","solar.jl"))
include(joinpath("assets","demand.jl"))
export Microgrid, Demand, Solar, Liion, ThermalStorage, H2Tank, FuelCell, Electrolyzer, Heater, Grid, GlobalParameters
export Electricity, Heat, Hydrogen
export add!
# Scenarios
include(joinpath("scenarios","scenarios.jl"))
include(joinpath("scenarios","reduction.jl"))
include(joinpath("scenarios","generation.jl"))
include(joinpath("scenarios","utils.jl"))
export Scenarios
export ManualReducer, SAAReducer, MeanValueReducer, FeatureBasedReducer
export UnitRangeTransform, ZScoreTransform
export PCAReduction, StatsReduction
export KmedoidsClustering
export MarkovGenerator, AnticipativeGenerator
export reduce, generate
# Optimization utils
include(joinpath("optimization","utils.jl"))
export Expectation, CVaR, WorstCase
# Operation optimization
include(joinpath("optimization","controller","dummy.jl"))
include(joinpath("optimization","controller","anticipative.jl"))
include(joinpath("optimization","controller","rb.jl"))
include(joinpath("optimization","controller","olfc.jl"))
export Dummy, RBC, Anticipative, OLFC
export RBCOptions, AnticipativeOptions, OLFCOptions
export initialize_controller!
# Investment optimization
include(joinpath("optimization","designer","manual.jl"))
include(joinpath("optimization","designer","milp.jl"))
include(joinpath("optimization","designer","metaheuristic.jl"))
export Manual, Metaheuristic, MILP
export MetaheuristicOptions, MILPOptions
export initialize_designer!
# Simulation
include(joinpath("simulation","informations.jl"))
include(joinpath("simulation","dynamics.jl"))
include(joinpath("simulation","power_balances.jl"))
include(joinpath("simulation","simulations.jl"))
export simulate!
# Utils
include(joinpath("utils","metrics.jl"))
include(joinpath("utils","plots.jl"))
include(joinpath("utils","saves.jl"))
export Metrics
export plot_operation, plot_metrics

end


text:
#include("benchmark.jl")

module FlatArraysTests

using Test, FlatArrays, Random
using Base:
    IndexStyle,
    has_offset_axes

for d in 1:3,
    (f, g) in ((firstindex, :firstindex),
               (lastindex,  :lastindex),
               (size,       :size),
               (axes,       :axis),
               (stride,     :stride))
    @eval $(Symbol(g,d))(A) = $f(A, $d)
end

#function memcpy!(A::DenseArray{T}, i::Integer,
#                 B::DenseArray{T}, j::Integer, n::Integer) where {T}
#    @assert firstindex(A) ≤ i ≤ i + n - 1 ≤ lastindex(A)
#    @assert firstindex(B) ≤ j ≤ j + n - 1 ≤ lastindex(B)
#    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
#          A, B, n*sizeof(T))
#end

function memcpy!(A::DenseArray{T}, B::DenseArray{T}, n::Integer) where {T}
    @assert firstindex(A) == 1 && 0 ≤ n ≤ lastindex(A)
    @assert firstindex(B) == 1 && 0 ≤ n ≤ lastindex(B)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
          A, B, n*sizeof(T))
    return A
end

function samevalues(A::AbstractArray{T}, B::AbstractArray{T},
                    exhaustive::Bool=false) where {T}
    @assert length(A) == length(B)
    j = 0
    for a in A
        j += 1
        B[j] == a || return false
    end
    if exhaustive
        i = 0
        for b in B
            i += 1
            A[i] == b || return false
        end
    end
    return true
end

@testset "FlatArrays" begin
    dims = (2,3,4,5)
    A = randn(dims)
    nelem = length(A)
    nrows = prod(dims[1:2])
    ncols = prod(dims[3:end])
    V = flatten(A, nelem)
    M = flatten(A, nrows, ncols)
    Q = ((firstindex,      1),
         (firstindex1,     1),
         (lastindex,       nelem),
         (lastindex1,      nelem),
         (length,          nelem),
         (ndims,           1),
         (eltype,          eltype(A)),
         (size,           (nelem,)),
         (size1,           nelem),
         (axes,           (1:nelem,)),
         (axis1,           1:nelem),
         (strides,        (1,)),
         (stride1,         1),
         (stride2,         nelem),
         (isflat,          true),
         (has_offset_axes, false),
         (IndexStyle,      IndexLinear()))
    for i in randperm(length(Q))
        (f, r) = Q[i]
        @test f(V) == r
    end
    Q = ((firstindex,      1),
         (firstindex1,     1),
         (firstindex2,     1),
         (lastindex,       nelem),
         (lastindex1,      nrows),
         (lastindex2,      ncols),
         (length,          nelem),
         (ndims,           2),
         (eltype,          eltype(A)),
         (size,           (nrows, ncols)),
         (size1,           nrows),
         (size2,           ncols),
         (axes,           (1:nrows, 1:ncols)),
         (axis1,           1:nrows),
         (axis2,           1:ncols),
         (strides,        (1, nrows)),
         (stride1,         1),
         (stride2,         nrows),
         (stride3,         nelem),
         (isflat,          true),
         (has_offset_axes, false),
         (IndexStyle,      IndexLinear()))
    for i in randperm(length(Q))
        (f, r) = Q[i]
        @test f(M) == r
    end
    @test isflat(A) == true
    @test isflat(A) == true
    V1 = view(A, :,  :, 1:3, 3)  # flat view
    @test isflat(V1) == true
    V2 = view(A, :,  :, 2:3, :) # non-flat view
    @test isflat(V2) == false
    @test reshape(A, nelem) == V
    @test reshape(A, nrows, ncols) == M
    B = similar(A, nrows, ncols)
    @test memcpy!(B, M, nelem) == reshape(A, nrows, ncols)
    @test samevalues(B, M, true)
    @test samevalues(B, V, true)
    @test samevalues(flatten(V1,length(V1)), V1, true)
    @test samevalues(flatten(V2,length(V2)), V2, true)
    @test_throws BoundsError V[-1]
    @test_throws BoundsError V[length(V)+1]
    @test_throws BoundsError M[1,0]
    @test_throws BoundsError M[nrows+1,0] == M[1,1]
end

end # module

max_stars_repo_name: emmt/FlatArrays

max_issues_repo_name: emmt/FlatArrays

max_forks_repo_name: emmt/FlatArrays


text:
# http://web.mit.edu/puzzle/www/2012/puzzles/william_s_bergman/behave/

@testset "behave" begin
    f1 = best_feature(["annieproulx", "commutative", "hugoweaving", "mountaindew", "mozambique", "sequoia"])
    @test f1.description == "has 5 unique vowels"

    f3 = best_feature(["almost", "biopsy", "chimp", "films", "ghost", "tux"])
    @test f3.description == "has at least 1 reverse alphabetical bigrams"

    f4 = best_feature(["balked", "barspoon", "highnoon", "klutzy", "onyx", "posted"])
    @test f4.description == "has at least 2 reverse sequential bigrams"

end


text:
export QBasedPolicy, get_prob

using Flux: softmax

"""
    QBasedPolicy(learner::Q, selector::S) -> QBasedPolicy{Q, S}

Use a Q-`learner` to generate the estimations of actions and use `selector` to get the action.
"""
Base.@kwdef struct QBasedPolicy{
    Q<:AbstractLearner,
    S<:AbstractActionSelector,
} <: AbstractPolicy
    learner::Q
    selector::S
end

(π::QBasedPolicy)(obs::Observation) = obs |> π.learner |> π.selector

get_prob(π::QBasedPolicy, s) = get_prob(π.selector, π.learner(s))

function update!(
    π::QBasedPolicy{<:Union{PrioritizedDQNLearner,RainbowLearner}},
    buffer::AbstractTurnBuffer,
)
    indexed_batch = extract_transitions(buffer, π)
    if !isnothing(indexed_batch)
        inds, batch = indexed_batch
        priorities = update!(π.learner, batch)
        isnothing(priorities) || (priority(buffer)[inds] .= priorities)
    end
end

update!(π::QBasedPolicy, model::AbstractEnvironmentModel; kw...) =
    update!(π.learner, model; kw...)

priority(transition, π::QBasedPolicy) = priority(transition, π.learner)

#####
# dispatches
#####

extract_transitions(buffer, π::QBasedPolicy) = extract_transitions(buffer, π.learner)

function extract_transitions(
    buffer::EpisodeTurnBuffer,
    π::QBasedPolicy{<:TDLearner{<:AbstractQApproximator,:ExpectedSARSA}},
)
    if length(buffer) > 0
        n = π.learner.n
        transitions = @views (
            states = state(buffer)[max(1, end - n - 1):end-1],
            actions = action(buffer)[max(1, end - n - 1):end-1],
            rewards = reward(buffer)[max(1, end - n):end],
            terminals = terminal(buffer)[max(1, end - n):end],
            next_states = state(buffer)[max(1, end - n):end],
        )
        s′ = transitions.next_states[end]
        merge(transitions, (probs_of_a′ = get_prob(π, s′),))
    else
        nothing
    end
end


text:
reload("SAMC")
#using SAMC
reload("MCBN")
#using MCBN

D = 8
template = zeros(Int,D,D)
cpds = MCBN.BayesNetDAI(D)

function factor_klds(bnd1::MCBN.BayesNetDAI, bnd2::MCBN.BayesNetDAI)
    0 #FIXME
end

function prior(bns::MCBN.BayesNetSampler)
    p_cpds = 0.2
    p_structural = 0.2

    #p_cpds::Float64
    #cpds::BayesNetDAI
    #p_structural::Float64
    #template::Matrix{Bool}
    #limparent::Int
    #limit number of parents
    d = size(bns.mat,1)
    tot = 0
    for i=1:d
        for j=1:d
            tot += bns.mat[i,j]
        end
        if tot > bns.limparent + 1
            return 1e10
        end
        tot = 0
    end

    adj_diff = sum(abs(bns.mat[bns.x,bns.x] - template))
    e_struct = p_structural * adj_diff

    # energy contribution from cpds
    cpd_diff = factor_klds(bns.bnd, cpds)
    e_cpd = p_cpds * cpd_diff
    return e_cpd + e_struct
end

preseed = rand(Uint)
srand(123)
#data = rand(1:2, D, 60)
gold = MCBN.random_net(D)
data = MCBN.draw_data(gold,10)
srand(preseed)

#bns = MCBN.BayesNetSampler(D, data, prior)
#samc1 = SAMC.set_energy_limits(bns, refden_power=0.0)
#samc1.thin = 1000
#samc1.stepscale = 20.0
#samc1.burn = 40_000
#@time SAMC.sample(samc1, 8_00_000, beta=0.6)

## MH
# bn = MCBN.BayesNetSampler(D, data, prior)
# mh = SAMC.MHRecord(bn, burn=burn, thin=thin)
# samptime = @elapsed SAMC.sample!(mh, iters)

## PopSAMC
# bngen = x->MCBN.BayesNetSampler(D, data, prior)
# samc2 = SAMC.set_energy_limits(bngen, 10, refden_power=0.0)
# samc2.thin = 100
# samc2.stepscale = 1_000.0
# samc2.burn = 10_000
cum_kld = Any[]

function edge_prob(bns::MCBN.BayesNetSampler)
    ord = sortperm(bns.x)
    return bns.mat[ord,ord]
end

kld_gold(bns::MCBN.BayesNetSampler) = MCBN.kld(gold, bns.bnd)

post_edge = SAMC.posterior_e(edge_prob,samc2)
post_kld = SAMC.posterior_e(kld_gold,samc2)


text:
module TestKernels
    using Base.Test
    using SmoothingKernels
    ks = SmoothingKernels.kernels

    for k in keys(ks)
        kernel = ks[k]
        for u in [0.0, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0]
            @test isapprox(kernel(-u), kernel(u))
        end
        integral, approxerr = quadgk(kernel, -Inf, Inf)
        @test abs(integral - 1.0) < approxerr
    end
end


text:
export VariableTimeModel

"""
    abstract type VariableTimeModel{I,J,O} <: Model{I,J,O}

A model that creates sfuncs based on the time delta between the parents and the current instance.
In general, the deltas can be different for different parents.
This type does not introduce any new functionality over Model.
Its purpose is to make explicit the fact that for this type of model, separate time deltas are possible.
Must implement 'make_initial', which takes the current time, and 'make_transition', which takes the current time and parent times.
"""
abstract type VariableTimeModel{I,J,O} <: Model{I,J,O} end


text:
function sparse_builder(fun::Function, ::Type{T}, M, N;
                        issymmetric=false, ishermitian=false) where T
    I = Int[]
    J = Int[]
    V = T[]

    set! = if issymmetric || ishermitian && T <: Real
        (i,j,v) -> begin
            iszero(v) && return
            push!(I, i)
            push!(J, j)
            push!(V, v)

            push!(I, j)
            push!(J, i)
            push!(V, v)
        end
    elseif ishermitian
        (i,j,v) -> begin
            iszero(v) && return
            push!(I, i)
            push!(J, j)
            push!(V, v)

            push!(I, j)
            push!(J, i)
            push!(V, conj(v))
        end
    else
        (i,j,v) -> begin
            iszero(v) && return
            push!(I, i)
            push!(J, j)
            push!(V, v)
        end
    end

    fun(set!)

    sparse(I, J, V, M, N)
end


text:
"""
    AbstractFockAddress

Supertype representing a Fock state.
"""
abstract type AbstractFockAddress{N,M} end

"""
    num_particles(::Type{<:AbstractFockAddress})
    num_particles(::AbstractFockAddress)

Number of particles represented by address.
"""
num_particles(a::AbstractFockAddress) = num_particles(typeof(a))
num_particles(::Type{<:AbstractFockAddress{N}}) where {N} = N

"""
    num_modes(::Type{<:AbstractFockAddress})
    num_modes(::AbstractFockAddress)

Number of modes represented by address.
"""
num_modes(a::AbstractFockAddress) = num_modes(typeof(a))
num_modes(::Type{<:AbstractFockAddress{<:Any,M}}) where {M} = M

"""
    num_components(::Type{<:AbstractFockAddress})
    num_components(::AbstractFockAddress)

Number of components in address.
"""
num_components(b::AbstractFockAddress) = num_components(typeof(b))

"""
    SingleComponentFockAddress{M}

A type representing a single component Fock state with `M` modes.

Implemented subtypes: [`BoseFS`](@ref), [`FermiFS`](@ref).

# Supported functionality

* [`find_mode`](@ref)
* [`find_occupied_mode`](@ref)
* [`is_occupied`](@ref)
* [`num_occupied_modes`](@ref)
* [`occupied_modes`](@ref): Lazy iterator.
* [`OccupiedModeMap`](@ref): `AbstractVector` with eager construction.
* [`excitation`](@ref): Create a new address.
* [`BoseFSIndex`](@ref) and [`FermiFSIndex`](@ref) for indexing.
"""
abstract type SingleComponentFockAddress{N,M} <: AbstractFockAddress{N,M} end

num_components(::Type{<:SingleComponentFockAddress}) = 1

"""
    find_mode(::SingleComponentFockAddress, i)

Find the `i`-th mode in address. Returns [`BoseFSIndex`](@ref) for [`BoseFS`](@ref), and
[`FermiFSIndex`](@ref) for [`FermiFS`](@ref). Can work on a tuple of modes. Does not check
bounds.

```jldoctest
julia> find_mode(BoseFS((1, 0, 2)), 2)
BoseFSIndex(occnum=0, mode=2, offset=2)

julia> find_mode(FermiFS((1, 1, 1, 0)), (2,3))
(FermiFSIndex(occnum=1, mode=2), FermiFSIndex(occnum=1, mode=3))
```

See [`SingleComponentFockAddress`](@ref).
"""
find_mode

"""
    find_occupied_mode(::SingleComponentFockAddress, k)
    find_occupied_mode(::BoseFS, k, [n])

Find the `k`-th occupied mode in address (with at least `n` particles).
Returns [`BoseFSIndex`](@ref) for [`BoseFS`](@ref), and [`FermiFSIndex`](@ref) for
[`FermiFS`](@ref). When unsuccessful it returns a zero index.

# Example

```jldoctest
julia> find_occupied_mode(FermiFS((1, 1, 1, 0)), 2)
FermiFSIndex(occnum=1, mode=2)

julia> find_occupied_mode(BoseFS((1, 0, 2)), 1)
BoseFSIndex(occnum=1, mode=1, offset=0)

julia> find_occupied_mode(BoseFS((1, 0, 2)), 1, 2)
BoseFSIndex(occnum=2, mode=3, offset=3)
```

See also [`occupied_modes`](@ref), [`OccupiedModeMap`](@ref),
[`SingleComponentFockAddress`](@ref).
"""
find_occupied_mode

"""
    is_occupied(::SingleComponentFockAddress, i)

Return `true` if index `i` points to an occupied mode.

See [`SingleComponentFockAddress`](@ref).
"""
is_occupied

"""
    num_occupied_modes(::SingleComponentFockAddress)

Get the number of occupied modes in address. Equivalent to
`length(`[`occupied_modes`](@ref)`(address))`, or the number of non-zeros in its ONR
representation.

# Example

```jldoctest
julia> num_occupied_modes(BoseFS((1, 0, 2)))
2
julia> num_occupied_modes(FermiFS((1, 1, 1, 0)))
3
```

See [`SingleComponentFockAddress`](@ref).
"""
num_occupied_modes

"""
    occupied_modes(::SingleComponentFockAddress)

Return a lazy iterator over all occupied modes in an address. Iterates over
[`BoseFSIndex`](@ref)s for [`BoseFS`](@ref), and over [`FermiFSIndex`](@ref)s for
[`FermiFS`](@ref). See [`OccupiedModeMap`](@ref) for an eager version.

# Example

```jldoctest
julia> b = BoseFS((1,5,0,4));

julia> foreach(println, occupied_modes(b))
BoseFSIndex(occnum=1, mode=1, offset=0)
BoseFSIndex(occnum=5, mode=2, offset=2)
BoseFSIndex(occnum=4, mode=4, offset=9)
```

```jldoctest
julia> f = FermiFS((1,1,0,1,0,0,1));

julia> foreach(println, occupied_modes(f))
FermiFSIndex(occnum=1, mode=1)
FermiFSIndex(occnum=1, mode=2)
FermiFSIndex(occnum=1, mode=4)
FermiFSIndex(occnum=1, mode=7)
```
See also [`find_occupied_mode`](@ref),
[`SingleComponentFockAddress`](@ref).
"""
occupied_modes

"""
    excitation(a::SingleComponentFockAddress, creations::NTuple{N}, destructions::NTuple{N})

Generate an excitation on address `a` by applying `creations` and `destructions`, which are
tuples of the appropriate address indices (i.e. [`BoseFSIndex`](@ref) for bosons, or
[`FermiFSIndex`](@ref) for fermions).

```math
a^†_{c_1} a^†_{c_2} \\ldots a_{d_1} a_{d_2} \\ldots |\\mathrm{a}\\rangle \\to
α|\\mathrm{nadd}\\rangle
```

Returns the new address `nadd` and the value `α`. If the excitation is illegal, returns an
arbitrary address and the value `0.0`.

# Example

```jldoctest
julia> f = FermiFS((1,1,0,0,1,1,1,1))
FermiFS{6,8}((1, 1, 0, 0, 1, 1, 1, 1))

julia> i, j, k, l = find_mode(f, (3,4,2,5))
(FermiFSIndex(occnum=0, mode=3), FermiFSIndex(occnum=0, mode=4), FermiFSIndex(occnum=1, mode=2), FermiFSIndex(occnum=1, mode=5))

julia> excitation(f, (i,j), (k,l))
(FermiFS{6,8}((1, 0, 1, 1, 0, 1, 1, 1)), -1.0)
```

See [`SingleComponentFockAddress`](@ref).
"""
excitation

"""
    OccupiedModeMap(add) <: AbstractVector

Get a map of occupied modes in address as an `AbstractVector` of indices compatible with
[`excitation`](@ref) - [`BoseFSIndex`](@ref) or [`FermiFSIndex`](@ref).

`OccupiedModeMap(add)[i]` contains the index for the `i`-th occupied mode.
This is useful because repeatedly looking for occupied modes with
[`find_occupied_mode`](@ref) can be time-consuming.
`OccupiedModeMap(add)` is an eager version of the iterator returned by
[`occupied_modes`](@ref). It is similar to [`onr`](@ref) but contains more information.

# Example

```jldoctest
julia> b = BoseFS((10, 0, 0, 0, 2, 0, 1))
BoseFS{13,7}((10, 0, 0, 0, 2, 0, 1))

julia> mb = OccupiedModeMap(b)
3-element OccupiedModeMap{7, BoseFSIndex}:
 BoseFSIndex(occnum=10, mode=1, offset=0)
 BoseFSIndex(occnum=2, mode=5, offset=14)
 BoseFSIndex(occnum=1, mode=7, offset=18)

julia> f = FermiFS((1,1,1,1,0,0,1,0,0))
FermiFS{5,9}((1, 1, 1, 1, 0, 0, 1, 0, 0))

julia> mf = OccupiedModeMap(f)
5-element OccupiedModeMap{5, FermiFSIndex}:
 FermiFSIndex(occnum=1, mode=1)
 FermiFSIndex(occnum=1, mode=2)
 FermiFSIndex(occnum=1, mode=3)
 FermiFSIndex(occnum=1, mode=4)
 FermiFSIndex(occnum=1, mode=7)

julia> mf == collect(occupied_modes(f))
true

julia> dot(mf, mb)
11

julia> dot(mf, 1:20)
17
```
See also [`dot`](@ref), [`SingleComponentFockAddress`](@ref).
"""
struct OccupiedModeMap{N,T} <: AbstractVector{T}
    indices::SVector{N,T} # N = min(N, M)
    length::Int
end

function OccupiedModeMap(add::SingleComponentFockAddress{N,M}) where {N,M}
    modes = occupied_modes(add)
    T = eltype(modes)
    # There are at most N occupied modes. This could be also @generated for cases where N ≫ M
    indices = MVector{min(N,M),T}(undef)
    i = 0
    for index in modes
        i += 1
        @inbounds indices[i] = index
    end
    return OccupiedModeMap(SVector(indices), i)
end

Base.size(om::OccupiedModeMap) = (om.length,)
function Base.getindex(om::OccupiedModeMap, i)
    @boundscheck 1 ≤ i ≤ om.length || throw(BoundsError(om, i))
    return om.indices[i]
end

"""
    dot(map::OccupiedModeMap, vec::AbstractVector)
    dot(map1::OccupiedModeMap, map2::OccupiedModeMap)
Dot product extracting mode occupation numbers from an [`OccupiedModeMap`](@ref) similar
to [`onr`](@ref).

```jldoctest
julia> b = BoseFS((10, 0, 0, 0, 2, 0, 1))
BoseFS{13,7}((10, 0, 0, 0, 2, 0, 1))

julia> mb = OccupiedModeMap(b)
3-element OccupiedModeMap{7, BoseFSIndex}:
 BoseFSIndex(occnum=10, mode=1, offset=0)
 BoseFSIndex(occnum=2, mode=5, offset=14)
 BoseFSIndex(occnum=1, mode=7, offset=18)

julia> dot(mb, 1:7)
27

julia> mb⋅(1:7) == onr(b)⋅(1:7)
true
```
See also [`SingleComponentFockAddress`](@ref).
"""
function LinearAlgebra.dot(map::OccupiedModeMap, vec::AbstractVector)
    value = zero(eltype(vec))
    for index in map
        value += vec[index.mode] * index.occnum
    end
    return value
end
LinearAlgebra.dot(vec::AbstractVector, map::OccupiedModeMap) = dot(map, vec)

# Defined for consistency. Could also be used to compute cross-component interactions in
# real space.
function LinearAlgebra.dot(map1::OccupiedModeMap, map2::OccupiedModeMap)
    i = j = 1
    value = 0
    while i ≤ length(map1) && j ≤ length(map2)
        index1 = map1[i]
        index2 = map2[j]
        if index1.mode == index2.mode
            value += index1.occnum * index2.occnum
            i += 1
            j += 1
        elseif index1.mode < index2.mode
            i += 1
        else
            j += 1
        end
    end
    return value
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Type{Schema},Vector{Symbol},Vector{DataType}})   # time: 0.0035495
    Base.precompile(Tuple{typeof(which(_columntable,(Tables.Schema{names, types},Any,)).generator.gen),Any,Any,Any,Any,Any})   # time: 0.0033783
end


text:
struct Normal{T}
    μ::T
    σ::T
    function Normal(m,s)
        μ,σ = promote(m,s)
        T = typeof(μ)
        new{T}(μ, σ)
    end
end
Normal(;μ,σ) = Normal(μ,σ)
Base.show(io::IO, o::Normal) = print(io, "Normal(μ=$(o.μ), σ=$(o.σ))")

struct NormalRegressor end

function (o::NormalRegressor)(arr::AbstractMatrix)
    @argcheck mod(size(arr,1), 2) == 0
    nc, nb = size(arr)
    b = Int(nc/2)
    cmu    = view(arr, 1:b)
    csigma = view(arr, b+1:2b)
    return Normal_from_2channels.(cmu, csigma)
end
function Normal_from_2channels(cmu, csigma)
    return Normal(cmu, softplus(csigma))
end

function Distributions.logpdf(o::Normal, x)
    μ = o.μ
    σ = o.σ
    T = typeof(μ)
    half = T(1) / T(2)
    B = -1/(2*σ^2) * (μ - x)^2
    A = -half * log2π - log(σ)
    return A + B
end

sampletype(::Type{Normal{T}}) where {T} = T
Distributions.mean(o::Normal) = o.μ
Distributions.std(o::Normal)  = o.σ
Distributions.Normal(o::Normal) = Distributions.Normal(o.μ, o.σ)


text:
"""
    Cell{T}

preserve 'fmtId' for Cell
currently there is no functionality for fmtId after DataType of cell is determined
"""
struct Cell{T}
    numFmtId::UInt16
    v::T
end
Cell(id, v::T) where T = Cell{T}(id, v)

function Cell(wb::WorkBook, x::EzXML.Node)
    # @show "cell data / " x["r"]
    fmt_id, T = find_cell_datatype(wb, x)

    cell_values = filter(iselement, nodes(x))
    # sometimes there can be no value, even through it has cell address
    value = if isempty(cell_values)
                missing
            else # Excel always store value `v` node at the end
                v = nodecontent.(cell_values)[end]
                v = get_cellvalue(T, wb, v)
            end
    Cell(fmt_id, value)
end

# fallback methods
Base.Symbol(x::Cell) = convert(Symbol, x)
Base.convert(::Type{Symbol}, x::Cell) = Symbol(x.v)

function find_cell_datatype(wb::WorkBook, x)
    fmt_id = haskey(x, "s") ? wb.xf_index_to_number_format[parse(Int, x["s"])+1] : 0
    T = if haskey(x, "t")
            ECMA_CELLTYPE[x["t"]]
        else
            fmt = haskey(ECMA_NUMBER_FORMAT, fmt_id) ? ECMA_NUMBER_FORMAT[fmt_id] :
                  wb.custom_number_format[fmt_id]
            fmt[1]
        end
    (fmt_id, T)
end

function get_cellvalue(::Type{Bool}, wb, x)::Bool
    if x == "1" || x == "0"
        x = parse(Int8, x)
    end
    Bool(x)
end
function get_cellvalue(::Type{DateTime}, wb, x)
    x
end
get_cellvalue(::Type{ExcelDateTime}, wb, x) = parse(ExcelDateTime, x)
get_cellvalue(::Type{Missing}, wb, x) = missing
get_cellvalue(::Type{String}, wb, x) = String(x)
get_cellvalue(::Type{Float64}, wb, x) = parse(Float64, x)
get_cellvalue(::Type{T}, wb, x) where T<:Integer= parse(T, x)
function get_cellvalue(::Type{SharedString}, wb, x)
    wb.shared_strings[parse(Int, x)+1]
end
get_cellvalue(::Type{GeneralFormat}, wb, x) = parse(GeneralFormat, x)
get_cellvalue(::Type{ScientificFormat}, wb, x) = parse(ScientificFormat, x)

# peel off 'Cell{T}' and returns T(value)
value(x::Cell) = x.v
value(x::Missing) = missing
function peeloff(ws::WorkSheet)
    x = ws.data
    isa(x, Array{Union{EzExcel.Cell, Missing}}) ? _peeloff(x) : x
end
function _peeloff(datas)
    v = value.(datas)
    T = unique(typeof.(v))
    if length(T) < 3
        convert(Array{Union{T...}}, v)
    else
        convert(Array{Any}, v)
    end
end


text:
module child

f1() = "one"
# f2() = 2 # method deleted
f3() = 3

end


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

JLLWrappers.@generate_main_file_header("SPIRV_LLVM_Translator")
JLLWrappers.@generate_main_file("SPIRV_LLVM_Translator", UUID("4a5d46fc-d8cf-5151-a261-86b458210efb"))
end  # module SPIRV_LLVM_Translator_jll


text:
#
# Formatting Julia expressions.
#

using PrettyPrinting:
    Layout,
    indent,
    list_layout,
    literal,
    pair_layout,
    pprint,
    tile

print_expr(io::IO, ex) =
    pprint(io, tile_expr(ex))

quoteof(obj) =
    obj

quoteof_inner(obj) =
    quoteof(obj)

quoteof(ref::Base.RefValue) =
    quoteof(ref.x)

quoteof(s::Symbol) =
    QuoteNode(s)

function quoteof(f::Function)
    s = nameof(f)
    if startswith(string(s), "#")
        ex = _reconstruct(f)
        if ex !== nothing
            s = ex
        end
    end
    s
end

quoteof(f::Union{Function,Type}, args::Vector{Any}) =
    Expr(:call, nameof(f), quoteof.(args)...)

quoteof(::typeof(broadcast), args::Vector{Any}) =
    if length(args) >= 1 && args[1] isa Function
        quoteof(broadcast, args[1], args[2:end])
    else
        Expr(:call, nameof(broadcast), quoteof.(args)...)
    end

function quoteof(::typeof(broadcast), f::Function, args::Vector{Any})
    ex = quoteof(f)
    if Meta.isexpr(ex, :(->), 2)
        names =
            if Meta.isexpr(ex.args[1], :tuple)
                ex.args[1].args
            else
                [ex.args[1]]
            end
        repl = Dict(zip(names, args))
        return _broadcast(ex.args[2], repl)
    end
    if ex isa Symbol && Base.operator_precedence(ex) > 0
        return Expr(:call, Symbol(".$ex"), args...)
    else
        return Expr(:(.), ex, Expr(:tuple, args...))
    end
end

function _broadcast(ex, repl)
    if ex isa Symbol
        if ex in keys(repl)
            ex = repl[ex]
        end
    elseif ex isa Expr
        ex = Expr(ex.head, (_broadcast(arg, repl) for arg in ex.args)...)
        if ex.head == :call && length(ex.args) >= 1 && ex.args[1] isa Symbol
            func = ex.args[1]
            ex =
                if Base.operator_precedence(func) > 0
                    func = Symbol(".$func")
                    Expr(ex.head, func, ex.args[2:end]...)
                else
                    Expr(:(.), func, Expr(:tuple, ex.args[2:end]...))
                end
        end
    end
    ex
end

quoteof(v::Vector) =
    Expr(:vect, quoteof.(v)...)

quoteof(p::Pair) =
    Expr(:call, :(=>), quoteof(p.first), quoteof(p.second))

tile_expr(obj; precedence=0) =
    tile(obj)

tile_expr(sym::Symbol; precedence=0) =
    literal(sym)

tile_expr(qn::QuoteNode; precedence=0) =
    literal(string(qn))

function tile_expr(ex::Expr; precedence=0)
    if ex.head == :call
        func = ex.args[1]
        if func isa Function
            func = nameof(func)
        end
        args = _flatten(func, ex.args[2:end])
        precedence′ = Base.operator_precedence(func)
        arg_lts = Layout[tile_expr(arg, precedence=precedence′) for arg in args]
        if func == :(=>) && length(arg_lts) == 2
            key_lt, val_lt = arg_lts
            pair_layout(key_lt, val_lt)
        elseif precedence′ > 0
            sep = func == :(:) ? "$func" : " $func "
            par =
                if precedence′ < precedence
                    ("(", ")")
                else
                    ("", "")
                end
            if length(arg_lts) == 2
                literal(par[1]) *
                pair_layout(arg_lts..., sep=sep, tab=0) *
                literal(par[2])
            else
                list_layout(arg_lts, par=par, sep=sep)
            end
        else
            par = ("$func(", ")")
            list_layout(arg_lts, par=par)
        end
    elseif ex.head == :tuple
        list_layout(Layout[tile_expr(arg) for arg in ex.args])
    elseif ex.head == :vect
        list_layout(Layout[tile_expr(arg) for arg in ex.args], par=("[", "]"))
    elseif ex.head == :ref && length(ex.args) >= 1
        tile_expr(ex.args[1]) *
        list_layout(Layout[tile_expr(arg) for arg in ex.args[2:end]], par=("[", "]"))
    elseif ex.head in (:(=), :(->), :(<:)) && length(ex.args) == 2
        ilt = tile_expr(ex.args[1])
        olt = tile_expr(ex.args[2])
        pair_layout(ilt, olt, sep=(ex.head == :(<:) ? "$(ex.head)" : " $(ex.head) "))
    elseif ex.head == :(...) && length(ex.args) == 1
        tile_expr(ex.args[1]) * literal("...")
    else
        literal(string(ex))
    end
end

function _flatten(func, args)
    if !(func === :(>>) || func === :(|>))
        return args
    end
    args′ = []
    for arg in args
        if Meta.isexpr(arg, :call) && arg.args[1] == func
            append!(args′, _flatten(func, arg.args[2:end]))
        else
            push!(args′, arg)
        end
    end
    args′
end

function _reconstruct(f::Function)
    infos = code_lowered(f)
    if length(infos) != 1
        return
    end
    info = infos[1]
    ssa = []
    for line in info.code
        if line === nothing
            continue
        end
        if Meta.isexpr(line, :call)
            ex = _reconstruct(line, ssa, info.slotnames)
            push!(ssa, ex)
        elseif Meta.isexpr(line, :return, 1)
            body = _reconstruct(line.args[1], ssa, info.slotnames)
            if body === nothing
                return
            end
            args = [(startswith(string(n), "#") ? Symbol("_$i") : n) for (i, n) in enumerate(info.slotnames[2:end])]
            return Expr(:(->), length(args) == 1 ? args[1] : Expr(:tuple, args...), body)
        else
            return
        end
    end
    return
end

function _reconstruct(ex, ssa, slots)
    if ex isa Core.SSAValue && checkbounds(Bool, ssa, ex.id)
        ex = ssa[ex.id]
    elseif ex isa Core.SlotNumber && checkbounds(Bool, slots, ex.id)
        s = slots[ex.id]
        if startswith(string(s), "#")
            s = Symbol("_$(ex.id-1)")
        end
        ex = s
    elseif ex isa GlobalRef
        ex = ex.name
    elseif ex isa Expr
        ex = Expr(ex.head, (_reconstruct(arg, ssa, slots) for arg in ex.args)...)
    end
    return ex
end


text:
"""
The base type used by all breeder sources.
"""
abstract BreederSource

"""
The base type used by all breeder source definitions.
"""
abstract BreederSourceDefinition

"""
Provides a definition of a selection breeder source.
"""
type SelectionBreederSourceDefinition <: BreederSourceDefinition
  label::AbstractString
  operator::SelectionDefinition

  SelectionBreederSourceDefinition(label::AbstractString, op::SelectionDefinition) =
    new(label, op)
end

type SelectionBreederSource <: BreederSource
  label::AbstractString
  operator::Selection
  eigen::Type

  SelectionBreederSource(label::AbstractString, s::Selection) =
    new(label, s)
end

"""
Composes a given selection breeder source.
"""
function compose!(
  def::SelectionBreederSourceDefinition,
  ::Species,
  ::Dict{AbstractString, BreederSource}
)
  op = compose!(def.operator)
  s = SelectionBreederSource(def.label, op)
  s.eigen = anonymous_type(breeder)
  s
end

"""
A selection source.
"""
selection_source(label::AbstractString, op::SelectionDefinition) =
  SelectionBreederSourceDefinition(label, op)

"""
Provides a definition of a variation breeder source.
"""
type VariationBreederSourceDefinition <: BreederSourceDefinition
  label::AbstractString
  source::AbstractString
  stage::AbstractString
  operator::VariationDefinition

  VariationBreederSourceDefinition(
    label::AbstractString,
    source::AbstractString,
    stage::AbstractString,
    operator::VariationDefinition
  ) =
    new(label, source, stage, operator)
end

type VariationBreederSource <: BreederSource
  label::AbstractString
  operator::Variation
  source_name::AbstractString
  stage_name::AbstractString
  eigen::Type
  source::BreederSource
  stage_getter::Function

  VariationBreederSource(
    label::AbstractString,
    v::Variation,
    source_name::AbstractString,
    stage_name::AbstractString
  ) =
    new(label, v, source_name, stage_name)
end

"""
Composes a variation breeder source.

**Parameters:**

* `def::VariationBreederSourceDefinition`, a definition of the breeder source
that is to be composed.
* `species::Species`, the species that this breeder belongs to.
* `sources::Dict{AbstractString, Source}`, a dictionary containing the sources
of the breeder that this source belongs to.
"""
function compose!(
  def::VariationBreederSourceDefinition,
  sp::Species,
  sources::Dict{AbstractString, BreederSource}
)
  # Find the stage and representation used by this source.
  stage = sp.stages[def.stage]
  representation = species.rep(stage) # TODO: FIX FIX FIX

  # Build the variation operator.
  op = compose!(def.operator, representation)

  # Build the breeder source.
  v = VariationBreederSource(def.label, op, def.source, def.stage)
  v.eigen = anonymous_type(breeder)
  v.source = sources[def.source]
  v.stage_getter =
    eval(Main, Base.parse("inds -> Wallace.IndividualStage{$(chromosome(representation))}[i.$(def.stage) for i in inds]"))
  v
end

"""
DOCUMENT: breeder.variation
"""
variation_source(
  label::AbstractString,
  source::AbstractString,
  stage::AbstractString,
  op::VariationDefinition
) =
  VariationBreederSourceDefinition(label, source, stage, op)

type MultiBreederSource <: BreederSource
  sources::Vector{BreederSource}
  MultiBreederSource(sources::Vector{BreederSource}) = new(sources)
end


text:
using Documenter, QuantEcon

makedocs(
    modules = [QuantEcon],
    format = :html,
    sitename = "QuantEcon.jl",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        # "User Guide" => "man/guide.md",
        "API" => Any[
            "QuantEcon" => "api/QuantEcon.md"
        ],
        "Contributing" => "man/contributing.md"
    ]
)

deploydocs(
    repo = "github.com/QuantEcon/QuantEcon.jl.git",
    branch = "gh-pages",
    target = "build",
    julia  = "0.6",
    make = nothing,
)


text:
#TODO: I feel a bit of a code smell for having to constantly repeat the name in calls like
# style=Style(...). Is there a better way is something like style=(:stroke=NC"red", ...)
# this could be a short-of-shorthand for the full naming? But then I introduce the many-ways
# to do the same thing.
type Style
    stroke::Union{RGB, Void}
    stroke_width::Float64
    fill::Union{RGB, Void}
    fill_opacity::Float64 # should check for being between 0-1
    function Style(;stroke=NC"black", stroke_width=1.0, fill=NC"white", fill_opacity=1.0)
        new(stroke, stroke_width, fill, fill_opacity)
    end
end

# Matplotlib requires the color settings using 'none' and not the python type None, which
# feels very odd. So a simple utility is needed when converting the color options to the
# Dict passed to python. Where nothing is set return :none for the correct conversion.

# This function is not type stable, but neither is the Color fields do to needed to accept
# nothing
function _colorkw(field::Union{RGB, Void})
    if field == nothing
        return :none
    else
        return field
    end
end

# keywords for matplotlib patches are uniform
function style2kw(s::Style)
    Dict(:edgecolor => _colorkw(s.stroke), :linewidth => s.stroke_width, :facecolor => _colorkw(s.fill), :alpha => s.fill_opacity)
end

# keyworks for Line2D are different than patches
function style2kw_l(s::Style)
    Dict(:color => _colorkw(s.stroke), :linewidth => s.stroke_width)
end

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
export dropout
export dropout!
export xdropout
export xdropout!

"""
    dropout(x::Variable{T}; p=0.1) -> y::Variable{T}

Randomly zeros some elements of the input tensor `x` with probability `p` using samples
from a Bernoulli distribution. This is an effective way to regularization and preventing
the co-adaptation of neurons. The output elements of `y` are scaled by a factor of `1/(1-p)`
during training. `dropout` should be removed at evaluation. Dropout is also viewed as a
mean of data augmentation.
"""
function dropout(x::Variable{T}; p=0.1) where T
    @assert 0.0<=p<1.0 "p is in [0,1), but got p=$p"
    τ = eltype(T)
    l = τ(1)
    p = τ(p)
    m = T(rand(τ, x.shape) .< (l - p)) .* (l/(l - p)) # weighted mask
    y = Variable{T}(ᵛ(x) .* m, x.backprop)
    if x.backprop
        function dropoutBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* m
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, dropoutBackward)
    end
    return y
end

"""
    dropout!(x::Variable{T}; p=0.1) -> y::Variable{T}

Randomly zeros some elements of the input tensor `x` with probability `p` using samples
from a Bernoulli distribution. This is an effective way to regularization and preventing
the co-adaptation of neurons. The output elements of `y` are scaled by a factor of `1/(1-p)`
during training. `dropout!` should be removed at evaluation. Dropout is also viewed as a
mean of data augmentation.
"""
function dropout!(x::Variable{T}; p=0.1) where T
    @assert 0.0<=p<1.0 "p is in [0,1), but got p=$p"
    τ = eltype(T)
    l = τ(1)
    p = τ(p)
    m = T(rand(τ, x.shape) .< (l - p)) .* (l/(l - p)) # weighted mask
    y = Variable{T}(dotmul!(ᵛ(x), m), x.backprop)
    if x.backprop
        function dropoutBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* m
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, dropoutBackward)
    end
    return y
end

"""
    xdropout(x::Variable; p=0.1, dim=1) -> y::Variable{T}

Randomly zeros some `slices` of the input tensor `x` with probability `p` using samples
from a Bernoulli distribution. This is an effective way to regularization and preventing
the co-adaptation of neurons. The output elements of `y` are scaled by a factor of `1/(1-p)`
during training. `xdropout` should be removed at evaluation. Dropout is also viewed as a
mean of data augmentation.
"""
function xdropout(x::Variable{T}; p=0.1, dim=1) where T
    @assert 0.0<=p<1.0 "p is in [0,1), but got p=$p"
    τ = eltype(T)
    l = τ(1)
    p = τ(p)
    N = size(x, dim)
    S = ntuple(i -> i≠dim ? 1 : N, ndims(x))
    m = T(rand(τ, S) .< (l - p)) .* (l/(l - p))  # weighted mask
    y = Variable{T}(ᵛ(x) .* m, x.backprop)
    if x.backprop
        function dropoutBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* m
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, dropoutBackward)
    end
    return y
end

"""
    xdropout!(x::Variable; p=0.1, dim=1) -> y::Variable{T}

Randomly zeros some `slices` of the input tensor `x` with probability `p` using samples
from a Bernoulli distribution. This is an effective way to regularization and preventing
the co-adaptation of neurons. The output elements of `y` are scaled by a factor of `1/(1-p)`
during training. `xdropout!` should be removed at evaluation. Dropout is also viewed as a
mean of data augmentation.
"""
function xdropout!(x::Variable{T}; p=0.1, dim=1) where T
    @assert 0.0<=p<1.0 "p is in [0,1), but got p=$p"
    τ = eltype(T)
    l = τ(1)
    p = τ(p)
    N = size(x, dim)
    S = ntuple(i -> i≠dim ? 1 : N, ndims(x))
    m = T(rand(τ, S) .< (l - p)) .* (l/(l - p))  # weighted mask
    y = Variable{T}(dotmul!(ᵛ(x), m), x.backprop)
    if x.backprop
        function dropoutBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* m
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, dropoutBackward)
    end
    return y
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
  # Vinicius Piccoli  - PhD student at M3G - UNICAMP

  #= Informations:

  This simple code reads protein pdb files and writes inputs for Packmol.

  The code only works for a simulation of a protein solvated by water and ionic liquids (e.g EMIMDCA, EMIMBF4, BMIMNO3, BMIMCl, etc.)

  Example:
  gmx pdb2gmx -f protein.pdb (crystallographic structure) -o processed.pdb
  measure_prot(processed.pdb,10) - this will return the dimensions of the protein in x,y and z axis plus 10 Angstrons.

  The conc function write an “inp” file (which will be the input for Packmol and, therefore, contains the information for the simulation box)
  and the topology file for the simulation

  conc(MM::Float64,Nome::String,Nomean::String, conc::Float64)
    MM - Ionic Liquid molar mass
    Nome - Name of the cation
    Nomean - Name of the anion
    conc - Ionic liquid concentration
  =#

  function measure_prot(nome,dist)

  ## Variables - Measurement of the protein(X, Y and Z axis)
    pdbfile   = nome
    segment   = "all"
    restype   = "all"
    atomtype  = "all"
    firstatom = "all"
    lastatom  = "all"

  ## Open the file -
    file = open("$pdbfile","r+")
    natoms = 0
    var = 0

    for line in eachline(file)

      data = split(line)
      consider = false
      if data[1] == "ATOM" || data[1] == "HEATOM"
        natoms = natoms + 1
        ss = data[11]
        rt = data[4]
        at = data[3]
        consider = true
      else
        continue
      end

      if firstatom != "all"
        if natoms < firstatom
          consider = false
        end
      end

      if lastatom != "all"
        if natoms > lastatom
          consider = false
        end
      end

      if segment != "all"
        if ss != segment
          consider = false
        end
      end

      if restype != "all"
        if rt != restype
          consider = false
        end
      end

      if atomtype != "all"
        if at != atomtype
          consider = false
        end
      end

      if consider == true
        global x = parse(Float64,data[7])
        global y = parse(Float64,data[8])
        global z = parse(Float64,data[9])
      end

      if var == 1
        if x < xmin
            xmin = x
        end
        if y < ymin
            ymin = y
        end
        if z < zmin
            zmin = z
        end
        if x > xmax
            xmax = x
        end
        if y > ymax
            ymax = y
        end
        if z > zmax
            zmax = z
        end
      else

        global   xmin = x
        global   ymin = y
        global   zmin = z
        global   xmax = x
        global   ymax = y
        global   zmax = z
        var = 1

      end

    end

  ## Box size - add 30 A in each size of box
    bx = round(Int,((xmax - xmin) + dist ))
    by = round(Int,((ymax - ymin) + dist ))
    bz = round(Int,((zmax - zmin) + dist ))

    return bx , by, bz

  end

  # Function to calculate number of components to put inside a box given a specific concentration

  function conc(MM::Float64,Nome::String,Nomean::String, conc::Float64)

      # number of particles calculationa

      # volume of the box (L)
      vol_box(a,b,c) = a*b*c*1e-27;

      # protein volume (L)
      vol_prot(m) = (m/(6.02e23))*1e-3;

      # Water volume (L)
      vol_wat(n,mw) = n * (mw/(6.02e23))*1e-3;

      # Solution volume
      vol_sol(vc,vp) = vc - vp;

      # number of ionic liquids molecules
      num_il(vs,cil) = round(Int128,(vs*cil*6.02e23));

      # Volume of ionic liquids moleculres
      v_il(nil,mil) = (nil*mil*1e-3)/(6.02e23);

      # each ionic liquid will occupy the half of the volume

      # Number of water molecules
      num_wat(vs,vil) = round(Int128,((vs - vil) * 6.02e23) / (18*1e-3));

      # Box dimensions
      lx,ly,lz = measure_prot("ubq.pdb")   # Ubiquitin pdb

      # Total concentration of ions
      c = conc;

      vs   = vol_box(2*lx,2*ly,2*lz) - vol_prot(8560); # Solution volume / Ubiquitin molar mass = 8560
      nil  = round(Int,(num_il(vs,c)));                # Number of ionic liquid molecules for a concentration "c"

      nwat1 = num_wat(vs, v_il(nil,MM));               # Number of water molecules
      nwat = nwat1

      io = open("box.inp","w")

      println(io,"tolerance 2.0")
      println(io,"output system.pdb")
      println(io,"add_box_sides 1.0")
      println(io,"filetype pdb")
      println(io,"seed -1")
      println(io,"                  ")
      println(io,"structure ubq.pdb")
      println(io," number 1")
      println(io," center")
      println(io," fixed 0. 0. 0. 0. 0. 0.")
      println(io,"end structure")
      println(io,"                  ")
      println(io,"structure WATER.pdb")
      println(io," number $nwat")
      println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
      println(io,"end structure")
      println(io,"                  ")
      println(io,"structure $(Nome)_VSIL.pdb")
      println(io," number $nil")
      println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
      println(io,"end structure")
      println(io,"                  ")
      println(io,"structure $(Nomean)_VSIL.pdb")
      println(io," number $nil")
      println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
      println(io,"end structure")

      close(io)

      topa = open("topola.top","r+")

      # number of lines !
      nlines = 0
      for line in eachline(topa)
        nlines += 1
      end

      close(topa)

      Nomean2="NOTHING"
      Nome2="NOTHING"

      topa = open("topola.top","r+")
      top = open("topol.top","w")
      ntline = 0

      for line in eachline(topa)

        ntline += 1

        if occursin("AN1",line) && ntline != nlines && ntline < nlines - 4
          println(top,replace(line,"AN1" => "$Nomean",count = 1))
        elseif occursin("CA1",line) && ntline < nlines - 4
          println(top,replace(line,"CA1" => "$Nome",count = 1))
        elseif occursin("AN2",line) && ntline < nlines - 4
          println(top,replace(line,"AN2" => "$Nomean2",count = 1))
        elseif occursin("CA2",line) && ntline < nlines - 4
          println(top,replace(line,"CA2" => "$Nome2",count = 1))
        elseif ntline == (nlines - 4)
          println(top,replace(line,line => "SOL               $(nwat+58)",count = 1)) # 58 cristalization water molecules

        # data for the ions

        elseif ntline == (nlines - 1)
          if Nomean == "DCA"
            println(top,replace(line,line => "NC                $nil",count = 1))
          elseif Nomean == "BF4"
            println(top,replace(line,line => "BF4                $nil",count = 1))
          elseif Nomean == "NO3"
            println(top,replace(line,line => "NO3                $nil",count = 1))
          elseif Nomean == "Cl"
            println(top,replace(line,line => "Cl                 $nil",count = 1))
          elseif Nomean == "SCN"
            println(top,replace(line,line => "SCN                $nil",count = 1))
          end
        elseif ntline == (nlines - 3)
          println(top,replace(line,line => "$(Nome[1:3])                $nil ",count = 1))
        else
          println(top,line)
        end
      end
      close(top)
      close(topa)

  end


text:
function main()
    N = parse(Int, readline())
    As = parse.(Int, split(readline()))

    flg = true
    cnt = -1
    while flg
        for i in 1:length(As)
            if As[i] % 2 == 0
                As[i] = As[i] / 2
            else
                flg = false
                break
            end
        end
        cnt += 1
    end
    println(cnt)
end

main()


text:
# Julia wrapper for header: /usr/local/include/libavutil/camellia.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_camellia_alloc,
    av_camellia_init,
    av_camellia_crypt

function av_camellia_alloc()
    ccall((:av_camellia_alloc,libavutil),Ptr{AVCAMELLIA},())
end

function av_camellia_init(ctx,key,key_bits::Integer)
    ccall((:av_camellia_init,libavutil),Cint,(Ptr{AVCAMELLIA},Ptr{UInt8},Cint),ctx,key,key_bits)
end

function av_camellia_crypt(ctx,dst,src,count::Integer,iv,decrypt::Integer)
    ccall((:av_camellia_crypt,libavutil),Void,(Ptr{AVCAMELLIA},Ptr{UInt8},Ptr{UInt8},Cint,Ptr{UInt8},Cint),ctx,dst,src,count,iv,decrypt)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
struct Reaction
    product_stoich_coeffs::SparseMatrixCSC{Float64, Int64}
    reactant_stoich_coeffs::SparseMatrixCSC{Float64, Int64}
    reactant_orders::SparseMatrixCSC{Float64, Int64}
    is_reversible::Array{Bool, 1}
    Arrhenius_coeffs::Array{Float64, 2}
    Arrhenius_0::Array{Float64, 2}
    Troe_::Array{Float64, 2}
    index_three_body::Array{Int64, 1}
    index_falloff::Array{Int64, 1}
    index_falloff_Troe::Array{Int64, 1}
    efficiencies_coeffs::SparseMatrixCSC{Float64, Int64}
    i_reactant::Array{Array{Int64, 1},1}
    i_product::Array{Array{Int64, 1},1}
    n_reactions::Int64
    vk::SparseMatrixCSC{Float64, Int64}
    vk_sum::Array{Float64, 1}
end

struct Thermo
    nasa_low::Array{Float64, 2}
    nasa_high::Array{Float64, 2}
    Trange::Array{Float64, 2}
    isTcommon::Bool
end

struct Solution
    n_species::Int64
    n_reactions::Int64
    MW::Vector{Float64}
    species_names::Vector{String}
    elements::Vector{String}
    ele_matrix::Array{Float64, 2}
    thermo::Thermo
    reaction::Reaction
end


text:
using RoME, IncrementalInference, TransformUtils
using Distributions
using KernelDensityEstimate
using Base.Test

tf = SE3([0.0;0.0;0.0], AngleAxis(pi/4,[0;0;1.0]))# Euler(pi/4,0.0,0.0) )

N = 1
initCov = 0.0001*eye(6)
[initCov[i,i] = 0.000001 for i in 4:6];
odoCov = deepcopy(initCov)
odo = Pose3Pose3( MvNormal(veeEuler(tf), odoCov) )

X = [0.01*randn(5,N);0*pi/4+0.01*randn(1,N)]

Y = X ⊕ odo

@test norm(Y[1:3]-zeros(3)) < 1.0
@test norm(Y[4:6]-[zeros(2);pi/4]) < 0.15

tf = SE3([0.0;0.0;0.0], AngleAxis(pi/4,[0;0;1.0]))# Euler(pi/4,0.0,0.0) )

N = 1
initCov = 0.01*eye(6)
[initCov[i,i] = 0.001 for i in 4:6];
odoCov = deepcopy(initCov)
odo = Pose3Pose3(  MvNormal(veeEuler(tf), odoCov) )

X = [0.01*randn(5,N);0*pi/4+0.01*randn(1,N)]

Y = X ⊕ odo

@test norm(Y[1:3]-zeros(3)) < 1.0
@test norm(Y[4:6]-[zeros(2);pi/4]) < 0.2

N = 75
fg = initfg()

initCov = eye(6)
[initCov[i,i] = 0.01 for i in 4:6];
odoCov = deepcopy(initCov)

println("Adding PriorPose3 to graph...")
v1 = addNode!(fg, :x1,  0.1*randn(6,N),  N=N)
initPosePrior = PriorPose3( MvNormal(zeros(6), initCov) )
f1  = addFactor!(fg,[v1], initPosePrior)

println("Ensure vertex initialized properly")
# start with to tight an initialization
muX1 = Base.mean(getVal(fg,:x1),2)
stdX1 = Base.std(getVal(fg,:x1),2)
@test sum(map(Int,abs.(muX1) .< 0.1)) == 6
@test sum(map(Int, 0.05 .< stdX1 .< 0.15)) == 6

println("Testing PriorPose3 evaluation...")
priorpts = evalFactor2(fg, fg.g.vertices[2], 1)
means = Base.mean(priorpts,2)
@test sum(map(Int,abs.(means[1:3]) .> 0.5)) == 0
@test sum(map(Int,abs.(means[4:6]) .> 0.05)) == 0

println("Adding Pose3Pose3 to graph...")
odo = SE3([10;0;0], Quaternion(0))
pts0X2 = projectParticles(getVal(fg,:x1), MvNormal(veeEuler(odo), odoCov) )
odoconstr = Pose3Pose3( MvNormal(veeEuler(odo), odoCov) )
v2 = addNode!(fg,:x2,  pts0X2, N=N)
addFactor!(fg,[v1;v2],odoconstr)

println("Testing Pose3Pose3 evaluation...")
X1pts = evalFactor2(fg, fg.g.vertices[4], 1)
X2pts = evalFactor2(fg, fg.g.vertices[4], 3)
X2ptsMean = Base.mean(X2pts,2)
X1ptsMean = Base.mean(X1pts,2)
@show X1ptsMean
@test  sum(map(Int, abs.(X1ptsMean) .< 1.0 )) == 6
@test  sum(map(Int, abs.(X2ptsMean - [10.0;0;0;0;0;0]) .< 1.0 )) == 6

println("Construct Bayes tree and perform inference...")
tree = prepBatchTree!(fg);
inferOverTree!(fg, tree, N=N)

println("Ensure basic parameters on x1,x2 after inference...")
# check mean and covariances after one up and down pass over the tree
muX1 = Base.mean(getVal(fg,:x1),2)
stdX1 = Base.std(getVal(fg,:x1),2)
@test sum(map(Int,abs.(muX1[1:3]) .< 1.0)) == 3
@test sum(map(Int,abs.(muX1[4:6]) .< 0.1)) == 3
@test sum(map(Int, 0.5 .< stdX1[1:3] .< 1.5)) == 3
@test sum(map(Int, 0.025 .< stdX1[4:6] .< 0.25)) == 3
muX2 = Base.mean(getVal(fg,:x2),2)
stdX2 = Base.std(getVal(fg,:x2),2)
@test sum(map(Int, abs.(muX2[1:3]-[10.0;0;0]) .< 1.0)) == 3
@test sum(map(Int, abs.(muX2[4:6]) .< 0.1)) == 3
@show println("previous test failure 0.75 .< $(round.(stdX2[1:3],2)) .< 2.25")
@test sum(map(Int, 0.75 .< stdX2[1:3] .< 2.25)) == 3
@test sum(map(Int, 0.05 .< stdX2[4:6] .< 0.25)) == 3

# println("Plot marginals to see what is happening")
# plotKDE(marginal(getVertKDE(fg,:x1),[1]))
# plotKDE(marginal(getVertKDE(fg,:x2),[1]))

#


text:
module ManifoldTests

using ..Manifolds
using ..ManifoldsBase

"""
    find_eps(x...)

Find an appropriate tolerance for given points or tangent vectors, or their types.
"""
find_eps(x...) = find_eps(Base.promote_type(map(number_eltype, x)...))
find_eps(x::Type{TN}) where {TN<:Number} = eps(real(TN))
find_eps(x) = find_eps(number_eltype(x))

function test_action end

function test_group end

function test_manifold end

function test_reversediff end

function test_forwarddiff end

export test_action, test_group, test_manifold

end # module


text:
println("please enter the password for WeChat")
str = "wozhenbang"
while true
    jk = readline()
    if jk == str || jk == "1"
        println("the password is right, we are working on starting the WeChat")
        for i in 1:100
            print(i," ")
        end
        run(`WeChat`)
    end
end

max_stars_repo_name: Basicconstruction/Easy-tools

max_issues_repo_name: Basicconstruction/Easy-tools

max_forks_repo_name: Basicconstruction/Easy-tools