text:
function remove_node!(integrator::DiffEqBase.DEIntegrator, I...)
    idxs = getindices(integrator.u, I...)
    for c in full_cache(integrator)
        remove_node!(c, I...)
    end
    if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
      for c in DiffEqBase.ratenoise_cache(integrator)
          remove_node!(c, I...)
      end
    end
    remove_node_non_user_cache!(integrator, idxs, I...) # required to do noise correctly
end

function add_node!(integrator::DiffEqBase.DEIntegrator, x, I...)
    cur_len = length(integrator.u)
    add_len = length(x)
    last_idx = length(integrator.u[I...].nodes)
    idx_start = getindices(integrator.u, last_idx)[end] + 1
    idxs = idx_start:idx_start+add_len-1
    for c in full_cache(integrator)
        add_node!(c, recursivecopy(x), I...)
    end
    if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
      for c in DiffEqBase.ratenoise_cache(integrator)
          add_node!(c, recursivecopy(x), I...)
      end
    end
    #addat_non_user_cache!(integrator, idxs)
    add_node_non_user_cache!(integrator, idxs, x, I...) # required to do noise correctly
end

function add_node!(integrator::DiffEqBase.DEIntegrator, x)
    cur_len = length(integrator.u)
    add_len = length(x)
    last_idx = length(integrator.u.nodes)
    idx_start = getindices(integrator.u, last_idx)[end] + 1
    idxs = idx_start:idx_start+add_len-1
    for c in full_cache(integrator)
        add_node!(c, recursivecopy(x))
    end
    if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
      for c in DiffEqBase.ratenoise_cache(integrator)
          add_node!(c, recursivecopy(x))
      end
    end
    add_node_non_user_cache!(integrator, idxs, fill!(similar(x, eltype(x)),0)) # required to do noise correctly
end

reshape(m::AbstractMultiScaleArray, i::Int...) = m

function remove_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,idxs,node...)
  remove_node_non_user_cache!(integrator,integrator.cache,node...)
end

remove_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
                         cache::OrdinaryDiffEq.OrdinaryDiffEqCache,idxs,node...) = nothing

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,idxs,x::AbstractArray)
  add_node_non_user_cache!(integrator,integrator.cache,x)
end

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,idxs,x::AbstractArray,node...)
  add_node_non_user_cache!(integrator,integrator.cache,x,node...)
end

add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
                        cache::OrdinaryDiffEq.OrdinaryDiffEqCache,
                        x::AbstractArray) = nothing
add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
                        cache::OrdinaryDiffEq.OrdinaryDiffEqCache,
                        x::AbstractArray,node...) = nothing

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
                      cache::OrdinaryDiffEq.RosenbrockMutableCache,x::AbstractArray)
  i = length(integrator.u)
  cache.J = similar(cache.J,i,i)
  cache.W = similar(cache.W,i,i)
  add_node_jac_config!(cache,cache.jac_config,i,x)
  add_node_grad_config!(cache,cache.grad_config,i,x)
  nothing
end

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
                      cache::OrdinaryDiffEq.RosenbrockMutableCache,x::AbstractArray,node...)
  i = length(integrator.u)
  cache.J = similar(cache.J,i,i)
  cache.W = similar(cache.W,i,i)
  add_node_jac_config!(cache,cache.jac_config,i,x,node...)
  add_node_grad_config!(cache,cache.grad_config,i,x,node...)
  nothing
end

function remove_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
                                     cache::OrdinaryDiffEq.RosenbrockMutableCache,
                                     node...)
  i = length(integrator.u)
  cache.J = similar(cache.J,i,i)
  cache.W = similar(cache.W,i,i)
  remove_node_jac_config!(cache,cache.jac_config,i,node...)
  remove_node_grad_config!(cache,cache.grad_config,i,node...)
  nothing
end

function add_node_jac_config!(cache,config::SparseDiffTools.ForwardColorJacCache,i,x,node...)
    @assert cache.jac_config.colorvec isa UnitRange
    cache.jac_config = SparseDiffTools.ForwardColorJacCache(cache.uf,cache.uprev,config.chunksize)
    nothing
end

function add_node_jac_config!(cache,config::SparseDiffTools.ForwardColorJacCache,i,x)
    @assert cache.jac_config.colorvec isa UnitRange
    cache.jac_config = SparseDiffTools.ForwardColorJacCache(cache.uf,cache.uprev,config.chunksize)
    nothing
end

function remove_node_jac_config!(cache,config::SparseDiffTools.ForwardColorJacCache,i,x,node...)
    @assert cache.jac_config.colorvec isa UnitRange
    cache.jac_config = SparseDiffTools.ForwardColorJacCache(cache.uf,cache.uprev,config.chunksize)
    nothing
end

function remove_node_jac_config!(cache,config::SparseDiffTools.ForwardColorJacCache,i,x)
    @assert cache.jac_config.colorvec isa UnitRange
    cache.jac_config = SparseDiffTools.ForwardColorJacCache(cache.uf,cache.uprev,config.chunksize)
    nothing
end

function add_node_jac_config!(cache,config::FiniteDiff.JacobianCache,i,x)
    #add_node!(cache.x1, fill!(similar(x, eltype(cache.x1)),0))
    add_node!(config.fx, recursivecopy(x))
    #cache.fx1 !== nothing && add_node!(cache.fx1, fill!(similar(x, eltype(cache.fx1)),0))
    config.colorvec = 1:i
    nothing
end

function add_node_jac_config!(cache,config::FiniteDiff.JacobianCache,i,x,I...)
    #add_node!(cache.x1, fill!(similar(x, eltype(cache.x1)),0), I...)
    add_node!(config.fx, recursivecopy(x), I...)
    #cache.fx1 !== nothing && add_node!(cache.fx1, fill!(similar(x, eltype(cache.fx1)),0), I...)
    config.colorvec = 1:i
    nothing
end

function remove_node_jac_config!(cache,config::FiniteDiff.JacobianCache,i,I...)
    #add_node!(cache.x1, fill!(similar(x, eltype(cache.x1)),0), I...)
    remove_node!(config.fx, I...)
    #cache.fx1 !== nothing && add_node!(cache.fx1, fill!(similar(x, eltype(cache.fx1)),0), I...)
    config.colorvec = 1:i
    nothing
end

function add_node_grad_config!(cache,grad_config::ForwardDiff.DerivativeConfig,i,x)
    cache.grad_config = ForwardDiff.DerivativeConfig(cache.tf,cache.du1,cache.uf.t)
    nothing
end

function add_node_grad_config!(cache,grad_config::ForwardDiff.DerivativeConfig,i,x,I...)
  cache.grad_config = ForwardDiff.DerivativeConfig(cache.tf,cache.du1,cache.uf.t)
  nothing
end

function remove_node_grad_config!(cache,grad_config::ForwardDiff.DerivativeConfig,i,x)
    cache.grad_config = ForwardDiff.DerivativeConfig(cache.tf,cache.du1,cache.uf.t)
    nothing
end

function remove_node_grad_config!(cache,grad_config::ForwardDiff.DerivativeConfig,i,x,I...)
  cache.grad_config = ForwardDiff.DerivativeConfig(cache.tf,cache.du1,cache.uf.t)
  nothing
end

function add_node_grad_config!(cache,grad_config::AbstractArray,i,x)
    cache.grad_config = ForwardDiff.Dual{typeof(ForwardDiff.Tag(cache.tf,eltype(cache.du1)))}.(cache.du1, cache.du1)
    nothing
end

function add_node_grad_config!(cache,grad_config::AbstractArray,i,x,I...)
  cache.grad_config = ForwardDiff.Dual{typeof(ForwardDiff.Tag(cache.tf,eltype(cache.du1)))}.(cache.du1, cache.du1)
  nothing
end

function remove_node_grad_config!(cache,grad_config::AbstractArray,i,x)
    cache.grad_config = ForwardDiff.Dual{typeof(ForwardDiff.Tag(cache.tf,eltype(cache.du1)))}.(cache.du1, cache.du1)
    nothing
end

function remove_node_grad_config!(cache,grad_config::AbstractArray,i,x,I...)
  cache.grad_config = ForwardDiff.Dual{typeof(ForwardDiff.Tag(cache.tf,eltype(cache.du1)))}.(cache.du1, cache.du1)
  nothing
end

function add_node_grad_config!(cache,grad_config::FiniteDiff.GradientCache,i,x,I...)
  grad_config.fx !== nothing && add_node!(grad_config.fx, recursivecopy(x), I...)
  grad_config.c1 !== nothing && add_node!(grad_config.c1, recursivecopy(x), I...)
  grad_config.c2 !== nothing && add_node!(grad_config.c2, recursivecopy(x), I...)
  grad_config
end

function add_node_grad_config!(cache,grad_config::FiniteDiff.GradientCache,i,x)
  grad_config.fx !== nothing && add_node!(grad_config.fx, recursivecopy(x))
  grad_config.c1 !== nothing && add_node!(grad_config.c1, recursivecopy(x))
  grad_config.c2 !== nothing && add_node!(grad_config.c2, recursivecopy(x))
  grad_config
end

function remove_node_grad_config!(cache,grad_config::FiniteDiff.GradientCache,i,I...)
  grad_config.fx !== nothing && remove_node!(grad_config.fx, I...)
  grad_config.c1 !== nothing && remove_node!(grad_config.c1, I...)
  grad_config.c2 !== nothing && remove_node!(grad_config.c2, I...)
  grad_config
end

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractSDEIntegrator,idxs,x,node...)
  #addat_non_user_cache!(integrator, idxs)
  if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
    add_node_noise!(integrator,idxs,x,node...)
    for c in rand_cache(integrator)
      add_node!(c,copy(x),node...)
    end
  end
end

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractSDEIntegrator,idxs,x)
  #addat_non_user_cache!(integrator, idxs)
  if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
    add_node_noise!(integrator,idxs,x)
    for c in rand_cache(integrator)
      add_node!(c,copy(x))
    end
  end
end

function remove_node_non_user_cache!(integrator::DiffEqBase.AbstractSDEIntegrator,idxs,node...)
  #deleteat_non_user_cache!(integrator, idxs)
  if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
    remove_node_noise!(integrator,node...)
    for c in rand_cache(integrator)
      remove_node!(c,node...)
    end
  end
end

### For if noise is an AMSA

function add_node_noise!(integrator,idxs,x,node...)
  for c in integrator.W.S₁
    add_node!(c[2],copy(x),node...)
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x),node...)
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end
  for c in integrator.W.S₂
    add_node!(c[2],copy(x),node...)
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x),node...)
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end

  add_node!(integrator.W.dW,copy(x),node...)
  integrator.W.dW[idxs] .= zero(eltype(integrator.u))
  add_node!(integrator.W.curW,copy(x),node...)
  integrator.W.curW[idxs] .= zero(eltype(integrator.u))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZ,copy(x),node...)
    integrator.W.dZ[idxs] .= zero(eltype(integrator.u))
    add_node!(integrator.W.curZ,copy(x),node...)
    integrator.W.curZ[idxs] .= zero(eltype(integrator.u))
  end

  i = length(integrator.u)
  add_node!(integrator.W.dWtilde,copy(x),node...)
  add_node!(integrator.W.dWtmp,copy(x),node...)
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZtmp,copy(x),node...)
    add_node!(integrator.W.dZtilde,copy(x),node...)
  end

  # fill in rands
  fill!(@view(integrator.W.curW[idxs]),zero(eltype(integrator.u)))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    fill!(@view(integrator.W.curZ[idxs]),zero(eltype(integrator.u)))
  end
end

function add_node_noise!(integrator,idxs,x)
  for c in integrator.W.S₁
    add_node!(c[2],copy(x))
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x))
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end
  for c in integrator.W.S₂
    add_node!(c[2],copy(x))
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x))
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end

  add_node!(integrator.W.dW,copy(x))
  integrator.W.dW[idxs] .= zero(eltype(integrator.u))
  add_node!(integrator.W.curW,copy(x))
  integrator.W.curW[idxs] .= zero(eltype(integrator.u))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZ,copy(x))
    integrator.W.dZ[idxs] .= zero(eltype(integrator.u))
    add_node!(integrator.W.curZ,copy(x))
    integrator.W.curZ[idxs] .= zero(eltype(integrator.u))
  end

  i = length(integrator.u)
  add_node!(integrator.W.dWtilde,copy(x))
  add_node!(integrator.W.dWtmp,copy(x))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZtmp,copy(x))
    add_node!(integrator.W.dZtilde,copy(x))
  end

  # fill in rands
  fill!(@view(integrator.W.curW[idxs]),zero(eltype(integrator.u)))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    fill!(@view(integrator.W.curZ[idxs]),zero(eltype(integrator.u)))
  end
end

function remove_node_noise!(integrator,node...)
  for c in integrator.W.S₁
    remove_node!(c[2],node...)
    if DiffEqBase.alg_needs_extra_process(integrator.alg)
      remove_node!(c[3],node...)
    end
  end
  for c in integrator.W.S₂
    remove_node!(c[2],node...)
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      remove_node!(c[3],node...)
    end
  end
  remove_node!(integrator.W.dW,node...)
  remove_node!(integrator.W.dWtilde,node...)
  remove_node!(integrator.W.dWtmp,node...)
  remove_node!(integrator.W.curW,node...)

  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    remove_node!(integrator.W.curZ,node...)
    remove_node!(integrator.W.dZtmp,node...)
    remove_node!(integrator.W.dZtilde,node...)
    remove_node!(integrator.W.dZ,node...)
  end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    fitdims(dims, D)

Fit tuple `dims` to a given length `D`
by repeating the last dimension.
"""
function fitdims(dims::Dims{N}, D) where {N}
  ntuple(i -> i ≤ N ? dims[i] : last(dims), D)
end

"""
    signarea(A, B, C)

Compute signed area of triangle formed
by points `A`, `B` and `C`.
"""
function signarea(A::Point{2}, B::Point{2}, C::Point{2})
  ((B - A) × (C - A)) / 2
end

"""
    iscollinear(A, B, C)

Tells whether or not the points
`A`, `B` and `C` are collinear.
"""
function iscollinear(A::Point{Dim,T}, B::Point{Dim,T}, C::Point{Dim,T}) where {Dim,T}
  # points A, B, C are collinear if and only if the
  # cross-products for segments AB and AC with respect
  # to all possible pairs of coordinates are zero
  AB, AC = B - A, C - A
  result = true
  for i in 1:Dim, j in (i+1):Dim
    u = Vec{2,T}(AB[i], AB[j])
    v = Vec{2,T}(AC[i], AC[j])
    if !isapprox(u × v, zero(T), atol=atol(T)^2)
      result = false
      break
    end
  end
  result
end

"""
    sideof(point, segment)

Determines on which side of the oriented `segment`
the `point` lies. Possible results are `:LEFT`,
`:RIGHT` or `:ON` the segment.
"""
function sideof(p::Point{2,T}, s::Segment{2,T}) where {T}
  a, b = vertices(s)
  area = signarea(p, a, b)
  ifelse(area > atol(T), :LEFT, ifelse(area < -atol(T), :RIGHT, :ON))
end

"""
    sideof(point, chain)

Determines on which side of the closed `chain` the
`point` lies. Possible results are `:INSIDE` or
`:OUTSIDE` the chain.
"""
function sideof(p::Point{2,T}, c::Chain{2,T}) where {T}
  w = windingnumber(p, c)
  ifelse(isapprox(w, zero(T), atol=atol(T)), :OUTSIDE, :INSIDE)
end

"""
    proj2D(points)

Convert a vector of 3D `points` into a vector of 2D
points living in a plane of maximum variance using SVD.
"""
function proj2D(points::AbstractVector{Point{3,T}}) where {T}
  # https://math.stackexchange.com/a/99317
  X = mapreduce(coordinates, hcat, points)
  μ = sum(X, dims=2) / size(X, 2)
  Z = X .- μ
  U = svd(Z).U
  u = U[:,1]
  v = U[:,2]
  [Point(z⋅u, z⋅v) for z in eachcol(Z)]
end

"""
    dropunits(T)

Returns the unitless type of a (unitful) type or value. See Unitful.jl.
i.e. `dropunits(1u"mm") == Int`
"""
dropunits(v) = typeof(one(v))

ext: jl


text:
n = rand(1:10)
A = matrixdepot("smallworld", n)
B = matrixdepot("smallworld", Int, n)

C = matrixdepot("erdrey", n)
C2 = MatrixDepot.shortcuts(C, 0.6)
@test nnz(C) <= nnz(C2)
println("'erdrey' passed test...")


text:
export rotateAndMoveRBF,rotateAndMoveRBFsimple,wrapRBFparamAndRotationsTranslations,splitRBFparamAndRotationsTranslations, getIpIdxs

function rotateAndMoveRBF(m::Vector{Float64},Mesh::RegularMesh,theta_phi_rad::Array{Float64,2},trans::Array{Float64,2};numParamOfRBF = 5, computeJacobian = 1)
nRBF = div(length(m),numParamOfRBF);
nMoves = size(theta_phi_rad,1);
if nMoves != size(trans,1)
	error("rotateAndMoveRBF: trans vec and theta_phi vec should be of the same length");
end

m_ans = zeros(length(m),nMoves);

for k=1:nMoves
	m_ans[1:numParamOfRBF:end,k] = m[1:numParamOfRBF:end];
end

if numParamOfRBF==5
	for k=1:nMoves
		m_ans[2:numParamOfRBF:end,k] = m[2:numParamOfRBF:end];
	end
else
	P = zeros(9,6);
	P[idx_trilA(),1:6] = getIdentityMat(6);
	P[4,2] = 1.0; P[7,3] = 1.0; P[8,5] = 1.0;
end

mid_domain = [(Mesh.domain[1]+Mesh.domain[2])/2;(Mesh.domain[3]+Mesh.domain[4])/2;(Mesh.domain[5]+Mesh.domain[6])/2];
y = zeros(3);
zt = zeros(3);
zp = zeros(3);
o6 = ones(Int64,6);
o3 = ones(Int64,3);

h = 0.0001;

## The jacobian
## Each move has 5 parameters (theta, phi, b1, b2, b3)
## Dims: J: (5*nRBFs + 5*nMoves) -> 5*nRBFs*nMoves
## The rotation affects only the centers x.
## x_new = R*(x - mid) + b + mid

# x_new + dx_new = (R+R_tag_theta*dtheta+R_tag_phi*dphi)*(x+dx-mid) + b + db
# dx_new = R_tag_theta*(x - mid)*dtheta + R_tag_phi*(x - mid)*dphi + R*dx + db

## for each dip, and basis function, we have:
# dalpha_r		[1  0	 0   0	0 	0]	dalpha
#  dbeta_r	=	[0	1	 0   0	0 	0]	dbeta
#     dx_r 		[0  0	 R   zt	zp  I]	dx
									  # dtheta
									  # dphi
									  # db
# So, for each dip and basis function, we have 1+1+9+3+3+3 = 20 nz.

## for each dip, and basis function, we have:
# dalpha_r		[1  0	 0   0	0 	0]	dalpha
#  dAi_r	=	[0	iRiR 0   Tt	Tp 	0]	dAi_r
#  dx_r 		[0  0	 R   zt	zp  I]	dx
									  # dtheta
									  # dphi
									  # db
# So, for each dip and basis function, we have 1+1+9+3+3+3 = 20 nz.

nnzInJacPerMovePerRBF = numParamOfRBF==5 ? 20 : 67;

I = zeros(Int32,nnzInJacPerMovePerRBF*nRBF*nMoves);
J = zeros(Int32,nnzInJacPerMovePerRBF*nRBF*nMoves);
V = zeros(Float64,nnzInJacPerMovePerRBF*nRBF*nMoves);

g_idx = 1;

for k=1:nMoves
	b = vec(trans[k,:]);
	theta_rad = theta_phi_rad[k,1];
	phi_rad = theta_phi_rad[k,2];
	R = getRotate3D(theta_rad,phi_rad);
	invR = inv(R);
	if computeJacobian == 1
		R_tag_phi = (getRotate3D(theta_rad,phi_rad + h) - getRotate3D(theta_rad,phi_rad - h))./(2*h);
		R_tag_theta = (getRotate3D(theta_rad+h,phi_rad) - getRotate3D(theta_rad-h,phi_rad))./(2*h);
		if numParamOfRBF==10
			iRiR = kron(invR',invR');
			iRiR = (iRiR*P)[idx_trilA(),:];
			iR_tag_phi = -invR*R_tag_phi*invR;
			iR_tag_theta = -invR*R_tag_theta*invR;
		end
	end
	for j=1:nRBF
		y[:] = b;
		xj = getX(j,m,numParamOfRBF);
		xj .-= mid_domain;
		BLAS.gemv!('N',1.0,R,xj,1.0,y); # y = y + R*xj;
		y .+= mid_domain;
		m_ans[numParamOfRBF*(j-1).+((numParamOfRBF-2):numParamOfRBF),k] = y;
		if numParamOfRBF==10
			(a1,a2,a3,a4,a5,a6) = getA(j,m);
			A = [a1 a2 a3 ; a2 a4 a5 ; a3 a5 a6];
			Arot = invR'*A*invR;
			m_ans[10*(j-1) .+ (2:7),k] = Arot[idx_trilA()];
		end
		if computeJacobian == 1
			offset = ((k-1)*(numParamOfRBF*nRBF) + (j-1)*numParamOfRBF);
			offsetRBF = numParamOfRBF*(j-1);
			offsetMoves = nRBF*numParamOfRBF + 5*(k-1);
			temp = collect(offset .+ (1:numParamOfRBF));
			BLAS.gemv!('N',1.0,R_tag_theta,xj,0.0,zt);
			BLAS.gemv!('N',1.0,R_tag_phi,xj,0.0,zp);
			y[:] .= 1.0;
			if numParamOfRBF==5
				temp35 = temp[3:5];
				I[g_idx:(g_idx+nnzInJacPerMovePerRBF-1)] = [temp;temp35;temp35;temp35;temp35;temp35];
				J[g_idx:(g_idx+nnzInJacPerMovePerRBF-1)] = [offsetRBF.+[1;2;3;3;3;4;4;4;5;5;5]; offsetMoves .+ [1;1;1;2;2;2;3;4;5]];
				V[g_idx:(g_idx+nnzInJacPerMovePerRBF-1)] = [1.0;1.0;R[:];zt;zp;y];
			else
				Tt = invR'*A*iR_tag_theta; Tt = Tt + Tt'; Tt = Tt[idx_trilA()];
				Tp = invR'*A*iR_tag_phi;   Tp = Tp + Tp'; Tp = Tp[idx_trilA()];
				temp27 = temp[2:7];
				temp810 = temp[8:10];
				I[g_idx:(g_idx+nnzInJacPerMovePerRBF-1)] =
					[temp[1];temp27;temp27;temp27;temp27;temp27;temp27;temp810;temp810;temp810;temp27;temp810;temp27;temp810;temp810];
				J[g_idx:(g_idx+nnzInJacPerMovePerRBF-1)] =
					[offsetRBF.+[1;2*o6;3*o6;4*o6;5*o6;6*o6;7*o6; 8*o3;9*o3;10*o3]; offsetMoves.+[1*o6;1*o3;2*o6;2*o3;3;4;5]];
				y[:] .= 1.0;
				V[g_idx:(g_idx+nnzInJacPerMovePerRBF-1)] = [1.0;iRiR[:];R[:];Tt;zt;Tp;zp;y];
			end
			g_idx += nnzInJacPerMovePerRBF;
		end
	end
end
Jac = spzeros(0,0);
if computeJacobian == 1
	Jac = sparse(I,J,V,numParamOfRBF*nRBF*nMoves,(numParamOfRBF*nRBF + 5*nMoves));
end
return m_ans,Jac
end

function rotateAndMoveRBFsimple(m::Vector{Float64},Mesh::RegularMesh,theta_phi_rad::Array,trans::Array;computeJacobian = 1,numParamOfRBF::Int64 = 5)
(m_ans,Jac) = rotateAndMoveRBF(m,Mesh,theta_phi_rad,trans;numParamOfRBF = numParamOfRBF,computeJacobian = computeJacobian);
if computeJacobian==1
	Jac = Jac[:,1:length(m)];
end
return (m_ans,Jac);
end

function wrapRBFparamAndRotationsTranslations(m::Vector{Float64},theta_phi_rad::Array,trans::Array)
m_new = zeros(length(m) + length(theta_phi_rad) + length(trans));
m_new[1:length(m)] = m;
m_new[(length(m)+1):end] = vec([theta_phi_rad'; trans']);
return m_new;
end

function splitRBFparamAndRotationsTranslations(m_wraped::Vector{Float64},nRBF::Int64,nMoves::Int64,numParamOfRBF::Int64 = 5)
len_m = (numParamOfRBF*nRBF);
m_RBF = m_wraped[1:len_m];
numParamRotTrans = 5;
theta_phi = [m_wraped[(len_m+1):numParamRotTrans:end] m_wraped[(len_m+2):numParamRotTrans:end]];
translation =  [m_wraped[(len_m+3):numParamRotTrans:end] m_wraped[(len_m+4):numParamRotTrans:end] m_wraped[(len_m+5):numParamRotTrans:end]];
return (m_RBF,theta_phi,translation)
end

function getIpIdxs(I_p, nRBF::Int64,nMoves::Int64,numParamOfRBF::Int64 = 5)
Idxs = zeros(Bool,nRBF*numParamOfRBF + nMoves*5);
for k = I_p
	offset = nRBF*numParamOfRBF + (k-1)*5;
	Idxs[(offset+1):(offset+5)] .= true;
end
Idxs[1:nRBF*numParamOfRBF] .= true;
return Idxs;
end

# R = randn(3,3);
# dR = 0.01*randn(3,3);
# u0 = inv(R);
# for k=1:5
	# dR.*=0.5;
	# u = inv(R+dR);
	# ut = u0 - u0*dR*u0;
	# println("norm(ut-u0): ",vecnorm(u-u0),", norm(ut - u0 - J0*dm): ",vecnorm(u - ut));
# end

# function idx_trilA()
	# return [1;2;3;5;6;9];
# end
# # A = randn(3,3); A = A'*A + eye(3);Avec = A[idx_trilA()];
# dA = 0.01*randn(3,3); dA = dA + dA';
# R = randn(3,3);
# dAvec = dA[idx_trilA()];
# P = zeros(9,6);
# P[idx_trilA(),1:6] = eye(6);
# P[4,2] = 1.0; P[7,3] = 1.0; P[8,5] = 1.0;

# U = R'*dA*R; Uvec = U[idx_trilA()];
# RR = (kron(R',R')*P); RR = RR[idx_trilA(),:];
# Vvec = RR*dAvec;
# norm(Vvec - Uvec)

# for k=1:5


text:
using Test
import AMRVW
const A = AMRVW
using LinearAlgebra

# test polynomials
p4 = [24.0, -50.0, 35.0, -10.0, 1.0]
p5 = [-120.0, 274.0, -225.0, 85.0, -15.0, 1.0]
p6 = [720.0, -1764.0, 1624.0, -735.0, 175.0, -21.0, 1.0]
pc = [0.0 + 1.0im, -1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 - 1.0im, 1.0 + 0.0im]

T = Float64
S = Complex{T}

csort(xs) = sort(sort(xs, by=x -> imag(x)), by=x -> real(x))

@testset "transformations" begin

    Ts = [BigFloat, Float64, Float32, Float16]

    # Real case first
    ## givensrot
    for T in Ts
        a,b = rand(T, 2)
        c,s,r = A.givensrot(a,b)
        M =  diagm(0 => ones(T, 2))
        Rx = (A.Rotator(c,s,1)*M) * [a,b]
        @test abs(Rx[2]) <= sqrt(eps(T))
    end

    ## turnover
    for T in Ts
        M = diagm(0 => ones(T, 4))
        for i in 1:10
            U,V,W  = A.random_rotator.(T, [2,3,2])
            U1, V1, W1 = A.turnover(U,V,W)
            @test all(isapprox.([U1,V1,W1]*M , [U,V,W]*M, atol=4eps(T)))
        end
    end

    ## fuse
    for T in Ts
        U,V  = A.random_rotator.(T, [2,2])
        M = diagm(0 => ones(T, 4))
        UV, Di = A.fuse(U,V)
        @test all(isapprox.(UV*M, [U,V]*M, atol=2eps(T)))
    end

    # ComplexReal case
    ## givensrot
    for T in Ts
        a,b = rand(Complex{T}, 2)
        c,s,r = A.givensrot(a,b)
        M =  diagm(0 => ones(Complex{T}, 2))
        Rx = (A.Rotator(c,s,1)*M) * [a,b]
        @test abs(Rx[2]) <= sqrt(eps(T))
    end

    ## turnover
    ### a special case that had caused issues
    T = Float64
    S = Complex{T}
    U = A.Rotator(-0.9990813749617048 + 0.042846941864777804im,
                                   -0.0007387675316362891, 2)
    V = A.Rotator(0.20398474677465908 - 0.08295034712193844im,
                                         -0.9754534653153004, 1)
    W = A.Rotator(0.03541646337113308 + 0.28843087042807297im,
                                   -0.9568454980331911, 2)

    M = diagm(0 => ones(Complex{T}, 4))
    U1, V1, W1 = A.turnover(U,V,W)
    @test all(isapprox.([U1,V1,W1]*M , [U,V,W]*M))

    for T in Ts
        M = diagm(0 => ones(Complex{T}, 4))
        for i in 1:10
            U,V,W  = A.random_rotator.(T, [2,3,2])
            U1, V1, W1 = A.turnover(U,V,W)
            @test norm([U,V,W]*M - [U1,V1,W1]*M) <= sqrt(eps(T))
        end
    end

    ## when terms are 0
    for T in Ts
        S = Complex{T}
        M = diagm(0 => ones(S, 3))

        V,W = A.random_rotator.(S, [2,1])
        a,b =  sincos(T(pi/4))
        U = A.Rotator(complex(a,b), zero(T), 1)

        A.turnover(U,V,W)
        U1, V1, W1 = A.turnover(U,V,W)
        @test all(isapprox.([U1,V1,W1]*M , [U,V,W]*M, atol=4eps(T)))

        U,W  = A.random_rotator.(S, [1,1])
        a,b =  sincos(T(pi/4))
        V = A.Rotator(complex(a,b), zero(T), 2)

        A.turnover(U,V,W)
        U1, V1, W1 = A.turnover(U,V,W)
        @test all(isapprox.([U1,V1,W1]*M , [U,V,W]*M, atol=4eps(T)))

        U,V = A.random_rotator.(S, [1,2])
        a,b =  sincos(T(pi/4))
        W = A.Rotator(complex(a,b), zero(T), 1)

        A.turnover(U,V,W)
        U1, V1, W1 = A.turnover(U,V,W)
        @test all(isapprox.([U1,V1,W1]*M , [U,V,W]*M, atol=4eps(T)))

    end

    ## fuse
    for T in Ts
        S = Complex{T}
        U,V  = A.random_rotator.(S, [2,2])
        M = diagm(0 => ones(S, 4))
        UV,Di = A.fuse(U,V)
        alpha = Di.c
        D = diagm(0 => [1, alpha, conj(alpha), 1])
        @test norm(UV*D*M - [U,V]*M) <= sqrt(eps(T))
    end

end

@testset "factorization" begin
    # RDS case first
    ps = [p4, p5, p6]
    Ts = [Float64]
    for p in ps
        for T in Ts
            d = length(p)-1
            state = A.amrvw(p)

            # Cxs = y
            C = state.RF.Ct'
            xs = A.basic_decompose(p)
            M = diagm(0 => ones(T, d+1))
            y = (C*M)*xs
            @test abs(y[1] - norm(xs)) <= sqrt(eps(T)) # check that y[1] is right
            @test norm(y[2:end]) <= sqrt(eps(T))

            # Matrix gives correct eigenvalues
            F = A.Matrix(state)
            rts = eigvals(F)
            @test norm(eigvals(F)  .- [1:d...] ) <= sqrt(eps(T))

            # state.Ct * state.B is I if we skip some
            M = diagm(0 => ones(T, d+1))
            @test norm(state.RF.Ct.x[2:end] * (state.RF.B.x[1:end-1] * M) - I) <= sqrt(eps(T))

            # state.Ct[1] * state.B[end] is twisted
            @test norm((state.RF.Ct.x[1] * (state.RF.B.x[end]*M))[end-1:end, end-1:end] - [0 -1; 1 0]) <= sqrt(eps(T))

            M = diagm(0=> ones(T,d+1))
            state = A.amrvw(p)
            e_1 = vcat(1, zeros(d))
            e_n = vcat(zeros(d-1),1,0)

            xs = A.basic_decompose(p)
            # check that R = Z + x*en'
            Z = state.RF.Ct * (state.RF.B * M) # tested above
            X = state.QF.Q  * (Z + xs*e_n')
            @test norm(sort(eigvals(X))  .- [0:d...] ) <= sqrt(eps(T))
        end
    end

    # CSS case
    ps = [pc]
    Ts = [Float64]
    for p in ps
        for T in Ts
            S = Complex{T}
            d = length(p)-1
            state = A.amrvw(p)

            # Cxs = y
            C = state.RF.Ct'
            xs = A.basic_decompose(p)
            M = diagm(0 => ones(S, d+1))
            y = (C*M)*xs
            @test abs(y[1]) ≈ norm(xs)  # check that y[1] is right
            @test norm(y[2:end]) <= sqrt(eps(T))

            # Matrix gives correct eigenvalues
            F = Matrix(state)
            rts = csort(eigvals(F))
            @test norm(csort(round2(eigvals(F)))  .- S[-1,-im,im,im,1] ) <= sqrt(eps(T))

        end
    end

end

@testset "bulge step" begin

    # RDS First
    ps = [p4, p5, p6]
    Ts = [Float64]

    # check that bulge step works
    for p in ps
        for T in Ts
            m = 2
            state = A.amrvw(p)
            storage =  A.make_storage(state.QF, m)
            ctr = A.make_counter(state)
            A.bulge_step(state.QF, state.RF,  storage, ctr, m)
            d = length(p)-1
            F = A.Matrix(state)
            rts = sort(eigvals(F))
            @test norm(sort(eigvals(F))  .- [1:d...] ) <= 100sqrt(eps(T))

            A.bulge_step(state.QF, state.RF,  storage, ctr, 2)
            allocs = @allocated A.bulge_step(state.QF, state.RF, storage, ctr, 2)
            @test_skip allocs == 0
        end
    end

    # for complex too
    for p in [pc]
        for S in [Complex{Float64}]
            state = A.amrvw(p)
            m =  1
            storage =  A.make_storage(state.QF, m)
            ctr = A.make_counter(state)

            A.bulge_step(state.QF, state.RF, storage, ctr, m)
            allocs = @allocated A.bulge_step(state.QF, state.RF, storage, ctr, m)
            @test_skip allocs == 0
        end
    end
end

@testset "deflation" begin

end

@testset "diagonal block" begin
    # RDS First
    ps = [p4, p5, p6]
    Ts = [Float64]

    for p in ps
        for T in Ts
            state = A.amrvw(p)
            m =  2
            storage =  A.make_storage(state.QF, m)
            ctr = A.make_counter(state)

            A.bulge_step(state.QF, state.RF, storage, ctr, m)
            ctr.it_count += 1
            A.bulge_step(state.QF, state.RF, storage, ctr, m)
            A.create_bulge(state.QF,  state.RF, storage, ctr)

            k = length(p) - 1
            @test Matrix(state.RF)[k,k] ≈ state.RF[k,k]
            @test Matrix(state.RF)[k-1,k] ≈ state.RF[k-1,k]
            @test Matrix(state.RF)[k-1,k-1] ≈ state.RF[k-1,k-1]

            A.diagonal_block!(storage.A, state.QF, state.RF, k-1, k-1)
            @test all(Matrix(state)[k-1:k, k-1:k] .≈ storage.A[1:2, 1:2])
        end
    end

    # Css code

end


text:
using LatestPkg
using Test

@testset "LatestPkg.jl" begin
    # Write your own tests here.
end


text:
using Test
using LocalScopeBenchmarks
using BenchmarkTools
using Statistics

function judge_loosely(t1, t2)
    judge(ratio(median(t1), median(t2)), time_tolerance=0.2)
end

global_x = 1.0

@testset "LocalScopeBenchmarks" begin
    @testset "Basic benchmarks" begin
        x = 1.0
        t1 = @benchmark(sin($x))
        t2 = @localbenchmark(sin(x))
        j = judge_loosely(t1, t2)
        @show j
        @test isinvariant(j)

        t1 = @benchmark($sin($x))
        t2 = @localbenchmark(sin(x))
        j = judge_loosely(t1, t2)
        @show j
        @test isinvariant(j)

        f = sin
        x = 1.0
        t1 = @benchmark($f($x))
        t2 = @localbenchmark(f(x))
        j = judge_loosely(t1, t2)
        @show j
        @test isinvariant(j)
    end

    @testset "Benchmarks with setup" begin
        @testset "Single setup" begin
            t1 = @benchmark sin(x) setup=(x = 2.0)
            t2 = @localbenchmark sin(x) setup=(x = 2.0)
            j = judge_loosely(t1, t2)
            @show j
            @test isinvariant(j)
        end

        @testset "Multiple setups" begin
            t1 = @benchmark atan(x, y) setup=(x = 2.0; y = 1.5)
            t2 = @localbenchmark atan(x, y) setup=(x = 2.0; y = 1.5)
            j = judge_loosely(t1, t2)
            @show j
            @test isinvariant(j)
        end

        @testset "Setups override local vars" begin
            x = 1.0
            @localbenchmark (@assert x == 2.0) setup=(x = 2.0) evals=1
        end

        @testset "Mixed setup and local vars" begin
            x = 1.0
            t1 = @benchmark atan($x, y) setup=(y = 2.0)
            t2 = @localbenchmark atan(x, y) setup=(y = 2.0)
            j = judge_loosely(t1, t2)
            @show j
            @test isinvariant(j)
        end
    end

    @testset "Additional kwargs" begin
        @testset "evals kwarg" begin
            x = 1.0
            t1 = @benchmark sin($x) evals=5
            t2 = @localbenchmark sin(x) evals=5
            j = judge_loosely(t1, t2)
            @show j
            @test isinvariant(j)
        end

        @testset "evals and setup kwargs" begin
            t1 = @benchmark sin(x) setup=(x = 2.0) evals=5
            t2 = @localbenchmark sin(x) setup=(x = 2.0) evals=5
            j = judge_loosely(t1, t2)
            @show j
            @test isinvariant(j)
        end
    end

    @testset "Test that local benchmarks are faster than globals" begin
        t1 = @benchmark sin(global_x) evals=5  # note the lack of $
        t2 = @localbenchmark sin(global_x) evals=5
        j = judge_loosely(t1, t2)
        @show j
        @test isregression(j)
    end

    @testset "Other macros" begin
        x = 1.0
        @localbtime sin(x)
        @localbelapsed sin(x)
    end

    @testset "Interpolated values" begin
        t1 = @benchmark sum($(rand(1000)))
        t2 = @localbenchmark sum($(rand(1000)))
        j = judge_loosely(t1, t2)
        @show j
        @test isinvariant(j)
    end
end


text:
using Gadfly, RDatasets, Distributions

set_default_plot_size(6inch, 3inch)

plot(dataset("ggplot2", "diamonds"), x=:Price, Geom.density)

# manual densities
dist = MixtureModel(Normal, [(0.5, 0.2), (1, 0.1)])
xs = rand(dist, 10^5)
plot(layer(x=xs, Geom.density, Theme(default_color=colorant"orange")),
layer(x=xs, Geom.density(bandwidth=0.003),
Theme(default_color=colorant"blue")),
layer(x=xs, Geom.density(bandwidth=0.25),
Theme(default_color=colorant"purple")),
Guide.manual_color_key("bandwidth", ["auto", "bw=0.003", "bw=0.25"],
["orange", "blue", "purple"]))


text:
# --- Scopes ---
function ctx_pushscope!(ctx, decl = nothing)
    push!(ctx.scopes, symtable_new(decl))
end

function ctx_popscope!(ctx)
    pop!(ctx.scopes)
end

# Returns the symbol for this variable id
# Can return nothing if no symbol with this name already declared
function ctx_fetchscope(ctx, id)
    for i in length(ctx.scopes):-1:1
        if haskey(ctx.scopes[i].syms, id)
            return ctx.scopes[i].syms[id]
        end
    end

    return nothing
end

# Fetches a global symbol
# Might return nothing
function ctx_fetchglobal(ctx, id)
    if haskey(ctx.scopes[1].syms, id)
        return ctx.scopes[1].syms[id]
    end

    return nothing
end

# --- Symbols ---
function ctx_newsymlocal!(ctx::Ctx, decl::Decl, location)
    type = decl.type
    id = decl.id

    @alphaassert length(ctx.scopes) >= 2 decl.location "Cannot declare local symbols outside of a function"

    func_scope = ctx.scopes[2]
    func_scope.decl.nlocals += 1
    position = func_scope.decl.nlocals

    scope = last(ctx.scopes)
    sym = sym_new(k_sym_local, type, id, position)

    @alphaassert !haskey(scope.syms, id) decl.location "Variable named '$id' already declared"

    scope.syms[id] = sym

    return sym
end

function ctx_newsymarg!(ctx::Ctx, type::DeclType, id::String, position::Int, location)
    scope = last(ctx.scopes)
    scope.nargs += 1
    position = scope.nargs
    sym = sym_new(k_sym_arg, type, id, position)

    @alphaassert !haskey(scope.syms, id) location "Argument named '$id' already declared"

    scope.syms[id] = sym
end

function ctx_newsymglobal!(ctx::Ctx, type::DeclType, id::String, location)
    global alphalib_funcs

    sym = sym_new(k_sym_global, type, id)
    scope = first(ctx.scopes)

    @alphaassert !haskey(scope.syms, id) location "Declaration named '$id' already declared"
    @alphaassert id ∉ alphalib_funcs location "Declaration named '$id' reserved for alpha interns, use another name"

    scope.syms[id] = sym
end

# --- Registers ---
scratchregs_new() = [
        r8,
        r9,
        r10,
        r11
    ]

arg_regs = [
        di,
        si,
        dx,
        cx,
        r8,
        r9
    ]

# Allocates a new scratch register
function ctx_newscratch!(ctx::Ctx)
    @assert length(ctx.scratch_regs) != 0 "Not enough scratch registers"

    reg = pop!(ctx.scratch_regs)

    return reg
end

function ctx_freescratch!(ctx::Ctx, reg)
    @assert reg ∉ ctx.scratch_regs "Cannot push reg $(regstr(reg)) to $(ctx.scratch_regs)"

    push!(ctx.scratch_regs, reg)
end

function pushinstr!(section, data::String; indent = true)
    if indent
        data = "    " * data
    end

    push!(section, data)
end

# Pushes code either in data section or current function's text section
function ctx_push!(ctx::Ctx, data::String; code = true, indent = true)
    section = if code ctx.code else ctx.data end

    pushinstr!(section, data, indent = indent)
end

function regstr(reg::Reg)
    global reg2str

    return reg2str[reg]
end

# --- Labels ---
function labelstr(label::Int)
    return ".L$label"
end

function ctx_newlabel!(ctx)
    ctx.n_labels += 1

    return ctx.n_labels
end

# --- Data ---
# Returns the name of the label
function ctx_newdata!(ctx, data)
    ctx.n_data += 1

    name = "D$(ctx.n_data)"

    instruction = "$name:\n$data"
    ctx_push!(ctx, instruction, code=false, indent=false)

    return name
end


text:
module VectorCalculus

export grad, divergence, curl

function gradhelper(f)
    s = size(f,1)
    if s == 1

        g = zeros(f)

    elseif s == 2

        f1 = f[1,:] # Extracts the first "row" from all dimensions.
        f2 = f[2,:]
        n = length(f1)
        f21 = reshape(f2-f1, 1, n)
        g = reshape([f21; f21], size(f))

    elseif s == 3

        f1 = f[1,:]
        f2 = f[2,:]
        f3 = f[3,:]
        n = length(f1)
        f21 = reshape(f2-f1, 1, n)
        f31 = reshape(f3-f1, 1, n)
        f32 = reshape(f3-f2, 1, n)
        g = reshape([f21 ; f31./2 ; f32], size(f))

    else

        f1 = f[1,:]
        f2 = f[2,:]
        fe = f[end,:]
        fm = f[end-1,:]
        n = length(f1)
        f21 = reshape(f2-f1, 1, n)
        fem = reshape(fe-fm, 1, n)
        g = reshape([f21 ; (f[3:end,:] - f[1:end-2,:])./2 ; fem], size(f))

    end

    return g

end

function gradhelper(f,h,n,k)

     s = size(f,1)

     if s == 1

         g = zeros(f)

     elseif s == 2

         hi = [h[2]-h[1] ; h[k] - h[k-1]]
         f1 = f[1,:]
         f2 = f[2,:]
         fk = f[k,:]
         fm = f[k-1,:]
         n = length(f1)
         f21 = reshape(f2-f1, 1, n)
         fkm = reshape(fk-fm, 1, n)
         g = reshape([f21 ; fkm]./hi[:,ones(Int, div(n, k), 1)], size(f))

     elseif s == 3 && k == 3

         hi = [h[2]-h[1] ; (h[3:k] - h[1:k-2]) ; h[k] - h[k-1]]
         f1 = f[1,:]
         f2 = f[2,:]
         f3 = f[3,:]
         n = length(f1)
         f21 = reshape(f2-f1, 1, n)
         f31 = reshape(f3-f1, 1, n)
         f32 = reshape(f3-f2, 1, n)
         g = reshape([f21; f31 ; f32] ./ hi[:,ones(Int, div(n, k), 1)], size(f))

     else

         hi = [h[2]-h[1] ; (h[3:k] - h[1:k-2]) ; h[k] - h[k-1]]

         f1 = f[1,:]
         f2 = f[2,:]
         f3 = f[3:k,:]
         f4 = f[1:k-2,:]
         fk = f[k,:]
         fm = f[k-1,:]
         n = length(f1)
         m = size(f3,1)
         f21 = reshape(f2-f1,1,n)
         fkm = reshape(fk-fm,1,n)

         if m == 3

             f34 = reshape(f3-f4,1,n)
             g = reshape([f21 ; f34 ; fkm] ./ hi[:,ones(Int, div(n, k), 1)], size(f))

         else

             g = reshape([f21 ; f3 - f4 ; fkm] ./ hi[:,ones(Int, div(n, k), 1)], size(f))

         end

     end

     return g

end

function grad(f)

    nd = ndims(f)
    grads = Array(Array{eltype(f), nd}, nd)

    for i=1:nd
        dims = [2; 1; collect(3:nd)]
        dims[1] = dims[i]
        dims[i] = 2
        grads[i] = ipermutedims(gradhelper(permutedims(f, dims)), dims)
    end

    return grads
end

function grad{T<:Number}(f, h::Vector{T}...)

    hn = length(h)
    nd = ndims(f)
    n = length(f)
    a = size(f)

    @assert nd == hn

    grads = Array(Array{eltype(f), nd}, nd)

    for i=1:nd
        dims = [2; 1; collect(3:nd)]
        dims[1] = dims[i]
        dims[i] = 2
        hi = h[i]
        k = a[dims[i]]
        if length(hi) == 1
            grads[i] = ipermutedims(gradhelper(permutedims(f, dims)), dims)./hi[1]
        elseif length(hi) == a[i]
            grads[i] = ipermutedims(gradhelper(permutedims(f, dims), hi, n, k), dims)
        end
    end

    return grads
end

function grad{T<:Number}(f, h::T...)

    hn = length(h)
    nd = ndims(f)
    n = length(f)
    a = size(f)

    @assert nd == hn

    grads = Array(Array{eltype(f), nd}, nd)

    for i=1:nd
        dims = [2; 1; collect(3:nd)]
        dims[1] = dims[i]
        dims[i] = 2
        hi = h[i]
        k = a[dims[i]]
        grads[i] = ipermutedims(subgrad(permutedims(f,dims)), dims)./hi
    end

    return grads
end

function divergence{T<:Number}(X::Vector{T}, Y::Vector{T}, Z::Vector{T}, U::Array{T,3}, V::Array{T,3}, W::Array{T,3})

    ux, uy, uz = grad(U, X, Y, Z)
    vx, vy, vz = grad(V, X, Y, Z)
    wx, wy, wz = grad(W, X, Y, Z)

    return ux + vy + wz
end

function divergence{T<:Number}(X::Array{T,3}, Y::Array{T,3}, Z::Array{T,3}, U::Array{T,3}, V::Array{T,3}, W::Array{T,3})

    x = X[1, :, 1]
    y = Y[:, 1, 1]
    z = Z[1, 1, :]

    ux, uy, uz = grad(U, x, y, z)
    vx, vy, vz = grad(V, x, y, z)
    wx, wy, wz = grad(W, x, y, z)

    return ux + vy + wz
end

function divergence{T<:Number}(X::Vector{T}, Y::Vector{T}, U::Matrix{T}, V::Matrix{T})

    ux, uy = grad(U, X, Y)
    vx, vy = grad(V, X, Y)

    return ux + vy
end

function divergence{T<:Number}(X::Matrix{T}, Y::Matrix{T}, U::Matrix{T}, V::Matrix{T})

    x = X[1, :]
    y = Y[:, 1]

    ux, uy = grad(U, x, y)
    vx, vy = grad(V, x, y)

    return ux + vy
end

function divergence{T<:Number}(U::Array{T,3}, V::Array{T,3}, W::Array{T,3})

    ux, uy, uz = grad(U)
    vx, vy, vz = grad(V)
    wx, wy, wz = grad(W)

    return ux + vy + wz
end

function divergence{T<:Number}(U::Matrix{T}, V::Matrix{T})

    ux, uy = grad(U)
    vx, vy = grad(V)

    return ux + vy
end

function curl{T<:Number}(X::Vector{T}, Y::Vector{T}, Z::Vector{T}, U::Array{T,3}, V::Array{T,3}, W::Array{T,3})

    ux, uy, uz = grad(U, X, Y, Z)
    vx, vy, vz = grad(V, X, Y, Z)
    wx, wy, wz = grad(W, X, Y, Z)

    cx = wy - vz
    cy = uz - wx
    cz = vx - uy
    ca = (cx.*U + cy.*V + cz.*W) ./ (2sqrt(U.^2 + V.^2 + W.^2))

    return cx, cy, cz, ca
end

function curl{T<:Number}(X::Array{T,3}, Y::Array{T,3}, Z::Array{T,3}, U::Array{T,3}, V::Array{T,3}, W::Array{T,3})

    x = X[1, :, 1]
    y = Y[:, 1, 1]
    z = Z[1, 1, :]

    ux, uy, uz = grad(U, x, y, z)
    vx, vy, vz = grad(V, x, y, z)
    wx, wy, wz = grad(W, x, y, z)

    cx = wy - vz
    cy = uz - wx
    cz = vx - uy
    ca = (cx.*U + cy.*V + cz.*W) ./ (2sqrt(U.^2 + V.^2 + W.^2))

    return cx, cy, cz, ca
end

function curl{T<:Number}(X::Vector{T}, Y::Vector{T}, U::Matrix{T}, V::Matrix{T})

    ux, uy = grad(U, X, Y)
    vx, vy = grad(V, X, Y)

    cz = vx - uy
    ca = cz./2

    return cz, ca
end

function curl{T<:Number}(X::Matrix{T}, Y::Matrix{T}, U::Matrix{T}, V::Matrix{T})

    x = X[1, :]
    y = Y[:, 1]

    ux, uy = grad(U, x, y)
    vx, vy = grad(V, x, y)

    cz = vx - uy
    ca = cz./2

    return cz, ca
end

function curl{T<:Number}(U::Array{T,3}, V::Array{T,3}, W::Array{T,3})

    ux, uy, uz = grad(U)
    vx, vy, vz = grad(V)
    wx, wy, wz = grad(W)

    cx = wy - vz
    cy = uz - wx
    cz = vx - uy
    ca = (cx.*U + cy.*V + cz.*W) ./ (2sqrt(U.^2 + V.^2 + W.^2))

    return cx, cy, cz, ca
end

function curl{T<:Number}(U::Matrix{T}, V::Matrix{T})

    ux, uy = grad(U)
    vx, vy = grad(V)

    cz = vx - uy
    ca = cz./2

    return cz, ca
end

end


text:
function getFFF()
    factorUrl =
        "http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_daily_TXT.zip"
    (data, dataNames, varnames) = readFamaFrenchRaw(factorUrl)

    newVarnames = [:MktRf, :SMB, :HML, :RF]
    rename!(data[1].vals, names(data[1].vals), newVarnames)

    fileName = joinpath(Pkg.dir("EconDatasets"), "data", "FFF.csv")
    writeTimedata(fileName, data[1])
end


text:
############################################################
# abstract types hierarchy
#
# joAbstractOperator
# * joAbstractLinearOperator
#   * joAbstractParallelableLinearOperator
# * joAbstractFosterLinearOperator
#   * joAbstractLinearOperatorInplace
# * joAbstractSAparallelOperator
#  * joAbstractSAparallelToggleOperator
#  * joAbstractSAparallelLinearOperator
# * joAbstractDAparallelOperator
#  * joAbstractDAparallelToggleOperator
#  * joAbstractDAparallelLinearOperator
#
############################################################

############################################################
# joAbstractOperator #######################################

export joAbstractOperator, joAbstractOperatorException

# type definition
abstract type joAbstractOperator end

# type exception
struct joAbstractOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractLinearOperator #################################

export joAbstractLinearOperator, joAbstractLinearOperatorException

# type definition
abstract type joAbstractLinearOperator{DDT<:Number,RDT<:Number} <: joAbstractOperator end

# type exception
struct joAbstractLinearOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractParallelableLinearOperator #####################

export joAbstractParallelableLinearOperator, joAbstractParallelableLinearOperatorException

# type definition
abstract type joAbstractParallelableLinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractLinearOperator{DDT,RDT} end

# type exception
struct joAbstractParallelableLinearOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractFosterLinearOperator ###########################

export joAbstractFosterLinearOperator, joAbstractFosterLinearOperatorException

# type definition
abstract type joAbstractFosterLinearOperator{DDT<:Number,RDT<:Number} <: joAbstractOperator end

# type exception
struct joAbstractFosterLinearOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractLinearOperatorInplace ##########################

export joAbstractLinearOperatorInplace, joAbstractLinearOperatorInplaceException

# type definition
abstract type joAbstractLinearOperatorInplace{DDT<:Number,RDT<:Number} <: joAbstractFosterLinearOperator{DDT,RDT} end

# type exception
struct joAbstractLinearOperatorInplaceException <: Exception
    msg :: String
end

############################################################
# joAbstractSAparallelOperator #############################

export joAbstractSAparallelOperator, joAbstractSAparallelOperatorException

# type definition
abstract type joAbstractSAparallelOperator{DDT<:Number,RDT<:Number,N} <: joAbstractOperator end

# type exception
struct joAbstractSAparallelOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractSAparallelToggleOperator #######################

export joAbstractSAparallelToggleOperator, joAbstractSAparallelToggleOperatorException

# type definition
abstract type joAbstractSAparallelToggleOperator{DDT<:Number,RDT<:Number,N} <: joAbstractSAparallelOperator{DDT,RDT,N} end

# type exception
struct joAbstractSAparallelToggleOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractSAparallelLinearOperator #######################

export joAbstractSAparallelLinearOperator, joAbstractSAparallelLinearOperatorException

# type definition
abstract type joAbstractSAparallelLinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractSAparallelOperator{DDT,RDT,N} end

# type exception
struct joAbstractSAparallelLinearOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractDAparallelOperator #############################

export joAbstractDAparallelOperator, joAbstractDAparallelOperatorException

# type definition
abstract type joAbstractDAparallelOperator{DDT<:Number,RDT<:Number,N} <: joAbstractOperator end

# type exception
struct joAbstractDAparallelOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractDAparallelToggleOperator #######################

export joAbstractDAparallelToggleOperator, joAbstractDAparallelToggleOperatorException

# type definition
abstract type joAbstractDAparallelToggleOperator{DDT<:Number,RDT<:Number,N} <: joAbstractDAparallelOperator{DDT,RDT,N} end

# type exception
struct joAbstractDAparallelToggleOperatorException <: Exception
    msg :: String
end

############################################################
# joAbstractDAparallelLinearOperator #######################

export joAbstractDAparallelLinearOperator, joAbstractDAparallelLinearOperatorException

# type definition
abstract type joAbstractDAparallelLinearOperator{DDT<:Number,RDT<:Number,N} <: joAbstractDAparallelOperator{DDT,RDT,N} end

# type exception
struct joAbstractDAparallelLinearOperatorException <: Exception
    msg :: String
end