diff --git a/src/cuckooSearch.jl b/src/cuckooSearch.jl
index 0fed3056a9fca939415471597755ac0a7b3b60f2..23751b5adc08a385f9ac62a9e33a3838e8c841eb 100644
--- a/src/cuckooSearch.jl
+++ b/src/cuckooSearch.jl
@@ -15,7 +15,8 @@ Based on implementation by
}
http://www.mathworks.com/matlabcentral/fileexchange/29809-cuckoo-search--cs--algorithm
"""
-function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)
+function cuckoo_search(f,X0, Lb,Ub;n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)
+ X00 = deepcopy(X0)
nd=size(X0,1);
X0t = X0'
Lb = Lb'
@@ -29,13 +30,18 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2
# Random initial solutions
nest = zeros(n,nd)
- nest[1,:] = X0
+ nest[1,:] = X0t
for i=2:n
nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb));
+ DEBUG && @assert !any(nest[i,:] .> Ub)
+ DEBUG && @assert !any(nest[i,:] .< Lb)
end
# Get the current best
+
fitness=10^20*ones(n,1);
fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness);
+ DEBUG && println("f(X0) = $(f(X00)), f(bestnest) = $(fmin)")
+ DEBUG && @assert X00 == X0
N_iter=0;
t0 = time()
## Starting iterations
@@ -72,6 +78,7 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2
## Post-optimization processing
## Display all the nests
println("Total number of iterations=",N_iter);
+ println("f(bestnest) = $(fmin)")
squeeze(bestnest',2),fmin
end
@@ -106,8 +113,10 @@ function get_cuckoos(nest,best,Lb,Ub)
s=s+stepsize.*randn(size(s));
# Apply simple bounds/limits
nest[j,:]=simplebounds(s,Lb,Ub);
+ DEBUG && @assert !any(nest[j,:] .> Ub)
+ DEBUG && @assert !any(nest[j,:] .< Lb)
end
- nest
+ return nest
end
## Find the current best nest
function get_best_nest(f,nest,newnest,fitness)
@@ -123,7 +132,7 @@ function get_best_nest(f,nest,newnest,fitness)
# Find the current best
(fmin,K) = findmin(fitness) ;
best=nest[K,:];
- fmin,best,nest,fitness
+ return fmin,best,nest,fitness
end
@@ -140,19 +149,21 @@ function empty_nests(nest,Lb,Ub,pa)
## New solution by biased/selective random walks
stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]);
new_nest=nest+stepsize.*K;
+ for j = 1:size(nest,1)
+ new_nest[j,:]=simplebounds(new_nest[j,:],Lb,Ub);
+ end
+ return new_nest
end
# Application of simple constraints
function simplebounds(s,Lb,Ub)
# Apply the lower bound
- ns_tmp=s;
- I=ns_tmp.<Lb;
- ns_tmp[I]=Lb[I];
+ I = s.<Lb;
+ s[I] = Lb[I];
# Apply the upper bounds
- J=ns_tmp.>Ub;
- ns_tmp[J]=Ub[J];
- # Update this new move
- s=ns_tmp;
+ J = s.>Ub;
+ s[J] = Ub[J];
+ return s
end
# # ## You can replace the following by your own functions
diff --git a/src/rbfmodels/trainRBF_ARX.jl b/src/rbfmodels/trainRBF_ARX.jl
index bba0cb4117f8cab5a6dec73335f2ddd3ee43481f..429c2e800ec021dece8b4ddb22ac36c67ea0164f 100644
--- a/src/rbfmodels/trainRBF_ARX.jl
+++ b/src/rbfmodels/trainRBF_ARX.jl
@@ -1,7 +1,8 @@
using Devectorize
using Clustering
-using Debug
-
+# using Debug
+include("levenberg_marquardt.jl")
+include("../cuckooSearch.jl")
type RbfNonlinearParameters
x::Vector{Float64}
@@ -9,14 +10,6 @@ type RbfNonlinearParameters
n_centers::Integer
end
-type RbfLinearParameters
- x::Vector{Float64}
-end
-
-type RbfParameters
- n::RbfNonlinearParameters
- l::RbfLinearParameters
-end
@@ -51,6 +44,7 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
params[si:si+n_state-1,iter] = clusterresult.centers[:,i]
C = cov(state[clusterresult.assignments .== i,:])
params[si+n_state:si+2n_state-1,iter] = diag(inv(C))
+ @assert !any(diag(inv(C)) .< 0)
end
errorvec[iter] = rms(predictionerror(params[:,iter]))
end
@@ -76,57 +70,12 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
newplot(X[:,1],X[:,2],"o"); title("Centers")
end
- function getΨ(Znl)
- RBF(x, Znl::VecOrMat,n_state::Integer) = exp(-(((x-Znl[1:n_state]).^2.*Znl[n_state+1:end])[1]))
- if normalized
- rowsum = ones(n_points)
- for (j,Zi) in enumerate(Znl)
- for i = n_state+1:2n_state
- Zi[i] = Zi[i] <= 0 ? 1.0 : Zi[i] # Reset to 1 if precision became negative
- end
- for i = 1:n_points
- statei = squeeze(state[i,:]',2)
- a = RBF(statei, Zi, n_state)
- Ψ[i,j] = a
- rowsum[i] += a
- end
- end
- for i = 1:n_points
- if rowsum[i] <= 1e-10
- continue
- end
- @devec Ψ[i,:] ./= rowsum[i]
- end
- else # Not normalized
- for (j,Zi) in enumerate(Znl)
- for i = n_state+1:2n_state
- Zi[i] = Zi[i] <= 0 ? 1.0 : Zi[i] # Reset to 1 if precision became negative
- end
- for i = 1:n_points
- statei = squeeze(state[i,:]',2)
- # statei = slice(state,i,:)
- Ψ[i,j] = RBF(statei, Zi, n_state)
- if DEBUG && !isfinite(Ψ[i,j])
- @show i,j,statei, Zi, n_state, Ψ[i,j]
- @show (statei-Zi[1:n_state]).^2
- @show Zi[n_state+1:end]
- # @show exp(-(((statei-Zi[1:n_state]).^2.*Zi[n_state+1:end])[1]))
- error("Stopping")
- end
- end
- end
- end
- if DEBUG && sum(!isfinite(Ψ)) > 0
- @show sum(!isfinite(Ψ))
- end
- return Ψ
- end
function fitlinear(V)
try
- assert(isa(V,Matrix))
- assert(isa(y,Vector))
+ DEBUG && assert(isa(V,Matrix))
+ DEBUG && assert(isa(y,Vector))
DEBUG && assert(!any(!isfinite(V)))
DEBUG && assert(!any(!isfinite(y)))
return V\y
@@ -170,42 +119,14 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
return J
end
- function getLinearRegressor(Ψ)
- if outputnet
- ii = 1
- for i = 1:na
- for k = 1:n_centers+1
- for j = 1:n_state
- for l = 1:n_points
- V[l,ii] = Ψ[l,k]*A[l,i]*state[l,j]
- end
- ii = ii+1
- end
- end
- end
- else
- ii = 1
- for i = 1:na
- for k = 1:n_centers+1
- for l = 1:n_points
- V[l,ii] = Ψ[l,k]*A[l,i]
- end
- ii = ii+1
- end
- end
- end
- if DEBUG && sum(!isfinite(V)) > 0
- @show sum(!isfinite(V))
- end
- return V
- end
+
function predictionerror(z)
znl = RbfNonlinearParameters(z,n_state,n_centers)
- getΨ(znl);
- getLinearRegressor(Ψ);
- zl = fitlinear(V);
- prediction = V*zl
+ psi = getΨ(deepcopy(Ψ), znl, state, n_points, n_state, normalized)
+ v = getLinearRegressor(deepcopy(V),psi,A,state,outputnet,na,n_state,n_centers,n_points)
+ zl = fitlinear(v);
+ prediction = v*zl
error = prediction-y
return error
end
@@ -243,9 +164,9 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
Znl = getcentersKmeans(state,nc); debug("gotcentersKmeans")
end
@ddshow Znl
- getΨ(Znl); debug("Got Ψ")
+ getΨ(Ψ, Znl, state, n_points, n_state, normalized); debug("Got Ψ")
@ddshow sum(!isfinite(Ψ))
- getLinearRegressor(Ψ); debug("Got linear regressor V")
+ getLinearRegressor(V,Ψ,A,state,outputnet,na,n_state,n_centers,n_points); debug("Got linear regressor V")
@ddshow size(V)
@ddshow sum(!isfinite(V))
Zl = fitlinear(V); debug("fitlinear")
@@ -254,55 +175,80 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
error = y - prediction
errors = zeros(liniters+1)
+ Lb = zeros(Znl.x)
+ Ub = zeros(Znl.x)
+ mas = maximum(state,1)'
+ mis = minimum(state,1)'
+ for i = 1:2n_state:n_centers*2n_state
+ Lb[i:i+n_state-1] = mis
+ Ub[i:i+n_state-1] = mas
+ Lb[i+n_state:i+2n_state-1] = 0.000001
+ Ub[i+n_state:i+2n_state-1] = 10*Znl.x[n_state+1:2n_state]
+ end
+ @show Lb
+ @show Ub
+
# ============= Main loop ================================================
debug("Calculating initial error")
- errors[1] = sse(predictionerror(Znl.x))
+ errors[1] = rms(predictionerror(Znl.x))
println("Training RBF_ARX Centers: $(Znl.n_centers), Nonlinear parameters: $(length(Znl.x)), Linear parameters: $(length(Zl))")
+ function g(z)
+ znl = RbfNonlinearParameters(z,n_state,n_centers)
+ w = fitlinear(V)
+ jacobian(znl,Ψ, w)
+ end
+ f(z) = predictionerror(z)
+ X0 = deepcopy(Znl.x)
for i = 1:liniters
- function g(z)
- znl = RbfNonlinearParameters(z,n_state,n_centers)
- w = fitlinear(V)
- jacobian(znl,Ψ, w)
- end
- f(z) = predictionerror(z)
-
- X0 = Znl.x
-
if i%2 == 1 || !cuckoosearch
@time res = levenberg_marquardt(f, g, X0,
maxIter = nonliniters,
tolG = 1e-7,
tolX = 1e-10,
show_trace=true,
- timeout = 60)
- Znl = RbfNonlinearParameters(saturatePrecision(res.minimum,n_state),n_state,n_centers)
+ timeout = 60,
+ n_state = n_state)
+ X0 = deepcopy(res.minimum)
+ assert(X0 == res.minimum)
+ @show rms(f(X0))
+ if DEBUG
+ _V = deepcopy(V)
+ _Ψ = deepcopy(Ψ)
+ end
+ @show rms(f(res.minimum))
+ if DEBUG
+ @show res.minimum == X0
+ @show _V == V
+ @show _Ψ == Ψ
+ end
+ assert(X0 == res.minimum)
+ # Znl = RbfNonlinearParameters(saturatePrecision(copy(res.minimum),n_state),n_state,n_centers)
+ Znl = RbfNonlinearParameters(deepcopy(res.minimum),n_state,n_centers)
errors[i+1] = res.f_minimum
# show(res.trace)
else
display("Using cuckoo search to escape local minimum")
- @time (bestnest,fmin) = cuckoo_search(x -> sum(f(x).^2),X0;
- n=30,
+ @time (bestnest,fmin) = cuckoo_search(x -> rms(f(x)),X0, Lb, Ub;
+ n=50,
pa=0.25,
Tol=1.0e-5,
max_iter = i < liniters-1 ? 80 : 200,
timeout = 120)
debug("cuckoo_search done")
- Znl = RbfNonlinearParameters(bestnest,n_state,n_centers)
+ X0 = deepcopy(bestnest)
+ @ddshow rms(f(X0))
+ Znl = RbfNonlinearParameters(deepcopy(bestnest),n_state,n_centers)
errors[i+1] = fmin
end
if abs(errors[i+1]-errors[i]) < 1e-10
display("No significant change in function value")
break
end
- getΨ(Znl)
- getLinearRegressor(Ψ)
- fitlinear(V)
- # Znl.x = res.minimum
end
# Test ===============================================
- getΨ(Znl)
- getLinearRegressor(Ψ)
+ getΨ(Ψ, Znl, state, n_points, n_state, normalized)
+ getLinearRegressor(V,Ψ,A,state,outputnet,na,n_state,n_centers,n_points)
Zl = fitlinear(V); debug("fitlinear")
prediction = V*Zl
error = y - prediction
@@ -361,4 +307,85 @@ end
+function getΨ(Ψ, Znl, state, n_points, n_state, normalized::Bool)
+ Ψ = normalized ? getΨnormalized(Ψ, Znl, state, n_points, n_state) : getΨnonnormalized(Ψ, Znl, state, n_points, n_state)
+ if DEBUG && sum(!isfinite(Ψ)) > 0
+ @show sum(!isfinite(Ψ))
+ end
+ return Ψ
+end
+
+function getΨnormalized(Ψ, Znl, state, n_points, n_state)
+ RBF(x, Znl::VecOrMat,n_state::Integer) = exp(-(((x-Znl[1:n_state]).^2.*Znl[n_state+1:end])[1]))
+ rowsum = ones(n_points)
+ for (j,Zin) in enumerate(Znl)
+ Zi = deepcopy(Zin)
+ # for i = n_state+1:2n_state
+ # Zi[i] = Zi[i] <= 0 ? 0.01 : Zi[i] # Reset to 1 if precision became negative
+ # end
+ for i = 1:n_points
+ statei = squeeze(state[i,:]',2)
+ a = RBF(statei, Zi, n_state)
+ Ψ[i,j] = a
+ rowsum[i] += a
+ end
+ end
+ for i = 1:n_points
+ if rowsum[i] <= 1e-10
+ continue
+ end
+ @devec Ψ[i,:] ./= rowsum[i]
+ end
+ return Ψ
+end
+
+function getΨnonnormalized(Ψ, Znl, state, n_points, n_state)
+ RBF(x, Znl::VecOrMat,n_state::Integer) = exp(-(((x-Znl[1:n_state]).^2.*Znl[n_state+1:end])[1]))
+ for (j,Zin) in enumerate(Znl)
+ Zi = deepcopy(Zin)
+ for i = 1:n_points
+ statei = squeeze(state[i,:]',2)
+ # statei = slice(state,i,:)
+ Ψ[i,j] = RBF(statei, Zi, n_state)
+ if DEBUG && !isfinite(Ψ[i,j])
+ @show i,j,statei, Zi, n_state, Ψ[i,j]
+ @show (statei-Zi[1:n_state]).^2
+ @show Zi[n_state+1:end]
+ # @show exp(-(((statei-Zi[1:n_state]).^2.*Zi[n_state+1:end])[1]))
+ error("Stopping")
+ end
+ end
+ end
+ return Ψ
+end
+
+function getLinearRegressor(V,Ψ,A,state,outputnet,na,n_state,n_centers,n_points)
+ if outputnet
+ ii = 1
+ for i = 1:na
+ for k = 1:n_centers+1
+ for j = 1:n_state
+ for l = 1:n_points
+ V[l,ii] = Ψ[l,k]*A[l,i]*state[l,j]
+ end
+ ii = ii+1
+ end
+ end
+ end
+ else
+ ii = 1
+ for i = 1:na
+ for k = 1:n_centers+1
+ for l = 1:n_points
+ V[l,ii] = Ψ[l,k]*A[l,i]
+ end
+ ii = ii+1
+ end
+ end
+ end
+ if DEBUG && sum(!isfinite(V)) > 0
+ @show sum(!isfinite(V))
+ end
+ return V
+end