From 94c05c3d6fbc198d1fe6e49b9602a5cf29cfa1fc Mon Sep 17 00:00:00 2001 From: Gaspard Jankowiak Date: Fri, 15 Mar 2024 13:08:33 +0100 Subject: [PATCH] import julia files --- Data/PAR_POPC-test.toml | 46 ++++++ PLUV.jl | 355 ++++++++++++++++++++++++++++++++++++++++ Utils.jl | 86 ++++++++++ test.jl | 121 ++++++++++++++ 4 files changed, 608 insertions(+) create mode 100644 Data/PAR_POPC-test.toml create mode 100644 PLUV.jl create mode 100644 Utils.jl create mode 100644 test.jl diff --git a/Data/PAR_POPC-test.toml b/Data/PAR_POPC-test.toml new file mode 100644 index 0000000..6a10738 --- /dev/null +++ b/Data/PAR_POPC-test.toml @@ -0,0 +1,46 @@ +[metadata] +data_file = "./Data/POPC_sub_logstep_3.dat" +title = "POPC-test" +save_folder = "POPC-test" +q_min = 0.05 # (nm^-1) +q_max = 7.0 # (nm^-1) +model = "SDP_POPC_RecBuf" +temperature = 500 # (none=0) +binning = 10 +error = 1 # (0..1) +plot = 0 # (0-1) + + +parameters = [ +{ symbol = "Norm", initial_value = 1e5, fixed = true, prior = nan, min = 1e-2, max = 1e+5 }, # Normalization +{ symbol = "nv", initial_value = 8, fixed = false, prior = nan, min = 3, max = 13 }, # Vesicle concencentration x10^6 (nm^-3) +# Vesicle's parameters +{ symbol = "Rm", initial_value = 23, fixed = false, prior = nan, min = 20, max = 28 }, # Radius of the vesicle (nm) +{ symbol = "Z", initial_value = 6, fixed = false, prior = nan, min = 5.1, max = 17 }, # Radius polydispersity +# Head-group Positions & Widths +{ symbol = "n_TR", initial_value = 0.0, fixed = true, prior = nan, min = 0.00, max = 0.14 }, # Tris fraction +{ symbol = "d_TR", initial_value = 1.0, fixed = true, prior = nan, min = 0.4, max = 1.6 }, # Tris width (nm) +{ symbol = "s_TR", initial_value = 0.29, fixed = true, prior = nan, min = 0.20, max = 0.35 }, # Tris position (nm) +{ symbol = "d_Chol", initial_value = 0.22, fixed = false, prior = 0.05, min = 0.164, max = 0.273 }, # CholCH3 position (nm) +{ symbol = "s_Chol", initial_value = 0.3, fixed = false, prior = 0.05, min = 0.225, max = 0.375 }, # CholCH3 width (nm) +{ symbol = "d_PCN", initial_value = 0.41, fixed = false, prior = 0.05, min = 0.308, max = 0.512 }, # PCN position (nm) +{ symbol = "s_PCN", initial_value = 0.25, fixed = false, prior = nan, min = 0.23, max = 0.37 }, # PCN width (nm) +{ symbol = "d_CG", initial_value = 0.22, fixed = false, prior = nan, min = 0.18, max = 0.26 }, # CG position (nm) +{ symbol = "s_CG", initial_value = 0.25, fixed = false, prior = nan, min = 0.23, max = 0.37 }, # CG width (nm) +# Acyl chain Positions & Widths +{ symbol = "A_L", initial_value = 0.654, fixed = false, prior = 0.02, min = 0.589, max = 0.719 }, # Area per lipid (nm²) +{ symbol = "s_D_C", initial_value = 0.08, fixed = false, prior = nan, min = 0.05, max = 0.15 }, # HC excess std (nm) +{ symbol = "s_CH2", initial_value = 0.244, fixed = false, prior = 0.05, min = 0.183, max = 0.305 }, # HC smearing (nm) +{ symbol = "d_CH", initial_value = 0.90, fixed = true, prior = nan, min = 0.5, max = 0.95 }, # CH position (nm) +{ symbol = "s_CH", initial_value = 0.305, fixed = true, prior = nan, min = 0.15, max = 0.4 }, # CH width (nm) +{ symbol = "s_CH3", initial_value = 0.30, fixed = false, prior = nan, min = 0.23, max = 0.37 }, # CH3 width (nm) +# Area & Volumes +{ symbol = "r_PCN", initial_value = 0.27, fixed = false, prior = 0.05, min = 0.203, max = 0.338 }, # V_PCN/V_HL +{ symbol = "r_CG", initial_value = 0.48, fixed = false, prior = 0.05, min = 0.360, max = 0.600 }, # V_GG/V_HL +{ symbol = "r12", initial_value = 0.81, fixed = true, prior = nan, min = 0.770, max = 0.851 }, # V_CH/V_CH2 +{ symbol = "r32", initial_value = 1.93, fixed = false, prior = 0.05, min = 1.448, max = 2.415 }, # V_CH3/V_CH2 +{ symbol = "T", initial_value = 37, fixed = true, prior = nan, min = 32, max = 42 }, # Temperature (°C) +{ symbol = "V_BW", initial_value = 0.030, fixed = false, prior = nan, min = 0.0288, max = 0.0303 }, # Bound-water volume (nm³) +# +{ symbol = "Con", initial_value = 2.2e-4, fixed = false, prior = nan, min = 0.7e-4, max = 3.7e-4 } # Constant (nm⁻¹) +] diff --git a/PLUV.jl b/PLUV.jl new file mode 100644 index 0000000..c10179b --- /dev/null +++ b/PLUV.jl @@ -0,0 +1,355 @@ +module PLUV + +import SpecialFunctions: erf + +###################################################################################################### +###################################################################################################### +###################################################################################################### +############### POPC/POPG 95/5 mol/mol LUVs ############### +############### Reconstitution buffer: TRIS 20 mM, EDTA 2 mM ############### +############### Buffer used for PLUV exp. in DESY in 2017: ############### +###################################################################################################### +###################################################################################################### +###################################################################################################### + +const parameter_names = [ + :norm, :nv, :rm, :z, :n_tr, :d_tr, :s_tr, :d_chol, :s_chol, :d_pcn, :s_pcn, :d_cg, :s_cg, + :a_l, :s_d_c, :s_ch2, :d_ch, :s_ch, :s_ch3, :r_pcn, :r_cg, :r12, :r32, :t, :v_bw, :con +] + +############### GLOBAL VARIABLES ############### + +# POPG molar ratio +const CONST_x_PG = 0.05 + +# ############## POPC and POPG ############### +# 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine +# 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol + +# number of chain groups +const CONST_n_CH = 2 +const CONST_n_CH2 = 28 +const CONST_n_CH3 = 2 + +#X-ray scattering length chain groups (nm) +const CONST_b_CH = 1.97256E-05 +const CONST_b_CH2 = 2.25435E-05 +const CONST_b_CH3 = 2.53615E-05 + +### POPC +# Lipid-head volume +const CONST_V_HL_PC = 0.331 # 0.320 +# X-ray scattering length of head groups (nm) +const CONST_b_PC = 2.73340E-04 +const CONST_b_CG = 1.88802E-04 +const CONST_b_PCN = 1.97256E-04 +const CONST_b_Chol = 7.60844E-05 +# Lipid-volume temperature-dependencies a0 + a1*T (nm^3) +const CONST_a0_V_POPC = 1.22810311835285 +const CONST_a1_V_POPC = 0.000934915795086395 + +### POPG +# Lipid-head volume +const CONST_V_HL_PG = 0.289 +# X-ray Scattering length of head groups (nm) +const CONST_b_PG = 2.47979E-04 +const CONST_b_PG1 = 1.32443E-04 +const CONST_b_PG2 = 1.15536E-04 + +# Lipid-volume temperature-dependencies a0 + a1*T (nm^3) +const CONST_a0_V_POPG = 1.17881068602663 +const CONST_a1_V_POPG = 0.00108364914520327 + + +############### Other variables ############### +### Water +# V_PW = 24.5e-3 #(nm^3) Perkins 2001 (Hydration shell in proteins) +# polynome coefficient for T-dependency of bulk-water-molecule volume (V_HW) +# Units in degree Celsius +const CONST_p0_VW = 0.0299218 +const CONST_p1_VW = -2.25941e-06 +const CONST_p2_VW = 2.5675e-07 +const CONST_p3_VW = -1.69661e-09 +const CONST_p4_VW = 6.52029e-12 +# polynome coefficient for T-dependency of bulk-water molar concentration (Cw) +#Units in degree Celsius +const CONST_p0_Cw = 55.5052 +const CONST_p1_Cw = 0.00131894 +const CONST_p2_Cw = -0.000334396 +const CONST_p3_Cw = 9.10861e-07 + +const CONST_b_HW = 2.8179E-05 +const CONST_d_shl = 0.31 # (nm) Perkins 2001 (Hydration shell in proteins) + +# Composition of the Reconstitution buffer (M) +const CONST_ctris = 0.02 +const CONST_cEDTA = 0.002 + +### Extra molecules + +# TRIS buffer +const CONST_b_tris = 1.860E-04 +const CONST_V_tris = 0.15147 # (nm^3) +# EDTA +const CONST_b_EDTA = 4.340E-04 +const CONST_V_EDTA = 0.56430 # (nm^3) + + +################################################################################################################## +################################################################################################################## + +######################################################### +#@njit(parallel=True) +function PDF_normal(x, mu, sig) + return exp(-(x - mu)^2 / (2 * sig^2)) / (sig * sqrt(2 * pi)) +end + +######################################################### +#@njit(parallel=True) +function PDF_lognormal(x, mu_x, sig_x) + # https://en.wikipedia.org/wiki/Log-normal_distribution + mu = log(mu_x^2 / sqrt(mu_x^2 + sig_x^2)) + sig = sqrt(log(1 + sig_x^2 / mu_x^2)) + return exp(-(log(x) - mu)^2 / (2 * sig^2)) / (x * sig * sqrt(2 * pi)) +end + +######################################################### +function lipid_volume(T) + return (1 - CONST_x_PG) * (CONST_a0_V_POPC + T * CONST_a1_V_POPC) + CONST_x_PG * (CONST_a0_V_POPG + T * CONST_a1_V_POPG) +end + +######################################################### +function water_volume(T) + return CONST_p0_VW + CONST_p1_VW * T + CONST_p2_VW * T^2 + CONST_p3_VW * T^3 + CONST_p4_VW * T^4 +end + +######################################################### +#@njit(parallel=True) +function FTreal_erf(q, mu, d, sig) + """ FTreal_erf(q, mu, d, sig) """ + return sin(q * d / 2.0) / (q * d / 2.0) * exp(-(q * sig)^2 / 2.0) * cos(q * mu) +end + +######################################################### +#@njit(parallel=True) +function FTreal_gauss(q, mu, sig) + """ FTreal_gauss(q, mu, sig) """ + return exp(-(q * sig)^2 / 2.0) * cos(q * mu) +end + +######################################################### +function Slab(x, mu, L, sig) + return 0.5 * (erf((x - (mu - L / 2.0)) / (sqrt(2) * sig)) - erf((x - (mu + L / 2)) / (sqrt(2) * sig))) +end + +######################################################### +function Gauss(x, V, mu, sig, A_L) + return V * PDF_normal(x, mu, sig) / A_L +end + +#################################### j0^2 MOMENTS (SCHULZ PDF) ######################################### +#################################### (checked!) ######################################### + +#@njit(parallel=True) +function mu4(q, Z, a) + return a^2 * (Z + 2) * (Z + 1) * (1 - (1 + 4 * q^2 * a^2)^(-(Z + 3) / 2.0) * cos((Z + 3) * atan(2 * q * a))) / (2 * q^2) +end + +################################################################################################################## +################################################################################################################## + +################################ SYMMETRIC VESICLE FOR X-RAY SLDs ######################################### + +######################################## SDP MODELLING ######################################### +#################################### SEPARATED FORM FACTOR ######################################### + +################################################################################################################## +################################################################################################################## + +######################################################### +######### Symmetric POPC bilayer ######################## +######### liposomes and proteoliposomes ################# +######################################################### + +################## +function intensity(P::NamedTuple, q) + return intensity(Float64[P...], q) +end + +function intensity(P::AbstractArray{Float64}, q) + + Norm, nv, Rm, Z, n_TR, d_TR, s_TR, d_Chol, s_Chol, d_PCN, s_PCN, d_CG, s_CG, + A_L, s_D_C, s_CH2, d_CH, s_CH, s_CH3, r_PCN, r_CG, r12, r32, T, V_BW, Con = P + + # [Example 1], fixed parameters + # Norm 1e5 # Normalization + # n_TR 0.0 # Tris fraction + # d_TR 1.0 # Tris width (nm) + # s_TR 0.29 # Tris position (nm) + # d_CH 0.90 # CH position (nm) + # s_CH 0.305 # CH width (nm) + # r12 0.81 # V_CH/V_CH2 + # T 37 # Temperature (°C) + + # [example 1] fixed + Cw = CONST_p0_Cw + CONST_p1_Cw * T + CONST_p2_Cw * T^2 + CONST_p3_Cw * T^3 + xtris = CONST_ctris / Cw # mole fraction of free TRIS in bulk + xEDTA = CONST_cEDTA / Cw # mole fraction of free EDTA in bulk + + # Volumes + # [example 1] fixed + V_L = lipid_volume(T) + V_HW = water_volume(T) + V_HC = V_L - ((1 - CONST_x_PG) * CONST_V_HL_PC + CONST_x_PG * CONST_V_HL_PG) + + # Calculation of mean D_C + D_C = V_HC / A_L + + # Quasi-molecular volumes + # [example 1] r12 fixed + V_CH2 = V_HC / (CONST_n_CH2 + CONST_n_CH * r12 + CONST_n_CH3 * r32) # Volume of CH2 groups + V_CH = V_CH2 * r12 # Volume of CH groups + V_CH3 = V_CH2 * r32 # Volume of CH3 groups + + V_CG = CONST_V_HL_PC * r_CG # Volume of CG group + V_PCN = CONST_V_HL_PC * r_PCN # Volume of PCN group + V_Chol = CONST_V_HL_PC * (1 - r_PCN - r_CG) # Volume of CholCH3 group + + # CONST + CONST_V_PG1 = CONST_V_HL_PG * 0.16 # Kucerka 2012 + CONST_V_PG2 = CONST_V_HL_PG * (1 - 0.51 - 0.16) # Kucerka 2012 + + ############### X-ray scattering lengths (nm) + # [example 1] rho_sol fixed + rho_sol = (CONST_b_HW + xtris * CONST_b_tris + xEDTA * CONST_b_EDTA) / V_HW + drho_Chol = ((1 - CONST_x_PG) * CONST_b_Chol / V_Chol + CONST_x_PG * CONST_b_PG2 / CONST_V_PG2) - rho_sol + drho_PCN = ((1 - CONST_x_PG) * CONST_b_PCN / V_PCN + CONST_x_PG * CONST_b_PG1 / CONST_V_PG1) - rho_sol + drho_CG = CONST_b_CG / V_CG - rho_sol + drho_TR = CONST_b_tris / CONST_V_tris - rho_sol + drho_CH = CONST_b_CH / V_CH - rho_sol + drho_CH2 = CONST_b_CH2 / V_CH2 - rho_sol + drho_CH3 = CONST_b_CH3 / V_CH3 - rho_sol + drho_HW = CONST_b_HW / V_BW - rho_sol + + ############### D_C polydispersity + # Number of integration points + N = 201 + + # Parameter dependent integration bounds + # HC_min, HC_max = D_C-3*s_D_C, D_C+3*s_D_C + # + # Hardcoded integration bounds for POPC, fixed T = 37, A_L in [0.598, 0.719] + HC_min, HC_max = V_HC / 0.719 - 0.5, V_HC / 0.589 + 0.5 + # Minimum and maximum samples (offset by 1/2 of the integration step w.r.t integration bounds for the midpoint rule) + integration_step = (HC_max - HC_min) / N + HC_first, HC_last = HC_min + 0.5 * integration_step, HC_max - 0.5 * integration_step + + HC_array = range(HC_first, HC_last, length=N) + Normal = PDF_normal.(HC_array, D_C, s_D_C) + + ############### c-prefactors + c_Chol = ((1 - CONST_x_PG) * V_Chol + CONST_x_PG * CONST_V_PG2) / A_L + c_PCN = ((1 - CONST_x_PG) * V_PCN + CONST_x_PG * CONST_V_PG1) / A_L + c_CG = V_CG / A_L + c_TR = CONST_V_tris * n_TR / A_L + + c_CH = V_CH * CONST_n_CH / V_HC * HC_array + c_CH3 = V_CH3 * CONST_n_CH3 / V_HC * HC_array + + d_comp = d_CG + d_PCN + d_Chol + CONST_d_shl + + ############### calculating scattering amplitude ----------------------------------------------- + Am = zeros(size(q, 1), N) + + for i in 1:N + + # Adding hydrocarbon-chain envelope + Am[:, i] += 2 * drho_CH2 * HC_array[i] * FTreal_erf.(q, 0, 2 * HC_array[i], s_CH2) + # Adding CH and CH3 groups + Am[:, i] += 2 * (drho_CH - drho_CH2) * c_CH[i] * FTreal_gauss.(q, d_CH, s_CH) + Am[:, i] += 2 * (drho_CH3 - drho_CH2) * c_CH3[i] * FTreal_gauss.(q, 0, s_CH3) + + # Adding hydration-water envelope + Am[:, i] += 4 * drho_HW * d_comp * FTreal_erf.(q, (HC_array[i] + d_comp / 2.0), d_comp, s_CH2) + # Adding CG, PCN and CholCH3 groups + Am[:, i] += 2 * (drho_TR - drho_HW) * c_TR * FTreal_gauss.(q, (HC_array[i] + d_TR / 2.0), s_TR) + Am[:, i] += 2 * (drho_CG - drho_HW) * c_CG * FTreal_gauss.(q, (HC_array[i] + d_CG / 2.0), s_CG) + Am[:, i] += 2 * (drho_PCN - drho_HW) * c_PCN * FTreal_gauss.(q, (HC_array[i] + d_CG + d_PCN / 2.0), s_PCN) + Am[:, i] += 2 * (drho_Chol - drho_HW) * c_Chol * FTreal_gauss.(q, (HC_array[i] + d_CG + d_PCN + d_Chol / 2.0), s_Chol) + + ############### Ensemble average + # multiply each columns of Am by Normal and sum along the columns, + # then multiply by integration step + I = vec(sum(Am .^ 2 .* Normal'; dims=2)) * integration_step + + alp = Rm / (Z + 1) + return @. (Norm * nv * 1e-6) * I * (16 * pi^2 * mu4(q, Z, alp)) + Con * (0.99 * (1.0 / (1 + exp(-8 * (q - 1.0)))) + 0.01) + end +end + +################## +function negative_water(P::NamedTuple) + return negative_water([P...]) +end + +function negative_water(P::AbstractArray{Float64}) + + _, _, _, _, n_TR, d_TR, s_TR, d_Chol, s_Chol, d_PCN, s_PCN, d_CG, s_CG, + A_L, _, s_CH2, _, _, _, r_PCN, r_CG, _, _, T, _, _ = P + + # Volumes + V_L = lipid_volume(T) + V_HC = V_L - ((1 - CONST_x_PG) * CONST_V_HL_PC + CONST_x_PG * CONST_V_HL_PG) + + # Calculation of mean D_C + D_C = V_HC / A_L + + # Quasi-molecular volumes + V_CG = CONST_V_HL_PC * r_CG # Volume of CG group + V_PCN = CONST_V_HL_PC * r_PCN # Volume of PCN group + V_Chol = CONST_V_HL_PC * (1 - r_PCN - r_CG) # Volume of CholCH3 group + + check = 0 + z_array = range(0.0, 4.0, length=81) + + CG = Gauss.(z_array, V_CG, D_C + d_CG / 2.0, s_CG, A_L) + PCN = Gauss.(z_array, V_PCN, D_C + d_CG + d_PCN / 2.0, s_PCN, A_L) + Chol = Gauss.(z_array, V_Chol, D_C + d_CG + d_PCN + d_Chol / 2.0, s_Chol, A_L) + TRIS = Gauss.(z_array, n_TR * CONST_V_tris, D_C + d_TR / 2.0, s_TR, A_L) + BW = Slab.(z_array, D_C + (d_CG + d_PCN + d_Chol + CONST_d_shl) / 2.0, d_CG + d_PCN + d_Chol + CONST_d_shl, s_CH2) - CG - PCN - Chol - TRIS + + for i in (BW) + if i < -1e-3 + check += 1 + end + end + + return check +end + +function get_free_parameters_idx(metadata::Dict) + m = metadata["parameters"] + return filter(i -> !m[i]["fixed"], eachindex(m)) +end + +function reduce_to_free_parameters(metadata::Dict, params_init::NamedTuple, lb, ub, q) + default_params = Float64[params_init...] + free_idx = get_free_parameters_idx(metadata) + P_reduced = view(default_params, free_idx) + lb_reduced = lb[free_idx] + ub_reduced = ub[free_idx] + function reduced(P::AbstractArray{Float64}) + P_reduced .= P + return intensity(default_params, q) + end + return reduced, P_reduced, lb_reduced, ub_reduced +end + +function reduce_to_free_parameters(metadata::Dict, params::Vector{Float64}) + free_idx = get_free_parameters_idx(metadata) + return params[free_idx] +end + +end # module PLUV + +# vim ts=4,sts=4,sw=4 diff --git a/Utils.jl b/Utils.jl new file mode 100644 index 0000000..7692bbd --- /dev/null +++ b/Utils.jl @@ -0,0 +1,86 @@ +module Utils + +import TOML: parsefile +import DelimitedFiles: readdlm +import DSP + +function load_init_params(filename::String) + meta = parsefile(filename)["metadata"] + params = NamedTuple(Dict(Symbol(p["symbol"]) => p["initial_value"] for p in meta["parameters"])) + lower_bounds = [p["min"] for p in meta["parameters"]] + upper_bounds = [p["max"] for p in meta["parameters"]] + return (meta=meta, params=params, lower_bounds=lower_bounds, upper_bounds=upper_bounds) +end + +function select_columns(QIE_array, qmin, qmax, bin, err_mul) + #""Select q / I(q) / err(q) data - Change I(q) value from mm^-1 to nm^-1' + + Q = Float64[] + IDATA = Float64[] + ERR = Float64[] + + for i in axes(QIE_array, 1) + # keep only the values of q within bounds + if qmin <= QIE_array[i, 1] <= qmax + # skip all but the (k*bin)-th values + if i % bin == 0 + if QIE_array[i, 2] > 0 + push!(Q, QIE_array[i, 1]) + push!(IDATA, 1e-6 * QIE_array[i, 2]) + if err_mul == 0 + push!(ERR, 0.01 * 1e-6 * QIE_array[i, 2]) + elseif QIE_array[i, 3] * err_mul > 0.01 * QIE_array[i, 2] + push!(ERR, 1e-6 * QIE_array[i, 3] * err_mul) + else + push!(ERR, 0.01 * 1e-6 * QIE_array[i, 2]) + end + end + end + end + end + + return (Q, IDATA, ERR) +end + +""" + compute_logscale_weights(q) + +Compute the weights W, which are proportional to the distance in the log space +""" +function compute_logscale_weights(q) + diff_log_q = diff(log10.(q)) + left_diff_log_q = view(diff_log_q, [firstindex(diff_log_q); eachindex(diff_log_q)]) + right_diff_log_q = view(diff_log_q, [eachindex(diff_log_q); lastindex(diff_log_q)]) + + return left_diff_log_q + right_diff_log_q +end + +function load_config(filename::String) + meta, params_init, lower_bounds, upper_bounds = load_init_params(filename) + f = open(meta["data_file"], "r") + qie_data = readdlm(f, '\t') + close(f) + return (meta=meta, + params_init=params_init, + lower_bounds=lower_bounds, + upper_bounds=upper_bounds, + qie_deta=qie_data) +end + +function lowpass_filter(i; σ=1) + half_w = 10 + x = range(-5, 5; length=2half_w + 1) + i_padded = view(i, [fill(firstindex(i), half_w); eachindex(i); fill(lastindex(i), half_w)]) + k = exp.(-0.5 * ((x / σ) .^ 2)) / (2 * σ * sqrt(2π)) + k ./= sum(k) + @show sum(k) + r = DSP.conv(i_padded, k)[2half_w+1:end-2half_w] + @show size(r) + return r +end + +function chi2(I_data, I_model, err) + return sum((I_data .- I_model ./ err) .^ 2) +end + +end # module Utils diff --git a/test.jl b/test.jl new file mode 100644 index 0000000..c6a8d4b --- /dev/null +++ b/test.jl @@ -0,0 +1,121 @@ +import GLMakie as M +import Metaheuristics as MH +import DelimitedFiles: readdlm + +include("PLUV.jl") +include("Utils.jl") + + +# filtered_idx = d[:, 2] .> 0.0 +# @views q, I, err = d[filtered_idx, 1], d[filtered_idx, 2], d[filtered_idx, 3] + +meta, params_init, lower_bounds, upper_bounds, qie_data = Utils.load_config("Data/PAR_POPC-test.toml") + +best_5k = begin + f_χ2 = open("POPC-test-5k/Results_collection-X2.dat", "r") + + d_χ2 = readdlm(f_χ2) + close(f_χ2) + + idx_min_χ2 = argmin(d_χ2)[1] + @show idx_min_χ2 + + f_params = open("POPC-test-5k/Results_collection.dat", "r") + for i in 1:(idx_min_χ2-1) + readline(f_params) + end + best_params = map(x -> parse(Float64, x), split(readline(f_params), ' ')) + @show best_params + + PLUV.reduce_to_free_parameters(meta, best_params) +end + +mean_5k = begin + f_params = open("POPC-test-5k/Results_collection.dat", "r") + params_pop = readdlm(f_params) + N_trials, _ = size(params_pop) + + mean_params = vec(sum(params_pop; dims=1) / N_trials) + + PLUV.reduce_to_free_parameters(meta, mean_params) +end + +# best_5k = PLUV.reduce_to_free_parameters(meta, [100000.0000000, 7.5004085, 23.5074449, 9.8664991, 0.0000000, 1.0000000, 0.2900000, 0.2203793, +# 0.2990402, 0.4122816, 0.3266636, 0.2276763, 0.2481895, 0.6642123, 0.1203572, 0.2629861, 0.9000000, 0.3050000, +# 0.2945315, 0.2762371, 0.4777401, 0.8100000, 2.0149998, 37.0000000, 0.0299139, 0.0002171]) +# +# best_5k = PLUV.reduce_to_free_parameters(meta, [100000.0000000, 7.6975592, 23.4041912, 9.7275630, 0.0000000, 1.0000000, 0.2900000, 0.2219937, 0.3000114, +# 0.4158804, 0.3278631, 0.2296156, 0.2475607, 0.6664143, 0.1191859, 0.2618609, 0.9000000, 0.3050000, 0.2963982, 0.2770345, 0.4762528, +# 0.8100000, 1.9706651, 37.0000000, 0.0299179, 0.0002167]) + + +q, I_data, err = Utils.select_columns(qie_data, meta["q_min"], meta["q_max"], meta["binning"], meta["binning"]) +q_all, I_all, err_all = Utils.select_columns(qie_data, meta["q_min"], meta["q_max"], 1, meta["binning"]) + +I_data_lp = Utils.lowpass_filter(I_data; σ=2) +I_all_lp = Utils.lowpass_filter(I_all; σ=2) + +w = Utils.compute_logscale_weights(q) +w_all = Utils.compute_logscale_weights(q_all) + +I_init = PLUV.intensity(params_init, q) + +intensity_reduced, P_reduced, lb_reduced, ub_reduced = PLUV.reduce_to_free_parameters(meta, params_init, lower_bounds, upper_bounds, q) + +bounds = MH.boxconstraints(lb=lb_reduced, ub=ub_reduced) + +function obj(P) + I_model = intensity_reduced(P) + return Utils.chi2(I_data, I_model, err) +end + +information = MH.Information(f_optimum=0.0) +options = MH.Options(f_calls_limit=1_000_000, f_tol=1e-5); +algorithm = MH.ECA(information=information, options=options) + +I_best_5k = intensity_reduced(best_5k) +I_mean_5k = intensity_reduced(mean_5k) + +if false + result = MH.optimize(obj, bounds, algorithm) + + P_best = MH.minimizer(result) + I_best = intensity_reduced(P_best) +end + +if true + fig = M.Figure() + ax = M.Axis(fig[1, 1]; xscale=log10, yscale=log10) + + M.lines!(ax, q, I_best, label="MH best (julia)") + M.scatter!(ax, q, I_data, label="data") + M.lines!(ax, q, I_best_5k, label="TSA best (5k)") + M.lines!(ax, q, I_mean_5k, label="TSA mean (5k)") + + M.axislegend() + + display(fig) +end + + +if false + fig = M.Figure() + ax = M.Axis(fig[1, 1]; xscale=log10, yscale=log10) + + M.scatter!(ax, q, I) + M.lines!(ax, q_all, I_all) + M.lines!(ax, q_all, I_all_lp) + #M.lines!(ax, q, I_init) + + display(fig) +end + +if false + fig = M.Figure() + ax = M.Axis(fig[1, 1]; xscale=log10) + + #M.scatter!(ax, q, w) + M.scatter!(ax, q_all, w_all) + + display(fig) +end