diff --git a/Project.toml b/Project.toml
index de415795e..19fcc4304 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,10 +4,12 @@ authors = ["ludoro <ludovicobessi@gmail.com>"]
 version = "1.1.2"
 
 [deps]
+Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
+GaussianMixtures = "cc18c42c-b769-54ff-9e2a-b28141a64aae"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LatinHypercubeSampling = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/Kriging.jl b/src/Kriging.jl
index bc948acac..d5bccba70 100644
--- a/src/Kriging.jl
+++ b/src/Kriging.jl
@@ -94,13 +94,12 @@ Constructor for type Kriging.
 -p: value between 0 and 2 modelling the
    smoothness of the function being approximated, 0-> rough  2-> C^infinity
 """
-function Kriging(x,y,lb::Number,ub::Number;p=1.0)
+function Kriging(x,y,lb::Number,ub::Number;p=1.0,theta=1.0)
     if length(x) != length(unique(x))
         println("There exists a repetion in the samples, cannot build Kriging.")
         return
     end
     mu,b,sigma,inverse_of_R = _calc_kriging_coeffs(x,y,p)
-    theta = 1.0
     Kriging(x,y,lb,ub,p,theta,mu,b,sigma,inverse_of_R)
 end
 
diff --git a/src/MOE.jl b/src/MOE.jl
new file mode 100644
index 000000000..b6a2e9f48
--- /dev/null
+++ b/src/MOE.jl
@@ -0,0 +1,268 @@
+using Clustering
+using GaussianMixtures
+using LinearAlgebra
+
+mutable struct MOE{X,Y,L,U,S,K,M,V,W} <: AbstractSurrogate
+    x::X
+    y::Y
+    lb::L
+    ub::U
+    local_surr::S
+    k::K
+    means::M
+    varcov::V
+    weights::W
+end
+
+
+#Radial structure:
+function RadialBasisStructure(;radial_function,scale_factor,sparse)
+    return (name = "RadialBasis", radial_function = radial_function, scale_factor = scale_factor, sparse = sparse)
+end
+
+#Kriging structure:
+function KrigingStructure(;p,theta)
+    return (name = "Kriging", p = p, theta = theta)
+end
+
+#Linear structure
+function LinearStructure()
+    return (name = "LinearSurrogate")
+end
+
+#InverseDistance structure
+function InverseDistanceStructure(;p)
+    return (name = "InverseDistanceSurrogate", p = p)
+end
+
+#Lobachesky structure
+function LobacheskyStructure(;alpha,n,sparse)
+    return (name = "LobacheskySurrogate", alpha = alpha, n = n, sparse = sparse)
+end
+
+function NeuralStructure(;model,loss,opt,n_echos)
+    return (name ="NeuralSurrogate", model = model ,loss = loss,opt = opt,n_echos = n_echos)
+end
+
+function RandomForestStructure(;num_round)
+    return (name = "RandomForestSurrogate", num_round = num_round)
+end
+
+function SecondOrderPolynomialStructure()
+    return (name = "SecondOrderPolynomialSurrogate")
+end
+
+function WendlandStructure(; eps, maxiters, tol)
+    return (name = "Wendland", eps = eps, maxiters = maxiters, tol = tol)
+end
+
+
+
+function MOE(x,y,lb::Number,ub::Number; k::Int = 2, local_kind = [RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0,sparse = false),RadialBasisStructure(radial_function = cubicRadial, scale_factor=1.0, sparse = false)])
+    if k != length(local_kind)
+        throw("Number of mixtures = $k is not equal to length of local surrogates")
+    end
+    n = length(x)
+    # find weight, mean and variance for each mixture
+    # For GaussianMixtures I need nxd matrix
+    X_G = reshape(x,(n,1))
+    moe_gmm = GaussianMixtures.GMM(k,X_G)
+    weights = GaussianMixtures.weights(moe_gmm)
+    means = GaussianMixtures.means(moe_gmm)
+    variances = moe_gmm.Σ
+
+
+    #cluster the points
+    #For clustering I need dxn matrix
+    X_C = reshape(x,(1,n))
+    KNN = kmeans(X_C, k)
+    x_c = [ Array{eltype(x)}(undef,0) for i = 1:k]
+    y_c = [ Array{eltype(y)}(undef,0) for i = 1:k]
+    a = assignments(KNN)
+    @inbounds for i = 1:n
+        pos = a[i]
+        append!(x_c[pos],x[i])
+        append!(y_c[pos],y[i])
+    end
+
+    local_surr = Dict()
+    for i = 1:k
+        if local_kind[i][1] == "RadialBasis"
+            #fit and append to local_surr
+            my_local_i = RadialBasis(x_c[i],y_c[i],lb,ub,rad = local_kind[i].radial_function, scale_factor = local_kind[i].scale_factor, sparse = local_kind[i].sparse)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "Kriging"
+            my_local_i = Kriging(x_c[i], y_c[i],lb,ub, p = local_kind[i].p, theta = local_kind[i].p)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i] == "LinearSurrogate"
+            my_local_i = LinearSurrogate(x_c[i], y_c[i],lb,ub)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "InverseDistanceSurrogate"
+            my_local_i = InverseDistanceSurrogate(x_c[i], y_c[i],lb,ub, local_kind[i].p)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "LobacheskySurrogate"
+            my_local_i = LobacheskySurrogate(x_c[i], y_c[i],lb,ub,alpha = local_kind[i].alpha , n = local_kind[i].n, sparse = local_kind[i].sparse)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "NeuralSurrogate"
+            my_local_i = NeuralSurrogate(x_c[i], y_c[i],lb,ub, model = local_kind[i].model , loss = local_kind[i].loss ,opt = local_kind[i].opt ,n_echos = local_kind[i].n_echos)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "RandomForestSurrogate"
+            my_local_i = RandomForestSurrogate(x_c[i], y_c[i],lb,ub, num_round = local_kind[i].num_round)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i] == "SecondOrderPolynomialSurrogate"
+            my_local_i = SecondOrderPolynomialSurrogate(x_c[i], y_c[i],lb,ub)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "Wendland"
+            my_local_i = Wendand(x_c[i], y_c[i],lb,ub, eps = local_kind[i].eps, maxiters = local_kind[i].maxiters, tol = local_kind[i].tol)
+            local_surr[i] = my_local_i
+        else
+            throw("A surrogate with name "* local_kind[i][1] *" does not exist or is not currently supported with MOE.")
+        end
+    end
+    return MOE(x,y,lb,ub,local_surr,k,means,variances,weights)
+end
+
+function MOE(x,y,lb,ub; k::Int = 2,
+            local_kind = [RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0, sparse = false),RadialBasisStructure(radial_function = cubicRadial, scale_factor=1.0, sparse = false)])
+
+    n = length(x)
+    d = length(lb)
+    #GMM parameters:
+    X_G = collect(reshape(collect(Base.Iterators.flatten(x)), (d,n))')
+    my_gmm = GMM(k,X_G,kind = :full)
+    weights = my_gmm.w
+    means = my_gmm.μ
+    varcov = my_gmm.Σ
+
+    #cluster the points
+    X_C = collect(reshape(collect(Base.Iterators.flatten(x)), (d,n)))
+    KNN = kmeans(X_C, k)
+    x_c = [ Array{eltype(x)}(undef,0) for i = 1:k]
+    y_c = [ Array{eltype(y)}(undef,0) for i = 1:k]
+    a = assignments(KNN)
+    @inbounds for i = 1:n
+        pos = a[i]
+        x_c[pos] = vcat(x_c[pos],x[i])
+        append!(y_c[pos],y[i])
+    end
+
+    local_surr = Dict()
+    for i = 1:k
+        if local_kind[i][1] == "RadialBasis"
+            #fit and append to local_surr
+            my_local_i = RadialBasis(x_c[i],y_c[i],lb,ub,rad = local_kind[i].radial_function, scale_factor = local_kind[i].scale_factor, sparse = local_kind[i].sparse)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "Kriging"
+            my_local_i = Kriging(x_c[i], y_c[i],lb,ub, p = local_kind[i].p, theta = local_kind[i].p)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i] == "LinearSurrogate"
+            my_local_i = LinearSurrogate(x_c[i], y_c[i],lb,ub)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "InverseDistanceSurrogate"
+            my_local_i = InverseDistanceSurrogate(x_c[i], y_c[i],lb,ub, local_kind[i].p)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "LobacheskySurrogate"
+            my_local_i = LobacheskySurrogate(x_c[i], y_c[i],lb,ub,alpha = local_kind[i].alpha , n = local_kind[i].n, sparse = local_kind[i].sparse)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "NeuralSurrogate"
+            my_local_i = NeuralSurrogate(x_c[i], y_c[i],lb,ub, model = local_kind[i].model , loss = local_kind[i].loss ,opt = local_kind[i].opt ,n_echos = local_kind[i].n_echos)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "RandomForestSurrogate"
+            my_local_i = RandomForestSurrogate(x_c[i], y_c[i],lb,ub, num_round = local_kind[i].num_round)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i] == "SecondOrderPolynomialSurrogate"
+            my_local_i = SecondOrderPolynomialSurrogate(x_c[i], y_c[i],lb,ub)
+            local_surr[i] = my_local_i
+
+        elseif local_kind[i][1] == "Wendland"
+            my_local_i = Wendand(x_c[i], y_c[i],lb,ub, eps = local_kind[i].eps, maxiters = local_kind[i].maxiters, tol = local_kind[i].tol)
+            local_surr[i] = my_local_i
+        else
+            throw("A surrogate with name "* local_kind[i][1] *" does not exist or is not currently supported with MOE.")
+        end
+    end
+    return MOE(x,y,lb,ub,local_surr,k,means,varcov,weights)
+end
+
+
+function _prob_x_in_i(x::Number,i,mu,varcov,alpha,k)
+    num = (1/sqrt(varcov[i]))*alpha[i]*exp(-0.5(x-mu[i])*(1/varcov[i])*(x-mu[i]))
+    den = sum([(1/sqrt(varcov[j]))*alpha[j]*exp(-0.5(x-mu[j])*(1/varcov[j])*(x-mu[j]))  for j = 1:k])
+    return num/den
+end
+
+function _prob_x_in_i(x,i,mu,varcov,alpha,k)
+    num = (1/sqrt(det(varcov[i])))*alpha[i]*exp(-0.5*(x .- mu[i,:])'*(inv(varcov[i]))*(x .- mu[i,:]))
+    den = sum([(1/sqrt(det(varcov[j])))*alpha[j]*exp(-0.5*(x .- mu[j,:])'*(inv(varcov[j]))*(x .- mu[j,:])) for j = 1:k])
+    return num/den
+end
+
+function (moe::MOE)(val)
+    return prod([moe.local_surr[i](val)*_prob_x_in_i(val,i,moe.means,moe.varcov,moe.weights,moe.k) for i = 1:moe.k])
+end
+
+
+function add_point!(moe::MOE,x_new,y_new)
+    if length(moe.x[1]) == 1
+        #1D
+        moe.x = vcat(moe.x,x_new)
+        moe.y = vcat(moe.y,y_new)
+        n = length(moe.x)
+
+        #New mixture parameters
+        X_G = reshape(moe.x,(n,1))
+        moe_gmm = GaussianMixtures.GMM(moe.k,X_G)
+        moe.weights = GaussianMixtures.weights(moe_gmm)
+        moe.means = GaussianMixtures.means(moe_gmm)
+        moe.varcov = moe_gmm.Σ
+
+        #Find cluster of new point(s):
+        n_added = length(x_new)
+        X_C = reshape(moe.x,(1,n))
+        KNN = kmeans(X_C, moe.k)
+        a = assignments(KNN)
+        #Recalculate only relevant surrogates
+        for i = 1:n_added
+            pos = a[n-n_added+i]
+            add_point!(moe.local_surr[i],moe.x[n-n_added+i],moe.y[n-n_added+i])
+        end
+    else
+        #ND
+        moe.x = vcat(moe.x,x_new)
+        moe.y = vcat(moe.y,y_new)
+        n = length(moe.x)
+        d = length(moe.lb)
+        #New mixture parameters
+        X_G = collect(reshape(collect(Base.Iterators.flatten(moe.x)), (d,n))')
+        my_gmm = GMM(moe.k,X_G,kind = :full)
+        moe.weights = my_gmm.w
+        moe.means = my_gmm.μ
+        moe.varcov = my_gmm.Σ
+
+        #cluster the points
+        X_C = collect(reshape(collect(Base.Iterators.flatten(moe.x)), (d,n)))
+        KNN = kmeans(X_C, moe.k)
+        a = assignments(KNN)
+        n_added = length(x_new)
+        for i = 1:n_added
+            pos = a[n-n_added+i]
+            add_point!(moe.local_surr[i],moe.x[n-n_added+i],moe.y[n-n_added+i])
+        end
+    end
+    nothing
+end
diff --git a/src/Surrogates.jl b/src/Surrogates.jl
index 8437f497c..6697693e4 100644
--- a/src/Surrogates.jl
+++ b/src/Surrogates.jl
@@ -16,7 +16,11 @@ include("SthenoKriging.jl")
 include("RandomForestSurrogate.jl")
 include("NeuralSurrogate.jl")
 include("Wendland.jl")
+include("MOE.jl")
 
+current_surrogates_MOE = ["Kriging","LinearSurrogate","LobacheskySurrogate","NeuralSurrogate",
+                     "RadialBasis","RandomForestSurrogate","SecondOrderPolynomialSurrogate","Wendland"]
+export current_surrogates_MOE
 export AbstractSurrogate, SamplingAlgorithm
 export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
 export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial
@@ -32,5 +36,9 @@ export InverseDistanceSurrogate
 export SecondOrderPolynomialSurrogate
 export SthenoKriging
 export Wendland
+export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure
+export LobacheskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure
+export WendlandStructure
+export MOE
 
 end
diff --git a/test/MOE.jl b/test/MOE.jl
new file mode 100644
index 000000000..b4757d094
--- /dev/null
+++ b/test/MOE.jl
@@ -0,0 +1,43 @@
+using Surrogates
+using Distributions
+#1D MOE
+
+n = 30
+lb = 0.0
+ub = 5.0
+x1 = Surrogates.sample(n,1,Normal(2.5,0.1))
+x2 = Surrogates.sample(n,1,Normal(1.0,0.4))
+x = vcat(x1,x2)
+f = x-> 2*x
+y = f.(x)
+#Standard definition
+my_moe = MOE(x,y,lb,ub)
+val = my_moe(3.0)
+add_point!(my_moe,3.0,6.0)
+add_point!(my_moe,[4.0,5.0],[8.0,10.0])
+
+#Local surrogates redefinition
+my_local_kind = [InverseDistanceStructure(p = 1.0),
+                 RadialBasisStructure(radial_function = cubicRadial, scale_factor=1.0,sparse = false)]
+my_moe = MOE(x,y,lb,ub,k = 2,local_kind = my_local_kind)
+
+
+
+#ND MOE
+n = 30
+lb = [0.0,0.0]
+ub = [5.0,5.0]
+x1 = Surrogates.sample(n,2,Normal(0.0,4.0))
+x2 = Surrogates.sample(n,2,Normal(3.0,5.0))
+x = vcat(x1,x2)
+f = x -> x[1]*x[2]
+y = f.(x)
+my_moe_ND = MOE(x,y,lb,ub)
+val = my_moe_ND((1.0,1.0))
+
+add_point!(my_moe_ND, (1.0,1.0), 1.0)
+
+#Local surr redefinition
+my_locals = [InverseDistanceStructure(p = 1.0),
+             RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0,sparse = false)]
+my_moe_redef = MOE(x,y,lb,ub,k = 2,local_kind = my_local_kind)
diff --git a/test/runtests.jl b/test/runtests.jl
index 47edb5218..41c535a8c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,3 +15,4 @@ using Surrogates
 @testset "AD_Compatibility" begin include("AD_compatibility.jl") end
 @testset "SthenoKriging.jl" begin include("SthenoKriging.jl") end
 @testset "Wendland" begin include("Wendland.jl") end
+@testset "MOE" begin include("MOE.jl") end