export lowrankvortexassim


# Create a function to perform the sequential assimilation for any sequential filter SeqFilter
function lowrankvortexassim(algo::SeqFilter, X, tspan::Tuple{S,S}, config::VortexConfig, data::SyntheticData;
	                        withfreestream::Bool = false,
							rxdefault::Int64 = 100, rydefault::Int64 = 100, israndomized::Bool=false, P::Parallel = serial) where {S<:Real}

	# Define the additive Inflation
	ϵX = config.ϵX
	ϵΓ = config.ϵΓ
	β = config.β
	ϵY = config.ϵY

	Ny = size(config.ss,1)

	ϵx = RecipeInflation([ϵX; ϵΓ])
	ϵmul = MultiplicativeInflation(β)
	# Set different times
	Δtobs = algo.Δtobs
	Δtdyn = algo.Δtdyn
	t0, tf = tspan
	step = ceil(Int, Δtobs/Δtdyn)

	n0 = ceil(Int64, t0/Δtobs) + 1
	J = (tf-t0)/Δtobs
	Acycle = n0:n0+J-1

	# Array dimensions
	Nypx, Ne = size(X)
	Nx = Nypx - Ny
	ystar = zeros(Ny)

	freestream = Freestream(0.0*im)

	cachevels = allocate_velocity(state_to_lagrange(X[Ny+1:Ny+Nx,1], config))

	h(x, t) = measure_state(x, t, config; withfreestream = true)
	press_itp = CubicSplineInterpolation((LinRange(real(config.ss[1]), real(config.ss[end]), length(config.ss)),
	                                   t0:config.Δt:tf), data.yt, extrapolation_bc =  Line())

	yt(t) = press_itp(real.(config.ss), t)
	Xf = Array{Float64,2}[]
	push!(Xf, copy(state(X, Ny, Nx)))

	Xa = Array{Float64,2}[]
	push!(Xa, copy(state(X, Ny, Nx)))

	# Run particle filter
	@showprogress for i=1:length(Acycle)
		# Forecast step
		@inbounds for j=1:step
			tj = t0+(i-1)*Δtobs+(j-1)*Δtdyn
			X, _ = vortex(X, tj, Ny, Nx, cachevels, config; withfreestream = true)
		end

		push!(Xf, deepcopy(state(X, Ny, Nx)))

		# Get real measurement
		ystar .= yt(t0+i*Δtobs)

		# Perform state inflation
		ϵmul(X, Ny+1, Ny+Nx)
		ϵx(X, Ny, Nx, config)

		# Filter state
		if algo.isfiltered == true
			@inbounds for i=1:Ne
				x = view(X, Ny+1:Ny+Nx, i)
				# x .= filter_state!(x, config)
			end
		end

		observe(h, X, t0+i*Δtobs, Ny, Nx; P = P)

		ϵ = algo.ϵy.σ*randn(Ny, Ne) .+ algo.ϵy.m

		# Low-rank basis
		Cx = zeros(Nx, Nx)
		Cy  = zeros(Ny, Ny)


		Dx = Diagonal(std(X[Ny+1:Ny+Nx, :]; dims = 2)[:,1])
		# Dx = I
		# Dϵ = Diagonal(std(ϵ; dims = 2)[:,1])
		Dϵ = config.ϵY*I
		# Dϵ = I

		@inbounds Threads.@threads for j=1:Ne
			# J_AD = AD_symmetric_jacobian_pressure(config.ss, vcat(state_to_lagrange(X[Ny+1:Ny+Nx,j], config)...), t0+i*Δtobs)
			J = analytical_jacobian_pressure(config.ss, vcat(state_to_lagrange(X[Ny+1:Ny+Nx,j], config)...), freestream, t0+i*Δtobs)
			Jj = J[:,1:3*config.Nv]
			Cx .+= 1/(Ne-1)*(inv(Dϵ)*Jj*Dx)'*(inv(Dϵ)*Jj*Dx)
			Cy .+= 1/(Ne-1)*(inv(Dϵ)*Jj*Dx)*(inv(Dϵ)*Jj*Dx)'
		end


		ry = min(Ny, rydefault)
		rx = min(Nx, rxdefault)

		# Lbx, V = pheig(Symmetric(Cx), rank = rx)
		# V = reverse(V, dims = 2)
		# Λx, V = eigen(Symmetric(Cx); sortby = λ -> -λ)
		V, Λx, _ = svd(Symmetric(Cx))
		V = V[:,1:rx]
		# @show norm(V*V'-I), norm(V'*V-I)

		# Lby, U = pheig(Symmetric(Cy), rank = ry)
		# U = reverse(U, dims = 2)
		# Λy, U = eigen(Symmetric(Cy); sortby = λ -> -λ)
		U, Λy, _ = svd(Symmetric(Cy))
		U = U[:,1:ry]
		# U = I

		Xbreve = V'*(inv(Dx)*(X[Ny+1:Ny+Nx,:] .- mean(X[Ny+1:Ny+Nx,:]; dims = 2)[:,1]))
		Xbrevepert = (1/sqrt(Ne-1))*(Xbreve .- mean(Xbreve; dims = 2)[:,1])

		HXbreve = U'*(inv(Dϵ)*(X[1:Ny,:] .- mean(X[1:Ny,:]; dims = 2)[:,1]))
		HXbrevepert = (1/sqrt(Ne-1))*(HXbreve .- mean(HXbreve; dims = 2)[:,1])

		ϵbreve = U'*(inv(Dϵ)*(ϵ .- mean(ϵ; dims =2)[:,1]))
		ϵbrevepert = (1/sqrt(Ne-1))*(ϵbreve .- mean(ϵbreve; dims = 2)[:,1])

		"Low-rank analysis step with representers, Evensen, Leeuwen et al. 1998"

		# K̆ = Xbrevepert*HXbrevepert'*inv(HXbrevepert*HXbrevepert' + ϵbrevepert*ϵbrevepert')
		# Kcxcy =  Dx*V*K̆*U'*inv(Dϵ)
		# Xpert = (1/sqrt(Ne-1))*(X[Ny+1:Ny+Nx,:] .- mean(X[Ny+1:Ny+Nx,:]; dims = 2)[:,1])
		# HXpert = (1/sqrt(Ne-1))*(X[1:Ny,:] .- mean(X[1:Ny,:]; dims = 2)[:,1])
		# ϵpert = (1/sqrt(Ne-1))*(ϵ .- mean(ϵ; dims = 2)[:,1])
		# Kenkf = Xpert*HXpert'*inv(HXpert*HXpert'+ϵpert*ϵpert')
		# K̆enkf = V'*inv(Dx)*Kenkf*Dϵ*U
		# @show norm(Kcxcy - Kenkf), norm(Kcxcy - Kenkf)/norm(Kenkf)
		# @show cumsum(svd(Kenkf).S) ./ sum(svd(Kenkf).S)
		# @show norm(K̆ - K̆enkf), norm(K̆ - K̆enkf)/norm(K̆enkf)

		b̆ = (HXbrevepert*HXbrevepert' + ϵbrevepert*ϵbrevepert')\(U'*(Dϵ\(ystar .- (X[1:Ny,:] + ϵ))))
		view(X,Ny+1:Ny+Nx,:) .+= Dx*V*(Xbrevepert*HXbrevepert')*b̆

		# X = algo(X, ystar, t0+i*Δtobs)

		# Filter state
		if algo.isfiltered == true
			@inbounds for i=1:Ne
				x = view(X, Ny+1:Ny+Nx, i)
				# x .= filter_state!(x, config)
			end
		end

		push!(Xa, deepcopy(state(X, Ny, Nx)))
		end

	return Xf, Xa
end
