text:
@doc raw"""
    adjoint_differential_geodesic_startpoint(M,p, q, t, X)

Compute the adjoint of $D_p γ(t; p, q)[X]$.

# See also

[`differential_geodesic_startpoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
adjoint_differential_geodesic_startpoint(M::MT, p, q, t, X) where {MT <: Manifold} = adjoint_Jacobi_field(M, p, q, t, X, βdifferential_geodesic_startpoint)

@doc raw"""
    adjoint_differential_geodesic_endpoint(M, p, q, t, X)

Compute the adjoint of $D_q γ(t; p, q)[X]$.

# See also

[`differential_geodesic_endpoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_geodesic_endpoint(M::MT, p, q, t, X) where {MT <: Manifold}
    return adjoint_Jacobi_field(M, q, p, 1-t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    adjoint_differential_exp_basepoint(M, p, X, Y)

Computes the adjoint of $D_p \exp_p X[Y]$.

# See also

[`differential_exp_basepoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_exp_basepoint(M::MT, p, X, Y) where {MT <: Manifold}
    return adjoint_Jacobi_field(M, p, exp(M, p, X), 1., Y, βdifferential_exp_basepoint)
end

@doc raw"""
    adjoint_differential_exp_argument(M, p, X, Y)

Compute the adjoint of $D_X\exp_p X[Y]$.
Note that $X ∈  T_p(T_p\mathcal M) = T_p\mathcal M$ is still a tangent vector.

# See also

[`differential_exp_argument`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_exp_argument(M::mT, p, X, Y) where {mT <: Manifold}
    return adjoint_Jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_argument)
end

@doc raw"""
    adjoint_differential_log_basepoint(M, p, q, X)

computes the adjoint of $D_p log_p q[X]$.

# See also
[`differential_log_basepoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_log_basepoint(M::Manifold, p, q, X)
    return adjoint_Jacobi_field(M, p, q, 0., X, βdifferential_log_basepoint)
end

@doc raw"""
    adjoint_differential_log_argument(M, p, q, X)

Compute the adjoint of $D_q log_p q[X]$.

# See also
[`differential_log_argument`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_log_argument(M::Manifold, p, q, X)
    return adjoint_Jacobi_field(M, p, q, 1., X, βdifferential_log_argument)
end

@doc raw"""
    Y = adjoint_differential_forward_logs(M, p, X)

Compute the adjoint differential of [`forward_logs`](@ref) $F$ orrucirng,
in the power manifold array `p`, the differential of the function

$F_i(p) = \sum_{j ∈ \mathcal I_i} \log_{p_i} p_j$

where $i$ runs over all indices of the `PowerManifold` manifold `M` and $\mathcal I_i$
denotes the forward neighbors of $i$
Let $n$ be the number dimensions of the `PowerManifold` manifold (i.e. `length(size(x)`)).
Then the input tangent vector lies on the manifold $\mathcal M' = \mathcal M^n$.

# Input

* `M`     – a `PowerManifold` manifold
* `p`     – an array of points on a manifold
* `X`     – a tangent vector to from the n-fold power of `p`, where n is the `ndims` of `p`

# Ouput

`Y` – resulting tangent vector in $T_p\mathcal M$ representing the adjoint
  differentials of the logs.
"""
function adjoint_differential_forward_logs(M::PowerManifold{𝔽,TM,TSize,TPR}, p, X) where {𝔽,TM,TSize,TPR}
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = [last(R).I...] # maxInd as Array
    N = PowerManifold(M.manifold, TPR(), power_size..., d)
    Y = zero_tangent_vector(M,p)
    for i in R # iterate over all pixel
        for k in 1:d # for all direction combinations
            I = [i.I...] # array of index
            J = I .+ 1 .* (1:d .== k) #i + e_k is j
            if all( J .<= maxInd ) # is this neighbor in range?
                j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
                # Y[i] += adjoint_differential_log_basepoint(M.manifold,p[i],p[j],X[i,k])
                set_component!(
                    M,
                    Y,
                    get_component(M,Y,i) .+ adjoint_differential_log_basepoint(
                        M.manifold,
                        get_component(M,p,i),
                        get_component(M,p,j),
                        get_component(N,X,i,k),
                    ),
                    i,
                )
                # Y[j] += adjoint_differential_log_argument(M.manifold,p[i],p[j],X[i,k])
                set_component!(
                    M,
                    Y,
                    get_component(M,Y,j) + adjoint_differential_log_basepoint(
                        M.manifold,
                        get_component(M,p,i),
                        get_component(M,p,j),
                        get_component(N,X,i,k),
                    ),
                    i,
                )
            end
        end # directions
    end # i in R
    return Y
end


text:
mutable struct polynomial_t <: LCMType
    timestamp::Int64
    num_coefficients::Int32
    coefficients::Vector{Float64}
end

@lcmtypesetup(polynomial_t,
    coefficients => (num_coefficients,),
)

ext: jl


text:
module runnertest

    using Test
    using DataFrames
    using CSVFiles
    using CSV
    using Query
    using Statistics
    #
    # local imports; see ~/.julia/config/startup.jl
    #
    using GlobalDecls
    using LAModelData
    using StayingPutSim
    using CareData
    using ONSCodes

    main_results =  CareData.makecareroutcomesframe(0)

    grantdata = CSV.File( DATADIR*"edited/GRANTS_2019.csv" ) |> DataFrame

    @testset "SP Sim" begin
        alldata = CareData.load_all()
        iterations = 200
        for iteration in 1:iterations
            yp_data = CSV.File( DATADIR*"created_doe/yp_data_$iteration.csv" ) |> DataFrame
            carer_data = CSV.File( DATADIR*"created_doe/carer_data_$iteration.csv" ) |> DataFrame

            rc = size( carer_data )[1]
            ry = size( yp_data )[1]
            @test rc == ry
            params = StayingPutSim.get_default_params()
            outcomes = CareData.CarerOutcomes(0.0,0.0,0.0)
            @time for r in 1:rc
                year = yp_data[r,:year]
                ccode = yp_data[r,:ccode]
                # this test isn't strictly needed since we only create for live councils
                if (! ONSCodes.is_aggregate( ccode )) && (! ( ccode in SKIPLIST ))
                    carer = CareData.carerfromrow( carer_data[r,:] )
                    yp = CareData.ypfromrow( yp_data[r,:])
                    @assert carer.id == yp_data[r,:carer]
                    @assert year == carer_data[r,:year]
                    which = alldata.ofdata.ccode .== ccode
                    ofdata = alldata.ofdata[which,:]
                    @assert size(ofdata)[1] == 1
                    ofdatar = ofdata[1,:]
                    # println( ofdata.council )
                    if ccode in keys( FUNCS )
                        outcomes = FUNCS[ccode](
                            ccode,
                            year,
                            carer,
                            yp,
                            ofdatar,
                            params )
                    else
                        outcomes = StayingPutSim.do_basic_calc(
                            ccode,
                            year,
                            carer,
                            yp,
                            ofdatar,
                            params )
                    end
                    CareData.addcareroutcomestoframe!(
                        main_results,
                        1,
                        iteration,
                        year,
                        ccode,
                        carer.id,
                        outcomes )
                    if r % 100 == 0
                        println( carer.id )
                    end
                end # good code check
            end # main loop
        end # iteration
        CSVFiles.save( DATADIR*"created_doe/main_results2.csv", main_results )

        by_la_sys_iteration_and_year = main_results |>
            @filter( _.year > 2018 && _.year < 2025 ) |>
            @groupby( [_.ccode,  _.year, _.sysno, _.iteration ] ) |>
            @map({
                index=key(_),
                year=first(_.year),
                ccode=first(_.ccode),
                sysno=first(_.sysno),
                iteration=first(_.iteration),
                cnt=length( _ ),
                income=sum( _.income_recieved )*52.0,
                contribs = sum( _.contributions_from_yp )*52.0,
                payments = sum( _.payments_from_la )*52.0
                }

                ) |>
            @orderby( [_.year,_.ccode,_.sysno, _.iteration ] ) |>
            DataFrame
        CSVFiles.save( DATADIR*"created_doe/by_la_sys_iteration_and_year.csv", by_la_sys_iteration_and_year, delim='\t' )
        show( by_la_sys_iteration_and_year )

        m1 = by_la_sys_iteration_and_year |> @join( alldata.grantdata, _.ccode, _.ccode, { _.year, __.ccode, _.payments, __.amount } ) |> DataFrame

        by_la_sys_iteration_and_year = CSV.File( DATADIR*"created_doe/by_la_sys_iteration_and_year.csv", delim='\t' ) |> DataFrame
        by_la_sys_and_year=by_la_sys_iteration_and_year |>
            @groupby( [_.ccode,  _.year, _.sysno] ) |>
            @map(
                    {
                       ccode  = first(_.ccode),
                       year   = first(_.year),
                       sysno  = first(_.sysno),

                       avg_cnt=mean( _.cnt  ),
                       min_cnt=minimum( _.cnt ),
                       max_cnt=maximum( _.cnt ),
                       pct_10_cnt=quantile( _.cnt, [0.10] )[1],
                       pct_25_cnt=quantile( _.cnt, [0.25] )[1],
                       pct_75_cnt=quantile( _.cnt, [0.75] )[1],
                       pct_90_cnt=quantile( _.cnt, [0.90] )[1],

                       avg_income=mean( _.income  ),
                       min_income=minimum( _.income ),
                       max_income=maximum( _.income ),
                       pct_10_income=quantile( _.income, [0.10] )[1],
                       pct_25_income=quantile( _.income, [0.25] )[1],
                       pct_75_income=quantile( _.income, [0.75] )[1],
                       pct_90_income=quantile( _.income, [0.90] )[1],

                       avg_contribs = mean( _.contribs ),
                       min_contribs = minimum( _.contribs ),
                       max_contribs = maximum( _.contribs ),
                       pct_10_contribs=quantile( _.contribs, [0.10] )[1],
                       pct_25_contribs=quantile( _.contribs, [0.25] )[1],
                       pct_75_contribs=quantile( _.contribs, [0.75] )[1],
                       pct_90_contribs=quantile( _.contribs, [0.90] )[1],

                       avg_payments = mean( _.payments  ),
                       min_payments = minimum( _.payments ),
                       max_payments = maximum( _.payments ),
                       pct_10_payments=quantile( _.payments, [0.10] )[1],
                       pct_25_payments=quantile( _.payments, [0.25] )[1],
                       pct_75_payments=quantile( _.payments, [0.75] )[1],
                       pct_90_payments=quantile( _.payments, [0.90] )[1]
                    }
                ) |>
            @orderby( [ _.ccode, _.year,_.sysno] ) |>
            DataFrame

            merged = join( by_la_sys_and_year, grantdata, on=:ccode, makeunique=true )

        CSVFiles.save( DATADIR*"created_doe/by_la_sys_and_year.csv", by_la_sys_and_year, delim='\t' )
        CSVFiles.save( DATADIR*"created_doe/by_la_sys_and_year_merged_with_grant.csv", merged, delim='\t' )
    end

    @testset "grant aggregation" begin
        grantdata = CSV.File( DATADIR*"edited/GRANTS_2019.csv" ) |> DataFrame
        mg = mergegrantstoregions( grantdata )
        print( mg )
    end

end


text:
import RowEchelon

#=
Homework 6
Section 2.3
Problems 1,2,3,4,5,6,7,8,9,10
=#

# Create function to handle these questions
function checkInvertable(matrix)
    try
        inv(matrix)
        print(matrix)
        println(" IS invertable.")
    catch
        print(matrix)
        println(" is NOT invertable.")
    end
end

# Problem 1
A = [5 7 ; -3 -6]
checkInvertable(A)

# Problem 2
A = [-4 2 ; 6 -3]
checkInvertable(A)

# Problem 3
A = [3 0 0 ; -3 -4 0 ; 8 5 -3]
checkInvertable(A)

# Problem 4
A = [-5 1 4 ; 0 0 0 ; 1 4 9]
checkInvertable(A)

# Problem 5
A = [3 0 -3 ; 2 0 4 ; -4 0 7]
checkInvertable(A)

# Problem 6
A = [1 -3 -6 ; 0 4 3 ; -3 6 0]
checkInvertable(A)

# Problem 7
A = [-1 -3 0 1 ; 3 5 8 -3 ; -2 -6 3 2 ; 0 -1 2 1]
checkInvertable(A)

# Problem 8
A = [3 4 7 4 ; 0 1 4 6 ; 0 0 2 8 ; 0 0 0 1]
checkInvertable(A)

# Problem 9
A = [4 0 -3 -7 ; -6 9 9 9 ; 7 -5 10 19 ; -1 2 4 -1]
checkInvertable(A)

# Problem 10
A = [5 3 1 7 9 ;
     6 4 2 8 -8 ;
     7 5 3 10 9 ;
     9 6 4 -9 -5 ;
     8 5 2 11 4]

checkInvertable(A)


text:
# 2D FWI on Overthrust model using minConf library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: May 2017
#

using PyCall, PyPlot, HDF5, opesciSLIM.TimeModeling, opesciSLIM.SLIM_optim, SeisIO

# Load starting model
~isfile("overthrust_2D_initial_model.h5") && run(`wget ftp://ftp.slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_2D_initial_model.h5`)
n,d,o,m0 = read(h5open("overthrust_2D_initial_model.h5","r"), "n", "d", "o", "m0")

# Set up model structure
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1./model0.m)
vmin = ones(Float32,model0.n) + 0.3f0
vmax = ones(Float32,model0.n) + 5.5f0
vmin[:,1:21] = v0[:,1:21]	# keep water column fixed
vmax[:,1:21] = v0[:,1:21]

# Slowness squared [s^2/km^2]
mmin = vec((1f0./vmax).^2)
mmax = vec((1f0./vmin).^2)

# Load data
block = segy_read("/scratch/slim/pwitte/overthrust2D/overthrust_2D.segy")
dobs = joData(block)

# Set up wavelet
src_geometry = Geometry(block; key="source", segy_depth_key="SourceDepth")
wavelet = source(src_geometry.t[1],src_geometry.dt[1],0.008f0)	# 8 Hz wavelet
q = joData(src_geometry,wavelet)

############################### FWI ###########################################

# Optimization parameters
srand(1)	# set seed of random number generator
niterations = 20
batchsize = 20
proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)],2),model0.n)	# bound projection operator
fTerm = 1f-2
gradTerm = 1f-4
num_feval = 0
max_feval = 100
fval = 0
fvals = zeros(Float32,niterations)
dsub = []
qsub = []

# evaluate FWI misfit function for current model
function misfit(x)
	global num_feval; num_feval += 1
	model_new = deepcopy(model0)
	model_new.m = proj(reshape(x,model0.n))
	return f,g = fwi_objective(model_new,qsub,dsub)
end

# Backtracking line search with Armijo-Goldstein condition
function backtracking_linesearch(objective, x, grad; f_curr=[], alpha=1f0, tau=.5f0, c=.5f0, maxiter=10)

	# Parameters
	p = -grad/norm(grad)	# normalized descent direction
	m = dot(p,grad)
	t = -c*m
	isempty(f_curr) && (f_curr = objective(x)[1])	# current function value
	f_new = objective(x + alpha*p)[1]	# function value with initial step size
	iteration = 0

	# Line search
	while(f_curr - f_new < alpha*t && iteration < maxiter)
		println("	Iter LS: ", iteration, "; ",f_curr - f_new," >= ", alpha*t, "; alpha: ",alpha)
		alpha *= tau
		f_new = objective(x + alpha*p)[1]
		iteration += 1
	end
	return alpha*p
end

# Line search with armijo rule and curvature
function linesearch_wolfe(objective, x, grad_curr; f_curr=[], alpha=1f-4, tau=10f0, c1=1f-4, c2=.1f0, maxiter=10)

	# Parameters
	p = -grad_curr/norm(grad_curr)	# normalized descent direction
	isempty(f_curr) && ((f_curr, grad_curr) = objective(x))	# current function value and gradient
	f_new, grad_new = objective(x + alpha*p)	# function value and gradient with initial step size
	iteration = 0

	# line search with wolfe conditions
	while(f_new > f_curr + c1*alpha*dot(p,grad_curr) || -dot(p,grad_new) > -c2*dot(p,grad_curr) && iteration < maxiter)

		println("	Iter LS: ", iteration, "; ",f_new," <= ", f_curr + c1*alpha*dot(p,grad_curr),"; ",
		        -dot(p,grad_new), " <= ", -c2*dot(p,grad_curr), "; alpha: ",alpha)
		alpha *= 5f0
		f_new, grad_new = objective(x + alpha*p)
		iteration += 1
	end
	return alpha*p
end

# Main loop
for j=1:niterations

	# select current batch
	idx = randperm(dobs.nsrc)[1:batchsize]
	dsub = subsample(dobs,idx)
	qsub = subsample(q,idx)

	# get fwi objective function value and gradient
	fval, gradient = fwi_objective(model0,qsub,dsub)
	println("FWI iteration no: ",j,"; function value: ",fval)

	# linesearch
	step = backtracking_linesearch(misfit, vec(model0.m), gradient; f_curr=fval, alpha=10.f0, tau=.1f0, c=1e-4)
	#step = linesearch_wolfe(misfit, vec(model0.m), gradient; f_curr=fval, alpha=1f-4, tau=5f0, c1=1f-4, c2=.9f0)

	# Update model and bound projection
	model0.m = proj(model0.m + reshape(step,model0.n))

	# termination criteria
	fvals[j] = fval
	#if fval <= fTerm || norm(gradient) <= gradTerm || num_feval > max_feval
	#	break
	#end

end

# Save results
write(h5open("results_linesearch_armijo.h5", "w"), "x", sqrt.(1./reshape(x,model0.n)), "fsave", fsave, "fhistory",fvals)
# write(h5open("results_linesearch_wolfe.h5", "w"), "x", sqrt.(1./reshape(x,model0.n)), "fsave", fsave, "fhistory",fvals)

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
module VLModelParametersDB

    # include -
    include("Include.jl")

    # export data store methods -
    export get_biophysical_parameters_table_for_system
    export get_biophysical_parameters_dictionary_for_system
    export transform

    export value
    export units
    export reference

end # module


text:
############################################################
# For advanced JOLI users who want to create their own types
############################################################
# joMyOperator #########################################
############################################################
# based in joLinearFunction
#   shows minimum type and the set of methods needed
#   to fully work with other joAbstractLinearOperator types
#
#   un-implemented methods are defined in
#       joAbstractLinearOperator/base_functions.jl
############################################################

export joMyOperator, joMyOperatorException

############################################################
## type definition

"""
joMyOperator type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function              # forward
- iop::Nullable{Function}    # inverse

"""
immutable joMyOperator{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function              # forward
    iop::Nullable{Function}    # inverse
end

############################################################
## type exceptions

type joMyOperatorException <: Exception
    msg :: String
end

############################################################
## joMyOperator - outer constructors

"""
joMyOperator outer constructor

    joMyOperator(m::Integer,n::Integer,
        fop::Function,...,
        iop::Function,...,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joMyOperator")

Look up argument names in help to joMyOperator type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joMyOperator(m::Integer,n::Integer,
    fop::Function,...,
    iop::Function,...,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joMyOperatorAll") =
        joMyOperator{DDT,RDT}(name,m,n,
            fop,...,
            iop,...,
            )

############################################################
## joMyOperator - overloaded Base functions

# conj(jo)
function conj{DDT,RDT}(A::joMyOperator{DDT,RDT})
...
end

# transpose(jo)
function transpose{DDT,RDT}(A::joMyOperator{DDT,RDT})
...
end

# adjoint(jo)
function adjoint{DDT,RDT}(A::joMyOperator{DDT,RDT})
...
end

############################################################
## overloaded Base *(...jo...)

# *(jo,mvec)
function *{ADDT,ARDT,mvDT<:Number}(A::joMyOperator{ADDT,ARDT},mv::AbstractMatrix{mvDT})
...
end

# *(jo,vec)
function *{ADDT,ARDT,vDT<:Number}(A::joMyOperator{ADDT,ARDT},v::AbstractVector{vDT})
...
end

############################################################
## overloaded Base -(...jo...)

# -(jo)
function -{DDT,RDT}(A::joMyOperator{DDT,RDT})
...
end


text:
"""
    BEA

A client for the U.S. Bureau of Economic Analysis API.
"""
module BEA
end


text:
using Test, BenchmarkTools# , DataStructures # redundant - for now
#module Heap end
#=
```inputs
n: whole positive natural number > 1 -autoincremented
|---| --- |
b    b     b
```
 #queue - Datastructures.Deque{Int}()
#=
   i=4 = b ;a=[]; ; n = 1;

if (i < 2 || i > b) # stop
    #crunch in numbers
    return
else
    #modify then call
    i = i - (n * h) # 4 -(1*1) = 3 <--- new b /a
    pt = (i, b) # (3,4); push!()
    bLookup!(a, b; h = 1)
    #push!() # ∘ ⚇

    n += 1
    setindex!(a, 1, i)
end
a=rand(10,1) ; size(a)
@test
bLookup!()

(1,2)
=#
=#
#= this is coded to work with DataStructures.Deque
"""
extractall!
a Tuple operation; gets back the
"""
Base.@propagate_inbounds function extractall!(tuple::Tuple{Any,Any}) #compiles
    dt1 = nothing
    dt2 = nothing
    count = 1
    _size = size(collect((tuple)))
    typeof(minimum(tuple))
    n = _size    # length(size(collect(minimum((tuple))))) #if DataStructure.Deque was enforced
    m = _size
    #println(typeof(tuple[n, m])) # for a Deque test
    @inbounds for i = 1:m, j = 1:n
        #typeof(tuple[j, i])
        arr[j, i] = tuple[j, i][1]
        tuple[j, i] = tuple[j, i][2]
        dt1[count] = zip(tuple[j, i][1], tuple[j, i][2]) #with zip
        dt2[count] = (tuple[j, i][1], tuple[j, i][2])    #no zip
        count += 1 #auto-Increment
    end

    return dt1, dt2 #debugging, still
end
=#
tup = ((1, 2); (2, 4))
typeof(tup)
#@btime dt1, dt2 = extractall!(tup) # vector, use []
#@time dt1, dt2 = extractall!(((1, 2); (2, 4))) # vector, use []
#extractall!()
#--------------------
#=vector {any} -> tuple
Input: Tuple length N must be pre-defined
Output:
=#

Vector2TupleFloat64(N) = Tuple(Float64(x) for x in N)

@btime Float64.(tuple($N...))
#  840.075 ns (4 allocations: 912 bytes)
@btime Float64.(Tuple($N))
#  895.109 ns (5 allocations: 1.22 KiB)

vector2TupleFloat64(N) = @btime Tuple(Float64(x) for x in N); # best loop Optimized for
# @btime Tuple(1.0N);   #this line is not optimized
#  20.139 μs (159 allocations: 4.16 KiB) and on 1.0

@btime Tuple(Float64(x) for x in N); # faster, most Optimal for allocation
#  6.836 μs (57 allocations: 1.42 KiB)

#@btime Tuple(1.0N); #not optimized
#  2.660 μs (90 allocations: 2.31 KiB) # not at all
#----------
#working
#Tuple2Vector
Tuple2Vector(x) = collect(Iterators.flatten(x))
#code starts here:
Base.@propagate_inbounds function bLookup!(a = 1, b = 4; h = 1)  #optimal #rapid: median 71.5 ns

    α = 1
    β = 1 #Whatever _b could be picked-up, inside loop it 'll get overwritten
    q = [] #DataStructures.Deque{Tuple{Int64,Int64}}()
    n = 1

    while α != b && β <= b # adding inbounds slightly fastens the total process
        #   2 + (2) = 4 = _b
        β = α + (n * h)  # = 3
        push!(q, (α, β)) #1 (1,2) #infer: sub-range [1,2]  #Deque (b=4 !isa b[1]=2  (now a[2] = b[1]=2 for next op ) a starts at 2
        α = β # fetch next value a = 2
        n += 1 # 2
    end

    return q
end

#@time
#-----
copiedHeap = []::Any
function linearSort(heap)
    for i in size(heap) #compiles for a vector
        if i > 1 #this is a must
            if heap[i] < heap[i-1] #check code
                heap[i], heap[i-1] = heap[i-1], heap[i] #swap
            end
            arr[i] = copiedHeap[i]
        end #arr[i] = copiedHeap[i] copy code
    end
    return heap
end

"""checks size (for the moment)
if length isa Even , do something
Otherwise, if it's odd, throws an errror
"""
function biSectSort(heap) #fixed (this only checks size)
    _size = copy(length(heap))
    _size -= 1 #decrease always a 1 from the 'theoretical' length (to fit in the `structure `)
    try
        if _size == length(heap) % 2 == 0 # ? #= add to s : =#   this checks size # why were you throwing error , irresponsibly?
        #do something useful #like, sor, linearly
        else #on the else, ok <---- runs this
            throw(error("Debugging: this throws an error Unexpected Error found "))
        end
    catch UnexpectedError # <---- runs this
        @error "ERROR: UnexpectedError Found-error thrown successfully " exception =
            (UnexpectedError, catch_backtrace())
    end
end
biSectSort(copiedHeap) #runs # biSectSort (generic function with 1 method) #if throws your error, then it compies file

try #requires t possible error: define t first

    # for i in size(heap) #compiles for a vector

    #=find middle
    -isEven (is there a generalization)#review: generalizationnot now
    -- check middle (Q1. How? ) #Q2. isEfficient?
    lsize(heap) % 2  === 0
    #check if size(heap)/2 is even → d =#
    _size = length(copiedHeap) - 1 # -1 #indicates #ERROR (unhandled )
    if _size % 2 == 0  # TODO: odd (store in a  proper `DataStructure`)
    #do something call function

    elseif _size % 2 != 0
        #throw an error exception
        throw(error("Size Error "))
    end
catch SizeError   #Name the Error: any arbitrary error name (memorizable)
    @error "ERROR: error was unhandeled" exception = (UnexpectedError, catch_backtrace())
end #fine!
#SizeError "ERROR: error was unhandeled" exception
#=catch error as an exception  (via  catch_backtrace)#
print error
end
end
turns out, it's more than a check of Evenity...=#
#-----
@benchmark heap = bLookup!() # get subranges #Returns Deque -> Dictionary  @test(E):0.000003 seconds (4 allocations: 192 bytes) @btime =  67.551 ns (4 allocations: 192 bytes)

"""

BenchmarkTools.Trial: 10000 samples with 977 evaluations.
 Range (min … max):   66.121 ns …  11.737 μs  ┊ GC (min … max):  0.00% … 98.61%
 Time  (median):      76.254 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   131.199 ns ± 420.001 ns  ┊ GC (mean ± σ):  19.74% ±  6.19%

  ▆█▆▄▄▃▃▃▃▂▁▂▁▂▂▃▄▃▃▂▂▁▁▁ ▁ ▁       ▁ ▁ ▁                      ▂
  ███████████████████████████████▇▇██████████▆▇▇▇▆▆▇▅▆▅▅▅▅▄▅▄▄▃ █
  66.1 ns       Histogram: log(frequency) by time        318 ns <

"""
heap = Blookup.bLookup!()
copiedHeap = deepcopy(heap) # 65.169 ns (4 allocations: 192 bytes)   #65.881 ns (4 allocations: 192 bytes) #0.000003 seconds (4 allocations: 192 bytes)
#=   memory estimate:  192 bytes
  allocs estimate:  4
  --------------
  minimum time:     65.138 ns (0.00% GC)
  median time:      72.885 ns (0.00% GC)
  mean time:        127.858 ns (19.01% GC)
  maximum time:     8.332 μs (98.91% GC) =#
arr = []
#---

@benchmark heap = linearSort(heap)
arr
typeof(heap)   # was a Deque(doublyqueue) -> Dict (not a  Tuple)

length(heap)
#TODO: Deque -> Heap
#popfirst!(heap)
#typeof(pop!(heap))
#N = [] #UncommentMe
#collect(pop!(heap)) #UncommentMe
#collect(pop!(heap))[1, 1] #collect creates an Array
#orderedArray = extract_all!(heap) #
t = maximum(heap)
t = minimum(heap)
#t = popfirst!(heap)
t[1][1]
unitrange = t[1][1]:t[2][1]
typeof((t[1][1], t[2][1]))
_tuple = (t[1][1], t[2][1])

#tuple -to-> vector

vector = collect(Iterators.flatten(_tuple))

collect(Iterators.flatten(x)) #OR

[x[j] for x in x for j in eachindex(x)]

#---
#for i in heap

dims = deepcopy(size(collect(minimum(heap)))); #heap = update(heap,1,popped) #<-----Tuple Dimensions
println(dims)
typeof(dims) #dims isa tuple (if `parent` was a Deque) #Type Int64 (if `it` were )
dims[1];
dims[1][2] #BoundsError (from [2])

_length = length(dims) #Tamas_Papp Tamas_Papp Oct 2017 (2,) is simply syntax for a tuple of 1 element #needed to avoid confusion with (2) == 2.

minHeap = typeof(minimum(heap)) #(1,2)#Tuple{{Int64, Int64}} # each inside value is a Tuple  minHeap isa Tuple

#for i k in enumerate(length)
#accessing tuple :  2 loops f_vp=
#dt1 = nothing  #UncommentMe 4 #Debugging
#dt2 = nothing
#funtion iterTuple(tuple=heap, m,n)

typeHeap = typeof(copiedHeap) # Heap isa Deque
#ω ={Union{Nothing,Rational}}# nothing  # [1, 2]Vector{Int64}[2, 3, 4]Vector{Int64}
#this yeilds [1,2,3,4]

pts = []::Vector{Any} #Array{Any,1}
#= ---credits to @trakofon
https://discourse.julialang.org/t/converting-vector-any-to-tuple-of-floats/13714/7 =#
#_first = popfirst!(copiedHeap)
copiedHeap = deepcopy(heap)
typeof(copiedHeap)
size(copiedHeap)
#for i in enumerate(size(copiedHeap)) # length(copiedHeap) length for Tuple alone  #UncommentMe
#copiedHeap = deepcopy(heap)
typeof(heap)  # Deque -> Dict -> Tuple
length(heap)

dims = deepcopy(size(collect(minimum(heap)))) # == size(copiedHeap)
ω = []
copiedHeap = deepcopy(heap)
try
    while true
        tmp = popfirst!(copiedHeap)
        ω = pushfirst!(ω, tmp)#, 1) #removes the 1's monotonicity
    end
catch #ErrorException
    #continue
    print(ω) # Any[(2, 4), 1, (1, 2), 1, (2, 4), 1, (1, 2), 1] frankly the weirdest (1's monotonicity should be removable, first)
end
#=======Testing Area Experimental beyond this line
#print(ω) #this is tuple
#tuple to vector
w = collect(Iterators.flatten(ω))

typeof(collect(_first))
#------------------------------------------------------------------------

# VectorArray = []
#TupleArray = Tuple(Float64(x) for x in VectorArray);
_firstTuple = Tuple(Float64(x) for x in _first); #Float64.(tuple("$_first)) Tuple(Float64(x) for x in N);
typeof(_firstTuple)
_firstTuple
insert!(TupleArray, _firstTuple, size(TupleArray))
typeof(pts)
push!(collect((_first)), pts)
α = size(collect(_first))

#------------------------------------------------------------------------
m = deepcopy(α)
length(m)
m[2][1]
m[0]
#m2 = deepcopy()
β[i] = (m[1]:m[2]) # get rangeInt

push!(pts, m[1])
push(pts, m[2]) #ok
ω[i] = push(ω, collect(β))
push(ω, collect(β))
#TODO: need for an OrderedSet
#----
s = OrderedSet() # <: Base.AbstractSet{T} # no method matching append!(::OrderedSet{Any}, ::Int64)
#Q1. Sort an OrderedSet (Ascending)
#Q2.

#push!(s, ) #what's its continuation
typeof(s) # OrderedSet{Any}

"""OrderedSet
MIT 6-006 LECTURE 3: https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-006-introduction-to-algorithms-fall-2011/lecture-videos/MIT6_006F11_lec03.pdf

Input: array A[1…n] of numbers.
Output: permutation B[1…n] of A s.t.
    B[1] ≤ B/[2] ≤ … ≤ B[n] .
    A = [7, 2, 5, 5, 9.6] → B = [2, 5, 5, 7, 9.6]
Q. How can we do it effeciently ?
#set sorting
1.Store giant array in my set  # Inefficient (Asymptotic effeciency matters)
 if onject don't exist we add
2. how do we implement insert(set)
sortset() # binary searh  ( O(n) = nlogn : n = sortall , logn()  build set , search once )
need for a needlessly confusing algorithm
3.Search :Iterate from begining of orderedSet , & return         [#worstcase n times (if found at last element)  ]
    OrderedSet{Any}
    all (A, n)()
Strategy:
3.1.1. BuildSet : Reserve n slts in memory  - order n times
3.1.2. Copyall into set
3.1.3. Update (insert/delete : `the Amortized Argument`:if `set` is not allowed, Grow Dynamically - Reserve memory )
3.1.4. Find minimum: only algorithm , list `things` isn't sorted:iterate over every single One
e.g. If sth I'm seeking has a smaller Id (Index) : Replace Histogram (you guys are More than `Qualified` to Implement !)
(while) Set Interface is Dynamic - keep adding stuff to it )
`the Amortized Argument` : on average: (How much) it will take order n times i.e. Expectation (assuming naive)
would be the mean E(x)

Note: even if it's not dynamic (i.e. to replace an existing key (Index))
    -2 numbers exist
    _"first,there is the existance of a solution, but paying attention to the Details is what really matters"_
takeaway: there is a difference between:
1. the existence of a Solution
2. actually paying attention to the Detail (inside of it)

situation :
real world semi-blackbox using only the following 6 functions:

sizehint!(s::OrderedSet, sz::Integer) = (sizehint!(s.dict, sz); s) #capacity reservation for improving Optimization
in(x, s::OrderedSet) = haskey(s.dict, x)

push!(s::OrderedSet, x) = (s.dict[x] = nothing; s)
pop!(s::OrderedSet, x) = (pop!(s.dict, x); x)
pop!(s::OrderedSet, x, deflt) = pop!(s.dict, x, deflt) == deflt ? deflt : x
delete!(s::OrderedSet, x) = (delete!(s.dict, x); s)

note: no push first nor pop first !
"""
#--- tuple -to-> vector

# Tuple -to -> Vect
[[Vector{Float64}(undef, 4) for _ = 1:5] for _ = 1:6]
#5-element Array{Array{Array{Float64,1},1},1}:

struct Run
    vector_result::Vector{Float64}
    matrix_result::Matrix{Float64}
end

struct Simulation
    #name::String
    parameters
    runs::Vector{Run}
end

s = Simulation(some_parameters, Vector{Run}())
#---- end
_size = sizehint!(s, 1);
println(_size);
last = pop!(s, 1)

#=
Base.@propagate_inbounds function insert!(s::OrderedSet{Any}, ::Any{T})

end
=#

Base.@propagate_inbounds function insert!(s::OrderedSet{Any}, ::Any{T}) where {T::Any}
    #1. check size
    _#size = sizehint!(s)
    #look for safety with try (no matter the size of s )
    canLah = true
    while canLah
        try
            last = pop!(s, 1) #reminiscent of the tower of hanoi
        #checkva
        catch #can't pop anymore
            canLah = false
            break
        end
    end

end

#---
_size = size(collect(minimum(copiedHeap)))  #<-- the tuple
#m = _size
#typeof(m) # m is a tuple
#m[1]m[1][1]

Size = length(copiedHeap) # the true number looking for

isequal = nothing
_size == Size ? isequal = true : isequal = false
println(isequal)
#TODO: construct a full loop  & iterate ( ort) =#

Base.@propagate_inbounds function getPoints(copiedHeap)
    _size = size(collect(minimum(copiedHeap)))  #<-- the tuple
    #m = _size
    #typeof(m) # m is a tuple
    #m[1]m[1][1]

    Size = length(copiedHeap)

    pts = []
    α = []
    β = []  #  Array{Int64}(Int64, size)
    ω = []
    i = 1 # must be defined
    γ = nothing #[]
    α = size(collect(minimum(copiedHeap)))  #<-- the tuple
    typeof(α)
    # the Strategy:
    for i in enumerate(Size) #  # \A tuple in Deque[Tuple{Int64, Int64}]'s   thing isa `Humongous` [Algorithm allows for an Arbitrary  number of k ]
        # print(typeof(k)) #Tuple{Int64, Int64}Tuple{Int64, Int64}
        #print(k)#Tuple{Int64, Int64}(1, 2)Tuple{Int64, Int64}(2, 4) # the correct Ideal subranges we want #Extravagant!
        #println(length(k))
        # Ideally,    β[i] = α[1]:α[2] #β is a rangeInt    #UncommentMe
        #1. find m
        #m = size(collect(minimum(copiedHeap)))
        #2. pop(first),  if possible
        α = size(collect(popfirst!(copiedHeap)))  #<-- the tuple
        typeof(α)
        #m = _size
        #typeof(m) # m i a
        #m[1]m[1][1] ok
        m = copy(α) #the intent #deepcopy(α)

        @inbounds β[i] = (m[1]:m[2]) # a tuple # gets rangeInt  #without Deque, errors out

        push!(pts, m[1]) #ok you mean β[i] ?
        push(pts, m[2]) #ok   or pushfirst!
        #  ω[i] = β[i] #collect(β))
        push(ω, collect(β[i]))
        #append!(γ, collect(β)) #no method matching append!(::OrderedSet{Any}, ::Vector{Int64})#sloved
        #print(typeof(α[1]))
        # push!(γ, m[1]) # push: method no method matching append!(::OrderedSet{Any}, ::Int64)
        # push!(γ, m[2]) #or pushfirst!
        #print(ω)
        #print(γ)
        #print(typeof(ω))
        #print(typeof(γ))
        #println()
        #println(α[1], α[2]) #UncommentMe
        #=
        for β in length() #access a desired subrange α (me: (1,2) or (2,4))
            println(β[1]);#println(typeof(β))
        end
        =#
        return ω, γ, β #, pts #, α[1]:α[2] #  debugging values
    end
end # no error
ω, γ, β = getPoints(copiedHeap) #TODO: re-view

function append!(s ::OrderedSet{Any}, a::Int64) #use ] OrderedCollections (for OrderedCollections.jl) #q. why Unordered Structure?  #OrderedSet not defined #review#1: use of outside code
return
end
#=
function append!(s ::OrderedSet{Any}, a::Int64) #q. why Unordered Structure?

@inbounds for i in enumerate length(s)
    @inbounds if i!= length(s)
        @inbounds if a > s[i] && a< s[i+1]
    #found the right place
                    return i # value at i
                    #insert a
                    #break  #if it returns, no need to break

elseif  a == s[i] # same as the one in set
    break; #skip
end
#---- Algorithms

end
"""merge sort (with Difive & Conquer Strategy) #the reason behind divideConquer
credits to @edubkendo
https://gist.github.com/edubkendo/528d40034fd7037c55ce
"""

function mergeSort(lhs::Array, rhs::Array)
    size = length(lhs) + length(rhs)
    retvals = [] #Array(Type(lhs), size)

@inbounds    for k = 1:size
@inbounds        if first(lhs) <= first(rhs)
        @inline    retvals[k] = shift!(lhs)
        else
        @inline    retvals[k] = shift!(rhs)
        end
        if isempty(lhs) #Divide & Conquer Strategy
            return transfer_tail(retvals, rhs, k)
        elseif isempty(rhs)
            return transfer_tail(retvals, lhs, k)
        end
    end
end

#ok
function transfer_tail(vals::Array, tail::Array, count::Int64) #TODO: #check_Utility
@inbounds  for k = (count + 1):(length(tail) + count)
@inbounds    vals[k] = shift!(tail)
  end
  vals
end

#sortmerge
function sortmerge(to_sort::Array)
  len = length(to_sort)
  len <= 1 && return to_sort
  head = to_sort[1:(ceil(len/2))]
  tail = to_sort[(length(head)+1):end]
  merge_sorted(sortmerge(head), sortmerge(tail))
end

@inline x = sortmerge([1,3,2,4]) #tryout
@inine y = sortmerge([1,2,3,4,5, 99, 54, 33, 21, 67, 88, 9, 39, 76, 88, 89, 100001, 45, 3000, 6,867,8,10]) #tryout

println("x: $x \n y: $y")

function heap2Range()

@inbounds for α in heap # α tuple in Deque[Tuple{Int64, Int64}]  a tuple `Gigantic` one tuple percieved number is 1  [Algorithm allows for an Arbitrary  number of k ]
    # print(typeof(k)) #Tuple{Int64, Int64}Tuple{Int64, Int64}
    #print(k)#Tuple{Int64, Int64}(1, 2)Tuple{Int64, Int64}(2, 4) # the correct Ideal subranges we want #Extravagant!
    #println(length(k))
    β = α[1]:α[2]
    ω = collect(β)
    print(ω)

    print(typeof(ω))
    #println()
    #println(α[1], α[2]) #UncommentMe
    #=
    for β in length(α) #access a desired subrange α (me: (1,2) or (2,4))
        println(β[1]);#println(typeof(β))
    end
    =#
end
#return w
#TODO: Tuple to array operation

_size = size(collect(minimum(copiedHeap)))  #<-- the tuple
m = _size

typeof(m) # m is a tuple

#m[1] m[1][1];
dt1 = []; dt2 = [];
_size =size(collect(minimum(tuple)))
typeof(minimum(tuple))
n = length(m)
count = 1
@inbounds for i = 1:m, j = 1:n
    @inline  dt1[count] =  collect(zip(copiedHeap[j, i][1], copiedHeap[j, i][2]))
    @inline  dt2[count] =  (copiedHeap[j, i][1], copiedHeap[j, i][2])
  #  f_value1[j, i] = copiedHeap[j, i][1]
   # f_value2[j, i] = copiedHeap[j, i][2]
 #  arr1 = dt1
#collect(dt2)
@inline count += 1 #auto-Increment
end

#typeof(dt1)
#typeof(dt2)
res2 = nothing
for i in enumerate(len)
    res2 =(collect(zip((popfirst!(heap))))[1, 1])

end
println(res2); println(typeof(res2));
size(A[1:1, 1])

#--------------------
#arrays indexing rewind

#N[[1],:] =typeof(collect(pop!(heap)))# boundsError #UncommentMe
N[[:],] = typeof(collect(pop!(heap))) #possible :Deque must me non-empty # invalid index Colon #LoadError: ArgumentError: invalid index: [Colon()] of type Vector{Colon}
N[:] =typeof(collect(pop!(heap))) # returns a 0-element vector
resultsize(idxs...) = tuple(Iterators.flatten(size.(idxs))...)

#0element vector{Any} at index [[1]]
# typeof(collect(pop!(heap)))
#typeof( collect(pop!(map(x->x, heap))))

containr = nothing
arr = nothing
@inbounds for i in enumerate(length(heap))
    @inline  containr += [arr collect(pop!(map(x[i]->x, heap)) ) ]
end
#=" how you define behavior is how
    julia will react "=#

ranges = []
a = 1; b= 4
#ret = nothing
A=[];B=[] #;I=[];
#q =(,)::Tuple{UnitRange{Int64}, Tuple{Int64, Int64}}, ::Type{Any} # Any #DataStructures.Deque{UnitRange{Int64}}()

#insert!()

A = collect( popfirst!(heap))

@inbounds for (n,i) in enumerate(tmp)
    @inline   b[n] = i+1
end

a,b = popfirst!(heap)

function extractfromHeap(heap)
    A = [] #1
    B = []
    #I = []
    p=1 ;a=1;b=1;
@inbounds  for p in enumerate(length(heap))
        a, b = popfirst!(heap)
        #ret = [ret, a:b]
        #push!(A, a)
        #push!(B, b)
    A = @inline InsertionSort(collect(zip(a,b,p)))
    B = @inline collect(zip(a:b,p))
    # insert!(A, a, p)        #=@inbounds=#
       # insert!(B, b, p)        #=@inbounds=#
       # insert!(I, p, p)        #=@inbounds=#
    end
    # A = [a] #1
    # B = [b]
    # I = [p]
    #   (A,B)[i]  = (a,b)
    #i += 1
    #(A[i],B[i]) = (a,b)
    #   (A[i],B[i],I[i]) = a,b,i
    # B[i] = b ;
    #  I[i] = i
    #   a[i],B[i],I[i] = a,b,i            #push!((a:b,i),q)
    return A, B #for Debugging
A,B = extractfromHeap(copiedHeap)

end

#ret[1]
#ret[2]
    #ranges  = [ ranges ,a:b]
    #range1 = ret1:ret2
    #ret3, ret3 = popfirst!(heap)

#range2 =
# ret2
#ret1:ret2

#--------------------
module Deque
function getIndex!(Q::Tuple{UnitRange{Int64}}, idx::Int64)
    ret = getQueue(Q, idx)
end

function getindex!(Q::Tuple{Tuple{Int64, Int64}},idx::Int64)
    ret = nothing
    for i in enumerate(size(Q))
        if i isa idx
        ret = pop!(Q)
        end
    end
    return ret
end
end
#=--- testing Area
getindex!(heap,1)
i = 4  #... i = 1 ; b[1] =4
a[1]= b[1] - (n[1]*h) =

b[2]

b[1]  # setindex!(b,1,)
b[1] =  1 + (1 * 1 )  = 2  #infer: sub-range [1,2]  (b=4 !isa b[1]=2  (now a[2] = b[1]=2 for next op ) a starts at 2
b[2] = a[2] + (n[2] * h) = 2 + (2 * 1)

b[2] = a[2] - (n[2] * h) = 2 -  (2 * 1)
#=
Base.setindex!
=# #Function
setindex!(collection, value, key...)
Store the given value at the given key or index within a collection. The syntax a[i,j,...] = x is converted by the compiler to (setindex!(a, x, i, j, ...); x).
=#
=#
"""trying to recursively find a subrage *writes down a while!* """
function  findSubrange2(a=1,b=4,n= 1) #requires revision #review#1: try to remove ; from the arguuments

    i=copy(n) #copies n

    #a=[]; _b=[]
    #_n = 1 #review#1:instead of writing here, you are passing it in the argument - as an optional parameter , thus if user forgets it , you expect it to run
    _a =1
    _b =1

    if _b ==b # isa b #reached the end #review#1:comment:1stif
    _a = a  #
    return
    #end #end-if
   #end #end-function

    if !(_b isa b) #review#1:comment:2ndif
        _a + (n[1] * h)
        i += 1

        while !(_b isa b)# use of while
            _b = _a + (_n * h)
            #update _a, _b
            i = 1:_b = _a + (i * h) # = 1 + (1 * 1) = 2  #infers: sub-range [1,2]  (b=4 !isa b[1]=2  (now a[2] = b[1]=2 for next op ) a starts at 2
            i += 1
            findSubrange(_a, _b, n = i)
        end #end-while

     #end#end-if
    end ##-end-if error:syntax "1" - not a `valid argument name` #TODO: check name
end #incomplete function #error

function whilefunction(a=_a,b=_b) #error at 706
    while _b  != b #isa b)#review#1:  use of while
         _b = _a + (_n * h)
        #update _a, _b
        i = 1:_b = _a + (i * h) # = 1 + (1 * 1) = 2  #infers: sub-range [1,2]  (b=4 !isa b[1]=2  (now a[2] = b[1]=2 for next op ) a starts at 2
        i += 1
        findSubrange(_a, _b, n = i)
    end #end-while
end
#--------------------=# #was erroring, because of `unterminatedError`!  -file formation complete

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
if abspath(PROGRAM_FILE) == @__FILE__
    l = [1, 2, 3]
    b = ["a", "b", "c"]
    x = 0
    @assert(l[x+2:end] == [2, 3])
    x = 1
    @assert(b[x+2:end] == ["c"])
    output = [1, 2, 3, 4, 5, 6]
    start = 1
    stop = 3
    @assert(output[begin:stop-start] == [1, 2])
    @assert(output[length(output):end] == [6])
    @assert(output[end] == 6)
    println("OK")
end


text:
@testset "set.jl" begin
    @testset "basic" begin
        s = AVLSet()
        items = ["anything", "anything2"]
        push!(s, items[1])
        r = collect(s)
        @test length(r) == 1
        @test items[1] in r
        push!(s, items[2])
        r = collect(s)
        @test length(r) == 2
        @test all(items .∈ Ref(s))
        @test all(items .∈ Ref(r))
        delete!(s, items[1])
        delete!(s, items[2])
        @test isempty(s)
    end

    @testset "constructor with AbstractVector input" begin
        items = rand(1_000)
        s = AVLSet(items)
        @test eltype(items) == eltype(s)
        @test all(items .∈ Ref(s))
        @test all(items .∈ Ref(collect(s)))
    end

    @testset "union" begin
        a = rand(1:1000, 800)
        b = rand(1:1000, 800)

        avl_a = AVLSet(a)
        avl_b = AVLSet(b)

        sa = Set(a)
        sb = Set(b)

        @test union(avl_a, avl_b) == union(sa, sb)
        @test union(avl_a, avl_a) == union(sa, sa)
        @test union(avl_a, avl_b, avl_b) == union(sa, sb)
        @test union(avl_a, sb, avl_b, avl_a, sb) == union(sa, sb)
        @test union!(avl_a, avl_b) == union(sa, sb)
        @test union!(avl_b, avl_a) == union(sa, sb)
    end

    @testset "setdiff" begin
        a = rand(1:1000, 800)
        b = rand(1:1000, 800)

        avl_a = AVLSet(a)
        avl_b = AVLSet(b)

        sa = Set(a)
        sb = Set(b)

        @test setdiff(avl_a, avl_b, avl_b) == setdiff(sa, sb)
        @test setdiff(avl_a, avl_a) == setdiff(sa, sa)
        @test setdiff(avl_a, sb) == setdiff(sa, sb)
    end

    @testset "intersect" begin
        a = rand(1:1000, 800)
        b = rand(1:1000, 800)

        avl_a = AVLSet(a)
        avl_b = AVLSet(b)

        sa = Set(a)
        sb = Set(b)

        @test intersect(avl_a, avl_b) == intersect(sa, sb)
        @test intersect(avl_a, avl_a) == intersect(sa, sa)

    end
end


text:
module JuliaKara_Base_GUI
using Gtk,Gtk.ShortNames, Graphics
using Cairo

export
    Grid,
    grid_draw,
    world_init,
    symbol_image,
    cover_field

res_path = joinpath(@__DIR__,"..","res","icons")

# Should be initialized at runtime
icons = Dict(
)

function __init__()
    icons[:kara] = read_from_png(joinpath(res_path,"bugnorth.png"))
    icons[:tree] = read_from_png(joinpath(res_path,"object_tree.png"))
    icons[:mushroom] = read_from_png(joinpath(res_path,"object_mushroom.png"))
    icons[:leaf] = read_from_png(joinpath(res_path,"object_leaf.png"))
end

struct Grid{Tx<:Real,Ty<:Real,Tw<:Real,Th<:Real}
    x::Tx
    y::Ty
    width::Tw
    height::Th
    xe::Int
    ye::Int
end

function grid_draw(gr::Grid,ctx::Gtk.CairoContext)
    new_path(ctx)
    move_to(ctx,gr.x+0.5,gr.y+0.5)
    set_source_rgb(ctx,0,0,0)
    for x in range(gr.x+0.5,stop=gr.x+0.5+gr.width,length=gr.xe+1)
        move_to(ctx,x,gr.y+0.5+gr.height)
        rel_line_to(ctx,0,-gr.height)
    end
    for y in range(gr.y+0.5,stop=gr.y+0.5+gr.height,length=gr.ye+1)
        move_to(ctx,gr.x+0.5,y)
        rel_line_to(ctx,gr.width,0)
    end
    stroke(ctx)
end

function grid_coordinate_real(gr::Grid,x::Int,y::Int)
    gr.x+0.5 + gr.width/gr.xe*(x-1),
    gr.y+0.5 + gr.height/gr.ye*(gr.ye-y)
end

function grid_coordinate_virt(gr::Grid,xr::Float64,yr::Float64)
    Int(floor((xr-gr.x+0.5)*gr.xe/gr.width+1)),
    Int(ceil(gr.ye - (yr-gr.y+0.5)*gr.ye/gr.height))
end

function world_init(title::AbstractString)
    b = GtkBuilder(filename=joinpath(@__DIR__,"layout.glade"))
    c = @GtkCanvas()
    fc = b["frame_canvas"]
    push!(fc,c)
    win = b["win_main"]
    showall(win)
    return b,win,c
end

function cover_field(gr::Grid,ctx::Gtk.CairoContext,x::Int,y::Int)
    wi = gr.width/gr.xe
    hi = gr.height/gr.ye
    xr,yr = grid_coordinate_real(gr,x,y)
    set_source_rgb(ctx,0.85,0.85,0.85)
    rectangle(ctx,xr,yr,wi,hi)
    fill(ctx)
    rectangle(ctx,xr,yr,wi,hi)
    stroke(ctx)
    set_source_rgb(ctx,0.0,0.0,0.0)
    rectangle(ctx,xr,yr,wi,hi)
    stroke(ctx)
end

function symbol_image(gr::Grid,ctx::Gtk.CairoContext,x::Int,y::Int,angle::T,image::Symbol) where T <: Real
    cairo_matrix = get_matrix(ctx)
    wi = gr.width/gr.xe
    hi = gr.height/gr.ye
    xr,yr = grid_coordinate_real(gr,x,y)
    img = icons[image]
    # Move to the center node
    translate(ctx,xr,yr)
    translate(ctx,wi/2,hi/2)
    # Rotate in the center node
    rotate(ctx,-angle)
    # Move back to the upper left node
    translate(ctx,-wi/2,-hi/2)
    # Scale the image to fit the grid
    scale(ctx,wi/img.width,hi/img.height)
    set_source_surface(ctx,img,0,0)
    paint(ctx)
    set_matrix(ctx,cairo_matrix)
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""VolumeAttachmentSpec is the specification of a VolumeAttachment request.

    IoK8sApiStorageV1beta1VolumeAttachmentSpec(;
        attacher=nothing,
        nodeName=nothing,
        source=nothing,
    )

    - attacher::String : Attacher indicates the name of the volume driver that MUST handle this request. This is the name returned by GetPluginName().
    - nodeName::String : The node that the volume should be attached to.
    - source::IoK8sApiStorageV1beta1VolumeAttachmentSource : Source represents the volume that should be attached.
"""
mutable struct IoK8sApiStorageV1beta1VolumeAttachmentSpec <: SwaggerModel
    attacher::Any # spec type: Union{ Nothing, String } # spec name: attacher
    nodeName::Any # spec type: Union{ Nothing, String } # spec name: nodeName
    source::Any # spec type: Union{ Nothing, IoK8sApiStorageV1beta1VolumeAttachmentSource } # spec name: source

    function IoK8sApiStorageV1beta1VolumeAttachmentSpec(;attacher=nothing, nodeName=nothing, source=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1beta1VolumeAttachmentSpec, Symbol("attacher"), attacher)
        setfield!(o, Symbol("attacher"), attacher)
        validate_property(IoK8sApiStorageV1beta1VolumeAttachmentSpec, Symbol("nodeName"), nodeName)
        setfield!(o, Symbol("nodeName"), nodeName)
        validate_property(IoK8sApiStorageV1beta1VolumeAttachmentSpec, Symbol("source"), source)
        setfield!(o, Symbol("source"), source)
        o
    end
end # type IoK8sApiStorageV1beta1VolumeAttachmentSpec

const _property_map_IoK8sApiStorageV1beta1VolumeAttachmentSpec = Dict{Symbol,Symbol}(Symbol("attacher")=>Symbol("attacher"), Symbol("nodeName")=>Symbol("nodeName"), Symbol("source")=>Symbol("source"))
const _property_types_IoK8sApiStorageV1beta1VolumeAttachmentSpec = Dict{Symbol,String}(Symbol("attacher")=>"String", Symbol("nodeName")=>"String", Symbol("source")=>"IoK8sApiStorageV1beta1VolumeAttachmentSource")
Base.propertynames(::Type{ IoK8sApiStorageV1beta1VolumeAttachmentSpec }) = collect(keys(_property_map_IoK8sApiStorageV1beta1VolumeAttachmentSpec))
Swagger.property_type(::Type{ IoK8sApiStorageV1beta1VolumeAttachmentSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1beta1VolumeAttachmentSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1beta1VolumeAttachmentSpec }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1beta1VolumeAttachmentSpec[property_name]

function check_required(o::IoK8sApiStorageV1beta1VolumeAttachmentSpec)
    (getproperty(o, Symbol("attacher")) === nothing) && (return false)
    (getproperty(o, Symbol("nodeName")) === nothing) && (return false)
    (getproperty(o, Symbol("source")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiStorageV1beta1VolumeAttachmentSpec }, name::Symbol, val)
end


text:
include("activate.jl")

using LinearAlgebra

n = Int(1e3);
k = 1e5

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

f(x) = 2 * norm(x - xp)^3 - norm(x)^2

# this is just for the example -> better explicitly define your gradient
grad!(storage, x) = ReverseDiff.gradient!(storage, f, x)

# pick feasible region
lmo = FrankWolfe.ProbabilitySimplexOracle(1.0); #radius needs to be float

# compute some initial vertex
x0 = collect(FrankWolfe.compute_extreme_point(lmo, zeros(n)))

# benchmarking Oracles
FrankWolfe.benchmark_oracles(f, grad!, () -> randn(n), lmo; k=100)

@time x, v, primal, dual_gap, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Nonconvex(),
    print_iter=k / 10,
    verbose=true,
);


text:
# ------------------------------------------------------------------------------------------
# # Multilayer Perceptron (MLP)
# (c) Deniz Yuret, 2019
# * Objectives: Stack multiple layers, demonstrate the power of nonlinear activation
# functions, use regularization and dropout to address overfitting.
# * Prerequisites: [Linear models](30.lin.ipynb), [MNIST](20.mnist.ipynb)
# * New functions:
# [relu](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.relu),
# [dropout](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.dropout),
# [param, param0](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.param),
# [xavier](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.xavier),
# [training](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.training),
# [gc](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.gc)
#
# ![image](https://github.com/denizyuret/Knet.jl/blob/master/docs/src/images/neural_net2.jpe
# g?raw=true)
# ([image source](http://cs231n.github.io/neural-networks-1))
#
# A multilayer perceptron, i.e. a fully connected feed-forward neural network, is basically
# a bunch of linear layers stuck together with non-linearities in between. In this notebook
# we will first chain together a couple of linear layers and see that it has no effect on
# performance (multiple linear layers still compute a linear function). However with a one
# line change (adding an elementwise nonlinear function between layers) the accuracy will go
# from 92% to 98%! It turns out the class of functions that can be represented by linear
# models is severely restricted and does not contain good MNIST classifiers. MLPs have no
# such restrictions, being **universal function approximators** they can represent any
# classifier given large enough layers and learn the correct classifier with sufficient
# data.
#
# Simple models with low capacity, e.g. the linear model, have an **underfitting** problem:
# the error does not drop to 0 even on the training set. Models with high capacity, e.g.
# MLP, have an **overfitting** problem: even when the training set error goes to 0, the test
# set error remains high. We will implement two methods to fight overfitting:
# **regularization** and **dropout**.
# ------------------------------------------------------------------------------------------

# Set display width, load packages, import symbols
ENV["COLUMNS"]=72
using Pkg; for p in ("Knet","Plots"); haskey(Pkg.installed(),p) || Pkg.add(p); end
using Knet: Knet, dir, zeroone, progress, sgd, load, save, gc, Param, KnetArray, gpu, Data, nll, relu, training, dropout # param, param0, xavier
using Statistics: mean
using Base.Iterators: flatten

# Load data (see mnist.ipynb)
include(Knet.dir("data","mnist.jl"))  # Load data
dtrn,dtst = mnistdata(xsize=(784,:)); # dtrn and dtst = [ (x1,y1), (x2,y2), ... ] where xi,yi are minibatches of 100

# For running experiments
function trainresults(file,model; o...)
    if (print("Train from scratch? "); readline()[1]=='y')
        takeevery(n,itr) = (x for (i,x) in enumerate(itr) if i % n == 1)
        r = ((model(dtrn), model(dtst), zeroone(model,dtrn), zeroone(model,dtst))
             for x in takeevery(length(dtrn), progress(sgd(model,repeat(dtrn,100)))))
        r = reshape(collect(Float32,flatten(r)),(4,:))
        Knet.save(file,"results",r)
        Knet.gc() # To save gpu memory
    else
        isfile(file) || download("http://people.csail.mit.edu/deniz/models/tutorial/$file",file)
        r = Knet.load(file,"results")
    end
    println(minimum(r,dims=2))
    return r
end

# Some utilities to make model definitions easier:
param(d...; init=xavier, atype=atype())=Param(atype(init(d...)))
param0(d...; atype=atype())=param(d...; init=zeros, atype=atype)
xavier(o,i) = (s = sqrt(2/(i+o)); 2s .* rand(o,i) .- s)
atype()=(gpu() >= 0 ? KnetArray{Float32} : Array{Float32})

# ------------------------------------------------------------------------------------------
# ## A generic multilayer model
# ------------------------------------------------------------------------------------------

# Let's define a chain of layers
struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)
(c::Chain)(d::Data) = mean(c(x,y) for (x,y) in d)

# ------------------------------------------------------------------------------------------
# ## Multiple linear layers do not improve over a single linear layer
# ------------------------------------------------------------------------------------------

# Define a linear layer (See lin.ipynb):
struct Layer0; w; b; end
Layer0(i::Int,o::Int) = Layer0(param(o,i),param0(o))
(l::Layer0)(x) = (l.w * x .+ l.b)

# Here is an example two layer model
model=Chain(Layer0(784,64), Layer0(64,10))
println.(summary.((l.w,l.b)) for l in model.layers);

# Train the two layer model
# 2.04e-01  100.00%┣████████████████┫ 60000/60000 [00:51/00:51, 1179.52i/s]
# [0.240742; 0.282535; 0.0691333; 0.0794]
mlp1 = trainresults("mlp113a.jld2", model);

# Here is a single layer (linear) model
model=Chain(Layer0(784,10))
println.(summary.((l.w,l.b)) for l in model.layers);

# Train the single layer (linear) model
# 2.30e-01  100.00%┣████████████████┫ 60000/60000 [00:40/00:40, 1502.26i/s]
# [0.242363; 0.267073; 0.0671333; 0.0746]
lin1 = trainresults("mlp113b.jld2", model);

using Plots; default(fmt=:png,ls=:auto)

# multilinear converges to a similar solution, not identical because problem is non-convex
plot([lin1[1,:], lin1[2,:], mlp1[1,:], mlp1[2,:]], ylim=(0.0,0.4),
    labels=[:trnLin :tstLin :trnMulti :tstMulti],xlabel="Epochs",ylabel="Loss")

# error results also close to the linear model
plot([lin1[3,:], lin1[4,:], mlp1[3,:], mlp1[4,:]], ylim=(0.0,0.1),
    labels=[:trnLin,:tstLin,:trnMulti,:tstMulti], xlabel="Epochs", ylabel="Error")

# ------------------------------------------------------------------------------------------
# ## Multiple linear layers are useless because they are equivalent to a single linear layer
# If we write down what is being computed and do some algebra, we can show that what is
# being computed is still an affine function of the input, i.e. stacking multiple linear
# layers does not increase the representational capacity of the model:
#
# \begin{align*}
# \hat{p} &= \mbox{soft}(W_2 (W_1 x + b_1) + b_2) \\
# &= \mbox{soft}((W_2 W_1)\, x + W_2 b_1 + b_2) \\
# &= \mbox{soft}(W x + b)
# \end{align*}
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Multi Layer Perceptron (MLP) adds non-linearities between layers
# ------------------------------------------------------------------------------------------

# Using nonlinearities between layers results in a model with higher capacity and helps underfitting
# relu(x)=max(0,x) is a popular function used for this purpose, it replaces all negative values with zeros.
struct Layer1; w; b; f; end
Layer1(i::Int,o::Int,f=relu) = Layer1(param(o,i),param0(o),f)
(l::Layer1)(x) = l.f.(l.w * x .+ l.b)

# We add a nonlinear activation function to all but the last layer
# 2.83e-02  100.00%┣████████████████┫ 60000/60000 [00:52/00:52, 1150.34i/s]
# [0.00636037; 0.084129; 0.000516667; 0.0243]
model = Chain(Layer1(784,64), Layer1(64,10,identity))
mlp2 = trainresults("mlp113c.jld2", model);

# ------------------------------------------------------------------------------------------
# ## MLP solves underfitting but still has an overfitting problem
# ------------------------------------------------------------------------------------------

# MLP solves the underfitting problem!  A more serious overfitting problem remains.
plot([lin1[1,:], lin1[2,:], mlp2[1,:], mlp2[2,:]], ylim=(0.0,0.4),
     labels=[:trnLin :tstLin :trnMLP :tstMLP],xlabel="Epochs",ylabel="Loss")

# Test error improves from 7.5% to 2.5%!
plot([lin1[3,:], lin1[4,:], mlp2[3,:], mlp2[4,:]], ylim=(0.0,0.1),
    labels=[:trnLin,:tstLin,:trnMLP,:tstMLP], xlabel="Epochs", ylabel="Error")

# ------------------------------------------------------------------------------------------
# ## MLP with L1/L2 regularization
# ------------------------------------------------------------------------------------------

# We add two new fields for L1 and L2 regularization
struct Chain2
    layers; λ1; λ2
    Chain2(layers...; λ1=0, λ2=0) = new(layers, λ1, λ2)
end

# The prediction and average loss do not change
(c::Chain2)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain2)(d::Data) = mean(c(x,y) for (x,y) in d)

# The loss function penalizes the L1 and/or L2 norms of parameters during training
function (c::Chain2)(x,y)
    loss = nll(c(x),y)
    if training() # Only apply regularization during training, only to weights, not biases.
        c.λ1 != 0 && (loss += c.λ1 * sum(sum(abs, l.w) for l in c.layers))
        c.λ2 != 0 && (loss += c.λ2 * sum(sum(abs2,l.w) for l in c.layers))
    end
    return loss
end

# 1.94e-01  100.00%┣████████████████┫ 60000/60000 [00:57/00:57, 1058.05i/s]
# [0.026522; 0.0802438; 0.00688333; 0.0231]
model = Chain2(Layer1(784,64), Layer1(64,10,identity); λ1=4f-5)
mlp3 = trainresults("mlp113d.jld2", model);

# overfitting less, test loss improves from 0.0841 to 0.0802
plot([mlp2[1,:], mlp2[2,:], mlp3[1,:], mlp3[2,:]], ylim=(0.0,0.15),
     labels=[:trnMLP :tstMLP :trnL1 :tstL1],xlabel="Epochs",ylabel="Loss")

# test error also improves: 0.0243 -> 0.0231
plot([mlp2[3,:], mlp2[4,:], mlp3[3,:], mlp3[4,:]], ylim=(0.0,0.04),
     labels=[:trnMLP :tstMLP :trnL1 :tstL1],xlabel="Epochs",ylabel="Error")

# ------------------------------------------------------------------------------------------
# ## MLP with dropout
# ------------------------------------------------------------------------------------------

@doc dropout

# Dropout is another way to address overfitting
struct Layer2; w; b; f; pdrop; end
Layer2(i::Int,o::Int,f=relu; pdrop=0) = Layer2(param(o,i),param0(o),f,pdrop)
(l::Layer2)(x) = l.f.(l.w * dropout(x,l.pdrop) .+ l.b)

# 1.06e-01  100.00%┣████████████████┫ 60000/60000 [00:53/00:53, 1142.58i/s]
# [0.0144843; 0.067272; 0.00376667; 0.0189]
model = Chain(Layer2(784,64,pdrop=0.2), Layer2(64,10,identity))
mlp4 = trainresults("mlp113e.jld2", model);

# overfitting less, loss results improve 0.0841 -> 0.0673
plot([mlp2[1,:], mlp2[2,:], mlp4[1,:], mlp4[2,:]], ylim=(0.0,0.15),
     labels=[:trnMLP :tstMLP :trnDrop :tstDrop],xlabel="Epochs",ylabel="Loss")

# this time error also improves 0.0243 -> 0.0189
plot([mlp2[3,:], mlp2[4,:], mlp4[3,:], mlp4[4,:]], ylim=(0.0,0.04),
     labels=[:trnMLP :tstMLP :trnDrop :tstDrop],xlabel="Epochs",ylabel="Error")

(mlperr=minimum(mlp2[4,:]),L1err=minimum(mlp3[4,:]),dropouterr=minimum(mlp4[4,:]))

(mlploss=minimum(mlp2[2,:]),L1loss=minimum(mlp3[2,:]),dropoutloss=minimum(mlp4[2,:]))

# ------------------------------------------------------------------------------------------
# ## MLP with larger hidden layer and dropout
# ------------------------------------------------------------------------------------------

# The current trend is to use models with higher capacity tempered with dropout
# 5.80e-02  100.00%┣████████████████┫ 60000/60000 [00:54/00:54, 1101.25i/s]
# [0.0035497; 0.0479644; 0.0003; 0.0149]
model = Chain(Layer2(784,256,pdrop=0.2), Layer2(256,10,identity))
mlp = trainresults("mlp113f.jld2", model);

# Both train and test loss is better with the larger model
plot([mlp4[1,:], mlp4[2,:], mlp[1,:], mlp[2,:]],ylim=(0,0.15),
    labels=[:trn64 :tst64 :trn256 :tst256],xlabel="Epochs",ylabel="Loss")

# We are down to 1.5% error.
plot([mlp4[3,:], mlp4[4,:], mlp[3,:], mlp[4,:]],ylim=(0,0.04),
    labels=[:trn64 :tst64 :trn256 :tst256],xlabel="Epochs",ylabel="Error")

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using BNP_WMReg_Joint
using Test

# @test hello("Julia") == "Hello, Julia"
# @test domath(2.0) ≈ 7.0

@test mean(InverseGamma(2.0, 2.0)) ≈ 1.0 # make sure they don't change parameterization
@test mean(Gamma(2.0, 2.0)) ≈ 4.0 # make sure they don't change parameterization

n = 500
K = 10
H = 50

y = randn(n)
X = randn(n, K)

prior = Prior_DPmRegJoint(1.0, # α_sh
1.0, # α_rate
randn(K+1), # β0star_ηy_mean
PDMat(Matrix(Diagonal(fill(10.0, K+1)))), # β0star_ηy_Cov
2.0*K, # Λ0star_ηy_df
PDMat(Matrix(Diagonal(fill(10.0, K+1)))), # Λ0star_ηy_S0
5.0, # s0_δy_df
1.0, # s0_δy_s0
randn(K), # μ0_μx_mean
PDMat(Matrix(Diagonal(fill(10.0, K)))), # μ0_μx_Cov
2.0*K, # Λ0_μx_df
PDMat(Matrix(Diagonal(fill(10.0, K)))), # Λ0_μx_S0
[randn(k) for k = (K-1):-1:1], # β0_βx_mean
[ PDMat(Matrix(Diagonal(fill(10.0, k)))) for k = (K-1):-1:1 ], # β0_βx_Cov
fill(2.0*K, K-1), # Λ0_βx_df
[ PDMat(Matrix(Diagonal(fill(10.0, k)))) for k = (K-1):-1:1 ], # Λ0_βx_S0
fill(5.0, K), # s0_δx_df
fill(1.0, K)) # s0_δx_s0

state = State_DPmRegJoint(randn(H), # μ_y,
randn(H, K), # β_y,
exp.(randn(H)), # δ_y,
randn(H, K), # μ_x,
[ randn(H, k) for k = (K-1):-1:1 ], # β_x,
exp.(randn(H, K)), # δ_x,
[ sample(Weights(ones(H))) for i = 1:n ], # S,
rDirichlet(ones(H), true), # lω,
1.0, # α,
randn(K+1), # β0star_ηy,
PDMat(Matrix(Diagonal(ones(K+1)))), # Λ0star_ηy,
5.0, # ν_δy,
1.0, # s0_δy,
randn(K), # μ0_μx,
PDMat(Matrix(Diagonal(ones(K)))), # Λ0_μx,
[ randn(k) for k = (K-1):-1:1 ], # β0_βx,
[ PDMat(Matrix(Diagonal(fill(10.0, k)))) for k = (K-1):-1:1 ], # Λ0_βx,
fill(5.0, K), # ν_δx,
ones(K), # s0_δx
[ PDMat(Matrix(Diagonal(fill(0.5, Int(K + K*(K+1)/2))))) for h = 1:H ], # cSig_ηx,
true # adapt
)

# several fuctions require this broadcasting behavior
broadcast(+, permutedims([1.0, 2.0]), zeros(3,2)) ≈ [1.0 2.0; 1.0 2.0; 1.0 2.0]

model = Model_DPmRegJoint(y, # y,
X, # X,
H, # H,
prior, # prior,
state # state
)

model.state.lNX = lNXmat(model.X, model.state.μ_x, model.state.β_x, model.state.δ_x)
@time lNXmat(model.X, model.state.μ_x, model.state.β_x, model.state.δ_x)
model.state.lωNX_vec = lωNXvec(model.state.lω, model.state.lNX)

@time update_alloc!(model, hcat(model.y, model.X))

ω_start = [0.2, 0.3, 0.25, 0.15, 0.10]
v = lω_to_v(log.(ω_start))

for i = 1:1000
    global lω = v_to_lω(v)
    global v = lω_to_v(lω)
end
exp.(lω) ≈ ω_start

@time update_vlω_mvSlice!(model)

Λβ0star_ηy = model.state.Λ0star_ηy * model.state.β0star_ηy
βΛβ0star_ηy = PDMats.quad(model.state.Λ0star_ηy, model.state.β0star_ηy)

h = 7
@time update_η_h_Met!(model, h, Λβ0star_ηy, βΛβ0star_ηy)

@time mcmc_DPmRegJoint!(model, 200, Monitor_DPmRegJoint(false, false, false), "out_progress.txt", 1, 10)
rm("out_progress.txt")

println(model.state.iter)
println(model.state.accpt / float(model.state.iter))
println(counts(model.state.S, 1:model.H))
println(exp.(model.state.lω))
println(model.state.α)

## Check sqfChol_to_Σ and lNX_sqfChol
mu = [1.0, 1.5, 2.0]
delta = [3.0, 4.0, 1.0]
betas = [ [0.0, 0.5], [2.5] ]

Sig = sqfChol_to_Σ(betas, delta)
Sig.mat

Sig.mat[3,2] ≈ -delta[3]*betas[2][1]
Sig.mat[3,1] ≈ delta[3]*(betas[1][1]*betas[2][1]-betas[1][2])
Sig.mat[2,2] ≈ delta[2] + delta[3]*betas[2][1]^2
Sig.mat[2,1] ≈ -delta[2]*betas[1][1] - delta[3]*betas[2][1]*(betas[1][1]*betas[2][1]-betas[1][2])
Sig.mat[1,1] ≈ delta[1] + delta[2]*betas[1][1]^2 + delta[3]*(betas[1][1]*betas[2][1]-betas[1][2])^2

X = [0.0, 2.0, 1.2]

lNX_sqfChol(X, mu, betas, delta) ≈ -0.5*length(mu)*log(2π) - 0.5*logdet(Sig) -0.5*( (X - mu)' * ( Sig \ (X - mu) ) )

XX = hcat(X, [-1.5, -1.0, 0.0])

lNX_sqfChol(XX, mu, betas, delta)


text:
"""
## sample!

This subroutine returns the local coordinates and weighting coefficients
of the integrating points.

### Function
```julia
sample!(fin_el, s, wt)
```

### Arguments
```julia
* fin_el::FiniteElement      : Finite element type
* s::Matrix{Float64}        : Local coordinates (Updated)
* wt::Vector{Float64}       : Weighting coefficients (Updated)
```
"""
function sample!(fin_el::Line, s::Matrix{Float64} , wt::Vector{Float64})
  #
  # This subroutine returns the local coordinates and weighting coefficients
  # of the integrating points.
  #
  nip = size(s,1)
  (s[:,1], wt[:]) = QuadGK.gauss(Float64, nip)
end


text:
# to run tests, type in Julia: Pkg.test("EmACQR")

include("../src/ACQR_test.jl")


text:
"""
$(TYPEDEF)

A hybrid noise model with both low and high frequency noise. The high frequency noise is characterized by Ohmic bath and the low frequence noise is characterized by the MRT width `W`.

$(FIELDS)
"""
struct HybridOhmicBath <: AbstractBath
    "MRT width (2π GHz)"
    W::Float64
    "low spectrum reorganization energy (2π GHz)"
    ϵl::Float64
    "strength of high frequency Ohmic bath"
    η::Float64
    "cutoff frequency"
    ωc::Float64
    "inverse temperature"
    β::Float64
end

function Base.show(io::IO, ::MIME"text/plain", m::HybridOhmicBath)
    print(
        io,
        "Hybrid Ohmic bath instance:\n",
        "W (mK): ",
        freq_2_temperature(m.W / 2 / pi),
        "\n",
        "ϵl (GHz): ",
        m.ϵl / 2 / pi,
        "\n",
        "η (unitless): ",
        m.η / 2 / π,
        "\n",
        "ωc (GHz): ",
        m.ωc / pi / 2,
        "\n",
        "T (mK): ",
        β_2_temperature(m.β)
    )
end

"""
    HybridOhmic(W, η, fc, T)

Construct HybridOhmicBath object with parameters in physical units. `W`: MRT width (mK); `η`: interaction strength (unitless); `fc`: Ohmic cutoff frequency (GHz); `T`: temperature (mK).
"""
function HybridOhmic(W, η, fc, T)
    # scaling of W to the unit of angular frequency
    W = 2 * pi * temperature_2_freq(W)
    # scaling of η because different definition of Ohmic spectrum
    η = 2 * π * η
    β = temperature_2_β(T)
    ωc = 2 * π * fc
    ϵl = W^2 * β / 2
    HybridOhmicBath(W, ϵl, η, ωc, β)
end

"""
    Gₕ(ω, bath::HybridOhmicBath)

High frequency noise spectrum of the HybridOhmicBath `bath`.
"""
function Gₕ(ω, bath::HybridOhmicBath)
    η = bath.η
    S0 = η / bath.β
    if isapprox(ω, 0, atol=1e-8)
        return 4 / S0
    else
        γ² = (S0 / 2)^2
        return η * ω * exp(-abs(ω) / bath.ωc) / (1 - exp(-bath.β * ω)) / (ω^2 + γ²)
    end
end

"""
    Gₗ(ω, bath::HybridOhmicBath)

Low frequency noise specturm of the HybridOhmicBath `bath`.
"""
function Gₗ(ω, bath::HybridOhmicBath)
    W² = bath.W^2
    ϵ = bath.ϵl
    sqrt(π / 2 / W²) * exp(-(ω - 4ϵ)^2 / 8 / W²)
end

function spectrum(ω, bath::HybridOhmicBath)
    Gl(x) = Gₗ(x, bath)
    Gh(x) = Gₕ(x, bath)
    integrand(x) = Gl(ω - x) * Gh(x)
    a, b = sort([0.0, bath.ϵl])
    res, err = quadgk(integrand, -Inf, a, b, Inf)
    res / 2 / π
end


text:
const idle_workers = Int[]
const workqueue = Channel(Inf)

function worker_loop()
    while true
        try
            id = myid()
            fetch(@spawnat 1 begin
                item = take!(FemtoCleaner.workqueue)
                println("Beginning work on worker ", id)
                item
            end)()
        catch e
            bt = catch_backtrace()
            Base.showerror(STDERR, e, bt)
        end
    end
end

queue(work) = nprocs() == 1 ? work() : push!(workqueue, work)


text:
function minimize!(objective, x0, scheme::QuasiNewton, B0=nothing, options=OptOptions(), cache=preallocate_qn_caches_inplace(x0))
    minimize!(objective, x0, (scheme, BackTracking()), B0, options, cache )
end
function minimize!(objective, x0, approach::Tuple{<:Any, <:LineSearch}, B0=nothing,
                   options::OptOptions=OptOptions(), # options
                   cache = preallocate_qn_caches_inplace(x0), # preallocate arrays for QN
                   )

    x, fx, ∇fx, z, fz, ∇fz, B = prepare_variables(objective, approach, x0, copy(x0), B0)

    # first iteration
    x, fx, ∇fx, z, fz, ∇fz, B, is_converged = iterate!(cache, x, fx, ∇fx, z, fz, ∇fz, B, approach, objective, options)

    iter = 0
    while iter <= options.maxiter && !is_converged
        iter += 1

        # take a step and update approximation
        x, fx, ∇fx, z, fz, ∇fz, B, is_converged = iterate!(cache, x, fx, ∇fx, z, fz, ∇fz, B, approach, objective, options, false)
    end
    return z, fz, ∇fz, iter
end

function iterate!(cache, x, fx, ∇fx, z, fz, ∇fz, B, approach::Tuple{<:Any, <:LineSearch}, objective, options, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch = approach

    y, d, s = cache.y, cache.d, cache.s

    # Move nexts into currs
    fx = fz
    copyto!(x, z)
    copyto!(∇fx, ∇fz)

    # Update current gradient and calculate the search direction
    d = find_direction!(d, B, ∇fx, scheme) # solve Bd = -∇f

    # Perform line search along d
    α, f_α, ls_success = linesearch(objective, d, x, fx, ∇fx, 1.0)

    # Calculate final step vector and update the state
    @. s = α * d
    @. z = x + s

    # Update approximation
    fz, ∇fz, B = update_obj!(objective, d, s, y, ∇fx, z, ∇fz, B, scheme, is_first)

    # Check for gradient convergence
    is_converged = converged(z, ∇fz, options.g_tol)

    return x, fx, ∇fx, z, fz, ∇fz, B, is_converged
end


text:
"""
`init_subspec!(m::Model990)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model990)
    if subspec(m) == "ss2"
        return
    elseif subspec(m) == "ss5"
        return ss5!(m)
    else
        error("This subspec is not defined.")
    end
end

"""
`ss5(m::Model990)`

Initializes subspecification 5 for Model990. Specifically, fixes ι_w
and ι_p to 0 (so that intermediate goods producers who do not readjust
prices and wages in a given period do not index to inflation.)
"""
function ss5!(m::Model990)

    m <= parameter(:ι_p, 0.0, fixed=true,
                   description= "ι_p: The persistence of last period's inflation in
                   the equation that describes the intertemporal
                   change in prices for intermediate goods producers
                   who cannot adjust prices. The change in prices is a
                   geometric average of steady-state inflation
                   (π_star, with weight (1-ι_p)) and last period's
                   inflation (π_{t-1})).",
                   tex_label="\\iota_p")

    m <= parameter(:ι_w,   0.0, fixed=true,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")
end


text:
__precompile__()
module Pidfile

export mkpidlock

using Base:
    IOError, UV_EEXIST, UV_ESRCH,
    Process

using Base.Filesystem:
    File, open, JL_O_CREAT, JL_O_RDWR, JL_O_RDONLY, JL_O_EXCL,
    samefile

using FileWatching: watch_file
using Base.Sys: iswindows

"""
    mkpidlock([f::Function], at::String, [pid::Cint, proc::Process]; kwopts...)

Create a pidfile lock for the path "at" for the current process
or the process identified by pid or proc. Can take a function to execute once locked,
for usage in `do` blocks, after which the lock will be automatically closed. If the lock fails
and `wait` is false, then an error is thrown.

Optional keyword arguments:
 - `mode`: file access mode (modified by the process umask). Defaults to world-readable.
 - `poll_interval`: Specify the maximum time to between attempts (if `watch_file` doesn't work)
 - `stale_age`: Delete an existing pidfile (ignoring the lock) if its mtime is older than this.
     The file won't be deleted until 25x longer than this if the pid in the file appears that it may be valid.
     By default this is disabled (`stale_age` = 0), but a typical recommended value would be about 3-5x an
     estimated normal completion time.
 - `wait`: If true, block until we get the lock, if false, raise error if lock fails.
"""
function mkpidlock end

mutable struct LockMonitor
    path::String
    fd::File

    global function mkpidlock(at::String, pid::Cint; kwopts...)
        local lock
        at = abspath(at)
        fd = open_exclusive(at; kwopts...)
        try
            write_pidfile(fd, pid)
            lock = new(at, fd)
            finalizer(close, lock)
        catch ex
            close(fd)
            rm(at)
            rethrow(ex)
        end
        return lock
    end
end

mkpidlock(at::String; kwopts...) = mkpidlock(at, getpid(); kwopts...)
mkpidlock(f::Function, at::String; kwopts...) = mkpidlock(f, at, getpid(); kwopts...)

function mkpidlock(f::Function, at::String, pid::Cint; kwopts...)
    lock = mkpidlock(at, pid; kwopts...)
    try
        return f()
    finally
        close(lock)
    end
end

# TODO: enable this when we update libuv
#Base.getpid(proc::Process) = ccall(:uv_process_get_pid, Cint, (Ptr{Void},), proc.handle)
#function mkpidlock(at::String, proc::Process; kwopts...)
#    lock = mkpidlock(at, getpid(proc))
#    @schedule begin
#        wait(proc)
#        close(lock)
#    end
#    return lock
#end

"""
    write_pidfile(io, pid)

Write our pidfile format to an open IO descriptor.
"""
function write_pidfile(io::IO, pid::Cint)
    print(io, "$pid $(gethostname())")
end

"""
    parse_pidfile(file::Union{IO, String}) => (pid, hostname, age)

Attempt to parse our pidfile format,
replaced an element with (0, "", 0.0), respectively, for any read that failed.
"""
function parse_pidfile(io::IO)
    fields = split(read(io, String), ' ', limit = 2)
    pid = tryparse(Cuint, fields[1])
    pid === nothing && (pid = Cuint(0))
    hostname = (length(fields) == 2) ? fields[2] : ""
    when = mtime(io)
    age = time() - when
    return (pid, hostname, age)
end

function parse_pidfile(path::String)
    try
        existing = open(path, JL_O_RDONLY)
        try
            return parse_pidfile(existing)
        finally
            close(existing)
        end
    catch ex
        isa(ex, EOFError) || isa(ex, IOError) || rethrow(ex)
        return (Cuint(0), "", 0.0)
    end
end

"""
    isvalidpid(hostname::String, pid::Cuint) :: Bool

Attempt to conservatively estimate whether pid is a valid process id.
"""
function isvalidpid(hostname::AbstractString, pid::Cuint)
    # can't inspect remote hosts
    (hostname == "" || hostname == gethostname()) || return true
    # pid < 0 is never valid (must be a parser error or different OS),
    # and would have a completely different meaning when passed to kill
    !iswindows() && pid > typemax(Cint) && return false
    # (similarly for pid 0)
    pid == 0 && return false
    # see if the process id exists by querying kill without sending a signal
    # and checking if it returned ESRCH (no such process)
    return ccall(:uv_kill, Cint, (Cuint, Cint), pid, 0) != UV_ESRCH
end

"""
    stale_pidfile(path::String, stale_age::Real) :: Bool

Helper function for `open_exclusive` for deciding if a pidfile is stale.
"""
function stale_pidfile(path::String, stale_age::Real)
    pid, hostname, age = parse_pidfile(path)
    age < -stale_age && @warn "filesystem time skew detected" path=path
    if age > stale_age
        if (age > stale_age * 25) || !isvalidpid(hostname, pid)
            return true
        end
    end
    return false
end

"""
    tryopen_exclusive(path::String, mode::Integer = 0o444) :: Union{Void, File}

Try to create a new file for read-write advisory-exclusive access,
return nothing if it already exists.
"""
function tryopen_exclusive(path::String, mode::Integer = 0o444)
    try
        return open(path, JL_O_RDWR | JL_O_CREAT | JL_O_EXCL, mode)
    catch ex
        (isa(ex, IOError) && ex.code == UV_EEXIST) || rethrow(ex)
    end
    return nothing
end

"""
    open_exclusive(path::String; mode, poll_interval, stale_age) :: File

Create a new a file for read-write advisory-exclusive access.
If `wait` is `false` then error out if the lock files exist
otherwise block until we get the lock.

For a description of the keyword arguments, see [`mkpidlock`](@ref).
"""
function open_exclusive(path::String;
                        mode::Integer = 0o444 #= read-only =#,
                        poll_interval::Real = 10 #= seconds =#,
                        wait::Bool = true #= return on failure if false =#,
                        stale_age::Real = 0 #= disabled =#)
    # fast-path: just try to open it
    file = tryopen_exclusive(path, mode)
    file === nothing || return file
    if !wait
        if file === nothing && stale_age > 0
            if stale_age > 0 && stale_pidfile(path, stale_age)
                @warn "attempting to remove probably stale pidfile" path=path
                try
                    rm(path)
                catch ex
                    isa(ex, IOError) || rethrow(ex)
                end
            end
            file = tryopen_exclusive(path, mode)
        end
        if file === nothing
            error("Failed to get pidfile lock for $(repr(path)).")
        else
            return file
        end
    end
    # fall-back: wait for the lock

    while true
        # start the file-watcher prior to checking for the pidfile existence
        t = @async try
            watch_file(path, poll_interval)
        catch ex
            isa(ex, IOError) || rethrow(ex)
            sleep(poll_interval) # if the watch failed, convert to just doing a sleep
        end
        # now try again to create it
        file = tryopen_exclusive(path, mode)
        file === nothing || return file
        Base.wait(t) # sleep for a bit before trying again
        if stale_age > 0 && stale_pidfile(path, stale_age)
            # if the file seems stale, try to remove it before attempting again
            # set stale_age to zero so we won't attempt again, even if the attempt fails
            stale_age -= stale_age
            @warn "attempting to remove probably stale pidfile" path=path
            try
                rm(path)
            catch ex
                isa(ex, IOError) || rethrow(ex)
            end
        end
    end
end

"""
    close(lock::LockMonitor)

Release a pidfile lock.
"""
function Base.close(lock::LockMonitor)
    isopen(lock.fd) || return false
    havelock = samefile(stat(lock.fd), stat(lock.path))
    close(lock.fd)
    if havelock # try not to delete someone else's lock
        rm(lock.path)
    end
    return havelock
end

end # module


text:
#!/usr/local/bin/julia
# ARGS=["REBC-AC8L-TP","REBC-AC8L-NT","tmp1.tsv","REBC-AC8L-TP-NT.M2_maflite.tsv"]

using DataFrames
tumor_id=ARGS[1]
normal_id=ARGS[2]
file1=ARGS[3]
maflite=ARGS[4]
file1a="vcf_1.tsv"

df = readtable(file1)
# size(df)
# describe(df)
# print(df)
delete!(df, [:FILTER,:QUAL, :NORMAL_PID,:x,:TUMOR_GT,:TUMOR_PL,:TUMOR_GQ,:TUMOR_DP,:TUMOR_PGT,:NORMAL_FOXOG,:NORMAL_QSS,:NORMAL_ALT_F2R1,:NORMAL_ALT_F1R2,:NORMAL_REF_F2R1,:NORMAL_REF_F1R2,:NORMAL_GT,:NORMAL_PGT,:NORMAL_GQ,:NORMAL_PL,:NORMAL_DP])
if ( :TUMOR_FOXOG in names(df) )  # FOXOG messed up by VariantAnnotator
    delete!(df, [:TUMOR_FOXOG])
end
# print(df)

for c in names(df)
    println(c)
    k=isna(df[c])
    if mean(k)==1.0
		delete!(df, c)
		continue
	end
    # if isa(df[c],DataArray{Float64,1})  & (mean(k)>0)
    #     v=df[c]
    #     v=map(x -> replace(x,NA,NaN), v)
    #     df[c]=v
    # end
    if isa(df[c],DataArray{String,1})
        df[k,c] = ""
    end
end
#rename!(df, [:normal_tumor_alt_count, :normal_tumor_ref_count], [:t_alt_count, :t_ref_count])
head(df)
a= df[:CHRO]
if !isa(a[1],Int)
    a=map(x -> replace(x,r"[X]", "23"), a)
    a=map(x -> replace(x,r"[Y]", "24"), a)
    a=map(x -> replace(x,r"[MT]", "25"), a)
    a=map(x -> replace(x,r"[M]", "25"), a)
    a=map(x -> parse(Int32,x), a)
end
df[:a] = a
# print(df)
sort!(df, cols = [:a, :POS])
delete!(df, [:a])
open(file1a, "w") do f
    writedlm(f, reshape(names(df), 1, length(names(df))), '\t')
    writedlm(f, convert(Array,df), '\t')
end

a=df[:ALT]
r=df[:REF]
n=length(a)
ar=fill("",n)
for i in 1:n
    #println(i)
    ar[i]=string(a[i],":",r[i])
end

l1=map(x -> length(x), ar)
df[:INDEL]=map(x -> x>3,l1)
df[:Variant_Type]=fill("SNV",n)
df[df[:INDEL],:Variant_Type]="INDEL"
dfindel=df[df[:INDEL],:]
dfsnv=df[~df[:INDEL],:]
#delete!(dfsnv, [:INDEL])
#delete!(dfindel, [:INDEL])
#delete!(df, [:INDEL])
df2=df

# set ALT0 to original ALT,shorten ALT to at most 5 bases, then copy back after merge

df2[:ALT0]=df2[:ALT]
alt0=df2[:ALT0]
alt=fill("",n)
n=length(alt0)
l0=map(x -> length(x), alt0)
for i in 1:n
    l1=min(l0[i],5)
    alt1=string(alt0[i])
    alt[i]=alt1[1:l1]
end
df2[:ALT]=alt

delete!(df, [:INDEL, :Variant_Type,:ALT0])

# label M1 with M1 flag
a= df[:CHRO]
if !isa(a[1],Int)
    a=map(x -> replace(x,r"[X]", "23"), a)
    a=map(x -> replace(x,r"[Y]", "24"), a)
    a=map(x -> replace(x,r"[MT]", "25"), a)
    a=map(x -> replace(x,r"[M]", "25"), a)
    a=map(x -> parse(Int32,x), a)
end
df[:a] = a
#print(df)
sort!(df, cols = [:a, :POS])
delete!(df, [:a])

isString(x::Number)=false
isString(x::DataArrays.NAtype)=false
isString(x::AbstractString)=true

for c in names(df)
     println(c)
    k=isna(df[c])
    if mean(k)==1.0
		delete!(df, c)
		continue
	end
    println(typeof(df[1,c]))
    if isString(df[1,c])
    	df[k,c] = ""
    end
end

# maflite fields
# SNP: build	chr	start	end	ref_allele	tum_allele1	tum_allele2	tumor_barcode	normal_barcode	tumor_f	init_t_lod	t_lod_fstar	t_alt_count	t_ref_count	judgement
# INDEL: build	contig	start_position	end_position	ref_allele	alt_allele1	alt_allele2	tumor_name	normal_name	tumor_f	n_ref_count	n_alt_count	t_ref_count	t_alt_count	init_n_lod	init_t_lod	judgement

df[:end]=df[:POS]

# start for indels is POS+1
kdel = map(x-> length(x)>1, df[:REF])
p=df[:POS]
p=p+kdel
df[:POS]=p
k = find(kdel)
ref=df[k,:REF]
ref=map(x->x[2:end],ref)
alt=df[k,:ALT]
alt=map(x->x[2:end],alt)
k1 = find(map(x-> length(x)<1, alt))
alt[k1]=map(x-> "-", alt[k1])
df[k,:REF]=ref
df[k,:ALT]=alt
dref = map(x-> length(x), ref)
dalt = map(x-> length(x), alt)
ddel = dref
df[k,:end]=df[k,:POS]+ddel

kins = map(x-> length(x)>1, df[:ALT])
# p=df[:POS]
# p=p+kins
# df[:POS]=p
# df[find(kins),:POS]
k = find(kins)
ref=df[k,:REF]
ref=map(x->x[2:end],ref)
k1 = find(map(x-> length(x)<1, ref))
ref[k1]=map(x-> "-", ref[k1])
alt=df[k,:ALT]
alt=map(x->x[2:end],alt)
df[k,:REF]=ref
df[k,:ALT]=alt
df[k,:end]=df[k,:POS]+1

# add build column
df[:build]=fill("37",size(df,1))
df[:tumor_barcode]=fill(tumor_id,size(df,1))
df[:normal_barcode]=fill(normal_id,size(df,1))
df[:judgement]=fill("KEEP",size(df,1))

# M1 counts
df[:n_alt_count]=df[:NORMAL_AD_ALT]
df[:n_ref_count]=df[:NORMAL_AD_REF]
df[:t_lod_fstar]=df[:TLOD]
df[:tumor_f]=df[:TUMOR_AF]
df[:t_alt_count]=df[:TUMOR_AD_ALT]
df[:t_ref_count]=df[:TUMOR_AD_REF]
rename!(df, [:CHRO,:POS, :REF, :ALT], [:chr, :start, :ref_allele,:alt_allele])#from = [:DB, :HCNT, :MAX_ED,:RPA,:STR, :RU, :NLOD,:TUMOR_ALT_F1R2,:MIN_ED,:ECNT, :TLOD, :TUMOR_PID,:TUMOR_AF,:TUMOR_FOXOG, :TUMOR_QSS, :TUMOR_ALT_F2R1,:TUMOR_AD_REF,:TUMOR_AD_ALT, :TUMOR_REF_F2R1, :TUMOR_REF_F1R2,:NORMAL_AF,:NORMAL_AD_REF,:NORMAL_AD_ALT]
#from = ["DB","HCNT","MAX_ED","RPA","STR","RU","NLOD","TUMOR_ALT_F1R2","MIN_ED","ECNT","TLOD","TUMOR_PID","TUMOR_AF","TUMOR_FOXOG","TUMOR_QSS","TUMOR_ALT_F2R1","TUMOR_AD_REF","TUMOR_AD_ALT","TUMOR_REF_F2R1","TUMOR_REF_F1R2","NORMAL_AF","NORMAL_AD_REF","NORMAL_AD_ALT"]
#println(from)

#for c in names(df)
#    c1 = string(c)
#    println(c1)
#    if length(find(Bool[contains(c1,i) for i in from]))>0
#        m2c = Symbol(string("M2_",c1))
#        rename!(df,c,m2c)
#    end
#end

maf=df
for c in names(maf)
    #println(c)
    if ~isString(maf[1,c])
        maf[c] = map(x -> string(x),maf[c])
    end
end

for c in names(maf)
    println(c)
    k=isna(maf[c])
    if mean(k)==1.0
        delete!(maf, c)
        continue
    end
    k=find(map(x->x=="NA",maf[c]))
    maf[k,c] = ""
    k=find(map(x->uppercase(x)=="TRUE",maf[c]))
    maf[k,c] = "1"
    k=find(map(x->uppercase(x)=="FALSE",maf[c]))
    maf[k,c] = "0"
end

open(maflite, "w") do f
    writedlm(f, reshape(names(maf), 1, length(names(maf))), '\t')
    writedlm(f, convert(Array,maf), '\t')
end


text:
# This file was generated, do not modify it. # hide
stand = machine(Standardizer(), Xs)
W = transform(stand, Xs)


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

@testset "mini-FL Tests" begin
   @testset "comp" begin
		@test Lar.comp([acos,atan,tan,sin])(pi/4) * 4 - pi < 0.1e-14
		@test Lar.comp([acos,atan,tan,sin])(1)==(acos∘atan∘tan∘sin)(1)
		@test Lar.comp([atan,tan])(pi/4)*4 > 3.14159265358979
   end

   @testset "cons Tests" begin
		@test ((Lar.cons([acos,atan,tan,sin])(pi/4) .<=
		[0.667458,0.665774,1.0,0.707107]))==[true,true,true,true]
		@test (Lar.cons([sin,cos])(pi/2).>=[1.0,0.0])==[true,true]
		@test Lar.cons([+,-])([2,3])==[[2,3],[-2,-3]]
		@test typeof(Lar.cons([+,-])([2,3]))==Array{Array{Int64,1},1}
   end

   @testset "k Tests" begin
		@test Lar.k(10)(100)==10
		@test Lar.k(10)(sin)==10
		@test Lar.k(sin)(cos)==sin
		@test Lar.k(sin)(1000)==sin
		@test Lar.k([1,2,3])(100)==[1,2,3]
		@test Lar.k([1,2,3])(sin)==[1,2,3]
   end

   @testset "aa Tests" begin
		@test (Lar.aa(sin)([0.,pi/2])==[0.0,1.0])
		@test typeof(Lar.aa(sin)([0.,pi/6,pi/2]))==Array{Float64,1}
		@test Lar.aa(sum)([[1,2],[3,4],[5,6],[7,8]])==[3,7,11,15]
   end
end


text:
# const SURFACE_BUFFER_FACTOR = 0/4 # abandoned for now

# """
#     Cavity(boundary,lattice,[dielectric,pump];align=false,[name]) -> domain
#     Cavity(lattice,boundary,[dielectric,pump];align=false,[name]) -> domain
# """
# Cavity(args...;kwargs...) = Domain(:Cavity,args...;kwargs...)
# """
#     Resonator(boundary,lattice,[dielectric,pump];align=false,[name]) -> domain
#     Resonator(lattice,boundary,[dielectric,pump];align=false,[name]) -> domain
# """
# Resonator(args...;kwargs...) = Domain(:Resonator,args...;kwargs...)
# """
#     Void(boundary,lattice;align=false,[name]) -> domain
#     Void(lattice,boundary;align=false,[name]) -> domain
# """
# Void(args...; align=false, kwargs...) = Domain(:Void,args...; align=align, kwargs...)
# """
#     Dielectric(lattice,boundary,[ε=1,F=0];align=false,[name]) -> domain
#     Dielectric(boundary,lattice,[ε=1,F=0];align=false,[name]) -> domain
# """
# Dielectric(args...;kwargs...) = Domain(:Dielectric,args...;kwargs...)
# Dielectric(bnd::Union{Boundary,Lattice},lat::Union{Boundary,Lattice},ε::Number,args...;kwargs...) = Domain(:Dielectric,bnd,lat,DielectricFunction(ε),args...;kwargs...)
# Dielectric(bnd::Union{Boundary,Lattice},lat::Union{Boundary,Lattice},ε,F::Number;kwargs...) = Domain(:Dielectric,bnd,lat,ε,PumpFunction(F);kwargs...)
# """
#     Resonator(boundary,lattice,[dielectric,pump];align=true,[name]) -> domain
#     Resonator(lattice,boundary,[dielectric,pump];align=true,[name]) -> domain
# """
# Waveguide(args...; kwargs...) = Domain(:Waveguide,args...; align=true, kwargs...)
# """
#     Lead(boundary,lattice,[dielectric,pump];align=true,[name]) -> domain
#     Lead(lattice,boundary,[dielectric,pump];align=true,[name]) -> domain
# """
# Lead(args...; kwargs...) = Domain(:Lead,args...; align=true, kwargs...)

# # polar spanning indices
# function generate_spanning_indices_polar(bnd::Boundary, lat::Lattice)
#     @assert typeof(bnd.shape)<:Union{Circle,Annulus} "polar coordinates only for circular outermost region"
#     @assert lat.x0==bnd.shape.x0 && lat.y0==bnd.shape.y0 "polar lattice origin $((lat.x0,lat.y0)) must be same as circle origin $((bnd.shape.x0,bnd.shape.y0))"
#     # i,j = get_lattice_index(lat,bnd.shape.x0,bnd.shape.y0)
#     # imin, imax = floor(Int,i), ceil(Int,i)
#     # jmin, jmax = floor(Int,j), ceil(Int,j)
#     # imin ≤ 0 ? imin = 0 : nothing
#     # jmin ≤ 0 ? jmin = 0 : nothing
#     imin = 0
#     imax = 1
#     jmin = 0
#     jmax = floor(Int,2π/lat.dθ)
#
#     flag_rmin = flag_rmax = flag_θmin = flag_θmax = true
#     init_flag = false
#     while flag_rmin || flag_rmax || flag_θmin|| flag_θmax
#         I = imin:imax
#         J = jmin:jmax
#         XYBool = BitArray(undef,2length(I)+2length(J))
#         count1 = 0
#         count2 = 0
#         temp = Array{Float64}(undef,2)
#         for i ∈ I
#             count2 += 1
#             X, Y = lat(i,J[1])
#             XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
#         end
#         (!init_flag || any(XYBool[(count1+1):count2])) && jmin≥1 ? (jmin -= 1; flag_θmin=true) : flag_θmin=false
#         !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
#         count1 = count2
#         for i ∈ I
#             count2 += 1
#             X,Y = lat(i,J[end])
#             XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
#         end
#         (!init_flag || any(XYBool[(count1+1):count2])) && (jmax+1)*lat.dθ≤2π ? (jmax += 1; flag_θmax=true) : flag_θmax=false
#         !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
#         count1 = count2
#         for j ∈ J
#             count2 += 1
#             X,Y = lat(I[1],j)
#             XYBool[1+count1+j-J[1]] = bnd.shape(X,Y)
#         end
#         (!init_flag || any(XYBool[(count1+1):count2])) && imin≥1 ? (imin -= 1; flag_rmin=true) : flag_rmin=false
#         !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
#         count1 = count2
#         for j ∈ J
#             count2 += 1
#             X,Y = lat(I[end],j)
#             XYBool[2length(I) + length(J) + j+1-J[1]] = bnd.shape(X,Y)
#         end
#         !init_flag || any(XYBool[(count1+1):count2]) ? (imax += 1; flag_rmax=true) : flag_rmax=false
#         !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
#     end
#     return imin, imax, jmin, jmax
# end

# function generate_xy(bnd::Boundary,lat::Lattice,imin,imax,jmin,jmax)
#     I, J = imin:imax, jmin:jmax
#     X = Array{Float64}(undef,length(I),length(J))
#     Y = Array{Float64}(undef,length(I),length(J))
#     for k ∈ CartesianIndices(X)
#         i,j = k[1],k[2]
#         X[k],Y[k] = lat(I[i],J[j])
#     end
#     return X, Y
# end

# function generate_interior(bnd::Boundary,lat::Lattice,x::Array,y::Array)
#     interior = BitArray(undef,size(x)...)
#     for k ∈ CartesianIndices(x)
#         interior[k] = bnd.shape(x[k],y[k])
#     end
#     return interior
# end

# function generate_translations(lat::Lattice,imin,imax,jmin,jmax)
#     I, J = imin:imax, jmin:jmax
#     x1 = Array{Float64}(undef,length(I),length(J))
#     x2 = Array{Float64}(undef,length(I),length(J))
#     x3 = Array{Float64}(undef,length(I),length(J))
#     x4 = Array{Float64}(undef,length(I),length(J))
#     y1 = Array{Float64}(undef,length(I),length(J))
#     y2 = Array{Float64}(undef,length(I),length(J))
#     y3 = Array{Float64}(undef,length(I),length(J))
#     y4 = Array{Float64}(undef,length(I),length(J))
#     for k ∈ CartesianIndices(x1)
#         i,j = k[1],k[2]
#         x1[k],y1[k] = lat(I[i]-1,J[j])
#         x2[k],y2[k] = lat(I[i]+1,J[j])
#         x3[k],y3[k] = lat(I[i],J[j]-1)
#         x4[k],y4[k] = lat(I[i],J[j]+1)
#     end
#     return (x1,y1),(x2,y2),(x3,y3),(x4,y4)
# end

function generate_surface!(interior::BitArray,lat::Lattice,bnd::Boundary,t1,t2,t3,t4,x,y)
    surface = falses(size(interior)...)
    surface_temp = BitArray(undef,size(interior)...)
    for k ∈ CartesianIndices(interior)
        surface_temp[k] = interior[k] &&
            (!bnd.shape(t1[1][k],t1[2][k]) ||
             !bnd.shape(t2[1][k],t2[2][k]) ||
             !bnd.shape(t3[1][k],t3[2][k]) ||
             !bnd.shape(t4[1][k],t4[2][k])
            )
        if surface_temp[k]
            # n,t,d = normal_distance(bnd.shape,x[k],y[k])
            # if d[1]<min(lat.dx*SURFACE_BUFFER_FACTOR,lat.dy*SURFACE_BUFFER_FACTOR)
                # interior[k] = false
                # surface[k] = false
                # surface[k+CartesianIndex( 1, 0)]=true
                # surface[k+CartesianIndex(-1, 0)]=true
                # surface[k+CartesianIndex( 0, 1)]=true
                # surface[k+CartesianIndex( 0,-1)]=true
            # else
                surface[k] = true
            # end
        end
    end
    return surface
end

function generate_neighbor_indices(interior,lat::Lattice)
    if lat.type==:Cartesian
        return generate_neighbor_indices_cartesian(interior)
    elseif lat.type==:Polar
        return generate_neighbor_indices_polar(interior)
    else
        throw(LatticeError(lat))
    end
end

function generate_neighbor_indices_cartesian(interior)
    NN = LinearIndices(interior) - reshape(cumsum(.!interior[:]),size(interior)...)
    nnxm = vcat(zeros(Int,1,size(interior,2)),NN[1:end-1,:])
    nnxp = vcat(NN[2:end,:],zeros(Int,1,size(interior,2)))
    nnym = hcat(zeros(Int,size(interior,1),1),NN[:,1:end-1])
    nnyp = hcat(NN[:,2:end],zeros(Int,size(interior,1),1))
    return nnxm,nnxp,nnym,nnyp
end

function generate_neighbor_indices_polar(interior)
    NN = LinearIndices(interior) - reshape(cumsum(.!interior[:]),size(interior)...)
    nnrm = vcat(zeros(Int,1,size(interior,2)),NN[1:end-1,:])
    nnrp = vcat(NN[2:end,:],zeros(Int,1,size(interior,2)))
    nnθm = hcat(NN[:,end],NN[:,1:end-1])
    nnθp = hcat(NN[:,2:end],NN[:,1])
    return nnrm,nnrp,nnθm,nnθp
end

function generate_corner(x,y,bnd::Boundary,lat::Lattice,surface,interior)
    if lat.type==:Cartesian
        return generate_corner_cartesian(x,y,bnd,lat,surface,interior)
    elseif lat.type==:Polar
        return generate_corner_polar(x,y,bnd,lat,surface,interior)
    else
        throw(LatticeError(lat))
    end
end

function generate_corner_cartesian(x,y,bnd::Boundary,lat::Lattice,surface,interior)
    corner = falses(size(surface)...)
    c = bnd.shape.corners
    dx = hypot(lat.dx,lat.dy)#*(1+SURFACE_BUFFER_FACTOR)
    if length(c)>0
        for k ∈ CartesianIndices(surface)
            if surface[k]
                d = map(z->hypot(x[k]-z[1],y[k]-z[2]),c)
                corner[k] = minimum(d) ≤ dx+eps(dx)
                X1, Y1 = x[k+CartesianIndex( 2, 0)], y[k+CartesianIndex( 2, 0)]
                X2, Y2 = x[k+CartesianIndex(-2, 0)], y[k+CartesianIndex(-2, 0)]
                X3, Y3 = x[k+CartesianIndex( 0, 2)], y[k+CartesianIndex( 0, 2)]
                X4, Y4 = x[k+CartesianIndex( 0,-2)], y[k+CartesianIndex( 0,-2)]
                count = sum(
                            (!bnd.shape(X1,Y1) && !interior[k+CartesianIndex( 1, 0)],
                             !bnd.shape(X2,Y2) && !interior[k+CartesianIndex(-1, 0)],
                             !bnd.shape(X3,Y3) && !interior[k+CartesianIndex( 0, 1)],
                             !bnd.shape(X4,Y4) && !interior[k+CartesianIndex( 0,-1)],)
                            )
                corner[k] = corner[k] && (count>1)
            end
        end
    end
    return corner
end

function generate_corner_polar(x,y,bnd,lat,surface,interior)
    return falses(size(surface)...)
end