Commit 970d0e81 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

Extracted jacobians and other functions from train_RBF_ARX

parent b0c08636
......@@ -17,6 +17,7 @@ using Optim
# returns: x, J
# x - least squares solution for x
# J - estimate of the Jacobian of f at x
# This code is copied from the package Optim.jl and modified to include bounds checking and a timeout feature
"""
function levenberg_marquardt(f::Function, g::Function, x0; tolX=1e-8, tolG=1e-12, maxIter=100, lambda=100.0, show_trace=false, timeout=Inf, n_state = 1)
......
......@@ -32,67 +32,11 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
n_points = length(y)
na = isa(A,Matrix) ? size(A,2) : 1
function saturatePrecision(x,n_state)
for i = 1:2n_state:length(x)
range = i+n_state:i+2n_state-1
x[range] = abs(x[range])
end
return x
end
function fitlinear(V)
try
DEBUG && assert(isa(V,Matrix))
DEBUG && assert(isa(y,Vector))
DEBUG && assert(!any(!isfinite(V)))
DEBUG && assert(!any(!isfinite(y)))
return V\y
catch ex
@show reshape(Znl.x,2n_state,n_centers)
display(ex)
error("Linear fitting failed")
end
end
function jacobian(Znl, Ψ, w)
J = Array(Float64,(n_points,length(Znl.x)))
ii = 1
for (k,Zi) in enumerate(Znl)
μ = Zi[1:n_state] # slice?
γ = Zi[n_state+1:end]
i1 = ii-1
for l = 1:n_points
Ψw = 1.0
if outputnet
for i = 1:na
for j = 1:n_state
ind = j + n_state*(k-1) + n_state*(n_centers+1)*(i-1)
Ψw += V[l,ind]*w[ind]
end
end
else
for i = 1:na
ind = k + (n_centers+1)*(i-1)
Ψw += V[l,ind]*w[ind]
end
end
for p = 1:n_state
x_μ = state[l,p]-μ[p]
J[l,i1+p] = 2*Ψw*x_μ*γ[p]
J[l,i1+n_state+p] = (-Ψw)*x_μ^2
end
end
ii += 2n_state
end
return J
end
function predictionerror(z)
znl = RbfNonlinearParameters(z,n_state,n_centers)
psi = getΨ(Ψ, znl, state, n_points, n_state, normalized)
getLinearRegressor(V,psi,A,state,outputnet,na,n_state,n_centers,n_points)
zl = fitlinear(V);
zl = fitlinear(V,y);
prediction = V*zl
error = prediction-y
return error
......@@ -135,7 +79,7 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
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")
Zl = fitlinear(V,y); debug("fitlinear")
@ddshow sum(!isfinite(Zl))
prediction = V*Zl
error = y - prediction
......@@ -148,8 +92,8 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
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)
w = fitlinear(V,y)
return outputnet ? jacobian_outputnet(znl,Ψ, w, V) : jacobian_no_outputnet(znl,Ψ, w, V)
end
f(z) = predictionerror(z)
X0 = deepcopy(Znl.x)
......@@ -205,7 +149,7 @@ function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=fal
# Test ===============================================
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")
Zl = fitlinear(V,y); debug("fitlinear")
prediction = V*Zl
error = y - prediction
if PYPLOT || WINSTON
......@@ -390,3 +334,81 @@ function getcentersKmeans(state, nc::Int, f::Function, n_state::Int)
ind = indmin(errorvec)
return RbfNonlinearParameters(params[:,ind],n_state, nc)
end
function jacobian_outputnet(Znl, Ψ, w, V)
n_points = size(Ψ,1)
n_state = Znl.n_state
n_centers = Znl.n_centers
J = Array(Float64,(n_points,length(Znl.x)))
ii = 1
for (k,Zi) in enumerate(Znl)
μ = Zi[1:n_state] # slice?
γ = Zi[n_state+1:end]
i1 = ii-1
for l = 1:n_points
Ψw = 1.0
for i = 1:na
for j = 1:n_state
ind = j + n_state*(k-1) + n_state*(n_centers+1)*(i-1)
Ψw += V[l,ind]*w[ind]
end
end
for p = 1:n_state
x_μ = state[l,p]-μ[p]
J[l,i1+p] = 2*Ψw*x_μ*γ[p]
J[l,i1+n_state+p] = (-Ψw)*x_μ^2
end
end
ii += 2n_state
end
return J
end
function jacobian_no_outputnet(Znl, Ψ, w,v)
n_points = size(Ψ,1)
n_state = Znl.n_state
n_centers = Znl.n_centers
J = Array(Float64,(n_points,length(Znl.x)))
ii = 1
for (k,Zi) in enumerate(Znl)
μ = Zi[1:n_state] # slice?
γ = Zi[n_state+1:end]
i1 = ii-1
for l = 1:n_points
Ψw = 1.0
for i = 1:na
ind = k + (n_centers+1)*(i-1)
Ψw += V[l,ind]*w[ind]
end
for p = 1:n_state
x_μ = state[l,p]-μ[p]
J[l,i1+p] = 2*Ψw*x_μ*γ[p]
J[l,i1+n_state+p] = (-Ψw)*x_μ^2
end
end
ii += 2n_state
end
return J
end
function saturatePrecision(x,n_state)
for i = 1:2n_state:length(x)
range = i+n_state:i+2n_state-1
x[range] = abs(x[range])
end
return x
end
function fitlinear(V,y)
try
DEBUG && assert(isa(V,Matrix))
DEBUG && assert(isa(y,Vector))
DEBUG && assert(!any(!isfinite(V)))
DEBUG && assert(!any(!isfinite(y)))
return V\y
catch ex
display(ex)
error("Linear fitting failed")
end
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment