text:
module PlasmoPipsNlpInterface

using MathProgBase.SolverInterface
using SparseArrays
using LinearAlgebra
import MPI
import JuMP
import Plasmo.PlasmoModelGraph
#import ..PlasmoModelGraph

include("PipsNlpSolver.jl")
using .PipsNlpSolver
export pipsnlp_solve

mutable struct ModelData
    d    #NLP evaluator
    n::Int
    m::Int
    local_m::Int
    jacnnz::Int
    hessnnz::Int
    firstIeq::Vector{Int}           #row index of equality constraint in 1st stage
    firstJeq::Vector{Int}           #column index of equality constraint in 1st stage
    firstVeq::Vector{Float64}       #coefficient of variable in equality constraint in 1st stage
    secondIeq::Vector{Int}          #row index of equality constraint in 2nd stage
    secondJeq::Vector{Int}          #column index of equality constraint in 2nd stage
    secondVeq::Vector{Float64}
    firstIineq::Vector{Int}
    firstJineq::Vector{Int}
    firstVineq::Vector{Float64}
    secondIineq::Vector{Int}
    secondJineq::Vector{Int}
    secondVineq::Vector{Float64}
    num_eqconnect::Int
    num_ineqconnect::Int
    eqconnect_lb::Vector{Float64}
    eqconnect_ub::Vector{Float64}
    ineqconnect_lb::Vector{Float64}
    ineqconnect_ub::Vector{Float64}
    eq_idx::Vector{Int}
    ineq_idx::Vector{Int}
    firstJeqmat
    secondJeqmat
    firstJineqmat
    secondJineqmat
    linkIeq::Vector{Int}
    linkJeq::Vector{Int}
    linkVeq::Vector{Float64}
    linkIineq::Vector{Int}
    linkJineq::Vector{Int}
    linkVineq::Vector{Float64}
    x_sol::Vector{Float64}
    coreid::Int
    loaded::Bool
    local_unsym_hessnnz::Int
end
ModelData() = ModelData(nothing,0,0,0,0,0,Int[],Int[], Float64[], Int[], Int[], Float64[],Int[],Int[],Float64[], Int[], Int[], Float64[], 0, 0, Float64[],Float64[],Float64[],Float64[], Int[], Int[],nothing, nothing, nothing, nothing, Int[],Int[], Float64[], Int[], Int[], Float64[],Float64[], 0, false, 0)

#Helper function
function getData(m::JuMP.Model)
    if haskey(m.ext, :Data)
        return m.ext[:Data]
    else
        error("This functionality is only available to model with extension data")
    end
end

function pipsnlp_solve(graph::PlasmoModelGraph.AbstractModelGraph,master_index::Int64,children_indices::Vector{Int64})

    if master_index == 0
        master_node = PlasmoModelGraph.ModelNode()
    else
        master_node = PlasmoModelGraph.getnode(graph,master_index)
    end
    children_nodes = [PlasmoModelGraph.getnode(graph,c_index) for c_index in children_indices]

    #need to check that the structure makes sense
    submodels = [PlasmoModelGraph.getmodel(child) for child in children_nodes]
    scen = length(children_nodes)
    #############################

    master = PlasmoModelGraph.getmodel(master_node)
    modelList = [master; submodels]

    #Add ModelData to each model
    for (idx,node) in enumerate(modelList)
        node.ext[:Data] = ModelData()
    end

    master_linear_lb = []
    master_linear_ub = []

    linkconstraints = PlasmoModelGraph.getlinkconstraints(graph) #get all of the link constraints in the graph
    #get arrays of lower and upper bounds for each link constraint
    for con in linkconstraints
        push!(master_linear_lb,con.lb)
        push!(master_linear_ub,con.ub)
    end

    nlinkeq = 0
    nlinkineq = 0
    eqlink_lb = Float64[]
    eqlink_ub = Float64[]
    ineqlink_lb = Float64[]
    ineqlink_ub = Float64[]

    #go through the link constraints
    for c in 1:length(linkconstraints)
        coeffs = linkconstraints[c].terms.coeffs
        vars   = linkconstraints[c].terms.vars
        connect = false     #the constraint connects a single node to the master
        allconnect = false  #the constraint connects subproblems to eachother
        node = nothing

		for (it,ind) in enumerate(coeffs)              #for each coefficient in the constraint
            #if the constraint has a variable from a sub-node
            if (!connect) && (vars[it].m) != master    #if the variable isn't in the master model
                connect = true
                node = vars[it].m                      #the submodel
		    end
            #if the constraint actually involved multiple nodes
		    if connect && (vars[it].m != master) && (vars[it].m != node)
		       allconnect = true
		       break  #don't need to check every variable if this statement runs
		    end
	    end

    	if (connect) && (!allconnect)   #if there's a master child connection only
            # println("Simple Links")
            # println(node.ext)
    	    local_data = getData(node)
    		if master_linear_lb[c] == master_linear_ub[c]  #if it's an equality constraint
                firstIeq = local_data.firstIeq
                firstJeq = local_data.firstJeq
                firstVeq = local_data.firstVeq
                secondIeq = local_data.secondIeq
                secondJeq = local_data.secondJeq
                secondVeq = local_data.secondVeq
                push!(local_data.eqconnect_lb, master_linear_lb[c])
                push!(local_data.eqconnect_ub, master_linear_ub[c])
                local_data.num_eqconnect += 1
                row = local_data.num_eqconnect
                for (it,ind) in enumerate(coeffs)
                    if (vars[it].m) == master
                           	  push!(firstIeq, row)
                          	  push!(firstJeq, vars[it].col)
                          	  push!(firstVeq, ind)
                    elseif (vars[it].m) == node
                           	  push!(secondIeq, row)
                          	  push!(secondJeq, vars[it].col)
                          	  push!(secondVeq, ind)
                    else
                      error("only supports connection between first stage variables and second stage variables from one specific scenario")
                    end
                end
    		else
                firstIineq = local_data.firstIineq
                firstJineq = local_data.firstJineq
                firstVineq = local_data.firstVineq
                secondIineq = local_data.secondIineq
                secondJineq = local_data.secondJineq
                secondVineq = local_data.secondVineq

                push!(local_data.ineqconnect_lb, master_linear_lb[c])
                push!(local_data.ineqconnect_ub, master_linear_ub[c])

                local_data.num_ineqconnect += 1
                row = local_data.num_ineqconnect
                for (it,ind) in enumerate(coeffs)
                    if (vars[it].m) == master
                        push!(firstIineq, row)
                        push!(firstJineq, vars[it].col)
                        push!(firstVineq, ind)
                    elseif (vars[it].m) == node
                        push!(secondIineq, row)
                        push!(secondJineq, vars[it].col)
                        push!(secondVineq, ind)
                    else
                        error("only supports connection between first stage variables and second stage variables from one specific scenario")
                    end
                end
	        end
	    end

	    if (allconnect)
            if master_linear_lb[c] == master_linear_ub[c]
                nlinkeq = nlinkeq + 1
                push!(eqlink_lb, master_linear_lb[c])
                push!(eqlink_ub, master_linear_ub[c])
                row = nlinkeq
                for (it,ind) in enumerate(coeffs)
                    node = vars[it].m
                    # println("Linking constraint")
                    # println(node.ext)
                    # println("Check model list again")
                    # println(modelList[1].ext)
                    local_data = getData(node)
                    linkIeq = local_data.linkIeq
                    linkJeq = local_data.linkJeq
                    linkVeq = local_data.linkVeq
                    push!(linkIeq, row)
                    push!(linkJeq, vars[it].col)
                    push!(linkVeq, ind)
                end
            else
                nlinkineq = nlinkineq + 1
                push!(ineqlink_lb, master_linear_lb[c])
                push!(ineqlink_ub, master_linear_ub[c])
                row = nlinkineq
                for (it,ind) in enumerate(coeffs)
                    node = vars[it].m
                    local_data = getData(node)
                    linkIineq = local_data.linkIineq
                    linkJineq = local_data.linkJineq
                    linkVineq = local_data.linkVineq
                    push!(linkIineq, row)
                    push!(linkJineq, vars[it].col)
                    push!(linkVineq, ind)
                end
            end
        end
    end
    #removeConnection(master, connectid)  #removes the linear constraints from the master model.  I don't add these to the master anyways
    master_data = getData(master)  #this isn't used anywhere?

    if haskey(master.ext, :nlinkeq)
        nlinkeq =  master.ext[:nlinkeq]
    end
    if haskey(master.ext, :nlinkineq)
        nlinkineq =  master.ext[:nlinkineq]
    end
    if haskey(master.ext, :eqlink_lb)
        eqlink_lb = copy(master.ext[:eqlink_lb])
        eqlink_ub = copy(master.ext[:eqlink_lb])
    end
    if haskey(master.ext, :eqlink_ub)
        eqlink_lb = copy(master.ext[:eqlink_ub])
        eqlink_ub = copy(master.ext[:eqlink_ub])
    end
    if haskey(master.ext, :ineqlink_lb)
        ineqlink_lb = copy(master.ext[:ineqlink_lb])
    end
    if haskey(master.ext, :ineqlink_ub)
        ineqlink_ub = copy(master.ext[:ineqlink_ub])
    end

    for (idx,node) in enumerate(modelList)
        local_data = getData(node)
        if haskey(node.ext, :linkIineq)
            local_data.linkIineq =  copy(node.ext[:linkIineq])
            node.ext[:linkIineq] = nothing
        end
        if haskey(node.ext, :linkJineq)
            local_data.linkJineq =  copy(node.ext[:linkJineq])
            node.ext[:linkJineq] = nothing
        end
        if haskey(node.ext, :linkVineq)
            local_data.linkVineq =  copy(node.ext[:linkVineq])
            node.ext[:linkVineq] = nothing
        end
        if haskey(node.ext, :linkIeq)
            local_data.linkIeq =  copy(node.ext[:linkIeq])
            node.ext[:linkIeq] = nothing
        end
        if haskey(node.ext, :linkJeq)
            local_data.linkJeq =  copy(node.ext[:linkJeq])
            node.ext[:linkJeq] = nothing
        end
        if haskey(node.ext, :linkVeq)
            local_data.linkVeq =  copy(node.ext[:linkVeq])
            node.ext[:linkVeq] = nothing
        end
    end

    function str_init_x0(nodeid, x0)
        node = modelList[nodeid+1]
        local_initval = copy(node.colVal)
        if any(isnan,node.colVal)
            local_initval[isnan.(node.colVal)] .= 0
            local_initval = min.(max.(node.colLower,local_initval),node.colUpper)
        end
        original_copy(local_initval,x0)
    end

    function str_prob_info(nodeid,flag, mode,col_lb,col_ub,row_lb,row_ub)
	 	#comm = MPI.COMM_WORLD
        if flag != 1
            node = modelList[nodeid+1]
            local_data = getData(node)
            if(!local_data.loaded)
                local_data.loaded = true
                if (nodeid > 0)
                    if haskey(node.ext, :warmStart)
    	                if node.ext[:warmStart] == true
              	      	   solve(node)
                   	   	end
                    end
            	end

    			nlp_lb, nlp_ub = JuMP.constraintbounds(node)
         		local_data.local_m  = length(nlp_lb)

    			newRowId = Array{Int}(undef,local_data.local_m)
    			eqId = 1
    			ineqId = 1
                for c in 1:local_data.local_m
                    if nlp_lb[c] == nlp_ub[c]
                        push!(local_data.eq_idx, c)
                        newRowId[c] = eqId
                        eqId +=  1
                    else
        				push!(local_data.ineq_idx, c)
        				newRowId[c] = ineqId
        				ineqId += 1
                    end
         		end
                local_data.m  = local_data.local_m + local_data.num_eqconnect + local_data.num_ineqconnect
                local_data.n = node.numCols
                local_data.x_sol = zeros(Float64,local_data.n)
                local_data.firstJeqmat = sparse(local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_data.num_eqconnect, master_data.n)
                local_data.secondJeqmat = sparse(local_data.secondIeq, local_data.secondJeq, local_data.secondVeq, local_data.num_eqconnect, local_data.n)
                local_data.firstJineqmat = sparse(local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_data.num_ineqconnect, master_data.n)
                local_data.secondJineqmat = sparse(local_data.secondIineq, local_data.secondJineq, local_data.secondVineq, local_data.num_ineqconnect, local_data.n)

                local_data.d = JuMP.NLPEvaluator(node)
                initialize(local_data.d, [:Grad,:Jac, :Hess])
                Ijac, Jjac = jac_structure(local_data.d)
                Ijaceq = Int[]
                Jjaceq = Int[]
                Ijacineq = Int[]
                Jjacineq = Int[]
                jac_eq_index = Int[]
                jac_ineq_index = Int[]
                for i in 1:length(Ijac)
                    c = Ijac[i]
                    if nlp_lb[c] == nlp_ub[c]
                        modifiedrow = newRowId[c]
                        push!(Ijaceq, modifiedrow)
                        push!(Jjaceq, Jjac[i])
                        push!(jac_eq_index, i)
                    else
                        modifiedrow = newRowId[c]
                        push!(Ijacineq, modifiedrow)
                        push!(Jjacineq, Jjac[i])
                        push!(jac_ineq_index,i)
                    end
                end

        		node.ext[:Ijaceq] = Ijaceq
         		node.ext[:Jjaceq] = Jjaceq
                node.ext[:Ijacineq] = Ijacineq
         		node.ext[:Jjacineq] = Jjacineq
         		node.ext[:jac_eq_index] = jac_eq_index
         		node.ext[:jac_ineq_index] = jac_ineq_index
         		Ihess, Jhess = hesslag_structure(local_data.d)

         		Hmap = Bool[]
         		node_Hrows = Int[]
         		node_Hcols = Int[]
         		for i in 1:length(Ihess)
            	    if Jhess[i] <= Ihess[i]
               	        push!(node_Hrows, Ihess[i])
                   		push!(node_Hcols, Jhess[i])
                   		push!(Hmap, true)
            	    else
                        push!(Hmap, false)
                    end
                end
         		val = ones(Float64,length(node_Hrows))
         		mat = sparse(node_Hrows,node_Hcols,val, local_data.n, local_data.n)
         		node.ext[:Hrows] = node_Hrows
         		node.ext[:Hcols] = node_Hcols
         		node.ext[:Hmap] = Hmap
         		local_hessnnz = length(mat.rowval)
         		local_data.local_unsym_hessnnz = length(Ihess)
         		local_data.hessnnz = local_hessnnz
		    end

	 	    if(mode == :Values)
               	nlp_lb, nlp_ub = JuMP.constraintbounds(node)
    			eq_lb=Float64[]
    			eq_ub=Float64[]
    			ineq_lb=Float64[]
    			ineq_ub=Float64[]
    			for i in 1: length(nlp_lb)
    			    if nlp_lb[i] == nlp_ub[i]
    			       push!(eq_lb, nlp_lb[i])
    			       push!(eq_ub, nlp_ub[i])
    			    else
    			       push!(ineq_lb, nlp_lb[i])
    			       push!(ineq_ub, nlp_ub[i])
    			    end
    			end

    			if nodeid !=  0
    			   eq_lb = [eq_lb; local_data.eqconnect_lb]
    			   eq_ub = [eq_ub; local_data.eqconnect_ub]
    			   ineq_lb = [ineq_lb; local_data.ineqconnect_lb]
    			   ineq_ub = [ineq_ub; local_data.ineqconnect_ub]
    			end

    			original_copy([eq_lb;ineq_lb], row_lb)
    			original_copy([eq_ub;ineq_ub], row_ub)
    			original_copy(node.colLower, col_lb)
    			original_copy(node.colUpper, col_ub)
            end
		    return (local_data.n,local_data.m)
		else
		    if(mode == :Values)
                original_copy([eqlink_lb; ineqlink_lb], row_lb)
                original_copy([eqlink_ub; ineqlink_ub], row_ub)
            end
            return (0, nlinkeq + nlinkineq)
		end
    end

    function str_eval_f(nodeid,x0,x1)
    	node = modelList[nodeid+1]
        local_data = getData(node)
        local_d = getData(node).d
        if nodeid ==  0
            local_x = x0
        else
            local_x = x1
        end
        local_scl = (node.objSense == :Min) ? 1.0 : -1.0
        f = local_scl*eval_f(local_d,local_x)
        return f
    end

    function str_eval_g(nodeid,x0,x1,new_eq_g, new_inq_g)
        node = modelList[nodeid+1]
        local_data = getData(node)
        local_d = getData(node).d
        if nodeid ==  0
            local_x = x0
        else
            local_x = x1
        end
        local_g = Array{Float64}(undef,local_data.local_m)
        eval_g(local_d, local_g, local_x)
        new_eq_g[1:end] = [local_g[local_data.eq_idx]; local_data.firstJeqmat*x0+local_data.secondJeqmat*x1]
        new_inq_g[1:end] = [local_g[local_data.ineq_idx]; local_data.firstJineqmat*x0+local_data.secondJineqmat*x1]
	    return Int32(1)
    end

    function str_write_solution(id::Integer, x::Vector{Float64}, y_eq::Vector{Float64}, y_ieq::Vector{Float64})
        node = modelList[id+1]
        local_data = getData(node)
        local_data.x_sol = copy(x)
        if id == 0
            node.objVal = str_eval_f(id, x, nothing)
        else
            node.objVal = str_eval_f(id, nothing, x)
        end
        r = MPI.Comm_rank(comm)
        local_data.coreid = r
    end

    function str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
    node = modelList[rowid+1]
    if rowid == colid
        local_data = getData(node)
        local_d = getData(node).d
        if colid ==  0
            local_x = x0
        else
            local_x = x1
        end
        local_grad_f = Array{Float64}(undef,local_data.n)
        eval_grad_f(local_d, local_grad_f, local_x)
        local_scl = (node.objSense == :Min) ? 1.0 : -1.0
        #scale!(local_grad_f,local_scl)
        rmul!(local_grad_f,local_scl)
        original_copy(local_grad_f, new_grad_f)
    elseif colid == 0
        new_grad_f[1:end] .= 0
    else
        @assert(false)
    end
    return Int32(1)
    end

    function array_copy(src,dest)
        @assert(length(src)==length(dest))
        for i in 1:length(src)
            dest[i] = src[i]-1
        end
    end

    function original_copy(src,dest)
        @assert(length(src)==length(dest))
        for i in 1:length(src)
            dest[i] = src[i]
        end
    end

    function str_eval_jac_g(rowid,colid,flag, x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
        if flag != 1
            node = modelList[rowid+1]
            local_data = getData(node)
            local_m_eq = length(local_data.eq_idx)
            local_m_ineq = length(local_data.ineq_idx)
            if mode == :Structure
                if (rowid == colid)
                    Ieq=[node.ext[:Ijaceq];local_data.secondIeq .+ local_m_eq]
                    Jeq=[node.ext[:Jjaceq];local_data.secondJeq]
                    Veq= ones(Float64, length(Ieq))
                    Iineq=[node.ext[:Ijacineq];local_data.secondIineq .+ local_m_ineq]
                    Jineq=[node.ext[:Jjacineq];local_data.secondJineq]
                    Vineq=ones(Float64, length(Iineq))
                    eqmat = sparse(Ieq, Jeq, Veq, local_m_eq + local_data.num_eqconnect, local_data.n)
                    ineqmat = sparse(Iineq, Jineq, Vineq, local_m_ineq + local_data.num_ineqconnect, local_data.n)
                else
                    eqmat = sparse(local_m_eq .+ local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_m_eq .+ local_data.num_eqconnect, master_data.n)
                    ineqmat = sparse(local_m_ineq .+ local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_m_ineq .+ local_data.num_ineqconnect, master_data.n)
                end
                return(length(eqmat.rowval), length(ineqmat.rowval))
            else
                if rowid == colid
                    if colid ==  0
                        local_x = x0
                    else
                        local_x = x1
                    end
                    local_values = Array{Float64}(undef,length(node.ext[:Ijaceq])+length(node.ext[:Ijacineq]))
                    eval_jac_g(local_data.d, local_values, local_x)
                    jac_eq_index = node.ext[:jac_eq_index]
                    jac_ineq_index = node.ext[:jac_ineq_index]
                    Ieq=[node.ext[:Ijaceq];local_data.secondIeq .+ local_m_eq]
                    Jeq=[node.ext[:Jjaceq];local_data.secondJeq]
                    Veq=[local_values[jac_eq_index];local_data.secondVeq]
                    Iineq=[node.ext[:Ijacineq];local_data.secondIineq .+ local_m_ineq]
                    Jineq=[node.ext[:Jjacineq];local_data.secondJineq]
                    Vineq=[local_values[jac_ineq_index];local_data.secondVineq]
                    eqmat = sparseKeepZero(Ieq, Jeq, Veq, local_m_eq + local_data.num_eqconnect, local_data.n)
                    ineqmat = sparseKeepZero(Iineq, Jineq, Vineq, local_m_ineq + local_data.num_ineqconnect, local_data.n)
                else
                    eqmat = sparseKeepZero(local_m_eq .+ local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_m_eq .+ local_data.num_eqconnect, master_data.n)
                    ineqmat = sparseKeepZero(local_m_ineq .+ local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_m_ineq .+ local_data.num_ineqconnect, master_data.n)
                end
                if length(eqmat.nzval) > 0
                    array_copy(eqmat.rowval,e_rowidx)
                    array_copy(eqmat.colptr,e_colptr)
                    original_copy(eqmat.nzval,e_values)
                end
                if length(ineqmat.nzval) > 0
                    array_copy(ineqmat.rowval,i_rowidx)
                    array_copy(ineqmat.colptr,i_colptr)
                    original_copy(ineqmat.nzval, i_values)
                end
            end
        else
            node = modelList[rowid+1]
            local_data = getData(node)
            linkIeq = local_data.linkIeq
            linkJeq = local_data.linkJeq
            linkVeq = local_data.linkVeq
            linkIineq = local_data.linkIineq
            linkJineq = local_data.linkJineq
            linkVineq = local_data.linkVineq

            if mode == :Structure
                return(length(linkVeq), length(linkVineq))
            else
                eqmat = sparse(linkIeq, linkJeq, linkVeq, nlinkeq, local_data.n)
                ineqmat = sparse(linkIineq, linkJineq, linkVineq, nlinkineq, local_data.n)
                if length(eqmat.nzval) > 0
                    array_copy(eqmat.rowval,e_rowidx)
                    array_copy(eqmat.colptr,e_colptr)
                    original_copy(eqmat.nzval,e_values)
                end
                if length(ineqmat.nzval) > 0
                    array_copy(ineqmat.rowval,i_rowidx)
                    array_copy(ineqmat.colptr,i_colptr)
                    original_copy(ineqmat.nzval, i_values)
                end
            end
        end
        return Int32(1)
    end

    function str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        node = modelList[colid+1]
        local_data = getData(node)
        if mode == :Structure
            if rowid == colid
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_val = ones(Float64,length(node_Hrows))
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                return length(mat.rowval)
            elseif colid == 0
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_Hrows,node_Hcols = exchange(node_Hrows,node_Hcols)
                node_val = ones(Float64,length(node_Hrows))
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                return length(mat.rowval)
            else
                return 0
            end
        else
            if rowid==colid
                node_Hmap = node.ext[:Hmap]
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                if colid ==  0
                    local_x = x0
                else
                    local_x = x1
                end
                local_unsym_values = Array{Float64}(undef,local_data.local_unsym_hessnnz)
                node_val = ones(Float64,length(node_Hrows))
                local_scl = (node.objSense == :Min) ? 1.0 : -1.0
                local_m_eq = length(local_data.eq_idx)
                local_m_ineq = length(local_data.ineq_idx)
                local_lambda = zeros(Float64, local_data.local_m)
                for i in 1:local_m_eq
                    local_lambda[local_data.eq_idx[i]] = lambda[i]
                end
                for i in 1:local_m_ineq
                    local_lambda[local_data.ineq_idx[i]] = lambda[i+local_m_eq]
                end

                eval_hesslag(local_data.d, local_unsym_values, local_x, obj_factor*local_scl, local_lambda)
                local_sym_index=1
                for i in 1:local_data.local_unsym_hessnnz
                    if node_Hmap[i]
                        node_val[local_sym_index] = local_unsym_values[i]
                        local_sym_index +=1
                    end
                end
                node_Hrows,node_Hcols = exchange(node_Hrows,node_Hcols)
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                array_copy(mat.rowval,rowidx)
                array_copy(mat.colptr,colptr)
                original_copy(mat.nzval,values)
            elseif colid ==	0
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_val = zeros(Float64,length(node_Hrows))
                node_Hrows,node_Hcols = exchange(node_Hrows,node_Hcols)
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                array_copy(mat.rowval,rowidx)
                array_copy(mat.colptr,colptr)
                original_copy(mat.nzval,values)
            else
            end
        end
        return Int32(1)
    end

    comm = MPI.COMM_WORLD
    if(MPI.Comm_rank(comm) == 0)
        t1 = time()
    end

    #Create FakeModel (The PIPS interface model) and pass all the functions it requires
    model = FakeModel(:Min,0, scen,
    str_init_x0, str_prob_info, str_eval_f, str_eval_g, str_eval_grad_f, str_eval_jac_g, str_eval_h,str_write_solution)
    prob = createProblemStruct(comm, model, true)
    ret = solveProblemStruct(prob)
    root = 0
    r = MPI.Comm_rank(comm)

    if(MPI.Comm_rank(comm) == 0)
        println("Timing Results:")
        println("init_x0  ",prob.t_jl_init_x0)
        println("str_init_x0  ", prob.t_jl_str_prob_info)
        println("eval_f  ", prob.t_jl_eval_f)
        println("eval_g0  ",prob.t_jl_eval_g)
        println("eval_grad_f  ", prob.t_jl_eval_grad_f)
        println("eval_jac  ", prob.t_jl_eval_jac_g)
        println("str_eval_jac  ", prob.t_jl_str_eval_jac_g)
        println("eval_h  ",  prob.t_jl_eval_h)
        println("str_eval_h ",  prob.t_jl_str_eval_h)
        println("eval_write_solution  ",  prob.t_jl_write_solution)
        println("PIPS-NLP time:   ",  time() - t1, " (s)")
    end

    for (idx,node) in enumerate(modelList)  #set solution values for each model
        local_data = getData(node)
        if idx != 1
            coreid = zeros(Int, 1)
            sc = MPI.Reduce(local_data.coreid, MPI.SUM, root, comm)
            if r == root
                coreid[1] = sc
            end
            MPI.Bcast!(coreid, length(coreid), root, comm)
            n = zeros(Int, 1)
            n[1] = local_data.n
            MPI.Bcast!(n,      length(n),      coreid[1], comm)
            if r != coreid[1]
                local_data.n = n[1]
                local_data.x_sol = zeros(Float64,local_data.n)
            end
            MPI.Bcast!(local_data.x_sol, local_data.n, coreid[1], comm)
            node.colVal = local_data.x_sol
        else
            node.colVal = local_data.x_sol
        end
    end

    status = :Unknown
    if ret == 0
        status = :Optimal
        #PlasmoModelGraph._setobjectivevalue(graph,prob.t_jl_eval_f)
    elseif ret == 1
        status = :Not_Finished
    elseif	ret == 2
        status = :Maximum_Iterations_Exceeded
    elseif	ret == 3
        stauts = :Infeasible_Problem_Detected
    elseif ret == 4
        status = :Restoration_needed
    else
        status = :UnknownStatus
    end

    return status
end  #end pips_nlp_solve

#############################################
# Helpers
#############################################
function exchange(a,b)
	 temp = a
         a=b
         b=temp
	 return (a,b)
end

function sparseKeepZero(I::AbstractVector{Ti},
    J::AbstractVector{Ti},
    V::AbstractVector{Tv},
    nrow::Integer, ncol::Integer) where {Tv,Ti<:Integer}
    N = length(I)
    if N != length(J) || N != length(V)
        throw(ArgumentError("triplet I,J,V vectors must be the same length"))
    end
    if N == 0
        return spzeros(eltype(V), Ti, nrow, ncol)
    end

    # Work array
    Wj = Array{Ti}(undef,max(nrow,ncol)+1)
    # Allocate sparse matrix data structure
    # Count entries in each row
    Rnz = zeros(Ti, nrow+1)
    Rnz[1] = 1
    nz = 0
    for k=1:N
        iind = I[k]
        iind > 0 || throw(ArgumentError("all I index values must be > 0"))
        iind <= nrow || throw(ArgumentError("all I index values must be ≤ the number of rows"))
        Rnz[iind+1] += 1
        nz += 1
    end
    Rp = cumsum(Rnz)
    Ri = Array{Ti}(undef,nz)
    Rx = Array{Tv}(undef,nz)

    # Construct row form
    # place triplet (i,j,x) in column i of R
    # Use work array for temporary row pointers
    @simd for i=1:nrow; @inbounds Wj[i] = Rp[i]; end
    @inbounds for k=1:N
        iind = I[k]
        jind = J[k]
        jind > 0 || throw(ArgumentError("all J index values must be > 0"))
        jind <= ncol || throw(ArgumentError("all J index values must be ≤ the number of columns"))
        p = Wj[iind]
        Vk = V[k]
        Wj[iind] += 1
        Rx[p] = Vk
        Ri[p] = jind
    end

    # Reset work array for use in counting duplicates
    @simd for j=1:ncol; @inbounds Wj[j] = 0; end

    # Sum up duplicates and squeeze
    anz = 0
    @inbounds for i=1:nrow
        p1 = Rp[i]
        p2 = Rp[i+1] - 1
        pdest = p1
        for p = p1:p2
            j = Ri[p]
            pj = Wj[j]
            if pj >= p1
                Rx[pj] = Rx[pj] + Rx[p]
            else
                Wj[j] = pdest
                if pdest != p
                    Ri[pdest] = j
                    Rx[pdest] = Rx[p]
                end
                pdest += one(Ti)
            end
        end
        Rnz[i] = pdest - p1
        anz += (pdest - p1)
    end

    # Transpose from row format to get the CSC format
    RiT = Array{Ti}(undef,anz)
    RxT = Array{Tv}(undef,anz)

    # Reset work array to build the final colptr
    Wj[1] = 1
    @simd for i=2:(ncol+1); @inbounds Wj[i] = 0; end
    @inbounds for j = 1:nrow
        p1 = Rp[j]
        p2 = p1 + Rnz[j] - 1
        for p = p1:p2
            Wj[Ri[p]+1] += 1
        end
    end
    RpT = cumsum(Wj[1:(ncol+1)])

    # Transpose
    @simd for i=1:length(RpT); @inbounds Wj[i] = RpT[i]; end
    @inbounds for j = 1:nrow
        p1 = Rp[j]
        p2 = p1 + Rnz[j] - 1
        for p = p1:p2
            ind = Ri[p]
            q = Wj[ind]
            Wj[ind] += 1
            RiT[q] = j
            RxT[q] = Rx[p]
        end
    end

    return SparseMatrixCSC(nrow, ncol, RpT, RiT, RxT)
end

function convert_to_c_idx(indicies)
    for i in 1:length(indicies)
        indicies[i] = indicies[i] - 1
    end
end

end #end module
# function removeConnection(master::JuMP.Model, deleterow)
#     deleteat!(master.linconstr,deleterow)
# end


text:
import Plots

"""Pre-execute hooks to set the plot size for the chunk """
function plots_set_size(chunk)
  w = chunk.options[:fig_width] * chunk.options[:dpi]
  h = chunk.options[:fig_height] * chunk.options[:dpi]
  Plots.default(size = (w,h))
  return chunk
end

push_preexecute_hook(plots_set_size)

ext: jl


text:
# common facilities

# scalar functions
"""
    xlogx(x::Number)

Compute `x * log(x)`, returning zero if `x` is zero.

```jldoctest
julia> StatsFuns.xlogx(0)
0.0
```
"""
function xlogx(x::Number)
    result = x * log(x)
    ifelse(iszero(x), zero(result), result)
end

"""
    xlogy(x::Number, y::Number)

Compute `x * log(y)`, returning zero if `x` is zero.

```jldoctest
julia> StatsFuns.xlogy(0, 0)
0.0
```
"""
function xlogy(x::Number, y::Number)
    result = x * log(y)
    ifelse(iszero(x) && !isnan(y), zero(result), result)
end

# The following bounds are precomputed versions of the following abstract
# function, but the implicit interface for AbstractFloat doesn't uniformly
# enforce that all floating point types implement nextfloat and prevfloat.
# @inline function _logistic_bounds(x::AbstractFloat)
#     (
#         logit(nextfloat(zero(float(x)))),
#         logit(prevfloat(one(float(x)))),
#     )
# end

@inline _logistic_bounds(x::Float16) = (Float16(-16.64), Float16(7.625))
@inline _logistic_bounds(x::Float32) = (-103.27893f0, 16.635532f0)
@inline _logistic_bounds(x::Float64) = (-744.4400719213812, 36.7368005696771)

"""
    logistic(x::Real)

The [logistic](https://en.wikipedia.org/wiki/Logistic_function) sigmoid function mapping a real number to a value in the interval [0,1],
```math
\\sigma(x) = \\frac{1}{e^{-x} + 1} = \\frac{e^x}{1+e^x}.
```

Its inverse is the [`logit`](@ref) function.
"""
logistic(x::Real) = inv(exp(-x) + one(x))

function logistic(x::Union{Float16, Float32, Float64})
    e = exp(x)
    lower, upper = _logistic_bounds(x)
    ifelse(
        x < lower,
        zero(x),
        ifelse(
            x > upper,
            one(x),
            e / (one(x) + e)
        )
    )
end

"""
    logit(p::Real)

The [logit](https://en.wikipedia.org/wiki/Logit) or log-odds transformation,
```math
\\log\\left(\\frac{x}{1-x}\\right), \\text{where} 0 < x < 1
```
Its inverse is the [`logistic`](@ref) function.
"""
logit(x::Real) = log(x / (one(x) - x))

"""
    log1psq(x::Real)

Return `log(1+x^2)` evaluated carefully for `abs(x)` very small or very large.
"""
log1psq(x::Real) = log1p(abs2(x))
function log1psq(x::Union{Float32,Float64})
    ax = abs(x)
    ax < maxintfloat(x) ? log1p(abs2(ax)) : 2 * log(ax)
end

"""
    log1pexp(x::Real)

Return `log(1+exp(x))` evaluated carefully for largish `x`.

This is also called the ["softplus"](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
transformation, being a smooth approximation to `max(0,x)`. Its inverse is [`logexpm1`](@ref).
"""
log1pexp(x::Real) = x < 18.0 ? log1p(exp(x)) : x < 33.3 ? x + exp(-x) : oftype(exp(-x), x)
log1pexp(x::Float32) = x < 9.0f0 ? log1p(exp(x)) : x < 16.0f0 ? x + exp(-x) : oftype(exp(-x), x)

"""
    log1mexp(x::Real)

Return `log(1 - exp(x))`

See:
 * Martin Maechler (2012) "Accurately Computing log(1 − exp(− |a|))",
   http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf

Note: different than Maechler (2012), no negation inside parentheses
"""
log1mexp(x::Real) = x < loghalf ? log1p(-exp(x)) : log(-expm1(x))

"""
    log2mexp(x::Real)

Return `log(2 - exp(x))` evaluated as `log1p(-expm1(x))`
"""
log2mexp(x::Real) = log1p(-expm1(x))

"""
    logexpm1(x::Real)

Return `log(exp(x) - 1)` or the "invsoftplus" function.
It is the inverse of [`log1pexp`](@ref) (aka "softplus").
"""
logexpm1(x::Real) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : oftype(exp(-x), x)
logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : oftype(exp(-x), x)

const softplus = log1pexp
const invsoftplus = logexpm1

"""
    log1pmx(x::Float64)

Return `log(1 + x) - x`

Use naive calculation or range reduction outside kernel range.  Accurate ~2ulps for all `x`.
"""
function log1pmx(x::Float64)
    if !(-0.7 < x < 0.9)
        return log1p(x) - x
    elseif x > 0.315
        u = (x-0.5)/1.5
        return _log1pmx_ker(u) - 9.45348918918356180e-2 - 0.5*u
    elseif x > -0.227
        return _log1pmx_ker(x)
    elseif x > -0.4
        u = (x+0.25)/0.75
        return _log1pmx_ker(u) - 3.76820724517809274e-2 + 0.25*u
    elseif x > -0.6
        u = (x+0.5)*2.0
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        u = (x+0.625)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    end
end

"""
    logmxp1(x::Float64)

Return `log(x) - x + 1` carefully evaluated.
"""
function logmxp1(x::Float64)
    if x <= 0.3
        return (log(x) + 1.0) - x
    elseif x <= 0.4
        u = (x-0.375)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    elseif x <= 0.6
        u = 2.0*(x-0.5)
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        return log1pmx(x - 1.0)
    end
end

# The kernel of log1pmx
# Accuracy within ~2ulps for -0.227 < x < 0.315
function _log1pmx_ker(x::Float64)
    r = x/(x+2.0)
    t = r*r
    w = @horner(t,
                6.66666666666666667e-1, # 2/3
                4.00000000000000000e-1, # 2/5
                2.85714285714285714e-1, # 2/7
                2.22222222222222222e-1, # 2/9
                1.81818181818181818e-1, # 2/11
                1.53846153846153846e-1, # 2/13
                1.33333333333333333e-1, # 2/15
                1.17647058823529412e-1) # 2/17
    hxsq = 0.5*x*x
    r*(hxsq+w*t)-hxsq
end

"""
    logaddexp(x::Real, y::Real)

Return `log(exp(x) + exp(y))`, avoiding intermediate overflow/undeflow, and handling non-finite values.
"""
function logaddexp(x::Real, y::Real)
    # ensure Δ = 0 if x = y = ± Inf
    Δ = ifelse(x == y, zero(x - y), abs(x - y))
    max(x, y) + log1pexp(-Δ)
end

Base.@deprecate logsumexp(x::Real, y::Real) logaddexp(x, y)

"""
    logsubexp(x, y)

Return `log(abs(e^x - e^y))`, preserving numerical accuracy.
"""
logsubexp(x::Real, y::Real) = max(x, y) + log1mexp(-abs(x - y))

"""
    logsumexp(X)

Compute `log(sum(exp, X))` in a numerically stable way that avoids intermediate over- and
underflow.

`X` should be an iterator of real numbers. The result is computed using a single pass over
the data.

# References

[Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)
"""
logsumexp(X) = _logsumexp_onepass(X)

"""
    logsumexp(X::AbstractArray{<:Real}; dims=:)

Compute `log.(sum(exp.(X); dims=dims))` in a numerically stable way that avoids
intermediate over- and underflow.

The result is computed using a single pass over the data.

# References

[Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)
"""
logsumexp(X::AbstractArray{<:Real}; dims=:) = _logsumexp(X, dims)

_logsumexp(X::AbstractArray{<:Real}, ::Colon) = _logsumexp_onepass(X)
function _logsumexp(X::AbstractArray{<:Real}, dims)
    # Do not use log(zero(eltype(X))) directly to avoid issues with ForwardDiff (#82)
    FT = float(eltype(X))
    xmax_r = reduce(_logsumexp_onepass_op, X; dims=dims, init=(FT(-Inf), zero(FT)))
    return @. first(xmax_r) + log1p(last(xmax_r))
end

function _logsumexp_onepass(X)
    # fallback for empty collections
    isempty(X) && return log(sum(X))
    return _logsumexp_onepass_result(_logsumexp_onepass_reduce(X, Base.IteratorEltype(X)))
end

# function barrier for reductions with single element and without initial element
_logsumexp_onepass_result(x) = float(x)
_logsumexp_onepass_result((xmax, r)::Tuple) = xmax + log1p(r)

# iterables with known element type
function _logsumexp_onepass_reduce(X, ::Base.HasEltype)
    # do not perform type computations if element type is abstract
    T = eltype(X)
    isconcretetype(T) || return _logsumexp_onepass_reduce(X, Base.EltypeUnknown())

    FT = float(T)
    return reduce(_logsumexp_onepass_op, X; init=(FT(-Inf), zero(FT)))
end

# iterables without known element type
_logsumexp_onepass_reduce(X, ::Base.EltypeUnknown) = reduce(_logsumexp_onepass_op, X)

## Reductions for one-pass algorithm: avoid expensive multiplications if numbers are reduced

# reduce two numbers
function _logsumexp_onepass_op(x1, x2)
    a = x1 == x2 ? zero(x1 - x2) : -abs(x1 - x2)
    xmax = x1 > x2 ? oftype(a, x1) : oftype(a, x2)
    r = exp(a)
    return xmax, r
end

# reduce a number and a partial sum
function _logsumexp_onepass_op(x, (xmax, r)::Tuple)
    a = x == xmax ? zero(x - xmax) : -abs(x - xmax)
    if x > xmax
        _xmax = oftype(a, x)
        _r = (r + one(r)) * exp(a)
    else
        _xmax = oftype(a, xmax)
        _r = r + exp(a)
    end
    return _xmax, _r
end
_logsumexp_onepass_op(xmax_r::Tuple, x) = _logsumexp_onepass_op(x, xmax_r)

# reduce two partial sums
function _logsumexp_onepass_op((xmax1, r1)::Tuple, (xmax2, r2)::Tuple)
    a = xmax1 == xmax2 ? zero(xmax1 - xmax2) : -abs(xmax1 - xmax2)
    if xmax1 > xmax2
        xmax = oftype(a, xmax1)
        r = r1 + (r2 + one(r2)) * exp(a)
    else
        xmax = oftype(a, xmax2)
        r = r2 + (r1 + one(r1)) * exp(a)
    end
    return xmax, r
end

"""
    softmax!(r::AbstractArray, x::AbstractArray)

Overwrite `r` with the `softmax` (or _normalized exponential_) transformation of `x`

That is, `r` is overwritten with `exp.(x)`, normalized to sum to 1.

See the [Wikipedia entry](https://en.wikipedia.org/wiki/Softmax_function)
"""
function softmax!(r::AbstractArray{R}, x::AbstractArray{T}) where {R<:AbstractFloat,T<:Real}
    n = length(x)
    length(r) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    u = maximum(x)
    s = 0.
    @inbounds for i = 1:n
        s += (r[i] = exp(x[i] - u))
    end
    invs = convert(R, inv(s))
    @inbounds for i = 1:n
        r[i] *= invs
    end
    r
end

"""
    softmax(x::AbstractArray{<:Real})

Return the [`softmax transformation`](https://en.wikipedia.org/wiki/Softmax_function) applied to `x`
"""
softmax!(x::AbstractArray{<:AbstractFloat}) = softmax!(x, x)
softmax(x::AbstractArray{<:Real}) = softmax!(similar(x, Float64), x)


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

mutable struct IoK8sApiCoreV1HTTPGetAction <: SwaggerModel
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    httpHeaders::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1HTTPHeader} } # spec name: httpHeaders
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    port::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: port
    scheme::Any # spec type: Union{ Nothing, String } # spec name: scheme

    function IoK8sApiCoreV1HTTPGetAction(;host=nothing, httpHeaders=nothing, path=nothing, port=nothing, scheme=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("httpHeaders"), httpHeaders)
        setfield!(o, Symbol("httpHeaders"), httpHeaders)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiCoreV1HTTPGetAction, Symbol("scheme"), scheme)
        setfield!(o, Symbol("scheme"), scheme)
        o
    end
end # type IoK8sApiCoreV1HTTPGetAction

const _property_map_IoK8sApiCoreV1HTTPGetAction = Dict{Symbol,Symbol}(Symbol("host")=>Symbol("host"), Symbol("httpHeaders")=>Symbol("httpHeaders"), Symbol("path")=>Symbol("path"), Symbol("port")=>Symbol("port"), Symbol("scheme")=>Symbol("scheme"))
const _property_types_IoK8sApiCoreV1HTTPGetAction = Dict{Symbol,String}(Symbol("host")=>"String", Symbol("httpHeaders")=>"Vector{IoK8sApiCoreV1HTTPHeader}", Symbol("path")=>"String", Symbol("port")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("scheme")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1HTTPGetAction }) = collect(keys(_property_map_IoK8sApiCoreV1HTTPGetAction))
Swagger.property_type(::Type{ IoK8sApiCoreV1HTTPGetAction }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1HTTPGetAction[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1HTTPGetAction }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1HTTPGetAction[property_name]

function check_required(o::IoK8sApiCoreV1HTTPGetAction)
    (getproperty(o, Symbol("port")) === nothing) && (return false)
    true
end

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


text:
"""
    Node(name, value, first, last, children, ruleType)
Node in the Abstract Syntax Tree generated by parsing according to a specified grammar.

`name` gives the name of the rule in the grammar.
`value` contains the complete String represented by this node and all its children.
`first` position of the first character matched by this node and its children within the parsed string
`last` position of the last character matched by this node and its children within the parsed string
`children`: Array of child-nodes
Creation adds the additional field `sym` containing the name as a Symbol.
"""
immutable Node
    name::AbstractString
    value::AbstractString
    first::Int
    last::Int
    children::Array #::Array{Node}
    ruleType::Type
end
Node(name::AbstractString, value::AbstractString, first::Int, last::Int, typ) =
    Node(name, value, first, last, [], typ)

########
# SHOW #
########

show{T}(io::IO, val::T, indent) = println(io, "$val ($(typeof(val)))")

# by default don't show anything
displayValue{T <: Rule}(value, ::Type{T}) = ""
# except for terminals and regex
displayValue(value, ::Type{Terminal}) = "'$value',"
displayValue(value, ::Type{RegexRule}) = "'$value',"
displayValue(value, ::Type{IntegerRule}) = "$value,"
displayValue(value, ::Type{FloatRule}) = "$value,"

function show(io::IO, node::Node, indent)
  println(io, "node($(node.name)) {$(displayValue(node.value, node.ruleType))$(node.ruleType)}")
  if isa(node.children, Array)
    for (i, child) in enumerate(node.children)
      print(io, "  "^indent)
      print(io, "$i: ")
      show(io, child, indent+1)
    end
  else
    print(io, "  "^(indent+1))
    show(io, node.children, indent+1)
  end
end

show(io::IO, node::Node) = show(io, node, 0)


text:
using StochasticDiffEq, DiffEqDevTools, Test
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_additivesystem

prob = prob_sde_additivesystem
prob = SDEProblem(prob.f,prob.g,prob.u0,(0.0,0.1),prob.p)

reltols = 1.0./10.0.^(1:4)
abstols = reltols#[0.0 for i in eachindex(reltols)]
setups = [Dict(:alg=>SRIW1())
          Dict(:alg=>EM(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>RKMil(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRIW1(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRA1(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRA1())
          ]
_names = ["SRIW1","EM","RKMil","SRIW1 Fixed","SRA1 Fixed","SRA1"]
test_dt = 0.1
wp = WorkPrecisionSet(prob,abstols,reltols,setups,test_dt;
                                     numruns=5,names=_names,error_estimate=:l2)

se = get_sample_errors(prob,setups[1],numruns=100,solution_runs=100)
se = get_sample_errors(prob,setups[1],numruns=[5,10,25,50,100,1000],solution_runs=100)

println("Now weak error without analytical solution")

prob2 = SDEProblem((du,u,p,t)->prob.f(du,u,p,t),prob.g,prob.u0,(0.0,0.1),prob.p)
test_dt = 1/10^4
appxsol_setup = Dict(:alg=>SRIW1(),:abstol=>1e-4,:reltol=>1e-4)
wp = WorkPrecisionSet(prob2,abstols,reltols,setups,test_dt;
                                     appxsol_setup = appxsol_setup,
                                     numruns=5,names=_names,error_estimate=:weak_final)

println("Get sample errors")

se2 = get_sample_errors(prob2,setups[1],test_dt,appxsol_setup = appxsol_setup,
                        numruns=[5,10,25,50,100],solution_runs=20)

@test all(se[1:5]-se2 .< 1e-1)

# Ensemble Problem with non-commutative noise process

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

u₀ = [1.0,1.0]
function f2!(du,u,p,t)
  du[1] = -273//512*u[1]
  du[2] = -1//160*u[1]-(-785//512+sqrt(2)/8)*u[2]
  return nothing
end
function g2!(du,u,p,t)
  du[1,1] = 1//4*u[1]
  du[1,2] = 1//16*u[1]
  du[2,1] = (1-2*sqrt(2))/4*u[1]
  du[2,2] = 1//10*u[1]+1//16*u[2]
  return nothing
end
dts = 1 .//2 .^(3:-1:0)
tspan = (0.0,3.0)

h2(z) = z^2 # but apply it only to u[1]

prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2))

numtraj = Int(1e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h2(sol[end][1]),false),
        prob_func = prob_func
        )

reltols = 1.0 ./ 4.0 .^ (1:4)
abstols = reltols#[0.0 for i in eachindex(reltols)]
setups = [
        Dict(:alg=>EM(),:dts=>dts),
        Dict(:alg=>SimplifiedEM(),:dts=>dts),
        Dict(:alg=>DRI1(),:dts=>dts, :adaptive=>false)
        ]
test_dt = 1//1000
appxsol_setup = Dict(:alg=>EM(), :dt=>test_dt)

# without analytical expectation value
wp1 = @time WorkPrecisionSet(ensemble_prob,abstols,reltols,setups,test_dt;
                            maxiters = 1e7,
                            verbose=false,save_everystep=false, save_start=false,
                            appxsol_setup = appxsol_setup,
                            trajectories=numtraj, error_estimate=:weak_final)

# with analytical expectation value
wp2 = @time WorkPrecisionSet(ensemble_prob,abstols,reltols,setups,test_dt;
                            maxiters = 1e7,
                            verbose=false,save_everystep=false, save_start=false,
                            appxsol_setup = appxsol_setup,expected_value=exp(-3.0),
                            trajectories=numtraj, error_estimate=:weak_final)

err1 = [wp1.wps[i].errors for i=1:length(setups)]
err2 = [wp2.wps[i].errors for i=1:length(setups)]

@test isapprox(err1, err2, atol=1e-4)


text:
#===============================================================================
                    New definitions of introduced updates
===============================================================================#

#===============================================================================
                            Imputation of paths
===============================================================================#
"""
    struct PathImputation{T} <: MCMCDiffusionImputation
        ρs::Vector{Vector{T}}
        aux_laws
    end

Update type. An indicator for imputation of unobserved path segments. `ρs` are
the memory parameters for the preconditioned Crank–Nicolson scheme and
`aux_laws` are the laws for the auxiliary diffusions.

    PathImputation(ρ::T, P) where T<:Number

Base constructor v1. Initialize all `ρs` with the same value `ρ` and set
auxiliary laws to `P`.

    PathImputation(ρ::T, P) where T<:Vector
Base constructor v2. Initialize each recording with its own value of `ρ` and set
auxiliary laws to `P`.
"""
struct PathImputation{T,A} <: MCMCDiffusionImputation
    ρs::Vector{Vector{T}}
    aux_laws
    adpt::A

    function PathImputation(ρ::T, P; adpt=NoAdaptation()) where T<:Number
        _adpt, A = init_path_imp_adpt(adpt)
        new{T,A}([[ρ]], P, _adpt)
    end

    function PathImputation(ρ::T, P; adpt=NoAdaptation()) where T<:Vector
        _adpt, A = init_path_imp_adpt(adpt)
        new{eltype(ρ),A}([ρ], P, _adpt)
    end

    function PathImputation(
            ρs::T, P; adpt=NoAdaptation()
        ) where T<:Vector{<:Vector}
        _adpt, A = init_path_imp_adpt(adpt)
        new{eltype(ρs[1]),A}(ρs, P, _adpt)
    end
end

init_update!(updt::eMCMC.MCMCUpdate, block_layout) = nothing

function init_update!(updt::PathImputation, block_layout)
    num_recordings = length(block_layout)
    Δ = num_recordings - length(updt.ρs)
    if Δ > 0
        for i in 1:Δ
            push!(updt.ρs, [updt.ρs[1][1]])
        end
    end
    for i in 1:length(updt.ρs)
        update_ρ!(updt.ρs[i], block_layout[i])
    end
    resize_adpt!(updt.adpt, num_recordings)
end

function update_ρ!(ρs::Vector, block_layout)
    N = length(block_layout[2])
    if length(ρs) < N
        for i in 1:(N-length(ρs))
            push!(ρs, ρs[1])
        end
    end
end

struct StartingPointsImputation <: MCMCDiffusionImputation end

struct StartingPointsLangevinImputation <: MCMCDiffusionImputation end

#===============================================================================
                        Conjugate Gaussian updates
===============================================================================#

struct DiffusionConjugGsnUpdate{K,TP} <: eMCMC.MCMCConjugateParamUpdate
    coords::K
    prior::TP
    adpt::NoAdaptation

    function DiffusionConjugGsnUpdate(coords, prior)
        if typeof(prior) <: Normal
            prior = Gaussian( [prior.μ], reshape(prior.σ, (1,1)) )
        elseif typeof(prior) <: MvNormal
            prior = Gaussian(prior.μ, prior.Σ.mat)
        else
            @assert typeof(prior) <: Gaussian
        end

        if typeof(coords) <: Number
            coords = [coords]
        end
        new{typeof(coords), typeof(prior)}(coords, prior, NoAdaptation())
    end
end

struct ConjugGsnHelper{T,N,K}
    relevant_recordings::Vector{Int64}
    sub_μ_idx::Vector{K}
    θ_names::Vector{Val}
    θᶜ_names::Vector{Val}

    function ConjugGsnHelper(
            updt::DiffusionConjugGsnUpdate,
            global_ws,
            PPP,
            lpn::LocalUpdtParamNames;
            computation_type=:outofplace
        )
        _glob_pnames = global_ws.pnames[updt.coords]
        relevant_recordings = Int64[]
        sub_μ_idx = Vector{Vector{Int64}}(undef, 0)
        θᶜ_names = Val[]
        θ_names = Val[]
        all_θ_names = Symbol[]

        for (i, PP) in enumerate(PPP)
            P = first(PP).P_target
            if length(lpn.θ_local_names) > 0
                push!(relevant_recordings, i)

                # sort by the local index of θ
                θ_loc_names = sort(collect(lpn.θ_local_names[i]), by=x->x[1])
                sub_idx = getindex.(θ_loc_names, 1)
                push!(sub_μ_idx, sub_idx)

                sub_name = getindex.(θ_loc_names, 2)
                push!(
                    θᶜ_names,
                    Val(
                        Tuple(
                            filter(
                                pn->(
                                    !(pn in sub_name) &&
                                    !DD.ignore_for_cu(Val(pn), P)
                                ),
                                DD.parameter_names(P)
                            )
                        )
                    )
                )
                push!( θ_names, Val(Tuple(sub_name)) )
                append!(all_θ_names, sub_name)
            else
                push!(sub_μ_idx, Int64[])
                push!(θᶜ_names, Val[])
                push!(θ_names, Val[])
            end
        end
        u_θ_names = unique(all_θ_names)
        sub_μ_idx = transform_μ_idx(Val(computation_type), sub_μ_idx)
        new{computation_type, length(u_θ_names), eltype(sub_μ_idx)}(
            relevant_recordings,
            sub_μ_idx,
            θ_names,
            θᶜ_names,
        )
    end

    ConjugGsnHelper() = new{:none, 0, Nothing}()
end

transform_μ_idx(::Val{:inplace}, sub_μ_idx) = sub_μ_idx

function transform_μ_idx(v::Val{:outofplace}, sub_μ_idx)
    map(i->transform_μ_idx(v, i, Val(length(i))), sub_μ_idx)
end

function transform_μ_idx(::Val{:outofplace}, sub_μ_idx, ::Val{N}) where N
    SVector{N,Int64}(sub_μ_idx)
end

function μ_init(cu_helper::ConjugGsnHelper{:outofplace,N}) where N
    zero(MVector{N,Float64})
end

function W_init(cu_helper::ConjugGsnHelper{:outofplace,N}) where N
    zero(MMatrix{N,N,Float64})
end

function μ_init(cu_helper::ConjugGsnHelper{:outofplace}, ::SVector{N}) where N
    zero(SVector{N,Float64})
end

μ_init(cu_helper::ConjugGsnHelper{:inplace,N}) where N = zeros(Float64, N)
W_init(cu_helper::ConjugGsnHelper{:inplace,N}) where N = zeros(Float64, (N,N))

function μ_init(cu_helper::ConjugGsnHelper{:inplace,N}, v) where N
    zeros(Float64, length(v))
end

@generated function conjug_θᶜ(
        P_target, pnames::Val{T}, ::ConjugGsnHelper{:outofplace}
    ) where T
    N = length(T)
    data = Expr(
        :tuple,
        1.0,
        (:(getfield(P_target, T[$i])) for i in 1:N)...,
    )
    vec = Expr(:call, MVector{N+1}, data)
end

@generated function conjug_θᶜ(
        P_target, pnames::Val{T}, ::ConjugGsnHelper{:inplace}
    ) where T
    N = length(T)
    data = Expr(
        :tuple,
        1.0,
        (:(getfield(P_target, T[$i])) for i in 1:N)...,
    )
    vec = Expr(:call, :collect, data)
end

function eMCMC.proposal!(
        updt::DiffusionConjugGsnUpdate,
        global_ws,
        local_ws::eMCMC.LocalWorkspace,
        step,
    )
    state°(local_ws) .= conjugate_draw(
        local_ws.cu_helper,
        local_ws.sub_ws_diff.P,
        updt.prior,
        local_ws.sub_ws_diff.XX,
    )
end

function conjugate_draw(cu_helper, PP, prior, XX)
    μ = μ_init(cu_helper)
    W = W_init(cu_helper)
    for i in cu_helper.relevant_recordings
        μi = μ_init(cu_helper, cu_helper.sub_μ_idx[i])
        Wi = μi*μi'
        P = first(PP[i]).P_target
        θᶜ = conjug_θᶜ(P, cu_helper.θᶜ_names[i], cu_helper)

        μi, Wi = compute_μ_and_W(
            μi, Wi, cu_helper.θ_names[i], cu_helper.θᶜ_names[i], θᶜ, P, XX[i],
            cu_helper
        )

        sidx = cu_helper.sub_μ_idx[i]
        view(μ, sidx) .+= μi
        view(W, sidx, sidx) .+= Wi
    end
    μ_post, Σ_post = compute_posterior_μ_and_Σ(μ, W, prior)
    rand(Gaussian(μ_post, Σ_post))
end

function compute_μ_and_W(
        μ, W, θ_names, θᶜ_names, θᶜ, P_target::S, XX, ::ConjugGsnHelper{:outofplace};
        nnh = DD.num_non_hypo(S)
    ) where S
    for X in XX
        num_pts = length(X)
        for i in 1:num_pts-1
            φₜ = φ(θ_names, X.t[i], X.x[i], P_target, nnh)
            φᶜₜ = φᶜ(θᶜ_names, θᶜ, X.t[i], X.x[i], P_target, nnh)
            Γ⁻¹ = DD.hypo_a_inv(X.t[i], X.x[i], P_target)
            dt = X.t[i+1] - X.t[i]
            dy = DD.nonhypo(X.x[i+1], P_target) - DD.nonhypo(X.x[i], P_target)
            μ = μ + φₜ'*Γ⁻¹*dy - φₜ'*Γ⁻¹*φᶜₜ*dt
            W = W + φₜ'*Γ⁻¹*φₜ*dt
        end
    end
    μ, W
end

function compute_μ_and_W(
        μ, W, θ_names, P_target, XX, θᶜ, ::ConjugGsnHelper{:inplace}
    )
    error("not implemented, need buffers...")
end

function compute_posterior_μ_and_Σ(μ, W, prior)
    Σ = Symmetric(inv(W + inv(prior.Σ)))
    Σ * (μ .+ prior.Σ\prior.μ), Σ
end

const _GENERATED = Val(:gen)

function φ(v::Val, t, x, P::S, num_non_smooth_coords=DD.num_non_hypo(S)) where S
    φ(_GENERATED, v, t, x, P, num_non_smooth_coords)
end

@generated function φ(::Val{:gen}, ::Val{Ts}, t, x, P::S, ::Val{N}) where {Ts,S,N}
	args = [Expr(:call, :Val, :(Ts[$i])) for i in 1:length(Ts)]
    data = Expr(
        :call,
        :tuplejoin,#TODO try to do tuple join not in the expression
        ( :(DD.phi($a, t, x, P)) for a in args )...
    )
    mat = Expr(:call, SMatrix{N,length(Ts)}, data)
    return mat
end

"""
    φᶜ(v, θᶜ, t, x, P::S, n=DD.num_non_hypo(S)) where S

Remainder term in the drift after removing φ'θ from it.
"""
function φᶜ(v, θᶜ, t, x, P::S, n=DD.num_non_hypo(S)) where S
    φᶜlinear(v, t, x, P, n) * θᶜ
end

function φᶜlinear(v, t, x, P::S, nnh=DD.num_non_hypo(S)) where S
    φᶜlinear(_GENERATED, v, t, x, P, nnh)
end

@generated function φᶜlinear(::Val{:gen}, ::Val{Ts}, t, x, P::S, ::Val{N}) where {Ts,S,N}
    args = [Expr(:call, :Val, :(Ts[$i])) for i in 1:length(Ts)]
    data = Expr(
        :call,
        :tuplejoin,
        :( DD.phi(Val(:intercept), t, x, P) ),
        ( :(DD.phi($a, t, x, P)) for a in args )...
    )
    mat = Expr(:call, SMatrix{N,length(Ts)+1}, data)
    return mat
end

@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = tuplejoin(tuplejoin(x, y), z...)
#TODO at the accept reject step it needs to accept a boolean not only a vec of booleans


text:
## SparseTensorDOK: dictionary of keys format

"""
    SparseTensorDOK{Tv,Ti<:Integer,N} <: AbstractSparseTensor{Tv,Ti,N}

`N`-dimensional sparse tensor stored in the **D**ictionary **O**f **K**eys format.
Elements are stored as a dictionary mapping indices (using type `Ti`) to values (of type `Tv`).

Fields:
+ `dims::Dims{N}`               : tuple of dimensions
+ `dict::Dict{NTuple{N,Ti},Tv}` : dictionary mapping indices to values
"""
struct SparseTensorDOK{Tv,Ti<:Integer,N} <: AbstractSparseTensor{Tv,Ti,N}
    dims::Dims{N}                   # Dimensions
    dict::Dict{NTuple{N,Ti},Tv}     # Dictionary

    function SparseTensorDOK{Tv,Ti,N}(dims::Dims{N},
                            dict::Dict{NTuple{N,Ti},Tv}) where {Tv,Ti<:Integer,N}
        check_Ti(dims, Ti)
        foreach(ind -> checkbounds_dims(dims, ind...), keys(dict))
        return new(dims, dict)
    end
end
SparseTensorDOK(dims::Dims{N}, dict::Dict{NTuple{N,Ti},Tv}) where {Tv,Ti<:Integer,N} =
    SparseTensorDOK{Tv,Ti,N}(dims, dict)

"""
    SparseTensorDOK{Tv,Ti<:Integer}(undef, dims)
    SparseTensorDOK{Tv,Ti<:Integer,N}(undef, dims)

Construct an uninitialized `N`-dimensional `SparseTensorDOK`
with indices using type `Ti` and elements of type `Tv`.
Here uninitialized means it has no stored entries.

Here `undef` is the `UndefInitializer`. If `N` is supplied,
then it must match the length of `dims`.

# Examples
```julia-repl
julia> A = SparseTensorDOK{Float64, Int8, 3}(undef, (2, 3, 4)) # N given explicitly
2×3×4 SparseTensorDOK{Float64, Int8, 3} with 0 stored entries

julia> B = SparseTensorDOK{Float64, Int8}(undef, (4,)) # N determined by the input
4-element SparseTensorDOK{Float64, Int8, 1} with 0 stored entries

julia> similar(B, 2, 4, 1) # use typeof(B), and the given size
2×4×1 SparseTensorDOK{Float64, Int8, 3} with 0 stored entries
```
"""
SparseTensorDOK{Tv,Ti,N}(::UndefInitializer, dims::Dims{N}) where {Tv,Ti<:Integer,N} =
    SparseTensorDOK(dims, Dict{NTuple{N,Ti},Tv}())
SparseTensorDOK{Tv,Ti}(::UndefInitializer, dims::Dims{N}) where {Tv,Ti<:Integer,N} =
    SparseTensorDOK{Tv,Ti,N}(undef, dims)

"""
    SparseTensorDOK(Ti, A::AbstractArray)

Convert an AbstractArray `A` into a `SparseTensorDOK`
with indices using type `Ti`.

# Examples
```julia-repl
julia> A = SparseTensorDOK(Int8, Float16[1.1 0.0 0.0; 2.1 0.0 2.3])
2×3 SparseTensorDOK{Float16, Int8, 2} with 3 stored entries:
  [1, 1]  =  1.1
  [2, 1]  =  2.1
  [2, 3]  =  2.3

julia> B = SparseTensorDOK(Int16, Float16[1.1 0.0 0.0; 2.1 0.0 2.3])
2×3 SparseTensorDOK{Float16, Int16, 2} with 3 stored entries:
  [1, 1]  =  1.1
  [2, 1]  =  2.1
  [2, 3]  =  2.3
```
"""
function SparseTensorDOK(Ti::Type{<:Integer}, A::AbstractArray)
    Tv, N = eltype(A), ndims(A)
    dims = size(A)
    nzidx = findall(!iszero, A)
    inds = convert(Vector{NTuple{N,Ti}}, CartesianIndices(A)[nzidx])
    vals = convert(Vector{Tv}, A[nzidx])
    SparseTensorDOK{Tv,Ti,N}(dims, Dict(inds .=> vals))
end

## Minimal AbstractArray interface

size(A::SparseTensorDOK) = A.dims

function getindex(A::SparseTensorDOK{Tv,<:Integer,N}, I::Vararg{Int,N}) where {Tv,N}
    @boundscheck checkbounds(A, I...)
    return get(A.dict, I, zero(Tv))
end

function setindex!(A::SparseTensorDOK{Tv,Ti,N}, v, I::Vararg{Int,N}) where {Tv,Ti<:Integer,N}
    @boundscheck checkbounds(A, I...)
    ind, val = convert(NTuple{N,Ti}, I), convert(Tv, v)
    if !iszero(val) || haskey(A.dict, ind)
        A.dict[ind] = val
    end
    return A
end

IndexStyle(::Type{<:SparseTensorDOK}) = IndexCartesian()

## Overloads for specializing outputs

similar(::SparseTensorDOK{<:Any,Ti}, ::Type{Tv}, dims::Dims{N}) where {Tv,Ti<:Integer,N} =
    SparseTensorDOK{Tv,Ti,N}(undef, dims)

## AbstractSparseTensor interface

function dropstored!(f::Function, A::SparseTensorDOK)
    filter!(p -> !f(p.second), A.dict)
    return A
end

numstored(A::SparseTensorDOK) = length(A.dict)
storedindices(A::SparseTensorDOK) = collect(keys(A.dict))
storedvalues(A::SparseTensorDOK) = collect(values(A.dict))
storedpairs(A::SparseTensorDOK) = pairs(A.dict)


text:
pat = r"(?:n(?:w|e)|e|w|s(?:w|e))"
tiles = open("tiles.txt") do f
    map(tile -> map(x -> x.match, tile), eachmatch.(pat, readlines(f)))
end

flipped = Dict()
for tile in tiles
    coord = (q = 0, r = 0)
    for d in tile
        q = 0
        r = 0
        if d == "nw"
            r = -1
        elseif d == "ne"
            q = 1
            r = -1
        elseif d == "e"
            q = 1
        elseif d == "se"
            r = 1
        elseif d == "sw"
            r = 1
            q = -1
        elseif d == "w"
            q = -1
        end
        coord = (q = coord.q + q, r = coord.r + r)
    end
    if haskey(flipped, coord)
        delete!(flipped, coord)
    else
        flipped[coord] = 1
    end
end

println(sum(values(flipped)))

function draw_map(flipped)
    min_q = 0
    max_q = 0
    min_r = 0
    max_r = 0
    for (q,r) in keys(flipped)
        if q > max_q
            max_q = q
        elseif q < min_q
            min_q = q
        end
        if r > max_r
            max_r = r
        elseif r < min_r
            min_r = r
        end
    end
    for r in min_r:max_r
        print(repeat(" ", r-min_r))
        for q in min_q:max_q
            if haskey(flipped, (q = q, r = r))
                print("$(flipped[(q = q, r = r)]) ")
            else
                print("o ")
            end
        end
        println()
    end
end

for i in 1:100
    global flipped
    # draw_map(flipped)
    # add adjacent white tiles
    for (coord, tile) in flipped
        if tile == 0
            continue
        end
        adj = [(0, -1), (+1, -1), (+1, 0), (0, +1), (-1, +1), (-1, 0)]
        for (q, r) in adj
            if !haskey(flipped, (q = coord.q + q, r = coord.r + r))
                flipped[(q = coord.q + q, r = coord.r + r)] = 0
            end
        end
    end

    # draw_map(flipped)
    result = deepcopy(flipped)
    for (coord, tile) in flipped
        blacks = 0
        adj = [(0, -1), (+1, -1), (+1, 0), (0, +1), (-1, +1), (-1, 0)]
        for (q, r) in adj
            if haskey(flipped, (q = coord.q + q, r = coord.r + r)) && flipped[(q = coord.q + q, r = coord.r + r)] == 1
                blacks += 1
            end
        end
        if tile == 1 && (blacks == 0 || blacks > 2)
            result[coord] = 0
        elseif tile == 0 && blacks == 2
            result[coord] = 1
        end
        # println("$(coord): $(tile) => $(result[coord]) ($blacks)")
    end

    flipped = result

    # remove empty
    for (coord, tile) in flipped
        if tile == 0
            delete!(flipped, coord)
        end
    end

    # if i < 10 || i % 10 == 0
    #     println(i,": ",sum(values(flipped)))
    # end
end

println(sum(values(flipped)))


text:
using Folds
using Folds.Testing: test_uses_threads
using Test

@testset "ThreadedEx(basesize = 1)" begin
    test_uses_threads(ThreadedEx(basesize = 1))
end

@testset "ThreadedEx()" begin
    test_uses_threads(ThreadedEx())
end


text:
using Interpolations

struct MultiDomain{TF <: Union{Float32, Float64}}
    enable_nudge::Bool
    nudge_time::TF

    u_nudge::Array{TF, 3}
    v_nudge::Array{TF, 3}
    w_nudge::Array{TF, 3}
    s_nudge::Array{TF, 3}
end

# Constructor for multidomain struct.
function MultiDomain(g::Grid, settings, TF)
    enable_nudge = settings["enable_nudge"]
    if enable_nudge
        nudge_time = settings["nudge_time"]
        return MultiDomain{TF}(
            enable_nudge,
            nudge_time,
            zeros(g.icells, g.jcells, g.kcells),
            zeros(g.icells, g.jcells, g.kcells),
            zeros(g.icells, g.jcells, g.kcells),
            zeros(g.icells, g.jcells, g.kcells))
    else
        nudge_time = 1e12
        return MultiDomain{TF}(
            enable_nudge,
            nudge_time,
            zeros(0, 0, 0),
            zeros(0, 0, 0),
            zeros(0, 0, 0),
            zeros(0, 0, 0))
    end
end

# Transfer the state
function transfer_state!(f_d::Fields, f_s::Fields, g_d::Grid, g_s::Grid)
    interp = interpolate((g_s.xh, g_s.y, g_s.z), f_s.u, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    f_d.u[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.xh[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])

    interp = interpolate((g_s.x, g_s.yh, g_s.z), f_s.v, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    f_d.v[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.x[g_d.is:g_d.ie], g_d.yh[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])

    interp = interpolate((g_s.x, g_s.y, g_s.zh), f_s.w, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    f_d.w[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.keh] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.zh[g_d.ks:g_d.keh])

    interp = interpolate((g_s.x, g_s.y, g_s.z), f_s.s, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    f_d.s[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])

    interp = interpolate((g_s.x, g_s.y, g_s.z), f_s.p, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    f_d.p[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])

    interp = interpolate((g_s.x, g_s.y), f_s.s_bot, (Gridded(Linear()), Gridded(Linear())))
    f_d.s_bot[g_d.is:g_d.ie, g_d.js:g_d.je] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je])

    interp = interpolate((g_s.x, g_s.y), f_s.s_gradbot, (Gridded(Linear()), Gridded(Linear())))
    f_d.s_gradbot[g_d.is:g_d.ie, g_d.js:g_d.je] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je])

    interp = interpolate((g_s.x, g_s.y), f_s.s_top, (Gridded(Linear()), Gridded(Linear())))
    f_d.s_top[g_d.is:g_d.ie, g_d.js:g_d.je] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je])

    interp = interpolate((g_s.x, g_s.y), f_s.s_gradtop, (Gridded(Linear()), Gridded(Linear())))
    f_d.s_gradtop[g_d.is:g_d.ie, g_d.js:g_d.je] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je])
end

## Define the upsample functions.
function make_inner_loop!(ex_list, ii, jj, kk, ifac, jfac, kfac, ioff, joff, koff, is_top)
    push!(ex_list, :(i_hi = $ifac*i+$ii+is_hi), :(j_hi = $jfac*j+$jj+js_hi))
    if is_top
        push!(ex_list, :(k_hi = ke_hi))
    else
        push!(ex_list, :(k_hi = $kfac*k+$kk+ks_hi))
    end

    i_pos = -1/2 + ioff + (1/2 - ioff + ii)/ifac
    j_pos = -1/2 + joff + (1/2 - joff + jj)/jfac
    k_pos = -1/2 + koff + (1/2 - koff + kk)/kfac

    i_lo_off = floor(Int, i_pos)
    j_lo_off = floor(Int, j_pos)
    k_lo_off = floor(Int, k_pos)

    if is_top
        push!(ex_list, :(i_lo = i + $i_lo_off + is_lo), :(j_lo = j + $j_lo_off + js_lo), :(k_lo = ke_lo))
    else
        push!(ex_list, :(i_lo = i + $i_lo_off + is_lo), :(j_lo = j + $j_lo_off + js_lo), :(k_lo = k + $k_lo_off + ks_lo))
    end

    fi = mod(i_pos, 1)
    fj = mod(j_pos, 1)
    fk = mod(k_pos, 1)

    coef = zeros(2, 2, 2)
    coef[1, 1, 1] = (((1-fi)*fj+fi-1)*fk+(fi-1)*fj-fi+1)
    coef[2, 1, 1] = ((fi*fj-fi)*fk-fi*fj+fi)
    coef[1, 2, 1] = ((fi-1)*fj*fk+(1-fi)*fj)
    coef[2, 2, 1] = (fi*fj-fi*fj*fk)
    coef[1, 1, 2] = ((fi-1)*fj-fi+1)*fk
    coef[2, 1, 2] = (fi-fi*fj)*fk
    coef[1, 2, 2] = (1-fi)*fj*fk
    coef[2, 2, 2] = fi*fj*fk

    kc_end = is_top ? 1 : 2

    rhs_list = []
    for kc in 1:kc_end, jc in 1:2, ic in 1:2
        if !(coef[ic, jc, kc] ≈ 0)
            push!(rhs_list, :( $(coef[ic, jc, kc])*lo[i_lo+$(ic-1), j_lo+$(jc-1), k_lo+$(kc-1)] ))
        end
    end

    rhs = Expr(:call, :+, rhs_list...)

    push!(ex_list, :(hi[i_hi, j_hi, k_hi] = $rhs))
end

macro upsample_lin_fast(suffix, ifac, jfac, kfac, ioff, joff, koff, add_top)
    ex_inner = []
    for kk in 0:kfac-1, jj in 0:jfac-1, ii in 0:ifac-1
        make_inner_loop!(ex_inner, ii, jj, kk, ifac, jfac, kfac, ioff, joff, koff, false)
    end
    ex_inner_block = quote
        Threads.@threads for k in 0:ktot-1
            for j in 0:jtot-1
                @inbounds @simd for i in 0:itot-1
                    $(Expr(:block, ex_inner...))
                end
            end
        end
    end

    if add_top
        ex_top = []
        kk = 0
        for jj in 0:jfac-1, ii in 0:ifac-1
            make_inner_loop!(ex_top, ii, jj, kk, ifac, jfac, kfac, ioff, joff, koff, true)
        end
        ex_top_block = quote
            Threads.@threads for j in 0:jtot-1
                @inbounds @simd for i in 0:itot-1
                    $(Expr(:block, ex_top...))
                end
            end
        end
    else
        ex_top_block = :()
    end

    name = Symbol(@sprintf "upsample_lin_%s!" suffix)
    ex = quote
        function $name(hi, lo, itot, jtot, ktot, is_hi, js_hi, ks_hi, ke_hi, is_lo, js_lo, ks_lo, ke_lo)
            $ex_inner_block
            $ex_top_block
        end
    end

    return esc(ex)
end

@upsample_lin_fast("222_u", 2, 2, 2, 0.5, 0  ,   0, false)
@upsample_lin_fast("222_v", 2, 2, 2,   0, 0.5,   0, false)
@upsample_lin_fast("222_w", 2, 2, 2,   0, 0  , 0.5, true )
@upsample_lin_fast("222_s", 2, 2, 2,   0, 0  ,   0, false)
@upsample_lin_fast("333_u", 3, 3, 3, 0.5, 0  ,   0, false)
@upsample_lin_fast("333_v", 3, 3, 3,   0, 0.5,   0, false)
@upsample_lin_fast("333_w", 3, 3, 3,   0, 0  , 0.5, true )
@upsample_lin_fast("333_s", 3, 3, 3,   0, 0  ,   0, false)
@upsample_lin_fast("444_u", 4, 4, 4, 0.5, 0  ,   0, false)
@upsample_lin_fast("444_v", 4, 4, 4,   0, 0.5,   0, false)
@upsample_lin_fast("444_w", 4, 4, 4,   0, 0  , 0.5, true )
@upsample_lin_fast("444_s", 4, 4, 4,   0, 0  ,   0, false)

function calc_nudge_fields!(md::MultiDomain, f_d::Fields, f_s::Fields, g_d::Grid, g_s::Grid)
    if !md.enable_nudge
        return
    end

    upsample_ratio = (g_d.itot, g_d.jtot, g_d.ktot) .÷ (g_s.itot, g_s.jtot, g_s.ktot)
    if upsample_ratio == (2, 2, 2)
        upsample_lin_222_u!(
            md.u_nudge, f_s.u, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

        upsample_lin_222_v!(
            md.v_nudge, f_s.v, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

        upsample_lin_222_w!(
            md.w_nudge, f_s.w, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.keh, g_s.is, g_s.js, g_s.ks, g_s.keh)

        upsample_lin_222_s!(
            md.s_nudge, f_s.s, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

    elseif upsample_ratio == (3, 3, 3)
        upsample_lin_333_u!(
            md.u_nudge, f_s.u, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

        upsample_lin_333_v!(
            md.v_nudge, f_s.v, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

        upsample_lin_333_w!(
            md.w_nudge, f_s.w, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.keh, g_s.is, g_s.js, g_s.ks, g_s.keh)

        upsample_lin_333_s!(
            md.s_nudge, f_s.s, g_s.itot, g_s.jtot, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

    elseif upsample_ratio == (4, 4, 4)
        upsample_lin_444_u!(
            md.u_nudge, f_s.u, g_s.imax, g_s.jmax, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

        upsample_lin_444_v!(
            md.v_nudge, f_s.v, g_s.imax, g_s.jmax, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

        upsample_lin_444_w!(
            md.w_nudge, f_s.w, g_s.imax, g_s.jmax, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.keh, g_s.is, g_s.js, g_s.ks, g_s.keh)

        upsample_lin_444_s!(
            md.s_nudge, f_s.s, g_s.imax, g_s.jmax, g_s.ktot,
            g_d.is, g_d.js, g_d.ks, g_d.ke, g_s.is, g_s.js, g_s.ks, g_s.ke)

    else
        println("WARNING: resorting to slow interpolations in nudging.")
        @sync begin
            Threads.@spawn begin
                interp = interpolate((g_s.xh, g_s.y, g_s.z), f_s.u, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
                md.u_nudge[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.xh[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])
            end

            Threads.@spawn begin
                interp = interpolate((g_s.x, g_s.yh, g_s.z), f_s.v, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
                md.v_nudge[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.x[g_d.is:g_d.ie], g_d.yh[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])
            end

            Threads.@spawn begin
                interp = interpolate((g_s.x, g_s.y, g_s.zh), f_s.w, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
                md.w_nudge[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.keh] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.zh[g_d.ks:g_d.keh])
            end

            Threads.@spawn begin
                interp = interpolate((g_s.x, g_s.y, g_s.z), f_s.s, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
                md.s_nudge[g_d.is:g_d.ie, g_d.js:g_d.je, g_d.ks:g_d.ke] .= interp(g_d.x[g_d.is:g_d.ie], g_d.y[g_d.js:g_d.je], g_d.z[g_d.ks:g_d.ke])
            end
        end
    end
end

function nudge_kernel!(a_tend, a, a_nudge, c_nudge, is, ie, js, je, ks, ke)
    @fast3d begin
        @fd (a_tend, a, a_nudge) a_tend -= c_nudge * (a - a_nudge)
    end
end

function calc_nudge_tend!(f::Fields, g::Grid, md::MultiDomain)
    if !md.enable_nudge
        return
    end

    c_nudge = 1 / md.nudge_time

    nudge_kernel!(f.u_tend, f.u, md.u_nudge, c_nudge, g.is, g.ie, g.js, g.je, g.ks, g.ke )
    nudge_kernel!(f.v_tend, f.v, md.v_nudge, c_nudge, g.is, g.ie, g.js, g.je, g.ks, g.ke )
    nudge_kernel!(f.w_tend, f.w, md.w_nudge, c_nudge, g.is, g.ie, g.js, g.je, g.ks, g.keh)
    nudge_kernel!(f.s_tend, f.s, md.s_nudge, c_nudge, g.is, g.ie, g.js, g.je, g.ks, g.ke )
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

dependencies = [
    "build_Zlib.v1.2.11.jl",
    "build_libpng.v1.0.0.jl",
    "build_libjpeg.v9.0.0-b.jl",
    "build_libtiff.v4.0.9.jl"
]

for elem in dependencies
    # it's a bit faster to run the build in an anonymous module instead of
    # starting a new julia process
    m = Module(:__anon__)
    Core.include(m, (joinpath(@__DIR__, elem)))
end

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libMagickWand"], :libwand),
]
# Download binaries from hosted location
version = v"6.9.10-10"
bin_prefix = "https://github.com/JuliaIO/ImageMagickBuilder/releases/download/v$version"

# Listing of files generated by BinaryBuilder:

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/ImageMagick.v6.9.10-10.aarch64-linux-gnu.tar.gz", "50809b9ad03f527998acb00d53ec647964c759282e0828d8e5123a3f30ad0a0f"),
    Linux(:aarch64, :musl) => ("$bin_prefix/ImageMagick.v6.9.10-10.aarch64-linux-musl.tar.gz", "24e6198bcb3ed639e815bb422c2681721052cb4212ceb2bfe509e155d9610c83"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/ImageMagick.v6.9.10-10.arm-linux-gnueabihf.tar.gz", "636872cd4a3c4078490164596258d13cddad52fef455e1fc45a8856893bd29d3"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/ImageMagick.v6.9.10-10.arm-linux-musleabihf.tar.gz", "d0c7c882fd4351b1f5f5a655047a193e07c1863dd9d6978b792c76c16d3b359c"),
    Linux(:i686, :glibc) => ("$bin_prefix/ImageMagick.v6.9.10-10.i686-linux-gnu.tar.gz", "9aa4fa4dd5a38cad9d0e99cb73dbbba839abd80222641ca997816ea4881a68cb"),
    Linux(:i686, :musl) => ("$bin_prefix/ImageMagick.v6.9.10-10.i686-linux-musl.tar.gz", "606dcdc8964db83c48ab669c1ac482fa726a728a73f9806ac88bdac3029703fa"),
    Windows(:i686) => ("$bin_prefix/ImageMagick.v6.9.10-10.i686-w64-mingw32.tar.gz", "cf28050759ff1fe56928449110206a23d5313a26bdab778737bc63640c1eefc0"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/ImageMagick.v6.9.10-10.powerpc64le-linux-gnu.tar.gz", "c83b2c9b5ab711d2a54643358017cbd0e29dd8ca9d6e6df81915d1ad82a1c400"),
    MacOS(:x86_64) => ("$bin_prefix/ImageMagick.v6.9.10-10.x86_64-apple-darwin14.tar.gz", "2299a8df97f91e60a8e070439276b8dfea926644bd93e3c9940560989e16e005"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/ImageMagick.v6.9.10-10.x86_64-linux-gnu.tar.gz", "3bcfb55798ef8ebdd291990febac393157c246163053b2504bc3637a162e05dd"),
    Linux(:x86_64, :musl) => ("$bin_prefix/ImageMagick.v6.9.10-10.x86_64-linux-musl.tar.gz", "8b2832b290aa2296fe78bca165fdd2a59c598d730f76227d82882a13d30511bd"),
    FreeBSD(:x86_64) => ("$bin_prefix/ImageMagick.v6.9.10-10.x86_64-unknown-freebsd11.1.tar.gz", "b31be069da6d3020ceeb4ce30e3e262bb784adfbcdfb67637f46b99ef9e396b6"),
    Windows(:x86_64) => ("$bin_prefix/ImageMagick.v6.9.10-10.x86_64-w64-mingw32.tar.gz", "c81df1c8adcd4facda669d8d1277af438f90bdcb61aa17da6d59e1566e9e0ae6"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_im.jl"), products)

function include_deps(name)
    """
    module $name
        import Libdl
        path = joinpath(@__DIR__, $(repr(string("deps_", name, ".jl"))))
        isfile(path) || error("$name wasn't build correctly. Please run Pkg.build(\\\"ImageMagick\\\")")
        include(path)
    end
    using .$name
    """
end

open("deps.jl", "w") do io
    for dep in (:zlib, :png, :jpeg, :tiff, :im)
        println(io, include_deps(dep))
    end
    println(io, """
    libversion() = $(repr(version))
    const libwand = im.libwand

    function check_deps()
        zlib.check_deps()
        png.check_deps()
        jpeg.check_deps()
        tiff.check_deps()
        im.check_deps()
    end
    """)
end


text:
function bleu(hyp, ref)
    N = 4
    hypothesis_length = 0
    references_length = 0
    clipped = Array(Float64, N)
    total = Array(Float64, N)

    # convert single translation to array
    # works for both single and corpus translation
    hcopy, rcopy = copy(hyp), copy(ref)
    if !isa(hyp, Array)
        hcopy = [hcopy]
        rcopy = isa(ref, Array) ? [rcopy] : [[rcopy]]
    end
    length(hcopy) == length(rcopy) || error("dimensions must match")
    translations = zip(hcopy, rcopy)

    for (hypothesis, references) in translations
        # word count dicts
        hypcounts = map(n -> ngram_counts(hypothesis,n), 1:N)
        refscounts = map(r -> map(n -> ngram_counts(r,n), 1:N), references)

        # corpus length
        hyplen = sum(values(hypcounts[1]))
        reflens = map(x->sum(values(x[1])), refscounts)
        differences = sort(map(x -> (abs(x - hyplen), x), reflens))
        hypothesis_length += hyplen
        references_length += differences[1][2]

        for n = 1:N
            c1 = collect(values(hypcounts[n]))
            c2 = map(
                x -> mapreduce(r -> get(r[n], x, 0), max, refscounts),
                keys(hypcounts[n]))
            if length(c1) >= 1
                clipped[n] += sum(map(k->min(c1[k],c2[k]), 1:length(c1)))
                total[n] += sum(c1)
            end
        end
    end

    scores = zeros(N)
    if length(translations) != 0
        scores = map(i -> clipped[i]/total[i], 1:N)
    end

    brevity_penalty = 1
    if hypothesis_length < references_length
        brevity_penalty = exp(1-references_length/hypothesis_length)
    end

    bleuN = Array(Float64, 4)
    for k = 1:4
        bleuN[k] = brevity_penalty * exp(mean(map(log, scores[1:k])))
    end

    (bleuN, brevity_penalty, hypothesis_length, references_length)
end

function countdict(xs)
    counts = Dict()
    for x in xs
        if haskey(counts, x)
            counts[x] += 1
        else
            counts[x] = 1
        end
    end
    return counts
end

get_words(sentence) = map(lowercase, split(sentence, " "))
build_ngram(words, n) = map(i -> join(words[i:i+n-1], " "), 1:length(words)-n+1)
ngram_counts(sentence, n) = countdict(build_ngram(get_words(sentence), n))

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
af = ArrayFire{CPU}()
arr1 = array(af, [[1,1] [1,1] [1,1]])
arr2 = array(af, [[1,1] [1,1] [1,1]])
b = PingPongBatcher(af, arr1, arr2)
println("ok")

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
#!/usr/bin/env julia

# Adapted from bin/cellfishing in the CellFishing.jl repo

using CellFishing
using HDF5
using DocOpt

using Dates
using Random
using SparseArrays
using StatsBase
using Statistics

args = docopt("""
Usage:
    cellfishing [--n-bits=<n-bits> --n-lshashes=<n-lshashes> --infer-stats=<infer-stats> --genes=<genes> --n-neighbors=<n-neighbors> --cutoff=<cutoff> --min-hits=<min-hits> --majority=<majority> --annotation=<annotation> --seed=<seed> --clean=<clean> --subsample-ref=<subsample-ref>] <reference> <query> <output>

CellFishing ― Cell search engine for scRNA-seq.

Options:
    -h --help                        Show this help.
    --version                        Show version.
    --n-bits=<n-bits>                The number of bits for hashing [default: 128].
    --n-lshashes=<n-lshashes>        The number of locality-sensitive hashes [default: 4].
    --infer-stats=<infer-stats>      Infer statistics from the query, the database, or both [default: both].
    --genes=<genes>                  Specify genes to use instead of selecting by CellFishing [default: cf_genes].
    --n-neighbors=<n-neighbors>      The number of nearest neighbors to search for [default: 10].
    --cutoff=<cutoff>...             Hamming distance cutoff [default: 120].
    --min-hits=<min-hits>            Minimal number of hits below which query will be rejected [default: 2].
    --majority=<majority>            Minimal fraction required to pass majority vote [default: 0.5].
    --annotation=<annotation>        Choose annotation column in meta table [default: cell_ontology_class].
    --clean=<clean>                  Clean data based on certain meta column [default: None].
    --subsample-ref=<subsample-ref>  Subsample reference data [default: None].
    --seed=<seed>                    Set random seed [default: None].
    --self                           Self-projection.
""")

function loadh5(datafile, clean=nothing, subsample=nothing, annotation=nothing, genes=nothing)
    HDF5.h5open(datafile) do file
        counts = read(file, "/exprs")
        if isa(counts, Dict)
            counts = SparseMatrixCSC(
                counts["shape"][2], counts["shape"][1],
                counts["indptr"] .+ 1, counts["indices"] .+ 1,
                counts["data"]
            )
        else
            counts = Matrix(counts)
        end
        featurenames = read(file, "/var_names")
        cellnames = read(file, "/obs_names")
        if annotation != nothing
            annotations = read(file, "/obs/$(annotation)")
        else
            annotations = nothing
        end
        if clean != nothing
            clean = read(file, "/obs/$(clean)")
            mask = clean .∈ [["" "na" "NA" "nan" "NaN"]]
            if count(mask) > 0
                print("Cleaning removed $(count(mask)) cells. ")
                counts = counts[:, .!mask]
                cellnames = cellnames[.!mask]
                if annotations != nothing
                    annotations = annotations[.!mask]
                end
            end
        end
        if subsample != nothing
            subsample = sample(1:size(counts, 2), subsample, replace=false)
            counts = counts[:, subsample]
            cellnames = cellnames[subsample]
            if annotations != nothing
                annotations = annotations[subsample]
            end
        end
        if genes != nothing
            genes = read(file, "/uns/$(genes)")
        else
            genes = nothing
        end
        return (counts=counts, featurenames=featurenames, cellnames=cellnames, annotations=annotations, genes=genes)
    end
end

macro message(msg, block)
    quote
        print(stderr, "  ", rpad(string($(esc(msg)), " "), 25, '\u2015'), " ")
        starttime = now()
        $(esc(block))
        println(stderr, Dates.canonicalize(Dates.CompoundPeriod(now() - starttime)))
    end
end

function writeresult(file, predictions, neighbors, time_per_cell)
    HDF5.h5open(file, "w") do file
        g = g_create(file, "prediction")
        for item in predictions
            k, v = item
            g[String(k)] = v
        end
        file["indexes"] = neighbors.indexes
        file["hammingdistances"] = neighbors.hammingdistances
        file["time"] = time_per_cell
    end
end

using Random: seed!
if args["--seed"] != "None"
    seed!(parse(Int, args["--seed"]))
end

if args["--clean"] == "None"
    args["--clean"] = nothing
end
if args["--subsample-ref"] == "None"
    args["--subsample-ref"] = nothing
else
    args["--subsample-ref"] = parse(Int, args["--subsample-ref"])
end
if args["--genes"] == "cf_genes"
    args["--genes"] = nothing
end

@message "Loading data" begin
    ref = loadh5(args["<reference>"], args["--clean"], args["--subsample-ref"], args["--annotation"], args["--genes"])
end
@message "Selecting features" begin
    if ref.genes == nothing
        features = CellFishing.selectfeatures(ref.counts, ref.featurenames)
    else
        nzv_features = ref.featurenames[dropdims(std(ref.counts, dims=2) .> 0, dims=2)]
        println("Removing $(size(setdiff(ref.genes, nzv_features), 1)) genes with zero variance.")
        ref_genes = intersect(nzv_features, ref.genes)
        features = CellFishing.Features(ref.featurenames, ref.featurenames .∈ [ref_genes])
    end
end
@message "Creating a database" begin
    n_bits = parse(Int, args["--n-bits"])
    n_lshashes = parse(Int, args["--n-lshashes"])
    database = CellFishing.CellIndex(ref.counts, features, metadata=ref.annotations, n_bits=n_bits, n_lshashes=n_lshashes)
end

k = parse(Float64, args["--n-neighbors"])
if k < 1
    k = round(k * size(ref.counts, 2))
    if args["--self"]
        k += 1
    end
end
k = Int64(k)

@message "Loading query data" begin
    query = loadh5(args["<query>"], args["--clean"])
end

@message "Searching the database" begin
    inferstats = args["--infer-stats"]
    inferstats = inferstats == "query"    ? :query :
                    inferstats == "database" ? :database :
                    inferstats == "both"     ? :both :
                    error("invalid --infer-stats option")
    start_time = now()
    neighbors = CellFishing.findneighbors(k, query.counts, query.featurenames, database, inferstats=inferstats)
end

time_per_cell = nothing
@message "Doing annotation" begin
    predictions = Dict()
    majority = parse(Float64, args["--majority"])
    min_hits = parse(Int, args["--min-hits"])
    cutoffs = split(args["--cutoff"], ",")
    for cutoff in cutoffs
        global time_per_cell
        predictions[cutoff] = Vector{String}()
        _cutoff = parse(Int, cutoff)
        for qid in 1:size(neighbors.indexes, 2)
            mask = neighbors.hammingdistances[:, qid] .≤ _cutoff
            annotations = database.metadata[neighbors.indexes[mask, qid]]
            if size(annotations, 1) < min_hits
                push!(predictions[cutoff], "rejected")
                continue
            end
            max_k, max_v = nothing, 0
            for item in countmap(annotations)
                k, v = item
                if v > max_v
                    max_k, max_v = k, v
                end
            end
            if max_v ≤ majority * size(annotations, 1)
                push!(predictions[cutoff], "ambiguous")
            else
                push!(predictions[cutoff], string(max_k))
            end
        end
        if time_per_cell == nothing
            time_per_cell = (now() - start_time).value / size(query.counts, 2)
        end
    end
end

@message "Saving results" begin
    writeresult(args["<output>"], predictions, neighbors, time_per_cell)
end


text:
# Threaded version of TaylorSeries.evaluate!
function evaluate_threads!(x::Array{Taylor1{T},1}, δt::T,
        x0::Union{Array{T,1},SubArray{T,1}}) where {T<:Number}
    # @assert length(x) == length(x0)
    Threads.@threads for i in eachindex(x, x0)
        x0[i] = evaluate( x[i], δt )
    end
    nothing
end

#### Threaded version of TaylorIntegration.stepsize
function stepsize_threads(q::AbstractArray{Taylor1{U},1}, epsilon::T) where
        {T<:Real, U<:Number}
    R = promote_type(typeof(norm(constant_term(q[1]), Inf)), T)
    h = convert(R, Inf)
    #= Threads.@threads =# for i in eachindex(q)
        @inbounds hi = TaylorIntegration.stepsize( q[i], epsilon )
        h = min( h, hi )
    end
    # If `isinf(h)==true`, we use the maximum (finite)
    # step-size obtained from all coefficients as above.
    # Note that the time step is independent from `epsilon`.
    if isinf(h)
        h = zero(R)
        #= Threads.@threads =# for i in eachindex(q)
            @inbounds hi = TaylorIntegration._second_stepsize(q[i], epsilon)
            h = max( h, hi )
        end
    end
    return h::R
end

# first step-size control from Jorba & Zou, 2005
function stepsize_jz05(q::AbstractArray{Taylor1{U}, N}, epsilon::T) where
        {T<:Real, U<:Number, N}
    nbodies = (length(q)-13)÷6
    q0_norminf = norm(constant_term.(q[1:6nbodies]), Inf)
    pred = epsilon*q0_norminf ≤ epsilon

    if pred
        p_jz05 = Int(ceil(-0.5log(epsilon)+1)) # atol
    else
        p_jz05 = Int(ceil(-0.5log(epsilon)+1)) # rtol
    end

    order = min(p_jz05, q[1].order) # q[1].order
    ordm1 = order-1
    invorder = 1/order
    invordm1 = 1/ordm1
    qordm1_norminf = norm(getcoeff.(q[1:6nbodies], ordm1), Inf)
    qorder_norminf = norm(getcoeff.(q[1:6nbodies], order), Inf)

    if pred
        ρ_ordm1 = ( 1/qordm1_norminf )^invordm1
        ρ_order = ( 1/qorder_norminf )^invorder
    else
        ρ_ordm1 = ( q0_norminf/qordm1_norminf )^invordm1
        ρ_order = ( q0_norminf/qorder_norminf )^invorder
    end
    ρ = min(ρ_ordm1, ρ_order)
    return ρ*exp(-2.0)
end

# ##### Constant timestep method: set timestep equal to 1 day
# function stepsize_threads(q::AbstractArray{Taylor1{U},1}, epsilon::T) where
#         {T<:Real, U<:Number}
#     R = promote_type(typeof(norm(constant_term(q[1]), Inf)), T)
#     h = convert(R, Inf)
#     h = TaylorIntegration.stepsize( q[1], epsilon )
#     return one(h)::R
# end

function taylorstep_threads!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
        dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params,
        parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Compute the Taylor coefficients
    TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params)
    # @time TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params)
    # @time TaylorIntegration.__jetcoeffs!(Val(false), f!, t, x, dx, xaux, params)
    # @time TaylorIntegration.__jetcoeffs!(Val(true), f!, t, x, dx, xaux, params)

    # Compute the step-size of the integration using `abstol`
    # δt = stepsize_threads(x, abstol)
    δt = stepsize_jz05(x, abstol)

    return δt
end

function taylorinteg_threads(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T,
        params = nothing; maxsteps::Int=500, parse_eqs::Bool=true, dense::Bool=false) where {T<:Real, U<:Number}

    # Allocation
    tv = Array{T}(undef, maxsteps+1)
    dof = length(q0)
    xv = Array{U}(undef, dof, maxsteps+1)
    if dense
        xv_interp = Array{Taylor1{U}}(undef, dof, maxsteps+1)
    end

    # Initialize the vector of Taylor1 expansions
    t = Taylor1(T, order)
    x = Array{Taylor1{U}}(undef, dof)
    dx = Array{Taylor1{U}}(undef, dof)
    xaux = Array{Taylor1{U}}(undef, dof)
    dx .= Taylor1.(zeros(U), order)

    # Initial conditions
    @inbounds t[0] = t0
    x .= Taylor1.(q0, order)
    x0 = deepcopy(q0)
    @inbounds tv[1] = t0
    @inbounds xv[:,1] .= q0
    sign_tstep = copysign(1, tmax-t0)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs = TaylorIntegration._determine_parsing!(parse_eqs, f!, t, x, dx, params)

    @show parse_eqs

    # Integration
    nsteps = 1
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep_threads!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive!
        # δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs) # δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        evaluate_threads!(x, δt, x0) # new initial condition
        # evaluate!(x, δt, x0) # new initial condition
        t0 += δt
        @inbounds t[0] = t0
        nsteps += 1
        @inbounds tv[nsteps] = t0
        if dense
            # @inbounds xv_interp[:,nsteps-1] .= deepcopy(x)
            Threads.@threads for i in eachindex(x0)
                @inbounds xv_interp[i,nsteps-1] = deepcopy(x[i])
            end
        else
            # @inbounds xv[:,nsteps] .= x0
            Threads.@threads for i in eachindex(x0)
                @inbounds xv[i,nsteps] = x0[i]
            end
        end
        Threads.@threads for i in eachindex(x0)
            @inbounds x[i][0] = x0[i]
            @inbounds dx[i] = Taylor1( zero(x0[i]), order )
        end
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    if dense
        return TaylorInterpolant(tv[1], view(tv.-tv[1],1:nsteps), view(transpose(view(xv_interp,:,1:nsteps-1)),1:nsteps-1,:))
    else
        return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:)
    end
end