module SystemIdentification
if !isdefined(:DEBUG); DEBUG = false; end
export
Model,LinearModel,NonLinearModel,
Network,
Polynom,PolynomMatrix,
TFdata,
AR,
ARX,
RBFARX,
FitStatistics,
FitResult,
IdData,
# Functions
ar,arx,getARregressor,getARXregressor,find_na,
toeplitz, kalman, kalman_smoother, forward_kalman, PCA, plotmodel

## Fit Methods =================
:LS
:LS_reg
:L1
:LM

## Types =======================
abstract Model
abstract LinearModel <: Model
abstract NonLinearModel <: Model
abstract Network <: NonLinearModel
typealias Polynom{T<:Real} Union{Array{T,1} , T}
typealias PolynomMatrix{T} Union{Array{Polynom{T},1},Polynom{T}, T}



"""
`a::Vector{Float64}`: The polynomial coeffs A(z) (not including the first 1)\n
`na::Int`\n
"""
type AR <: LinearModel
    a::Polynom
    na::Int
end
function Base.show(m::AR)
    print("A(z) = 1")
    for i = 1:length(m.a)
        print(m.a[i], " + z^-$i")
    end
    print("\n")
end

"""
`a::Vector{Float64}`: The polynomial coeffs A(z) (not including the first 1)\n
`b::VecOrMat{Float64}`: The polynomial coeffs B(z)\n
`na::Int`\n
`nb::Vector{Int}`\n
"""
type ARX <: LinearModel
    a::Polynom
    b::PolynomMatrix
    na::Int
    nb::Polynom{Int}
end

type TFdata <: LinearModel
    P
    F
end

type RBFARX <: Network
    na::Int
    n_centers::Int
    n_state::Int
    outputnet::Bool
    normalized::Bool
end

type FitStatistics
    RMS::Float64
    FIT::Float64
    AIC::Float64
    FitStatistics(y::Vector{Float64},yh::Vector{Float64},d) = new(rms(y-yh),fit(y,yh),aic(y-yh,d))
end

type FitResult
    error::VecOrMat{Float64}
    fit::FitStatistics
    method::Symbol
    λ::Float64
    function FitResult(y,yh,d::Int,method::Symbol, λ = 0)
        error = y-yh
        new(error, FitStatistics(y,yh, d),method, λ)
    end
end

type IdData
    y::VecOrMat{Float64}
    u::VecOrMat{Float64}
    Ts::Float64
    IdData(y::VecOrMat{Float64}) = new(y,[],0)
    IdData(y::VecOrMat{Float64} ,Ts::Float64) = new(y,[],Ts)
    IdData(y::VecOrMat{Float64}, u::VecOrMat{Float64}) = new(y,u,0)
end

## Helper functions
rms(x) = sqrt(mean(x.^2))
sse(x) = sum(x.^2)
fit(y,yh) = 100 * (1-rms(y-yh)./rms(y-mean(y)));
aic(x,d) = log(sse(x)) + 2d/length(x)

##
Base.show(fit::FitStatistics) = println("Fit RMS:$(fit.RMS), FIT:$(fit.FIT), AIC:$(fit.AIC)")


include("armax.jl")
include("kalman.jl")
include("PCA.jl")
# include("kernelPCA.jl")
include("toeplitz.jl")

end