From 970d0e8126a0b95b9998be1a22f2524b4eb2f1f0 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlsson <cont-frb@ulund.org>
Date: Thu, 10 Sep 2015 12:52:10 +0200
Subject: [PATCH] Extracted jacobians and other functions from train_RBF_ARX

---
 src/rbfmodels/levenberg_marquardt.jl |   1 +
 src/rbfmodels/trainRBF_ARX.jl        | 144 +++++++++++++++------------
 2 files changed, 84 insertions(+), 61 deletions(-)

diff --git a/src/rbfmodels/levenberg_marquardt.jl b/src/rbfmodels/levenberg_marquardt.jl
index 030d5cb..3f30a3f 100644
--- a/src/rbfmodels/levenberg_marquardt.jl
+++ b/src/rbfmodels/levenberg_marquardt.jl
@@ -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)
 
diff --git a/src/rbfmodels/trainRBF_ARX.jl b/src/rbfmodels/trainRBF_ARX.jl
index b2db6ef..1605083 100644
--- a/src/rbfmodels/trainRBF_ARX.jl
+++ b/src/rbfmodels/trainRBF_ARX.jl
@@ -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
-- 
GitLab