using LocalAnisotropies
using GeoStats
using Test
import LocalAnisotropies: rotmat

@testset "LocalAnisotropies.jl" begin
    # convert data to LocalAnisotropy
    dummy = georef((az=1:10, r1=1:10, r2=1:10), PointSet(rand(2,10)))
    pars  = localanisotropies(dummy, [:az], [:r1,:r2], :GSLIB)

    # convert between different rotation conventions
    angs1 = convertangles([30,30,30], :GSLIB, :Datamine)
    angs2 = convertangles(angs1, :Datamine, :GSLIB)
	rmat  = rotmat([1,1,1], [30,30,30], :GSLIB)
    @test all(rmat .≈ rotmat([1,1,1], angs1, :Datamine))
    @test all(rmat .≈ rotmat([1,1,1], angs2, :GSLIB))
    angs_ = convertangles.(pars.rotation, :GSLIB)

    # pars_ = convertpars(pars, convention=:Datamine)
    # exportpars("C:\\Users\\test.csv", pars, :GSLIB)

	# local anisotropies from pointset coordinates
	data  = rmat[1] * vcat(rand(2,100),zeros(1,100))
	pset  = PointSet(data)
	nhood = KNearestSearch(pset, 10)
    gpars = localanisotropies(Geometric, nhood, simplify=true)
    gpars = localanisotropies(Geometric, nhood, simplify=false)

    grid2d = (10,10)
    grid3d = (10,10,5)

    for dims in (grid2d, grid3d)
        # reference scenario for tests
		R = length(dims)
		D = if R==2
    		georef((P=[sin(i)+j for i in 1:dims[1], j in 1:dims[2]],))
    	else
    		georef((P=[sin(i)+j+k for i in 1:dims[1], j in 1:dims[2], k in 1:dims[3]],))
    	end

        n = round(Int, 0.2*prod(dims))
        S = sample(D, n, replace=false)
        S = georef(values(S), PointSet(centroid.(domain(S))))
        G = CartesianGrid(dims...)
        searcher = KNearestSearch(G, 10)

        # get local anisotropies
        lpars = localanisotropies(Gradients, D, :P, 3)

        # rescale magnitude and interpolate local anisotropies
        lpars = rescale_magnitude(lpars, r2=(0.2,1.0))
        lpars = reference_magnitude(lpars, :Y)
        lpars = smoothpars(lpars, searcher)

		spars = nnpars(lpars, D, S)

        # interpolate in a coarser grid
        G_ = if R==2
			CartesianGrid((5,5),(0.5,0.5),(2.0,2.0))
		else
			CartesianGrid((5,5,3),(0.5,0.5,0.5),(2.0,2.0,2.0))
		end
        lpars_ = idwpars(lpars, searcher, G_)

        # Estimation problem
        P = EstimationProblem(S, G, :P)
        γ = NuggetEffect(0.1) + 0.9*ExponentialVariogram(range=60.0)

        # LocalKriging (MW)
        MW = LocalKriging(:P => (variogram=γ, localaniso=lpars, method=:MovingWindows))
        s1 = solve(P, MW)

        # LocalKriging (KC)
        KC = LocalKriging(:P => (variogram=γ, localaniso=lpars, method=:KernelConvolution))
        s2 = solve(P, KC)

        # Spatial deformation: anisotropic distances
        Sd1, Dd1 = deformspace(S, G, lpars, AnisoDistance, anchors=1500)
        Pd1 = EstimationProblem(Sd1, Dd1, :P)
        s3 = solve(Pd1, Kriging(:P => (variogram=γ,)))
		x3 = to_3d(s3)

        # Spatial deformation: anisotropic variogram distances
        Sd2, Dd2 = deformspace(S, G, lpars, KernelVariogram, γ, anchors=1500)
        Pd2 = EstimationProblem(Sd2, Dd2, :P)
        s4 = solve(Pd2, Kriging(:P => (variogram=γ,)))
		x4 = to_3d(s4)

        # Spatial deformation: geodesic anisotropic distances
        LDa = graph(S, G, lpars, AnisoDistance, searcher)
        Sd3, Dd3 = deformspace(LDa, GraphDistance, anchors=1500)
        Pd3 = EstimationProblem(Sd3, Dd3, :P)
        s5 = solve(Pd3, Kriging(:P => (variogram=γ,)))
		x5 = to_3d(s5)

        # Spatial deformation: geodesic anisotropic variogram distances
        LDv = graph(S, G, lpars, KernelVariogram, γ, searcher)
        Sd4, Dd4 = deformspace(LDv, GraphDistance, anchors=1500)
        Pd4 = EstimationProblem(Sd4, Dd4, :P)
        s6 = solve(Pd4, Kriging(:P => (variogram=γ,)))
		x6 = to_3d(s6)
    end
end
