Skip to content
Snippets Groups Projects
Commit 32b22926 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

Solved the issue of negative precision

parent 59d2e3ac
Branches
No related tags found
1 merge request!1Dev
Show whitespace changes
Inline Side-by-side
src/cuckooSearch.jl
+ 22
− 11
View file @ 32b22926
...@@ -15,7 +15,8 @@ Based on implementation by ...@@ -15,7 +15,8 @@ Based on implementation by
} }
http://www.mathworks.com/matlabcentral/fileexchange/29809-cuckoo-search--cs--algorithm 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); nd=size(X0,1);
X0t = X0' X0t = X0'
Lb = Lb' Lb = Lb'
...@@ -29,13 +30,18 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2 ...@@ -29,13 +30,18 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2
# Random initial solutions # Random initial solutions
nest = zeros(n,nd) nest = zeros(n,nd)
nest[1,:] = X0 nest[1,:] = X0t
for i=2:n for i=2:n
nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb)); nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb));
DEBUG && @assert !any(nest[i,:] .> Ub)
DEBUG && @assert !any(nest[i,:] .< Lb)
end end
# Get the current best # Get the current best
fitness=10^20*ones(n,1); fitness=10^20*ones(n,1);
fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness); 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; N_iter=0;
t0 = time() t0 = time()
## Starting iterations ## Starting iterations
...@@ -72,6 +78,7 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2 ...@@ -72,6 +78,7 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2
## Post-optimization processing ## Post-optimization processing
## Display all the nests ## Display all the nests
println("Total number of iterations=",N_iter); println("Total number of iterations=",N_iter);
println("f(bestnest) = $(fmin)")
squeeze(bestnest',2),fmin squeeze(bestnest',2),fmin
end end
...@@ -106,8 +113,10 @@ function get_cuckoos(nest,best,Lb,Ub) ...@@ -106,8 +113,10 @@ function get_cuckoos(nest,best,Lb,Ub)
s=s+stepsize.*randn(size(s)); s=s+stepsize.*randn(size(s));
# Apply simple bounds/limits # Apply simple bounds/limits
nest[j,:]=simplebounds(s,Lb,Ub); nest[j,:]=simplebounds(s,Lb,Ub);
DEBUG && @assert !any(nest[j,:] .> Ub)
DEBUG && @assert !any(nest[j,:] .< Lb)
end end
nest return nest
end end
## Find the current best nest ## Find the current best nest
function get_best_nest(f,nest,newnest,fitness) function get_best_nest(f,nest,newnest,fitness)
...@@ -123,7 +132,7 @@ function get_best_nest(f,nest,newnest,fitness) ...@@ -123,7 +132,7 @@ function get_best_nest(f,nest,newnest,fitness)
# Find the current best # Find the current best
(fmin,K) = findmin(fitness) ; (fmin,K) = findmin(fitness) ;
best=nest[K,:]; best=nest[K,:];
fmin,best,nest,fitness return fmin,best,nest,fitness
end end
...@@ -140,19 +149,21 @@ function empty_nests(nest,Lb,Ub,pa) ...@@ -140,19 +149,21 @@ function empty_nests(nest,Lb,Ub,pa)
## New solution by biased/selective random walks ## New solution by biased/selective random walks
stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]); stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]);
new_nest=nest+stepsize.*K; 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 end
# Application of simple constraints # Application of simple constraints
function simplebounds(s,Lb,Ub) function simplebounds(s,Lb,Ub)
# Apply the lower bound # Apply the lower bound
ns_tmp=s; I = s.<Lb;
I=ns_tmp.<Lb; s[I] = Lb[I];
ns_tmp[I]=Lb[I];
# Apply the upper bounds # Apply the upper bounds
J=ns_tmp.>Ub; J = s.>Ub;
ns_tmp[J]=Ub[J]; s[J] = Ub[J];
# Update this new move return s
s=ns_tmp;
end end
# # ## You can replace the following by your own functions # # ## You can replace the following by your own functions
......
src/rbfmodels/trainRBF_ARX.jl
+ 140
− 113
View file @ 32b22926
using Devectorize using Devectorize
using Clustering using Clustering
using Debug # using Debug
include("levenberg_marquardt.jl")
include("../cuckooSearch.jl")
type RbfNonlinearParameters type RbfNonlinearParameters
x::Vector{Float64} x::Vector{Float64}
...@@ -9,14 +10,6 @@ type RbfNonlinearParameters ...@@ -9,14 +10,6 @@ type RbfNonlinearParameters
n_centers::Integer n_centers::Integer
end 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 ...@@ -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] params[si:si+n_state-1,iter] = clusterresult.centers[:,i]
C = cov(state[clusterresult.assignments .== i,:]) C = cov(state[clusterresult.assignments .== i,:])
params[si+n_state:si+2n_state-1,iter] = diag(inv(C)) params[si+n_state:si+2n_state-1,iter] = diag(inv(C))
@assert !any(diag(inv(C)) .< 0)
end end
errorvec[iter] = rms(predictionerror(params[:,iter])) errorvec[iter] = rms(predictionerror(params[:,iter]))
end end
...@@ -76,57 +70,12 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal ...@@ -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") newplot(X[:,1],X[:,2],"o"); title("Centers")
end 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) function fitlinear(V)
try try
assert(isa(V,Matrix)) DEBUG && assert(isa(V,Matrix))
assert(isa(y,Vector)) DEBUG && assert(isa(y,Vector))
DEBUG && assert(!any(!isfinite(V))) DEBUG && assert(!any(!isfinite(V)))
DEBUG && assert(!any(!isfinite(y))) DEBUG && assert(!any(!isfinite(y)))
return V\y return V\y
...@@ -170,42 +119,14 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal ...@@ -170,42 +119,14 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
return J return J
end 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) function predictionerror(z)
znl = RbfNonlinearParameters(z,n_state,n_centers) znl = RbfNonlinearParameters(z,n_state,n_centers)
getΨ(znl); psi = getΨ(deepcopy(Ψ), znl, state, n_points, n_state, normalized)
getLinearRegressor(Ψ); v = getLinearRegressor(deepcopy(V),psi,A,state,outputnet,na,n_state,n_centers,n_points)
zl = fitlinear(V); zl = fitlinear(v);
prediction = V*zl prediction = v*zl
error = prediction-y error = prediction-y
return error return error
end end
...@@ -243,9 +164,9 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal ...@@ -243,9 +164,9 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
Znl = getcentersKmeans(state,nc); debug("gotcentersKmeans") Znl = getcentersKmeans(state,nc); debug("gotcentersKmeans")
end end
@ddshow Znl @ddshow Znl
getΨ(Znl); debug("Got Ψ") getΨ(Ψ, Znl, state, n_points, n_state, normalized); debug("Got Ψ")
@ddshow sum(!isfinite(Ψ)) @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 size(V)
@ddshow sum(!isfinite(V)) @ddshow sum(!isfinite(V))
Zl = fitlinear(V); debug("fitlinear") Zl = fitlinear(V); debug("fitlinear")
...@@ -254,55 +175,80 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal ...@@ -254,55 +175,80 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
error = y - prediction error = y - prediction
errors = zeros(liniters+1) 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 ================================================ # ============= Main loop ================================================
debug("Calculating initial error") 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))") println("Training RBF_ARX Centers: $(Znl.n_centers), Nonlinear parameters: $(length(Znl.x)), Linear parameters: $(length(Zl))")
for i = 1:liniters
function g(z) function g(z)
znl = RbfNonlinearParameters(z,n_state,n_centers) znl = RbfNonlinearParameters(z,n_state,n_centers)
w = fitlinear(V) w = fitlinear(V)
jacobian(znl,Ψ, w) jacobian(znl,Ψ, w)
end end
f(z) = predictionerror(z) f(z) = predictionerror(z)
X0 = deepcopy(Znl.x)
X0 = Znl.x for i = 1:liniters
if i%2 == 1 || !cuckoosearch if i%2 == 1 || !cuckoosearch
@time res = levenberg_marquardt(f, g, X0, @time res = levenberg_marquardt(f, g, X0,
maxIter = nonliniters, maxIter = nonliniters,
tolG = 1e-7, tolG = 1e-7,
tolX = 1e-10, tolX = 1e-10,
show_trace=true, show_trace=true,
timeout = 60) timeout = 60,
Znl = RbfNonlinearParameters(saturatePrecision(res.minimum,n_state),n_state,n_centers) 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 errors[i+1] = res.f_minimum
# show(res.trace) # show(res.trace)
else else
display("Using cuckoo search to escape local minimum") display("Using cuckoo search to escape local minimum")
@time (bestnest,fmin) = cuckoo_search(x -> sum(f(x).^2),X0; @time (bestnest,fmin) = cuckoo_search(x -> rms(f(x)),X0, Lb, Ub;
n=30, n=50,
pa=0.25, pa=0.25,
Tol=1.0e-5, Tol=1.0e-5,
max_iter = i < liniters-1 ? 80 : 200, max_iter = i < liniters-1 ? 80 : 200,
timeout = 120) timeout = 120)
debug("cuckoo_search done") 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 errors[i+1] = fmin
end end
if abs(errors[i+1]-errors[i]) < 1e-10 if abs(errors[i+1]-errors[i]) < 1e-10
display("No significant change in function value") display("No significant change in function value")
break break
end end
getΨ(Znl)
getLinearRegressor(Ψ)
fitlinear(V)
# Znl.x = res.minimum
end end
# Test =============================================== # Test ===============================================
getΨ(Znl) getΨ(Ψ, Znl, state, n_points, n_state, normalized)
getLinearRegressor(Ψ) getLinearRegressor(V,Ψ,A,state,outputnet,na,n_state,n_centers,n_points)
Zl = fitlinear(V); debug("fitlinear") Zl = fitlinear(V); debug("fitlinear")
prediction = V*Zl prediction = V*Zl
error = y - prediction error = y - prediction
...@@ -361,4 +307,85 @@ end ...@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment