text:
# Spawn a thread (using pthreads on Unix/OSX and Windows threads on Windows)
# to implement the "heartbeat" message channel, which is just a ZMQ
# socket that echoes all messages.  This is implemented with the zmq_device
# call in libzmq, which simply blocks forever, so the usual lack of
# thread safety in Julia should not be an issue here.

# entry point for new thread
function heartbeat_thread(sock::Ptr{Void})
    ccall((:zmq_device,ZMQ.zmq), Cint, (Cint, Ptr{Void}, Ptr{Void}),
          2, sock, sock)
    nothing # not correct on Windows, but irrelevant since we never return
end
const heartbeat_c = cfunction(heartbeat_thread, Void, (Ptr{Void},))

if @windows? false : true
    const threadid = Array(Int, 128) # sizeof(pthread_t) is <= 8 on Linux & OSX
end

function start_heartbeat(sock)
    @windows? begin
        ccall(:_beginthread, Int, (Ptr{Void}, Cuint, Ptr{Void}),
              heartbeat_c, 0, sock.data)
    end : begin
        ccall((:pthread_create, :libpthread), Cint,
              (Ptr{Int}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
              threadid, C_NULL, heartbeat_c, sock.data)
    end
end


text:
export DynamicMatLawLoss, preprocessing, LSfittingStress, LSfittingStressHelper

@doc """
    - 'domain': Domain, finite element domain, for data structure
    - 'E_all':  Float64[NT+1, neles*nGauss, nstrains], all strains for the whole simulation from training data
    - 'w∂E∂u_all': Float64[NT+1, neles*nGauss, ndofs_per_element, nstrains],
                multiplication of the Gaussian weight and ∂E∂u^T for the whole simulation from training data
    - 'F_tot': Float64[NT, neqs] approximated internal force for the whole simulation, from time n=1 to time n=NT

    form the loss function, based on dynamic equilibrium
        (Mddu + fint(NN, E, DE) + MIDddu_bc = fext

    we have
        F_tot  = fext - MIDddu_bc - Mddu
    loss = ∑ ||fint(NN, E, DE) - F_tot||^2

    todo
    We assume the time step is constant for the training data
"""->
function DynamicMatLawLoss(domain::Domain, E_all::Array{Float64}, w∂E∂u_all::Array{Float64},
     F_tot::Array{Float64}, nn::Function)
    # you can use fint instead of computed F_tot to debug
    # F_tot =  hcat(domain.history["fint"]...)'
    # define variables
    neles = domain.neles
    nGauss = length(domain.elements[1].weights)
    nstrains = size(E_all,3)

    NT = size(E_all,1)-1
    @assert size(E_all)==(NT+1, neles*nGauss, nstrains)
    @assert size(F_tot)==(NT, domain.neqs)

    #convert to tensors (Tensorflow variables)
    E_all = constant(E_all)
    F_tot = constant(F_tot)
    w∂E∂u_all = constant(w∂E∂u_all)

    #instead using for loop for time step update, while loop is used for efficiency (specifically for Tensorflow)
    #time step is 1, 2, 3, ....NT+1, 1 is the initial condition
    function cond0(i, ta_loss, ta_σ)
        i<=NT+1
    end

    function body(i, ta_loss, ta_σ)
        E = E_all[i]
        DE = E_all[i-1]
        w∂E∂u = w∂E∂u_all[i]
        σ0 = read(ta_σ, i-1)
        fint, σ = tfAssembleInternalForce(domain,nn,E,DE,w∂E∂u,σ0)
        ta_σ = write(ta_σ, i, σ)
        ta_loss = write(ta_loss, i, sum((fint-F_tot[i-1])^2))
        i+1, ta_loss, ta_σ
    end

    σ0 = constant(zeros(neles*nGauss, nstrains))
    # tensorflow stress collection at each time step
    ta_σ = TensorArray(NT+1); ta_σ = write(ta_σ, 1, σ0)
    # tensorflow loss collection at each time step
    ta_loss = TensorArray(NT+1); ta_loss = write(ta_loss, 1, constant(0.0))
    i = constant(2, dtype=Int32)
    _, out, _ = while_loop(cond0, body, [i,ta_loss, ta_σ]; parallel_iterations=20)

    total_loss = sum(stack(out)[2:NT])
    return total_loss
end

# function DynamicMatLawLossWithTailLoss(domain::Domain, E_all::Array{Float64}, w∂E∂u_all::Array{Float64},
#     F_tot::Array{Float64}, nn::Function, H0::Array{Float64}, n_tail::Int64, stress_scale::Float64)
#    # todo, use fint instead of computed F_tot
#    # F_tot =  hcat(domain.history["fint"]...)'
#    # define variables
#    neles = domain.neles
#    nGauss = length(domain.elements[1].weights)
#    nstrains = size(E_all,3)

#    NT = size(E_all,1)-1
#    @assert size(E_all)==(NT+1, neles*nGauss, nstrains)
#    @assert size(F_tot)==(NT, domain.neqs)
#    E_all = constant(E_all)
#    F_tot = constant(F_tot)
#    w∂E∂u_all = constant(w∂E∂u_all)

#    function cond0(i, ta_loss, ta_σ, ta_tail)
#        i<=NT+1
#    end

#    function body(i, ta_loss, ta_σ, ta_tail)
#        E = E_all[i]
#        DE = E_all[i-1]
#        w∂E∂u = w∂E∂u_all[i]
#        σ0 = read(ta_σ, i-1)
#        fint, σ = tfAssembleInternalForce(domain,nn,E,DE,w∂E∂u,σ0)
#        ta_σ = write(ta_σ, i, σ)
#        ta_loss = write(ta_loss, i, sum((fint-F_tot[i-1])^2))
#        tail_loss = tf.cond(i<=NT+1-n_tail,
#             ()->constant(0.0),
#             ()->begin
#                 σ_all = nn(E, DE, σ0)
#                 norm( (E-DE)*H0 + σ0/stress_scale - σ_all/stress_scale )
#             end
#        )
#        ta_tail = write(ta_tail, i, tail_loss)
#        i+1, ta_loss, ta_σ, ta_tail
#    end

#    σ0 = constant(zeros(neles*nGauss, nstrains))
#    ta_σ = TensorArray(NT+1); ta_σ = write(ta_σ, 1, σ0)
#    ta_tail = TensorArray(NT+1); ta_tail = write(ta_tail, 1, constant(0.0))
#    ta_loss = TensorArray(NT+1); ta_loss = write(ta_loss, 1, constant(0.0))
#    i = constant(2, dtype=Int32)
#    _, out, _, tout = while_loop(cond0, body, [i,ta_loss, ta_σ, ta_tail]; parallel_iterations=20)

#    total_loss = sum(stack(out)[2:NT])
#    tail_loss = sum(stack(tout)[2:NT])
#    return total_loss, tail_loss
# end

# function DynamicMatLawLossInternalVariable(domain::Domain, E_all::Array{Float64},
#     w∂E∂u_all::Array{Float64}, F_tot::Array{Float64}, nn::Function, n_internal::Int64=128)
#    # todo, use fint instead of computed F_tot
#    # F_tot =  hcat(domain.history["fint"]...)'
#    # define variables
#    neles = domain.neles
#    nGauss = length(domain.elements[1].weights)
#    nstrains = size(E_all,3)

#    NT = size(E_all,1)-1
#    @assert size(E_all)==(NT+1, neles*nGauss, nstrains)
#    @assert size(F_tot)==(NT, domain.neqs)
#    E_all = constant(E_all)
#    F_tot = constant(F_tot)
#    w∂E∂u_all = constant(w∂E∂u_all)

#    function cond0(i, ta_loss, ta_σ, ta_α)
#        i<=NT+1
#    end

#    function body(i, ta_loss, ta_σ, ta_α)
#        E = E_all[i]
#        DE = E_all[i-1]
#        w∂E∂u = w∂E∂u_all[i]
#        σ0 = read(ta_σ, i-1)
#        α0 = read(ta_α, i-1)
#        fint, σ, α = tfAssembleInternalForce(domain,nn,E,DE,w∂E∂u,σ0, α0)
#        ta_σ = write(ta_σ, i, σ)
#        ta_loss = write(ta_loss, i, sum((fint-F_tot[i-1])^2))
#        ta_α = write(ta_α, i, α)
#        i+1, ta_loss, ta_σ, ta_α
#    end

#    σ0 = constant(zeros(neles*nGauss, nstrains))
#    ta_σ = TensorArray(NT+1); ta_σ = write(ta_σ, 1, σ0)
#    ta_loss = TensorArray(NT+1); ta_loss = write(ta_loss, 1, constant(0.0))
#    ta_α = TensorArray(NT+1);  ta_α = write(ta_α, 1, constant(zeros(neles*nGauss, n_internal)))
#    i = constant(2, dtype=Int32)
#    _, out, _ = while_loop(cond0, body, [i,ta_loss, ta_σ, ta_α]; parallel_iterations=20)

#    total_loss = sum(stack(out)[2:NT])
#    return total_loss
# end

@doc """
    - 'domain': Domain,  finite element domain
    - 'globdat': GlobalData,  finite element data structure
    - 'state_history': Float64[NT+1][ndofs], displace history of all time steps and all nodes, including time 0
                    hcat(state_history...) gives a matrix of Float64[NT+1, ndofs]
    - 'fext_history': Float64[NT+1][neqns] or Float64[NT][neqns]: external force load of all time steps and all equations,
                    with or without the initial step
                    hcat(fext_history...) gives a matrix of size(neqns, NT+1)
    nn: Neural network
    Δt: time step size

    compute loss function from state and external force history
"""->
function DynamicMatLawLoss(domain::Domain, globdat::GlobalData, state_history::Array{T}, fext_history::Array{S}, nn::Function, Δt::Float64) where {T, S}
    # todo convert to E_all, Ftot
    domain.history["state"] = state_history
    F_tot, E_all, w∂E∂u_all = preprocessing(domain, globdat, hcat(fext_history...), Δt)
    DynamicMatLawLoss(domain, E_all, w∂E∂u_all, F_tot, nn)
end

# function DynamicMatLawLoss(domain::Domain, globdat::GlobalData, state_history::Array{T},
#     fext_history::Array{S}, nn::Function, Δt::Float64, H0::Array{Float64}, n_tail::Int64, stress_scale::Float64) where {T, S}
#     # todo convert to E_all, Ftot
#     domain.history["state"] = state_history
#     F_tot, E_all, w∂E∂u_all = preprocessing(domain, globdat, hcat(fext_history...), Δt)
#     DynamicMatLawLossWithTailLoss(domain, E_all, w∂E∂u_all, F_tot, nn, H0, n_tail, stress_scale)
# end

# function DynamicMatLawLoss(domain::Domain, globdat::GlobalData, state_history::Array{T}, fext_history::Array{S}, nn::Function, Δt::Float64, n::Int64) where {T, S}
#     domain.history["state"] = state_history
#     F_tot, E_all, w∂E∂u_all = preprocessing(domain, globdat, hcat(fext_history...), Δt, n)
#     DynamicMatLawLoss(domain, E_all, w∂E∂u_all, F_tot, nn)
# end

# function DynamicMatLawLossInternalVariable(domain::Domain, globdat::GlobalData, state_history::Array{T}, fext_history::Array{S}, nn::Function, Δt::Float64, n_internal::Int64) where {T, S}
#     domain.history["state"] = state_history
#     F_tot, E_all, w∂E∂u_all = preprocessing(domain, globdat, hcat(fext_history...), Δt)
#     DynamicMatLawLossInternalVariable(domain, E_all, w∂E∂u_all, F_tot, nn, n_internal)
# end

@doc """
    Postprocess the data from domain.history and F_ext
    - 'domain': Domain, domain.history["state"] is Float64[NT+1][ndofs]
    - 'globdat': GlobalData, for mass matrix
    - 'F_ext': Float64[neqns, NT+1] or Float64[neqns, NT], external force history with/without the initial step
    - 'Δt': Float64, time step size

    Return:
    - 'F_tot': Float64[NT, neqns], approxmiated internal force
    - 'E_all': Float64[NT+1, neles*nGauss, nstrains], all strains for the whole simulation from training data
    - 'w∂E∂u_all': Float64[NT+1, neles*nGauss, neqns_per_elem, nstrains], multiplication of the Gaussian weight and ∂E∂u^T
    for the whole simulation from training data

    todo
    We assume the time step is constant for the training data
"""->
function preprocessing(domain::Domain, globdat::GlobalData, F_ext::Array{Float64},Δt::Float64)
    U = hcat(domain.history["state"]...)

    M = globdat.M
    MID = globdat.MID

    NT = size(U,2)-1

    #Acceleration of the time-dependent Dirichlet boundary nodes
    bc_acc = zeros(sum(domain.EBC.==-2),NT)
    if !(globdat.EBC_func===nothing)
        for i = 1:NT
            _, bc_acc[:,i]  = globdat.EBC_func(Δt*i)
        end
    end

    ∂∂U = zeros(size(U,1), NT+1)
    ∂∂U[:,2:NT] = (U[:,1:NT-1]+U[:,3:NT+1]-2U[:,2:NT])/Δt^2

    if size(F_ext,2)==NT+1
        F_tot = F_ext[:,2:end] - M*∂∂U[domain.dof_to_eq,2:end] - MID*bc_acc
    elseif size(F_ext,2)==NT

        F_tot = F_ext - M*∂∂U[domain.dof_to_eq,2:end] - MID*bc_acc
    else
        error("F size is not valid")
    end

    neles = domain.neles
    nGauss = length(domain.elements[1].weights)
    neqns_per_elem = length(getEqns(domain,1))

    nstrains = div((domain.elements[1].eledim + 1)*domain.elements[1].eledim, 2)

    E_all = zeros(NT+1, neles*nGauss, nstrains)
    w∂E∂u_all = zeros(NT+1, neles*nGauss, neqns_per_elem, nstrains)

    for i = 1:NT+1
        domain.state = U[:, i]
        # Loop over the elements in the elementGroup to construct strain and geo-matrix E_all and w∂E∂u_all
        for iele  = 1:neles
            element = domain.elements[iele]

            # Get the element nodes
            el_nodes = getNodes(element)

            # Get the element nodes
            el_eqns = getEqns(domain,iele)

            el_dofs = getDofs(domain,iele)

            el_state  = getState(domain,el_dofs)

            # Get the element contribution by calling the specified action
            E, w∂E∂u = getStrain(element, el_state)

            E_all[i, (iele-1)*nGauss+1:iele*nGauss, :] = E

            w∂E∂u_all[i, (iele-1)*nGauss+1:iele*nGauss,:,:] = w∂E∂u
        end
    end
    @info "preprocessing end..."
    return F_tot'|>Array, E_all, w∂E∂u_all
end

# @doc """
# preprocessing(domain::Domain, globdat::GlobalData, F_ext::Array{Float64},Δt::Float64, n::Int64)

# Same as `preprocessing`, except that only the first `n` steps are considered
# """->
# function preprocessing(domain::Domain, globdat::GlobalData, F_ext::Array{Float64},Δt::Float64, n::Int64)
#     U = hcat(domain.history["state"]...)
#     # @info " U ", size(U),  U'
#     M = globdat.M
#     MID = globdat.MID

#     NT = size(U,2)-1

#     #Acceleration of Dirichlet nodes
#     bc_acc = zeros(sum(domain.EBC.==-2),NT)
#     if !(globdat.EBC_func===nothing)
#         for i = 1:NT
#             _, bc_acc[:,i]  = globdat.EBC_func(Δt*i)
#         end
#     end

#     ∂∂U = zeros(size(U,1), NT+1)
#     ∂∂U[:,2:NT] = (U[:,1:NT-1]+U[:,3:NT+1]-2U[:,2:NT])/Δt^2

#     if size(F_ext,2)==NT+1
#         F_tot = F_ext[:,2:end] - M*∂∂U[domain.dof_to_eq,2:end] - MID*bc_acc
#     elseif size(F_ext,2)==NT
#         F_tot = F_ext - M*∂∂U[domain.dof_to_eq,2:end] - MID*bc_acc
#     else
#         error("F size is not valid")
#     end

#     neles = domain.neles
#     nGauss = length(domain.elements[1].weights)
#     neqns_per_elem = length(getEqns(domain,1))

#     nstrains = div((domain.elements[1].eledim + 1)*domain.elements[1].eledim, 2)

#     E_all = zeros(NT+1, neles*nGauss, nstrains)
#     w∂E∂u_all = zeros(NT+1, neles*nGauss, neqns_per_elem, nstrains)

#     for i = 1:NT+1
#         domain.state = U[:, i]
#         # @info "domain state", domain.state
#         # Loop over the elements in the elementGroup to construct strain and geo-matrix E_all and w∂E∂u_all
#         for iele  = 1:neles
#             element = domain.elements[iele]

#             # Get the element nodes
#             el_nodes = getNodes(element)

#             # Get the element nodes
#             el_eqns = getEqns(domain,iele)

#             el_dofs = getDofs(domain,iele)

#             el_state  = getState(domain,el_dofs)

#             # Get the element contribution by calling the specified action
#             E, w∂E∂u = getStrain(element, el_state)
#             # if i==2
#             #     @info (el_state, E)
#             # end

#             # @show E, nGauss
#             E_all[i, (iele-1)*nGauss+1:iele*nGauss, :] = E

#             w∂E∂u_all[i, (iele-1)*nGauss+1:iele*nGauss,:,:] = w∂E∂u
#         end
#     end
#     @info "preprocessing end..."
#     # # DEBUG
#     # fext = hcat(domain.history["fint"]...)
#     F_tot = F_tot'|>Array

#     F_tot[1:n,:], E_all[1:n+1, :, :], w∂E∂u_all[1:n+1, :, :, :]
# end

##########################################################
@doc """

Compute the Fint, dFint_dS_comp, E_all, S_all and dfint_dstress, based on the state in domain and the stress in S_comp
- 'domain': Domain
- 'S_comp': Float64[neles, nstrain] or Float64[npoints, nstrain], the stresses in each element(Constant) or at each geometric point(Linear)
- 'method': String "Constant" or "Linear" for Constant approach or Linear approach

Return
    Fint: Float64[neqns], constructed from state in domain
    dFint_dS_comp: Float64[neqs,  size(S_comp,1)*nstrain], sparse matrix representation
    E_all: Float64[neles*ngps_per_elem, nstrain], strains at each Gaussian quadrature point ,constructed from state in domain
    S_all: Float64[neles*ngps_per_elem, nstrain], stresses at each Gaussian quadrature point ,constructed from S_comp
"""->
function LSfittingStressHelper(domain, S_comp::Array{Float64}, method::String)
    neles = domain.neles
    nnodes = domain.nnodes
    npoints = domain.npoints
    eledim = domain.elements[1].eledim
    nstrain = div((eledim + 1)*eledim, 2)
    ngps_per_elem = length(domain.elements[1].weights)
    neqs = domain.neqs
    eledim = domain.elements[1].eledim
    elepoints = 2^eledim  #number of anchor points of an element

    node_to_point = domain.node_to_point

    #check
    @assert(method=="Constant" || method=="Linear")
    method=="Constant" ? @assert(size(S_comp,1) == neles) : @assert(size(S_comp,1) == npoints)

    Fint = zeros(Float64, neqs)
    ii = Int64[]; jj=Int64[]; vv=Float64[]
    E_all = zeros(neles*ngps_per_elem, nstrain)
    S_all = similar(E_all)

    # each element involves ele_per_S_comp number of compact Stress points,
    # for Constant the stress is assume to be constant ele_per_S_comp = 1
    # for Linear the stress is assume to be bilinear ele_per_S_comp = 4
    # ele_per_S_comp = method=="Constant" ? 1 : 4;

    ele_per_S_comp = method=="Constant" ? 1 : elepoints;

    # stress array at each Gaussian point
    S_ele = zeros(Float64, ngps_per_elem, nstrain)
    # dstress_dcompact_Stress,
    # for Constant, compact stress contains 1 in each element
    # for Linear, compact stress contains 4 in each element
    dS_ele_all_dS_ele_comp =  zeros(Float64, ngps_per_elem*nstrain, ele_per_S_comp * nstrain)

    # Loop over the elements in the elementGroup, to construct S_ele, and dS_ele_all_dS_ele_comp
    for iele  = 1:neles
      element = domain.elements[iele]

      # Get the element nodes
      el_eqns = getEqns(domain,iele)

      el_dofs = getDofs(domain,iele)

      el_state  = getState(domain, el_dofs)

      el_coords_comp = element.coords[1:elepoints, :]

      #the node number is the element point number
      el_nodes_comp = node_to_point[element.elnodes[1:elepoints]]

        if method == "Constant"
            dS_ele_all_dS_ele_comp[:,:]  .= 0.0
            for igp = 1:ngps_per_elem
                S_ele[igp, :] = S_comp[iele, :]

                for i = 1:nstrain
                    dS_ele_all_dS_ele_comp[(igp-1)*nstrain + i, i] = 1.0
                end
            end

        else
            @assert(method == "Linear")
            @assert(eledim <= 2)
            ngp = Int64((ngps_per_elem)^(1.0/eledim))

            hs = (eledim==1) ? get1DElemShapeData(el_coords_comp , ngp)[end] : get2DElemShapeData(el_coords_comp , ngp)[end]

            dS_ele_all_dS_ele_comp[:,:]  .= 0.0
            @assert(length(hs[1]) == ele_per_S_comp)

            for igp = 1:ngps_per_elem
                S_ele[igp, :] = S_comp[el_nodes_comp, :]' * hs[igp]

                for i = 1:nstrain
                    for j = 1:ele_per_S_comp
                        dS_ele_all_dS_ele_comp[(igp-1)*nstrain + i, (j-1)*nstrain+i] = hs[igp][j]
                    end
                end
            end
        end

      # Get the element contribution by calling the specified action
      fint, dfint_dS_ele_all = getForceAndDforceDstress(element, el_state, S_ele)
      dfint_dS_ele_comp = dfint_dS_ele_all * dS_ele_all_dS_ele_comp

      el_eqns_active = el_eqns .>= 1
      el_eqns_active_idx = el_eqns[el_eqns_active]
      Fint[el_eqns_active_idx] += fint[el_eqns_active]

      # K[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += stiff[el_eqns_active,el_eqns_active]
      dfint_dS_ele_comp_active = dfint_dS_ele_comp[el_eqns_active,:]

      for i = 1:length(el_eqns_active_idx)
        for j = 1:ele_per_S_comp
            for k = 1:nstrain
                push!(ii, el_eqns_active_idx[i])
                if method == "Constant"
                    push!(jj, (iele-1)*nstrain + k)
                else
                    push!(jj, (el_nodes_comp[j] - 1)*nstrain + k)
                end

                push!(vv, dfint_dS_ele_comp_active[i,(j-1)*nstrain + k])
            end
        end
      end

      E_ele, _ = getStrain(element, el_state)

      E_all[(iele-1)*ngps_per_elem+1:iele*ngps_per_elem, :] = E_ele
      S_all[(iele-1)*ngps_per_elem+1:iele*ngps_per_elem, :] = S_ele

    end

    dFint_dS_comp =  sparse(ii, jj, vv, neqs,  size(S_comp,1)*nstrain)

    return Fint, dFint_dS_comp, E_all, S_all
  end

@doc """
    Fit the stress from the strain and external load conditions,
    For quadratic element, nx, ny, ... elements in each direction
    The number of equations are neqs ≈ 2*(2nx + 1)*(2ny + 1)
    Linear approach : Assume the stress in each element is constant, the number of unknows are 3*nx*ny
    Constant approach: Assume the stress in each element is linear, the number of unknows are 3*(nx+1)*(ny+1) on each nodes
    Least Square problem is solving for stress

    - 'domain': Domain
    - 'globdat': GlobalData
    - 'state_history': Float64[NT+1, ndofs], displacement fields for each freedoms
    - 'F_ext': Float64[neqns, NT+1] or Float64[neqns, NT+1], external force vector for each equations
    - 'Δt': Float64, time step size
    - 'method': String "Constant" or "Linear" for Constant approach or Linear approach

    Return:
    - 'E_all':  Float64[NT+1, neles*ngps_per_elem, nstrain], strain fields
    - 'S_all':  Float64[NT+1, neles*ngps_per_elem, nstrain], fitted stress fields at each Gaussian quadratic point

    Remark: Turss 1D requires method == "Constant", otherwise the system is under constraint
"""->
function LSfittingStress(domain::Domain, globdat::GlobalData, state_history::Array{Float64}, F_ext::Array{Float64},Δt::Float64, method::String)
    nnodes = domain.nnodes
    neles = domain.neles
    eledim = domain.elements[1].eledim
    nstrain = div((eledim + 1)*eledim, 2)

    neles = domain.neles
    ngps_per_elem = length(domain.elements[1].weights)
    neqs = domain.neqs

    U = state_history
    M = globdat.M
    MID = globdat.MID

    NT = size(U,2)-1

    #Acceleration of Dirichlet nodes
    bc_acc = zeros(sum(domain.EBC.==-2),NT)
    if !(globdat.EBC_func===nothing)
        for i = 1:NT
            _, bc_acc[:,i]  = globdat.EBC_func(Δt*i)
        end
    end

    ∂∂U = zeros(size(U,1), NT+1)
    ∂∂U[:,2:NT] = (U[:,1:NT-1]+U[:,3:NT+1]-2U[:,2:NT])/Δt^2
    ∂∂U[:,NT+1] = (U[:,NT+1]-U[:,NT])/Δt

    @show "compute estimated F_int"
    if size(F_ext,2)==NT+1
        F_tot = F_ext[:,2:end] - M*∂∂U[domain.dof_to_eq,2:end] - MID*bc_acc
    elseif size(F_ext,2)==NT
        F_tot = F_ext - M*∂∂U[domain.dof_to_eq,2:end] - MID*bc_acc
    else
        error("F size is not valid", size(F_ext,2) , NT)
    end

    E_all = zeros(NT+1, neles*ngps_per_elem, nstrain)
    S_all = zeros(NT+1, neles*ngps_per_elem, nstrain)

    @assert(domain.npoints > 0)
    nS_comp = (method == "Constant") ? neles : domain.npoints;

    @assert(nS_comp*nstrain < neqs)

    S_comp_all = zeros(Float64, NT+1, nS_comp, nstrain)

    @assert(nS_comp*nstrain <= domain.neqs)

    for it = 2:NT+1

        domain.state = U[:, it-1]

        Fint, dFint_dS_comp, E, S = LSfittingStressHelper(domain, S_comp_all[it, :, :], method)

        #min ||F_tot -  F_int(E, S)||^2 = min ||F_tot -  dF_int/dS_comp * S_comp||^2

        S_comp = dFint_dS_comp\F_tot[:,it-1]

        S_comp = reshape(S_comp, nstrain, nS_comp)'|>Array
        Fint, dFint_dS_comp, E, S = LSfittingStressHelper(domain, S_comp, method)

        @show "LS loss is ", norm(Fint - F_tot[:,it-1])

        E_all[it, :, :] = E
        S_all[it, :, :] = S

        S_comp_all[it,:,:] = S_comp

    end

    E_all, S_all
end

ext: jl


text:
#=

  Archery Puzzle in Julia ConstraintSolver.jl

  http://www.eonhq.com/m/puzzles/images/archery-puzzle.jpg
  Archery puzzle by Sam Loyd:
  """
  How close can the young archer come to scoring a total of
  100 - using as many arrows as she please.
  [The targets are: 16, 17, 23, 24, 39, 40]
  """
  Via: The Aperiodical: "Manchester MathsJam June 2012 Recap"
  http://aperiodical.com/2012/06/manchester-mathsjam-june-2012-recap/

  Solution (unique):
  x = [2,4,0,0,0,0]
  z = 100

  i.e. 2 arrows on score 16 and 4 arrows on score 17:
    2*16 + 4*17 = 100

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function archery_puzzle(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit" => 6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 6
    targets =  [16, 17, 23, 24, 39, 40]
    target = 100
    @variable(model, 0 <= x[1:n] <= 10, Int)
    @variable(model, 0 <= z <= n*maximum(targets), Int)
    @variable(model, 0 <= diff_z <= n*maximum(targets), Int)

    scalar_product(model,x,targets,:(==),z)
    my_abs(model,target,z,diff_z)

    @objective(model,Min,diff_z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                diff_z_val = convert.(Integer,JuMP.value.(diff_z; result=sol))
                println("x:$x_val")
                println("z:$z_val diff_z:$diff_z_val")
                println(["target:$(targets[i]) x:$(x_val[i])  " for i in 1:n if x_val[i] > 0]|>join)
                println()

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time archery_puzzle(true,true)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Iterating
# ---------

mutable struct SegmentBTreeIteratorState{K, V, B}
    leaf::Union{Nothing, LeafNode{K, V, B}}
    i::Int
end

function Base.iterate(
        t::SegmentBTree{K, V, B},
        state::SegmentBTreeIteratorState{K, V, B}=iterinitstate(t)) where {K, V, B}
    if state.leaf === nothing || isempty(state.leaf)
        return nothing
    end
    leaf = notnothing(state.leaf)
    entry = leaf.entries[state.i]
    if state.i < length(leaf)
        state.i += 1
    else
        state.leaf = leaf.right
        state.i = 1
    end
    return entry, state
end

function iterinitstate(t::SegmentBTree{K,V,B}) where {K, V, B}
    # traverse to the first leaf node
    node = t.root
    while !isa(node, LeafNode{K, V, B})
        node = node.children[1]
    end
    return SegmentBTreeIteratorState(node, 1)
end

# Iterate from a given starting from the first Segment that intersects a given point.
struct SegmentFromIterator{K, V, B}
    t::SegmentBTree{K, V, B}
    p::K
end

Base.eltype(::Type{SegmentFromIterator{K, V, B}}) where {K, V, B} = V
Base.IteratorSize(::Type{SegmentFromIterator{K, V, B}}) where {K, V, B} = Base.SizeUnknown()

function from(t::SegmentBTree{K, V, B}, p) where {K, V, B}
    return SegmentFromIterator{K, V, B}(t, convert(K, p))
end

function Base.iterate(
        it::SegmentFromIterator{K, V, B},
        state::SegmentBTreeIteratorState{K, V, B}=iterinitstate(it)) where {K, V, B}
    return iterate(it.t, state)
end

function iterinitstate(it::SegmentFromIterator{K, V, B}) where {K, V, B}
    node, i = firstfrom(it.t, it.p)
    if i == 0
        return SegmentBTreeIteratorState{K, V, B}(nothing, i)
    else
        return SegmentBTreeIteratorState{K, V, B}(node, i)
    end
end


text:
# This is the backbone of the VMC code.  Both determinants and Pfaffians are available.

abstract VMCMatrix{T}

#####################################################################
#################### Determinantal wavefunctions ####################
#####################################################################

type CeperleyMatrix{T} <: VMCMatrix{T}
    # FIXME: make current_state an Enum, i.e., enforce that it's one of a certain set of states
    current_state::Symbol

    detmat::Matrix{T}
    inverse::Matrix{T}

    det::T
    old_det::T

    detrat::T
    Q_inv::Matrix{T}

    # pending indices
    cs::Vector{Int}
    rs::Vector{Int}

    # these apply equally well to single column or row updates
    pending_cols::Matrix{T}
    pending_rows::Matrix{T}

    # these are only necessary for combined column and row updates,
    # and are determined from pending_cols and pending_rows
    cols_tilde::Matrix{T}
    rows_tilde::Matrix{T}

    CeperleyMatrix() = new(:UNINITIALIZED)

    function CeperleyMatrix(mat::Matrix{T})
        @assert size(mat)[1] == size(mat)[2]
        LUmat = lufact(mat)
        detval = det(LUmat)
        return new(:READY_FOR_UPDATE, mat, inv(LUmat), detval, detval)
    end
end

CeperleyMatrix(mat::Matrix) = CeperleyMatrix{eltype(mat)}(mat)
CeperleyMatrix() = CeperleyMatrix{Float64}() # default to Float64

# this should be unnecessary to call anywhere below,
# unless there was some rogue action on the detmat and inverse fields from the outside
function consistent_Ceperley_dims(cmat::CeperleyMatrix)
    return size(cmat.detmat)[1] == size(cmat.detmat)[2] &&
           size(cmat.inverse)[1] == size(cmat.inverse)[2] &&
           size(cmat.detmat)[1] == size(cmat.inverse)[1]
end

# check inverse for numerical error due to SMW
function check_inverse{T}(cmat::CeperleyMatrix{T}, tol::Float64 = 1e-10)
    @assert cmat.current_state == :READY_FOR_UPDATE

    error = sum(abs(cmat.detmat * cmat.inverse - eye(T, size(cmat.detmat)[1]))) / size(cmat.detmat)[1]
    is_good = error < tol

    return is_good, error
end

function fix_inverse!(cmat::CeperleyMatrix)
    LU = lufact(cmat.detmat)
    cmat.inverse = inv(LU)
    cmat.det = cmat.old_det = det(LU)
    return nothing
end

# FIXME: should be able to get name of passed argument from argument itself..?
function check_and_fix_inverse!(cmat::CeperleyMatrix, name::ASCIIString)
    cmat_error = check_inverse(cmat)
    if !cmat_error[1]
        printlnf("Recomputing $(name).inverse from scratch!  Error = $(cmat_error[2])")
        fix_inverse!(cmat)
    end
    # FIXME: have this return something so that we can keep track of how often errors occurred
    return nothing
end

function check_and_reset_inverse_and_det!(cmat::CeperleyMatrix, name::ASCIIString,
                                          inv_tol::Float64 = 1e-10, det_tol::Float64 = 1e-10)
    @assert cmat.current_state == :READY_FOR_UPDATE

    LU = lufact(cmat.detmat)
    inv_exact = inv(LU)
    det_exact = det(LU)

    inv_error = sum(abs(cmat.inverse - inv_exact)) / size(inv_exact)[1]
    det_error = abs((cmat.det - det_exact) / det_exact)

    if inv_error > inv_tol
        printlnf("Inverse error for $(name) = $(inv_error)")
    end
    if det_error > det_tol
        printlnf("Determinant error for $(name) = $(det_error)")
    end

    # reset inverse and determinant to exact values
    cmat.inverse = inv_exact
    cmat.det = cmat.old_det = det_exact

    # do we want/need to return this?
    return inv_error, det_error
end

function get_detrat(cmat::CeperleyMatrix)
    @assert cmat.current_state == :COLUMNS_UPDATE_PENDING ||
            cmat.current_state == :ROWS_UPDATE_PENDING ||
            cmat.current_state == :COLUMNS_AND_ROWS_UPDATE_PENDING
    return cmat.detrat
end

# I don't think it's appropriate for another class to get det and detmat during an update (see asserts)
function get_det(cmat::CeperleyMatrix)
    @assert cmat.current_state == :READY_FOR_UPDATE
    return cmat.det
end

function get_detmat(cmat::CeperleyMatrix)
    @assert cmat.current_state == :READY_FOR_UPDATE
    return cmat.detmat
end

function get_det_exact(cmat::CeperleyMatrix)
    @assert cmat.current_state == :READY_FOR_UPDATE
    return det(cmat.detmat)
end

function relative_determinant_error(cmat::CeperleyMatrix)
    @assert cmat.current_state == :READY_FOR_UPDATE
    det_exact = det(cmat.detmat)
    return abs((cmat.det - det_exact) / det_exact)
end

###########################################################
#################### Swap columns/rows ####################
###########################################################

function swap_columns!{T}(A::Matrix{T}, c1::Int, c2::Int)
    m, n = size(A)
    if (1 <= c1 <= n) && (1 <= c2 <= n)
        @inbounds begin
            for i in 1:m
                A[i, c1], A[i, c2] = A[i, c2], A[i, c1]
            end
        end
        return A
    else
        throw(BoundsError())
    end
end

function swap_columns{T}(A::Matrix{T}, c1::Int, c2::Int)
    m, n = size(A)
    Anew = copy(A)
    if (1 <= c1 <= n) && (1 <= c2 <= n)
        @inbounds begin
            for i in 1:m
                Anew[i, c1], Anew[i, c2] = A[i, c2], A[i, c1]
            end
        end
        return Anew
    else
        throw(BoundsError())
    end
end

function swap_rows!{T}(A::Matrix{T}, r1::Int, r2::Int)
    m, n = size(A)
    if (1 <= r1 <= m) && (1 <= r2 <= m)
        @inbounds begin
            for j in 1:n
                A[r1, j], A[r2, j] = A[r2, j], A[r1, j]
            end
        end
        return A
    else
        throw(BoundsError())
    end
end

function swap_rows{T}(A::Matrix{T}, r1::Int, r2::Int)
    m, n = size(A)
    Anew = copy(A)
    if (1 <= r1 <= m) && (1 <= r2 <= m)
        @inbounds begin
            for j in 1:n
                Anew[r1, j], Anew[r2, j] = A[r2, j], A[r1, j]
            end
        end
        return Anew
    else
        throw(BoundsError())
    end
end

function swap_columns!(cmat::CeperleyMatrix, c1::Int, c2::Int)
    @assert c1 <= size(cmat.detmat)[2]
    @assert c2 <= size(cmat.detmat)[2]
    @assert c1 != c2
    @assert cmat.current_state == :READY_FOR_UPDATE

    swap_columns!(cmat.detmat, c1, c2)
    swap_rows!(cmat.inverse, c1, c2)

    cmat.det *= -1
    # no need to do anything to detrat since it's only relevant while an update is pending

    return cmat
end

function swap_rows!(cmat::CeperleyMatrix, r1::Int, r2::Int)
    @assert r1 <= size(cmat.detmat)[1]
    @assert r2 <= size(cmat.detmat)[1]
    @assert r1 != r2
    @assert cmat.current_state == :READY_FOR_UPDATE

    swap_rows!(cmat.detmat, r1, r2)
    swap_columns!(cmat.inverse, r1, r2)

    cmat.det *= -1
    # no need to do anything to detrat since it's only relevant while an update is pending

    return cmat
end

#########################################################################
#################### Start updates (compute detrats) ####################
#########################################################################

# NB: In the detrat_from_* functions, we keep the type parameter in the definition to assert that the
# type of the input CeperleyMatrix object and the cols / rows Matrix objects have the same type T.

# single column
function detrat_from_columns_update!{T}(cmat::CeperleyMatrix{T}, c::Int, col::Matrix{T})
    @assert size(col)[2] == 1
    @assert c <= size(cmat.detmat)[2]
    @assert size(col)[1] == size(cmat.detmat)[1]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    cmat.detrat = (cmat.inverse[c, :] * col)[1]
    cmat.Q_inv = [1. / cmat.detrat].'

    cmat.old_det = cmat.det
    cmat.det *= cmat.detrat

    cmat.cs = [c]
    cmat.pending_cols = col
    cmat.current_state = :COLUMNS_UPDATE_PENDING

    return cmat
end

# single column (no bang)
function detrat_from_columns_update{T}(cmat::CeperleyMatrix{T}, c::Int, col::Matrix{T})
    @assert size(col)[2] == 1
    @assert c <= size(cmat.detmat)[2]
    @assert size(col)[1] == size(cmat.detmat)[1]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    return (cmat.inverse[c, :] * col)[1]
end

# multiple columns
function detrat_from_columns_update!{T}(cmat::CeperleyMatrix{T}, cs::Vector{Int}, cols::Matrix{T})
    @assert length(cs) == size(cols)[2]
    @assert length(cs) <= size(cmat.detmat)[2]
    # these asserts are somewhat costly (the first more so)
    # @assert all(cs .<= size(cmat.detmat)[2])
    # @assert cs == unique(cs)
    @assert size(cols)[1] == size(cmat.detmat)[1]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # generalized matrix determinant lemma
    Q = zeros(T, (length(cs), length(cs))) # should we even allocate here?
    Q = cmat.inverse[cs, :] * cols
    LUQ = lufact(Q) # = lufact(cmat.inverse[cs, :] * cols)
    cmat.detrat = det(LUQ)
    cmat.Q_inv = inv(LUQ)

    cmat.old_det = cmat.det
    cmat.det *= cmat.detrat

    cmat.cs = cs
    cmat.pending_cols = cols
    cmat.current_state = :COLUMNS_UPDATE_PENDING

    return cmat
end

# multiple columns (no bang)
function detrat_from_columns_update{T}(cmat::CeperleyMatrix{T}, cs::Vector{Int}, cols::Matrix{T})
    @assert length(cs) == size(cols)[2]
    @assert length(cs) <= size(cmat.detmat)[2]
    # these asserts are somewhat costly (the first more so)
    # @assert all(cs .<= size(cmat.detmat)[2])
    # @assert cs == unique(cs)
    @assert size(cols)[1] == size(cmat.detmat)[1]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # generalized matrix determinant lemma
    return det(cmat.inverse[cs, :] * cols)
end

# single row
function detrat_from_rows_update!{T}(cmat::CeperleyMatrix{T}, r::Int, row::Matrix{T})
    @assert size(row)[1] == 1
    @assert r <= size(cmat.detmat)[1]
    @assert size(row)[2] == size(cmat.detmat)[2]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    cmat.detrat = (row * cmat.inverse[:, r])[1]
    cmat.Q_inv = [1. / cmat.detrat].'

    cmat.old_det = cmat.det
    cmat.det *= cmat.detrat

    cmat.rs = [r]
    cmat.pending_rows = row
    cmat.current_state = :ROWS_UPDATE_PENDING

    return cmat
end

# single row (no bang)
function detrat_from_rows_update{T}(cmat::CeperleyMatrix{T}, r::Int, row::Matrix{T})
    @assert size(row)[1] == 1
    @assert r <= size(cmat.detmat)[1]
    @assert size(row)[2] == size(cmat.detmat)[2]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    return (row * cmat.inverse[:, r])[1]
end

# multiple rows
function detrat_from_rows_update!{T}(cmat::CeperleyMatrix{T}, rs::Vector{Int}, rows::Matrix{T})
    @assert length(rs) == size(rows)[1]
    @assert length(rs) <= size(cmat.detmat)[1]
    # these asserts are somewhat costly (the first more so)
    # @assert all(rs .<= size(cmat.detmat)[1])
    # @assert rs == unique(rs)
    @assert size(rows)[2] == size(cmat.detmat)[2]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # generalized matrix determinant lemma
    Q = zeros(T, (length(rs), length(rs))) # should we even allocate here?
    Q = rows * cmat.inverse[:, rs]
    LUQ = lufact(Q) # = lufact(rows * cmat.inverse[:, rs])
    cmat.detrat = det(LUQ)
    cmat.Q_inv = inv(LUQ)

    cmat.old_det = cmat.det
    cmat.det *= cmat.detrat

    cmat.rs = rs
    cmat.pending_rows = rows
    cmat.current_state = :ROWS_UPDATE_PENDING

    return cmat
end

# multiple rows (no bang)
function detrat_from_rows_update{T}(cmat::CeperleyMatrix{T}, rs::Vector{Int}, rows::Matrix{T})
    @assert length(rs) == size(rows)[1]
    @assert length(rs) <= size(cmat.detmat)[1]
    # these asserts are somewhat costly (the first more so)
    # @assert all(rs .<= size(cmat.detmat)[1])
    # @assert rs == unique(rs)
    @assert size(rows)[2] == size(cmat.detmat)[2]
    @assert cmat.current_state == :READY_FOR_UPDATE

    # generalized matrix determinant lemma
    return det(rows * cmat.inverse[:, rs])
end

# multiple (in general) rows and columns
function detrat_from_columns_and_rows_update!{T}(cmat::CeperleyMatrix{T},
                                                 cs::Vector{Int}, cols::Matrix{T},
                                                 rs::Vector{Int}, rows::Matrix{T})
    @assert length(cs) == size(cols)[2]
    @assert length(rs) == size(rows)[1]
    @assert length(cs) <= size(cmat.detmat)[2]
    @assert length(rs) <= size(cmat.detmat)[1]
    @assert all(cs .<= size(cmat.detmat)[2])
    @assert all(rs .<= size(cmat.detmat)[1])
    @assert cs == unique(cs)
    @assert rs == unique(rs)
    @assert size(cols)[1] == size(cmat.detmat)[1]
    @assert size(rows)[2] == size(cmat.detmat)[2]
    @assert cmat.current_state == :READY_FOR_UPDATE

    if cols[rs, :] != rows[:, cs]
        warn("cols[rs, :] != rows[:, cs] in detrat_from_columns_and_rows_update.
Proceeding anyway, putting the columns on top..") # see Q[nc+1:end, 1:nc] assignment below
    end

    cols_tilde = copy(cols)
    rows_tilde = copy(rows)
    cols_tilde[rs, :] = 0.5(cols[rs, :] + cmat.detmat[rs, cs])
    rows_tilde[:, cs] = 0.5(rows[:, cs] + cmat.detmat[rs, cs])

    # generalized matrix determinant lemma
    nc = length(cs); nr = length(rs)
    Q = zeros(T, (nc + nr, nc + nr))
    Q[1:nc, 1:nc] = cmat.inverse[cs, :] * cols_tilde
    Q[1:nc, nc+1:end] = cmat.inverse[cs, rs]
    Q[nc+1:end, nc+1:end] = rows_tilde * cmat.inverse[:, rs]
    Q[nc+1:end, 1:nc] = rows_tilde * cmat.inverse * cols_tilde - 2cols_tilde[rs, :] + cmat.detmat[rs, cs]
    # Q[nc+1:end, 1:nc] = rows_tilde * cmat.inverse * cols_tilde - 2rows_tilde[:, cs] + cmat.detmat[rs, cs]

    LUQ = lufact(Q)
    cmat.detrat = det(LUQ)
    cmat.Q_inv = inv(LUQ)

    cmat.old_det = cmat.det
    cmat.det *= cmat.detrat

    cmat.cs = cs
    cmat.rs = rs
    cmat.pending_cols = cols
    cmat.pending_rows = rows
    cmat.cols_tilde = cols_tilde
    cmat.rows_tilde = rows_tilde
    cmat.current_state = :COLUMNS_AND_ROWS_UPDATE_PENDING

    return cmat
end

# multiple (in general) rows and columns (no bang)
function detrat_from_columns_and_rows_update{T}(cmat::CeperleyMatrix{T},
                                                cs::Vector{Int}, cols::Matrix{T},
                                                rs::Vector{Int}, rows::Matrix{T})
    @assert length(cs) == size(cols)[2]
    @assert length(rs) == size(rows)[1]
    @assert length(cs) <= size(cmat.detmat)[2]
    @assert length(rs) <= size(cmat.detmat)[1]
    @assert all(cs .<= size(cmat.detmat)[2])
    @assert all(rs .<= size(cmat.detmat)[1])
    @assert cs == unique(cs)
    @assert rs == unique(rs)
    @assert size(cols)[1] == size(cmat.detmat)[1]
    @assert size(rows)[2] == size(cmat.detmat)[2]
    @assert cmat.current_state == :READY_FOR_UPDATE

    if cols[rs, :] != rows[:, cs]
        warn("cols[rs, :] != rows[:, cs] in detrat_from_columns_and_rows_update.
Proceeding anyway, putting the columns on top..") # see Q[nc+1:end, 1:nc] assignment below
    end

    cols_tilde = copy(cols)
    rows_tilde = copy(rows)
    cols_tilde[rs, :] = 0.5(cols[rs, :] + cmat.detmat[rs, cs])
    rows_tilde[:, cs] = 0.5(rows[:, cs] + cmat.detmat[rs, cs])

    # generalized matrix determinant lemma
    nc = length(cs); nr = length(rs)
    Q = zeros(T, (nc + nr, nc + nr))
    Q[1:nc, 1:nc] = cmat.inverse[cs, :] * cols_tilde
    Q[1:nc, nc+1:end] = cmat.inverse[cs, rs]
    Q[nc+1:end, nc+1:end] = rows_tilde * cmat.inverse[:, rs]
    Q[nc+1:end, 1:nc] = rows_tilde * cmat.inverse * cols_tilde - 2cols_tilde[rs, :] + cmat.detmat[rs, cs]
    # Q[nc+1:end, 1:nc] = rows_tilde * cmat.inverse * cols_tilde - 2rows_tilde[:, cs] + cmat.detmat[rs, cs]

    return det(Q)
end

# single column, multiple rows (in general)
detrat_from_columns_and_rows_update!{T}(cmat::CeperleyMatrix{T},
                                        c::Int, cols::Matrix{T},
                                        rs::Vector{Int}, rows::Matrix{T}) = (
    detrat_from_columns_and_rows_update!(cmat, [c], cols, rs, rows) )
detrat_from_columns_and_rows_update{T}(cmat::CeperleyMatrix{T},
                                       c::Int, cols::Matrix{T},
                                       rs::Vector{Int}, rows::Matrix{T}) = (
    detrat_from_columns_and_rows_update(cmat, [c], cols, rs, rows) )

# multiple columns (in general), single row
detrat_from_columns_and_rows_update!{T}(cmat::CeperleyMatrix{T},
                                        cs::Vector{Int}, cols::Matrix{T},
                                        r::Int, rows::Matrix{T}) = (
    detrat_from_columns_and_rows_update!(cmat, cs, cols, [r], rows) )
detrat_from_columns_and_rows_update{T}(cmat::CeperleyMatrix{T},
                                       cs::Vector{Int}, cols::Matrix{T},
                                       r::Int, rows::Matrix{T}) = (
    detrat_from_columns_and_rows_update(cmat, cs, cols, [r], rows) )

# single column, single row
detrat_from_columns_and_rows_update!{T}(cmat::CeperleyMatrix{T},
                                        c::Int, cols::Matrix{T},
                                        r::Int, rows::Matrix{T}) = (
    detrat_from_columns_and_rows_update!(cmat, [c], cols, [r], rows) )
detrat_from_columns_and_rows_update{T}(cmat::CeperleyMatrix{T},
                                        c::Int, cols::Matrix{T},
                                        r::Int, rows::Matrix{T}) = (
    detrat_from_columns_and_rows_update(cmat, [c], cols, [r], rows) )

###########################################################################
#################### Finish updates (compute inverses) ####################
###########################################################################

function finish_columns_update!(cmat::CeperleyMatrix)
    @assert cmat.current_state == :COLUMNS_UPDATE_PENDING

    cmat.detmat[:, cmat.cs] = cmat.pending_cols

    # Sherman-Morrison-Woodbury
    old_inverse_rows = cmat.inverse[cmat.cs, :]
    cmat.inverse -= (cmat.inverse * cmat.pending_cols) * cmat.Q_inv * cmat.inverse[cmat.cs, :]
    cmat.inverse[cmat.cs, :] = cmat.Q_inv * old_inverse_rows

    cmat.old_det = cmat.det

    cmat.current_state = :READY_FOR_UPDATE

    return cmat
end

function finish_rows_update!(cmat::CeperleyMatrix)
    @assert cmat.current_state == :ROWS_UPDATE_PENDING

    cmat.detmat[cmat.rs, :] = cmat.pending_rows

    # Sherman-Morrison-Woodbury
    old_inverse_cols = cmat.inverse[:, cmat.rs]
    cmat.inverse -= cmat.inverse[:, cmat.rs] * cmat.Q_inv * (cmat.pending_rows * cmat.inverse)
    cmat.inverse[:, cmat.rs] = old_inverse_cols * cmat.Q_inv

    cmat.old_det = cmat.det

    cmat.current_state = :READY_FOR_UPDATE

    return cmat
end

function finish_columns_and_rows_update!(cmat::CeperleyMatrix)
    @assert cmat.current_state == :COLUMNS_AND_ROWS_UPDATE_PENDING

    # put the columns down after the rows (see detrat_from_columns_and_rows_update)
    cmat.detmat[cmat.rs, :] = cmat.pending_rows
    cmat.detmat[:, cmat.cs] = cmat.pending_cols

    # Sherman-Morrison-Woodbury
    nc = length(cmat.cs);
    old_inverse = copy(cmat.inverse)
    cmat.inverse -= ( (cmat.inverse * cmat.cols_tilde) * cmat.Q_inv[1:nc, 1:nc] * cmat.inverse[cmat.cs, :]
                     + cmat.inverse[:, cmat.rs] * cmat.Q_inv[nc+1:end, nc+1:end] * (cmat.rows_tilde * cmat.inverse)
                     + cmat.inverse[:, cmat.rs] * cmat.Q_inv[nc+1:end, 1:nc] * cmat.inverse[cmat.cs, :]
                     + (cmat.inverse * cmat.cols_tilde) * cmat.Q_inv[1:nc, nc+1:end] * (cmat.rows_tilde * cmat.inverse) )
    cmat.inverse[cmat.cs, :] = ( cmat.Q_inv[1:nc, 1:nc] * old_inverse[cmat.cs, :]
                               + cmat.Q_inv[1:nc, nc+1:end] * (cmat.rows_tilde * old_inverse) )
    cmat.inverse[:, cmat.rs] = ( old_inverse[:, cmat.rs] * cmat.Q_inv[nc+1:end, nc+1:end]
                               + (old_inverse * cmat.cols_tilde) * cmat.Q_inv[1:nc, nc+1:end] )
    cmat.inverse[cmat.cs, cmat.rs] = -cmat.Q_inv[1:nc, nc+1:end]

    cmat.old_det = cmat.det

    cmat.current_state = :READY_FOR_UPDATE

    return cmat
end

########################################################
#################### Cancel updates ####################
########################################################

function cancel_columns_update!(cmat::CeperleyMatrix)
    @assert cmat.current_state == :COLUMNS_UPDATE_PENDING
    cmat.det = cmat.old_det
    cmat.current_state = :READY_FOR_UPDATE
    return cmat
end

function cancel_rows_update!(cmat::CeperleyMatrix)
    @assert cmat.current_state == :ROWS_UPDATE_PENDING
    cmat.det = cmat.old_det
    cmat.current_state = :READY_FOR_UPDATE
    return cmat
end

function cancel_columns_and_rows_update!(cmat::CeperleyMatrix)
    @assert cmat.current_state == :COLUMNS_AND_ROWS_UPDATE_PENDING
    cmat.det = cmat.old_det
    cmat.current_state = :READY_FOR_UPDATE
    return cmat
end

########################################################################
#################### Householder --> exact Pfaffian ####################
########################################################################

is_antisymmetric(A, tol::Float64 = 1e-15) = maximum(abs(A + A.')) < tol

function Householder(x::Vector{Float64})
    absx = norm(x)
    u = copy(x)
    u[1] = x[1] + (x[1] > 0 ? absx : -absx)
    beta = 2 / norm(u)^2
    return u, beta
end

function Householder(x::Vector{Complex128})
    absx = norm(x)
    u = copy(x)
    u1p = x[1] + exp(im*angle(x[1])) * absx
    u1m = x[1] - exp(im*angle(x[1])) * absx
    u[1] = abs(u1p) > abs(u1m) ? u1p : u1m
    beta = 2 / norm(u)^2
    return u, beta
end

function Householder{T}(x::Matrix{T})
    @assert size(x)[2] == 1 || size(x)[1] == 1
    return Householder(vec(x))
end

# FIXME: this implementation is very slow/naive -- use Pfaffian.jl instead
function Pf{T}(A::Matrix{T})
    @assert size(A)[1] == size(A)[2]
    @assert is_antisymmetric(A)
    N = size(A)[1]

    if mod(N, 2) == 1
        return 0.
    end

    for i in 1:2:N-2 # only need to do every other row/col to get actual Pfaffian (see sign in return statement)
        u, beta = Householder(A[i, i+1:end])
        u = [zeros(T, i); u]
        p = beta * A * conj(u)
        K = 0.5beta * (u' * p)[1]
        A = A - p * u.' + u * p.' + 2K * u * u.'
        # @assert is_antisymmetric(A) # this is too harsh
    end

    # = sign from Householder determinants * product of upper diagonal, with only every other element included
    return (-1.)^((N-2)/2) * prod(diag(A, 1)[1:2:end])
end

################################################################
#################### Pfaffian wavefunctions ####################
################################################################

# Pfaffian <--> Ceperley
# pfaff <--> det
# pf <--> c

type PfaffianMatrix{T} #<: VMCMatrix{T}
    current_state::Symbol

    Amat::Matrix{T}
    inverse::Matrix{T}

    # pfaff::T
    pfaffrat::T
    # Q_inv::Matrix{T}

    # pending index (only have single row/column complete updates currently available)
    ks::Vector{Int}

    # pending row/col (only have single row/column complete updates currently available)
    pending_rows::Matrix{T}

    PfaffianMatrix() = new(:UNINITIALIZED)

    function PfaffianMatrix(A::Matrix{T})
        @assert size(A)[1] == size(A)[2]
        @assert is_antisymmetric(A)
        return new(:READY_FOR_UPDATE, A, inv(A))
    end
end

PfaffianMatrix(A::Matrix) = PfaffianMatrix{eltype(A)}(A)
PfaffianMatrix() = PfaffianMatrix{Float64}() # default to Float64

# FIXME: Can we have something like a VMCMatrix abstract base class,
# so we don't need to redefine the following functions from their Ceperley counterparts?

# check inverse for numerical error due to SMW
function check_inverse{T}(pfmat::PfaffianMatrix{T}, tol::Float64 = 1e-10)
    @assert pfmat.current_state == :READY_FOR_UPDATE

    error = sum(abs(pfmat.Amat * pfmat.inverse - eye(T, size(pfmat.Amat)[1]))) / size(pfmat.Amat)[1]
    is_good = error < tol

    return is_good, error
end

function fix_inverse!(pfmat::PfaffianMatrix)
    pfmat.inverse = inv(pfmat.Amat)
    return pfmat
end

function get_pfaffrat(pfmat::PfaffianMatrix)
    @assert pfmat.current_state == :UPDATE_PENDING
    return pfmat.pfaffrat
end

###########################################################################
#################### Start updates (compute pfaffrats) ####################
###########################################################################

# NB: In the pfaffrat_from_* functions, we keep the type parameter in the definition to assert that the
# type of the input PfaffianMatrix object and the cols / rows Matrix objects have the same type T.
# NB: Below, only ever actually use row updates, since antisymmetry makes row and column updates equivalent.

# single row
function pfaffrat_from_rows_update!{T}(pfmat::PfaffianMatrix{T}, k::Int, row::Matrix{T})
    @assert size(row)[1] == 1
    @assert row[1, k] == 0 # so the new A can be antisymmetric
    @assert k <= size(pfmat.Amat)[1]
    @assert size(row)[2] == size(pfmat.Amat)[2]
    @assert pfmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    pfmat.pfaffrat = (row * pfmat.inverse[:, k])[1]

    pfmat.ks = [k]
    pfmat.pending_rows = row
    pfmat.current_state = :UPDATE_PENDING

    return pfmat
end

# single row (no bang)
function pfaffrat_from_rows_update{T}(pfmat::PfaffianMatrix{T}, k::Int, row::Matrix{T})
    @assert size(row)[1] == 1
    @assert row[1, k] == 0 # so the new A can be antisymmetric
    @assert k <= size(pfmat.Amat)[1]
    @assert size(row)[2] == size(pfmat.Amat)[2]
    @assert pfmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    return (row * pfmat.inverse[:, k])[1]
end

# 2 rows (no bang)
function pfaffrat_from_rows_update{T}(pfmat::PfaffianMatrix{T}, ks::Vector{Int}, rows::Matrix{T})
    @assert length(ks) == 2
    @assert size(rows)[1] == 2
    @assert rows[1, ks[1]] == rows[2, ks[2]] == 0 # so the new A can be antisymmetric
    @assert maximum(ks) <= size(pfmat.Amat)[1]
    @assert size(rows)[2] == size(pfmat.Amat)[2]
    @assert pfmat.current_state == :READY_FOR_UPDATE

    # matrix determinant lemma
    return ( rows[2, ks[1]] * pfmat.inverse[ks[1], ks[2]]
          + (rows[2, :] * pfmat.inverse[:, ks[2]])[1] * (rows[1, :] * pfmat.inverse[:, ks[1]])[1]
          - (rows[2, :] * pfmat.inverse[:, ks[1]])[1] * (rows[1, :] * pfmat.inverse[:, ks[2]])[1]
          + pfmat.inverse[ks[1], ks[2]] * (rows[2, :] * pfmat.inverse * rows[1, :].')[1] )
          # note sign/transpose on last factor of last term (see notes)
end

# single column
function pfaffrat_from_columns_update!{T}(pfmat::PfaffianMatrix{T}, k::Int, col::Matrix{T})
    return pfaffrat_from_rows_update!(pfmat, k, -col.')
end

# single column (no bang)
function pfaffrat_from_columns_update{T}(pfmat::PfaffianMatrix{T}, k::Int, col::Matrix{T})
    return pfaffrat_from_rows_update(pfmat, k, -col.')
end

# 2 columns (no bang)
function pfaffrat_from_columns_update{T}(pfmat::PfaffianMatrix{T}, ks::Vector{Int}, cols::Matrix{T})
    return pfaffrat_from_rows_update(pfmat, ks, -cols.')
end

###########################################################################
#################### Finish updates (compute inverses) ####################
###########################################################################

# only single row/col currently available
function finish_update!(pfmat::PfaffianMatrix)
    @assert pfmat.current_state == :UPDATE_PENDING

    pfmat.Amat[pfmat.ks, :] = pfmat.pending_rows
    pfmat.Amat[:, pfmat.ks] = -pfmat.pending_rows.'

    # Sherman-Morrison-Woodbury
    old_inverse_rows = pfmat.inverse[pfmat.ks, :]
    pfmat.inverse -= ( (pfmat.inverse * -pfmat.pending_rows.') * pfmat.inverse[pfmat.ks, :]
                      + pfmat.inverse[:, ks] * (pfmat.pending_rows * pfmat.inverse) ) / pfmat.pfaffrat
    pfmat.inverse[ks, :] = old_inverse_rows / pfmat.pfaffrat
    pfmat.inverse[:, ks] = -pfmat.inverse[ks, :].'
    pfmat.inverse[ks, ks] = 0

    pfmat.current_state = :READY_FOR_UPDATE

    return pfmat
end

########################################################
#################### Cancel updates ####################
########################################################

function cancel_update!(pfmat::PfaffianMatrix)
    @assert pfmat.current_state == :UPDATE_PENDING
    pfmat.current_state = :READY_FOR_UPDATE
    return pfmat
end


text:
# Packages used in this chapter
# To check for installed packages in V1.0 we now need to use the PKG API

using Pkg
ipks = Pkg.installed()
pkgs = ["BenchmarkTools","MacroTools"]
for p in pkgs
    !haskey(ipks,p) && Pkg.add(p)
end

# We can define these here or as required
using BenchmarkTools, MacroTools, PyPlot

# We need these as some functions have moved from Base to Stdlib
#
using Printf, Statistics

recip(x::Number) =
(x == zero(typeof(x))) ? error("Invalid reciprocal") : one(typeof(x)) / x

recip(2)

recip(recip(2))  # recip does not ALWAYS return the same type

recip(11//17)

recip(11 + 17im)

aa = rand(3)
recip(aa)

## recip(a::Array) = [a[i] = recip(a[i]) for i = 1:length(a)]
## OR
recip(a::Array) = map(recip,a)
recip(aa)

map(sin,recip(aa)) # and can map functions to the array

bb = [2.1 3.2 4.3; 9.8 8.7 7.6]

recip(bb) # our definition works but returns a 1-D array

cc = recip(aa)'.*recip(bb)  # because they are 1-D arrays we can do this.

incr(x) = x + 1

@code_native incr(2)

@code_native incr(2.7)

@code_native incr(2//7)

@code_native incr(2.0 + 7.0im)

# Look at the intermediat stages

dump(:(incr(2.7)))

@code_lowered incr(2.7)

@code_typed incr(2.7)

@code_llvm(incr(2.7))

ex1 = :((x^2 + y^2 - 2*x*y)^0.5)

ex2 = quote
    (x^2 + y^2 - 2*x*y)^0.5
end

# If we instantiate the variables x,y,z then we can evaluate the expression
# This is NOT a function call, so values need to be known before hand
#
x = 1.1; y = 2.5;
eval(ex1)

# Note that the @eval macro does NOT behave as the function call
@eval ex1

# Unless the expression is prefixed with '$'
@eval $ex1

# i.e.  @eval $(ex) === eval(ex)

# Although ex1 and ex2 are equivalent the AST of the expression
# is slightly different
# Here it is for ex1, we leave ex2 to the reader
#
dump(ex1)

# Dump is very verbose, so using show_expr creates an S-expression version
#
Meta.show_sexpr(ex1)

# Traversing a tree
#
function traverse!(ex, symbols) end

function traverse!(ex::Symbol, symbols)
    push!(symbols, ex)
end

function traverse!(ex::Expr, symbols)
    if ex.head == :call  # function call
        for arg in ex.args[2:end]
            traverse!(arg, symbols)  # recursive
        end
    else
        for arg in ex.args
            traverse!(arg, symbols)  # recursive
        end
    end
end

function traverse(ex::Expr)
    symbols = Symbol[]
    traverse!(ex, symbols)
    return unique(symbols)  # Don't output duplicate
end

traverse(ex1)

macro pout(n)
    if typeof(n) == Expr
       println(n.args)
    end
    return n
end

@pout x

@pout (x^2 + y^2 - 2*x*y)^0.5

macro dotimes(n, body)
    quote
        for i = 1:$(esc(n))
            $(esc(body))
        end
    end
end

@dotimes 3 print("Hi")

i = 0; @dotimes 3 [global i += 1; println(i*i)]

import Base.macroexpand
macroexpand(ex::Expr) = macroexpand(Main,ex)

macroexpand(:(@assert n > 0))

# Expand our @dotimes function, whic is somewhat simpler
#
macroexpand(Main,:(@dotimes 3 [global i += 1; println(i*i)]))

# Not all macro expansions produce short boiler plate
#
using Printf, Statistics
aa = [rand() for i = 1:100000]
@printf "The average value is %f.3 over %d trials"  mean(aa) length(aa)

macroexpand(:(@printf "The average value is %f.3 over %d trials"  mean(aa) length(aa)))

# A more useful version of dotimes macro is until
# Creates a loop and breaks out when a condition fails
#
macro until(condition, block)
    quote
        while true
            $(esc(block))
            if $(esc(condition))
                break
            end
        end
    end
end

i = 0; @until (i >= 3) [global i += 1; println(i*i)]

# The until macro can be used to implement a simple if-then-else
#
macro trigger(cond, body1, body2)
    :(if !$cond
        $(esc(body1))
    else
        $(esc(body2))
    end)
end

fac(n::Integer) = (n == 1) ? 1 : n*fac(n-1)

n = 10
@trigger (n < 1) fac(n) ArgumentError("$n not positive")

n = -1; @trigger (n < 1) fac(n) ArgumentError("$n not positive")

# Evaluate a polynomial with a function
# This is NOT a macro version

function poly_native(x, a...)
  p=zero(x)
  for i = 1:length(a)
    p = p + a[i] *  x^(i-1)
  end
  return p
end

f_native(x) = poly_native(x,1,2,3,4,5)
f_native(2.1)

# Neither IS this!

function poly_horner(x, a...)
  b=zero(x)
  for i = length(a):-1:1
    b = a[i] + b * x
  end
  return b
end

# f -> (((5*x + 4)*x + 3)*x + 2)*x + 1
#
f_horner(x) = poly_horner(x,1,2,3,4,5)
f_horner(2.1)

# for the macro version we need a help function
# Define mad(x,a,b)
# [Julia has this function too :- muladd(x,a,b)]

mad(x,a,b) = a*x + b
mad(2.1,5,4)

# And NOW use horner's method in a macro
# [ p... is a variable list of arguments passed as an array]

macro horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(mad(t, $ex, $(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

@horner 2.1 1 2 3 4 5

# The saving of eliminating the loop for a 5th order polynomial is not great
# However for larger arrays and complex calculations this can be substantial
#
macroexpand(Main,:(@horner 2.1 1 2 3 4 5))

# Mike Innes' Macrotools has a useful set of macros and utility functions
# We will discuss the package in more detail later

using MacroTools:postwalk

# The postwalk function splits an expression into symbols and then
# reconstructs it, so we can apply different operations to each symbol

ex = :(1 + (2 + 3) + 4)
p = postwalk(ex) do x
    x isa Integer ? fac(x) : x
end

# Evaluate the expression
eval(p)

map(x -> @show(x), [1,2,3,4]);

postwalk(ex) do x
         @show x
end

@capture(ex, a_ + b_ + c_)

b

a*eval(b) + c  #=> 1*5 + 4

reduce(+, 1:10)

plus(a, b) = :($a + $b)

p = reduce(plus, 1:10)

eval(p)

# Do something useful
# This is the series expansion for the SINE function
# We need the factorial function and will use the version from the STDLIB
#
using SpecialFunctions
k = 2
pp = [:($((-1)^k) * x^$(1+2k) / $(factorial(1+2k))) for k = 0:5]

# We can reduce this to a single expression
reduce(plus,pp)

# ... and evaluate it for a specific value of x
x = 2.1; eval(reduce(plus,pp))

# A simple generated function to execute a multiply + add
# The resukt returned is a symbol
#
@generated function mad(a,b,c)
    Core.println("Calculating: a*b + c")
    return :(a * b + c)
end

# Call the function
mad(2.3,1.7,1.1)

# And again with the same arguments
# This time it is not reevaluated
#
mad(2.3,1.7,1.1)

# With different arguments the function is reevaluated
mad(2.3,1.7,1)

# To illustrate the use consider the following function ...
# ... which multiplies the size of dimensions of a n-D array
#
# First a 'normal' function definition

function pdims(x::Array{T,N}) where {T,N}
  s = 1
  for i = 1:N
    s = s * size(x, i)
  end
  return s
end

# And then a 'generated function' version
@generated function gpdims(x::Array{T,N}) where {T,N}
  ex = :(1)
  for i = 1:N
     ex = :(size(x, $i) * $ex)
  end
  return ex
end

# We need an array to test the function
aa = reshape([rand() for i = 1:1000],10,5,5,4); size(aa)

# Both functions provide the same result
#
pdims(aa) == gpdims(aa)

# But if we look at the lowered code
@code_lowered pdims(aa)

@code_lowered gpdims(aa)

# Which results in very different generated native code
@code_native pdims(aa)

@code_native gpdims(aa)

# Previous versions of Julia had an example section, absent in version v1.0
# An interesting one was modular integer, which required a little re-editing
#
module ModInts
export ModInt

import Base: +, -, *, /, inv

struct ModInt{n} <: Integer
  k::Int
  ModInt{n}(k) where n = new(mod(k,n))
end

Base.show(io::IO, k::ModInt{n}) where n =
    print(io, get(io, :compact, false) ? k.k : "$(k.k) mod $n")

(a::ModInt{n} + b::ModInt{n}) where n = ModInt{n}(a.k+b.k)
(a::ModInt{n} - b::ModInt{n}) where n = ModInt{n}(a.k-b.k)
(a::ModInt{n} * b::ModInt{n}) where n = ModInt{n}(a.k*b.k)
-(a::ModInt{n}) where n = ModInt{n}(-a.k)

inv(a::ModInt{n}) where n = ModInt{n}(invmod(a.k, n))
(a::ModInt{n} / b::ModInt{n}) where n = a*inv(b)

Base.promote_rule(::Type{ModInt{n}}, ::Type{Int}) where n = ModInt{n}
Base.convert(::Type{ModInt{n}}, i::Int) where n = ModInt{n}(i)

end

# Test the module
#
using Main.ModInts
m1 = ModInt{11}(2)
m2 = ModInt{11}(7)
m3 = 3*m1 + m2  # => 13 mod 11

# Because of multiple dispatch we can do the following
#
mm = reshape([ModInt{11}(rand(0:10)) for i = 1:100],10,10)
ma = [ModInt{11}(rand(0:10)) for i = 1:10]

mm.*ma'

## An Ordered Pair Module
#
module OrdPairs

import Base: +,-,*,/,==,!=,>,<,>=,<=
import Base: abs,conj,inv,zero,one,show
import LinearAlgebra: transpose,adjoint,norm

export OrdPair

struct OrdPair{T<:Number}
    a::T
    b::T
end

OrdPair(x::Number) = OrdPair(x, zero(T))

value(u::OrdPair)   = u.a
epsilon(u::OrdPair) = u.b

zero(::Type{OrdPairs.OrdPair}) = OrdPair(zero(T),zero(T))
one(::Type{OrdPairs.OrdPair})  = OrdPair(one(T),zero(T))

abs(u::OrdPair)  = abs(value(u))
norm(u::OrdPair) = norm(value(u))

+(u::OrdPair, v::OrdPair) = OrdPair(value(u) + value(v), epsilon(u) + epsilon(v))
-(u::OrdPair, v::OrdPair) = OrdPair(value(u) - value(v), epsilon(u) - epsilon(v))
*(u::OrdPair, v::OrdPair) = OrdPair(value(u)*value(v), epsilon(u)*value(v) + value(u)*epsilon(v))
/(u::OrdPair, v::OrdPair) = OrdPair(value(u)/value(v),(epsilon(u)*value(v) - value(u)*epsilon(v))/(value(v)*value(v)))

==(u::OrdPair, v::OrdPair) = norm(u) == norm(v)
!=(u::OrdPair, v::OrdPair) = norm(u) != norm(v)
>(u::OrdPair, v::OrdPair)  = norm(u) > norm(v)
>=(u::OrdPair, v::OrdPair) = norm(u) >= norm(v)
<(u::OrdPair, v::OrdPair)  = norm(u) < norm(v)
<=(u::OrdPair, v::OrdPair) = norm(u) <= norm(v)

+(x::Number, u::OrdPair) = OrdPair(value(u) + x, epsilon(u))
+(u::OrdPair, x::Number) = x + u

-(x::Number, u::OrdPair) = OrdPair(x - value(u), epsilon(u))
-(u::OrdPair, x::Number) = OrdPair(value(u) - x, epsilon(u))

*(x::Number, u::OrdPair) = OrdPair(x*value(u), x*epsilon(u))
*(u::OrdPair, x::Number) = x*u

/(u::OrdPair, x::Number) = (1.0/x)*u

conj(u::OrdPair)  = OrdPair(value(u),-epsilon(u))
inv(u::OrdPair)   = one(OrdPair)/u

transpose(u::OrdPair) = u
transpose(uu::Array{OrdPairs.OrdPair,2}) = [uu[j,i] for i=1:size(uu)[1],j=1:size(uu)[2]]
adjoint(u::OrdPair) = u

convert(::Type{OrdPair}, x::Number) = OrdPair(x,zero(x))
promote_rule(::Type{OrdPair}, ::Type{<:Number}) = OrdPair

## show(io::IO,u::OrdPair) = print(io,value(u)," + (",epsilon(u),")ϵ")

function show(io::IO,u::OrdPair)
 op::String = (epsilon(u) < 0.0) ? " - " : " + ";
 print(io,value(u),op,abs(epsilon(u)),"ϵ")
end

end

# Exercise the module
# Notice how the show() routine provides pretty-print output
using Main.OrdPairs
p1 = OrdPair(2.3,-1.7)
p2 = OrdPair(4.4,0.9)
p1 * p2

p2/p1

# Again we can operate on array of Ordered Pairs
using Statistics
pp = [OrdPair(rand(),rand()) for i in 1:100]
mean(pp)

# And promote a rational in a mixed OP
p3 = OrdPair(2.3, 11/7)

# The module is not a full implementation
p4 = OrdPair(2.3, 11.0 + 7.2im)

# Also a number of other arithmetic functions need to be imported from Base
std(p1)


text:
# This file was generated, do not modify it. # hide
@load DecisionTreeClassifier
tree_model = DecisionTreeClassifier()


text:
function energy(v, S)
    1/2*(dot(v,inv(S)*v) + log(det(S)))
end

#to do: add prior
function energy(estimates)
    foldxl(+, estimates |> Map(x -> x[end-1:end]) |> Map(x -> energy(x...)))
end


text:
module GenerativePreTrain

using Flux
using Requires
using Requires: @init
using BSON

using ..Transformers: Abstract3DTensor
using ..Basic
using ..Pretrain: isbson, iszip, isnpbson, zipname, zipfile, findfile
export Gpt, load_gpt_pretrain, lmloss

include("./gpt.jl")
include("./tokenizer.jl")
include("./npy2bson.jl")
include("./load_pretrain.jl")

end


text:
using SymPy

@syms px1 py1 lambda1 vx1 vy1 px2 py2 lambda2 vx2 vy2

lx1 = px1 + lambda1*vx1
ly1 = py1 + lambda1*vy1

lx2 = px2 + lambda2*vx2
ly2 = py2 + lambda2*vy2

lx12 = solve(lx1-lx2, lambda2)[1]
l2i = subs(ly2, lambda2, lx12)

sol = solve(ly1-l2i, lambda1)[1]

xi = simplify(subs(lx1, lambda1, sol))
yi = simplify(subs(ly1, lambda1, sol))