text:
#edge pusing algorithm for Hessian reverse AD

function reset_hess{I,V}(tape::Tape{I,V})
    for i = 1:length(tape.bhi)
        @inbounds tape.bhi[i] = Vector{Int}()
        @inbounds tape.bhv[i] = Vector{Float64}()
    end
    fill!(tape.bh_idxes,zero(I))
    tape.nzh = -one(I)     #hess indicator
    nothing
end

@inline function push_edge{I}(bhi::Vector{Vector{I}},to::I,from::I, nvar::I)
    # @show "push_edge - ",to," <--- ", from
    @inbounds push!(bhi[to],from)
    if to <= nvar && from > nvar
        push_edge(bhi,from,to,nvar)
    end
end

function hess_struct{I}(ttstart::I, ttend::I, tt::Vector{I}, trend::I, tr::Vector{I}, nvar::I, bhi::Vector{Vector{I}})
    #@timing tape.enable_timing_stats tic()
    idx = ttend
    trlen = trend - 1

    @inbounds while idx > ttstart
        @inbounds ntype = tt[idx]
        # @show ntype, trlen, idx
        if ntype == TYPE_P
            #@assert false
        elseif ntype == TYPE_V
            idx -= 3
        elseif ntype == TYPE_O
            @inbounds n = tt[idx-1]
            @inbounds oc = tt[idx-2]

            #pushing
            @inbounds i_id = tt[idx-3]
            @inbounds lvi = bhi[i_id]
            for j = 1:length(lvi)
                # @inbounds p_id = lvi[j].i
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    for j0=trlen-n+1:trlen
                        @inbounds ci_id = tr[j0]
                        push_edge(bhi,ci_id, ci_id, nvar)
                        for j1=j0+1:trlen
                            @inbounds cii_id = tr[j1]
                            push_edge(bhi,cii_id, ci_id, nvar)
                        end
                    end
                else  #when i_id != p_id
                    for j0 = trlen -n + 1:trlen
                        @inbounds ci_id = tr[j0]
                        push_edge(bhi,ci_id, p_id, nvar)
                    end
                end
            end

            #creating
            @inbounds op_sym = OP[oc]
            # @show op_sym
            if (op_sym ==:+ || op_sym == :-)
                #zeros
            elseif (n == 1 && op_sym != :-) #1-ary operator
                @inbounds ci_id = tr[trlen]
                push_edge(bhi,ci_id, ci_id, nvar)
            elseif (op_sym == :*)
                # @show "times ", n
                for j0 = trlen -n + 1:trlen
                    @inbounds ci_id = tr[j0]
                    for j1=j0+1:trlen
                        @inbounds cii_id = tr[j1]
                        push_edge(bhi,cii_id, ci_id, nvar)
                    end
                end
            elseif op_sym == :/ # binary operator /
                assert(n==2)
                @inbounds ri_id = tr[trlen]
                @inbounds li_id = tr[trlen-1]
                push_edge(bhi,ri_id,li_id, nvar)
                push_edge(bhi,ri_id,ri_id, nvar)
            else # other binary
                assert(n==2)
                @inbounds ri_id = tr[trlen]
                @inbounds li_id = tr[trlen-1]
                push_edge(bhi,li_id, li_id, nvar)
                push_edge(bhi,ri_id, li_id, nvar)
                push_edge(bhi,ri_id, ri_id, nvar)
            end
            trlen -= n
            idx -= 5
        elseif ntype == TYPE_O1
            @inbounds oc = tt[idx-2]
            #pushing
            @inbounds i_id = tt[idx-3]
            @inbounds lvi = bhi[i_id]
            @inbounds c_id = tr[trlen]

            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi,c_id, c_id, nvar)
                else
                    push_edge(bhi,c_id,p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :/ || op_sym == :^
                push_edge(bhi,c_id,c_id, nvar)
            end
            trlen -= 1
            idx -= 5
        elseif ntype == TYPE_O2
            @inbounds oc = tt[idx-2]
            #pushing
            @inbounds i_id = tt[idx-3]
            @inbounds lvi = bhi[i_id]
            @inbounds c_id = tr[trlen]

            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi,c_id, c_id, nvar)
                else
                    push_edge(bhi,c_id,p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :^
                push_edge(bhi,c_id,c_id, nvar)
            end
            trlen -= 1
            idx -= 5
        elseif ntype == TYPE_O3
            @inbounds oc = tt[idx-3]
            @inbounds i_id = tt[idx-4]
            @inbounds vidx = tt[idx-1]

            @inbounds lvi = bhi[i_id]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi,vidx, vidx, nvar)
                else
                    push_edge(bhi,vidx,p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :/ || op_sym == :^
                push_edge(bhi,vidx,vidx, nvar)
            end
            idx -= 6
        elseif ntype == TYPE_O4
            @inbounds oc = tt[idx-3]
            @inbounds i_id = tt[idx-4]
            @inbounds vidx = tt[idx-1]

            @inbounds lvi = bhi[i_id]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi,vidx, vidx, nvar)
                else
                    push_edge(bhi,vidx,p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :^
                push_edge(bhi,vidx,vidx, nvar)
            end
            idx -= 6
        elseif ntype == TYPE_O5
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]
            @inbounds vidx = tt[idx-1]

            @inbounds lvi = bhi[i_id]
            @inbounds c_id = tr[trlen]
            #pusing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi, vidx, vidx, nvar)
                    push_edge(bhi, c_id, vidx, nvar)
                    push_edge(bhi, c_id, c_id, nvar)
                else
                    push_edge(bhi, vidx, p_id, nvar)
                    push_edge(bhi, c_id, p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym ==:*
                push_edge(bhi,c_id,vidx, nvar)
            elseif op_sym ==:/
                push_edge(bhi,c_id,vidx, nvar)
                push_edge(bhi,c_id,c_id, nvar)
            elseif op_sym ==:^
                push_edge(bhi,vidx,vidx, nvar)
                push_edge(bhi,c_id,vidx, nvar)
                push_edge(bhi,c_id,c_id, nvar)
            else
                #@assert op_sym == :+ || op_sym ==:-
            end
            trlen -= 1
            idx -= 5
        elseif ntype == TYPE_O6
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]
            @inbounds vidx = tt[idx-1]

            @inbounds lvi = bhi[i_id]
            @inbounds c_id = tr[trlen]
            #pusing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi, c_id, c_id, nvar)
                    push_edge(bhi, vidx, c_id, nvar)
                    push_edge(bhi, vidx, vidx, nvar)
                else
                    push_edge(bhi, c_id, p_id, nvar)
                    push_edge(bhi, vidx, p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym ==:*
                push_edge(bhi,vidx,c_id, nvar)
            elseif op_sym ==:/
                push_edge(bhi,vidx,c_id, nvar)
                push_edge(bhi,vidx,vidx, nvar)
            elseif op_sym ==:^
                push_edge(bhi,c_id,c_id, nvar)
                push_edge(bhi,vidx,c_id, nvar)
                push_edge(bhi,vidx,vidx, nvar)
            else
                #@assert op_sym == :+ || op_sym ==:-
            end
            trlen -= 1
            idx -= 5

        elseif ntype == TYPE_O7
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]
            @inbounds vidx = tt[idx-1]

            @inbounds lvi = bhi[i_id]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                if p_id == i_id
                    push_edge(bhi,vidx,vidx, nvar)
                else
                    push_edge(bhi,vidx,p_id, nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym != :-
                push_edge(bhi,vidx,vidx, nvar)
            end
            idx -= 5
        else
            #@assert false
        end
    end #end while loop

    @assert trlen == 0 && idx == (ttstart-1)

    #@timing tape.enable_timing_stats tape.ep_structure += toq()
    nothing
end

@inline function hess_findnz{I,V}(bhi::Vector{Vector{I}}, idxmap::Vector{I}, HH::Vector{V}, nvar::I)
    #@timing tape.enable_timing_stats tic()
    nnz = zero(I)
    for i = 1:nvar
        @inbounds lvi = bhi[i]
        for j = 1:length(lvi)
            @inbounds v_id = lvi[j]
            if v_id <= nvar
                nnz += 1
            end
        end
    end
    resize!(idxmap, nnz)
    fill!(idxmap,0)
    resize!(HH, nnz)
    #@timing tape.enable_timing_stats tape.ep_findnz += toq()
    return nnz
end

function hess_struct_recovery{I}(bhi::Vector{Vector{I}}, idxmap::Vector{I}, nvar::I, nnz_original::I, h_I::Vector{I}, h_J::Vector{I})
    #@timing tape.enable_timing_stats tic()
    @assert length(h_I) == length(h_J)
    @assert length(idxmap) == nnz_original

    hh_I = zeros(I,nnz_original)
    hh_J = zeros(I,nnz_original)
    start = one(I)
    for i=1:nvar
        @inbounds lvi = bhi[i]
        for j = 1:length(lvi)
            @inbounds v_id = lvi[j]
            if v_id <= nvar
                if v_id <= i  #assert lower trangular
                    @inbounds hh_I[start] = i
                    @inbounds hh_J[start] = v_id
                else
                    @inbounds hh_I[start] = v_id
                    @inbounds hh_J[start] = i
                end
                start += 1
            end
        end
    end
    @assert (start - 1) == nnz_original "$start $(nnz_original)"

    #postprocess non-repeated
    mergedindices = zeros(I,length(hh_I))
    @inbounds mergednnz = [0]
    mergedmap = zeros(I,length(hh_I))
    function combine(idx1,idx2)
        @inbounds if mergedmap[idx1] == 0 && mergedmap[idx2] != 0
            @inbounds mergednnz[1] += 1
            @inbounds mergedmap[idx1] = idx2
            @inbounds mergedindices[mergednnz[1]] = idx1
            return idx2
        else
            @inbounds @assert mergedmap[idx2] == 0
            @inbounds mergednnz[1] += 1
            @inbounds mergedmap[idx2] = idx1
            @inbounds mergedindices[mergednnz[1]] = idx2
            return idx1
        end
    end
    Hmat = sparse(hh_I, hh_J, [i for i in 1:length(hh_I)], nvar, nvar, combine)
    nnz_non_repeat = length(Hmat.nzval)
    pre = length(h_I)
    start = pre + 1
    append!(h_I, zeros(nnz_non_repeat))
    append!(h_J, zeros(nnz_non_repeat))

    for row in 1:nvar
        @inbounds for pos in Hmat.colptr[row]:(Hmat.colptr[row+1]-1)
            @inbounds col = Hmat.rowval[pos]
            @inbounds origidx = Hmat.nzval[pos] # this is the original index of this element
            @inbounds idxmap[origidx] = pos
            @inbounds h_I[start] = row
            @inbounds h_J[start] = col
            start += 1
        end
    end

    @inbounds for k in 1:mergednnz[1]
        @inbounds origidx = mergedindices[k]
        @inbounds mergedwith = mergedmap[origidx]
        @inbounds @assert idxmap[origidx] == 0
        @inbounds @assert idxmap[mergedwith] != 0
        @inbounds idxmap[origidx] = idxmap[mergedwith]
    end

    for i in 1:nnz_original
         @inbounds @assert idxmap[i] != 0
    end

    # @assert start - pre  - 1 == nz "$start $pre $(nz)"
    #@timing tape.enable_timing_stats tape.ep_structure_recovery += toq()

    nothing
end

function forward_pass_2ord{I,V}(ttstart::I, ttend::I, tt::Vector{I}, vvals::Array{V,1}, pvals::Array{V,1}, stk::Vector{V}, vals::Vector{V})
    #@timing tape.enable_timing_stats tic()

    idx = ttstart
    vallen = zero(I)
    stklen = zero(I)

    @inbounds while idx <= ttend
        ntype = tt[idx]
        # @show ntype
        # if ntype == TYPE_P
        #     #@assert false
        # else
        if ntype == TYPE_V
            stklen += 1
            @inbounds stk[stklen] = vvals[tt[idx+1]]
            idx += 3

        elseif ntype == TYPE_O
            @inbounds oc = tt[idx+2]
            @inbounds n = tt[idx+3]

            if n == 1 # 1-argument functions
                @inbounds val = eval_0ord(OP[oc],stk[stklen])
                vallen += 1
                @inbounds vals[vallen] = stk[stklen]
                @inbounds stk[stklen] = val
            elseif n == 2
                # @show OP[oc],stk
                @inbounds val = eval_0ord(OP[oc],stk[stklen-1],stk[stklen])
                @inbounds vals[vallen+1] = stk[stklen-1]
                @inbounds vals[vallen+2] = stk[stklen]
                vallen += 2
                stklen -= 1
                @inbounds stk[stklen] = val
            else
                @inbounds val = eval_0ord(OP[oc],stk,stklen-n+1,stklen)
                @simd for i=1:n
                    @inbounds vals[vallen+i] = stk[stklen-n+i]
                end
                vallen += n
                stklen -= n-1
                @inbounds stk[stklen] = val
            end
            idx += 5

        elseif ntype == TYPE_O1
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds val = eval_0ord(OP[oc],pval,stk[stklen])
            vallen += 1
            @inbounds vals[vallen] = stk[stklen]
            @inbounds stk[stklen] = val
            idx += 5

        elseif ntype == TYPE_O2
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds val = eval_0ord(OP[oc],stk[stklen],pval)
            vallen += 1
            vals[vallen] = stk[stklen]
            @inbounds stk[stklen] = val
            idx += 5

        elseif ntype == TYPE_O3
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds vval = vvals[tt[idx+4]]
            stklen += 1
            @inbounds stk[stklen] = eval_0ord(OP[oc],pval,vval)
            idx += 6

        elseif ntype == TYPE_O4
            @inbounds oc = tt[idx+2]
            @inbounds pval = pvals[tt[idx+3]]
            @inbounds vval = vvals[tt[idx+4]]
            stklen += 1
            @inbounds stk[stklen] = eval_0ord(OP[oc],vval,pval)
            idx += 6

        elseif ntype == TYPE_O5
            @inbounds oc = tt[idx+2]
            @inbounds vval = vvals[tt[idx+3]]
            @inbounds val = eval_0ord(OP[oc],vval,stk[stklen])
            vallen += 1
            @inbounds vals[vallen] = stk[stklen]
            @inbounds stk[stklen] = val
            idx += 5

        elseif ntype == TYPE_O6
            @inbounds oc = tt[idx+2]
            @inbounds vval = vvals[tt[idx+3]]
            @inbounds val = eval_0ord(OP[oc],stk[stklen],vval)
            vallen += 1
            @inbounds vals[vallen] = stk[stklen]
            @inbounds stk[stklen] = val
            idx += 5

        elseif ntype == TYPE_O7
            @inbounds oc = tt[idx+2]
            @inbounds vval = vvals[tt[idx+3]]
            stklen += 1
            @inbounds stk[stklen] = eval_0ord(OP[oc],vval)
            idx += 5

        # else
        #     #@assert false
        end
    end
    @assert stklen == 1 && idx == (ttend + 1)
    vallen += 1
    @inbounds vals[vallen] = stk[stklen]
    #@timing tape.enable_timing_stats tape.ep_forward_time += toq()
    return vallen
end

@inline function update(bhv,to,from,w,bh_idxes,nvar)
    # @show "update - ", to, "<-- ", from, w
    @inbounds bh_idxes[to] += 1
    # @assert tape.bhi[to][tape.bh_idxes[to]] == from
    @inbounds bhv[to][bh_idxes[to]] = w

    if to <= nvar && from > nvar
        update(bhv,from,to,w,bh_idxes,nvar)
    end
end

function reverse_pass_2ord{I,V}(ttstart::I, ttend::I, tt::Vector{I}, trend::I, tr::Vector{I},
    bhi::Vector{Vector{I}}, bhv::Vector{Vector{V}}, bh_idxes::Vector{I}, nvar::I,
    vvals::Vector{V}, pvals::Vector{V}, vallen::I, vals::Vector{V}, stk::Vector{V}, imm::Vector{V})

    idx = ttend
    trlen = trend - 1
    vallen -= 1

    adjs = stk
    adjlen = 1
    @inbounds adjs[adjlen] = one(V)  #initialize adjoint = 1

    while idx > ttstart
        #@timing tape.enable_timing_stats tic()   #timing
        #@timing tape.enable_timing_stats node_id = 0  #timing

        @inbounds ntype = tt[idx]
        @inbounds adj = adjs[adjlen]

        # @show ntype,trlen,idx

        # if ntype == TYPE_P
        #     #@assert false
        # else
        if ntype == TYPE_V
            #@timing tape.enable_timing_stats node_id = tt[idx-1]  #timing
            adjlen -= 1
            idx -= 3 #skip TYPE_V

        elseif ntype == TYPE_O
            @inbounds n = tt[idx-1]
            @inbounds oc = tt[idx-2]

            @inbounds i_id = tt[idx-3]
            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds op_sym = OP[oc]

            if n == 1
                @inbounds eval_2ord(OP[oc],vals[vallen],imm)
                @inbounds ci_id = tr[trlen]
                #pushing
                for j = 1:length(lvi)
                    @inbounds p_id = lvi[j]
                    @inbounds w = lvv[j]
                    if p_id == i_id
                        @inbounds w_bar = imm[1]*imm[1]*w
                        update(bhv,ci_id, ci_id, w_bar,bh_idxes,nvar)
                    else
                        @inbounds w_bar = imm[1]*w
                        update(bhv,ci_id, p_id, ci_id == p_id? 2.0*w_bar:w_bar,bh_idxes,nvar)
                    end
                end
                #creating
                if op_sym!=:-
                    @inbounds w_bar = adj*imm[2]
                    update(bhv,ci_id, ci_id, w_bar,bh_idxes,nvar)
                end
            elseif n == 2
                @inbounds eval_2ord(OP[oc],vals[vallen-1],vals[vallen],imm)
                @inbounds dl = imm[1]
                @inbounds dr = imm[2]
                @inbounds li_id = tr[trlen-1]
                @inbounds ri_id = tr[trlen]

                #pushing
                for j = 1:length(lvi)
                    @inbounds p_id = lvi[j]
                    @inbounds w = lvv[j]
                    dlw = dl*w
                    drw = dr*w
                    if p_id == i_id
                        update(bhv,li_id, li_id, dl*dlw,bh_idxes,nvar)
                        dlrw = dl*drw
                        update(bhv,ri_id, li_id, ri_id == li_id? 2.0*dlrw:dlrw,bh_idxes,nvar)
                        update(bhv,ri_id, ri_id, dr*drw,bh_idxes,nvar)
                    else
                        update(bhv,li_id, p_id, li_id == p_id? 2.0*dlw:dlw,bh_idxes,nvar)
                        update(bhv,ri_id, p_id, ri_id == p_id? 2.0*drw:drw,bh_idxes,nvar)
                    end
                end

                #creating
                if op_sym == :+ || op_sym==:-
                    #nothing
                elseif op_sym == :/
                    #@assert n==2
                    @inbounds dlr = imm[4]
                    @inbounds drr = imm[5]
                    w_bar = adj*dlr
                    update(bhv,ri_id,li_id, ri_id == li_id?2.0*w_bar:w_bar,bh_idxes,nvar)
                    update(bhv,ri_id,ri_id,adj*drr,bh_idxes,nvar)

                elseif op_sym == :*
                    #@assert n==2
                    @inbounds dlr = imm[4]
                    update(bhv,ri_id, li_id,ri_id == li_id? 2.0*adj:adj,bh_idxes,nvar)

                else
                    @assert op_sym == :^
                    @inbounds dll = imm[3]
                    @inbounds dlr = imm[4]
                    @inbounds drr = imm[5]
                    update(bhv,li_id,li_id,adj*dll,bh_idxes,nvar)
                    w_bar = adj*dlr
                    update(bhv,ri_id,li_id,ri_id==li_id?2.0*w_bar:w_bar,bh_idxes,nvar)
                    update(bhv,ri_id,ri_id,adj*drr,bh_idxes,nvar)

                end
            else
                @assert n>2
                @inbounds eval_2ord(OP[oc],vals,vallen-n+1,vallen,imm)
                #pushing
                if op_sym == :+
                    for j = 1:length(lvi)
                        @inbounds p_id = lvi[j]
                        @inbounds w = lvv[j]
                        if p_id == i_id
                            for j0=trlen-n+1:trlen
                                @inbounds ci_id = tr[j0]
                                w_bar = w
                                update(bhv, ci_id, ci_id, w_bar,bh_idxes,nvar)
                                for j1=j0+1:trlen
                                    @inbounds cii_id = tr[j1]
                                    update(bhv,cii_id, ci_id,w_bar,bh_idxes,nvar)
                                end
                            end
                        else
                            for j0 = trlen -n + 1:trlen
                                @inbounds ci_id = tr[j0]
                                w_bar = w
                                update(bhv,ci_id, p_id, ci_id == p_id? 2.0*w_bar:w_bar,bh_idxes,nvar)
                            end
                        end
                    end
                    #creating
                        #nothing
                else
                    @assert op_sym == :*
                    for j = 1:length(lvi)
                        @inbounds p_id = lvi[j]
                        @inbounds w = lvv[j]
                        if p_id == i_id
                            k=1
                            for j0=trlen-n+1:trlen
                                @inbounds ci_id = tr[j0]
                                @inbounds t0 = imm[k]
                                t1 = t0 * w
                                w_bar = t0 * t1
                                update(bhv, ci_id, ci_id, w_bar,bh_idxes,nvar)
                                k0 = k + 1
                                for j1=j0+1:trlen
                                    @inbounds cii_id = tr[j1]
                                    @inbounds w_bar1 = t1*imm[k0]
                                    update(bhv,cii_id, ci_id,cii_id==ci_id? 2.0*w_bar1:w_bar1,bh_idxes,nvar)
                                    k0 += 1
                                end
                                k += 1
                            end
                        else
                            k=1
                            for j0 = trlen -n + 1:trlen
                                @inbounds ci_id = tr[j0]
                                @inbounds w_bar = imm[k]*w
                                update(bhv,ci_id, p_id, ci_id == p_id? 2.0*w_bar:w_bar,bh_idxes,nvar)
                                k += 1
                            end
                        end
                    end
                    #creating
                    # @show "creating", adj
                    k=n+1
                    @simd for j0 = trlen -n + 1:trlen
                        @inbounds ci_id = tr[j0]
                        for j1=j0+1:trlen
                            @inbounds cii_id = tr[j1]
                            @inbounds w_bar = adj*imm[k]
                            update(bhv,cii_id, ci_id, cii_id == ci_id? 2.0*w_bar:w_bar,bh_idxes,nvar)
                            k += 1
                        end
                    end
                end
            end

            #updating
            if op_sym == :+
                @simd for i=1:n
                    @inbounds adjs[adjlen+i-1] = adj
                end
            else
                @simd for i=1:n
                    @inbounds adjs[adjlen+i-1] = imm[i]*adj
                end
            end
            adjlen += n-1
            trlen -= n
            vallen -= n
            idx -= 5
        elseif ntype == TYPE_O1
            @inbounds pval = pvals[tt[idx-1]]
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]
            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds eval_2ord(OP[oc],pval,vals[vallen],imm)
            @inbounds c_id = tr[trlen]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds t0 = imm[2]*w
                if p_id == i_id
                    @inbounds w_bar = imm[2]*t0
                    update(bhv,c_id,c_id,w_bar,bh_idxes,nvar)
                else
                    w_bar = t0
                    update(bhv,c_id,p_id,c_id==p_id?2.0*w_bar:w_bar,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :/ || op_sym == :^
                @inbounds update(bhv,c_id,c_id,adj*imm[5],bh_idxes,nvar)
            end

            #updating
            @inbounds adjs[adjlen] = imm[2]*adj
            trlen -= 1
            vallen -= 1
            idx -= 5

        elseif ntype == TYPE_O2
            @inbounds pval = pvals[tt[idx-1]]
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]
            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds eval_2ord(OP[oc],vals[vallen],pval,imm)
            @inbounds c_id = tr[trlen]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds t0 = imm[1]*w
                if p_id == i_id
                    @inbounds w_bar = imm[1]*t0
                    update(bhv,c_id,c_id,w_bar,bh_idxes,nvar)
                else
                    w_bar = t0
                    update(bhv,c_id,p_id,c_id==p_id?2.0*w_bar:w_bar,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :^
                @inbounds update(bhv,c_id,c_id,adj*imm[3],bh_idxes,nvar)
            end

            #updating
            @inbounds adjs[adjlen] = imm[1]*adj
            trlen -= 1
            vallen -= 1
            idx -= 5

        elseif ntype == TYPE_O3
            @inbounds vidx = tt[idx-1]
            @inbounds vval = vvals[vidx]
            @inbounds pval = pvals[tt[idx-2]]
            @inbounds oc = tt[idx-3]
            @inbounds i_id = tt[idx-4]
            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds eval_2ord(OP[oc],pval,vval,imm)
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds t0 = imm[2]*w
                if p_id == i_id
                    @inbounds w_bar = imm[2]*t0
                    update(bhv,vidx,vidx,w_bar,bh_idxes,nvar)
                else
                    w_bar = t0
                    update(bhv,vidx,p_id,vidx==p_id?2.0*w_bar:w_bar,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :/ || op_sym == :^
                @inbounds update(bhv,vidx,vidx,adj*imm[5],bh_idxes,nvar)
            end

            #updating
            adjlen -= 1
            idx -= 6

        elseif ntype == TYPE_O4
            @inbounds vidx = tt[idx-1]
            @inbounds vval = vvals[vidx]
            @inbounds pval = pvals[tt[idx-2]]
            @inbounds oc = tt[idx-3]
            @inbounds i_id = tt[idx-4]
            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds eval_2ord(OP[oc],vval,pval,imm)
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds t0 = imm[1]*w
                if p_id == i_id
                    @inbounds w_bar = imm[1]*t0
                    update(bhv,vidx,vidx,w_bar,bh_idxes,nvar)
                else
                    w_bar = t0
                    update(bhv,vidx,p_id,vidx==p_id?2.0*w_bar:w_bar,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym == :^
                @inbounds update(bhv,vidx,vidx,adj*imm[3],bh_idxes,nvar)
            end

            #updating
            adjlen -= 1
            idx -= 6

        elseif ntype == TYPE_O5
            @inbounds vidx = tt[idx-1]
            @inbounds vval = vvals[vidx]
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]

            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds eval_2ord(OP[oc],vval,vals[vallen],imm)
            @inbounds c_id = tr[trlen]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds dlw = imm[1]*w
                @inbounds drw = imm[2]*w
                if p_id == i_id
                    @inbounds update(bhv,vidx,vidx,imm[1]*dlw,bh_idxes,nvar)
                    @inbounds update(bhv,c_id,vidx,imm[1]*drw,bh_idxes,nvar)
                    @inbounds update(bhv,c_id,c_id,imm[2]*drw,bh_idxes,nvar)
                else
                    update(bhv,vidx,p_id,vidx==p_id?2.0*dlw:dlw,bh_idxes,nvar)
                    update(bhv,c_id,p_id,c_id==p_id?2.0*drw:drw,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            @inbounds dll = imm[3]
            @inbounds dlr = imm[4]
            @inbounds drr = imm[5]
            adlr = adj*dlr
            if op_sym == :*
                update(bhv,c_id,vidx,c_id == vidx?2.0*adlr:adlr,bh_idxes,nvar)
            elseif op_sym == :/
                update(bhv,c_id,vidx,c_id == vidx?2.0*adlr:adlr,bh_idxes,nvar)
                update(bhv,c_id,c_id,adj*drr)
            elseif op_sym == :^
                update(bhv,vidx,vidx,adj*dll)
                update(bhv,c_id,vidx,c_id == vidx?2.0*adlr:adlr,bh_idxes,nvar)
                update(bhv,c_id,c_id,adj*drr,bh_idxes,nvar)
            else
                #@assert op_sym == :+ || op_sym ==:-
            end

            #updating
            # @inbounds adjs[adjlen] = imm[1]*adj
            @inbounds adjs[adjlen] = imm[2]*adj
            trlen -= 1
            vallen -= 1
            idx -= 5

        elseif ntype == TYPE_O6
            @inbounds vidx = tt[idx-1]
            @inbounds vval = vvals[vidx]
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]

            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]

            @inbounds eval_2ord(OP[oc],vals[vallen],vval,imm)
            @inbounds c_id = tr[trlen]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds dlw = imm[1]*w
                @inbounds drw = imm[2]*w
                if p_id == i_id
                    @inbounds update(bhv,c_id,c_id,imm[1]*dlw,bh_idxes,nvar)
                    @inbounds update(bhv,vidx,c_id,imm[1]*drw,bh_idxes,nvar)
                    @inbounds update(bhv,vidx,vidx,imm[2]*drw,bh_idxes,nvar)
                else
                    update(bhv,c_id,p_id,c_id==p_id?2.0*dlw:dlw,bh_idxes,nvar)
                    update(bhv,vidx,p_id,vidx==p_id?2.0*drw:drw,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            @inbounds dll = imm[3]
            @inbounds dlr = imm[4]
            @inbounds drr = imm[5]
            adlr = adj*dlr
            if op_sym == :*
                update(bhv,vidx,c_id,vidx == c_id?2.0*adlr:adlr,bh_idxes,nvar)
            elseif op_sym == :/
                update(bhv,vidx,c_id,vidx == c_id?2.0*adlr:adlr,bh_idxes,nvar)
                update(bhv,vidx,vidx,adj*drr,bh_idxes,nvar)
            elseif op_sym == :^
                update(bhv,c_id,c_id,adj*dll,bh_idxes,nvar)
                update(bhv,vidx,c_id,vidx == c_id?2.0*adlr:adlr,bh_idxes,nvar)
                update(bhv,vidx,vidx,adj*drr,bh_idxes,nvar)
            else
                #@assert op_sym == :+ || op_sym ==:-
            end

            #updating
            @inbounds adjs[adjlen] = imm[1]*adj
            trlen -= 1
            vallen -= 1
            idx -= 5

        else
            @assert ntype == TYPE_O7
            @inbounds vidx = tt[idx-1]
            @inbounds oc = tt[idx-2]
            @inbounds i_id = tt[idx-3]
            @inbounds vval = vvals[vidx]

            @inbounds eval_2ord(OP[oc],vval,imm)

            @inbounds lvi = bhi[i_id]
            @inbounds lvv = bhv[i_id]
            #pushing
            for j = 1:length(lvi)
                @inbounds p_id = lvi[j]
                @inbounds w = lvv[j]
                @inbounds t1 = imm[1]*w
                if p_id == i_id
                    @inbounds w_bar = imm[1]*t1
                    update(bhv,vidx,vidx,w_bar,bh_idxes,nvar)
                else
                    update(bhv,vidx,p_id,vidx==p_id?2.0*t1:t1,bh_idxes,nvar)
                end
            end
            #creating
            @inbounds op_sym = OP[oc]
            if op_sym != :-
                @inbounds w_bar = adj*imm[2]
                update(bhv,vidx,vidx,w_bar,bh_idxes,nvar)
            end
            adjlen -= 1
            idx -= 5
        end

        #@timing tape.enable_timing_stats tape.ep_reverse_times[node_id] += toq()
        # println("++++++++++++++++++++++++++++++++++++")
    end  #end while

    # @show idx, vallen, adjlen, trlen
    @assert idx == (ttstart-1) && vallen == 0 && adjlen == 0 && trlen == 0
    nothing
end

@inline function hess_recovery{I,V}(bhi::Vector{Vector{I}},bhv::Vector{Vector{V}}, idxmap::Vector{I}, HH::Vector{V}, nvar::I, start::Int, H::Vector{V})
    #@timing tape.enable_timing_stats tic()
    pre = start
    for i = 1:nvar
        @inbounds lvi = bhi[i]
        @inbounds lvv = bhv[i]
        for j = 1:length(lvi)
            @inbounds v_id = lvi[j]
            if v_id <= nvar
                @inbounds HH[start] = lvv[j]
                start += 1
            end
        end
    end
    nnz =  start - pre
    @assert nnz == length(HH)

    fill!(H,0)
    for i = 1:nnz
        @inbounds H[idxmap[i]] += HH[i]
    end

    #@timing tape.enable_timing_stats tape.ep_recovery_time += toq()
    nothing
end

#Interface function
function hess_structure{I,V}(ts::I, te::I,tt::Vector{I},
    re::I,tr::Vector{I},
    hs::HessStorage{I,V}, h_I::Vector{I}, h_J::Vector{I})

    hess_struct(ts,te,tt,re,tr,hs.gnvar,hs.bhi)
    @assert length(hs.bhi) == length(hs.bhv) == length(hs.bh_idxes)
    for i=1:length(hs.bhi)
        @inbounds resize!(hs.bhv[i],length(hs.bhi[i]))
    end
    nnz_repeat = hess_findnz(hs.bhi, hs.idxmap, hs.HH, hs.gnvar)
    hess_struct_recovery(hs.bhi, hs.idxmap, hs.gnvar, nnz_repeat, h_I, h_J)
    return length(h_I)
end

function hess_reverse{I,V}(ts::I, te::I, tt::Vector{I},
    re::I,tr::Vector{I},
    vvals::Vector{V},pvals::Vector{V}, stk::Vector{V}, vals::Vector{V}, imm::Vector{V}, hs::HessStorage{I,V}, start::I, H::Vector{V})
    #@timing tape.enable_timing_stats tape.ep_n += 1

    vallen = forward_pass_2ord(ts, te, tt, vvals, pvals, stk, vals)
    # @assert vallen == tape.nnode

    fill!(hs.bh_idxes,zero(I))
    reverse_pass_2ord(ts, te, tt, re, tr, hs.bhi, hs.bhv, hs.bh_idxes, hs.gnvar
        ,vvals,pvals, vallen, vals, stk, imm)
    # @time recovery2(tape, tape.bhi,tape.bhv,tape.hess,tape.nvar, factor)
    hess_recovery(hs.bhi, hs.bhv, hs.idxmap, hs.HH, hs.gnvar, start, H)
    nothing
end


text:
export inverselines, invorder_firstdim
"""
    inverselines(nbit::Int; n_reg::Int=nbit) -> ChainBlock

inverse first `n_reg` lines

TODO:
deprecate this function, it is not used.
"""
function inverselines(nbit::Int; n_reg::Int=nbit)
    c = chain(nbit)
    for i = 1:(n_reg ÷ 2)
        push!(c, swap(i,(n_reg-i+1)))
    end
    c
end

"""
    invorder_firstdim(v::VecOrMat) -> VecOrMat

inverse the bit order of first dimension.
"""
function invorder_firstdim(v::Matrix)
    w = similar(v)
    n = size(v, 1) |> log2i
    n_2 = n ÷ 2
    mask = [bmask(i, n-i+1) for i in 1:n_2]
    @simd for b in basis(n)
        @inbounds w[breflect(b, mask; nbits=n)+1,:] = v[b+1,:]
    end
    w
end

function invorder_firstdim(v::Vector)
    n = length(v) |> log2i
    n_2 = n ÷ 2
    w = similar(v)
    #mask = SVector{n_2, Int}([bmask(i, n-i+1)::Int for i in 1:n_2])
    mask = [bmask(i, n-i+1)::Int for i in 1:n_2]
    @simd for b in basis(n)
        @inbounds w[breflect(b, mask; nbits=n)+1] = v[b+1]
    end
    w
end

ext: jl


text:
#= src/snorm.jl

References:

  J. Dixon. Estimating extremal eigenvalues and condition numbers of matrices.
    SIAM J. Numer. Anal. 20 (4): 812-814, 1983.

  J. Kuczynski, H. Wozniakowski. Estimating the largest eigenvalue by the power
    and Lanczos algorithms with a random start. SIAM J. Matrix Anal. Appl. 13
    (4), 1094-1122, 1992.
=#

# spectral norm
  function snorm(A::AbstractLinOp{T}, opts::LRAOptions=LRAOptions(T); args...) where T
    opts = isempty(args) ? opts : copy(opts; args...)
    m, n   = size(A)
    isherm = ishermitian(A)
    xn     = crandn(T, n)
    xm     = Array{T}(undef, m)
    xnrm   = norm(xn)
    s      = one(real(T))
    t      = 0
    niter  = 0
    while s > 0 && abs(s - t) > max(opts.atol, t*opts.rtol)
      if niter == opts.snorm_niter
        opts.verb && warn("iteration limit ($(opts.snorm_niter)) reached")
        break
      end
      niter += 1
      rmul!(xn, 1/xnrm)
      if isherm
        mul!(xm, A, xn)
        copyto!(xn, xm)
      else
         mul!(xm, A, xn)
        mul!(xn, A', xm)
      end
      xnrm = norm(xn)
      t = s
      s = isherm ? xnrm : sqrt(xnrm)
    end
    s
  end

snorm(A, args...; kwargs...) = snorm(LinOp(A), args...; kwargs...)

# spectral norm difference
snormdiff(A::AbstractLinOp{T}, B::AbstractLinOp{T}, args...; kwargs...) where {T} =
  snorm(A - B, args...; kwargs...)
snormdiff(A, B, args...; kwargs...) =
  snormdiff(LinOp(A), LinOp(B), args...; kwargs...)


text:
using CoolProp
using Compat
using Base.Test
const branchname = begin
  if (isdefined(:LibGit2))
    LibGit2.branch(LibGit2.GitRepo(abspath(@__FILE__, "..", "..")));
  else
    Base.Git.branch(dir = abspath(@__FILE__, "..", ".."));
  end
end
info("On $branchname");

include("testThrows.jl");

dl = CoolProp.get_debug_level();
CoolProp.set_debug_level(10);
@test CoolProp.get_debug_level()==10;
CoolProp.set_debug_level(Int(dl));
const coolpropfluids = map(Compat.String, split(get_global_param_string("FluidsList"),','));
const coolpropparameters = map(Compat.String, split(get_global_param_string("parameter_list"),','));
const coolpropmix = map(Compat.String, split(get_global_param_string("predefined_mixtures"), ','));
#all trivials taken from http://www.coolprop.org/coolprop/HighLevelAPI.html#table-of-string-inputs-to-propssi-function
const coolproptrivialparameters = ["ACENTRIC", "DIPOLE_MOMENT", "FH", "FRACTION_MAX", "FRACTION_MIN",
  "GAS_CONSTANT", "GWP100", "GWP20", "GWP500", "HH", "M", "ODP", "PCRIT", "PH", "PMAX", "PMIN", "PTRIPLE",
  "P_REDUCING", "RHOCRIT", "RHOMASS_REDUCING", "RHOMOLAR_CRITICAL", "RHOMOLAR_REDUCING", "TCRIT", "TMAX",
  "TMIN", "TTRIPLE", "T_FREEZE", "T_REDUCING"];
const trivalwithnumval = ["FH","GWP100","PMIN","TMIN","P_REDUCING","PCRIT","GWP20","GAS_CONSTANT","PMAX","RHOCRIT","TCRIT","T_REDUCING","ACENTRIC","GWP500","RHOMOLAR_REDUCING","TMAX","TTRIPLE","PH","M","PTRIPLE","RHOMOLAR_CRITICAL","ODP","HH"];
const fails_any_props_trivals = ["DIPOLE_MOMENT","FRACTION_MAX","FRACTION_MIN","RHOMASS_REDUCING","T_FREEZE"];
const fails_critical_point = ["DiethylEther","R134a","R116","SulfurDioxide","n-Pentane","R11","CycloPropane","MDM","n-Nonane","Oxygen","DimethylCarbonate","R41","R227EA","R245fa","trans-2-Butene","n-Propane","MM","Air","R236FA","Neon","SES36","Fluorine","n-Undecane","Isohexane","MD4M","IsoButane","D5"];
const fails_tcrit_eq_treducing = ["R134a","R116","n-Pentane","R11","n-Nonane","MDM","Oxygen","R41","MM","Neon","Fluorine","n-Undecane","Isohexane","Helium","IsoButane","D5"];
#high
info("********* High Level Api *********");
#get_global_param_string
for param in ["version", "gitrevision", "errstring", "warnstring", "FluidsList", "incompressible_list_pure", "incompressible_list_solution", "mixture_binary_pairs_list", "parameter_list", "predefined_mixtures", "HOME", "cubic_fluids_schema"]
  try
    get_global_param_string(param);
  catch err
    warn("get_global_param_string($param) fails with $err");
  end
end
#PropsSI
@test (PropsSI("T","P",101325.0,"Q",0.0,"Water")-373.1242958476879)<1e-5
#get_fluid_param_string PhaseSI PropsSI
include("testFluids.jl");
#get_parameter_information_string PropsSI
include("testParameters.jl");
#PropsSI
include("testConstants.jl");
#AbstractState_factory get_input_pair_index AbstractState_update get_param_index AbstractState_keyed_output AbstractState_free AbstractState_set_binary_interaction_double AbstractState_set_fractions
include("testLow.jl");
#set_reference_stateS
#F2K K2F
@test K2F(F2K(100)) ≈ 100
#HAPropsSI
dt=1e-3;
@test_approx_eq_eps (HAPropsSI("H", "T", 300+dt, "P", 100000, "Y", 0)-HAPropsSI("H", "T", 300-dt, "P", 100000, "Y", 0))/2/dt HAPropsSI("C", "T", 300, "P", 100000, "Y", 0) 2e-9
#PropsSI
h0 = -15870000.0; # J/kg
s0 = 3887.0; #J/kg
rho0 = 997.1;
T0 = 298.15;
M = PropsSI("molemass", "Water");
#set_reference_stateD("Water", T0, rho0/M, h0*M, s0*M);
#@test PropsSI("H", "T", T0, "P", 101325, "Water") ≈ -1.5870107493843542e7
set_reference_stateS("Water", "DEF");
@test PropsSI("H", "T", T0, "P", 101325, "Water") ≈ 104920.1198093371

if (haskey(ENV, "includelocalwrapper") && ENV["includelocalwrapper"]=="on")
  #saturation_ancillary
  @test saturation_ancillary("R410A","I",1,"T", 300) ≈ 0.004877519938463293
  #cair_sat
  @test cair_sat(400) ≈ 202.53656679999574
end
#config
if (haskey(ENV, "includelocalwrapper") && ENV["includelocalwrapper"]=="on")
  set_config_string("ALTERNATIVE_TABLES_DIRECTORY", "")
end
try
  set_config("MAXIMUM_TABLE_DIRECTORY_SIZE_IN_GB", 1.0)
catch err
  warn("set_config for double fails with: $err")
end
try
  set_config("NORMALIZE_GAS_CONSTANTS", true)
catch err
  warn("set_config for bool fails with: $err")
end


text:
using StatsBase
using Base.Test

fnecdf = ecdf(randn(10000000))
@test_approx_eq_eps fnecdf([-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]) [0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975] 1e-3
@test_approx_eq_eps fnecdf(1.96) 0.975 1e-3
@test_approx_eq fnecdf([-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]) map(fnecdf, [-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96])

fnecdf = ecdf([0.5])
@test fnecdf([zeros(5000); ones(5000)]) == [zeros(5000); ones(5000)]


text:
@test_nowarn GibbsTypePriors.Vnk_NGG(950,50, 0.5, 0.2)
@test_nowarn GibbsTypePriors.Vnk_2PD(7, 6, 0.5, 0.01)

@test_nowarn GibbsTypePriors.Vnk_NGG_rec(100,1, 0.5, 0.2)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
module KittyTerminalImages

using Base64: base64encode
using Rsvg
using Cairo: FORMAT_ARGB32, CairoImageSurface, CairoContext
import Cairo
using Base.Multimedia: xdisplayable
using ImageCore: RGBA, channelview
using CodecZlib: ZlibCompressor, ZlibCompressorStream
using Requires: @require
using Interpolations: interpolate, BSpline, Linear
using PNGFiles

import Base: display

export pushKittyDisplay!, forceKittyDisplay!, set_kitty_config!, get_kitty_config

struct KittyDisplay <: AbstractDisplay end

include("configuration.jl")
include("images.jl")

function __init__()
    # TODO verify that we are actually using kitty
    pushKittyDisplay!()
    @require Compose="a81c6b42-2e10-5240-aca2-a61377ecd94b" include("integration/Compose.jl")
end

function draw_temp_file(img)

    # TODO ensure that there is no racing condition with these tempfiles
    path, io = mktemp()
    PNGFiles.save(io, img)
    close(io)

    payload = transcode(UInt8, base64encode(path))
    write_kitty_image_escape_sequence(stdout, payload, f=100, t='t', X=1, Y=1, a='T')
end

function draw_direct(img)

    # TODO this adds some unnecessary channels for alpha and colors that are not always necessary
    # TODO might be easier to write to a png then we have some compression, then we can also maybe remove zlib
    img_rgba = permutedims(channelview(RGBA.(img)), (1, 3, 2))
    img_encoded = base64encode(ZlibCompressorStream(IOBuffer(vec(reinterpret(UInt8, img_rgba)))))

    (_, width, height) = size(img_rgba)

    buff = IOBuffer()
    partitions = Iterators.partition(transcode(UInt8, img_encoded), 4096)
    for (i, payload) in enumerate(partitions)

        m = (i == length(partitions)) ? 0 : 1 # 0 if this is the last data chunk

        if i == 1
            write_kitty_image_escape_sequence(buff, payload, f=32, s=width, v=height, X=1, Y=1, a='T', o='z', m=m)
        else
            write_kitty_image_escape_sequence(buff, payload, m=m)
        end

    end

    write(stdout, take!(buff))

    return
end

# allows to define custom behaviour for special cases of show within this package
show_custom(io::IO, m::MIME, x) = show(io, m , x)

# values for control data: https://sw.kovidgoyal.net/kitty/graphics-protocol.html#control-data-reference
function write_kitty_image_escape_sequence(io::IO, payload::AbstractVector{UInt8}; controll_data...)

    cmd_prefix = "\033_G"
    first_iteration = true
    for (key, value) in controll_data
        if first_iteration
            first_iteration = false
        else
            cmd_prefix *= ','
        end
        cmd_prefix *= string(key)
        cmd_prefix *= '='
        cmd_prefix*= string(value)
    end
    cmd_prefix *= ';'
    cmd_postfix = "\033\\"
    cmd = [transcode(UInt8, cmd_prefix); payload; transcode(UInt8, cmd_postfix)]
    write(io, cmd)

    return
end

# This is basically just what is defined in Base.Multimedia.display(x)
# but tries to display with KittyDisplay before trying the rest of
# the stack.
function _display(@nospecialize x)
    displays = Base.Multimedia.displays
    for d in reverse(vcat(displays, KittyDisplay()) )
        if xdisplayable(d, x)
            try
                display(d, x)
                return
            catch e
                if !(isa(e, MethodError) && e.f in (display, show))
                    rethrow()
                end
            end
        end
    end
    throw(MethodError(display, (x,)))
end

function forceKittyDisplay!()
    @eval display(@nospecialize x) = _display(x)
    return
end

function pushKittyDisplay!()
    d = Base.Multimedia.displays
    if !isempty(d) && !isa(d[end], KittyDisplay)
        Base.Multimedia.pushdisplay(KittyDisplay())
    end
    return
end

# Supported mime types, they are tried in order that they appear of the list returned
# svg should be preferred png so that we can apply scaling to a vector graphics instead of
# pixels if both formats are supported but because of a bug, (https://github.com/simonschoelly/KittyTerminalImages.jl/issues/4),
# that has not been solved yet, some svg's are not rendered correctly
function kitty_mime_types()
    if get_kitty_config(:prefer_png_to_svg)
        [MIME"image/png"(), MIME"image/svg+xml"()]
    else
        [MIME"image/svg+xml"(), MIME"image/png"()]
    end
end

function display(d::KittyDisplay, x)
    for m in kitty_mime_types()
        if showable(m, x)
            display(d, m, x)
            return
        end
    end
    throw(MethodError(display, (x,)))
end

function display(d::KittyDisplay,
                 m::MIME"image/png", x; scale=get_kitty_config(:scale, 1.0))
    buff = IOBuffer()
    show_custom(buff, m, x)
    seekstart(buff)  # we need to reset the IOBuffer to it's start
    img = PNGFiles.load(buff)
    img = imresize(img; ratio=scale)

    if get_kitty_config(:transfer_mode) == :direct
        draw_direct(img)
    else
        draw_temp_file(img)
    end

    return
end

function display(d::KittyDisplay,
                 m::MIME"image/svg+xml", x; scale=get_kitty_config(:scale, 1.0))
    # Write x to a cairo buffer a and the use the png display method
    buff = IOBuffer()
    show_custom(buff, m, x)
    svg_data = String(take!(buff))
    handle = Rsvg.handle_new_from_data(svg_data)
    dims = Rsvg.handle_get_dimensions(handle)
    width = round(Int, dims.width * scale)
    height = round(Int, dims.height * scale)
    surface = CairoImageSurface(width, height, FORMAT_ARGB32)
    context = CairoContext(surface)
    Cairo.scale(context, scale, scale)
    Rsvg.handle_render_cairo(context, handle)
    # Rsvg.handle_free(handle) # this leads to error messages
    display(d, MIME"image/png"(), surface; scale=1.0) # scaling already happened to svg
    return
end

randkitty() = '\U1f600' + rand([-463, -504, 63, 57, 64, 60, 62, 59, 58, 56, 61, 897, -465])

end # module


text:
using Match

input = joinpath(@__DIR__, "input")
input = readlines(input)

INACTIVE = "."
ACTIVE = "#"

initial = map(line -> split(line, ""), input)
initial = reshape(vcat(initial...), :, length(initial))

function locate(initial, dims)
    (x, y) = size(initial)
    locs = []
    for i in 1:x
        for j in 1:y
            if initial[CartesianIndex(i, j)] == ACTIVE
                array = zeros(Int, dims)
                array[1] = i; array[2] = j
                push!(locs, CartesianIndex(array...))
            end
        end
    end
    return locs
end

function adjacent(dims)
    arr = zeros(Int, dims)
    origin = CartesianIndex(arr...)
    return filter(coord -> coord != origin, CartesianIndices((repeat([-1:1], dims)...,)))
end

function iter(sources, dims, cycles)
    next = Set(copy(sources))
    for _ in 1:cycles
        field = next
        neighbors = Dict{CartesianIndex{dims}, Int}()
        for source in field
            for adj in map(coord -> coord + source, adjacent(dims))
                haskey(neighbors, adj) ? neighbors[adj] = neighbors[adj] + 1 : neighbors[adj] = 1
            end
        end
        next = Set{CartesianIndex{dims}}()
        for (coord, count) in neighbors
            count == 3 && push!(next, coord)
            coord in field && count in 2:3 && push!(next, coord)
        end
    end
    return collect(next)
end

dims = 3
locs = locate(initial, dims)
locs = iter(locs, dims, 6)
p1 = length(locs)

dims = 4
locs = locate(initial, dims)
locs = iter(locs, dims, 6)
p2 = length(locs)

@assert(p1 == 257)
@assert(p2 == 2532)

println("-----------------------------------------------------------------------")
println("conway cubes -- part one :: $p1")
println("conway cubes -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
# I had some problems to generate plot during CI: here is a script to
# generate them before committing...  This is ugly and not scalable,
# but for the moment this is the only solution I found.
#

using GnuplotScripting
figdir = joinpath(@__DIR__,"figures")
figfile(f::AbstractString)=joinpath(figdir,f)

# Script 1
#
gp = GnuplotScript()

X=[-pi:0.1:pi;];
Ys =sin.(X);
Yc =cos.(X);

id=register_data(gp,hcat(X,Ys,Yc))
free_form(gp,"replot '$id' u 1:3 w l t 'cos'")
free_form(gp,"replot '$id' u 1:2 w l t 'sin'")

export_png(gp, figfile("trig.png"))

# Script 2
#
X = -2:0.1:2
Y = -2:0.1:2
M = [exp(-x_i*x_i-y_j*y_j) for x_i=X, y_j=Y]

gp = GnuplotScript()

id = register_data(gp, M)

free_form(gp,"set autoscale fix")
free_form(gp,"plot '$id' matrix using 1:2:3 with image")

export_png(gp, figfile("2D.png"))

# vlines
#
using GnuplotScripting

gp = GnuplotScript()

add_vertical_line(gp,-5.0,name="left")
add_vertical_line(gp,+5.0,name="right")

free_form(gp,"plot exp(-x*x/25) with line t 'Gaussian'")

export_png(gp, figfile("vlines.png"))


text:
using YaoBase, TupleTools

export focus!, relax!, partial_tr, exchange_sysenv

"""
    contiguous_shape_orders(shape, orders)

Merge the shape and orders if the orders are contiguous. Returns the
new merged shape and order.

# Example

```jldoctest; setup=:(using YaoArrayRegister)
julia> YaoArrayRegister.contiguous_shape_orders((2, 3, 4), (1, 2, 3))
([24], [1])
```
"""
function contiguous_shape_orders(shape, orders)
    new_shape, new_orders = Int[], Int[]
    prv = -1
    for cur in orders
        if cur == prv + 1
            new_shape[end] *= shape[cur]
        else
            push!(new_orders, cur)
            push!(new_shape, shape[cur])
        end
        prv = cur
    end
    # NOTE: some of the orders are merged above
    #       we use sortperm to retrieve correct
    #       orders
    inv_orders = sortperm(new_orders)
    return new_shape[inv_orders], invperm(inv_orders)
end

# NOTE: don't use Vector for move_ahead, it's way slower!
# Before:
# julia> @benchmark move_ahead($c, $o)
# BenchmarkTools.Trial:
#   memory estimate:  14.87 KiB
#   allocs estimate:  392
#   --------------
#   minimum time:     29.062 μs (0.00% GC)
#   median time:      31.316 μs (0.00% GC)
#   mean time:        39.724 μs (14.37% GC)
#   maximum time:     39.013 ms (99.86% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1
# After:
# julia> @benchmark move_ahead($c, $o)
# BenchmarkTools.Trial:
#   memory estimate:  4.04 KiB
#   allocs estimate:  24
#   --------------
#   minimum time:     2.848 μs (0.00% GC)
#   median time:      3.045 μs (0.00% GC)
#   mean time:        4.013 μs (18.20% GC)
#   maximum time:     4.494 ms (99.86% GC)
#   --------------
#   samples:          10000
#   evals/sample:     9

"""
    move_ahead(collection, orders)

Move `orders` to the beginning of `collection`.
"""
move_ahead(collection, orders) = (orders..., setdiff(collection, orders)...)
move_ahead(ndim::Int, orders) = (orders..., setdiff(1:ndim, orders)...)

function group_permutedims(A::AbstractArray, orders)
    @assert length(orders) == ndims(A) "number of orders does not match number of dimensions"
    return unsafe_group_permutedims(A, orders)
end

# forward directly if the length is the same with ndims
function group_permutedims(A::AbstractArray{T,N}, orders::NTuple{N,Int}) where {T,N}
    return unsafe_group_permutedims(A, orders)
end

function unsafe_group_permutedims(A::AbstractArray, orders)
    s, o = contiguous_shape_orders(size(A), orders)
    return permutedims(reshape(A, s...), o)
end

"""
    is_order_same(locs) -> Bool

Check if the order specified by `locs` is the same as current order.
"""
is_order_same(locs) = all(a == b for (a, b) in zip(locs, 1:length(locs)))

# NOTE: locations is not the same with orders
# locations: some location of the wire
# orders: includes all the location of the wire in some order
function YaoBase.focus!(r::AbstractArrayReg{D}, locs) where {D}
    if is_order_same(locs)
        arr = r.state
    else
        new_orders = move_ahead(nactive(r) + 1, locs)
        arr = group_permutedims(hypercubic(r), new_orders)
    end
    r.state = reshape(arr, D^length(locs), :)
    return r
end

function YaoBase.relax!(r::AbstractArrayReg{D}, locs; to_nactive::Int = nqudits(r)) where {D}
    r.state = reshape(state(r), D^to_nactive, :)
    if !is_order_same(locs)
        new_orders = TupleTools.invperm(move_ahead(to_nactive + 1, locs))
        r.state = reshape(group_permutedims(hypercubic(r), new_orders), D^to_nactive, :)
    end
    return r
end

function YaoBase.partial_tr(r::AbstractArrayReg{D}, locs) where D
    orders = setdiff(1:nqudits(r), locs)
    r2 = focus!(copy(r), orders)
    state = sum(rank3(r2); dims = 2)
    return normalize!(arrayreg(reshape(state, :, size(state, 3)); nbatch=nbatch(r), nlevel=D))
end

"""
    exchange_sysenv(reg::AbstractArrayReg) -> AbstractRegister

Exchange system (focused qubits) and environment (remaining qubits).
"""
function exchange_sysenv(reg::AbstractArrayReg)
    r3 = rank3(reg)
    M, N, B = size(r3)
    arrayreg(reshape(permutedims(r3, (2, 1, 3)), :, M*B); nbatch=nbatch(reg), nlevel=nlevel(reg))
end


text:
module TestShow

using Test
using Try
using TryExperimental: @tryable, NotImplementedError

@tryable dummy

function test_tryable()
    @test string(dummy) == "dummy"
    @test sprint(show, dummy) == "dummy"
    @test sprint(show, "text/plain", dummy) == "dummy (tryable function with 1 method)"
end

function test_notimplementederror()
    ex = NotImplementedError(identity, (1, 2, 3), (a = 4, b = 5))
    msg = sprint(showerror, ex)
    @test occursin("identity(1, 2, 3; a = 4, b = 5)", msg)
end

end  # module


text:
using Documenter, Interact, Literate
src = joinpath(@__DIR__, "src")
Literate.markdown(joinpath(src, "tutorial.jl"), src, codefence = "```julia" => "```")

makedocs(
    sitename = "Interact",
    authors = "JuliaGizmos",
    pages = [
        "Introduction" => "index.md",
        "Observables" => "observables.md",
        "Widgets" => "widgets.md",
        "Custom widgets" => "custom_widgets.md",
        "Modifiers" => "modifiers.md",
        "Layout" => "layout.md",
        "Deploying the web app" => "deploying.md",
        "Tutorial" => "tutorial.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaGizmos/Interact.jl.git",
)


text:
## -------------  Model

# Make random model parameters for testing and intialization
# All in [1, 2]
# Changes GLOBAL_RNG
function make_random_parameters(mdl :: Model, randSeed :: Integer)
    Random.seed!(randSeed);
    return ModelParams(rand(n_params(mdl)) .+ 1.0)
end

# Objective function to be minimized
function objective(mdl :: Model,  tgStats :: ModelStats,  mp :: ModelParams)
    mStats = solve(mdl, mp);

    # Deviation
    dev = deviation(mdl, tgStats, mStats);
    return dev
end

# Guess to model params and reverse
function guess_to_mp(mdl :: Model, guessV :: Vector{Float64})
    @assert length(guessV) == n_params(mdl)
    return ModelParams(guessV)
end

function mp_to_guess(mdl :: Model, mp :: ModelParams)
    return mp.alphaV
end

## -----------  Solve

function solve(mdl :: Model, mp :: ModelParams)
    # endow_jV = make_endowments(mdl,  mp);
    # optValue_mV = make_option_values(mdl, mp);

    prob_mjM = decision_probabilities(mdl, mp);
    choice_sjM = simulate_choices(mdl, prob_mjM, 320);
    mStats = compute_stats(mdl, choice_sjM);

    return mStats
end

## Make endowments of households
function make_endowments(mdl :: Model)

    # Random.seed!(583);
    # outV = rand(mdl.N, mdl.K + 1) * vcat(1.0, mp.alphaV);
    outV = collect(1 : mdl.N) ./ mdl.N;
    return outV :: Vector{Float64}
    # endow_kjM = zeros(mdl.K, nInd);
    # endowV = collect(1 : nInd) ./ nInd .* 2;
    # for k = 1 : mdl.K
    #     endow_kjM[k,:] = (endowV .^ 0.6) .* (alphaV[k] .+ betaV[k] .* endowV);
    # end
    # return endow_kjM
end

## Option value = fixed matrix * alphaV
function make_option_values(mdl :: Model, mp :: ModelParams)
    Random.seed!(439);
    outV = rand(mdl.M, mdl.K + 1) * vcat(1.0, mp.alphaV);
    return outV :: Vector{Float64}
end

## Decision probabilities
function decision_probabilities(mdl :: Model, mp :: ModelParams)
    v_mV = make_option_values(mdl, mp);
    value_jV = make_endowments(mdl);
    # The trick is to find a bounded function with many parameters
    util_mjM = exp.(sin.(5.0 .* v_mV * transpose(value_jV)));
    prob_mjM = util_mjM ./ sum(util_mjM, dims = 1);
    # prob_mjM = zeros(mdl.M, mdl.N);
    # for j = 1 : mdl.N
    #     prob_mjM[:, j] = decision_prob(mdl, endow_jV[j],  optValue_mV);
    # end

    @assert all(prob_mjM .> 0.0001)
    return prob_mjM
end

## Decision probabilties for one household
# Inputs are endowments and option characteristics
# Must be nonlinear. Otherwise some parameters do not affect decision probs (intercepts)
# function decision_prob(mdl :: Model, endow :: Float64, optValue_mV :: Vector{Float64})
#     # Utility of alternatives
#     utilV = zeros(mdl.M);
#     for m = 1 : mdl.M
#         # Ensures positive utility
#         utilV[m] = sum((endow_kV .* optValue_kmM[:, m]) .^ 0.5);
#     end
#     @assert all(utilV .> 0.0)
#     @assert all(utilV .< 1e5)
#     probV = utilV ./ sum(utilV);
#     @assert all(probV .> 0.0)
#     @assert all(probV .< 1.0)
#     return probV
# end

# This changes GLOBAL_RNG
function simulate_choices(mdl :: Model, prob_mjM, randSeed :: Integer)
    Random.seed!(randSeed);
    choice_sjM = zeros(Int64, mdl.nSim, mdl.N);
    for j = 1 : mdl.N
        d = Categorical(prob_mjM[:, j]);
        choice_sjM[:, j] = rand(d, mdl.nSim);
    end
    return choice_sjM
end

"""
	compute_stats

Compute model statistics
"""
function compute_stats(mdl :: Model,  choice_sjM :: Matrix{T1})  where T1 <: Integer

    choiceCl_mV, Nc = choice_classes(mdl);
    endowCl_jV, Nendow = endowment_classes(mdl);

    # Count by [choice class, endowment class]
    cnt_ceM = zeros(Int64,  Nc, Nendow);
    for j = 1 : mdl.N
        # Endowment class for this agent
        iEndow = endowCl_jV[j];
        # Choice class by simulated instance of hh j
        choiceClV = choiceCl_mV[choice_sjM[:, j]];
        for ic = 1 : Nc
            # Count how many simulated agents choose this choice class
            cnt_ceM[ic, iEndow] += round(Int64, sum(choiceClV .== ic));
        end
    end
    frac_ceM = cnt_ceM ./ sum(cnt_ceM);
    mStats = ModelStats(frac_ceM);
    return mStats
end

function choice_classes(mdl :: Model)
    Nendow = min(mdl.nChoiceCl, mdl.M - 2);
    pctUbV = collect(range(1.0 ./ Nendow, 1.0, length = Nendow));
    choiceCl_kV = discretize_given_percentiles(collect(1.0 : mdl.M), pctUbV, false);
    return choiceCl_kV, Nendow
end

function endowment_classes(mdl :: Model);
    Nendow = min(mdl.nEndowCl, mdl.N - 2);
    endow_jV = make_endowments(mdl);
    pctUbV = collect(range(1.0 ./ Nendow, 1.0, length = Nendow));
    endowCl_jV = discretize_given_percentiles(endow_jV, pctUbV, false);

    return endowCl_jV, Nendow
end

## ------------  Deviation

function deviation(mdl :: Model,  tgStats :: ModelStats,  mStats :: ModelStats)
    dev_ceM = mStats.frac_ceM .- tgStats.frac_ceM;
    dev = sum(abs.(dev_ceM));
    mdl.currentIter += 1;

    if mdl.currentIter < 100
        showFreq = 10;
    else
        showFreq = mdl.showFreq;
    end
    if rem(mdl.currentIter, showFreq) == 0
        println("Iter $(mdl.currentIter).   Dev: $dev");
    end
    return dev :: Float64
end

# Test the objective function with random guesses
# Because NLopt returns FORCED_STOP if an error occurs
function calc_many_deviations()
    mdl = Model(showFreq = 20)
    nRuns = 400;
    # showFreq = 20;

    mp = make_random_parameters(mdl, 343);
    tgStats = solve(mdl, mp);
    dev0 = objective(mdl, tgStats, mp);

    Random.seed!(449);
    seedV = 123 .+ round.(Int,  1000 .* rand(nRuns));

    for i1 = 1 : nRuns
        mpGuess = make_random_parameters(mdl, seedV[i1]);
        devGuess = objective(mdl, tgStats, mpGuess);
        @assert abs(devGuess - dev0) > 0.001
        # guessV = mp_to_guess(mdl, mpGuess);
        # dev0 = dev_fct(guessV, []);
        # @assert abs(dev0 - dev00) < 0.001
        # if rem(i1, showFreq) == 0
        #     println("Iteration $i1.  Deviation $dev0")
        #     # println(guessV)
        # end
    end

    return true
end

## Check that changing guesses changes solution
function check_changing_guesses(mdl)
    mp = make_random_parameters(mdl, 932);
    tgStats = solve(mdl, mp);
    guessV = mp_to_guess(mdl, mp);
    dev0 = objective(mdl, tgStats, mp);

    equalV = fill(false, n_params(mdl));
    for i1 = 1 : n_params(mdl)
        guess2V = deepcopy(guessV);
        guess2V[i1] += 0.1;
        mp2 = OptimizationLH.guess_to_mp(mdl, guess2V);
        dev2 = OptimizationLH.objective(mdl, tgStats, mp2);
        equalV[i1] = abs(dev2 - dev0) < 0.001;
    end

    success = all(equalV .== false);
    if success
        println("All guesses affect objective");
    else
        @warn "Not all guesses affect objective"
        println("Guesses that do not:");
        println(findall(equalV));
    end

    return success
end

# --------------


text:
"""
    sequentialise(data)

Convert an array of arrays to a concatenated array throught the first dimension.

It would be useful when plotting figures using data obtained from Transducers.

# Examples
```jldoctest
julia> using LazyFym

julia> sequentialise([[1, 2], [3, 4], [5, 6]])
3×2 Array{Int64,2}:
 1  2
 3  4
 5  6
```
"""
function sequentialise(data)
    return vcat([reshape(data[i], (1, size(data[i])...)) for i in 1:length(data)]...)
end


text:
using PyPlot, AcousticWave
# =============generate synthetic data set================================
nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
v  = 3000. * ones(nz, nx);
# initialize source term
pos_src = [101 75; 101 125]; f0 = 30.0; ot = [0.0, 0.0];
src = InitSources(pos_src, f0, ot, dt);
vmax = maximum(v); vmin = minimum(v);
fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
pos_rec = hcat(iz, ix);
shot = MultiStepForward(pos_rec, src, fidMtx, tmax=0.5);
wt = Ricker(f0, dt); ntw = length(wt);
mu = 0.1; lambda = 120.0; tmax = 0.5; maxit = 10;

(dis, obj) =CG_dis(shot, fidMtx, pos_rec, wt, maxit=10, mu=0.001, tmax=0.5)
dis = reshape(dis, nz, nx)
SeisPlot(dis, name="/Users/wenlei/Desktop/CG.pdf", dx=10, dy=10, xlabel="Length", xunits="(m)", ylabel="Depth", yunits="(m)")
dis = zeros(Float64,nz,nx)
dis[pos_src[1,1], pos_src[1,2]] = 1.; dis[pos_src[2,1], pos_src[2,2]] = 1.;
dis = vec(dis)
ntw = 5*ntw;
(w, obj) = CG_w(ntw, shot, fidMtx, pos_rec, dis, maxit=10, mu=0.001, tmax=0.5)
plot(w);
plot(wt);
xlable = "Time (s)"
ylable = "Amplitude"
savefig("/Users/wenlei/Desktop/wavelet.pdf")

dis = fista(shot, wt, mu, lambda, fidMtx, pos_rec, tmax, maxit)

# (dis, w) = alterMini(w, shot, fidMtx, pos_rec, mu, lambda, outer_it=6, fistait=10, CGit=5, tmax=0.5)
# (dis, w) = alternateCG(w, shot, fidMtx, pos_rec, outer_it=6)

# ============================================================================================================
# using PyPlot, AcousticWave
# nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3; f0 = 30.0
# v  = 2500. * ones(nz, nx);
# v[101:end, 1:end] = 3000.;
# vmax = maximum(v); vmin = minimum(v);
#
# Nz = nz + 2*ext; Nx = nx + 2*ext
# w = Ricker(f0, dt); dis = randn(nz*nx); it = floor(Int64, length(w)/2);
# spt = dis2spt(w, dis, nz, nx, ext, iflag, dt, it);
# # SeisPlot(reshape(spt.pz, Nz, Nx));
# # SeisPlot(reshape(spt.px, Nz, Nx));
#
# tmp = reshape(spt.pz, Nz, Nx) + reshape(spt.px, Nz, Nx);
# tmp = tmp[ext+1:ext+nz, ext+1:ext+nx];
# norm(w[it]*dis-vec(tmp))
#
# spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, it)
# spt1.vz = randn(length(spt1.vz));
# spt1.vx = randn(length(spt1.vx));
# spt1.pz = randn(length(spt1.pz));
# spt1.px = randn(length(spt1.px));
# dis1 = spt2dis(w, spt1);
# tmp_spt = dot(spt.vz, spt1.vz) + dot(spt.vx, spt1.vx) + dot(spt.pz, spt1.pz) + dot(spt.px, spt1.px)
# tmp_dis = dot(dis, dis1)

# ===========================================================================================================
# using PyPlot, AcousticWave
# nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3; f0 = 30.0
# v  = 2500. * ones(nz, nx);
# v[101:end, 1:end] = 3000.;
# vmax = maximum(v); vmin = minimum(v);
# w = Ricker(f0, dt); it = 33;
# spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, it)
# dis  = zeros(nz, nx); dis[101,75] = 1.0; dis[101,125] = 1.0; dis = vec(dis);
# AddSources!(spt2, w, dis)
# spt3 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, it)
# f0 = 30.0; ot = [0.0,0.0]; pos_src = [101 75;101 125]
# src = InitSources(pos_src, f0, ot, dt);
# AddSources!(spt3, src)
# norm(spt2.vx-spt3.vx)
# norm(spt2.vz-spt3.vz)
# norm(spt2.pz-spt3.pz)
# norm(spt2.px-spt3.px)

# =============================================================================================================
using PyPlot, AcousticWave
nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
v  = 3000. * ones(nz, nx);
# initialize source term
pos_src = [101 75; 101 125]; f0 = 30.0; ot = [0.0, 0.0];
src = InitSources(pos_src, f0, ot, dt);
vmax = maximum(v); vmin = minimum(v);
fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
pos_rec = hcat(iz, ix);
w = Ricker(f0, dt);
dis = zeros(nz, nx); dis[pos_src[1,1], pos_src[1,2]] = 1.0;
dis[pos_src[2,1], pos_src[2,2]] = 1.0; dis = vec(dis);
shot = MultiStepForward(pos_rec, w, dis, fidMtx, tmax=0.5);
dis1 = MultiStepAdjoint(w, shot, fidMtx)

# ============================================================================================================
# dot product test
# using PyPlot, AcousticWave
# nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
# v  = 3000. * ones(nz, nx);
# # initialize source term
# pos_src = [101 75; 101 125]; f0 = 30.0; ot = [0.0, 0.0];
# src = InitSources(pos_src, f0, ot, dt);
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
# iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
# pos_rec = hcat(iz, ix);
# w = Ricker(f0, dt);
# dis = zeros(nz, nz);
# for ix = 1: nx
#     for iz = 1 : nz
#         if rand() <= 0.1
#            dis[iz, ix] = randn()
#         end
#     end
# end
# dis[pos_src[2,1], pos_src[2,2]] = 1.0; dis = vec(dis);
# shot = MultiStepForward(pos_rec, w, dis, fidMtx, tmax=0.5);
#
# nt = shot.nt
# shot1 = InitShot(0, 0, pos_rec, 0.0, nt, dt)
# nrec  = size(pos_rec, 1)
# ls    = floor(Int64, length(w)/2)
# for irec = 1 : nrec
#     tmp  = conv(rand(nt), w)
#     shot1.d[:, irec] = tmp[ls+1:end-ls]
# end
# dis1 = MultiStepAdjoint(w, shot1, fidMtx)
#
# tmp  = dot(dis, dis1)
# tmp1 = dot(vec(shot.d), vec(shot1.d))

# ============================================================================================================
# using PyPlot, AcousticWave
# nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3; f0 = 30.0
# v  = 2500. * ones(nz, nx);
# v[101:end, 1:end] = 3000.;
# vmax = maximum(v); vmin = minimum(v);
# Nz = nz + 2*ext; Nx = nx + 2*ext
# w = Ricker(f0, dt); dis = randn(nz*nx);
# nt = length(w)
# spts = Array(SnapShot, nt)
# for it = 1 : nt
#     spts[it] = dis2spt(w, dis, nz, nx, ext, iflag, dt, it)
# end
# it = 36
# a = reshape(spts[it].pz, 261, 261);
# b = a[ext+1:ext+nz, ext+1:ext+nx];
# tmp = b ./ reshape(dis, nz, nx)
#
#
# spts1 = Array(SnapShot, nt)
# for it = 1 : nt
#     spts1[it] = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, it)
#     spts1[it].vz = randn(length(spts1[it].vz))
#     spts1[it].vx = randn(length(spts1[it].vx))
#     spts1[it].pz = randn(length(spts1[it].pz))
#     spts1[it].px = randn(length(spts1[it].px))
# end
# w1 = zeros(nt)
# for it = 1 : nt
#     spt2wlet(w1, dis, spts1[it])
# end
#
# tmp = 0.0
# for it = 1 : nt
#     tmp = tmp + dot(spts[it].vz, spts1[it].vz)
#     tmp = tmp + dot(spts[it].vx, spts1[it].vx)
#     tmp = tmp + dot(spts[it].pz, spts1[it].pz)
#     tmp = tmp + dot(spts[it].px, spts1[it].px)
# end
# tmp1 = dot(w, w1)

# ============================================================================================================
# dot product test
# using PyPlot, AcousticWave
# nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
# v  = 3000. * ones(nz, nx);
# # initialize source term
# pos_src = [101 75; 101 125]; f0 = 30.0; ot = [0.0, 0.0];
# src = InitSources(pos_src, f0, ot, dt);
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
# iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
# pos_rec = hcat(iz, ix);
# w = Ricker(f0, dt);
# dis = zeros(nz, nz);
# for ix = 1: nx
#     for iz = 1 : nz
#         if rand() <= 0.1
#            dis[iz, ix] = randn()
#         end
#     end
# end
# dis = vec(dis)
# shot = MultiStepForward(pos_rec, w, dis, fidMtx, tmax=0.5);
#
#
# nt = shot.nt
# shot1 = InitShot(0, 0, pos_rec, 0.0, nt, dt)
# nrec  = size(pos_rec, 1)
# ls    = floor(Int64, length(w)/2)
# for irec = 1 : nrec
#     tmp  = conv(rand(nt), w)
#     shot1.d[:, irec] = tmp[ls+1:end-ls]
# end
# ntw = length(w)
# w1  = MultiStepAdjoint(ntw, dis, shot1, fidMtx)
#
# tmp  = dot(w, w1)
# tmp1 = dot(vec(shot.d), vec(shot1.d))

# ===========================================================================================================
using PyPlot, AcousticWave
# =============generate synthetic data set================================
nz = 201; nx = 201; ext = 30; iflag = 1; dz = 5.; dx = 5.; dt = 1e-3;
v  = 3000. * ones(nz, nx);
# initialize source term
pos_src = [101 75; 101 125]; f0 = 30.0; ot = [0.0, 0.0];
src = InitSources(pos_src, f0, ot, dt);
vmax = maximum(v); vmin = minimum(v);
fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
iz = collect(5:5:201); ix = 20 * ones(Int64, length(iz));
pos_rec = hcat(iz, ix); tmax = 0.5; maxit=20;
w = Ricker(f0, dt)
lambda = power(maxit, fidMtx, w, pos_rec, tmax)

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
"""
Representation of a state for the SingleOIcludedIntersection POMDP,
depends on ADM VehicleState type
"""
mutable struct SingleOIState
    crash ::Bool
    ego::VehicleState
    car::VehicleState
end

# copy b to a
function Base.copyto!(a::SingleOIState, b::SingleOIState)
    a.crash = b.crash
    a.ego = b.ego
    a.car = b.car
end

function Base.hash(s::SingleOIState, h::UInt64 = zero(UInt64))
    return hash(s.crash, hash(s.ego, hash(s.car, h)))
end

function Base.:(==)(a::SingleOIState, b::SingleOIState)
    return a.crash == b.crash && a.ego == b.ego && a.car == b.car
end

#### Observaton type

const SingleOIObs = SingleOIState

#### Action type

mutable struct SingleOIAction
    acc::Float64
end

function Base.copyto!(a::SingleOIAction, b::SingleOIAction)
    a.acc = b.acc
end

function Base.hash(a::SingleOIAction, h::UInt64 = zero(UInt64))
    return hash(a.acc, h)
end

function Base.:(==)(a::SingleOIAction, b::SingleOIAction)
    return a.acc == b.acc
end

#### POMDP type

mutable struct SingleOIPOMDP <: POMDP{SingleOIState, SingleOIAction, SingleOIObs}
    env::SimpleInterEnv
    ego_type::VehicleDef
    car_type::VehicleDef
    max_acc::Float64
    pos_res::Float64
    vel_res::Float64
    s_start::Float64
    s_goal::Float64
    off_grid::VecSE2
    ΔT::Float64 # decision frequency
    p_birth::Float64
    a_noise::Float64
    pos_obs_noise::Float64
    vel_obs_noise::Float64
    collision_cost::Float64
    action_cost::Float64
    goal_reward::Float64
    γ::Float64 # discount factor
end

function SingleOIPOMDP(; env::SimpleInterEnv = SimpleInterEnv(),
                   ego_type::VehicleDef = VehicleDef(),
                   car_type::VehicleDef = VehicleDef(),
                   max_acc::Float64 = 2.0,
                   pos_res::Float64 = 2.,
                   vel_res::Float64 = 2.,
                   s_start::Float64 = env.params.stop_line - ego_type.length/2,
                   s_goal::Float64 = 30.,
                   off_grid::VecSE2 = VecSE2(14., -5, pi/2),
                   ΔT::Float64  = 0.5,
                   p_birth::Float64 = 0.1,
                   a_noise::Float64 = 1.0,
                   pos_obs_noise::Float64 = 1.0,
                   vel_obs_noise::Float64 = 1.0,
                   collision_cost::Float64 = -1.,
                   action_cost::Float64 = 0.0,
                   goal_reward::Float64 = 1.,
                   γ::Float64  = 0.95)
    return SingleOIPOMDP(env,
                   ego_type,
                   car_type,
                   max_acc,
                   pos_res,
                   vel_res,
                   s_start,
                   s_goal,
                   off_grid,
                   ΔT,
                   p_birth,
                   a_noise,
                   pos_obs_noise,
                   vel_obs_noise,
                   collision_cost,
                   action_cost,
                   goal_reward,
                   γ)
end

### REWARD MODEL ##################################################################################

function POMDPs.reward(pomdp::SingleOIPOMDP, s::SingleOIState, a::SingleOIAction, sp::SingleOIState)
    r = 0.
    if sp.crash
        r += pomdp.collision_cost
    end
    if sp.ego.posF.s >= pomdp.s_goal
        r += pomdp.goal_reward
    elseif a.acc > 0.
        r += pomdp.action_cost
    else
        r += pomdp.action_cost
    end
    return r
end

# other method for SARSOP
function POMDPs.reward(pomdp::SingleOIPOMDP, s::SingleOIState, a::SingleOIAction)
    return reward(pomdp, s, a, s)
end

function POMDPs.isterminal(pomdp::SingleOIPOMDP, s::SingleOIState)
    return s.crash || s.ego.posF.s >= pomdp.s_goal
end

######## DISTRIBUTION ############################################################################

"""
Concrete type to represent a distribution over state for the SingleOIcludedCrosswalk POMDP Problem
"""
mutable struct SingleOIDistribution
    p::Vector{Float64}
    it::Vector{SingleOIState}
end

SingleOIDistribution() = SingleOIDistribution(Float64[], SingleOIState[])

# transition and observation pdf
function POMDPs.pdf(d::SingleOIDistribution, s::SingleOIState)
    for (i, sp) in enumerate(d.it)
        if sp==s
            return d.p[i]
        end
    end
    return 0.
end

function POMDPs.rand(rng::AbstractRNG, d::SingleOIDistribution)
    ns = sample(d.it, Weights(d.p)) # sample a neighbor state according to the distribution c
    return ns
end

"""
    most_likely_state(d::SingleOIDistribution)
returns the most likely state given distribution d
"""
function most_likely_state(d::SingleOIDistribution)
    val, ind = fargmax(d.p)
    return d.it[ind]
end

### HELPERS

function POMDPs.discount(pomdp::SingleOIPOMDP)
    return pomdp.γ
end

function state_to_scene(pomdp::SingleOIPOMDP, s::SingleOIState)
    scene = Scene()
    ego = Vehicle(s.ego, pomdp.ego_type, 1)
    car = Vehicle(s.car, pomdp.car_type, 2)
    push!(scene, ego)
    push!(scene, car)
    return scene
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# Calculate the thermal expansivity αᵥ of an EOS

using Calculus: derivative, gradient  # to evaluate derivatives

"Directional derivative d(ln ρ)/dT"
function partiallnT(eos::EOS, P, T)
    derivative(δT -> log(eos(P, T + δT)), 0)::Float64
end
function partiallnT(s::StitchedEOS, P, T)
    i = findfirst(eos -> (P, T) in eos, s.eoses)
    i == 0 ? NaN : partiallnT(s.eoses[i], P, T)
end

"Gradient ∇(ln ρ)"
function gradientlog(eos::EOS, P, T)
    gradient(PT -> log(eos(PT[1], PT[2])), [P, T])::Vector{Float64}
end
function gradientlog(eos::StitchedEOS, P, T)
    i = findfirst(eos -> (P, T) in eos, s.eoses)
    i == 0 ? [NaN, NaN] : gradient(s.eoses[i], P, T)
end

"Thermal expansivity α of an EOS at pressure `P` and temperature `T`"
function thermalexpansivity(eos::EOS, P, T)
    -1 * partiallnT(eos, P, T)
end


text:
module FSM

#using DataAssim

const fsm = joinpath(dirname(@__FILE__), "..", "deps", "FSM")

export FsmType, FsmInput

export run_fsm, run_fsm!

#export q_noise, perturb_input

export precip_liquid, precip_solid

include("run_fsm.jl")
#include("utils_pfilter.jl")
include("utils_meteo.jl")

end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# The Bellman Ford algorithm for single-source shortest path

###################################################################
#
#   The type that capsulates the state of Bellman Ford algorithm
#
###################################################################

type NegativeCycleError <: Exception end

# AbstractPathState is defined in core
type BellmanFordState{T<:Number}<:AbstractPathState
    parents::Vector{Int}
    dists::Vector{T}
end

function bellman_ford_shortest_paths!{R<:Real}(
    graph::SimpleGraph,
    sources::AbstractVector{Int},
    distmx::AbstractArray{R, 2},
    state::BellmanFordState
)

    active = Set{Int}()
    for v in sources
        state.dists[v] = 0
        state.parents[v] = 0
        push!(active, v)
    end
    no_changes = false
    for i in 1:nv(graph)
        no_changes = true
        new_active = Set{Int}()
        for u in active
            for v in fadj(graph, u)
                edist = distmx[u, v]
                if state.dists[v] > state.dists[u] + edist
                    state.dists[v] = state.dists[u] + edist
                    state.parents[v] = u
                    no_changes = false
                    push!(new_active, v)
                end
            end
        end
        if no_changes
            break
        end
        active = new_active
    end
    no_changes || throw(NegativeCycleError())
    return state
end

"""Uses the [Bellman-Ford algorithm](http://en.wikipedia.org/wiki/Bellman–Ford_algorithm)
to compute shortest paths between a source vertex `s` or a set of source
vertices `ss`. Returns a `BellmanFordState` with relevant traversal information
(see below).
"""
function bellman_ford_shortest_paths{T}(
    graph::SimpleGraph,

    sources::AbstractVector{Int},
    distmx::AbstractArray{T, 2} = DefaultDistance()
    )
    nvg = nv(graph)
    state = BellmanFordState(zeros(Int,nvg), fill(typemax(T), nvg))
    bellman_ford_shortest_paths!(graph, sources, distmx, state)
end

bellman_ford_shortest_paths{T}(
    graph::SimpleGraph,
    v::Int,
    distmx::AbstractArray{T, 2} = DefaultDistance()
) = bellman_ford_shortest_paths(graph, [v], distmx)

function has_negative_edge_cycle(graph::SimpleGraph)
    try
        bellman_ford_shortest_paths(graph, vertices(graph))
    catch e
        isa(e, NegativeCycleError) && return true
    end
    return false
end

function enumerate_paths(state::AbstractPathState, dest::Vector{Int})
    parents = state.parents

    num_dest = length(dest)
    all_paths = Array(Vector{Int},num_dest)
    for i=1:num_dest
        all_paths[i] = Vector{Int}()
        index = dest[i]
        if parents[index] != 0 || parents[index] == index
            while parents[index] != 0
                push!(all_paths[i], index)
                index = parents[index]
            end
            push!(all_paths[i], index)
            reverse!(all_paths[i])
        end
    end
    all_paths
end

enumerate_paths(state::AbstractPathState, dest) = enumerate_paths(state, [dest])[1]
enumerate_paths(state::AbstractPathState) = enumerate_paths(state, [1:length(state.parents);])

"""Given a path state `state` of type `AbstractPathState` (see below), returns a
vector (indexed by vertex) of the paths between the source vertex used to
compute the path state and a destination vertex `v`, a set of destination
vertices `vs`, or the entire graph. For multiple destination vertices, each
path is represented by a vector of vertices on the path between the source and
the destination. Nonexistent paths will be indicated by an empty vector. For
single destinations, the path is represented by a single vector of vertices,
and will be length 0 if the path does not exist.

For Floyd-Warshall path states, please note that the output is a bit different,
since this algorithm calculates all shortest paths for all pairs of vertices:
`enumerate_paths(state)` will return a vector (indexed by source vertex) of
vectors (indexed by destination vertex) of paths. `enumerate_paths(state, v)`
will return a vector (indexed by destination vertex) of paths from source `v`
to all other vertices. In addition, `enumerate_paths(state, v, d)` will return
a vector representing the path from vertex `v` to vertex `d`.
"""
enumerate_paths


text:
const _GW_TASK_DAT_FILE_NAME = "dat.jls"
_taskdat_file(taskdir::String) = joinpath(taskdir, _GW_TASK_DAT_FILE_NAME)
_taskdat_file(gwt::GWTask) = _taskdat_file(task_dir(gwt))

const _GW_TASK_EXPR_KEY = "expr"
const _GW_TASK_SRC_KEY = "src"
const _GW_TASK_DESC_KEY = "desc"
const _GW_TASK_LANG_KEY = "lang"

_task_dat(gwt::GWTask) = get!(gwt, _GW_TASK_DAT_FILE_NAME) do
    Dict{String, Any}()
end

function _write_task_dat!(gwt::GWTask; ex, src, desc, lang)

    dfile = _taskdat_file(gwt)
    dat = _task_dat(gwt)
    dat[_GW_TASK_EXPR_KEY] = ex
    dat[_GW_TASK_SRC_KEY] = src
    dat[_GW_TASK_DESC_KEY] = desc
    dat[_GW_TASK_LANG_KEY] = lang
    serialize(dfile, dat)

    return dat
end

function _read_task_dat!(gwt::GWTask)

    tfile = _taskdat_file(gwt)
    dat = deserialize(tfile)
    merge!(_task_dat(gwt), dat)

    return dat
end

function _task_desc(gwt::GWTask)
    dat = _task_dat(gwt)
    get!(dat, _GW_TASK_DESC_KEY, "")
end

function _task_expr(gwt::GWTask)
    dat = _task_dat(gwt)
    get!(dat, _GW_TASK_EXPR_KEY, nothing)
end

function _task_src(gwt::GWTask)
    dat = _task_dat(gwt)
    get!(() -> _expr_src(gwt), dat, _GW_TASK_SRC_KEY)
end

function _task_lang(gwt::GWTask)
    dat = _task_dat(gwt)
    get!(dat, _GW_TASK_LANG_KEY, "")
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
module Analysis

ENV["GKSwstype"]="nul"
using JuliaDB,GLM,StatsPlots,Dates

function exp_regress(cases)
    aggregate_cases = groupreduce((Confirmed=+,), cases, :ObservationDate, select=:Confirmed)
    case_data = select(aggregate_cases,
                       (:ObservationDate => d -> (d - Date(2020,1,1)).value, :Confirmed => n -> log(2,n)))
    regression = lm(@formula(Confirmed ~ 1 + ObservationDate), case_data)
    (regression,
     plot(select(aggregate_cases,:ObservationDate), [exp2.(select(case_data,:Confirmed)),exp2.(predict(regression,case_data))], legend=:bottomright, label=["data" "regression"],yaxis=:log))
end

function filter_date(s,e,cases)
    filter(:ObservationDate => d -> s <= d <= e,cases)
end

function filter_country(c,cases)
    filter(Symbol("Country/Region") => r -> r == c,cases)
end

function filter_state(postal,fullname,cases)
    filter(Symbol("Province/State") => ps -> length(ps) >= 2 && (occursin(postal,ps) || ps == fullname), cases)
end

function last_n_days(n,cases)
    filter_date(today() - Day(n), today(),cases)
end

function process_day_file(path,day)
    fullpath = path * "/csse_covid_19_data/csse_covid_19_daily_reports/" * day * ".csv"
    raw = loadtable(fullpath)
    date = Date(day,"m-d-y")
    insertcols(select(raw, (:Confirmed,:Deaths,:Recovered,Symbol("Country/Region"),Symbol("Province/State"))), 1, :ObservationDate => repeat([date],length(raw)))
end

function load_csse_data(path)
    fullpath = path * "/csse_covid_19_data/csse_covid_19_daily_reports/"
    function parse_date(dt)
        try
            Date(dt)
        catch
            Date(2020,01,01)
        end
    end
    report_data = [process_day_file(fullpath * file) for file in filter(f -> f[end-2:end] == "csv",readdir(fullpath))]
    dropmissing(reduce(merge, report_data))
end

function add_new_day(path,day,t)
    new_data = process_day_file(path)
    merge(new_data,t)
end

end


text:
"""
    fire!(actr)

Selects a production rule and registers conflict resolution and new events for selected production rule

# Arguments

- `actr`: an ACT-R model object
"""
function fire!(actr::AbstractACTR)
    actr.procedural.state.busy ? (return nothing) : nothing
    rules = select_rule(actr)
    isempty(rules) ? (return nothing) : nothing
    rule = rules[1]
    description = "Selected "*rule.name
    type = "model"
    id = get_name(actr)
    tΔ = rnd_time(.05)
    resolving(actr, true)
    f(r, a, v) = (resolving(a, v), r.action())
    register!(actr, f, after, tΔ, rule, actr, false; description, type, id)
    return nothing
end

resolving(actr, v) = actr.procedural.state.busy = v

function select_rule(actr)
    rules = get_rule_set(actr)
    isempty(rules) ? (return rules) : nothing
    compute_utilities!(actr, rules)
    _,idx = findmax(x -> x.utility, rules)
    return rules[idx:idx]
end

function get_rule_set(actr)
    @unpack mmp_utility = actr.parms
    rules = get_rules(actr)
    if mmp_utility
        # use @code_string to filter out conditions that do not contain
        # slots. Consider caching a boolean that indicates whether a production
        # rule is eligible for partial matching
        return rules
    end
    return filter(r->match(actr, r), rules)
end

get_rules(actr) = actr.procedural.rules

function match(actr, rule)
    return all_match(actr, rule.conditions)
end

"""
    all_match(actr, conditions)

Checks whether all conditions of a production rule are satisfied.

# Arguments

- `actr`: an ACT-R model object
- `conditions`: a tuple of functions representing production rule conditions.
"""
function all_match(actr, conditions)
    for c in conditions
        !c(actr) ? (return false) : nothing
    end
    return true
end

function compute_utilities!(actr, rules)
    for rule in rules
        compute_utility!(actr, rule)
    end
    return nothing
end

function compute_utility!(actr, rule)
    @unpack utility_noise, mmp_utility = actr.parms
    mmp_utility ? compute_penalties!(actr, rule) : nothing
    utility_noise ? add_noise!(actr, rule) : nothing
    total_utility!(rule)
end

function compute_penalties!(actr, rule)
    @unpack model_trace = actr.scheduler
    @unpack δu,util_mmp_fun = actr.parms
    model_trace ? println("rule: ", rule.name) : nothing
    penalty = 0.0
    for c in rule.conditions
        penalty += util_mmp_fun(actr, c)
    end
    model_trace ? println("") : nothing
    rule.utility_penalty = penalty * δu
    return nothing
end

function utility_match(actr::ACTR, condition)
    @unpack model_trace = actr.scheduler
    penalty = condition(actr) ? 0.0 : 1.0
    model_trace ? condition_trace(actr, condition, penalty) : nothing
    return penalty
end

function condition_trace(actr, condition, penalty)
    str = @code_string condition(actr)
    println(str)
    println("penalty: $penalty")
    return nothing
end

function add_noise!(actr, rule)
    @unpack σu = actr.parms
    rule.utility_noise = rand(Normal(0, σu))
end

function total_utility!(r::Rule)
    r.utility_mean = r.initial_utility + r.utility_penalty
    r.utility = r.utility_mean + r.utility_noise
    return nothing
end

function Rule(;
    utility = 0.0,
    initial_utility = 0.0,
    utility_mean = 0.0,
    utility_penalty = 0.0,
    utility_noise = 0.0,
    conditions,
    name = "",
    actr,
    task,
    action,
    args = (),
    kwargs...
    )
    funs = conditions(actr, args...; kwargs...)
    can_pm = can_partial_match(funs, actr)
    Rule(
        utility,
        initial_utility,
        utility_mean,
        utility_penalty,
        utility_noise,
        conditions(actr, args...; kwargs...),
        ()->action(actr, task; kwargs...),
        can_pm,
        name
     )
end

function can_partial_match(conditions, actr)
    for c in conditions
        str = @code_string c(actr)
        !occursin("slots", str) ? (return false) : nothing
    end
    return true
end