Fixed bug when creating periodic subLHC

MrUrq · Mar 30, 2020 · 0e1f7a8 · 0e1f7a8 · MrUrq · Apr 3, 2020
1 parent 90d758f
commit 0e1f7a8
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 16 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "LatinHypercubeSampling"
 uuid = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
 authors = ["Magnus Urquhart <[email protected]>"]
-version = "1.6.3"
+version = "1.6.4"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

diff --git a/src/AudzeEglaisObjective.jl b/src/AudzeEglaisObjective.jl
@@ -1,4 +1,4 @@
-@inline function _AudzeEglaisDist(LHC,periodic_ae,ae_power)
+@inline function _AudzeEglaisDist(LHC,periodic_ae,ae_power,periodic_n)
     n,d = size(LHC)
     dist = 0.0
 
@@ -9,7 +9,7 @@
             for k = 1:d
                 if periodic_ae
                     @inbounds dist_comp = abs(LHC[i,k]-LHC[j,k])
-                    dist_comp = min(dist_comp,n-dist_comp)
+                    dist_comp = min(dist_comp,periodic_n-dist_comp)
                 else
                     @inbounds dist_comp = LHC[i,k]-LHC[j,k]
                 end
@@ -28,13 +28,13 @@
     return output
 end
 
-function _AudzeEglaisObjective(dim::Continuous,LHC,periodic_ae,ae_power)
-    output = _AudzeEglaisDist(LHC,periodic_ae,ae_power)    
+function _AudzeEglaisObjective(dim::Continuous,LHC,periodic_ae,ae_power,periodic_n)
+    output = _AudzeEglaisDist(LHC,periodic_ae,ae_power,periodic_n)    
     return output
 end
 
-function _AudzeEglaisObjective(dim::Categorical,LHC,periodic_ae,ae_power)
-    output = _AudzeEglaisDist(LHC,periodic_ae,ae_power)     
+function _AudzeEglaisObjective(dim::Categorical,LHC,periodic_ae,ae_power,periodic_n)
+    output = _AudzeEglaisDist(LHC,periodic_ae,ae_power,periodic_n)     
     output == Inf ? 0 : output    
 end
 
@@ -45,20 +45,21 @@ distance as the objective function. Note this is the inverse of the typical
 Audze-Eglais distance which normally is minimized.
 """
 function AudzeEglaisObjective(LHC::T; dims::Array{V,1} =[Continuous() for i in 1:size(LHC,2)],
-                                            interSampleWeight::Float64=1.0,periodic_ae::Bool=false,ae_power::Union{Int,Float64}=2
+                                            interSampleWeight::Float64=1.0,periodic_ae::Bool=false,ae_power::Union{Int,Float64}=2,
+                                            periodic_n::Int = size(LHC,1)
                                             ) where T <: AbstractArray where V <: LHCDimension
 
     out = 0.0
 
     #Compute the objective function among all points
-    out += _AudzeEglaisObjective(Continuous(),LHC,periodic_ae,ae_power)*interSampleWeight
+    out += _AudzeEglaisObjective(Continuous(),LHC,periodic_ae,ae_power,periodic_n)*interSampleWeight
 
     #Compute the objective function within each categorical dimension
     categoricalDimInds = findall(x->typeof(x)==Categorical,dims)
     for i in categoricalDimInds
         for j = 1:dims[i].levels
             subLHC = @view LHC[LHC[:,i] .== j,:] 
-            out += _AudzeEglaisObjective(dims[i],subLHC,periodic_ae,ae_power)*dims[i].weight
+            out += _AudzeEglaisObjective(dims[i],subLHC,periodic_ae,ae_power,periodic_n)*dims[i].weight
         end
     end
 
@@ -67,8 +68,9 @@ end
 
 # Remove depwarning in release 2.x.x
 function AudzeEglaisObjective!(dist,LHC::T; dims::Array{V,1} =[Continuous() for i in 1:size(LHC,2)],
-                                            interSampleWeight::Float64=1.0,periodic_ae::Bool=false,ae_power::Union{Int,Float64}=2
+                                            interSampleWeight::Float64=1.0,periodic_ae::Bool=false,ae_power::Union{Int,Float64}=2,
+                                            periodic_n::Int = size(LHC,1)
                                             ) where T <: AbstractArray where V <: LHCDimension
     @warn "AudzeEglaisObjective!(dist,LHC) is deprecated and does not differ from AudzeEglaisObjective(LHC)"
-    AudzeEglaisObjective(LHC; dims = dims, interSampleWeight = interSampleWeight, periodic_ae=periodic_ae, ae_power=ae_power)
+    AudzeEglaisObjective(LHC; dims = dims, interSampleWeight = interSampleWeight, periodic_ae=periodic_ae, ae_power=ae_power, periodic_n=periodic_n)
 end
diff --git a/src/LatinHypercubeSampling.jl b/src/LatinHypercubeSampling.jl
@@ -266,12 +266,20 @@ function LHCoptim!(X::Array{Int,2},gens;    popsize::Int=100,
 end
 
 """
-    function subLHCoptim(X,n::Int,gens;popsize::Int=100,ntour::Int=2,ptour::Float64=0.8,periodic_ae::Bool=false,ae_power::Union{Int,Float64}=2)
+    function subLHCoptim(X,n::Int,gens; popsize::Int=100,
+                                        ntour::Int=2,
+                                        ptour::Float64=0.8,
+                                        periodic_ae::Bool=false,
+                                        ae_power::Union{Int,Float64}=2)
 Produce an optimized Latin Hyper Cube with `n` sample points from a subset of points in
 `X`. Optimization is run for `gens` generations. Returns a tuple of the sample plan and
 the optimization fitness history.
 """
-function subLHCoptim(X,n::Int,gens;popsize::Int=100,ntour::Int=2,ptour::Float64=0.8,periodic_ae::Bool=false,ae_power::Union{Int,Float64}=2)
+function subLHCoptim(X,n::Int,gens; popsize::Int=100,
+                                    ntour::Int=2,
+                                    ptour::Float64=0.8,
+                                    periodic_ae::Bool=false,
+                                    ae_power::Union{Int,Float64}=2)
 
     #preallocate memory
     nLarge, d = size(X)
@@ -298,7 +306,10 @@ function subLHCoptim(X,n::Int,gens;popsize::Int=100,ntour::Int=2,ptour::Float64=
     for i = 1:popsize+1
         subInds = sample(1:nLarge, n, replace = false)
         pop[i] = X[subInds,:]
-        fitness[i] = AudzeEglaisObjective(pop[i];periodic_ae=periodic_ae,ae_power=ae_power)
+        fitness[i] = AudzeEglaisObjective(pop[i];
+                                            periodic_ae=periodic_ae,
+                                            ae_power=ae_power,
+                                            periodic_n=nLarge)
     end
 
 
@@ -334,7 +345,10 @@ function subLHCoptim(X,n::Int,gens;popsize::Int=100,ntour::Int=2,ptour::Float64=
 
         #evaluate fitness
         for i = 1:popsize+1
-            fitness[i] = AudzeEglaisObjective(nextpop[i];periodic_ae=periodic_ae,ae_power=ae_power)
+            fitness[i] = AudzeEglaisObjective(nextpop[i];
+                                                periodic_ae=periodic_ae,
+                                                ae_power=ae_power,
+                                                periodic_n=nLarge)
         end
 
         #set the first individual to the best and save the fitness