text:
#include(".../Cloudarray_time_test.jl")
#
#   DESCRIPTION: Tests the CloudArray construtor method in order to get it's cloud management execution time
#		 It will apply the function call parallel_reduce(+,darray) to the created carray
# 		 For a sequencial array (eg. 1:100) it should return an arithmetic progression ---> (n[1]+n[last])*n/2 which results 5050 in the 1:100 case
#
#
#       OPTIONS: ---
#  DEPENDENCIES: Infra.jl, DistributedArrays.jl, CloudArray.jl, PyPlot.jl, sshpass
#          BUGS: ---
#         NOTES: ---
#        AUTHOR: Naelson Douglas C. Oliveira
#  ORGANIZATION: GSD-UFAL
#       CREATED: 2015-10-09 05:32
###=============================================================================

#To use simple plots: https://github.com/nolta/Winston.jl

#include("CloudArray.jl")
#run(`source cloud_setup.sh`)
#include("../Infra.jl")
#include("../CloudArray.jl")

const _BYTE_ = 1
const _KBYTE_ = 1024_BYTE_
const _MBYTE_ = 1024_KBYTE_
const _GBYTE_ = 1024_MBYTE_
const _PBYTE_ = 1024_GBYTE_
###=========================User interface=======================================
#How many times you want the code to run on the same input?

#TODO passar essas variáveis para o usuário como parâmetro
#TODO criar um help

executions= 1
chunks = 1
#Input file. You also can use or generate an Array, just like this ---> input = [1:100]
input = "big_input.txt"

#The max memory size used by the DArray in every containner
chunk_size = 10_KBYTE_
#first test = 200KB
#Second test = 50KB

#The folder where the output will be store ---> format output_folder/input/'data_execution_'[git commit ID]_[current date/time]
output_folder = "tests_output"
###=============================================================================

#Case study 1: parallel map
#TODO: parallel map
#Case study 2: parallel reduce
parallel_reduce(f,darray) = reduce(f, map(fetch, { @spawnat p reduce(f, localpart(darray)) for p in workers()} ))

darray = dzeros(1)

#Output data treatment
commit = chomp(readall(`git log --format="%H" -n 1 `))
mkpath(string(output_folder,"/",input))
path  = string(output_folder,"/",input,"/","data_execution_",commit,"_",Dates.now(),".txt")
touch(path)
output_file = open(path,"a")

#Get the input file size and calculate the chunks usage
input_file = open(input)
input_size = stat(input_file).size
close(input_file)
chunk_size = Int64(round(input_size*2))

#It loops 'execution' times and in each iteration it uses as many chunks as the iteration number.
#eg in iteration 4 it will create a CloudArray wit 4 chunks and in iteration 2 it will create a carray with 2 chunks
#Note: for some values the number may be shifted by +-1
for i=1:chunks
	eval_time = Array(Float64,1)
	for j=1:executions

		println("___________________")
		println("Test loop with ",i," cores. Evaluation: ",j," out of ",executions,"")
		#println("\nTest call ",j," out of ",executions,"\n Starting carray using ",i," chunks")
		println("___________________")

		#creates a cloudarray with the given input
		carray  = DArray(input,chunk_size)
		#Get the amount of chunks in the new carray
		n_chunks = length(carray.chunks)
		println("CArray created with ",n_chunks," chunk(s) of ",chunk_size," for bytes each chunk\n")

		#Here is what the test happens: it runs the parallel_reduce function and calculates its execution time with tic() toc()
		tic()
		parallel_reduce(+,carray)
		#Does the evaluation and stores the time
		push!(eval_time,toc())
		#clear the variable
		carray = dzeros(1)
		@everywhere gc()

		if j==executions

			#Lower the chunk_size for each iteration to increase the chunks. The input is always the same(i.e. same size)

			if n_chunks > 1
				chunk_size = Int64(trunc(chunk_size - chunk_size/n_chunks))
			else
				chunk_size = Int64(trunc(chunk_size/2))
			end

			eval_mean = mean(eval_time)
			write(output_file,string(n_chunks,",",eval_mean,"\n"))
			flush(output_file)
			#we remove all containers and all data in order to free resources and start a new step in the loop
			delete_containers(all)
		end

	end
end

#flush and close the IO link to the log output
flush(output_file)
close(output_file)


text:
@testset "rook.jl" begin
    @printf("> rook\n")

	@test Bb.orthoSlide(p2u["a1"]) == Bb.MASK_RANKS[1] | Bb.MASK_FILES[1] ⊻ p2u["a1"] ⊻ p2u["a8"] ⊻ p2u["h1"]
    @test Bb.orthoSlide(p2u["h8"]) == Bb.MASK_RANKS[8] | Bb.MASK_FILES[8] ⊻ p2u["h8"] ⊻ p2u["a8"] ⊻ p2u["h1"]
    @test Bb.orthoSlide(p2u["e4"]) == (Bb.MASK_RANKS[4] | Bb.MASK_FILES[5] ⊻ p2u["e4"] ⊻
                                       p2u["e8"] ⊻ p2u["e1"] ⊻ p2u["a4"] ⊻ p2u["h4"])
    @test Bb.ROOK_MASKS[p2u["a1"]] == Bb.MASK_RANKS[1] | Bb.MASK_FILES[1] ⊻ p2u["a1"] ⊻ p2u["a8"] ⊻ p2u["h1"]
    @test Bb.ROOK_MASKS[p2u["h8"]] == Bb.MASK_RANKS[8] | Bb.MASK_FILES[8] ⊻ p2u["h8"] ⊻ p2u["a8"] ⊻ p2u["h1"]
    @test Bb.rookMasksGen()[p2u["e4"]] == (Bb.MASK_RANKS[4] | Bb.MASK_FILES[5] ⊻ p2u["e4"] ⊻
                                           p2u["e8"] ⊻ p2u["e1"] ⊻ p2u["a4"] ⊻ p2u["h4"])

    @test Bb.orthoAttack(p2u["a1"], p2u["a2"]) == Bb.MASK_RANKS[1] ⊻ p2u["a1"] | p2u["a2"]
    @test Bb.orthoAttack(p2u["a1"], Bb.EMPTY) == Bb.MASK_RANKS[1] | Bb.MASK_FILES[1] ⊻ p2u["a1"]
    @test Bb.orthoAttack(p2u["e3"], Bb.EMPTY) == Bb.MASK_RANKS[3] | Bb.MASK_FILES[5] ⊻ p2u["e3"]
end

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "distributions" begin
    include("test_distribution_functions.jl")
    include("test_standard_uniform.jl")
    include("test_standard_normal.jl")
    include("test_hierarchical_distribution.jl")
    include("test_log_uniform.jl")

    Test.@testset "multimodal student-t distribution" begin
        @test_throws ArgumentError BAT.MultimodalStudentT(n=1)
        for i in 2:6
            μ = randn()
            σ = rand(0.01:0.1:2)
            mmc = BAT.MultimodalStudentT(μ=μ,σ=σ, ν=1, n=i)
            @test isequal(mean(mmc), fill(NaN, i))
            @test isequal(var(mmc), fill(NaN, i))
            @test isequal(diag(cov(mmc)), fill(NaN, i))
            ps = StatsBase.params(mmc)
            @test sort(ps[1]) == sort(vcat([μ, -μ], zeros(i-2)))
            @test ps[2] == [σ for j in 1:i]
            @test ps[3] == i
        end
    end

    Test.@testset "funnel distribution" begin
        for i in 2:6
            a = randn()
            b = randn()
            funnel = BAT.FunnelDistribution(a=a, b=b, n=i)
            @test mean(funnel) == zeros(i)
            @test StatsBase.params(funnel) == (a, b, i)
        end
    end
end


text:
module KnetLayers

using Knet
import Knet: save, load, gc, rnnforw, rnninit
export gpu, KnetArray,
    relu, sigm, elu, invx, mat, bmm,
    logp, logsumexp, nll, bce, logistic, accuracy, zeroone, dropout, softmax,
    Data, minibatch,
    param, param0, params,
    gaussian, xavier, bilinear, setseed,
    SGD, Momentum, Nesterov, Adam, Adagrad, Adadelta, Rmsprop, update!, optimizers,
    train!, converge!, adam!, sgd!, nesterov!, rmsprop!, progress!,
    adam, converge, train, sgd, nesterov, rmsprop, progress,
    hyperband, goldensection, cpucopy, gpucopy,
    value, grad, cat1d, Param, @diff, @zerograd

"""
    KnetLayers.dir(path...)
Construct a path relative to KnetLayers root.
# Example
```julia
julia> KnetLayers.dir("src")
"/Users/ekin/git/KnetLayers/src"
```
"""
dir(path...) = joinpath(dirname(@__DIR__),path...)
seed! = Knet.seed!

arrtype = Array{Float32}
function __init__()
    global arrtype = gpu()>=0 ? KnetArray{Float32} : Array{Float32}
end
#Setters for atype
settype!(t::T) where T<:Type{KnetArray{V}} where V <: AbstractFloat = gpu()>=0 ? (global arrtype=t) : error("No GPU available")
settype!(t::T) where T<:Type{Array{V}} where V <: AbstractFloat = (global arrtype=t)
settype!(t::Union{Type{KnetArray},Type{Array}}) = settype!(t{Float32})
#Param init function
Prm(x)=Param(atype(x));   export Prm

include("core.jl");
include("primitive.jl");   export Multiply, Embed, Linear, Dense, BatchNorm
include("nonlinear.jl");   export ReLU,Sigm,Tanh,LeakyReLU,ELU,Dropout,LogSoftMax,SoftMax,LogSumExp
include("loss.jl");        export CrossEntropyLoss, BCELoss, LogisticLoss
include("cnn.jl");         export Pool,UnPool,DeConv,Conv
include("special.jl");     export MLP
include("rnn.jl");         export RNN,SRNN,LSTM,GRU,RNNOutput,PadRNNOutput,PadSequenceArray
include("chain.jl");       export Chain
include("../data/IndexedDict.jl");

end # module


text:
using Test
using QSimulator
import QSimulator: hamiltonian

function test_1()

  for i=1:100
      r1 = Resonator("r1", 1.0, 10)

      system = CompositeQSystem()

      system += r1

      state = randn(10) + 1im*randn(10)
      state = state/LinearAlgebra.norm(state,2)
      t = 1.0+10*abs(randn())
      #t = abs(randn())

      ans = state .* exp.(-1im*2*pi*t*Float64[k for k in 0:9])
      ev_state = unitary_evolution(state,system,.1,0.,t)

      @test isapprox(dot(ans,ev_state), 1.0, atol = 1e-12)
      #@test isapprox(norm(state, 2) 1.0 1e-12)
  end
end

function test_drive()

    # create the qubits
    q1 = Qubit("q1", 4.0)
    q2 = Qubit("q2", 4.5)

    # create a time dependant drive
    drive = MicrowaveControl("drive1", 4.0, timeStep = 1.0)

    # create the Hamiltonian with 'objects' and 'interactions'
    sys = q1 + q2 + FlipFlop(q1, q2, 0.005) + SemiClassicalDipole(Field(drive), q1, 0.01);

    seq = ones(200)
    load_sequence!(drive, seq);

    # test starting eigenvals
    E_test = LinearAlgebra.eigvals(hamiltonian(sys, 0.0))
    E = [-2.4999878473470395e-5,
            3.999975000433993,
            4.500024999566007,
            8.500024999878473]
    @test isapprox(E, E_test, atol = 1e-12)

    # define the initial state
    rhoIn = zeros(ComplexF64, dim(sys), dim(sys))
    rhoIn[1] = 1;
    #define the measurement projector
    evecs = LinearAlgebra.eigvecs(hamiltonian(sys))
    measOp = evecs[:,2] * evecs[:,2]'
    real(measOp)

    function run_simulation(sys, rhoIn, measOp, times)
    populations = zeros(length(times)-1)
    U = QSimulator.eye(ComplexF64, dim(sys))
    dt = 0.01
        for ct = 1:length(times)-1
            U = unitary_propagator(sys, dt, 0.0, times[ct+1])
            populations[ct] = real(LinearAlgebra.tr(U*rhoIn*U'*measOp))
        end
    return populations
    end

    times = 0.0:20.0
    populations_test = run_simulation(sys, rhoIn, measOp, times);
    populations = [0.001007229516897708,
                     0.003966998572908745,
                     0.008873938406129277,
                     0.01570878588050283,
                     0.024444712566132513,
                     0.03504742273249326,
                     0.047475299645262155,
                     0.06167955185968114,
                     0.0776044282644603,
                     0.09518740596921914,
                     0.11435947513329113,
                     0.1350453610888861,
                     0.15716387889776454,
                     0.1806281818148819,
                     0.20534618442196084,
                     0.23122082912331615,
                     0.2581505756295072,
                     0.28602967689427616,
                     0.31474873168055584,
                     0.3441949611994505]

    @test isapprox(populations, populations_test, atol = 1e-12)
end

@testset "unitary" begin
    test_1()
end

@testset "driven" begin
    test_drive()
end


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

@time begin
    import SparseArrays
    import Base.Threads: @threads, threadid

    # -------------------------------------------------------------------
    # Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in
    # the Julia Pkg REPL to install the package, then you must activate
    # the package enviroment (see README)
    import Chemostat_Heerden2013
    const ChHd = Chemostat_Heerden2013
    const Hd  = ChHd.HeerdenData
    const BD  = ChHd.BegData
    const iJR = ChHd.iJR904

    # -------------------------------------------------------------------
    # run add "https://github.com/josePereiro/Chemostat" in the
    # julia Pkg REPL for installing the package
    import Chemostat
    import Chemostat.LP.MathProgBase

    const Ch = Chemostat

    const ChU = Ch.Utils
    const ChSS = Ch.SteadyState
    const ChLP = Ch.LP
    const ChEP = Ch.MaxEntEP
    const ChSU = Ch.SimulationUtils

    import ChemostatPlots
    const ChP = ChemostatPlots

    import UtilsJL
    const UJL = UtilsJL

    using Statistics
    using ProgressMeter
    using Base.Threads
    using Serialization

    # -------------------------------------------------------------------
    using Plots, FileIO
    import GR
    GR.inline("png")
end

## ----------------------------------------------------------------------------
# ME data index
ME_INDEX = ChU.load_data(iJR.MAXENT_VARIANTS_ME_INDEX_FILE; verbose = false);

## ----------------------------------------------------------------------------
# LP data
LPDAT = ChU.load_data(iJR.LP_DAT_FILE; verbose = false)

# -------------------------------------------------------------------
const ME_Z_OPEN_G_OPEN          = :ME_Z_OPEN_G_OPEN           # Do not use extra constraints
const ME_MAX_POL                = :ME_MAX_POL                 # Match ME and Dy biom average and constraint av_ug
const ME_Z_EXPECTED_G_BOUNDED   = :ME_Z_EXPECTED_G_BOUNDED    # Match ME and Dy biom average and constraint av_ug
const ME_Z_EXPECTED_G_MOVING    = :ME_Z_EXPECTED_G_MOVING     #
const ME_Z_FIXXED_G_BOUNDED     = :ME_Z_FIXXED_G_BOUNDED      # Fix biom around observed

# -------------------------------------------------------------------
fileid = "2.1"
function mysavefig(p, pname; params...)
    fname = UJL.mysavefig(p, string(fileid, "_", pname), iJR.MODEL_FIGURES_DIR; params...)
    @info "Plotting" fname
end
myminmax(a::Vector) = (minimum(a), maximum(a))
CONC_IDERS = ["GLC", "SUCC", "AC", "FORM"]
FLX_IDERS = ["GLC", "SUCC", "AC", "FORM"]

ALL_METHODS = [
    # ME_Z_OPEN_G_OPEN,
    ME_MAX_POL,
    # ME_Z_FIXXED_G_BOUNDED,
    # ME_Z_EXPECTED_G_BOUNDED,
    # ME_Z_EXPECTED_G_MOVING
]

Hd_rxns_map = iJR.load_rxns_map()
Hd_mets_map = iJR.load_mets_map()

## ----------------------------------------------------------------------------
# collect DAT
include("2.2.1_collect_DAT.jl")

## ----------------------------------------------------------------------------
EXPS = DAT[:EXPS]

exp_colors = let
    colors = Plots.distinguishable_colors(length(EXPS))
    Dict(exp => color for (exp, color) in zip(EXPS, colors))
end

ider_colors = Dict(
    "GLC" => :red, "SUCC" => :yellow,
    "AC" => :orange, "FORM" => :blue,
    "D" => :black,
)

method_colors = Dict(
    ME_Z_OPEN_G_OPEN => :red,
    ME_MAX_POL => :blue,
    ME_Z_EXPECTED_G_BOUNDED => :orange,
    ME_Z_EXPECTED_G_MOVING => :purple,
    ME_Z_FIXXED_G_BOUNDED => :green,
)

## -------------------------------------------------------------------
# PLOTS
# -------------------------------------------------------------------

## -------------------------------------------------------------------
# polytope box volume
include("2.2.2_pol_vox_volume.jl")

## -------------------------------------------------------------------
# Var study
include("2.2.3_var_study.jl")

## -------------------------------------------------------------------
# MSE study
include("2.2.4_MSE_study.jl")

## -------------------------------------------------------------------
# proj 2D
include("2.2.5_proj2d_study.jl")

## -------------------------------------------------------------------
# Beta tendencies
include("2.2.6_Beta_tendencies.jl")

## -------------------------------------------------------------------
# flx correlations
include("2.2.7_flx_correlations.jl")

## -------------------------------------------------------------------
# inbounds study
include("2.2.7_flx_correlations.jl")

## -------------------------------------------------------------------
# marginals
include("2.2.9_marginals.jl")

## -------------------------------------------------------------------
# legends
include("2.2.10_leyends.jl")


text:
module API
using HDF5
import HPAT

operators = [:data_source_HDF5,:data_source_TXT,:Kmeans,:LinearRegression,:NaiveBayes]

for op in operators
    @eval export $op
end

@noinline function data_source_HDF5(T::DataType, var::AbstractString, file_name::AbstractString)
    arr::T = h5read(file_name, var)
    return arr
end

@noinline function data_sink_HDF5(y, var::AbstractString, file_name::AbstractString)
    h5write(file_name, var,y)
    return nothing
end

@noinline function data_source_TXT(T::DataType, file_name::AbstractString)
    arr::T = readdlm(file_name, ' ', eltype(T))
    return arr
end

@noinline function Kmeans{T}(points::Matrix{T}, numCenter::Int, iterNum::Int)
    # naive backup sequential implementation
    D = size(points,1) # number of features
    N = size(points,2) # number of instances
    centroids = rand(T, D, numCenter)
    for l in 1:iterNum
        dist = [ T[sqrt(sum((points[:,i].-centroids[:,j]).^2)) for j in 1:numCenter] for i in 1:N]
        labels = Int[indmin(dist[i]) for i in 1:N]
        centroids = T[ sum(points[j,labels.==i])/sum(labels.==i) for j in 1:D, i in 1:numCenter]
    end
    return centroids
end

@noinline function LinearRegression{T}(points::Matrix{T}, responses::Vector{T})
    # TODO: write actual sequential implementation
    # return random result to enable type inference
    return rand(T, size(points,1),2)
end

@noinline function NaiveBayes{T}(points::Matrix{T}, responses::Vector{T})
    # TODO: write actual sequential implementation
    # return random result to enable type inference
    return rand(T, size(points,1),2)
end

@noinline function join(t1::Vector{Vector},t2::Vector{Vector})
    t1_cols = length(t1)
    t2_cols = length(t2)
    # 1 key column for output
    out_size = t1_cols+t2_cols-1;
    out = [ [] for i in 1:out_size]
    # for each key in table 1
    for i in 1:length(t1[1])
        # for each key in table 2
        for j in 1:length(t2[1])
            if t1[1][i]==t2[1][j]
                # save key to output
                push!(out[1],t1[1][i])
                # save table 1 values
                for ii in 2:length(t1)
                    push!(out[ii], t1[ii][i])
                end
                # save table 2 values
                for jj in 2:length(t2)
                    push!(out[t1_cols+jj-1], t2[jj][j])
                end
            end
        end
    end
    return out
end

@noinline function aggregate{T}(t1k::Vector{T}, new_cols::Vector{ANY})
    new_cols_len = length(new_cols)
    # group values for each key into a separate array
    # key => [ [...] , [...] ...]
    groups = Dict{T,Vector{Vector}}()
    # for each key
    for i in 1:length(t1k)
        # initialize key if seen for first time
        if !haskey(groups,t1k[i])
            # one array for each output
            groups[t1k[i]] = Array(Vector,new_cols_len)
            for j in 1:new_cols_len
                groups[t1k[i]][j] = []
            end
        end
        # save group values
        for j in 1:new_cols_len
            push!(groups[t1k[i]][j], new_cols[j][1][i])
        end
    end
    #dprintln(groups)
    # output is array of column vectors
    out = Array(Vector,new_cols_len+1)
    # first array is keys
    out[1] = collect(keys(groups))
    for j in 1:new_cols_len
        out[j+1] = []
    end
    for j in 1:new_cols_len
        func = new_cols[j][2]
        #println(func)
        for key in out[1]
            push!(out[j+1], func(groups[key][j]))
        end
    end
    return out
end

@noinline function table_filter!(cond::BitArray{1}, columns::Vector{Vector})
    for i in 1:length(columns)
        columns[i] = columns[i][cond]
    end
end

#=
@doc """
function join{T1,T12,T22}(t1c1::Vector{T1}, t1c2::Vector{T12}, t2c1::Vector{T1}, t2c2::Vector{T22})
       out1 = t1c1 .+ t2c1
       out2 = t1c2
       return out1,out2
end

"""
# join up to 50 column tables
const MAX_COLUMNS = 50

for k in 1:MAX_COLUMNS
    for l in 1:MAX_COLUMNS
        # create type symbol names
        # T11 T12 ... T21 T22...
        typ_names1 = Symbol[Symbol("T1$j") for j in 1:k]
        typ_names2 = Symbol[Symbol("T2$j") for j in 1:l]
        typ_names = [typ_names1;typ_names2]
        # join{T1,T11,T12...T21,T22...}
        fexpr = Expr(:curly,:join,:T1,typ_names...)
        # t1c1::Vector{T11}, t1c2::Vector{T12} ....
        arrs1 = Expr[ Expr(:(::),Symbol("t1c$j"), Expr(:curly,:Vector,typ_names1[j])) for j in 1:k]
        arrs2 = Expr[ Expr(:(::),Symbol("t2c$j"), Expr(:curly,:Vector,typ_names2[j])) for j in 1:l]
        # key array: t1k::Vector{T1}
        arrk1 = Expr(:(::),Symbol("t1k"), Expr(:curly,:Vector,:T1))
        arrk2 = Expr(:(::),Symbol("t2k"), Expr(:curly,:Vector,:T1))
        arrs = [arrk1; arrs1;arrk2; arrs2]
        # fcall = Expr(:call, fexpr,arrs...)
        # out1=T1[], out2=T2[]
        out_decls1 = Expr[ Expr(:(=),Symbol("out1$j"), Expr(:ref,Symbol("T1$j"))) for j in 1:k]
        out_decls2 = Expr[ Expr(:(=),Symbol("out2$j"), Expr(:ref,Symbol("T2$j"))) for j in 1:l]
        outk_decl = Expr(:(=),:outk, Expr(:ref,:T1))
        out_decls = [outk_decl; out_decls1; out_decls2]
        out_block = Expr(:block, out_decls...)
        # list of output symbols: outk out11 out12... out21 out22...
        outs1 = Symbol[Symbol("out1$j") for j in 1:k]
        outs2 = Symbol[Symbol("out1$j") for j in 1:k]
        outs = [:outk;outs1;outs2]
        # push!(out11,t1c1[i])
        push_calls1 = [ Expr(:call,:push!,Symbol("out1$j"),Symbol("t1c$j[i]"))  for j in 1:k]
        # push!(out21,t2c1[j])
        push_calls2 = [ Expr(:call,:push!,Symbol("out2$j"),Symbol("t2c$j[j]"))  for j in 1:l]
        push_calls = Expr(:block, push_calls1..., push_calls2...)
        @eval begin
            @noinline function ($fexpr)($(arrs...))
                $out_block
                for i in 1:length(t1k)
                    for j in length(t2k)
                        if t1k[i]==t2k[j]
                            push!(outk,t1k[i])
                            $push_calls
                        end
                    end
                end
                return $(Expr(:tuple,outs...))
            end
        end

    end # for k
end # for l
=#

end # module


text:
using OrderedCollections, Test

@testset "OrderedSet" begin

    @testset "Constructors" begin
        @test isa(OrderedSet(), OrderedSet{Any})
        @test isa(OrderedSet([1,2,3]), OrderedSet{Int})
        @test isa(OrderedSet{Int}([3]), OrderedSet{Int})
        data_in = (1, "banana", ())
        s = OrderedSet(data_in)
        data_out = collect(s)
        @test isa(data_out, Array{Any,1})
        @test tuple(data_out...) === data_in
        @test tuple(data_in...) === tuple(s...)
        @test length(data_out) == length(data_in)
    end

    @testset "hash" begin
        s1 = OrderedSet{String}(["bar", "foo"])
        s2 = OrderedSet{String}(["foo", "bar"])
        s3 = OrderedSet{String}(["baz"])
        @test hash(s1) != hash(s2)
        @test hash(s1) != hash(s3)
    end

    @testset "isequal" begin
        @test  isequal(OrderedSet(), OrderedSet())
        @test !isequal(OrderedSet(), OrderedSet([1]))
        @test  isequal(OrderedSet{Any}(Any[1,2]), OrderedSet{Int}([1,2]))
        @test !isequal(OrderedSet{Any}(Any[1,2]), OrderedSet{Int}([1,2,3]))

        @test  isequal(OrderedSet{Int}(), OrderedSet{AbstractString}())
        @test !isequal(OrderedSet{Int}(), OrderedSet{AbstractString}([""]))
        @test !isequal(OrderedSet{AbstractString}(), OrderedSet{Int}([0]))
        @test !isequal(OrderedSet{Int}([1]), OrderedSet{AbstractString}())
        @test  isequal(OrderedSet{Any}([1,2,3]), OrderedSet{Int}([1,2,3]))
        @test  isequal(OrderedSet{Int}([1,2,3]), OrderedSet{Any}([1,2,3]))
        @test !isequal(OrderedSet{Any}([1,2,3]), OrderedSet{Int}([1,2,3,4]))
        @test !isequal(OrderedSet{Int}([1,2,3]), OrderedSet{Any}([1,2,3,4]))
        @test !isequal(OrderedSet{Any}([1,2,3,4]), OrderedSet{Int}([1,2,3]))
        @test !isequal(OrderedSet{Int}([1,2,3,4]), OrderedSet{Any}([1,2,3]))
    end

    @testset "eltype, empty" begin
        s1 = empty(OrderedSet([1,"hello"]))
        @test isequal(s1, OrderedSet())
        @test eltype(s1) === Any
        s2 = empty(OrderedSet{Float32}([2.0f0,3.0f0,4.0f0]))
        @test isequal(s2, OrderedSet())
        @test eltype(s2) === Float32
    end

    @testset "show" begin
        @test endswith(sprint(show, OrderedSet()), "OrderedSet{Any}()")
        @test endswith(sprint(show, OrderedSet(['a'])), "OrderedSet{Char}(['a'])")
    end

    @testset "Core Functionality" begin
        s = OrderedSet(); push!(s,1); push!(s,2); push!(s,3)
        @test !isempty(s)
        @test in(1,s)
        @test in(2,s)
        @test length(s) == 3
        push!(s,1); push!(s,2); push!(s,3)
        @test length(s) == 3
        @test pop!(s,1) == 1
        @test !in(1,s)
        @test in(2,s)
        @test length(s) == 2
        @test_throws KeyError pop!(s,1)
        @test pop!(s,1,:foo) == :foo
        @test length(delete!(s,2)) == 1
        @test !in(1,s)
        @test !in(2,s)
        @test pop!(s) == 3
        @test length(s) == 0
        @test isempty(s)
    end

    @testset "copy" begin
        data_in = (1,2,9,8,4)
        s = OrderedSet(data_in)
        c = copy(s)
        @test isequal(s,c)
        v = pop!(s)
        @test !in(v,s)
        @test  in(v,c)
        push!(s,100)
        push!(c,200)
        @test !in(100,c)
        @test !in(200,s)
    end

    @testset "sizehint!, empty" begin
        s = OrderedSet([1])
        @test isequal(sizehint!(s, 10), OrderedSet([1]))
        @test isequal(empty!(s), OrderedSet())
        # TODO: rehash
    end

    @testset "iterate" begin
        for data_in in ((7,8,4,5),
                        ("hello", 23, 2.7, (), [], (1,8)))
            s = OrderedSet(data_in)

            s_new = OrderedSet()
            for el in s
                push!(s_new, el)
            end
            @test isequal(s, s_new)

            t = tuple(s...)

            @test t === data_in
            @test length(t) == length(s)
            for (e,f) in zip(t,s)
                @test e === f
            end
        end
    end

    @testset "union" begin
        @test isequal(union(OrderedSet([1])),OrderedSet([1]))
        s = ∪(OrderedSet([1,2]), OrderedSet([3,4]))
        @test isequal(s, OrderedSet([1,2,3,4]))
        s = union(OrderedSet([5,6,7,8]), OrderedSet([7,8,9]))
        @test isequal(s, OrderedSet([5,6,7,8,9]))
        s = OrderedSet([1,3,5,7])
        union!(s,(2,3,4,5))
        # TODO: order is not the same, so isequal should return false...
        @test isequal(s,OrderedSet([1,2,3,4,5,7]))
    end

    @testset "intersect" begin
        @test isequal(intersect(OrderedSet([1])),OrderedSet([1]))
        s = ∩(OrderedSet([1,2]), OrderedSet([3,4]))
        @test isequal(s, OrderedSet())
        s = intersect(OrderedSet([5,6,7,8]), OrderedSet([7,8,9]))
        @test isequal(s, OrderedSet([7,8]))
        @test isequal(intersect(OrderedSet([2,3,1]), OrderedSet([4,2,3]), OrderedSet([5,4,3,2])), OrderedSet([2,3]))
    end

    @testset "setdiff" begin
        @test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet()),        OrderedSet([1,2,3]))
        @test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([1])),     OrderedSet([2,3]))
        @test isequal(setdiff(OrderedSet([1,2,3]), Set([1])),            OrderedSet([2,3]))
        @test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([1,2])),   OrderedSet([3]))
        @test isequal(setdiff(OrderedSet([1,2,3]), Set([1,2])),          OrderedSet([3]))
        @test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([1,2,3])), OrderedSet())
        @test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([4])),     OrderedSet([1,2,3]))
        @test isequal(setdiff(OrderedSet([1,2,3]), OrderedSet([4,1])),   OrderedSet([2,3]))
        s = OrderedSet([1,3,5,7])
        setdiff!(s,(3,5))
        @test isequal(s,OrderedSet([1,7]))
        s = OrderedSet([1,2,3,4])
        setdiff!(s, OrderedSet([2,4,5,6]))
        @test isequal(s,OrderedSet([1,3]))
    end

    @testset "ordering" begin
        @test OrderedSet() < OrderedSet([1])
        @test OrderedSet([1]) < OrderedSet([1,2])
        @test !(OrderedSet([3]) < OrderedSet([1,2]))
        @test !(OrderedSet([3]) > OrderedSet([1,2]))
        @test OrderedSet([1,2,3]) > OrderedSet([1,2])
        @test !(OrderedSet([3]) <= OrderedSet([1,2]))
        @test !(OrderedSet([3]) >= OrderedSet([1,2]))
        @test OrderedSet([1]) <= OrderedSet([1,2])
        @test OrderedSet([1,2]) <= OrderedSet([1,2])
        @test OrderedSet([1,2]) >= OrderedSet([1,2])
        @test OrderedSet([1,2,3]) >= OrderedSet([1,2])
        @test !(OrderedSet([1,2,3]) >= OrderedSet([1,2,4]))
        @test !(OrderedSet([1,2,3]) <= OrderedSet([1,2,4]))
    end

    @testset "issubset, symdiff" begin
        for (l,r) in ((OrderedSet([1,2]),     OrderedSet([3,4])),
                    (OrderedSet([5,6,7,8]), OrderedSet([7,8,9])),
                    (OrderedSet([1,2]),     OrderedSet([3,4])),
                    (OrderedSet([5,6,7,8]), OrderedSet([7,8,9])),
                    (OrderedSet([1,2,3]),   OrderedSet()),
                    (OrderedSet([1,2,3]),   OrderedSet([1])),
                    (OrderedSet([1,2,3]),   OrderedSet([1,2])),
                    (OrderedSet([1,2,3]),   OrderedSet([1,2,3])),
                    (OrderedSet([1,2,3]),   OrderedSet([4])),
                    (OrderedSet([1,2,3]),   OrderedSet([4,1])))
            @test issubset(intersect(l,r), l)
            @test issubset(intersect(l,r), r)
            @test issubset(l, union(l,r))
            @test issubset(r, union(l,r))
            @test isequal(union(intersect(l,r),symdiff(l,r)), union(l,r))
        end
        @test ⊆(OrderedSet([1]), OrderedSet([1,2]))

        @test ⊊(OrderedSet([1]), OrderedSet([1,2]))
        @test !⊊(OrderedSet([1]), OrderedSet([1]))
        @test ⊈(OrderedSet([1]), OrderedSet([2]))

        @test symdiff(OrderedSet([1,2,3,4]), OrderedSet([2,4,5,6])) == OrderedSet([1,3,5,6])

        @test isequal(symdiff(OrderedSet([1,2,3,4]), OrderedSet([2,4,5,6])), OrderedSet([1,3,5,6]))

    end

    @testset "filter" begin
        s = OrderedSet([1,2,3,4])
        @test isequal(filter(isodd,s), OrderedSet([1,3]))
        filter!(isodd, s)
        @test isequal(s, OrderedSet([1,3]))
    end

    @testset "first" begin
        @test_throws ArgumentError first(OrderedSet())
        @test first(OrderedSet([2])) == 2
    end

    @testset "empty set" begin
        d = OrderedSet{Char}()
        @test length(d) == 0
        @test isempty(d)
        @test !('c' in d)
        push!(d, 'c')
        @test !isempty(d)
        empty!(d)
        @test isempty(d)
    end

    @testset "access, modification" begin
        d = OrderedSet{Char}()

        for c in 'a':'z'
            push!(d, c)
        end

        for c in 'a':'z'
            @test c in d
        end

        @test collect(d) == collect('a':'z')
    end

end # @testset OrderedSet


text:
#!/home/skwarkmj/sw/julia/julia

using ArgParse
using GaussDCA
using PlmDCA

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--gaps", "-g"
            help = "Fraction of gaps to use"
        "--di"
            help = "an option without argument, i.e. a flag"
            action = :store_true
        "arg1"
            help = "input alignment"
            required = true
        "arg2"
            help = "lambdaJ 0"
            required = true
        "arg3"
            help = "Output file"
            required = true
    end

    return parse_args(s)
end

function main()
   parsed_args = parse_commandline()
   x = PlmDCA.plmdca(parsed_args["arg1"], lambdaJ=float(parsed_args["arg2"]))
   score=x.score
   Jtensor=x.Jtensor
   hvector=x.hvector

   io = open(parsed_args["arg3"], "w")
   for s in score
        @printf(io, "%d,%d,%f\n",s[1],s[2],s[3])
   end
   close(io)

   io = open(string(parsed_args["arg3"], ".Jtensor"), "w")
   N = size(Jtensor, 3)
   q = size(Jtensor, 1)
   for i = 1:N, j = 1:N
	@printf(io, "%d %d", i, j)
	for q1=1:q, q2=1:q
		@printf(io, " %f", Jtensor[q1,q2,i,j])
	end
	@printf(io, "\n")
   end
   close(io)
   io = open(string(parsed_args["arg3"], ".hvector"), "w")
   N = size(hvector, 2)
   q = size(Jtensor, 1)
   for i = 1:N
	@printf(io, "%d", i)
	for q1=1:q
		@printf(io, " %f", hvector[q1,i])
	end
	@printf(io, "\n")
   end
   close(io)
end

main()


text:
"""
	open(f::FormattedFilename, args...; kwargs...) -> ::FormattedIO
	open(fn::Function, f::FormattedFilename, args...; kwargs...)

Similar to `open(filename::AbstractString, ...)`, but operates on the filename
wrapped by `f`.

Format/coding information associated with `f` is preserved in the returned
`FormattedIO`. To specify format/coding or infer from the IO stream, use
`specify` or `infer` on the returned object.

# Examples

```julia
julia> f1 = infer("insulin.dat")
FormattedFilename: insulin.dat
 unknown format

julia> f2 = open(f1)
FormattedIO: IOStream(<file insulin.dat>)
 unknown format

julia> f3 = infer(f2)
FormattedIO: IOStream(<file insulin.dat>)
 inferred format: structure/x-pdb
```
"""
function Base.open(f::FormattedFilename, args...; kwargs...)
	io = open(f.resource, args...; kwargs...)
	FormattedIO(io, f.format, f.coding, f.format_guesses, f.coding_guesses)
end

function Base.open(fn::Function, f::FormattedFilename, args...; kwargs...)
	f2 = open(f, args...; kwargs...)
	results = fn(f2)
	close(f2)
	results
end

"""
	close(f::FormattedIO)

Similar to `close(io::IO)`, closes the wrapped IO stream.
"""
Base.close(f::FormattedIO) = close(f.resource)

"""
	read(f::Formatted, args...; kwargs...) -> value

Similar to `read(filename|io, ...)`, but operates on the resource wrapped by
`f`, reading a single object.

The format and coding are resolved, and appropriate reader and decoder are
selected.

# Examples

```julia
julia> f1 = infer("insulin.pdb")
FormattedFilename: insulin.pdb
 inferred format: structure/x-pdb

julia> mol1 = read(f1)
...

julia> f2 = infer("myoglobin.gro")
FormattedFilename: myoglobin.gro
 inferred format: structure/x-gro

julia> f3 = open(f2)
FormattedIO: IOStream(<file myoglobin.gro>)
 inferred format: structure/x-gro

julia> mol2 = read(f3)
```
"""
function Base.read(f::Formatted, args...; kwargs...)
	format = resolveformat(f)
    coding = resolvecoding(f)
	read(f.resource, format, coding, args...; kwargs...)
end

function Base.read(filename::AbstractString, format::MIME,
		coding::Union{MIME,Nothing}, args...; kwargs...)
	open(filename) do io
		read(io, format, coding, args...; kwargs...)
	end
end

function Base.read(io::IO, format::MIME, coding::Union{MIME,Nothing}, args...;
		kwargs...)
	reader = resolvereader(format)
	if coding == nothing
		read(reader, format, io, args...; kwargs...)
	else
		decoder = resolvedecoder(coding)
		io2 = TranscodingStream(decoder(), io)
		read(reader, format, io2, args...; kwargs...)
	end
end

"""
	read(reader, mime, io, args...; kwargs...) -> value

Read `value` from `io` using a specific `reader` for the `mime` format.

This version should be specialised for readers that wish to handle a specific
format. Arguments `args` and `kwargs` make it possible to customize reader
behavior.

# Arguments

	reader::FormatHandler
	mime::MIME
	io::IO

# Example

```julia
struct MyIO <: Formats.FormatHandler end

function Base.read(::MyIO, ::MIME"image/png", io::IO)
	# Read data from io to construct value x...
	x
end
```
"""
Base.read

"""
	read!(f::Formatted, output, args...; kwargs...) -> value

Similar to `read!(filename|io, array, ...)`, but operates on the resource
wrapped by `f`, reading a single object into pre-allocated `output`.

The format and coding are resolved, and appropriate reader and decoder are
selected.

# Examples

```julia
julia> f1 = infer("insulin.pdb")
FormattedFilename: insulin.pdb
 inferred format: structure/x-pdb

julia> mol1 = MolecularModel()
0-particle 0-key MolecularModel

julia> read!(f1, mol1)
[...]

julia> f2 = infer("myoglobin.gro")
FormattedFilename: myoglobin.gro
 inferred format: structure/x-gro

julia> f3 = open(f2)
FormattedIO: IOStream(<file myoglobin.gro>)
 inferred format: structure/x-gro

julia> mol2 = MolecularModel()
0-particle 0-key MolecularModel

julia> read!(f3, mol2)
```
"""
function Base.read!(f::FormattedFilename, output, args...; kwargs...)
	open(f) do f2
		read!(f2, output, args...; kwargs...)
	end
end

function Base.read!(f::FormattedIO, output, args...; kwargs...)
	format = resolveformat(f)
	reader = resolvereader(format)
    coding = resolvecoding(f)
	if coding == nothing
		io = f.resource
	else
		decoder = resolvedecoder(coding)
		io = TranscodingStream(decoder(), f.resource)
	end
	read!(reader, format, io, output, args...; kwargs...)
end

"""
	read!(reader, mime, io, output, args...; kwargs...) -> output

Read `value` from `io` into pre-allocated `output` using a specific `reader`
for the `mime` format.

This version should be specialised for readers that wish to handle a specific
format and reading into pre-allocated output. Arguments `args` and `kwargs`
make it possible to customize reader behavior.
"""
Base.read!

"""
	write(f::Formatted, x, args...; kwargs...)

Similar to `write(filename|io, x, ...)`, but operates on the resource wrapped by
`f`.

The format and coding are resolved, and appropriate writer and encoder are
selected.

# Examples

```julia
julia> mol = read(infer("insulin.pdb"))
[...]

julia> write(infer("insulin.gro"), mol)
[...]

julia> f = open(infer("insulin.gro.gz"), "w")
FormattedIO: IOStream(<file insulin.gro.gz>)
 inferred format: structure/x-gro
 inferred coding: application/gzip

julia> write(f, mol)
[...]
```
"""
function Base.write(f::Formatted, x, args...; kwargs...)
	format = resolveformat(f)
	coding = resolvecoding(f)
	write(f.resource, format, coding, x, args...; kwargs...)
end

function Base.write(filename::AbstractString, format::MIME,
		coding::Union{MIME,Nothing}, x, args...; kwargs...)
	open(filename, "w") do io
		write(io, format, coding, x, args...; kwargs...)
	end
end

function Base.write(io::IO, format::MIME, coding::Union{MIME,Nothing}, x,
		args...; kwargs...)
	writer = resolvewriter(format)
	pos0 = position(io)
	if coding == nothing
		write(writer, format, io, x, args...; kwargs...)
	else
		encoder = resolveencoder(coding)
		io2 = TranscodingStream(encoder(), io)
		write(writer, format, io2, x, args...; kwargs...)
	end
	position(io) - pos0
end

"""
	write(writer, mime, io, x, args...; kwargs...)

Write `x` to `io` using a specific `writer` for the `mime` format.

This version should be specialised for writers that wish to handle a specific
format. Arguments `args` and `kwargs` make it possible to customize writer
behavior.

# Arguments

	writer::FormatHandler
	mime::Type{<:MIME}
	io::IO

# Example

```julia
struct MyIO <: Formats.FormatHandler end

function Base.write(::MyIO, ::MIME"image/png", io::IO, x::MyT)
	# Write x to the io stream...
end
```
"""
Base.write

"""
	openf(filename|formatted, args...; kwargs...) -> ::FormattedIO
    openf(fn::Function, filename|formatted, args...; kwargs...)

Similar to `open(filename, ...)`, but format/coding are inferred automatically
if missing.

# Arguments

	filename::AbstractString
    formatted::Formatted

This is a convenience function; `openf(filename)` is equivalent to
`open(formatted(filename))`.
"""
openf(filename::Union{AbstractString,FormattedFilename}, args...; kwargs...) =
        open(formatted(filename), args...; kwargs...)

function openf(fn::Function, filename::Union{AbstractString,FormattedFilename},
        args...; kwargs...)
    f = open(filename, args...; kwargs...)
    results = fn(f)
    close(f)
	results
end

"""
	readf(filename|io|formatted, args...; kwargs...) -> value

Similar to `read(filename|io, ...)`, but format/coding are inferred
automatically if missing.

# Arguments

    filename::AbstractString
    io::IO
    formatted::Formatted

This is a convenience function; `readf(resource)` is equivalent to
`read(formatted(resource))`.
"""
readf(resource::Union{AbstractString,IO,Formatted}, args...; kwargs...) =
        read(formatted(resource), args...; kwargs...)

"""
	readf!(filename|io|formatted, output, args...; kwargs...) -> output

Similar to `read!(filename|io, output, ...)`, but format/coding are inferred
automatically if missing.

# Arguments

    filename::AbstractString
    io::IO
    formatted::Formatted

This is a convenience function; `readf!(resource, output)` is equivalent to
`read!(formatted(resource), output)`.
"""
readf!(resource::Union{AbstractString,IO,Formatted}, output, args...;
		kwargs...) = read!(formatted(resource), output, args...; kwargs...)

"""
	writef(filename|io|formatted, x, args...; kwargs...)

Similar to `write(filename|io, x, ...)`, but format/coding are inferred
automatically if unspecified.

# Arguments

    filename::AbstractString
    io::IO
    formatted::Formatted

This is a convenience function; `writef(resource, x)` is equivalent to
`write(formatted(resource), x)`.
"""
writef(resource::Union{AbstractString,IO,Formatted}, x, args...; kwargs...) =
        write(formatted(resource), x, args...; kwargs...)


text:
#=
QuantumIsingModel
initialize!
getMeasureOperator
getStates

=#

"""
	QuantumIsingModel
store information for quantum ising model

```julia
externalfield = 1.0e-13
temperature = 0.0
spininfo = SpinInfo("quantum_ising_2","asym",externalfield,temperature)
trotterparameter = 0.01
trotteriteration = 100
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
testModel = QuantumIsingModel(spininfo, trotterinfo)
```
"""
type QuantumIsingModel{T} <:QuantumSpinModel{T}
	spininfo::SpinInfo{T}
	trotterinfo::TrotterInfo{T}

	factorW::Array{T,2}
	factorWp::Array{T,2}

	function QuantumIsingModel{T}(spininfo::SpinInfo{T}, trotterinfo::TrotterInfo{T})
		this = new{T}()
		this.spininfo = spininfo
		this.trotterinfo = trotterinfo
		initialize!(this)
	return this
	end
end

QuantumIsingModel{T}(spininfo::SpinInfo{T}, trotterinfo::TrotterInfo{T}) =
QuantumIsingModel{T}(spininfo, trotterinfo)

#---
# functions for constructor
function initialize!{T}(qising::QuantumIsingModel{T})
	checkInfo!(qising.spininfo, qising.trotterinfo)
	factorW = getFactorW(qising)
	factorWp = getFactorWp(qising)
	qising.factorW = factorW
	qising.factorWp = factorWp
end

function getFactorW{T}(qising::QuantumIsingModel{T})
	tau = getTrotterparameter(qising)
	coshterm = sqrt(cosh(tau))
	sinhterm = sqrt(sinh(tau))
	return factorW = [coshterm sinhterm; coshterm -sinhterm]
end

function getFactorWp{T}(qising::QuantumIsingModel{T})
	tau = getTrotterparameter(qising)
	h = getExternalfield(qising)
	expterm = exp(tau * h / 2.0)
	factorWp = [expterm one(T)/expterm; expterm -one(T)/expterm]./sqrt(2.0)
	return factorWp
end

#---
# functions for IsingModel
function getMeasureOperator{T}(qising::QuantumIsingModel{T})
	return [-one(T), one(T)]
end

function getStates(qising::QuantumIsingModel)
	return 2
end

function isSymmetricFactorization(qising::QuantumIsingModel)
	return false
end

# set functions
function setTemperature!{T}(qising::QuantumIsingModel{T}, temperature::T)
	if !isZeroTemperature(temperature)
		iteration = -log(2.0, getTrotterparameter(qising) * temperature)
		iteration = try
			Int(iteration)
			setTrotterIteration!(qising.trotterinfo, iteration)
		catch
			error("cannot set the temperature")
		end
	end
	setTemperature!(qising.spininfo, temperature)
	initialize!(qising)
end

function setEnvParameters!{T}(qising::QuantumIsingModel{T}, temperature::T, externalfield::T)
	setExternalfield!(qising, externalfield)
	setTemperature!(qising,temperature)
end

#---
# test to make a example SpinModel
function testQuantumIsingModel()
	externalfield = 1.0e-13
	temperature = 0.0
	spininfo = SpinInfo("quantum_ising_2","asym",externalfield,temperature)
	trotterparameter = 0.01
	trotteriteration = 100
	trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
	testModel = QuantumIsingModel(spininfo, trotterinfo)
	return testModel
end


text:
@testset "Treebanks" begin

    datadir = joinpath(@__DIR__, "data")

    @testset "Treebank Readers" begin
        @testset "Typed" begin
            T = TypedDependency
            tb1 = TreebankReader{T}(joinpath(datadir, "wsj_0001.dp"), add_id=true)
            tb2 = TreebankReader{T}(joinpath(datadir, "wsj_0001.ids.dp"), add_id=false)
            @test deptype(tb1) == deptype(tb2) == T
            @test collect(tb1) == collect(tb2)
        end

        @testset "Untyped" begin
            T = UntypedDependency
            tb1 = TreebankReader{T}(joinpath(datadir, "wsj_0001.untyped.dp"), add_id=true)
            tb2 = TreebankReader{T}(joinpath(datadir, "wsj_0001.untyped.ids.dp"), add_id=false)
            @test deptype(tb1) == deptype(tb2) == T
            @test collect(tb1) == collect(tb2)
        end

        @testset "CoNLLU Multiword Tokens" begin
            tree = collect(TreebankReader{CoNLLU}(joinpath(datadir, "multiwordtoken.conllu")))[1]
            @test length(tree) == 24
            @test length(tree.mwts) == 1
            @test length(tree.emptytokens) == 0
        end

        @testset "CoNLLU Empty Tokens" begin
            tree = collect(TreebankReader{CoNLLU}(joinpath(datadir, "emptytokens.conllu")))[1]
            @test length(tree) == 6
            @test length(tree.mwts) == 0
            @test length(tree.emptytokens) == 1
        end
    end

    @testset "Treebanks" begin
        files = [joinpath(datadir, file) for file in readdir(datadir)
                 if endswith(file, ".conllu")]
        @test length(files) == 3

        @test_throws Exception Treebank{CoNLLU}("not a corpus")

        treebank = Treebank{CoNLLU}(files)
        @test showstr(treebank) == "Treebank{CoNLLU} of 3 file(s)"
        trees = collect(treebank)
        @test length(treebank) == length(trees) == 4
        @test deptype(treebank) == CoNLLU

        treebank2 = Treebank{CoNLLU}(datadir, pattern=r".conllu$")
        trees2 = collect(treebank2)
        @test length(trees2) == 4
        @test trees2 == trees

        for TS in (ArcEager, ArcHybrid, ArcStandard, ArcSwift)
            oracle = StaticOracle(TS())
            @test collect(xys(oracle, trees)) == collect(xys(oracle, trees2))
        end

        np = joinpath(datadir, "nonprojective.conll")
        ptreebank = Treebank{CoNLLU}(np, allow_nonprojective=false)
        nptreebank = Treebank{CoNLLU}(np, allow_nonprojective=true)
        @test length(collect(ptreebank)) == 2
        @test length(collect(nptreebank)) == 3
    end

    @testset "Oracles & Projectivity" begin
        tb = Treebank{CoNLLU}(joinpath(datadir, "nonprojective1.conll"))
        np_oracle = StaticOracle(ListBasedNonProjective())
        p_oracle  = StaticOracle(ArcEager())
        np_data = collect(xys(np_oracle, tb))
        p_data = collect(xys(p_oracle, tb))
        @test length(np_data) > 1
        @test length(p_data) == 0
        @test length(xys(p_oracle, first(tb))) == 0

        count = 0
        for (cfg, t) in xys(p_oracle, first(tb))
            count += 1
        end
        @test count == 0
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -7.107954588851326 <= q <= 8.439913276188518)

                     add_NL_constraint(m, :(gelu(-0.4953615570310643 + -0.5197659840811282*gelu(0.7941458217392698 + 0.49653685598476827*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.07122185861106178*gelu(-0.06420985369672083 + 0.9131512287471262*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.004105284291951339*gelu(-0.8124442669202567 + 0.4437121866105662*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) + gelu(0.7344840166166842 + -0.5219551557665656*gelu(0.7941458217392698 + 0.49653685598476827*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.12237894168399821*gelu(-0.06420985369672083 + 0.9131512287471262*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.5945370485807278*gelu(-0.8124442669202567 + 0.4437121866105662*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) + gelu(-0.7546740103129919 + -0.5743271753296177*gelu(0.7941458217392698 + 0.49653685598476827*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.8533482485058403*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + -0.13525586788633204*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + 0.959961363428667*gelu(-0.06420985369672083 + 0.9131512287471262*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.03263883937508805*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.0015134736925910275*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))) + -0.0685664132974031*gelu(-0.8124442669202567 + 0.4437121866105662*gelu(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3])) + 0.26832550078458484*gelu(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3])) + 0.2804882536336004*gelu(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# this script invokes the solver for any physics

using PDESolver
using EulerEquationMod
using AdvectionEquationMod
using SimpleODEMod
using EllipticEquationMod

mesh, sbp, eqn, opts = run_solver(ARGS[1])


text:
# convenience copnstructors for constant / linear / cubic spline interpolations
# 1D version
ConstantInterpolation(range::AbstractRange, vs::AbstractVector; extrapolation_bc = Throw()) =
    extrapolate(scale(interpolate(vs, BSpline(Constant())), range), extrapolation_bc)
ConstantInterpolation(range::AbstractVector, vs::AbstractVector; extrapolation_bc = Throw()) =
    extrapolate(interpolate((range, ), vs, Gridded(Constant())), extrapolation_bc)
LinearInterpolation(range::AbstractRange, vs::AbstractVector; extrapolation_bc = Throw()) =
    extrapolate(scale(interpolate(vs, BSpline(Linear())), range), extrapolation_bc)
LinearInterpolation(range::AbstractVector, vs::AbstractVector; extrapolation_bc = Throw()) =
    extrapolate(interpolate((range, ), vs, Gridded(Linear())), extrapolation_bc)
CubicSplineInterpolation(range::AbstractRange, vs::AbstractVector;
                         bc = Line(OnGrid()), extrapolation_bc = Throw()) =
    extrapolate(scale(interpolate(vs, BSpline(Cubic(bc))), range), extrapolation_bc)

# multivariate versions
ConstantInterpolation(ranges::NTuple{N,AbstractRange}, vs::AbstractArray{T,N};
                    extrapolation_bc = Throw()) where {N,T} =
    extrapolate(scale(interpolate(vs, BSpline(Constant())), ranges...), extrapolation_bc)
ConstantInterpolation(ranges::NTuple{N,AbstractVector}, vs::AbstractArray{T,N};
                    extrapolation_bc = Throw()) where {N,T} =
    extrapolate(interpolate(ranges, vs, Gridded(Constant())), extrapolation_bc)
LinearInterpolation(ranges::NTuple{N,AbstractRange}, vs::AbstractArray{T,N};
                    extrapolation_bc = Throw()) where {N,T} =
    extrapolate(scale(interpolate(vs, BSpline(Linear())), ranges...), extrapolation_bc)
LinearInterpolation(ranges::NTuple{N,AbstractVector}, vs::AbstractArray{T,N};
                    extrapolation_bc = Throw()) where {N,T} =
    extrapolate(interpolate(ranges, vs, Gridded(Linear())), extrapolation_bc)
CubicSplineInterpolation(ranges::NTuple{N,AbstractRange}, vs::AbstractArray{T,N};
                         bc = Line(OnGrid()), extrapolation_bc = Throw()) where {N,T} =
    extrapolate(scale(interpolate(vs, BSpline(Cubic(bc))), ranges...), extrapolation_bc)

"""
    etp = LinearInterpolation(knots, A; extrapolation_bc=Throw())

A shorthand for `extrapolate(interpolate(knots, A, scheme), extrapolation_bc)`,
where `scheme` is either `BSpline(Linear())` or `Gridded(Linear())` depending on whether
`knots` are ranges or vectors.
"""
LinearInterpolation

"""
    etp = ConstantInterpolation(knots, A; extrapolation_bc=Throw())

A shorthand for `extrapolate(interpolate(knots, A, scheme), extrapolation_bc)`,
where `scheme` is either `BSpline(Constant())` or `Gridded(Constant())` depending on whether
`knots` are ranges or vectors.
"""
ConstantInterpolation

"""
    etp = CubicSplineInterpolation(knots, A; bc=Line(OnGrid()), extrapolation_bc=Throw())

A shorthand for `extrapolate(interpolate(knots, A, BSpline(Cubic(bc))), extrapolation_bc)`.
"""
CubicSplineInterpolation


text:
# ----------
# Vector Calculus Module, Only work on peroideric boundary!
# ----------

function Curl(B1::Array,B2::Array,B3::Array;
              Lx = 2π, Ly = Lx, Lz = Lx,T = Float32)
    nx,ny,nz = size(B1);
    grid = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
    cB1,cB2,cB3 = Curl(B1,B2,B3,grid;Lx = Lx, Ly = Ly, Lz = Lz,T = T)
    return cB1,cB2,cB3;
end

function Curl(B1::Array,B2::Array,B3::Array,grid;
              Lx = 2π, Ly = Lx, Lz = Lx,T = Float32)
    #funtion of computing ∇×Vector using the fourier method
    # fft(∇×Vector) -> im * k × V
    #| i j k  |
    #| x y z  |
    #|B1 B2 B3|
    #
    nx,ny,nz = size(B1);
    B1h = zeros(Complex{T},(div(nx,2)+1,ny,nz));
    B2h = zeros(Complex{T},(div(nx,2)+1,ny,nz));
    B3h = zeros(Complex{T},(div(nx,2)+1,ny,nz));
    cBxh = copy(B1h);
    cByh = copy(B2h);
    cBzh = copy(B3h);
    mul!(B1h, grid.rfftplan, B1);
    mul!(B2h, grid.rfftplan, B2);
    mul!(B3h, grid.rfftplan, B3);

    for k in 1:nz, j in 1:ny,i in 1:div(nx,2)+1
       x,y,z = grid.kr[i],grid.l[j],grid.m[k];
       cBxh[i,j,k] = im*(y*B3h[i,j,k] - z*B2h[i,j,k]);
       cByh[i,j,k] = im*(z*B1h[i,j,k] - x*B3h[i,j,k]);
       cBzh[i,j,k] = im*(x*B2h[i,j,k] - y*B1h[i,j,k]);
    end

    cB1,cB2,cB3 = zeros(T,size(B1)),zeros(T,size(B1)),zeros(T,size(B1));
    ldiv!(cB1, grid.rfftplan, deepcopy(cBxh));
    ldiv!(cB2, grid.rfftplan, deepcopy(cByh));
    ldiv!(cB3, grid.rfftplan, deepcopy(cBzh));
    return cB1,cB2,cB3
end

function Div(B1::Array,B2::Array,B3::Array;
             Lx = 2π, Ly = Lx, Lz = Lx,T = Float32)
    nx,ny,nz = size(B1);
    grid = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
    cB1 = Div(B1,B2,B3,grid;Lx = Lx, Ly = Ly, Lz = Lz,T = T);

    return cB1
end

function Div(B1::Array,B2::Array,B3::Array,grid;
             Lx = 2π, Ly = Lx, Lz = Lx,T = Float32)
    #funtion of computing ∇̇ ⋅ Vector using the fourier method
    # fft(∇·Vector) -> im * k ⋅ V
    # = im* x*B1 + y*B2 + z*B3
    nx,ny,nz = size(B1);
    B1h = zeros(Complex{T},(div(nx,2)+1,ny,nz));
    B2h = zeros(Complex{T},(div(nx,2)+1,ny,nz));
    B3h = zeros(Complex{T},(div(nx,2)+1,ny,nz));
    Dot = copy(B1h);
    mul!(B1h, grid.rfftplan, B1);
    mul!(B2h, grid.rfftplan, B2);
    mul!(B3h, grid.rfftplan, B3);

    for k in 1:nz, j in 1:ny,i in 1:div(nx,2)+1
       x,y,z = grid.kr[i],grid.l[j],grid.m[k];
       Dot[i,j,k] = x*B1h[i,j,k] + y*B2h[i,j,k] + z*B3h[i,j,k];
    end

    cB1 = zeros(T,size(B1))
    ldiv!(cB1, grid.rfftplan, deepcopy(Dot));

    return cB1
end

function LaplaceSolver(B::Array; Lx=2π, Ly = Lx, Lz = Lz, T = Float32)
    nx,ny,nz = size(B);
    grid = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
    Φ   = LaplaceSolver(B,grid; Lx=2π, Ly = Lx, Lz = Lz, T = Float32);
    return Φ
end

function LaplaceSolver(B::Array,grid; Lx=2π, Ly = Lx, Lz = Lz, T = Float32)
    #=
    funtion of computing ΔΦ = B using the fourier method, must be peroidic condition
    Considering in k-space, k² Φ' = B', we would get Φ = F(B'/k²)
    =#
    nx,ny,nz = size(B);
    Φ    = zeros(T,nx,ny,nz);
    Bh   = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    mul!(Bh, grid.rfftplan, B);
    for k in 1:nz, j in 1:ny, i in 1:div(nx,2)+1
       x,y,z = grid.kr[i],grid.l[j],grid.m[k];
       k² = x^2 + y^2 + z^2;
       Bh[i,j,k] = Bh[i,j,k]/k²;
       if k² == 0; Bh[i,j,k] = 0; end
    end
    ldiv!(Φ, grid.rfftplan, deepcopy(Bh));
    return Φ;
end

function Crossproduct(A1,A2,A3,B1,B2,B3)
    C1 = @.  (A2*B3 - A3*B2);
    C2 = @. -(A1*B3 - A3*B1);
    C3 = @.  (A1*B2 - A2*B1);
    return C1,C2,C3
end

function Dotproduct(A1,A2,A3,B1,B2,B3)
    return A1.*B1 + A2.*B2 + A3.*B3
end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
export PermSparse
export randpermsparse

immutable PermSparse{T<:Real} <: AbstractPerm{T}
    data::Dict{T}
    plen::T
end

PermSparse{T<:Real}(p::AbstractVector{T}) =  PermSparse(PermPlain.listtosparse(p)...)
PermSparse{T<:Real}(cycs::AbstractArray{Array{T,1},1}) = PermSparse(PermPlain.cycstosparse(cycs)...)

PermSparse(a::Tuple) = PermSparse(PermPlain.tupcollect(a))
PermSparse{T<:Real}(::Type{T}) = PermSparse(T[])
PermSparse{T<:Real}(::Type{T}, a::Tuple) = PermSparse(PermPlain.tupcollect(T,a))
PermSparse() = PermSparse(Int[])

function getindex{T}(ps::PermSparse{T}, k::Real)
    res = get(ps.data,k,zero(T))
    res == zero(T) ? convert(T,k) : res
end

pmap(ps::PermSparse, k::Real) = getindex(ps,k)
copy(p::PermSparse) = PermSparse(copy(p.data),p.plen)
setindex!(p::PermSparse, i::Int, k::Integer) = (p.data)[k] = i
length(ps::PermSparse) = ps.plen
plength(ps::PermSparse) = ps.plen
isperm(ps::PermSparse) = PermPlain.isperm(ps.data)
isid(ps::PermSparse) = length(ps.data) == 0
order(ps::PermSparse) = PermPlain.permorder(ps.data)
supportsize(ps::PermSparse) = length(ps.data)
support(ps::PermSparse) = collect(keys(ps.data))
sign(ps::PermSparse) = PermPlain.permsgn_from_lengths(PermPlain.cyclelengths(ps.data))
cyclelengths(ps::PermSparse) = PermPlain.cyclelengths(ps.data)

one{T}(::Type{PermSparse{T}}) = PermSparse(T)
one(::Type{PermSparse}) = PermSparse(Int)
one{T}(p::PermSparse{T}) = PermSparse(T)
one(p::PermSparse) = PermSparse(eltype(p))

*(ps::PermSparse, k::Real) = ps[k]
*{T<:String}(p::PermSparse, v::T) = PermPlain.permapply(p.data,v)
*{T<:AbstractVector}(p::PermSparse, v::T) = PermPlain.permapply(p.data,v)
/(k::Real, p::PermSparse) = PermPlain.preimage(p.data,k)
*(p1::PermSparse, p2::PermSparse) = PermSparse(PermPlain.permcompose(p1.data,p2.data)...)
^(ps::PermSparse, n::Integer) = PermSparse(PermPlain.permpower(ps.data,n)...)
^(ps::PermSparse, n::FloatingPoint) = PermSparse(PermPlain.permpower(ps.data,n)...)
==(p1::PermSparse, p2::PermSparse) = p1.data == p2.data
function kron(p1::PermSparse, p2::PermSparse)
    PermSparse(PermPlain.permkron(p1.data,p2.data)...)
end


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/elastodynamics/demo_elastodynamics.py.html

module elastodynamics

using FenicsPy

# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = true
parameters["form_compiler"]["optimize"] = true

# Define mesh
_mesh = BoxMesh(Point(0., 0., 0.), Point(1., 0.1, 0.04), 60, 10, 5)

# Sub domain for clamp at left end
left = "near(x[0], 0.) && on_boundary"

# Sub domain for rotation at right end
right = "near(x[0], 1.) && on_boundary"

# Elastic parameters
E  = 1000.0
nu = 0.3
mu    = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

# Mass density
rho = 1.0

# Rayleigh damping coefficients
eta_m = 0.
eta_k = 0.

# Generalized-alpha method parameters
alpha_m = 0.2
alpha_f = 0.4
gamma   = 0.5+alpha_f-alpha_m
beta    = (gamma+0.5)^2/4.

# Time-stepping parameters
T       = 4.0
Nsteps  = 50
dt = T/Nsteps

p0 = 1.
cutoff_Tc = T/5
# Define the loading as an expression depending on t
p = Expression(("0", "t <= tc ? p0*t/tc : 0", "0"), t=0, tc=cutoff_Tc, p0=p0, degree=0)

# Define function space for displacement, velocity and acceleration
V = VectorFunctionSpace(_mesh, "CG", 1)
# Define function space for stresses
Vsig = TensorFunctionSpace(_mesh, "DG", 0)

# Test and trial functions
du = TrialFunction(V)
u_ = TestFunction(V)
# Current (unknown) displacement
u = FeFunction(V, name="Displacement")
# Fields from previous time step (displacement, velocity, acceleration)
u_old = FeFunction(V)
v_old = FeFunction(V)
a_old = FeFunction(V)

# Create mesh function over the cell facets
boundary_subdomains = MeshFunction("size_t", _mesh, _mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
# force_boundary = AutoSubDomain(right)
force_boundary = CompiledSubDomain(right)
force_boundary.mark(boundary_subdomains, 3)

# Define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)

# Set up boundary condition at left end
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, left)

# Stress tensor
function sigma(r)
    return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r))
end

# Mass form
function m(u, u_)
    return rho*inner(u, u_)*dx
end

# Elastic stiffness form
function k(u, u_)
    return inner(sigma(u), sym(grad(u_)))*dx
end

# Rayleigh damping form
function c(u, u_)
    return eta_m*m(u, u_) + eta_k*k(u, u_)
end

# Work of external forces
function Wext(u_)
    return dot(u_, p)*dss(3)
end

# Update formula for acceleration
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
function update_a(u, u_old, v_old, a_old)
    return (u-u_old-dt*v_old)/beta/dt^2 - (1-2*beta)/2/beta*a_old
end

# Update formula for velocity
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
function update_v(a, u_old, v_old, a_old)
    return v_old + dt*((1-gamma)*a_old + gamma*a)
end

function update_fields(u, u_old, v_old, a_old)
    """Update fields at the end of each time step."""

    # Get vectors (references)
    u_vec, u0_vec  = u.vector(), u_old.vector()
    v0_vec, a0_vec = v_old.vector(), a_old.vector()

    # use update functions using vector arguments
    a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec)
    v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec)

    # Update (u_old <- u)
    v_old.vector()[:] = v_vec
    a_old.vector()[:] = a_vec
    u_old.assign(u)
end

function avg(x_old, x_new, alpha)
    return alpha*x_old + (1-alpha)*x_new
end

# Residual
a_new = update_a(du, u_old, v_old, a_old)
v_new = update_v(a_new, u_old, v_old, a_old)
res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_) +
      k(avg(u_old, du, alpha_f), u_) - Wext(u_)
a_form = lhs(res)
L_form = rhs(res)

# Define solver for reusing factorization
K, res = assemble_system(a_form, L_form, bc)
solver = LUSolver(K, "mumps")
solver.parameters["symmetric"] = true

time = [i*T/Nsteps for i in 0:Nsteps]
u_tip = zeros(Nsteps+1)
energies = zeros(Nsteps+1, 4)
E_damp = 0
E_ext = 0
sig = FeFunction(Vsig, name="sigma")
xdmf_file = XDMFFile("elastodynamics/elastodynamics-results.xdmf")
xdmf_file.parameters["flush_output"] = true
xdmf_file.parameters["functions_share_mesh"] = true
xdmf_file.parameters["rewrite_function_mesh"] = false

function local_project(v, V, u=nothing)
    """Element-wise projection using LocalSolver"""
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u === nothing
        u = FeFunction(V)
        solver.solve_local_rhs(u)
        return u
    else
        solver.solve_local_rhs(u)
        return u
    end
end

for i in 1:Nsteps
    dt = time[i+1]-time[i]
    t = time[i+1]

    # Forces are evaluated at t_{n+1-alpha_f}=t_{n+1}-alpha_f*dt
    p.t = t-float(alpha_f*dt)

    # Solve for new displacement
    res = assemble(L_form)
    bc.apply(res)
    solver.solve(K, u.vector(), res)

    # Update old fields with new quantities
    update_fields(u, u_old, v_old, a_old)

    # Save solution to XDMF format
    xdmf_file.write(u, t)

    # Compute stresses and save to file
    local_project(sigma(u), Vsig, sig)
    xdmf_file.write(sig, t)

    p.t = t
    # Record tip displacement and compute energies
    # Note: Only works in serial
    if dolfin.MPI.comm_world.size == 1
        u_tip[i+1] = u(1., 0.05, 0.)[2]
    end
    E_elas = assemble(0.5*k(u_old, u_old))
    E_kin = assemble(0.5*m(v_old, v_old))
    global E_damp += dt*assemble(c(v_old, v_old))
    # E_ext += assemble(Wext(u-u_old))
    E_tot = E_elas+E_kin+E_damp #-E_ext
    energies[i+1, :] = [E_elas, E_kin, E_damp, E_tot]
    println("Time: ", t, "\tu_tip: ", u_tip[i+1], "\tenergies : ", energies[i+1, :] )
end

# 确保文件不被损坏
xdmf_file.close()

#using PyPlot

#if dolfin.MPI.comm_world.size == 1
#    plot(time, u_tip)
#    xlabel("Time")
#    ylabel("Tip displacement")
#    ylim(-0.5, 0.5)
#end

#if dolfin.MPI.comm_world.rank == 0
#    plot(time, energies)
#    legend(("elastic", "kinetic", "damping", "total"))
#    xlabel("Time")
#    ylabel("Energies")
#    ylim(0, 0.0011)
#end

end # module elastodynamics