trainRBF_ARX.jl

using Devectorize
using Clustering
# using Debug
include("levenberg_marquardt.jl")
include("../cuckooSearch.jl")

type RbfNonlinearParameters
    x::Vector{Float64}
    n_state::Integer
    n_centers::Integer
end


# for op = (:+, :*, :\, :/)
#     @eval ($op)(a::RbfNonlinearParameters,b) = ($op)(a.x,b)
#     @eval ($op)(b,a::RbfNonlinearParameters) = ($op)(b,a.x)
# end

Base.display(p::RbfNonlinearParameters) = println("RbfNonlinearParameters: Parameters = $(length(p.x)), n_state(x) = $(p.n_state), n_centers(x) = $(p.n_centers)")
Base.start(p::RbfNonlinearParameters) =  1
Base.done(p::RbfNonlinearParameters, state) = state > p.n_centers
Base.next(p::RbfNonlinearParameters, state) = (p.x[1+(state-1)*2*p.n_state:state*2p.n_state])::VecOrMat, state + 1
"""Train your RBF-ARX network.`trainRBF_ARX(y, na, nc; state = :A, liniters=3,nonliniters=50, normalized=false, initialcenters="equidistant", inputpca=false, outputnet = true, cuckoosearch = false)`\n
The number of centers is equal to `nc` if Kmeans is used to get initial centers, otherwise the number of centers is `nc^n_state`\n
`n_state` is equal to the state dimension, possibly reduced to `inputpca` if so desired.\n
The number of nonlinear parameters is `n_centers × n_state`\n
The number of linear parameters is `outputnet ? n_state × (n_centers+1) × (na)+1) : (na)×(n_centers+1)+1)`"""
function trainRBF_ARX(y, A, state, nc; liniters=3,nonliniters=50, normalized=false, initialcenters="equidistant", inputpca=false, outputnet = true, cuckoosearch = false, cuckooiter=100, timeout=60)
    n_points = length(y)
    na = isa(A,Matrix) ? size(A,2) : 1

    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,y);
        prediction = V*zl
        error = prediction-y
        return error
    end

    # Get initial centers ================================
    Znl::RbfNonlinearParameters
    state,n_state, TT = inputtransform(state, inputpca)

    if lowercase(initialcenters) == "equidistant"
        initialcenters = :equidistant
        n_centers = nc^n_state
    else
        initialcenters = :kmeans
        n_centers = nc
    end

    Ψ = Array(Float64,(n_points,n_centers+1))
    Ψ[:,end] = 1.0
    V = outputnet ? V = Array(Float64,(n_points, n_state* (n_centers+1)* (na)+1)) : V = Array(Float64,(n_points, (na)*(n_centers+1)+1))
    V[:,end] = 1.0
    if initialcenters == :equidistant
        Znl = getcentersEq(state,nc); debug("gotcentersEq")
    else
        Znl = getcentersKmeans(state, nc, predictionerror, n_state); debug("gotcentersKmeans")
    end
    @show Znl
    @show size(state)
    getΨ(Ψ, Znl, state, n_points, n_state, normalized); debug("Got Ψ")
    @ddshow sum(!isfinite(Ψ))
    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,y); debug("fitlinear")