Merge pull request #58 from acfr/57-polar-parameterisation-for-ren-is…

…-incorrect

Fixed polar parameterisation bug

src/Base/direct_params.jl

@@ -80,8 +80,6 @@ function DirectParams{T}(
         B2  = glorot_normal(nx, nu; T=T, rng=rng)
         D12 = glorot_normal(nv, nu; T=T, rng=rng)
 
-        ρ = zeros(1)
-
         # Make orthogonal X
         X = glorot_normal(2nx + nv, 2nx + nv; T=T, rng=rng)
         X = Matrix(qr(X).Q)
@@ -107,13 +105,15 @@ function DirectParams{T}(
                  -C1 H22 B1'
                  F  B1  P] + ϵ * I
 
-        ρ = zeros(T, 1)
         X = Matrix{T}(cholesky(Htild).U) # H = X'*X
 
     else
         error("Undefined initialisation method ", init)
     end
 
+    # Polar parameter
+    ρ = [norm(X)]
+
     # Free parameter for E
     Y1 = glorot_normal(nx, nx; T=T, rng=rng)
 

src/ParameterTypes/contracting_ren.jl

@@ -169,9 +169,10 @@ end
 function direct_to_explicit(ps::ContractingRENParams{T}, return_h::Bool=false) where T
 
     ϵ = ps.direct.ϵ
-    ρ = ps.direct.ρ
+    ρ = ps.direct.ρ[1]
     X = ps.direct.X
-    H = ps.direct.polar_param ? exp(ρ[1])*(X'*X + ϵ*I) / norm(X)^2 : X'*X + ϵ*I
+    polar_param = ps.direct.polar_param
+    H = x_to_h(X, ϵ, polar_param, ρ)
 
     !return_h && (return hmatrix_to_explicit(ps, H))
     return H

src/ParameterTypes/general_ren.jl

@@ -112,8 +112,9 @@ function direct_to_explicit(ps::GeneralRENParams{T}, return_h=false) where T
 
     # Implicit parameters
     ϵ = ps.direct.ϵ
-    ρ = ps.direct.ρ
+    ρ = ps.direct.ρ[1]
     X = ps.direct.X
+    polar_param = ps.direct.polar_param
 
     X3 = ps.direct.X3
     Y3 = ps.direct.Y3
@@ -150,11 +151,7 @@ function direct_to_explicit(ps::GeneralRENParams{T}, return_h=false) where T
     Γ1 = [C2'; D21'; zeros(nx, ny)] * Q * [C2 D21 zeros(ny, nx)]
     Γ2 = [C2_imp'; D21_imp'; B2_imp] * (𝑅 \ [C2_imp D21_imp B2_imp'])
 
-    if ps.direct.polar_param 
-        H = exp(ρ[1])*(X'*X + ϵ*I) / norm(X)^2 + Γ2 - Γ1
-    else
-        H = X'*X + ϵ*I + Γ2 - Γ1
-    end
+    H = x_to_h(X, ϵ, polar_param, ρ) + Γ2 - Γ1
 
     # Get explicit parameterisation
     !return_h && (return hmatrix_to_explicit(ps, H, D22))

src/ParameterTypes/lipschitz_ren.jl

@@ -90,8 +90,9 @@ function direct_to_explicit(ps::LipschitzRENParams{T}, return_h=false) where T
 
     # Implicit parameters
     ϵ = ps.direct.ϵ
-    ρ = ps.direct.ρ
+    ρ = ps.direct.ρ[1]
     X = ps.direct.X
+    polar_param = ps.direct.polar_param
 
     X3 = ps.direct.X3
     Y3 = ps.direct.Y3
@@ -123,11 +124,7 @@ function direct_to_explicit(ps::LipschitzRENParams{T}, return_h=false) where T
     Γ1 = [C2'; D21'; zeros(nx, ny)] * [C2 D21 zeros(ny, nx)] * (-1/γ)
     Γ2 = [C2_imp'; D21_imp'; B2_imp] * (𝑅 \ [C2_imp D21_imp B2_imp'])
 
-    if ps.direct.polar_param 
-        H = exp(ρ[1])*(X'*X + ϵ*I) / norm(X)^2 + Γ2 - Γ1
-    else
-        H = X'*X + ϵ*I + Γ2 - Γ1
-    end
+    H = x_to_h(X, ϵ, polar_param, ρ) + Γ2 - Γ1
 
     # Get explicit parameterisation
     !return_h && (return hmatrix_to_explicit(ps, H, D22))

src/ParameterTypes/passive_ren.jl

@@ -96,8 +96,9 @@ function direct_to_explicit(ps::PassiveRENParams{T}, return_h=false) where T
 
     # Implicit parameters
     ϵ = ps.direct.ϵ
-    ρ = ps.direct.ρ
+    ρ = ps.direct.ρ[1]
     X = ps.direct.X
+    polar_param = ps.direct.polar_param
 
     X3 = ps.direct.X3
     Y3 = ps.direct.Y3
@@ -121,12 +122,7 @@ function direct_to_explicit(ps::PassiveRENParams{T}, return_h=false) where T
 
     Γ2 = [C2'; D21_imp'; B2_imp] * (𝑅 \ [C2 D21_imp B2_imp'])
 
-    if ps.direct.polar_param 
-        # See Eqns 29 of TAC paper 
-        H = exp(ρ[1])*X'*X / norm(X)^2 + Γ2
-    else
-        H = X'*X + ϵ*I + Γ2
-    end
+    H = x_to_h(X, ϵ, polar_param, ρ) + Γ2
 
     # Get explicit parameterisation
     !return_h && (return hmatrix_to_explicit(ps, H, D22))

src/ParameterTypes/utils.jl

@@ -75,4 +75,8 @@ function hmatrix_to_explicit(ps::AbstractRENParams, H::AbstractMatrix{T}, D22::A
 
     return ExplicitParams{T}(A, B1, B2, C1, C2, D11, D12, D21, D22, bx, bv, by)
 
+end
+
+function x_to_h(X::AbstractMatrix{T}, ϵ::T, polar_param::Bool, ρ::T) where T
+    polar_param ? (ρ^2)*(X'*X) / norm(X)^2 + ϵ*I : X'*X + ϵ*I
 end