text:
# default behavior, dummy behavior, and convenient functions for all widgets

widgetStaggerPosx = 0
widgetStaggerPosy = 0

# only use these inside widget constructor, when their states are
# not yet fully formed.
function link_parent_child( p::TwObj{TwScreenData}, c::TwObj, height::Real, width::Real, posy::Any, posx::Any )
    registerTwObj( p, c )
    alignxy!( c, height,width,posx,posy, parent= p )
    configure_newwinpanel!( c )
    log( "Screen-"*string(objtype(c))*": x=" * string(c.xpos) * " y=" * string(c.ypos) )
end

function link_parent_child( p::TwObj{TwListData}, c::TwObj, height::Real, width::Real, posy::Any, posx::Any )
    @lintpragma( "Ignore unused height")
    @lintpragma( "Ignore unused width")
    @lintpragma( "Ignore unused posy")
    @lintpragma( "Ignore unused posx")
    update_list_canvas( p )
    begx = 0
    begy = 0
    if p.data.horizontal
        for sw in p.data.widgets
            begx += sw.width
        end
    else
        for sw in p.data.widgets
            begy += sw.height
        end
    end

    @assert c.screen.value == nothing
    @assert c.window == nothing

    # by the time a list is being added, its contents must be fully populated
    if objtype( c ) == :List
        @assert c.data.pad == nothing
        update_list_canvas( c )
        height = c.data.canvasheight
        width = c.data.canvaswidth
    end
    push!( p.data.widgets, c )
    update_list_canvas(p)
    alignxy!( c, height,width,begx,begy, parent= p )
    c.hasFocus = false
    c.window = TwWindow( WeakRef( p ), c.ypos, c.xpos, c.height, c.width )
    log( "List-"*string(objtype(c))*": x=" * string(c.xpos) * " y=" * string(c.ypos) )
    log( " orig begxy: x="*string(begx) * " y=" * string(begy) )
    log( " new geom: h=" *string(p.height)*" w=" * string(p.width) )
end

function link_parent_child( p::TwObj, c::TwObj, height::Real, width::Real, posy::Any, posx::Any )
    alignxy!( c, height,width,posx,posy, parent= p )
    c.window = TwWindow( WeakRef( p ), c.ypos, c.xpos, c.height, c.width )
    log( string(objtype(p))*"-"*string(objtype(c))*": x=" * string(c.xpos) * " y=" * string(c.ypos) )
end

function configure_newwinpanel!( obj::TwObj )
    obj.window = newwin( obj.height,obj.width,obj.ypos,obj.xpos )
    obj.panel = new_panel( obj.window )
    cbreak()
    noecho()
    keypad( obj.window, true )
    nodelay( obj.window, true )
    wtimeout( obj.window, 100 )
    curs_set( 0 )
end

function alignxy!( o::TwObj, h::Real, w::Real, x::Any, y::Any;
        relative::Bool=false, # if true, o.xpos = parent.x + x
        parent = o.screen.value )
    global widgetStaggerPosx, widgetStaggerPosy
    if typeof( parent ) <: TwScreen
        parentwin = parent.window
        ( parbegy, parbegx ) = getwinbegyx( parentwin )
        ( parmaxy, parmaxx ) = getwinmaxyx( parentwin )
    else
        tmppar = parent
        while( typeof( tmppar.window) <: TwWindow )
            tmppar = tmppar.window.parent.value
        end
        if objtype( tmppar ) == :List
            parmaxy = tmppar.data.canvasheight
            parmaxx = tmppar.data.canvaswidth
        else
            parmaxy = tmppar.height
            parmaxx = tmppar.width
        end
        log( "parmaxy=" * string(parmaxy) * "parmaxx=" * string( parmaxx ) )
        parbegx = parbegy = 0
    end
    if typeof( h ) <: Integer
        o.height = min( h, parmaxy )
    elseif typeof( h ) <: AbstractFloat&& 0.0 < h <= 1.0
        o.height = (@compat round(Int, parmaxy * h ))
        if o.height == 0
            throw( "height is too small")
        end
    else
        throw( "Illegal ysize " * string( h ) )
    end

    if typeof( w ) <: Integer
        o.width = min( w, parmaxx )
    elseif typeof( w ) <: AbstractFloat&& 0.0 < w <= 1.0
        o.width = (@compat round(Int, parmaxx * w ))
        if o.width == 0
            throw( "width is too small")
        end
    else
        throw( "Illegal xsize " * string( w ) )
    end

    if relative
        if typeof( x ) <: Integer && typeof( y ) <: Integer
            (begy, begx) = getwinbegyx( o.window )
            xpos = x+begx
            ypos = y+begy
        else
            throw( "Illegal relative position" )
        end
    else
        xpos = x
        ypos = y
    end

    gapx = max( 0, parmaxx - o.width )
    gapy = max( 0, parmaxy - o.height )
    lastx = parbegx + gapx
    lasty = parbegy + gapy
    if x == :left
        xpos = parbegx
    elseif x == :right
        xpos = parbegx + gapx
    elseif x == :center
        xpos = @compat round(Int, parbegx + gapx / 2 )
    elseif x == :random
        xpos = @compat round(Int, parbegx + gapx * rand() )
    elseif x == :staggered
        if widgetStaggerPosx > gapx
            widgetStaggerPosx = 0
        end
        xpos = parbegx + widgetStaggerPosx
        widgetStaggerPosx += 4
    elseif typeof( x ) <: AbstractFloat&& 0.0 <= x <= 1.0
        xpos = @compat round(Int, parbegx + gapx * x )
    end
    xpos = max( min( xpos, lastx ), parbegx )

    if y == :top
        ypos = parbegy
    elseif y == :bottom
        ypos = parbegy + gapy
    elseif y == :center
        ypos = @compat round(Int, parbegy + gapy / 2 )
    elseif y == :random
        ypos = @compat round(Int, parbegy + gapy * rand() )
    elseif y == :staggered
        if widgetStaggerPosy > gapy
            widgetStaggerPosy = 0
        end
        ypos = parbegy + widgetStaggerPosy
        widgetStaggerPosy += 2
    elseif typeof( y ) <: AbstractFloat&& 0.0 <= y <= 1.0
        ypos = @compat round(Int, parbegy + gapy * y )
    end
    ypos = max( min( ypos, lasty ), parbegy )
    o.xpos = xpos
    o.ypos = ypos
end

#=
function unsetFocus( o::TwObj )
    curs_set( 0 )
    obj.hasfocus = false
    unfocus(obj)
end

function setFocus( o::TwObj )
    obj.hasfocus = true
    focus(obj)
    curs_set(1)
end

function switchFocus( newobj, oldobj )
    if oldobj != newobj
        unsetFocus( oldobj )
        setFocus( newobj )
    end
end
=#

# a general blocking API to make a widget a dialogue
function activateTwObj( o::TwObj, tokens::Any=nothing )
    maxy, maxx = getwinmaxyx( o.window )
    werase( o.window )
    mvwprintw( o.window, maxy>>1, maxx>>1, "%s", "..." )
    wrefresh( o.window )

    draw(o)
    if tokens == nothing #just wait for input
        while true
            update_panels()
            doupdate()
            token = readtoken( o.window )
            status = inject( o, token ) # note that it could be :nochar
            if status == :exit_ok
                return o.value
            elseif status == :exit_nothing # most likely a cancel
                return nothing
            end # default is to continue
        end
    else
        for token in tokens
            update_panels()
            doupdate()
            status = inject( o, token )
            if status == :exit_ok
                return o.value
            elseif status == :exit_nothing # most likely a cancel
                return nothing
            end
        end
        # exhausted all the tokens, no obvious response
        unregisterTwObj( o.screen.value, o )
        return nothing
    end
end

function inject( _::TwObj, k::Any )
    if k == :esc
        return :exit_nothing
    else
        return :pass
    end
end

erase( o::TwObj ) = werase( o.window )

function move( o::TwObj, x, y, relative::Bool, refresh::Bool=false )
    begy, begx = getwinbegyx( o.window )
    alignxy!( o, o.height, o.width, x, y, relative=relative )

    xdiff = o.xpos - begx
    ydiff = o.ypos - begy

    if xdiff == 0 && ydiff == 0
        return
    end
    move_panel( o.panel, o.ypos, o.xpos )
    touchwin( o.screen.value )
    if refresh
        draw(o)
    end
end

function focus( _::TwObj )
end

function unfocus( _::TwObj )
end

refresh( o::TwObj ) = (erase(o);draw(o))

helptext( _::TwObj ) = utf8("")


text:
# The sea floor is getting steeper. Maybe the sleigh keys got carried this way?
#
# A massive school of glowing lanternfish swims past. They must spawn quickly
# to reach such large numbers - maybe exponentially quickly? You should model
# their growth rate to be sure.
#
# Although you know nothing about this specific species of lanternfish, you
# make some guesses about their attributes. Surely, each lanternfish creates a
# new lanternfish once every 7 days.
#
# However, this process isn't necessarily synchronized between every
# lanternfish - one lanternfish might have 2 days left until it creates another
# lanternfish, while another might have 4. So, you can model each fish as a
# single number that represents the number of days until it creates a new
# lanternfish.
#
# Furthermore, you reason, a new lanternfish would surely need slightly longer
# before it's capable of producing more lanternfish: two more days for its
# first cycle.
#
# So, suppose you have a lanternfish with an internal timer value of 3:
#
# - After one day, its internal timer would become 2.
# - After another day, its internal timer would become 1.
# - After another day, its internal timer would become 0.
# - After another day, its internal timer would reset to 6, and it would create
#   a new lanternfish with an internal timer of 8.
# - After another day, the first lanternfish would have an internal timer of 5,
#   and the second lanternfish would have an internal timer of 7.
#
# A lanternfish that creates a new fish resets its timer to 6, not 7 (because 0
# is included as a valid timer value). The new lanternfish starts with an
# internal timer of 8 and does not start counting down until the next day.
#
# Realizing what you're trying to do, the submarine automatically produces a
# list of the ages of several hundred nearby lanternfish (your puzzle input).
# For example, suppose you were given the following list:
#
# 3,4,3,1,2
#
# This list means that the first fish has an internal timer of 3, the second
# fish has an internal timer of 4, and so on until the fifth fish, which has an
# internal timer of 2. Simulating these fish over several days would proceed as
# follows:
#
# Initial state: 3,4,3,1,2
# After  1 day:  2,3,2,0,1
# After  2 days: 1,2,1,6,0,8
# After  3 days: 0,1,0,5,6,7,8
# After  4 days: 6,0,6,4,5,6,7,8,8
# After  5 days: 5,6,5,3,4,5,6,7,7,8
# After  6 days: 4,5,4,2,3,4,5,6,6,7
# After  7 days: 3,4,3,1,2,3,4,5,5,6
# After  8 days: 2,3,2,0,1,2,3,4,4,5
# After  9 days: 1,2,1,6,0,1,2,3,3,4,8
# After 10 days: 0,1,0,5,6,0,1,2,2,3,7,8
# After 11 days: 6,0,6,4,5,6,0,1,1,2,6,7,8,8,8
# After 12 days: 5,6,5,3,4,5,6,0,0,1,5,6,7,7,7,8,8
# After 13 days: 4,5,4,2,3,4,5,6,6,0,4,5,6,6,6,7,7,8,8
# After 14 days: 3,4,3,1,2,3,4,5,5,6,3,4,5,5,5,6,6,7,7,8
# After 15 days: 2,3,2,0,1,2,3,4,4,5,2,3,4,4,4,5,5,6,6,7
# After 16 days: 1,2,1,6,0,1,2,3,3,4,1,2,3,3,3,4,4,5,5,6,8
# After 17 days: 0,1,0,5,6,0,1,2,2,3,0,1,2,2,2,3,3,4,4,5,7,8
# After 18 days: 6,0,6,4,5,6,0,1,1,2,6,0,1,1,1,2,2,3,3,4,6,7,8,8,8,8
#
# Each day, a 0 becomes a 6 and adds a new 8 to the end of the list, while each
# other number decreases by 1 if it was present at the start of the day.
#
# In this example, after 18 days, there are a total of 26 fish. After 80 days,
# there would be a total of 5934.
#
# Find a way to simulate lanternfish. How many lanternfish would there be after
# 80 days?

# read input
fish = parse.(Int8, split(readline("input.txt"), ','))

for _ in 1:80
  fish .-= 1
  replicating = findall(isequal(-1), fish)
  if length(replicating) > 0
    fish[replicating] .= 6
    append!(fish, repeat([8], length(replicating)))
  end
end

println("Number of fish: ", length(fish))

ext: jl


text:
function cartesianmoveindex(direction)
    directiondict = Dict("UP" => (-1, 0), "DOWN" => (1, 0), "LEFT" => (0, -1), "RIGHT" => (0, 1))

    return CartesianIndex(directiondict[direction])
end

function swaptiles!(board, x, y)
    temp = board[x]
    board[x] = board[y]
    board[y] = temp
end

function moves(board)
    first = findfirst(iszero, board)
    movedict = Dict("UP" => "DOWN", "DOWN" => "UP", "LEFT" => "RIGHT", "RIGHT" => "LEFT")

    directions = []
    for i in ("UP", "DOWN", "LEFT", "RIGHT")
        index = first + cartesianmoveindex(i)

        if checkbounds(Bool, board, index)
            action = (index, movedict[i])
            push!(directions, action)
        end
    end

    return directions
end

function copyboard(action, board)
    newboard = copy(board)
    index1 = action[1]
    index2 = index1 + cartesianmoveindex(action[2])
    swaptiles!(newboard, index1, index2)

    return newboard
end

## determine all options for move states
function expansion(board)
    states = []
    possacts = moves(board)

    for act in possacts
        push!(states, copyboard(act, board))
    end

    return states
end

max_stars_repo_name: brice-allen/npuz

max_issues_repo_name: brice-allen/npuz

max_forks_repo_name: brice-allen/npuz


text:
# Generic solve methods
"iLQR solve method"
function solve!(prob::Problem{T,Discrete}, solver::iLQRSolver{T}) where T<:AbstractFloat
    reset!(solver)
    to = solver.stats[:timer]

    n,m,N = size(prob)
    J = Inf

    logger = default_logger(solver)

    # Initial rollout
    rollout!(prob)
    live_plotting(prob,solver)

    J_prev = cost(prob.obj, prob.X, prob.U, get_dt_traj(prob))

    with_logger(logger) do
        record_iteration!(prob, solver, J_prev, Inf)
        for i = 1:solver.opts.iterations
            J = step!(prob, solver, J_prev)

            # check for cost blow up
            if J > solver.opts.max_cost_value
                @warn "Cost exceeded maximum cost"
                return solver
            end

            @timeit to "copy" begin
                copyto!(prob.X, solver.X̄)
                copyto!(prob.U, solver.Ū)
            end

            dJ = abs(J - J_prev)
            J_prev = copy(J)
            record_iteration!(prob, solver, J, dJ)
            live_plotting(prob,solver)

            println(logger, InnerLoop)
            @timeit to "convergence" evaluate_convergence(solver) ? break : nothing
        end
    end
    return solver
end

function step!(prob::Problem{T}, solver::iLQRSolver{T}, J::T) where T
    to = solver.stats[:timer]
    @timeit to "jacobians" jacobian!(prob,solver)
    @timeit to "cost expansion" cost_expansion!(prob,solver)
    @timeit to "backward pass" ΔV = backwardpass!(prob,solver)
    @timeit to "forward pass" forwardpass!(prob,solver,ΔV,J)
end

function cost_expansion!(prob::Problem{T},solver::iLQRSolver{T}) where T
    reset!(solver.Q)
    if solver.opts.square_root
        cost_expansion_sqrt!(solver.Q, prob.obj, prob.X, prob.U, get_dt_traj(prob))
    else
        cost_expansion!(solver.Q, prob.obj, prob.X, prob.U, get_dt_traj(prob))
    end
end

"Plot state, control trajectories"
function live_plotting(prob::Problem{T},solver::iLQRSolver{T}) where T
    if solver.opts.live_plotting == :state
        p = plot(prob.X,title="State trajectory")
        display(p)
    elseif solver.opts.live_plotting == :control
        p = plot(prob.U,title="Control trajectory")
        display(p)
    else
        nothing
    end
end

function record_iteration!(prob::Problem{T}, solver::iLQRSolver{T}, J::T, dJ::T) where T
    solver.stats[:iterations] += 1
    push!(solver.stats[:cost], J)
    push!(solver.stats[:dJ], dJ)
    push!(solver.stats[:gradient],calculate_gradient(prob,solver))
    dJ == 0.0 ? solver.stats[:dJ_zero_counter] += 1 : solver.stats[:dJ_zero_counter] = 0

    @logmsg InnerLoop :iter value=solver.stats[:iterations]
    @logmsg InnerLoop :cost value=J
    @logmsg InnerLoop :dJ   value=dJ
    @logmsg InnerLoop :grad value=solver.stats[:gradient][end]
    @logmsg InnerLoop :zero_count value=solver.stats[:dJ_zero_counter][end]
end

function calculate_gradient(prob::Problem,solver::iLQRSolver)
    if solver.opts.gradient_type == :todorov
        gradient = gradient_todorov(prob,solver)
    elseif solver.opts.gradient_type == :feedforward
        gradient = gradient_feedforward(solver)
    elseif solver.opts.gradient_type == :ℓ2
        gradient = norm(compute_gradient(prob,solver))
    elseif solver.opts.gradient_type == :ℓinf
        gradient = norm(compute_gradient(prob,solver),Inf)
    end
    return gradient
end

function compute_gradient(prob::Problem,solver::iLQRSolver)
    N = prob.N
    reset!(solver.Q)
    cost_expansion!(solver.Q,prob.obj,prob.X,prob.U,get_dt_traj(prob))
    ∇ℓ = []
    for k = 1:N-1
        push!(∇ℓ,solver.Q[k].x)
        push!(∇ℓ,solver.Q[k].u)
    end
    push!(∇ℓ,solver.Q[N].x)

    vcat(∇ℓ...)
end

"""
$(SIGNATURES)
    Calculate the problem gradient using heuristic from iLQG (Todorov) solver
"""
function gradient_todorov(prob::Problem,solver::iLQRSolver)
    N = prob.N
    maxes = zeros(N)
    for k = 1:N-1
        maxes[k] = maximum(abs.(solver.d[k])./(abs.(prob.U[k]).+1))
    end
    mean(maxes)
end

"""
$(SIGNATURES)
    Calculate the infinity norm of the gradient using feedforward term d (from δu = Kδx + d)
"""
function gradient_feedforward(solver::iLQRSolver)
    norm(solver.d,Inf)
end

function evaluate_convergence(solver::iLQRSolver)
    # Check for cost convergence
    # note the  dJ > 0 criteria exists to prevent loop exit when forward pass makes no improvement
    if 0.0 < solver.stats[:dJ][end] < solver.opts.cost_tolerance
        return true
    end

    # Check for gradient convergence
    if solver.stats[:gradient][end] < solver.opts.gradient_norm_tolerance
        return true
    end

    # Check total iterations
    if solver.stats[:iterations] >= solver.opts.iterations
        return true
    end

    # Outer loop update if forward pass is repeatedly unsuccessful
    if solver.stats[:dJ_zero_counter] > solver.opts.dJ_counter_limit
        return true
    end

    return false
end

function regularization_update!(solver::iLQRSolver,status::Symbol=:increase)
    if status == :increase # increase regularization
        # @logmsg InnerLoop "Regularization Increased"
        solver.dρ[1] = max(solver.dρ[1]*solver.opts.bp_reg_increase_factor, solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = max(solver.ρ[1]*solver.dρ[1], solver.opts.bp_reg_min)
        if solver.ρ[1] > solver.opts.bp_reg_max
            @warn "Max regularization exceeded"
        end
    elseif status == :decrease # decrease regularization
        solver.dρ[1] = min(solver.dρ[1]/solver.opts.bp_reg_increase_factor, 1.0/solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = solver.ρ[1]*solver.dρ[1]*(solver.ρ[1]*solver.dρ[1]>solver.opts.bp_reg_min)
    end
end

"Project dynamically infeasible state trajectory into feasible space using TVLQR"
function projection!(prob::Problem{T},opts::iLQRSolverOptions{T}=iLQRSolverOptions{T}()) where T
    # backward pass - project infeasible trajectory into feasible space using time varying lqr
    solver_ilqr = AbstractSolver(prob, opts)
    backwardpass!(prob, solver_ilqr)

    # rollout
    rollout!(prob,solver_ilqr,0.0)

    # update trajectories
    copyto!(prob.X, solver_ilqr.X̄)
    copyto!(prob.U, solver_ilqr.Ū)
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
# Dissipation based homotopy predictor-corrector globalization for solving
# steady problems using Newton's method
# based on Brown and Zingg, "A monolithic Homotopy Continuation Algorithm
# with applications to Computational Fluid Dynamics"
# Journal of Computational Physics 321, (2016), 55-75
# specifically, Algorithm 2

mutable struct HomotopyData{Tsol, Tjac}

  physics_func::Function  # the function that will be solved at lambda = 0
  g_func::Function        # the function to be blended with physics_func

  time::Timings
  lambda::Float64 # homotopy parameter
  myrank::Int

  # parameters
  lambda_min::Float64
  lambda_cutoff::Float64
  itermax::Int   # maximum number of homotopy iterations
  res_reltol::Float64  # residual tolerance for the final solve
  res_abstol::Float64  # residual tolerance for the final solve
  krylov_reltol0::Float64  # Krylov tolerance that should be set before
                           # entering Newton's method
  orig_newton_globalize_euler::Bool  # value to reset opts afterwards
  use_pc::Bool  # true if PC is not PCNone
  homotopy_globalize_euler::Bool
  tighten_early::Bool  # tighten nonlinear tolerance when lambda reaches
                       # lambda_cutoff

  # working variables
  iter::Int  # current major iteration
  homotopy_tol::Float64  # the residual tolerance for homotopy iterations
                         # when lambda > 0
  delta_max::Float64  # step size limit
  psi_max::Float64  # max angle between tangent vectors (radians)
  h_max::Float64  # maximum h value
  psi::Float64  # current angle between tangent vectors (radians)
  tan_norm_1::Float64  # previous tangent vector norm
  res_norm::Float64   # current norm of physics (not homotopy) residual
  h::Float64  # step size

  # arrays
  q_vec0::Array{Tsol, 1}
  delta_q::Array{Tsol, 1}
  tan_vec::Array{Tjac, 1}
  tan_vec_1::Array{Tjac, 1}
  dHdLambda_real::Array{Tjac, 1}
  rhs_vec::Array{Tsol, 1}

  # composite objects
  recalc_policy::RecalculationPolicy
  ls::LinearSolver
  newton_data::NewtonData
  fconv::IO

  # homotopy shock capturing
  sensor_easy::AbstractShockSensor  # a shock sensor that is easy to solve
  opts_orig::Dict{String, Bool}     # options dictionary used for resetting
                                    # the real options dictionary to is
                                    # original state

  function HomotopyData{Tsol, Tjac}(mesh, sbp, eqn, opts,
                           physics_func::Function, g_func::Function
                          ) where {Tsol, Tjac}

    #println("\nConstructing HomotopyData")
    time = eqn.params.time
    lambda = 1.0  # homotopy parameter
    myrank = mesh.myrank

    # some parameters
    lambda_min = 0.0
    lambda_cutoff = 0.000  # was 0.005
    itermax = opts["itermax"]::Int
    res_reltol=opts["res_reltol"]::Float64
    res_abstol=opts["res_abstol"]::Float64
    krylov_reltol0 = opts["krylov_reltol"]::Float64
    orig_newton_globalize_euler = opts["newton_globalize_euler"]  # reset value
    homotopy_globalize_euler = opts["homotopy_globalize_euler"]
    tighten_early = opts["homotopy_tighten_early"]::Bool

    # counters/loop variables
    iter = 1
    homotopy_tol = 1e-2
    delta_max = 1.0  # step size limit, set to 1 initially,
    psi_max = opts["homotopy_psi_max"]
    h_max = 0.1  # maximum h
    psi = 0.0  # angle between tangent vectors
    tan_norm_1 = 0.0  # previous tangent vector norm
    res_norm = 0.0  # norm of residual (not homotopy)
    h = 0.05  # step size
    lambda -= h  # the jacobian is ill-conditioned at lambda=1, so skip it
    #h = lambda/20

    recalc_policy = getRecalculationPolicy(opts, "homotopy")
    # log file
    if myrank == 0
      fconv = BufferedIO("convergence.dat", "a+")
    else
      fconv = DevNull
    end

    # needed arrays
    q_vec0 = zeros(eqn.q_vec)
    delta_q = zeros(eqn.q_vec)
    tan_vec = zeros(Tjac, length(eqn.q_vec))  # tangent vector
    tan_vec_1 = zeros(tan_vec)  # previous tangent vector
    dHdLambda_real = zeros(Tjac, length(eqn.q_vec))

    # stuff for newtonInner
    # because the homotopy function is a pseudo-physics, we can reuse the
    # Newton PC and LO stuff, supplying homotopyPhysics in ctx_residual
    pc, lo = getHomotopyPCandLO(mesh, sbp, eqn, opts)
    use_pc = !(typeof(pc) <: PCNone)

    ls = StandardLinearSolver(pc, lo, eqn.comm, opts)

    # configure NewtonData
    newton_data, rhs_vec = setupNewton(mesh, mesh, sbp, eqn, opts, ls)
    newton_data.itermax = 30

    # homotopy shock capturing
    sensor_easy = createShockSensor(mesh, sbp, eqn, opts,
                                    opts["homotopy_shock_sensor"])
    opts_orig = Dict{String, Bool}()

    obj = new(physics_func, g_func, time, lambda, myrank,
              # parameters
              lambda_min, lambda_cutoff, itermax, res_reltol, res_abstol,
              krylov_reltol0, orig_newton_globalize_euler, use_pc,
              homotopy_globalize_euler, tighten_early,
              # working variables
              iter, homotopy_tol, delta_max, psi_max, h_max, psi,
              tan_norm_1, res_norm, h,
              # arrays
              q_vec0, delta_q, tan_vec, tan_vec_1, dHdLambda_real, rhs_vec,
              # composed objects
              recalc_policy, ls, newton_data, fconv,
              # homotopy shock capturing
              sensor_easy, opts_orig)

    lo.hdata = obj
    setLambda(obj, lambda)

    return obj
  end
end

function setLambda(hdata::HomotopyData, lambda::Number)

  hdata.lambda = lambda

  return nothing
end

"""
  This function solves steady problems using a dissipation-based
  predictor-correcor homotopy globalization for Newtons method.

  Inputs:
    physics_func: the function to solve, ie. func(q) = 0  mathematically.
                  func must have the signature func(mesh, sbp, eqn, opts)

    g_func: the function that evalutes G(q), the dissipation.
            If explicit jacobian calculation is used, then
            [`evalHomotopyJacobian`](@ref) must evaluate the jacobian of this
            function wrt. q
    sbp: an SBP operator
    eqn: a AbstractSolutionData.  On entry, eqn.q_vec must be the
         initial condition.  On exit, eqn.q_vec will be the solution to
         func(q) = 0
    opts: options dictionary

  This function uses Newtons method internally, and supports all the
  different jacobian types and jacobian calculation methods that
  Newton does.

  On entry, eqn.q_vec should contain the initial guess for q.  On exit
  it will contain the solution for func(q) = 0.

  This function is reentrant.

  **Options Keys**

   * calc_jac_explicit
   * itermax
   * res_reltol
   * res_abstol
   * krylov_reltol
   * homotopy_globalize_euler: activate implicit euler globalization when
                               lambda = 0
   * homotopy_tighten_early: tighten the newton solve tolerance and active
                             implicit Euler globalization when
                             `lambda < lambda_cutoff`, typically 0.005.
                             This may slow down the Newton solve significantly
   * homotopy_shock_capturing: see below
   * homotopy_shock_sensor: shock sensor for `R_easy` with
                            `homotopy_shock_capturing`

  This function supports jacobian/preconditioner freezing using the
  prefix "newton".  Note that this recalcuation policy only affects
  this function, and not newtonInner, and defaults to never recalculating
  (and letting newtonInner update the jacobian/preconditioner according to
  its recalculationPolicy).

  **Homotopy Shock Capturing**

  Does homotopy from one shock capturing scheme to another.
  The homotopy function is:

  R_euler (1-lambda)*R_hard + lambda*R_easy

  where R_euler is the physics residual, excluding shock capturing, R_hard
  is the shock sensor that is difficult to converge, and R_easy is the shock
  sensor that is easy to converge.  It is recommended that the initial
  condition satisfy R_euler + R_easy.

  In this mode, `g_func` is ignored.

"""
function predictorCorrectorHomotopy(physics_func::Function,
                  g_func::Function,
                  mesh::AbstractMesh{Tmsh},
                  sbp::AbstractOperator,
                  eqn::AbstractSolutionData{Tsol, Tres},
                  opts) where {Tsol, Tres, Tmsh}

  #----------------------------------------------------------------------------
  # setup

  Tjac = real(Tres)
  hdata = HomotopyData{Tsol, Tjac}(mesh, sbp, eqn, opts, physics_func, g_func)
  myrank = hdata.myrank

  # because the homotopy function is a pseudo-physics, we can reuse the
  # Newton PC and LO stuff, supplying homotopyPhysics in ctx_residual
  homotopyPhysics = (mesh, sbp, eqn, opts, t) -> _homotopyPhysics(hdata, mesh, sbp, eqn, opts, t)
  rhs_func = physicsRhs
  ctx_residual = (homotopyPhysics,)

  # calculate physics residual
  res_norm = real(physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (physics_func,)))
  res_norm_0 = res_norm

  # print to log file
  @mpi_master println(hdata.fconv, 0, " ", res_norm, " ", 0.0)
  @mpi_master flush(hdata.fconv)

  eqn.majorIterationCallback(0, mesh, sbp, eqn, opts, BSTDOUT)

  #----------------------------------------------------------------------------
  # main loop
  while res_norm > res_norm_0*hdata.res_reltol && res_norm > hdata.res_abstol &&
        hdata.iter < hdata.itermax  # predictor loop

    @unpack hdata h lambda iter
    @mpi_master begin
      println(BSTDOUT, "\npredictor iteration ", iter, ", lambda = ", lambda)
      println(BSTDOUT, "res_norm = ", res_norm)
      println(BSTDOUT, "res_norm/res_norm_0 = ", res_norm/res_norm_0)
    end

    # calculate homotopy residual
    #homotopy_norm = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (homotopyPhysics,))

    # update tolerances (linear and nonlinear)
    updateTolerances(hdata, opts)

   # calculate the PC and LO if needed
    doRecalculation(hdata.recalc_policy, hdata.iter, hdata.newton_data.ls,
                    mesh, sbp, eqn, opts, ctx_residual, 0.0)

    # do corrector steps
    copy!(hdata.q_vec0, eqn.q_vec)  # save initial q to calculate delta_q later
    newtonInner(hdata.newton_data, mesh, sbp, eqn, opts, rhs_func, hdata.ls,
                hdata.rhs_vec, ctx_residual)

    # predictor step calculation
    if abs(lambda - hdata.lambda_min) > eps()

      # compute tangent vector
      calcTangentVector(hdata, mesh, sbp, eqn, opts)

      takePredictorStep(hdata, mesh, sbp, eqn, opts)

    end  # end if lambda too large

    # calculate physics residual at new state q
    res_norm = real(physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (physics_func,),))

    # print to log file
    @mpi_master println(hdata.fconv, iter, " ", res_norm, " ", h )
    @mpi_master flush(hdata.fconv)

    eqn.majorIterationCallback(iter, mesh, sbp, eqn, opts, BSTDOUT)

    hdata.iter += 1
  end  # end while loop

  print(BSTDOUT, "\n")

  # inform user of final status
  @mpi_master if hdata.iter >= hdata.itermax
    println(BSTDERR, "Warning: predictor-corrector did not converge in $(hdata.iter) iterations")

  elseif res_norm <= hdata.res_abstol
    println(BSTDOUT, "predictor-corrector converged with absolute residual norm $res_norm")
  elseif res_norm/res_norm_0 <= hdata.res_reltol
    tmp = res_norm/res_norm_0
    println(BSTDOUT, "predictor-corrector converged with relative residual norm $tmp")
  end

  # reset options dictionary
  opts["newton_globalize_euler"] = hdata.orig_newton_globalize_euler

  free(hdata.newton_data)

  flush(BSTDOUT)
  flush(BSTDERR)

  return nothing
end

"""
  This function can be wrapped inside a lambda function to make it look like
  evalResidual.
"""
function _homotopyPhysics(hdata::HomotopyData, mesh, sbp, eqn, opts, t)

  # this function is only for use with Newton's method, where parallel
  # communication is started outside the physics
  # q_vec -> q

  if opts["homotopy_shock_capturing"]
    evalHomotopyFunction_sc(hdata, mesh, sbp, eqn, opts, t)

  else
    fill!(eqn.res, 0.0)
    res_homotopy = zeros(eqn.res)

    # calculate physics residual
    # call this function before g_func, to receive parallel communication
    hdata.physics_func(mesh, sbp, eqn, opts, t)

    # calculate homotopy function
    hdata.g_func(mesh, sbp, eqn, opts, res_homotopy)

    # combine (use lambda from outer function)
    lambda_c = 1 - hdata.lambda # complement of lambda
    for i=1:length(eqn.res)
      eqn.res[i] =  lambda_c*eqn.res[i] + hdata.lambda*res_homotopy[i]
    end

    #println("homotopy physics exiting with residual norm ", norm(vec(eqn.res)))
  end

  return nothing
end

"""
  Updates the linear and nonlinear tolerances for the Newton solve.
  The tolerances usually only change then lambda is close to or at 0.
  This also has the effect of resetting the tolerances if Newton changed
  them (for example, doing inexact Newton-Krylov)

  **Inputs**

   * hdata: `HomotopyData`
   * opts: options dictonary

"""
function updateTolerances(hdata::HomotopyData, opts)

  @unpack hdata res_reltol res_abstol homotopy_globalize_euler tighten_early
  @unpack hdata lambda_cutoff lambda_min homotopy_tol krylov_reltol0 lambda
  @unpack hdata newton_data myrank

  # if we have finished traversing the homotopy path, solve the
  # homotopy problem (= the physics problem because lambda is zero)
  # tightly
  changed_tols = false
  if abs(lambda - lambda_min) <= eps()
    @mpi_master println(BSTDOUT, "tightening homotopy tolerance at lambda = 0")
    changed_tols = true
    homotopy_tol = res_reltol

    reltol = krylov_reltol0 # smaller than newton tolerance
    abstol = 1e-50

    # enable globalization if required
    opts["newton_globalize_euler"] = opts["homotopy_globalize_euler"]
    newton_data.itermax = opts["itermax"]
  elseif lambda < lambda_cutoff && tighten_early

    # turn on implicit euler and tighten tolerances
    # This helps for shock problems when the algorithm takes very small steps
    # in lambda near the end, because lambda may become too small to prevent
    # Newton from diverging
    @mpi_master println(BSTDOUT, "tightening homotopy tolerance at lambda cutoff")
    changed_tols = true
    homotopy_tol = 1e-4

    reltol = krylov_reltol0  # smaller than newton tolerance
    abstol = 1e-50

    # enable globalization if required
    opts["newton_globalize_euler"] = opts["homotopy_globalize_euler"]
    newton_data.itermax = opts["itermax"]
  else  # normal iteration
    reltol = krylov_reltol0  # reset in case Newton did inexact Newton-Krylov
    abstol = -1.0
  end

  if changed_tols
    @mpi_master begin
      println(BSTDOUT, "setting homotopy tolerance to ", homotopy_tol)
      println(BSTDOUT, "ksp reltol = ", reltol)
      println(BSTDOUT, "ksp abstol = ", abstol)
    end
  end

  # record updated tolerances
  newton_data.res_reltol = homotopy_tol
  hdata.homotopy_tol = homotopy_tol

  # reset tolerances in case newton is doing inexact-NK and thus has changed
  # the tolerances
  setTolerances(newton_data.ls, reltol, abstol, -1, -1)

  return nothing
end

"""
  Solve for the tangent vector, overwriting `hdata.tan_vec`

  **Inputs**

   * hdata: `HomotopyData`
   * mesh
   * sbp
   * eqn
   * opts
"""
function calcTangentVector(hdata::HomotopyData, mesh, sbp, eqn, opts)

  @unpack hdata rhs_vec dHdLambda_real tan_vec myrank ls

  # calculate dHdLambda at new q value
  calcdHdLambda(hdata, mesh, sbp, eqn, opts, rhs_vec)
  for i=1:length(rhs_vec)
    dHdLambda_real[i] = real(rhs_vec[i])
  end

  # calculate tangent vector dH/dq * t = dH/dLambda
  @mpi_master println(BSTDOUT, "solving for tangent vector")
  # the accuracy of the tangent vector is not critical, so re-use the
  # Jacobian from the last Newton iteration
  #calcPCandLO(ls, mesh, sbp, eqn, opts, ctx_residual, 0.0)

  tsolve = @elapsed linearSolve(ls, dHdLambda_real, tan_vec)
  eqn.params.time.t_solve += tsolve

  return nothing
end

"""
  This function calculates dH/dLambda, where H is the homotopy function
  calculated by homotopyPhysics.  The differentiation is done analytically.

  **Inputs**

   * hdata: `HomotopyData`
   * mesh
   * sbp
   * eqn: eqn.res and eqn.res_vec are overwritten
   * opts

  Inputs/Outputs
    res_vec: vector to store dH/dLambda in

  Aliasing restrictions: res_vec and eqn.res_vec may not alias
"""
function calcdHdLambda(hdata::HomotopyData, mesh, sbp, eqn, opts, res_vec)

  @unpack hdata physics_func g_func

  # it appears this only gets called after parallel communication is done
  # so no need to start communication here

  if opts["homotopy_shock_capturing"]
    calcdHdLambda_sc(hdata, mesh, sbp, eqn, opts, 0.0, res_vec)
  else
    res_homotopy = zeros(eqn.res)

    # calculate physics residual
    physics_func(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

    # calculate homotopy function
    g_func(mesh, sbp, eqn, opts, res_homotopy)
    array3DTo1D(mesh, sbp, eqn, opts, res_homotopy, res_vec)

    # combine them
    for i=1:length(res_vec)
      res_vec[i] -= eqn.res_vec[i]
    end
  end

  return nothing
end

"""
  This function takes the predictor step.  This requires a few actions:

  * normalize tangent vector (overwrites `hdata.tan_vec`)
  * compute `delta_q` (overwrites `hdata.delta_q` and uses `hdata.q_vec0`
  * applies the various heuristics to update the step size (step length, angle,
    etc.)

"""
function takePredictorStep(hdata::HomotopyData, mesh, sbp, eqn, opts)

  @unpack hdata tan_vec tan_vec_1 tan_norm_1 delta_q q_vec0 h lambda
  @unpack hdata psi_max delta_max h_max lambda_min

  # normalize tangent vector
  # t = [z, -1], where z is a delta_q sized vector, so the norm is
  # sqrt(calcNorm(z)^2 + (1)^2).  In the code tan_vec = z, so what is
  # really happenening is the z component of t is being normalize.  The
  # -1 component will be handled below
  tan_norm = calcNorm(eqn, tan_vec)
  tan_norm = sqrt(tan_norm*tan_norm + 1)
  scale!(tan_vec, 1/tan_norm)

  # compute delta_q
  for i=1:length(eqn.q)
    delta_q[i] = eqn.q_vec[i] - q_vec0[i]
  end
  delta = calcNorm(eqn, delta_q)
  #println("L2 norm delta = ", delta)

  psi = psi_max
  if hdata.iter > 1
    # compute phi = acos(tangent_i_1 dot tangent_i)
    # however, use the L2 norm for the z part of the tangent vector
    tan_term = calcL2InnerProduct(eqn, tan_vec_1, tan_vec)
    # now add the normalized -1 components of t_i and t_i_1 to complete
    # the inner product between t_i and t_i_1 = inner_product(z_i, z_i_1)
    # + -1*-1/(||z_i||*||z_i_1||)
    tan_norm_term = (1/tan_norm)*(1/tan_norm_1)
    arg = tan_term + tan_norm_term
    arg = clamp(arg, -1.0, 1.0)
    psi = acos( arg )
  end

  # calculate step size
  fac = max(real(psi/psi_max), sqrt(delta/delta_max))
  #println("psi/psi_max = ", psi/psi_max)
  #println("delta/delta_max = ", delta/delta_max)
  #println("fac = ", fac)
  h /= fac
  h = min(h, h_max)
  lambda = max(lambda_min, lambda - h)
  if lambda < eps()
    lambda = 0.0
  end
  #println("new h = ", h)

  # take predictor step
  for i=1:length(eqn.q_vec)
    eqn.q_vec[i] += h*tan_vec[i]
  end

  # save variables
  setLambda(hdata, lambda)
  copy!(tan_vec_1, tan_vec)
  hdata.tan_norm_1 = tan_norm
  hdata.h = h

  return nothing
end

function getHomotopyPCandLO(mesh, sbp, eqn, opts)

  # get PC
  if opts["jac_type"] <= 2
    pc = PCNone(mesh, sbp, eqn, opts)
  else
    pc = HomotopyMatPC(mesh, sbp, eqn, opts)
  end

  jactype = opts["jac_type"]
  if jactype == 1
    lo = HomotopyDenseLO(pc, mesh, sbp, eqn, opts)
  elseif jactype == 2
    lo = HomotopySparseDirectLO(pc, mesh, sbp, eqn, opts)
  elseif jactype == 3
    lo = HomotopyPetscMatLO(pc, mesh, sbp, eqn, opts)
  elseif jactype == 4
    lo = HomotopyPetscMatFreeLO(pc, mesh, sbp, eqn, opts)
  end

  return pc, lo
end

# Define PC and LO objects
#------------------------------------------------------------------------------
# PC:

mutable struct HomotopyMatPC <: AbstractPetscMatPC
  pc_inner::NewtonMatPC
  lambda::Float64  # homotopy parameter
end

function HomotopyMatPC(mesh::AbstractMesh, sbp::AbstractOperator,
                    eqn::AbstractSolutionData, opts::Dict)

  pc_inner = NewtonMatPC(mesh, sbp, eqn, opts)
  lambda = 1.0

  return HomotopyMatPC(pc_inner, lambda)
end

function calcPC(pc::HomotopyMatPC, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  # compute the Jacobian of the Newton PC
  calcPC(pc.pc_inner, mesh, sbp, eqn, opts, ctx_residual, t)

  if opts["calc_jac_explicit"]
    A = getBasePC(pc).A
    assembly_begin(A, MAT_FINAL_ASSEMBLY)
    assembly_end(A, MAT_FINAL_ASSEMBLY)

    lambda_c = 1 - lambda # complement of lambda
    scale!(A, lambda_c)

    # compute the homotopy contribution to the Jacobian
    assembler = _AssembleElementData(A, mesh, sbp, eqn, opts)
    evalHomotopyJacobian(mesh, sbp, eqn, opts, assembler, lambda)
  end

  return nothing
end

#------------------------------------------------------------------------------
# LO:

mutable struct HomotopyDenseLO <: AbstractDenseLO
  lo_inner::NewtonDenseLO
  hdata::HomotopyData

  function HomotopyDenseLO(lo_inner::NewtonDenseLO)
    return new(lo_inner)
  end
end

function HomotopyDenseLO(pc::PCNone, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonDenseLO(pc, mesh, sbp, eqn, opts)
  return HomotopyDenseLO(lo_inner)
end

mutable struct HomotopySparseDirectLO <: AbstractSparseDirectLO
  lo_inner::NewtonSparseDirectLO
  hdata::HomotopyData

  function HomotopySparseDirectLO(lo_inner::NewtonSparseDirectLO)
    return new(lo_inner)
  end

end

function HomotopySparseDirectLO(pc::PCNone, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonSparseDirectLO(pc, mesh, sbp, eqn, opts)

  return HomotopySparseDirectLO(lo_inner)
end

mutable struct HomotopyPetscMatLO <: AbstractPetscMatLO
  lo_inner::NewtonPetscMatLO
  hdata::HomotopyData

  function HomotopyPetscMatLO(lo_inner::NewtonPetscMatLO)
    return new(lo_inner)
  end

end

function HomotopyPetscMatLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonPetscMatLO(pc, mesh, sbp, eqn, opts)

  return HomotopyPetscMatLO(lo_inner)
end

mutable struct HomotopyPetscMatFreeLO <: AbstractPetscMatFreeLO
  lo_inner::NewtonPetscMatFreeLO
  hdata::HomotopyData

  function HomotopyPetscMatFreeLO(lo_inner::NewtonPetscMatFreeLO)
    return new(lo_inner)
  end

end

"""
  Homotopy mat-free linear operator constructor

  **Inputs**

   * pc
   * mesh
   * sbp
   * eqn
   * opts
   * rhs_func: rhs_func from [`newtonInner`](@ref)
"""
function HomotopyPetscMatFreeLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonPetscMatFreeLO(pc, mesh, sbp, eqn, opts)

  return HomotopyPetscMatFreeLO(lo_inner)
end

#------------------------------------------------------------------------------
# Required functions

"""
  Homotopy matrix-explicit linear operators
"""
const HomotopyMatLO = Union{HomotopyDenseLO, HomotopySparseDirectLO, HomotopyPetscMatLO}

function calcLinearOperator(lo::HomotopyMatLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  # ctx_residual[1] = is the *homotopy* function H, so for the non-explict
  # case, calling calcLinearOperator(lo.lo_inner) is fine because it will
  # do coloring

  if opts["homotopy_shock_capturing"] && opts["calc_jac_explicit"]
    evalJacobian_homotopy_sc(lo, mesh, sbp, eqn, opts, ctx_residual, t)
  else

    calcLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t)

    if opts["calc_jac_explicit"]
      A = getBaseLO(lo).A
      assembly_begin(A, MAT_FINAL_ASSEMBLY)
      assembly_end(A, MAT_FINAL_ASSEMBLY)

      lambda_c = 1 - lo.hdata.lambda # complement of lambda
      scale!(A, lambda_c)

      # compute the homotopy contribution to the Jacobian
      assembler = _AssembleElementData(A, mesh, sbp, eqn, opts)
      evalHomotopyJacobian(mesh, sbp, eqn, opts, assembler, lo.hdata.lambda)
    end
  end

  return nothing
end

function calcLinearOperator(lo::HomotopyPetscMatFreeLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  if opts["homotopy_shock_capturing"]
    error("homotopy shock capturing not supported for matrix-free")
  else
    calcLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t)

    setLOCtx(lo, mesh, sbp, eqn, opts, ctx_residual, 0.0)

    # nothing to do here
  end

  return nothing
end

function applyLinearOperator(lo::HomotopyPetscMatFreeLO, mesh::AbstractMesh,
                       sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol},
                       opts::Dict, ctx_residual, t, x::AbstractVector,
                       b::AbstractVector) where Tsol

  @assert !(Tsol <: AbstractFloat)  # complex step only!

  if opts["homotopy_shock_capturing"]
    error("homotopy shock capturing not supported for matrix-free")
  else
    # ctx_residual[1] = homotopyPhysics, so this computes both the physics
    # and homotopy contribution
    applyLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t, x, b)
  end

  return nothing
end

function applyLinearOperatorTranspose(lo::HomotopyPetscMatFreeLO,
                             mesh::AbstractMesh,
                             sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol},
                             opts::Dict, ctx_residual, t, x::AbstractVector,
                             b::AbstractVector) where Tsol

  error("applyLinearOperatorTranspose() not supported by HomotopyPetscMatFreeLO")

end

include("shock_homotopy.jl")


text:
using Distributions
using DataStructures
using DataFrames
using Test

import GlobalSensitivityAnalysis: ishigami, split_output, first_order, total_order, sample
#
# STEP 1. try running confidence intervals a few times and looking at the means etc.
# and then compare that to doing the same for Python SALib
#

#
# STEP 1A. Uniform Case
#

# define problem and some constants
N = 1_000
num_resamples = 10_000
conf_level = 0.95
calc_second_order = false
iter = 10

data = SobolData(
    params = OrderedDict(:x1 => Uniform(-3.14159265359, 3.14159265359),
        :x2 => Uniform(-3.14159265359, 3.14159265359),
        :x3 => Uniform(-3.14159265359, 3.14159265359)),
    N = N,
    calc_second_order = calc_second_order
)

D = length(data.params)

# run sampling since this is deterministic
samples = sample(data) # generate samples using Sobol sequence
Y = ishigami(samples) # run model (example)

# iteratively perform Sobol Analysis and save confidence intervals
firstorder_all = zeros(iter, D)
totalorder_all = zeros(iter, D)

for i = 1:iter
    results = analyze(data, Y; num_resamples = num_resamples)
    firstorder_all[i,:] = results[:firstorder_conf]
    totalorder_all[i,:] = results[:totalorder_conf]
end

@test mean(firstorder_all, dims = 1) ≈ [0.06369217 0.05347996 0.05821841] atol = 1e-3
@test mean(totalorder_all, dims = 1) ≈ [0.08792231 0.04021016 0.02634918] atol = 1e-3

#
# STEP 1B. Non-Uniform Case
#

# define problem and some constants
N = 1_000
num_resamples = 10_000
conf_level = 0.95
calc_second_order = false
iter = 10

data = SobolData(
    params = OrderedDict(:x1 => Normal(1, 0.2),
        :x2 => Uniform(0.75, 1.25),
        :x3 => LogNormal(0, 0.5)),
    N = N,
    calc_second_order = calc_second_order
)

D = length(data.params)

# run sampling since this is deterministic
samples = sample(data) # generate samples using Sobol sequence
Y = ishigami(samples) # run model (example)

# iteratively perform Sobol Analysis and save confidence intervals
firstorder_all = zeros(iter, D)
totalorder_all = zeros(iter, D)

for i = 1:iter
    results = analyze(data, Y; num_resamples = num_resamples)
    firstorder_all[i,:] = results[:firstorder_conf]
    totalorder_all[i,:] = results[:totalorder_conf]
end

@test mean(firstorder_all, dims = 1) ≈ [0.00638154 0.0948373  0.67803243] atol = 1e-1
@test mean(totalorder_all, dims = 1) ≈ [0.01245481 0.09457124 0.10595523] atol = 1e-1

#
# STEP 2. Deterministic example
#

# define problem and constants
N = 5
conf_level = 0.95
calc_second_order = false

data = SobolData(
    params = OrderedDict(:x1 => Uniform(-3.14159265359, 3.14159265359),
        :x2 => Uniform(-3.14159265359, 3.14159265359),
        :x3 => Uniform(-3.14159265359, 3.14159265359)),
    N = N,
    calc_second_order = calc_second_order
)
D = length(data.params)

# sampling
samples = sample(data)
Y = ishigami(samples)

# analysis
A, B, AB, BA = split_output(Matrix(Y), N, D, calc_second_order)

# confidence intervals  - constants
r = [3 2 3;
    4 1 1;
    3 2 0;
    1 2 2;
    1 1 1]
r = r .+ 1
Z = quantile(Normal(0.0, 1.0),1 - (1 - conf_level)/2) # calculate z* for CI

# confidence intervals - calculate
firstorder_conf = Array{Float64}(undef, D)
totalorder_conf = Array{Float64}(undef, D)

for i in 1:D
    firstorder_conf[i] = Z * std(first_order(A[r], AB[r, i], B[r]))
    totalorder_conf[i] = Z * std(total_order(A[r], AB[r, i], B[r]))
end

@test firstorder_conf ≈ [1.57157685, 0.61660013, 0.81296287] atol = 1e-8
@test totalorder_conf ≈ [2.26661086, 0.36204958, 0.46523933] atol = 1e-8

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

include("make_schema_docs.jl")
include("find_attr_groups.jl")
AttrGroups.main()


text:
using Documenter, FastTransforms, Literate, Plots

plotlyjs()

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

examples = [
    "chebyshev.jl",
    "disk.jl",
    "nonlocaldiffusion.jl",
    "padua.jl",
    "sphere.jl",
    "spinweighted.jl",
    "subspaceangles.jl",
    "triangle.jl",
]

function uncomment_objects(str)
    str = replace(str, "###```@raw" => "```\n\n```@raw")
    str = replace(str, "###<object" => "<object")
    str = replace(str, "###```\n```" => "```")
    str
end

for example in examples
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR; execute=true, postprocess = uncomment_objects)
end

makedocs(
            doctest = false,
            format = Documenter.HTML(),
            sitename = "FastTransforms.jl",
            authors = "Richard Mikael Slevinsky",
            pages = Any[
                    "Home" => "index.md",
                    "Examples" => [
                        "generated/chebyshev.md",
                        "generated/disk.md",
                        "generated/nonlocaldiffusion.md",
                        "generated/padua.md",
                        "generated/sphere.md",
                        "generated/spinweighted.md",
                        "generated/subspaceangles.md",
                        "generated/triangle.md",
                        ],
                    ]
        )

deploydocs(
    repo   = "github.com/JuliaApproximation/FastTransforms.jl.git",
    )


text:
"""
Construct vectorfields and derivatives from bases and coefficients.
"""

function construct_vectorfield(coeff, B)
    @assert length(coeff) == length(B)
    a = copy(coeff)
    function vectorfield(x)
        out = [0., 0.]
        for i in 1:length(B)
            out += a[i]*B[i](x)
        end
        return out
    end
    return vectorfield
end

function construct_jacobian_matrix(coeff, DB)
    @assert length(coeff) == length(DB)
    a = copy(coeff)
    function jacobian_matrix(x)
        out = [ [0. 0.]; [0. 0.]]
        for i in 1:length(DB)
            out += a[i]*DB[i](x)
        end
        return out
    end
    return jacobian_matrix
end

function jacobian_determinant(Df)
    function J(x)
        return det(Df(x))
    end
end

function gradient_jacobian_determinant(Df)
    J = jacobian_determinant(Df)
    function DJ(x)
        return ForwardDiff.gradient(J, x)
    end
    return DJ
end


text:
using StaticArrays, Interpolations
using Images, CoordinateTransformations, Rotations
using RegisterMismatch
using RegisterQD
using Random

#import BlockRegistration, RegisterOptimize
#using RegisterCore, RegisterPenalty, RegisterDeformation, RegisterMismatch, RegisterFit

using Test, TestImages

@testset "QuadDIRECT tests with random images" begin
    Random.seed!(879)

    ##### Translations
    #2D
    moving = rand(50,50)
    tfm0 = Translation(-4.7, 5.1) #ground truth
    newfixed = warp(moving, tfm0)
    itp = interpolate(newfixed, BSpline(Linear()))
    etp = extrapolate(itp, NaN)
    fixed = etp(Base.axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation
    thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))]))
    mxshift = (10,10)

    tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, rtol=0)

    @test sum(abs.(tfm0.translation - tfm.translation)) < 1e-3

    #again with cropping
    tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, crop=true, rtol=0)
    @test sum(abs.(tfm0.translation - tfm.translation)) < 1e-3

    #3D
    moving = rand(30,30,30)
    tfm0 = Translation(-0.9, 2.1,1.2) #ground truth
    newfixed = warp(moving, tfm0)
    itp = interpolate(newfixed, BSpline(Linear()))
    etp = extrapolate(itp, NaN)
    fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation
    thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))]))
    mxshift = (5,5,5)

    tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, rtol=0)

    @test sum(abs.(tfm0.translation - tfm.translation)) < 0.1

    ######Rotations + Translations
    #2D
    moving = centered(rand(50,50))
    tfm0 = Translation(-4.0, 5.0) ∘ LinearMap(RotMatrix(pi/360)) #ground truth
    newfixed = warp(moving, tfm0)
    itp = interpolate(newfixed, BSpline(Linear()))
    etp = extrapolate(itp, NaN)
    fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation
    thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))]))
    mxshift = (10,10)
    mxrot = pi/90
    minwidth_rot = [0.0002]
    SD = SDiagonal(@SVector(ones(ndims(fixed))))

    tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, thresh=thresh, maxevals=1000, rtol=0, fvalue=1e-8)
    tfm = arrayscale(tfm, SD)
    @test sum(abs.(tfm0.linear - tfm.linear)) < 1e-3

    #3D
    moving = centered(rand(30,30,30))
    tfm0 = Translation(-1.0, 2.1,1.2) ∘ LinearMap(RotXYZ(pi/360, pi/180, pi/220)) #ground truth
    newfixed = warp(moving, tfm0)
    itp = interpolate(newfixed, BSpline(Linear()))
    etp = extrapolate(itp, NaN)
    fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation
    thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))]))
    mxshift = (5,5,5)
    mxrot = [pi/90; pi/90; pi/90]
    minwidth_rot = fill(0.0002, 3)
    SD = SDiagonal(@SVector(ones(ndims(fixed))))

    tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, thresh=thresh, maxevals=1000, rtol=0)
    tfm = arrayscale(tfm, SD)
    @test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(RotXYZ(tfm.linear)[:], tfm.translation))) < 0.1

#NOTE: the 2D test below fails rarely and the 3D test fails often, apparently because full affine is too difficult with these images
    #####General Affine Transformations
    #2D
    moving = centered(rand(50,50))
    shft = SArray{Tuple{2}}(rand(2).+2.0)
    #random displacement from the identity matrix
    mat = SArray{Tuple{2,2}}([1 0; 0 1] + rand(2,2)./40 + -rand(2,2)./40)
    tfm0 = AffineMap(mat, shft) #ground truth
    newfixed = warp(moving, tfm0)
    itp = interpolate(newfixed, BSpline(Linear()))
    etp = extrapolate(itp, NaN)
    fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation
    thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))]))
    mxshift = (5,5)
    SD = SDiagonal(@SVector(ones(ndims(fixed))))

    tfm, mm = qd_affine(fixed, moving, mxshift; SD=SD, thresh=thresh, rtol=0, fvalue=1e-3, maxevals=10^4)
    tfm = arrayscale(tfm, SD)
    @test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1

        #The tests below fail.  Probably two factors contributing to failure:
        # 1) Full affine 3D is a lot of parameters so it's just difficuilt (12 parameters)
        # 2) The current way of computing mismatch may be flawed for scaling transformations
        #    because the algorithm output shows it tends to compress/expand the image in order to remove
        #    pixels from the denominator in the mismatch calculation (especially bad for small images).
        #3D
        #moving = centered(rand(20,20,20));
        #shft = SArray{Tuple{3}}(2.6, 0.1, -3.3);
        ##random displacement from the identity matrix
        #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./30 + -rand(3,3)./30);
        #tfm0 = AffineMap(mat, shft); #ground truth
        #newfixed = warp(moving, tfm0);
        #inds = intersect.(axes(moving), axes(newfixed))
        #fixed = newfixed[inds...]
        #moving = moving[inds...]
        #thresh = 0.1 * (sum(abs2.(fixed[.!(isnan.(fixed))]))+sum(abs2.(moving[.!(isnan.(moving))])));
        #mxshift = (10,10,10)
        #SD = eye(ndims(fixed));

        #tfm, mm = qd_affine(centered(fixed), centered(moving), mxshift, SD; thresh=thresh, rtol=0, fvalue=1e-4);

        #@test mm <= 1e-4
        #@test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1

        #not random
        #moving = zeros(10,10,10);
        #moving[5:7, 5:7, 5:7] = 1.0
        #moving = centered(moving)
        #shft = SArray{Tuple{3}}(0.6, 0.1, -0.3);
        ##random displacement from the identity matrix
        #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./50 + -rand(3,3)./50);
        #tfm00 = AffineMap(mat, shft);
        ##tfm0 = recenter(tfm00, center(moving)); #ground truth
        #tfm0 = tfm00 #ground truth
        #newfixed = warp(moving, tfm0);
        #inds = intersect.(axes(moving), axes(newfixed))
        #fixed = newfixed[inds...]
        #moving = moving[inds...]
        #thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))]));
        #mxshift = (5,5,5)
        #SD = eye(ndims(fixed));
        #@test RegisterOptimize.aff(vcat(tfm00.translation[:], tfm00.linear[:]), fixed, SD) == tfm0

        #tfm, mm = qd_affine(centered(fixed), centered(moving), mxshift, SD; thresh=thresh, rtol=0, fvalue=1e-4);

        #@test mm <= 1e-4
        #@test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1
    #TODO this test passes if run individually, but breaks when included in the testset.
end #tests with random images

@testset "Suppression of printing" begin
    a, b = rand(5, 5), rand(5, 5)
    ca, cb = centered(a), centered(b)
    mxshift = (2,2)
    mxrot = 0.5
    # minwidth_rot = 1e-3
    mktemp() do path, io
        redirect_stdout(io) do
            qd_translate(a, b, (2,2); print_interval=typemax(Int))
            qd_rigid(ca, cb, mxshift, mxrot; print_interval=typemax(Int))
            qd_affine(ca, cb, mxshift; print_interval=typemax(Int))
        end
        flush(io)
        str = read(path, String)
        @test isempty(str)
    end
end


text:
"""
**SeisPlotCoordinates**

*Plot shot and receiver coordinates or fold map. Optionally can print to a file. *

**IN**
* h: headers (a vector where each element is of type Header)

* style="sxsygxgy" or "fold"
* cmap="Greys"
* aspect="auto"
* pclip=98
* vmin="NULL"
* vmax="NULL"
* title=" "
* xlabel=" "
* xunits=" "
* ylabel=" "
* yunits=" "
* ox=0
* dx=1
* oy=0
* dy=1
* dpi=100
* wbox=6
* hbox=6
* name="NULL"
* interpolation="none"
* titlesize=20
* labelsize=15
* ticksize=10
* fignum="NULL"

**OUT**

*Credits: Aaron Stanton, 2015*

"""
function SeisPlotCoordinates(headers;style="sxsygxgy",cmap="Greys",aspect="auto",pclip=98,vmin="NULL",vmax="NULL",title=" ",xlabel=" ",xunits=" ",ylabel=" ",yunits=" ",ox=0,dx=1,oy=0,dy=1,dpi=100,wbox=6,hbox=6,name="NULL",interpolation="none",titlesize=20,labelsize=15,ticksize=10,fignum="NULL")

	plt.ion()
	if (fignum == "NULL")
		fig = plt.figure(figsize=(wbox, hbox), dpi=dpi, facecolor="w", edgecolor="k")
	else
		fig = plt.figure(num=fignum,figsize=(wbox, hbox), dpi=dpi, facecolor="w", edgecolor="k")
	end
	if (style == "sxsygxgy")
		sx = Seismic.ExtractHeader(headers,"sx")
		sy = Seismic.ExtractHeader(headers,"sy")
		gx = Seismic.ExtractHeader(headers,"gx")
		gy = Seismic.ExtractHeader(headers,"gy")
		im = plt.scatter(gx,gy,marker="^",s=5,color="b");
		im = plt.scatter(sx,sy,marker="*",s=8,color="r");
		#plt.axis([ox,ox + (size(in,2)-1)*dx,oy + (size(in,1)-1)*dy,oy])
	elseif (style == "fold")
		imx = Seismic.ExtractHeader(headers,"imx")
		imy = Seismic.ExtractHeader(headers,"imy")
		imx_min = minimum(imx)
		imx_max = maximum(imx)
		nmx = imx_max - imx_min + 1
		imy_min = minimum(imy)
		imy_max = maximum(imy)
		nmy = imy_max - imy_min + 1
		fold = zeros(nmy,nmx)
		for itrace = 1 : length(imx)
			imx_index = imx[itrace] - imx_min + 1
			imy_index = imy[itrace] - imy_min + 1
			fold[imy_index,imx_index] += 1
		end
		if (vmin=="NULL" || vmax=="NULL")
			if (pclip<=100)
				a = -quantile(abs(fold[:]),(pclip/100))
			else
				a = -quantile(abs(fold[:]),1)*pclip/100
			end
			b = -a
		else
			a = vmin
			b = vmax
		end
		im = plt.imshow(fold,cmap=cmap,vmin=a,vmax=b,extent=[ox,ox + (size(fold,2)-1)*dx,oy + (size(fold,1)-1)*dy,oy],aspect=aspect)
		ax = plt.gca()
		ax[:invert_yaxis]();
		plt.colorbar()
	else
		println("no other plotting style defined yet")
	end

	plt.title(title)
	plt.xlabel(join([xlabel " " xunits]))
	plt.ylabel(join([ylabel " " yunits]))

	if (name == "NULL")
		plt.show()
	else
		plt.savefig(name,dpi=dpi)
		plt.close()
	end
	return im
end


text:
struct TA_Gaugefields_4D_mpi{NC,NumofBasis} <: TA_Gaugefields_4D{NC}
    a::Array{Float64,5}
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NC::Int64
    NumofBasis::Int64
    generators::Union{Nothing,Generator}

    PEs::NTuple{4,Int64}
    PN::NTuple{4,Int64}
    mpiinit::Bool
    myrank::Int64
    nprocs::Int64
    myrank_xyzt::NTuple{4,Int64}

    function TA_Gaugefields_4D_mpi(u::Gaugefields_4D{NC}) where NC
        NumofBasis = ifelse(NC == 1,1,NC^2-1)
        if NC <= 3
            generators = nothing
        else
            generators = Generator(NC)
        end

        return new{NC,NumofBasis}(zeros(Float64,NumofBasis,u.PN[1],u.PN[2],u.PN[3],u.PN[4]),u.NX,u.NY,u.NZ,u.NT,u.NC,NumofBasis,generators,
            u.PEs,u.PN,true,u.myrank,u.nprocs,u.myrank_xyzt
        )
    end
end

function Base.setindex!(x::T,v,i...)  where T<: TA_Gaugefields_4D_mpi
    @inbounds x.a[i...] = v
end

function Base.getindex(x::T,i...) where T<: TA_Gaugefields_4D_mpi
    @inbounds return x.a[i...]
end

function barrier(x::TA_Gaugefields_4D_mpi)
    MPI.Barrier(comm)
end

function add_U!(c::TA_Gaugefields_4D_mpi{NC,NumofBasis},α::N,a::TA_Gaugefields_4D_mpi{NC,NumofBasis}) where {NC, N<:Number,NumofBasis}
    NT = c.NT
    NZ = c.NZ
    NY = c.NY
    NX = c.NX
    #NumofBasis = c.NumofBasis
    for it=1:c.PN[4]
        for iz=1:c.PN[3]
            for iy=1:c.PN[2]
                for ix=1:c.PN[1]
                    for k=1:NumofBasis
                        c[k,ix,iy,iz,it] = c[k,ix,iy,iz,it] + α*a[k,ix,iy,iz,it]
                    end
                end
            end
        end
    end

    barrier(c)
    #error("add_U! is not implemented in type $(typeof(c)) ")
end

function clear_U!(U::TA_Gaugefields_4D_mpi{NC,NumofBasis}) where {NC,NumofBasis}
    #NumofBasis = Uμ.NumofBasis
    for it=1:U.PN[4]
        for iz=1:U.PN[3]
            for iy=1:U.PN[2]
                for ix=1:U.PN[1]
                    for k=1:NumofBasis
                        @inbounds U[k,ix,iy,iz,it] = 0
                    end
                end
            end
        end
    end
    barrier(U)
end

function gauss_distribution!(p::TA_Gaugefields_4D_mpi{NC,NumofBasis};σ=1.0) where {NC,NumofBasis}
    d = Normal(0.0, σ)
    #NumofBasis = Uμ.NumofBasis
    pwork = rand(d,prod(p.PN)*NumofBasis)
    icount = 0
    @inbounds for it=1:p.PN[4]
        for iz=1:p.PN[3]
            for iy=1:p.PN[2]
                @simd for ix=1:p.PN[1]
                    for k=1:NumofBasis
                        icount += 1
                        p[k,ix,iy,iz,it] = pwork[icount]
                    end
                end
            end
        end
    end
    barrier(p)
end

function Base.:*(x::TA_Gaugefields_4D_mpi{NC,NumofBasis},y::TA_Gaugefields_4D_mpi{NC,NumofBasis}) where {NC,NumofBasis}
    #NumofBasis = Uμ.NumofBasis
    s = 0.0
    @inbounds for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                @simd for ix=1:x.PN[1]
                    for k=1:NumofBasis
                        s += x[k,ix,iy,iz,it]*y[k,ix,iy,iz,it]
                    end
                end
            end
        end
    end
    s = MPI.Allreduce(s,MPI.SUM,comm)
    barrier(x)

    return s
end

function Traceless_antihermitian_add!(c::TA_Gaugefields_4D_mpi{3,NumofBasis},factor,vin::Union{Gaugefields_4D_wing_mpi{3},Gaugefields_4D_nowing_mpi{3}}) where NumofBasis
    #error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
    fac13 = 1/3
    NX = vin.NX
    NY = vin.NY
    NZ = vin.NZ
    NT = vin.NT

    for it=1:vin.PN[4]
        for iz=1:vin.PN[3]
            for iy=1:vin.PN[2]
                @simd for ix=1:vin.PN[1]
                    v11 = getvalue(vin,1,1,ix,iy,iz,it)
                    v22 = getvalue(vin,2,2,ix,iy,iz,it)
                    v33 = getvalue(vin,3,3,ix,iy,iz,it)

                    tri = fac13*(imag(v11)+imag(v22)+imag(v33))

                    #=
                    vout[1,1,ix,iy,iz,it] = (imag(v11)-tri)*im
                    vout[2,2,ix,iy,iz,it] = (imag(v22)-tri)*im
                    vout[3,3,ix,iy,iz,it] = (imag(v33)-tri)*im
                    =#
                    y11 = (imag(v11)-tri)*im
                    y22 = (imag(v22)-tri)*im
                    y33 = (imag(v33)-tri)*im

                    v12 = getvalue(vin,1,2,ix,iy,iz,it)
                    v13 = getvalue(vin,1,3,ix,iy,iz,it)
                    v21 = getvalue(vin,2,1,ix,iy,iz,it)
                    v23 = getvalue(vin,2,3,ix,iy,iz,it)
                    v31 = getvalue(vin,3,1,ix,iy,iz,it)
                    v32 = getvalue(vin,3,2,ix,iy,iz,it)

                    x12 = v12 - conj(v21)
                    x13 = v13 - conj(v31)
                    x23 = v23 - conj(v32)

                    x21 = - conj(x12)
                    x31 = - conj(x13)
                    x32 = - conj(x23)

                    #=
                    vout[1,2,ix,iy,iz,it) = 0.5  * x12
                    vout[1,3,ix,iy,iz,it) = 0.5  * x13
                    vout[2,1,ix,iy,iz,it) = 0.5  * x21
                    vout[2,3,ix,iy,iz,it) = 0.5  * x23
                    vout[3,1,ix,iy,iz,it) = 0.5  * x31
                    vout[3,2,ix,iy,iz,it) = 0.5  * x32
                    =#
                    y12 = 0.5  * x12
                    y13 = 0.5  * x13
                    y21 = 0.5  * x21
                    y23 = 0.5  * x23
                    y31 = 0.5  * x31
                    y32 =  0.5  * x32

                    c[1,ix,iy,iz,it] = ( imag(y12) + imag(y21) )*factor + c[1,ix,iy,iz,it]
                    c[2,ix,iy,iz,it] = ( real(y12) - real(y21) )*factor + c[2,ix,iy,iz,it]
                    c[3,ix,iy,iz,it] = ( imag(y11) - imag(y22) )*factor + c[3,ix,iy,iz,it]
                    c[4,ix,iy,iz,it] = ( imag(y13) + imag(y31) )*factor  + c[4,ix,iy,iz,it]
                    c[5,ix,iy,iz,it] = ( real(y13) - real(y31) )*factor  + c[5,ix,iy,iz,it]

                    c[6,ix,iy,iz,it] = ( imag(y23) + imag(y32) )*factor  + c[6,ix,iy,iz,it]
                    c[7,ix,iy,iz,it] = ( real(y23) - real(y32) )*factor + c[7,ix,iy,iz,it]
                    c[8,ix,iy,iz,it] = sr3i *
                            ( imag(y11) + imag(y22) -
                                    2*imag(y33) )*factor  + c[8,ix,iy,iz,it]

                end
            end
        end
    end
    barrier(c)

end

function Traceless_antihermitian_add!(c::TA_Gaugefields_4D_mpi{2,NumofBasis},factor,vin::Union{Gaugefields_4D_wing_mpi{2},Gaugefields_4D_nowing_mpi{2}}) where NumofBasis
    #error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
    fac12 = 1/2
    NX = vin.NX
    NY = vin.NY
    NZ = vin.NZ
    NT = vin.NT

    for it=1:vin.PN[4]
        for iz=1:vin.PN[3]
            for iy=1:vin.PN[2]
                @simd for ix=1:vin.PN[1]
                    v11 = getvalue(vin,1,1,ix,iy,iz,it)
                    v22 = getvalue(vin,2,2,ix,iy,iz,it)

                    tri = fac12*(imag(v11)+imag(v22))

                    v12 = getvalue(vin,1,2,ix,iy,iz,it)
                    #v13 = getvalue(vin,1,3,ix,iy,iz,it)
                    v21 = getvalue(vin,2,1,ix,iy,iz,it)

                    x12 = v12 - conj(v21)

                    x21 = - conj(x12)

                    y11 = (imag(v11)-tri)*im
                    y12 = 0.5  * x12
                    y21 = 0.5  * x21
                    y22 = (imag(v22)-tri)*im

                    c[1,ix,iy,iz,it] = (imag(y12)+imag(y21))*factor  + c[1,ix,iy,iz,it]
                    c[2,ix,iy,iz,it] = (real(y12)-real(y21))*factor  + c[2,ix,iy,iz,it]
                    c[3,ix,iy,iz,it] = (imag(y11)-imag(y22))*factor  + c[3,ix,iy,iz,it]

                end
            end
        end
    end
    barrier(c)

end

"""
-----------------------------------------------------c
     !!!!!   vin and vout should be different vectors

     Projectin of the etraceless antiermite part
     vout = x/2 - Tr(x)/6
     wher   x = vin - Conjg(vin)
-----------------------------------------------------c
    """
function Traceless_antihermitian!(c::TA_Gaugefields_4D_mpi{3,NumofBasis},vin::Union{Gaugefields_4D_wing_mpi{3},Gaugefields_4D_nowing_mpi{3}}) where NumofBasis
    #error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
    fac13 = 1/3
    NX = vin.NX
    NY = vin.NY
    NZ = vin.NZ
    NT = vin.NT

    for it=1:vin.PN[4]
        for iz=1:vin.PN[3]
            for iy=1:vin.PN[2]
                @simd for ix=1:vin.PN[1]
                    v11 = getvalue(vin,1,1,ix,iy,iz,it)
                    v22 = getvalue(vin,2,2,ix,iy,iz,it)
                    v33 = getvalue(vin,3,3,ix,iy,iz,it)

                    tri = fac13*(imag(v11)+imag(v22)+imag(v33))

                    #=
                    vout[1,1,ix,iy,iz,it) = (imag(v11)-tri)*im
                    vout[2,2,ix,iy,iz,it) = (imag(v22)-tri)*im
                    vout[3,3,ix,iy,iz,it) = (imag(v33)-tri)*im
                    =#
                    y11 = (imag(v11)-tri)*im
                    y22 = (imag(v22)-tri)*im
                    y33 = (imag(v33)-tri)*im

                    v12 = getvalue(vin,1,2,ix,iy,iz,it)
                    v13 = getvalue(vin,1,3,ix,iy,iz,it)
                    v21 = getvalue(vin,2,1,ix,iy,iz,it)
                    v23 = getvalue(vin,2,3,ix,iy,iz,it)
                    v31 = getvalue(vin,3,1,ix,iy,iz,it)
                    v32 = getvalue(vin,3,2,ix,iy,iz,it)

                    x12 = v12 - conj(v21)
                    x13 = v13 - conj(v31)
                    x23 = v23 - conj(v32)

                    x21 = - conj(x12)
                    x31 = - conj(x13)
                    x32 = - conj(x23)

                    #=
                    vout[1,2,ix,iy,iz,it] = 0.5  * x12
                    vout[1,3,ix,iy,iz,it] = 0.5  * x13
                    vout[2,1,ix,iy,iz,it] = 0.5  * x21
                    vout[2,3,ix,iy,iz,it] = 0.5  * x23
                    vout[3,1,ix,iy,iz,it] = 0.5  * x31
                    vout[3,2,ix,iy,iz,it] = 0.5  * x32
                    =#
                    y12 = 0.5  * x12
                    y13 = 0.5  * x13
                    y21 = 0.5  * x21
                    y23 = 0.5  * x23
                    y31 = 0.5  * x31
                    y32 =  0.5  * x32

                    c[1,ix,iy,iz,it] = ( imag(y12) + imag(y21) )
                    c[2,ix,iy,iz,it] = ( real(y12) - real(y21) )
                    c[3,ix,iy,iz,it] = ( imag(y11) - imag(y22) )
                    c[4,ix,iy,iz,it] = ( imag(y13) + imag(y31) )
                    c[5,ix,iy,iz,it] = ( real(y13) - real(y31) )

                    c[6,ix,iy,iz,it] = ( imag(y23) + imag(y32) )
                    c[7,ix,iy,iz,it] = ( real(y23) - real(y32) )
                    c[8,ix,iy,iz,it] = sr3i *
                            ( imag(y11) + imag(y22) -
                                    2*imag(y33) )
                end
            end
        end
    end
    barrier(c)

end

function Traceless_antihermitian!(c::TA_Gaugefields_4D_mpi{NC,NumofBasis},vin::Union{Gaugefields_4D_wing_mpi{NC},Gaugefields_4D_nowing_mpi{NC}}) where {NC,NumofBasis}
    @assert NC != 3 && NC != 2
    #NC = vout.NC
    fac1N = 1/NC
    nv = vin.NV

    NX = vin.NX
    NY = vin.NY
    NZ = vin.NZ
    NT = vin.NT
    g = c.generators
    matrix = zeros(ComplexF64,NC,NC)
    a = zeros(ComplexF64,length(g))

    for it=1:vin.PN[4]
        for iz=1:vin.PN[3]
            for iy=1:vin.PN[2]
                @simd for ix=1:vin.PN[1]
                    tri = 0.0
                    @simd for k=1:NC
                        tri += imag(getvalue(vin,k,k,ix,iy,iz,it))
                    end
                    tri *= fac1N
                    @simd for k=1:NC
                        #vout[k,k,ix,iy,iz,it] = (imag(getvalue(vin,k,k,ix,iy,iz,it))-tri)*im
                        matrix[k,k] = (imag(getvalue(vin,k,k,ix,iy,iz,it))-tri)*im
                    end

                    @simd for k2=k1+1:NC
                        vv = 0.5*(getvalue(vin,k1,k2,ix,iy,iz,it) - conj(getvalue(vin,k2,k1,ix,iy,iz,it)))
                        #vout[k1,k2,ix,iy,iz,it] = vv
                        #vout[k2,k1,ix,iy,iz,it] = -conj(vv)
                        matrix[k1,k2] = vv
                        matrix[k2,k1] = -conj(vv)
                    end

                    matrix2lie!(a,g,matrix)
                    for k = 1:length(g)
                        c[k,ix,iy,iz,it] = 2*imag(a[k])
                    end

                end
            end
        end
    end
    barrier(c)

end

function Traceless_antihermitian_add!(c::TA_Gaugefields_4D_mpi{NC,NumofBasis},factor,vin::Union{Gaugefields_4D_wing_mpi{NC},Gaugefields_4D_nowing_mpi{NC}}) where {NC,NumofBasis}
    @assert NC != 3 && NC != 2 "NC should be NC >4! in this function"
    #NC = vout.NC
    fac1N = 1/NC
    nv = vin.NV

    NX = vin.NX
    NY = vin.NY
    NZ = vin.NZ
    NT = vin.NT
    g = c.generators
    matrix = zeros(ComplexF64,NC,NC)
    a = zeros(ComplexF64,length(g))

    for it=1:vin.PN[4]
        for iz=1:vin.PN[3]
            for iy=1:vin.PN[2]
                @simd for ix=1:vin.PN[1]
                    tri = 0.0
                    @simd for k=1:NC
                        tri += imag(getvalue(vin,k,k,ix,iy,iz,it))
                    end
                    tri *= fac1N
                    @simd for k=1:NC
                        #vout[k,k,ix,iy,iz,it] = (imag(getvalue(vin,k,k,ix,iy,iz,it))-tri)*im
                        matrix[k,k] = (imag(getvalue(vin,k,k,ix,iy,iz,it))-tri)*im
                    end

                    for k1=1:NC
                        @simd for k2=k1+1:NC
                            vv = 0.5*(getvalue(vin,k1,k2,ix,iy,iz,it) - conj(getvalue(vin,k2,k1,ix,iy,iz,it)))
                            #vout[k1,k2,ix,iy,iz,it] = vv
                            #vout[k2,k1,ix,iy,iz,it] = -conj(vv)
                            matrix[k1,k2] = vv
                            matrix[k2,k1] = -conj(vv)
                        end
                    end

                    matrix2lie!(a,g,matrix)
                    for k = 1:length(g)
                        c[k,ix,iy,iz,it] = 2*imag(a[k])*factor + c[k,ix,iy,iz,it]
                    end

                end
            end
        end
    end
    barrier(c)

end

function exptU!(uout::T,t::N,u::TA_Gaugefields_4D_mpi{NC,NumofBasis},temps::Array{T,1}) where {N <: Number, T <: Gaugefields_4D, NC,NumofBasis} #uout = exp(t*u)
    @assert NC != 3 && NC != 2 "This function is for NC != 2,3, now, NC = $NC, and NumofBasis = $NumofBasis"
    g = u.generators
    NT = u.NT
    NZ = u.NZ
    NY = u.NY
    NX = u.NX
    V = zeros(ComplexF64,NC,NC)

    u0 = zeros(ComplexF64,NC,NC)
    a = zeros(Float64,length(g))
    for it=1:vin.PN[4]
        for iz=1:vin.PN[3]
            for iy=1:vin.PN[2]
                @simd for ix=1:vin.PN[1]
                    for k=1:length(a)
                        a[k] = u[k,ix,iy,iz,it]
                    end

                    lie2matrix!(u0,g,a)
                    V[:,:] = exp(t*(im/2)*u0)
                    for k2=1:NC
                        for k1=1:NC
                            v = V[k1,k2]
                            setvalue!(uout,v,k1,k2,ix,iy,iz,it)
                        end
                    end

                    #uout[:,:,ix,iy,iz,it] = exp(t*(im/2)*u0)

                end
            end
        end
    end
    barrier(uout)
    #error("exptU! is not implemented in type $(typeof(u)) ")
end

function exptU!(uout::T,t::N,u::TA_Gaugefields_4D_mpi{3,NumofBasis},temps::Array{T,1}) where {N <: Number, T <: Gaugefields_4D,NumofBasis} #uout = exp(t*u)
    ww = temps[1]
    w = temps[2]
    NT = u.NT
    NZ = u.NZ
    NY = u.NY
    NX = u.NX

    for it=1:u.PN[4]
        for iz=1:u.PN[3]
            for iy=1:u.PN[2]
                for ix=1:u.PN[1]
                    c1 = t*u[1,ix,iy,iz,it] * 0.5
                    c2 = t*u[2,ix,iy,iz,it] * 0.5
                    c3 = t*u[3,ix,iy,iz,it] * 0.5
                    c4 = t*u[4,ix,iy,iz,it] * 0.5
                    c5 = t*u[5,ix,iy,iz,it] * 0.5
                    c6 = t*u[6,ix,iy,iz,it] * 0.5
                    c7 = t*u[7,ix,iy,iz,it] * 0.5
                    c8 = t*u[8,ix,iy,iz,it] * 0.5
                    csum = c1+c2+c3+c4+c5+c6+c7+c8
                    if csum == 0
                        v = 1
                        setvalue!(w,v,1,1,ix,iy,iz,it)# =   1
                        v = 0
                        setvalue!(w,v,1,2,ix,iy,iz,it)
                        setvalue!(w,v,1,3,ix,iy,iz,it)# =    0
                        setvalue!(w,v,2,1,ix,iy,iz,it)# =  0
                        v = 1
                        setvalue!(w,v,2,2,ix,iy,iz,it)# =   1
                        v = 0
                        setvalue!(w,v,2,3,ix,iy,iz,it)# =   0
                        setvalue!(w,v,3,1,ix,iy,iz,it)# =   0
                        setvalue!(w,v,3,2,ix,iy,iz,it)# =   0
                        v = 1
                        setvalue!(w,v,3,3,ix,iy,iz,it)# =   1

                        v = 1
                        setvalue!(ww,v,1,1,ix,iy,iz,it)# =   1
                        v = 0
                        setvalue!(ww,v,1,2,ix,iy,iz,it)# =   0
                        setvalue!(ww,v,1,3,ix,iy,iz,it)# =   0
                        setvalue!(ww,v,2,1,ix,iy,iz,it)# =   0
                        v = 1
                        setvalue!(ww,v,2,2,ix,iy,iz,it)# =   1
                        v = 0
                        setvalue!(ww,v,2,3,ix,iy,iz,it)# =   0
                        setvalue!(ww,v,3,1,ix,iy,iz,it)# =   0
                        setvalue!(ww,v,3,2,ix,iy,iz,it)# =   0
                        v = 1
                        setvalue!(ww,v,3,3,ix,iy,iz,it)# =   1
                        continue
                    end

                    #x[1,1,icum] =  c3+sr3i*c8 +im*(  0.0 )
                    v1 = c3+sr3i*c8
                    v2 = 0.0
                    #x[1,2,icum] =  c1         +im*( -c2   )
                    v3 = c1
                    v4 = -c2
                    #x[1,3,icum] =  c4         +im*(-c5   )
                    v5 = c4
                    v6 = -c5

                    #x[2,1,icum] =  c1         +im*(  c2   )
                    v7 = c1
                    v8 = c2

                    #x[2,2,icum] =  -c3+sr3i*c8+im*(  0.0 )
                    v9 =-c3+sr3i*c8
                    v10 = 0.0

                    #x[2,3,icum] =  c6         +im*( -c7   )
                    v11 = c6
                    v12 = -c7

                    #x[3,1,icum] =  c4         +im*(  c5   )
                    v13 = c4
                    v14 = c5

                    #x[3,2,icum] =  c6         +im*(  c7   )
                    v15 = c6
                    v16 = c7
                    #x[3,3,icum] =  -sr3i2*c8  +im*(  0.0 )
                    v17 = -sr3i2*c8
                    v18 = 0.0

                    #c find eigenvalues of v
                    trv3 = (v1 + v9 + v17) / 3.0
                    cofac = v1 *  v9  - v3^2  - v4^2 +
                            v1 * v17 -  v5^2  - v6^2 +
                            v9 * v17 - v11^2 - v12^2
                    det = v1 * v9 * v17 -
                            v1 * (v11^2 + v12^2) -
                            v9 * (v5^2 + v6^2) -
                            v17 * (v3^2 + v4^2) +
                            (v5 * (v3 * v11 - v4 * v12) +
                            v6 * (v3 * v12 + v4 * v11)) * 2.0

                    p3 = cofac / 3.0 - trv3^2
                    q = trv3 * cofac - det - 2.0 * trv3^3
                    x = sqrt(-4.0 * p3) + tinyvalue

                    arg = q / (x * p3)

                    arg = min(1, max(-1, arg))
                    theta = acos(arg) / 3.0
                    e1 = x * cos(theta) + trv3
                    theta = theta + pi23
                    e2 = x * cos(theta) + trv3
            #       theta = theta + pi23
            #       e3 = x * cos(theta) + trv3
                    e3 = 3.0 * trv3 - e1 - e2

            # solve for eigenvectors

                    w1 =   v5 * (v9 - e1) - v3 * v11 + v4 * v12
                    w2 = - v6 * (v9 - e1) + v4 * v11 + v3 * v12
                    w3 =   (v1 - e1) * v11 - v3 * v5 - v4 * v6
                    w4 = - (v1 - e1) * v12 - v4 * v5 + v3 * v6
                    w5 = - (v1 - e1) * (v9 - e1) +  v3^2 + v4^2
                    w6 = 0.0
                    #println("1c $w1 $w2 $w3 $w4 $w5 $w6 ",)
                    #coeffv = sqrt(w1^2 + w2^2 + w3^2 + w4^2 + w5^2)

                    #coeff = ifelse(coeffv == zero(coeffv),0,coeffv)
                    coeff = 1.0 / sqrt(w1^2 + w2^2 + w3^2 + w4^2 + w5^2)
                    #println("1 ",coeff)

                    w1 = w1 * coeff
                    w2 = w2 * coeff
                    w3 = w3 * coeff
                    w4 = w4 * coeff
                    w5 = w5 * coeff

                    w7 =   v5 * (v9 - e2) - v3 * v11 + v4 * v12
                    w8 = - v6 * (v9 - e2) + v4 * v11 + v3 * v12
                    w9 = (v1 - e2) * v11 - v3 * v5 - v4 * v6
                    w10 = - (v1 - e2) * v12 - v4 * v5 + v3 * v6
                    w11 = - (v1 - e2) * (v9 - e2) +  v3^2 + v4^2
                    w12 = 0.0

                    coeff = 1.0 / sqrt(w7^2  + w8^2 + w9^2+ w10^2 + w11^2)

                    w7  = w7  * coeff
                    w8  = w8  * coeff
                    w9  = w9  * coeff
                    w10 = w10 * coeff
                    w11 = w11 * coeff

                    w13 =   v5 * (v9 - e3) - v3 * v11 + v4 * v12
                    w14 = - v6 * (v9 - e3) + v4 * v11 + v3 * v12
                    w15 =   (v1 - e3) * v11 - v3 * v5 - v4 * v6
                    w16 = - (v1 - e3) * v12 - v4 * v5 + v3 * v6
                    w17 = - (v1 - e3) * (v9 - e3) +  v3^2 + v4^2
                    w18 = 0.0

                    coeff = 1.0 / sqrt(w13^2 + w14^2 + w15^2+ w16^2 + w17^2)
                    w13 = w13 * coeff
                    w14 = w14 * coeff
                    w15 = w15 * coeff
                    w16 = w16 * coeff
                    w17 = w17 * coeff

            # construct the projection v
                    c1 = cos(e1)
                    s1 = sin(e1)
                    ww1  = w1  * c1 - w2 * s1
                    ww2  = w2  * c1 + w1 * s1
                    ww3  = w3  * c1 - w4 * s1
                    ww4  = w4  * c1 + w3 * s1
                    ww5  = w5  * c1 - w6 * s1
                    ww6  = w6  * c1 + w5 * s1

                    c2 = cos(e2)
                    s2 = sin(e2)
                    ww7  = w7  * c2 - w8 * s2
                    ww8  = w8  * c2 + w7 * s2
                    ww9  = w9  * c2 - w10 * s2
                    ww10 = w10 * c2 + w9 * s2
                    ww11 = w11 * c2 - w12 * s2
                    ww12 = w12 * c2 + w11 * s2

                    c3 = cos(e3)
                    s3 = sin(e3)
                    ww13 = w13 * c3 - w14 * s3
                    ww14 = w14 * c3 + w13 * s3
                    ww15 = w15 * c3 - w16 * s3
                    ww16 = w16 * c3 + w15 * s3
                    ww17 = w17 * c3 - w18 * s3
                    ww18 = w18 * c3 + w17 * s3

                    v = w1+im*w2
                    setvalue!(w,v,1,1,ix,iy,iz,it)# =   1
                    v = w3+im*w4
                    setvalue!(w,v,1,2,ix,iy,iz,it)
                    v = w5+im*w6
                    setvalue!(w,v,1,3,ix,iy,iz,it)# =    0
                    v = w7+im*w8
                    setvalue!(w,v,2,1,ix,iy,iz,it)# =  0
                    v = w9+im*w10

                    setvalue!(w,v,2,2,ix,iy,iz,it)# =   1
                    v = w11+im*w12
                    setvalue!(w,v,2,3,ix,iy,iz,it)# =   0
                    v = w13+im*w14
                    setvalue!(w,v,3,1,ix,iy,iz,it)# =   0
                    v = w15+im*w16
                    setvalue!(w,v,3,2,ix,iy,iz,it)# =   0
                    v = 1\w17+im*w18
                    setvalue!(w,v,3,3,ix,iy,iz,it)# =   1

                    v = ww1+im*ww2
                    setvalue!(ww,v,1,1,ix,iy,iz,it)# =   1
                    v = ww3+im*ww4
                    setvalue!(ww,v,1,2,ix,iy,iz,it)# =   0
                    v = ww5+im*ww6
                    setvalue!(ww,v,1,3,ix,iy,iz,it)# =   0
                    v = ww7+im*ww8
                    setvalue!(ww,v,2,1,ix,iy,iz,it)# =   0
                    v = ww9+im*ww10
                    setvalue!(ww,v,2,2,ix,iy,iz,it)# =   1
                    v =  ww11+im*ww12
                    setvalue!(ww,v,2,3,ix,iy,iz,it)# =   0
                    v = ww13+im*ww14
                    setvalue!(ww,v,3,1,ix,iy,iz,it)# =   0
                    v = ww15+im*ww16
                    setvalue!(ww,v,3,2,ix,iy,iz,it)# =   0
                    v =  ww17+im*ww18
                    setvalue!(ww,v,3,3,ix,iy,iz,it)# =   1

                end
            end
        end
    end
    barrier(uout)

    mul!(uout,w',ww)
    barrier(uout)

end
const tinyvalue =1e-100

function exptU!(uout::T,t::N,u::TA_Gaugefields_4D_mpi{2,NumofBasis},temps::Array{T,1}) where {N <: Number, T <: Gaugefields_4D,NumofBasis} #uout = exp(t*u)
    NT = u.NT
    NZ = u.NZ
    NY = u.NY
    NX = u.NX

    for it=1:u.PN[4]
        for iz=1:u.PN[3]
            for iy=1:u.PN[2]
                @simd for ix=1:u.PN[1]
                    #icum = (((it-1)*NX+iz-1)*NY+iy-1)*NX+ix
                    u1 = t*u[1,ix,iy,iz,it]/2
                    u2 = t*u[2,ix,iy,iz,it]/2
                    u3 = t*u[3,ix,iy,iz,it]/2
                    R = sqrt(u1^2+u2^2+u3^2) +  tinyvalue
                    sR = sin(R)/R
                    #sR = ifelse(R == 0,1,sR)
                    a0 = cos(R)
                    a1 = u1*sR
                    a2 = u2*sR
                    a3 = u3*sR

                    v = cos(R) + im*a3
                    setvalue!(uout,v,1,1,ix,iy,iz,it)
                    v = im*a1 + a2
                    setvalue!(uout,v,1,2,ix,iy,iz,it)
                    #uout[2,1,ix,iy,iz,it) =
                    v = im*a1 - a2
                    setvalue!(uout,v,2,1,ix,iy,iz,it)
                    #uout[2,2,ix,iy,iz,it)=
                    v = cos(R) - im*a3
                    setvalue!(uout,v,2,2,ix,iy,iz,it)

                end
            end
        end
    end

    barrier(uout)

end


text:
using OrthogonalPolynomialsQuasi, ContinuumArrays, QuasiArrays, BandedMatrices,
        SemiseparableMatrices, LazyArrays, ArrayLayouts, Test

import QuasiArrays: MulQuasiMatrix
import ContinuumArrays: MappedBasisLayout, MappedWeightedBasisLayout
import LazyArrays: arguments
import SemiseparableMatrices: VcatAlmostBandedLayout

@testset "ODEs" begin
    @testset "p-FEM" begin
        S = Jacobi(true,true)
        w = JacobiWeight(true,true)
        D = Derivative(axes(S,1))
        P = Legendre()

        @test w.*S isa QuasiArrays.BroadcastQuasiMatrix

        M = P\(D*(w.*S))
        @test M isa BandedMatrix
        @test M[1:10,1:10] == diagm(-1 => -2.0:-2:-18.0)

        N = 10
        A = D* (w.*S)[:,1:N]
        @test A.args[1] == P
        @test P\(D*(w.*S)[:,1:N]) isa MulMatrix

        L = D*(w.*S)
        Δ = L'L
        @test Δ isa MulMatrix
        @test Δ[1:3,1:3] isa BandedMatrix
        @test bandwidths(Δ) == (0,0)

        L = D*(w.*S)[:,1:N]

        A  = *((L').args..., L.args...)
        @test A isa MulMatrix

        Δ = L'L
        @test Δ isa MulMatrix
        @test bandwidths(Δ) == (0,0)
        @test BandedMatrix(Δ) == Δ
        @test BandedMatrix(Δ) isa BandedMatrix
    end

    @testset "∞-FEM" begin
        S = Jacobi(true,true)
        w = JacobiWeight(true,true)
        D = Derivative(axes(w,1))
        WS = w.*S
        L = D* WS
        Δ = L'L
        P = Legendre()

        f = P * Vcat(randn(10), Zeros(∞))
        (P\WS)'*(P'P)*(P\WS)
        B = BroadcastArray(+, Δ, (P\WS)'*(P'P)*(P\WS))
        @test colsupport(B,1) == 1:3

        @test axes(B.args[2].args[1]) == (Base.OneTo(∞),Base.OneTo(∞))
        @test axes(B.args[2]) == (Base.OneTo(∞),Base.OneTo(∞))
        @test axes(B) == (Base.OneTo(∞),Base.OneTo(∞))

        @test BandedMatrix(view(B,1:10,13:20)) == zeros(10,8)

        F = qr(B);
        b = Vcat(randn(10), Zeros(∞))
        @test B*(F \ b) ≈ b
    end

    @testset "Collocation" begin
        P = Chebyshev()
        D = Derivative(axes(P,1))
        n = 300
        x = cos.((0:n-2) .* π ./ (n-2))
        cfs = [P[-1,1:n]'; (D*P)[x,1:n] - P[x,1:n]] \ [exp(-1); zeros(n-1)]
        u = P[:,1:n]*cfs
        @test u[0.1] ≈ exp(0.1)

        P = Chebyshev()
        D = Derivative(axes(P,1))
        D2 = D^2 * P # could be D^2*P in the future
        n = 300
        x = cos.((1:n-2) .* π ./ (n-1)) # interior Chebyshev points
        C = [P[-1,1:n]';
            D2[x,1:n] + P[x,1:n];
            P[1,1:n]']
        cfs = C \ [1; zeros(n-2); 2] # Chebyshev coefficients
        u = P[:,1:n]*cfs  # interpret in basis
        @test u[0.1] ≈ (3cos(0.1)sec(1) + csc(1)sin(0.1))/2
    end

    @testset "∞-dimensional Dirichlet" begin
        S = Jacobi(true,true)
        w = JacobiWeight(true,true)
        D = Derivative(axes(S,1))
        X = Diagonal(Inclusion(axes(S,1)))

        @test_broken (Legendre() \ S)*(S\(w.*S))
        @test (Ultraspherical(3/2)\(D^2*(w.*S)))[1:10,1:10] ≈ diagm(0 => -(2:2:20))
    end

    @testset "rescaled" begin
        x = Inclusion(0..1)
        S = Jacobi(1.0,1.0)[2x.-1,:]
        D = Derivative(x)
        f = S*[[1,2,3]; zeros(∞)]
        g = Jacobi(1.0,1.0)*[[1,2,3]; zeros(∞)]
        @test f[0.1] ≈ g[2*0.1-1]
        h = 0.0000001
        @test (D*f)[0.1] ≈ (f[0.1+h]-f[0.1])/h atol=100h
        @test Jacobi(2.0,2.0)[2x.-1,:] \ (D*S).args[1] isa BandedMatrix
        @test (Jacobi(2.0,2.0)[2x.-1,:] \ (D*S))[1:10,1:10] == diagm(1 => 4:12)

        P = Legendre()[2x.-1,:]
        w = JacobiWeight(1.0,1.0)
        wS = (w .* Jacobi(1.0,1.0))[2x.-1,:]
        @test MemoryLayout(wS) isa MappedWeightedBasisLayout
        f = wS*[[1,2,3]; zeros(∞)]
        g = (w .* Jacobi(1.0,1.0))*[[1,2,3]; zeros(∞)]
        @test f[0.1] ≈ g[2*0.1-1]
        h = 0.0000001
        @test (D*f)[0.1] ≈ (f[0.1+h]-f[0.1])/h atol=100h
        @test P == P

        @test P.parent == (D*wS).args[1].parent
        DwS = apply(*,D,wS)
        A,B = P,arguments(DwS)[1];
        @test (A.parent\B.parent) == Eye(∞)
        @test (A \ B)[1:10,1:10] == diagm(-1 => ones(9))
        @test (P \ DwS)[1:10,1:10] == diagm(-1 => -4:-4:-36)
    end

    @testset "Beam" begin
        P = JacobiWeight(0.0,0.0) .* Jacobi(0.0,0.0)
        x = axes(P,1)
        D = Derivative(x)
        @test (D*P).args[1] == Jacobi{Float64}(1,1)
        @test (Jacobi(1,1)\(D*P))[1:10,1:10] ≈ (Jacobi(1,1) \ (D*Legendre()))[1:10,1:10]

        S = JacobiWeight(2.0,2.0) .* Jacobi(2.0,2.0)
        @test (Legendre() \ S)[1,1] ≈ 0.533333333333333333
        Δ² = (D^2*S)'*(D^2*S)
        M = S'S
    end

    @testset "Ultraspherical spectral method" begin
        T = ChebyshevT()
        U = ChebyshevU()
        x = axes(T,1)
        D = Derivative(x)
        A = U\(D*T) - U\T
        @test copyto!(BandedMatrix{Float64}(undef, (10,10), (0,2)), view(A,1:10,1:10)) == A[1:10,1:10]
        L = Vcat(T[1:1,:], A)
        @test L[1:10,1:10] isa AlmostBandedMatrix
        @test MemoryLayout(L) isa VcatAlmostBandedLayout
        u = L \ [ℯ; zeros(∞)]
        @test T[0.1,:]'u ≈ (T*u)[0.1] ≈ exp(0.1)

        C = Ultraspherical(2)
        A = C \ (D^2 * T) - C\(x .* T)
        L = Vcat(T[[-1,1],:], A)
        @test qr(L).factors[1:10,1:10] ≈ qr(L[1:13,1:10]).factors[1:10,1:10]
        u = L \ [airyai(-1); airyai(1); Zeros(∞)]
        @test T[0.1,:]'u ≈ airyai(0.1)

        ε = 0.0001
        A = ε^2 * (C \ (D^2 * T)) - C\(x .* T)
        L = Vcat(T[[-1,1],:], A)
        u = L \ [airyai(-ε^(-2/3)); airyai(ε^(2/3)); zeros(∞)]
        @test T[-0.1,:]'u ≈ airyai(-0.1*ε^(-2/3))
    end
end


text:
using LinearElasticity

# Define problem
nels = (60,20,4)
sizes = (1.0,1.0,1.0)
E = 1.0;
ν = 0.3;
force = -1.0;
problem = PointLoadCantilever(nels, sizes, E, ν, force)

# Build element stiffness matrices and force vectors
einfo = ElementFEAInfo(problem);

# Assemble global stiffness matrix and force vector
ginfo = assemble(problem, einfo);

# Solve for node displacements
u = ginfo.K \ ginfo.f

mesh = problem.ch.dh.grid
node_dofs = problem.metadata.node_dofs

node_displacements = reshape(u[node_dofs], 3, JuAFEM.getnnodes(mesh))

using Makie
using GeometryTypes

# Overloads to get the fem nodes working!

function AbstractPlotting.to_vertices(cells::AbstractVector{<: JuAFEM.Node{N, T}}) where {N, T}
    convert_arguments(nothing, cells)[1]
end
function AbstractPlotting.to_gl_indices(cells::AbstractVector{<: JuAFEM.Cell})
    tris = GLTriangle[]
    for cell in cells
        to_triangle(tris, cell)
    end
    tris
end

function to_triangle(tris, cell::JuAFEM.Cell{3, 8, 6})
    nodes = cell.nodes
    push!(tris, GLTriangle(nodes[1], nodes[2], nodes[5]))
    push!(tris, GLTriangle(nodes[5], nodes[2], nodes[6]))

    push!(tris, GLTriangle(nodes[6], nodes[2], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[6], nodes[7]))

    push!(tris, GLTriangle(nodes[7], nodes[8], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[8], nodes[4]))

    push!(tris, GLTriangle(nodes[4], nodes[8], nodes[5]))
    push!(tris, GLTriangle(nodes[5], nodes[4], nodes[1]))

    push!(tris, GLTriangle(nodes[1], nodes[2], nodes[3]))
    push!(tris, GLTriangle(nodes[3], nodes[1], nodes[4]))
end

function AbstractPlotting.convert_arguments(P, x::AbstractVector{<: JuAFEM.Node{N, T}}) where {N, T}
    convert_arguments(P, reinterpret(Point{N, T}, x))
end
#TODO make this work without creating a Node
cnode = Node(zeros(Float32, length(mesh.nodes)))
scene = Makie.mesh(mesh.nodes, mesh.cells, color = cnode, colorrange = (0.0, 33.0), shading = false)
mplot = scene[end]
displacevec = reinterpret(Vec{3, Float64}, node_displacements, (size(node_displacements, 2),))
displace = norm.(displacevec)

new_nodes = broadcast(1:length(mesh.nodes), mesh.nodes) do i, node
    JuAFEM.Node(ntuple(Val(3)) do j
        node.x[j] + node_displacements[j, i]
    end)
end

mesh!(mesh.nodes, mesh.cells, color = (:gray, 0.4))

scatter!(Point3f0.(getfield.(new_nodes, :x)), markersize = 0.1)
# TODO make mplot[1] = new_nodes work

mplot.input_args[1][] = new_nodes
# TODO make mplot[:color] = displace work
push!(cnode, displace)

points = reinterpret(Point{3, Float64}, mesh.nodes)
arrows!(points, displacevec, linecolor = (:black, 0.3))
scene