Spectral estimation with Control Variate Variance Reduction

Whitening with various data and whiteners
In [1]:
using PyPlot

using Statistics

se = include("SpecEst.jl")
at = include("../../AnalysisToolbox.jl")
whf = include("../../WhiteningFilters.jl")

ut = include("../Autocov_Spec_Est_utils.jl")

Unable to init server: Could not connect: Connection refused
Unable to init server: Could not connect: Connection refused

(.:230074): Gdk-CRITICAL **: 16:56:24.840: gdk_cursor_new_for_display: assertion 'GDK_IS_DISPLAY (display)' failed

(.:230074): Gdk-CRITICAL **: 16:56:24.844: gdk_cursor_new_for_display: assertion 'GDK_IS_DISPLAY (display)' failed
Out[1]:
Main.Autocov_Spec_Est_utils
In [28]:
function Get_Specs(data_type; N = 10^5, Nen = 25, L = 1500, Nex = 2L+2, flags...)

    O = data_type(N; flags...)
    X = O.X
    spec= O.spec

    labels = ["Mesa", "DSP_spec"]

    S = Array{Any,1}(undef,2)
    cv = false

    @time S[1] = se.spec_est(X; estimator = at.spec_GB, Nex, cv, Nen)

    @time S[2] = se.spec_est(X; estimator = se.Mesa, Nex, cv, Nen)

    SS = copy(S);

    S = Array{Any,1}(undef,2)
    cv = true

    @time S[1] = se.spec_est(X; estimator = at.spec_GB, Nex, cv, Nen)

    @time S[2] = se.spec_est(X; estimator = se.Mesa, Nex, cv, Nen)

    SS_cv = copy(S);


    S = Array{Any,1}(undef,2)
    cv = true
    whole = true

    @time S[1] = se.spec_est(X; estimator = at.spec_GB, Nex, cv, Nen, whole)

    @time S[2] = se.spec_est(X; estimator = se.Mesa, Nex, cv, Nen, whole)

    SS_cv_whole = copy(S);

    return SS, SS_cv, SS_cv_whole, labels, spec
end

function DisplayPlots(SS, SS_cv, SS_cv_whole, labels, spec, Run = 1;
        Nex = 2000,
        zoombox = [pi - .1,pi + .1,1e-7, 1e-2],
        smo = at._smoother(4,5; ty = "bin", two_sided = false)
    )

    wgrid = at.Θ(Nex)
    Spec = spec.(wgrid);
    figure(figsize = (12,12))
    suptitle( "Run $Run")
    for i = 1:2
        subplot(2,2,2i-1)
        title("estimated spectra using "*labels[i])
        semilogy(SS[i].wgrid,SS[i].S[:], "k", label = "Sans CV")
        semilogy(SS_cv[i].wgrid,SS_cv[i].S[:], "--g", label = "CV averaged")
        semilogy(SS_cv[i].wgrid,at.smooth_spec(SS_cv[i].S[:],smo), ":g", label = "CV ave smo")
        semilogy(SS_cv_whole[i].wgrid,SS_cv_whole[i].S[:], "--r", label = "CV whole")
        semilogy(SS_cv_whole[i].wgrid,at.smooth_spec(SS_cv_whole[i].S[:],smo), ":r", label = "With CV whl smo")

        semilogy(wgrid,Spec,":k",label = "truth")
        ylabel("power")
        xlabel(L"\omega")
        legend()

        subplot(2,2,2i)
            title("estimated spectra using "*labels[i])
        semilogy(SS[i].wgrid,SS[i].S[:], "k", label = "Sans CV")
        semilogy(SS_cv[i].wgrid,SS_cv[i].S[:], "--g", label = "CV averaged")
        semilogy(SS_cv[i].wgrid,at.smooth_spec(SS_cv[i].S[:],smo), ":g", label = "CV ave smo")
        semilogy(SS_cv_whole[i].wgrid,SS_cv_whole[i].S[:], "--r", label = "CV whole")
        semilogy(SS_cv_whole[i].wgrid,at.smooth_spec(SS_cv_whole[i].S[:],smo), ":r", label = "With CV whl smo")

        semilogy(wgrid,Spec,":k",label = "truth")
        axis(zoombox)
        ylabel("power")
        xlabel(L"\omega"*" (near "*L"\pi)")
        legend()
    end
end
Out[28]:
DisplayPlots (generic function with 3 methods)

Guassian

In [20]:
sig =.5
K = 1000
f(x) = x < pi ? K/sqrt(2pi*sig)*exp(-x^2/sig) : K/sqrt(2pi*sig)*exp(-(2pi-x)^2/sig)
spec = x -> f(x) # divide by two if we require X real

N = 10^5

L = 1500
Nex = 2L+2
Nen = 25

SS, SS_cv, SS_cv_whole, labels, spec = Get_Specs(se.Guassian; N, Nen, L, Nex, spec, rl = true);

DisplayPlots(SS, SS_cv, SS_cv_whole, labels, spec, 1; Nex)

  0.008367 seconds (22.68 k allocations: 11.848 MiB)
  0.033011 seconds (1.76 k allocations: 972.906 KiB)
  7.248784 seconds (30.54 M allocations: 18.513 GiB, 7.70% gc time)
  4.983639 seconds (25.93 M allocations: 10.337 GiB, 8.22% gc time)
  7.623371 seconds (31.70 M allocations: 19.149 GiB, 7.48% gc time)
  5.362695 seconds (27.01 M allocations: 10.845 GiB, 7.94% gc time)
In [22]:
SS, SS_cv, SS_cv_whole, labels, spec = Get_Specs(se.Guassian; N, Nen, L, Nex, spec, rl = true);

DisplayPlots(SS, SS_cv, SS_cv_whole, labels, spec, 2; Nex)

  0.011168 seconds (22.68 k allocations: 11.848 MiB)
  0.029362 seconds (1.76 k allocations: 974.094 KiB)
  8.114507 seconds (30.59 M allocations: 18.646 GiB, 8.52% gc time)
  5.212410 seconds (26.09 M allocations: 10.703 GiB, 8.29% gc time)
  8.504349 seconds (31.76 M allocations: 19.287 GiB, 8.46% gc time)
  5.532595 seconds (27.13 M allocations: 11.125 GiB, 7.98% gc time)
In [24]:
SS, SS_cv, SS_cv_whole, labels, spec = Get_Specs(se.Guassian; N, Nen, L, Nex, spec, rl = true);

DisplayPlots(SS, SS_cv, SS_cv_whole, labels, spec, 3; Nex)

  0.011273 seconds (22.68 k allocations: 11.848 MiB)
  0.029580 seconds (1.76 k allocations: 970.125 KiB)
  8.092203 seconds (30.69 M allocations: 18.866 GiB, 8.36% gc time)
  5.223843 seconds (25.96 M allocations: 10.385 GiB, 11.18% gc time, 0.12% compilation time)
  7.900454 seconds (31.86 M allocations: 19.521 GiB, 8.78% gc time)
  5.452400 seconds (27.01 M allocations: 10.847 GiB, 8.93% gc time)

ARMA

In [34]:
L = 1500
Nex = 2L+2
Nen = 25
zoombox = [3pi/2*.9 - .1,3pi/2*.9 + .1,1e-9, 1e1]

SS, SS_cv, SS_cv_whole, labels, spec = Get_Specs(se.ARMA; N, Nen, L, Nex, rl = true);

DisplayPlots(SS,SS_cv,SS_cv_whole, labels, spec, 1; Nex, zoombox)

  0.008681 seconds (22.68 k allocations: 11.848 MiB)
  0.039064 seconds (1.76 k allocations: 974.281 KiB)
  7.993912 seconds (31.23 M allocations: 20.098 GiB, 9.92% gc time)
  6.566516 seconds (27.37 M allocations: 13.617 GiB, 8.55% gc time)
  8.252586 seconds (32.40 M allocations: 20.757 GiB, 8.39% gc time)
  6.982337 seconds (28.45 M allocations: 14.130 GiB, 8.39% gc time)
In [35]:
SS, SS_cv, SS_cv_whole, labels, spec = Get_Specs(se.ARMA; N, Nen, L, Nex, rl = true);

DisplayPlots(SS,SS_cv,SS_cv_whole, labels, spec, 2; Nex, zoombox)

  0.008659 seconds (22.68 k allocations: 11.848 MiB)
  0.038740 seconds (1.76 k allocations: 974.219 KiB)
  8.764456 seconds (31.37 M allocations: 20.419 GiB, 16.59% gc time)
  6.598355 seconds (27.42 M allocations: 13.715 GiB, 8.51% gc time)
  8.771727 seconds (32.55 M allocations: 21.087 GiB, 12.61% gc time)
  7.087696 seconds (28.51 M allocations: 14.259 GiB, 8.95% gc time)
In [36]:
SS, SS_cv, SS_cv_whole, labels, spec = Get_Specs(se.ARMA; N, Nen, L, Nex, rl = true);

DisplayPlots(SS,SS_cv,SS_cv_whole, labels, spec, 3; Nex, zoombox)

  0.034773 seconds (22.69 k allocations: 11.848 MiB, 64.41% gc time)
  0.032983 seconds (1.76 k allocations: 973.219 KiB)
  7.944414 seconds (31.14 M allocations: 19.895 GiB, 9.88% gc time)
  6.580645 seconds (27.43 M allocations: 13.750 GiB, 8.53% gc time)
  8.164479 seconds (32.31 M allocations: 20.537 GiB, 8.42% gc time)
  7.015068 seconds (28.53 M allocations: 14.303 GiB, 8.32% gc time)

KSE (Second mode)

In [307]:
N = 10^5

O = se.KSE(N; rl = true)
X = O.X

Sol load location: ../../../../data/KSE_Data/ks_pred_lin1e5_r1.jld
Out[307]:
1×100000 Matrix{ComplexF64}:
 -0.0622791-0.0665026im  -0.0602835-0.0630604im  …  -0.191963+0.0385498im
In [308]:
L = 1500
Nex = 2L + 1 + 1

S = Array{Any,1}(undef,2)
cv = false


labels = ["Mesa", "DSP_spec"]

@time S[1] = se.spec_est(X; estimator = at.spec_GB, Nex, cv, Nen, disc)

@time S[2] = se.spec_est(X; estimator = se.Mesa, Nex, cv, Nen, disc)

SS = copy(S);

  0.015303 seconds (22.67 k allocations: 16.493 MiB)
  0.052928 seconds (1.76 k allocations: 1.719 MiB)
In [309]:
S = Array{Any,1}(undef,2)
cv = true
Nen = 25 # div(N,Nen) = 4000
disc = 0

labels = ["Mesa", "DSP_spec"]

@time S[1] = se.spec_est(X; estimator = at.spec_GB, Nex, cv, Nen, disc)

# @time S[2] = se.spec_est(X; estimator = se.Mesa, Nex, cv, Nen, disc)

SS_cv = copy(S);

  6.877931 seconds (29.10 M allocations: 15.238 GiB, 11.12% gc time)
In [310]:
S = Array{Any,1}(undef,2)
cv = true
Nen = 25 # div(N,Nen) = 4000
whole = true

labels = ["Mesa", "DSP_spec"]

@time S[1] = se.spec_est(X; estimator = at.spec_GB, Nex, cv, Nen, whole)

# @time S[2] = se.spec_est(X; estimator = se.Mesa, Nex, cv, Nen, whole)

SS_cv_whole = copy(S);

  7.203617 seconds (30.23 M allocations: 15.806 GiB, 11.05% gc time)
In [314]:
labels = ["spec_GB", "Mesa"]
wgrid = at.Θ(Nex)
Spec = spec.(wgrid);
figure(figsize = (12,12))
suptitle( "KSE Run 1")
for i = 1:1
    subplot(2,2,2i-1)
    title("estimated spectra using "*labels[i])
    loglog(SS[i].wgrid,SS[i].S[:], "k", label = "Sans CV")
    semilogy(SS_cv[i].wgrid,SS_cv[i].S[:], "--g", label = "CV averaged")
    semilogy(SS_cv[i].wgrid,at.smooth_spec(SS_cv[i].S[:],smo), ":b", label = "CV ave smo")
    semilogy(SS_cv_whole[i].wgrid,SS_cv_whole[i].S[:], "--r", label = "CV whole")
    semilogy(SS_cv_whole[i].wgrid,at.smooth_spec(SS_cv_whole[i].S[:],smo), ":r", label = "With CV whl smo")

    ylabel("power")
    xlabel(L"\omega")
    legend()

#     subplot(2,2,2i)
#     title("estimated spectra using "*labels[i])
#     loglog(SS[i].wgrid,SS[i].S[:], "k", label = "Sans CV")
#     semilogy(SS_cv[i].wgrid,SS_cv[i].S[:], "--g", label = "CV averaged")
#     semilogy(SS_cv[i].wgrid,at.smooth_spec(SS_cv[i].S[:],smo), ":b", label = "CV ave smo")    
#     semilogy(SS_cv_whole[i].wgrid,SS_cv_whole[i].S[:], "--r", label = "CV whole")
#     semilogy(SS_cv_whole[i].wgrid,at.smooth_spec(SS_cv_whole[i].S[:],smo), ":r", label = "With CV whl smo")

#     axis([3pi/2*.9 - .2,3pi/2*.9 + .2,1e-8, 1e1])
#     ylabel("power")
#     xlabel(L"\omega"*" (near "*L"\pi)")
#     legend()

end
i = 2
subplot(2,2,2i-1)
title("estimated spectra using "*labels[i])
loglog(SS[i].wgrid,SS[i].S[:], "k", label = "Sans CV")
# semilogy(SS_cv[i].wgrid,SS_cv[i].S[:], "--g", label = "CV averaged")
# semilogy(SS_cv[i].wgrid,at.smooth_spec(SS_cv[i].S[:],smo), ":b", label = "CV ave smo")
# semilogy(SS_cv_whole[i].wgrid,SS_cv_whole[i].S[:], "--r", label = "CV whole")
# semilogy(SS_cv_whole[i].wgrid,at.smooth_spec(SS_cv_whole[i].S[:],smo), ":r", label = "With CV whl smo")

ylabel("power")
xlabel(L"\omega")
legend()
Out[314]:
PyObject <matplotlib.legend.Legend object at 0x7f674c49fee0>