
function plotresonancegroups(projects_dir)
    #import Clustering

    PyPlot.close("all")
    fig_num = 1

    Random.seed!(25)
    PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])




    #projects_dir = save_folder_path
    #projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/Nam2022_Serine"
    #projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/L-Serine-700"
    #projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/final/D-(+)-Glucose-NRC-600"


    plots_save_folder = joinpath(projects_dir, "plots")
    isdir(plots_save_folder) || mkdir(plots_save_folder)
    project_title = "resonance groupings"




    T = Float64

    ### load block.
    #load_path = joinpath(joinpath(projects_dir, project_name), "results_full.bson")
    load_path = joinpath(projects_dir, "results_full.bson")
    dict = BSON.load(load_path)
    As = collect( dict[:As][i] for i = 1:length(dict[:As]) )
    Es = collect( NMRSpectraSimulator.κCompoundFIDType(As[i]) for i = 1:length(As) )

    Δsys_cs = convert(Vector{Vector{Float64}}, dict[:Δsys_cs])
    y = convert(Vector{Complex{Float64}}, dict[:y])
    U_y = convert(Vector{Float64}, dict[:U_y])
    SW = dict[:SW]
    fs = dict[:fs]
    ν_0ppm = dict[:ν_0ppm]

    hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
    ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)
    P_y = hz2ppmfunc.(U_y)


    ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(As, hz2ppmfunc)
    ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))
    u_offset = 0.5
    u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
    u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)


    P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
    U = ppm2hzfunc.(P)

    p_star_set = dict[:p_star_set]
    κ_lb_default = dict[:κ_lb_default]
    κ_ub_default = dict[:κ_ub_default]
    κ_star_set = dict[:κ_star_set]
    d_star_set = dict[:d_star_set]
    β_star_set = dict[:β_star_set]
    λ_star_set = dict[:λ_star_set]
    cost_inds_set = dict[:cost_inds_set]
    w = dict[:w]

    As = collect( Es[n].core for n = 1:length(Es))
    ### end load block.

    hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
    ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)




    ### plot.


    f = uu->NMRSpectraSimulator.evalmixture(uu, As)

    # test params.
    ΩS_ppm = collect( hz2ppmfunc.( NMRSpectraSimulator.combinevectors(A.Ωs) ./ (2*π) ) for A in As )
    #ΩS_ppm_flat = NMRSpectraSimulator.combinevectors(ΩS_ppm)




    P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
    U = ppm2hzfunc.(P)

    ## parameters that affect qs.
    # A.ss_params.d, A.ss_params.κs_λ, A.ss_params.κs_β
    # A.d_singlets, A.αs_singlets, A.Ωs_singlets, A.β_singlets, A.λ0, A.κs_λ_singlets
    q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As)

    #g = uu->evalitpproxycompoundresonance(uu, As[1])
    A = As[1]
    g1 = uu->evalitpproxysysresonance(A.qs, uu, A.ss_params.d, A.ss_params.κs_λ, A.ss_params.κs_β)
    g2 = uu->evalsingletsresonance(uu, A.d_singlets, A.αs_singlets, A.Ωs_singlets,
    A.β_singlets, A.λ0, A.κs_λ_singlets)

    f_U = f.(U)
    q_U = q.(U)
    g1_U = g1.(U)
    g2_U = g2.(U)

    #g_qs_U, g_singlets_U = convertresonancetimeseries(g_U)


    z1 = Destruct.destruct(g1_U);
    q_singlets_evals = Destruct.destruct(g2_U);


    # collect( length(z1[i]) for i = 1:length(z1))
    # collect( length(z2[i]) for i = 1:length(z2))
    qs_evals = assembleqsevals(z1, A.qs)


    sanity_check = zeros(Complex{Float64}, length(q_U))
    if !isempty(z1)
        sanity_check += sum(z1)
    end
    if !isempty(q_singlets_evals)
        sanity_check += sum(q_singlets_evals)
    end
    @assert norm(sanity_check - q_U) < 1e-14
    println("sanity: ", norm(sanity_check - q_U))
    println()



    discrepancy = abs.(f_U-q_U)
    max_val, ind = findmax(discrepancy)
    println("relative discrepancy = ", norm(discrepancy)/norm(f_U))
    println("max discrepancy: ", max_val)
    println()

    ## visualize.
    PyPlot.plot(P, real.(f_U), label = "f")
    PyPlot.plot(P, real.(q_U), label = "q")

    #visualizesubgroups(P, qs_evals, q_singlets_evals, real)

    PyPlot.legend()
    PyPlot.xlabel("ppm")
    PyPlot.ylabel("real")
    PyPlot.title("f vs q")

    #@assert 1==2

    plots_save_path = joinpath(plots_save_folder, "groups_real.html")
    title_string = "$(project_title), real"
    savefiggroups(plots_save_path, title_string,
        P, q_U, qs_evals, q_singlets_evals, real)

    #
    plots_save_path = joinpath(plots_save_folder, "groups_imag.html")
    title_string = "$(project_title), imag"
    savefiggroups(plots_save_path, title_string,
        P, q_U, qs_evals, q_singlets_evals, real)
    #
    plots_save_path = joinpath(plots_save_folder, "groups_modulus.html")
    title_string = "$(project_title), modulus"
    savefiggroups(plots_save_path, title_string,
        P, q_U, qs_evals, q_singlets_evals, real)

end
