From d90732ab342fe6d2bfd4ce5e478fe5b5dd560e87 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson <cont-frb@ulund.org> Date: Tue, 13 Sep 2016 16:02:37 +0200 Subject: [PATCH] Clean up --- src/MCMC/testStan.jl | 109 ------------------------------------ src/PCA.jl | 14 +++++ src/SystemIdentification.jl | 2 +- src/armax.jl | 21 ++++--- src/cuckooSearch.jl | 1 + src/idinput.jl | 6 +- src/kernelPCA.jl | 12 ---- src/least_squares.jl | 8 ++- src/sobol.jl | 13 ++--- src/toeplitz.jl | 11 ---- src/transfer_functions.jl | 14 ++--- 11 files changed, 50 insertions(+), 161 deletions(-) delete mode 100644 src/MCMC/testStan.jl delete mode 100644 src/kernelPCA.jl delete mode 100644 src/toeplitz.jl diff --git a/src/MCMC/testStan.jl b/src/MCMC/testStan.jl deleted file mode 100644 index f2e8cbc..0000000 --- a/src/MCMC/testStan.jl +++ /dev/null @@ -1,109 +0,0 @@ -######### Stan program example ########### - -# module Tmp -const σw0 = 1.0 -const σw = 1 -const σv = 1.0 -const theta0 = [0.5, 25, 8] -s2piσv = log(sqrt(2*pi) * σv) -ProjDir = pwd() - -function f_sample(x::Vector, t::Int64) - c = 8*cos(1.2*(t-1)) - @inbounds for i = 1:length(x) - x[i] = 0.5*x[i] + 25*x[i]./(1+x[i]^2) + c + σw*randn() - end - return x -end -f_sample(x::Float64, t::Int64) = 0.5*x + 25*x/(1+x^2) + 8*cos(1.2*(t-1)) + σw*randn() - - - -T = 100 -M = 1 -x = Array(Float64,T) -y = Array(Float64,T) -x0 = 0 - -x[1] = σw*randn() -y[1] = σv*randn() -for t = 1:T-1 - x[t+1] = f_sample(x[t],t) - y[t+1] = 0.05x[t+1]^2 + σv*randn() -end # t - - - - -using Stan, Mamba - - -odemodel =" -data { - int<lower=1> T; - real y[T]; -} - -parameters { - real theta[3]; - real sigma[2]; - real x[T]; -} - -model { - sigma[1] ~ cauchy(1,1); - sigma[2] ~ cauchy(1,1); - theta[1] ~ cauchy(0.5,0.5); - theta[2] ~ cauchy(25,25); - theta[3] ~ cauchy(8,8); - x[1] ~ normal(0,1); - y[1] ~ normal(0.05*x[1]*x[1], sigma); - for (t in 1:(T-1)){ - x[t+1] ~ normal(theta[1]*x[t] + theta[2]*x[t]/(1+x[t]*x[t]) + theta[3]*cos(1.2*(t-1)),sigma[1]); - y[t+1] ~ normal(0.05*x[t+1]*x[t+1], sigma[2]); - } -} -" - -odedict = - Dict( - "T" => T, - "y" => y) - -stanmodel = Stanmodel(name="ode", model=odemodel, nchains=4, update=10000) -@time sim1 = stan(stanmodel, [odedict], ProjDir, diagnostics=false, CmdStanDir=CMDSTAN_HOME) - -## Subset Sampler Output to variables suitable for describe(). -monitor = ["lp__", "accept_stat__"] -sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)] - -describe(sim) -println() - -p = plot(sim, [:trace, :mean, :density, :autocor], legend=true) -draw(p, ncol=4, filename="$(stanmodel.name)-infoplot", fmt=:pdf) - - - -## Subset Sampler Output to variables suitable for describe(). -monitor = ["theta.1","theta.2","theta.3"] -sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)] - -describe(sim) -println() - -p = plot(sim, [:trace, :mean, :density, :autocor], legend=true) -draw(p, ncol=4, filename="$(stanmodel.name)-thetaplot", fmt=:pdf) - - - -## Subset Sampler Output to variables suitable for describe(). -monitor = ["sigma.1", "sigma.2"] -sim = sim1[:, monitor, :] -describe(sim) -println() - -p = plot(sim, [:trace, :mean, :density, :autocor], legend=true) -draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-sigmaplot", fmt=:pdf) - - diff --git a/src/PCA.jl b/src/PCA.jl index c852e7d..4505982 100644 --- a/src/PCA.jl +++ b/src/PCA.jl @@ -6,3 +6,17 @@ function PCA(W) score = score*diagm(latent) C,score,latent,W0 end + + +using MLKernels +function kernelPCA(X; α=1.0) + κ = GaussianKernel(α) + K = kernelmatrix(κ,X) + N = size(K)[1] + In = fill(1/N,(N,N)) + K = K-In*K - K*In + In*K*In # Make sure mean is zero + + (D,V) = eig(K) + Kpc = K*V + Kpc,D,V +end diff --git a/src/SystemIdentification.jl b/src/SystemIdentification.jl index 7841796..053144e 100644 --- a/src/SystemIdentification.jl +++ b/src/SystemIdentification.jl @@ -13,7 +13,7 @@ FitResult, IdData, # Functions ar,arx,getARregressor,getARXregressor,find_na, -toeplitz, kalman, kalman_smoother, forward_kalman, PCA +toeplitz, kalman, kalman_smoother, forward_kalman, PCA, plotmodel ## Fit Methods ================= :LS diff --git a/src/armax.jl b/src/armax.jl index 255572a..f716403 100644 --- a/src/armax.jl +++ b/src/armax.jl @@ -93,7 +93,7 @@ getARXregressor(iddata::IdData, na, nb) = getARXregressor(iddata.y,iddata.u, na, """Plots the RMSE and AIC for model orders up to `n`. Useful for model selection""" -function find_na(y::Vector{Float64},n::Int) +function find_na(y::Vector,n::Int) error = zeros(n,2) for i = 1:n w,e = ar(y,i) @@ -102,18 +102,21 @@ function find_na(y::Vector{Float64},n::Int) print(i,", ") end println("Done") - plotsub(error,"-o") - show() + scatter(error, show=true) end -import PyPlot -function PyPlot.plot(y,m::AR) +""" + plotmodel(y,m::AR) + +Plots a signal `y` and the output of the model `m` +""" +function plotmodel(y,m::AR) na = length(m.a) y,A = getARregressor(y,na) yh = A*m.a error = y-yh - newplot(y,"k") - plot(yh,"b") - plot(error,"r") - title("Fitresult, AR, n_a: $na, RMSE = $(rms(error)) Fit = $(fit(y,yh))") + plot(y,c=:black) + plot(yh,c=:b) + plot(error,c=:r, title="Fitresult, AR, n_a: $na, RMSE = $(rms(error)) Fit = $(fit(y,yh))") + end diff --git a/src/cuckooSearch.jl b/src/cuckooSearch.jl index 858249d..f813651 100644 --- a/src/cuckooSearch.jl +++ b/src/cuckooSearch.jl @@ -1,6 +1,7 @@ using Devectorize """ `cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)`\n + `n` = Number of nests (or different solutions) `pa=0.25` Discovery rate of alien eggs/solutions Change this if you want to get better results diff --git a/src/idinput.jl b/src/idinput.jl index 714497e..74f21bc 100644 --- a/src/idinput.jl +++ b/src/idinput.jl @@ -1,3 +1,8 @@ +""" + idinput(N, class = "prbs"; band = 1) + +Create a input signal for a sys.id. experiment, currently only supports prbs signals. +""" function idinput(N, class = "prbs"; band = 1) output = zeros(N) if lowercase(class) == "prbs" @@ -18,4 +23,3 @@ function idinput(N, class = "prbs"; band = 1) return output end - diff --git a/src/kernelPCA.jl b/src/kernelPCA.jl deleted file mode 100644 index f6d727e..0000000 --- a/src/kernelPCA.jl +++ /dev/null @@ -1,12 +0,0 @@ -using MLKernels -function kernelPCA(X; α=1.0) - κ = GaussianKernel(α) - K = kernelmatrix(κ,X) - N = size(K)[1] - In = fill(1/N,(N,N)) - K = K-In*K - K*In + In*K*In # Make sure mean is zero - - (D,V) = eig(K) - Kpc = K*V - Kpc,D,V -end diff --git a/src/least_squares.jl b/src/least_squares.jl index 4a0395d..9eec77b 100644 --- a/src/least_squares.jl +++ b/src/least_squares.jl @@ -36,9 +36,13 @@ function dsvd(A,b,λ) dsvd(A,U,S,V,b,λ) end +""" + Lcurve(normE, normX, λ) + +Plots the L-curve. normE and normX are obtained from, e.g., `dsvd` +""" function Lcurve(normE, normX, λ) - plot(normE,normX,xscale=:log10,yscale=:log10,m=:o) + plot(normE,normX,xscale=:log10,yscale=:log10,m=:o, xlabel="RMSE", ylabel="||k||", title="L-curve") annotations = [(normE[i],normX[i],"λ=$(round(λ[i],8))") for i in 1:length(λ)] annotate!(annotations) - xlabel!("RMSE"); ylabel!("||k||"); title!("L-curve") end diff --git a/src/sobol.jl b/src/sobol.jl index 5ab9290..360ed24 100644 --- a/src/sobol.jl +++ b/src/sobol.jl @@ -194,21 +194,18 @@ function InitSobol(dim::Int64) end """ -`using Winston` Plots the first 4 outputs of sobel(2,512) """ function test_sobol() X, nextseed, MeM = sobol(2,512) - figure(); pp = plot(X[:,1],X[:,2],".b"); hold(true) + plot(X[:,1],X[:,2],c=:b) X, nextseed, MeM = sobol(X, nextseed, MeM ) - plot(X[:,1],X[:,2],".g") + plot(X[:,1],X[:,2],c=:g) X, nextseed, MeM = sobol(X, nextseed, MeM ) - plot(X[:,1],X[:,2],".r") + plot(X[:,1],X[:,2],c=:r) X, nextseed, MeM = sobol(X, nextseed, MeM ) - plot(X[:,1],X[:,2],".m") - - display(pp) -end \ No newline at end of file + plot(X[:,1],X[:,2],c=:m) +end diff --git a/src/toeplitz.jl b/src/toeplitz.jl deleted file mode 100644 index c4cfa70..0000000 --- a/src/toeplitz.jl +++ /dev/null @@ -1,11 +0,0 @@ -function toeplitz{T}(c::Array{T},r::Array{T}) - nc = length(c) - nr = length(r) - A = zeros(T, nc, nr) - A[:,1] = c - A[1,:] = r - for i in 2:nr - A[2:end,i] = A[1:end-1,i-1] - end - A -end diff --git a/src/transfer_functions.jl b/src/transfer_functions.jl index eb0814b..3d66fbf 100644 --- a/src/transfer_functions.jl +++ b/src/transfer_functions.jl @@ -1,5 +1,11 @@ using DSP +""" + tfest(y,u) + +Estimate a transfer function model using the Correlogram approach +H = Syu/Suu +""" function tfest(y,u) Cyu = xcorr(y,u) Cuu = xcorr(u,u) @@ -9,11 +15,3 @@ function tfest(y,u) end tfest(iddata::IdData) = tfest(iddata.y,iddata.u) - - -N = 200000; -u = randn(N); -y = filt(ones(5),5,u); - -H = tfest(y,u); -semilogy(H.F,abs(H.P)) -- GitLab