initial_forcing.jl

initial_forcing = [64.79861184519945, -0.03956030217266119, -0.054641255219961, 0.1438437598358666, -0.006979979733325994, -0.018870338457863758, 0.0018286864505904798, 0.04255465840503827, 0.3659851712160052, 1.1678154520629573, -0.09920684852949448, -37.04190904240955, 0.134931019548738, -5.808707055633873, 0.04197231403115209, -0.191235897674358, -0.0005651295600417024, -0.09171877227387193, 8.1351216950623, 0.000826587745230315, 4.979592748986883, 0.3292794160236103, 0.06762453753371207, 0.168831120814038, 0.004048558625848306, 0.0008675782119757018, -0.01928763037706797, 0.5574567887828206, -0.015534993721929277, 0.011287179334352642, -0.058563893453183365, 0.013350864610056202, -0.04818672484320412, 0.04327800296616796, 0.14114202015747723]

initial_lte.jl

initial_lte = [0.05472690701449043, 0.3692548806960122, -0.24753874917918603, -0.05000001427918294, 21.26883249027012, 1.9376360856723416, 1.364428314879506, -0.049999596066355685, 60.93788942755487, -0.04995210980031579, 1.0791580265623102, 0.08257963892580181, -2.5125481948555675, -0.05362717322096584, -0.04156533028625943, -0.042382255866398606, -0.050107117571746476, -0.03201558873669222, -51.768232204714636, -0.04999536267299161, -96.41537111287732, -0.05004485070293295, 94.0821418760349, -0.04979042015456535, 61.54318252237016, -0.04170931204366888, -0.04700811996589601, 0.050148228745085684, -30.206098882787778, -0.039679524919648, 46.520297973555664, 34.87907310257046, -0.12327563607325141, -0.04976907875797756, 0.4947908911801988, 1.8906774195314668, -1.3817423458669187]

opt_forcing.jl

# Optimize the parameters to minimize the error criterion
#initial_guess = [-0.025264182634587424, -0.029140067289270374, 0.33335443215989946, 0.6580020559602734, -0.10725608828547605, -1.1847105470519264, 0.13433345333775906, -2.1008202881783085, 0.021125467765861355, 0.4120658850558322, -0.016158641018002635, 0.03236285202717148, 0.8506732807107917, 0.011878524326867833, 2.4180032507037317, -0.020304516265993628, 0.011736152444470115, 0.15985870492287657, -0.015854657416780382, -0.002761708443297949, 0.012780057850313886, 0.14306010895396556, -0.02120487697681831]
#initial_guess = [ -0.07622441209264813, 0.05848405851236095, 0.3154700648067066, 0.6900993037194362, -0.013263359735130646, -1.3195081448882422, 0.017732880476693386, -2.1320749066400753, -0.008743393329515431, 0.3841386393424862, -0.038823183104708024, 0.33246435744947994, 0.8933862535216588, 0.01794191693270748, 2.6220168157993915, -0.0011282294376351103, 0.013716914621117383, 0.16589705547853292, 0.000226579086536897, -0.0008205875136117957, 0.0014706566043327335, 0.14799841798195218, -0.07719460242860718]
#initial_forcing =  [3.3138702285653516, -0.0005082513862306119, -1.1233230985759448, 22.446546519178696, -0.0018808514536766403, -0.0012518794495253703, -0.08026210440826871, 0.05313076196575556, 0.2920241927616061, 0.6901630011277564, -0.012519053465143296, -1.3196114545648108, 0.016732870685059794, -2.1320110899573437, -0.009069252050063285, 0.19613827079408658, -0.04174761284433637, 0.3080150731797063, 0.8934406851223046, 0.016816080489430788, 2.6222303423744426, -0.001287350786093868, 0.04735239959249578, 0.1665924205086596, 0.0002442766282887312, -0.0009420366619222261, 0.0013891058915139448, 0.14808648792629242, -0.08084173012369256]
#initial_forcing =  [6.545043709719594, -0.0022739160901968894, -0.5582716989403155, 12.214713803700972, 2.9599031076309643e-5, 0.002628668845652561, -0.08666761661481144, 0.07187335892297696, 0.3153535736336279, 0.727474010362788, -0.01175799164351175, -1.319438249012814, 0.015594992836841057, -2.2814143222027314, -0.02654239764052265, 0.17578722026856125, -0.03882913506352988, 0.3301325843842333, 0.8530158760458952, 0.015958664952036083, 2.883191405352314, -0.0012359493764822135, 0.0055731112690743665, 0.20229464547324186, 0.00019893005497487252, -0.0008152076247390968, -0.0012926978914918758, 0.16657943369583172, -0.08816609826277817]
#initial_forcing = [3.601115326395929, -0.01918853417117629, -1.0184800385259067, 2.6948508914911233, 4.121147731228225e-5, 0.0026240014336916004, 
# -0.006588475562514759, -0.024653465194151288, 0.300633961089336, 0.7377827340472846, -0.11094314878255129, -1.3234273407769426, 0.13206530520960433, 
# -2.2814069755743644, 0.024131756682530106, -0.06286126009668094, -0.006802017915351862, 0.057862530603353945, 0.8550600405134318, 0.01751701820792908, 
#  2.930411598891345, -0.030327974927040573, 0.013501193155745973, 0.23168868493465336, -0.00642036502064293, -0.002769618153148527, -0.011059271254404495, 
#  0.16652726401807688, 0.1586057719202488,  0.000001, 0.1, 0.000001, 0.1, 0.000001, 0.1]
#initial_forcing = [4.131547161815728, -0.032640724491523485, -0.06372389140101843, 6.399064799077626, 0.003143403581994039, 0.061752042818841515, -0.02538475677338493, -0.0320578848806179, 0.33393490212871996, 0.7383377110781572, -0.10680503986362515, -1.480948239870044, 0.1349511117480185, -2.2819422026648732, 0.0211846202832046, -0.07945837951504822, -0.016243252042947496, 0.031149083716982495, 0.8566777942201425, 0.014853890757343275, 3.208148751798296, -0.01580105829594853, 0.011551449095390067, 0.17038473579807464, -0.014581447747530002, -0.0029989103406074503, -0.012937674366122769, 0.16723577137992315, -0.0017373030358161124, -0.00029283328285881975, 0.1367505729892078, -0.00014172721355304113, 0.012601125677971784, 0.0008200501010441476, 0.033945372043258544]
# initial_forcing =  [5.338397311453116, -0.01971407762179162, -0.10491364363194991, 2.490908601513402, 0.47046822826492063, -0.7723060802389159, -0.10820152188820532, -0.12392258559542828, 0.13101610777279893, 0.7381091913328977, -0.02626345578743629, -1.6470011311282193, 0.05279553793391339, -2.2073480244856425, 0.0093775733099673, -0.02162721657184183, -0.07011962024498493, 0.009385154244105236, 0.8957438595650802, 0.007154424762958634, 3.9034327241567977, -0.005652351831627007, 0.004740717872063474, 0.4239253061627341, -0.01932393000257608, -0.008426374295803235, -0.003990990447879322, 0.2409624533505573, 0.16694944967611428, 0.0010500462456156476, 3.8158797289926527, 0.00013943463398619657, -0.43945813870264083, 0.0001211076387880486, -0.18439380428658714]
#initial_forcing = [17.7044810559976, -0.04032847478014285, -0.05894671887687778, 1.1493691404279447, 0.042592922698210835, 0.0065443972675049, -0.03283329379997429, -0.034388307371705396, 0.2875630134548254, 0.7382710493446801, -0.09768736974821803, -1.4063165755048628, 0.11658223743872476, -2.244501849215369, 0.014854792428515208, -0.07078175024509024, -0.02034475904183493, 0.028442021659134856, 0.9782392856930763, 0.01322101048148294, 3.122323674520081, -0.01931771893094451, 0.009823557166019854, 0.1802852102145696, -0.021677177321619272, -0.010014220316173144, -0.008912065882859681, 0.20387685125298072, 0.001962941554343735, 0.0012154390117943463, 0.10379225990170687, 0.0017382400777123082, 0.061163924019539524, 2.4382437707112215e-5, 0.022682778745075783]
initial_forcing = [17.452261519342578, -0.04014697502663168, -0.058491124069511014, 1.1446873857837834, 0.04922782709737959, 0.008468794118487784, -0.032695644728911656, -0.03427277223462172, 0.2887534816161097, 0.738270918322002, -0.09808846684804419, -1.4063177895231558, 0.11706080721830646, -2.2445010287341227, 0.014913541326633457, -0.07088709743543452, -0.02026164062112959, 0.02857082756534435, 0.978239454746619, 0.013246829878291183, 3.122325276764329, -0.019270263697322286, 0.009865160261669653, 0.1803713160265151, -0.02160582386846471, -0.00997677138545945, -0.008946668474894068, 0.203878206121115, 0.0016902376597589918, 0.0012588053531657054, 0.2099903925180267, 0.0017226539985551955, 0.06397988267431463, 2.5002316827213288e-5, 0.02185189405733515]

include("initial_forcing.jl");
#initial_forcing = [19.027743793182978, -0.0374322509268958, -0.06188211606201738, 0.5388015907658139, 0.07183094512341788, 0.028569415519486864, -0.03205776209313545, -0.02760489955619465, 0.32658769039211644, 0.8701442751366372, -0.09944802821096452, -1.4291509751451843, 0.113868322863351, -1.8887381781469434, 0.03720173549635413, -0.07043161753825136, -0.013674746223632852, -0.00251961797435271, 0.936974913861008, 0.020048362129241683, 3.488234133678723, -0.02338344344447066, -0.0017817839139574352, 0.2039928481144155, -0.022352798814591724, -0.01578743271770451, -0.024080476246585897, 0.21775344202893634, 0.005783123435734556, -0.008696683505101195, 0.20319581865062344, -0.0033184795382876567, 0.07151546349205856, 0.005199936632817127, 0.022379734451424098]

using Optim, CSV, DataFrames
using Distributions, Plots

options = Optim.Options(
    iterations=200,  # maximum number of iterations
    show_trace=true,
    show_every=2
)
options1 = Optim.Options(
    iterations=8000,  # maximum number of iterations
    show_trace=true,
    show_every=100
)


df=CSV.read("dlod_compare1.csv", DataFrame; delim=',')
#df=CSV.read("dLOD3original.dat", DataFrame; delim='\t')

X = Vector{Float64}(df[2:6000, 1])'
Y = Vector{Float64}(df[8:6006, 2])';  # 8 6

Yr = 365.24864; 
Yr = 365.2476154;  # + 365.241237718675-0.0012622739+0.00764
Drac = 2*pi*Yr/27.212220815;
Trop = 2*pi*Yr/27.321661554;
Anom = 2*pi*Yr/27.554549886;
Nodal = Drac-Trop;
Peri = (Trop-Anom)/2.0;

Ds(x, A, Drac, p) = A*sin(2*Drac*(x+p))
Ts(x, A, Trop, p) = A*sin(2*Trop*(x+p))
T2s(x, A, Trop, Drac, p) = A*sin((Trop+Drac)*(x+p))
As(x, A, Anom, p) = A*sin(Anom*(x+p))
Semi(x, A, p) = A*sin(4*pi*(x+p))
Evect(x, A, Trop, Anom, p) = A*sin((2*Trop-Anom-4*pi)*(x+p))
Syn(x, A, Trop, p) = A*sin((2*Trop-4*pi)*(x+p))
Annual(x, A, p) = A*cos(2*pi*(x+p))
N(x, A, Nodal, p) = A*sin(Nodal*(x+p))
N2(x, A, Nodal, p) = A*sin(Nodal/2.0*(x+p))
P(x, A, Peri, p) = A*sin(Peri*(x+p))

function integrate(x)
    y = zeros(length(x))  # Initialize the output vector to zeros
    offset = sum(x)/length(x);
    y[1] = -offset; # x[1]  # Compute the first output sample
    for n = 2:length(x)
        y[n] =x[n-1] + y[n-1] 
    end
    return y
end

YI = integrate(Y);

# Define the model
function model(x, damp1, damp2, initial, Value, bg, trend,
               K1, K2, TA, TP, AA, AP, T2A, T2P, SA, SP, A2A, SYA, SYP, EA, EP, E2A, YA, YP, Q1, Q2, DA, DP, K3,   NA, NP, PA, PP, N2A, N2P) 
   Mm = As(x,AA, Anom, AP);
   Msm = Evect(x, EA, Trop, Anom, EP);
   Mfp = T2s(x + A2A*Mm + E2A*Msm, T2A, Trop, Drac, T2P);
   Mf =  Ts(x + A2A*Mm + E2A*Msm,  TA,  Trop,       TP);
   Mfd = Ds(x + A2A*Mm + E2A*Msm,  DA,  Drac,       DP);
   Msf = Syn(x + Q1*Mf + Q2*Mfp, SYA, Trop, SYP);
   First = Mf + Mfp + Mfd + Semi(x, SA, SP) + Annual(x, YA, YP) + Msm + Msf + N(x, NA, Nodal, NP) +  N2(x, N2A, Nodal, N2P) +  P(x, PA, Peri, PP) +
           As(x + K1*(Mf+Mfp+Mfd) + K2*Mm + K3*Msf, AA, Anom, AP); 
   return (First) 
end

# Define the error criterion as the sum of squared errors
function error_criterion(params)
    model_data = model.(X, params...)
    return sum((model_data - Y).^2)
end

result = optimize(error_criterion, initial_forcing, NelderMead(), options)

initial_forcing = Optim.minimizer(result)
# println("Optimal parameters: ", initial_forcing)

Z=model.(X, initial_forcing...); 
println("CC=", cor(Z',Y'));
display(plot(X[1:600],Y[1:600])); display(plot!(X[1:600],Z[1:600]))

# Define the vectors
# names = ["K1", "K2", "TA", "TP", "AA", "AP", "T2A", "T2P", "SA", "SP", "A2A", "SYA", "SYP", "EA", "EP", "E2A", "YA", "YP", "Q1", "Q2", "DA", "DP", "K3"]

# Zip the vectors and print each pair
#for (name, value) in zip(names, initial_guess)
#    println("$name = $value")
#end

#data = [X' Z']
#writedlm("dlod3julia.dat", data, '\t')

######################################################################################################################################################
Month=1

function delta_semiannual(i, N, M)
  i = i % (2*N)
  if i == M
     return -1
  elseif i == N+M      
     return 1
  else
     return 0
  end
end

#function apply_delta(vec, N, M, Value)
#    return [model(vec[i],  Optim.minimizer(result)...) * Value*delta_semiannual(i, N, M) for i in 1:length(vec)]
#end

function iir_filter(x, b0, a1)
    y = zeros(length(x))  # Initialize the output vector to zeros
    y[1] = b0 * x[1]  # Compute the first output sample
    for n = 2:length(x)
        y[n] = b0 * x[n] + a1 * y[n-1]  # Apply the IIR filter equation
    end
    return y
end


function forced(vec, damp1, damp2, initial, Value, bg, trend, 
    K1, K2, TA, TP, AA, AP, T2A, T2P, SA, SP, A2A, SYA, SYP, EA, EP, E2A, YA, YP, Q1, Q2, DA, DP, K3,  NA, NP, PA, PP, N2A, N2P)
    # Initialize output with zeros
    out_vec = zeros(Float64, length(vec))
    #in_vec = apply_delta(vec, 6, Month, Value, tdelta)
    in_vec = [model(vec[i], damp1, damp2, initial, Value, bg, trend,
              K1, K2, TA, TP, AA, AP, T2A, T2P, SA, SP, A2A, SYA, SYP, EA, EP, E2A, YA, YP, Q1, Q2, DA, DP, K3,  NA, NP, PA, PP, N2A, N2P) * Value*delta_semiannual(i, 6, Month) for i in 1:length(vec)]
#    in_vec = [model(vec[i], damp1, damp2, initial, Value, bg, trend,
#              K1, K2, TA, TP, AA, AP, T2A, T2P, SA, SP, A2A, SYA, SYP, EA, EP, E2A, YA, YP, Q1, Q2, DA, DP, K3,  NA, NP, PA, PP, N2A, N2P) for i in 1:length(vec)]
#    integrated = integrate(in_vec)
    out_vec[1] = initial+in_vec[1]
    for i in 2:length(vec)
        out_vec[i] = out_vec[i-1] + in_vec[i]+0.001*trend; #   + bg*integrated[i] + trend;
    end
    out_vec = iir_filter(out_vec, damp1, damp2)
    return out_vec
#     return in_vec
end


df1=CSV.read("lte_results.csv", DataFrame; delim=',');
df=CSV.read("forcing.dat", DataFrame; delim=',');
N0 = nrow(df);
timeline = Vector{Float64}(df[:, 1]);
ltemodel = Vector{Float64}(df1[:, 2]);
#data     = Vector{Float64}(df[:, 3]);
forcing  = Vector{Float64}(df[:, 4]);


# println("CC=", cor(forcing,F_out))

function forcing_error_criterion(params)
    F_out = forced(timeline, params...); 
    return sum((F_out[12:N0] - forcing[12:N0]).^2)                     # + error_criterion(params)
end

# initial_forcing =   [[1.0639602216003332, -0.027065046139161643, -3.1406298700873267, 8.795958049170562, 0.00038524642975911654, 0.000001]; initial_guess]
# initial_forcing =   [[3.3138702285653516, -0.0005082513862306119, -1.1233230985759448, 22.446546519178696, -0.0018808514536766403, -0.0012518794495253703]; initial_guess]

res = optimize(forcing_error_criterion, initial_forcing, NelderMead(), options1)
initial_forcing = Optim.minimizer(res); println("initial_forcing = ", initial_forcing); F_out = forced(timeline, initial_forcing... );
open("initial_forcing.jl", "w") do f
    println(f, "initial_forcing = ", initial_forcing)
end

display(plot(timeline, forcing)); display(plot!(timeline, F_out)); println("CC=", cor(forcing,F_out))


Z=model.(X, initial_forcing...); println("CC ZY=", cor(Z',Y'))
#display(plot(X[1:600],Y[1:600])); display(plot!(X[1:600],Z[1:600]))


function lte(fv,
    K0, Level, Fund, A1, P1, A2, P2, A3, P3, A4, P4, A5, P5, A6, P6, A7, P7, A8, P8, A9, P9, A10, P10, A11, P11, A12, P12,
    AnnualAmp, AnnualPhase, SemiAmp, SemiPhase, Freq, Amp,  Phase, Freq1, Amp1,  Phase1)
    vec=deepcopy(fv);
    for i in 1:length(fv)
        vec[i] = Level + K0*fv[i] + A1*sin(2.0*pi*Fund*1*fv[i] + P1) + A2*sin(2.0*pi*Fund*2*fv[i] + P2) + A3*sin(2.0*pi*Fund*3*fv[i] + P3) + A4*sin(2.0*pi*Fund*4*fv[i] + P4) +
                                    A5*sin(2.0*pi*Fund*5*fv[i] + P5) + A6*sin(2.0*pi*Fund*6*fv[i] + P6) + A7*sin(2.0*pi*Fund*7*fv[i] + P7) + A8*sin(2.0*pi*Fund*8*fv[i] + P8) +
                                    A9*sin(2.0*pi*Fund*9*fv[i] + P9) + A10*sin(2.0*pi*Fund*10*fv[i] + P10) + A11*sin(2.0*pi*Fund*11*fv[i] + P11) + A12*sin(2.0*pi*Fund*12*fv[i] + P12) +
                    AnnualAmp*cos(2*pi*X[i] + AnnualPhase) + SemiAmp*cos(4*pi*X[i] + SemiPhase)  + Amp*sin(2.0*pi*Freq*fv[i] + Phase)  + Amp1*sin(2.0*pi*Freq1*fv[i] + Phase1);
    end
    return vec
end

function lte_error_criterion(params)
    L_out = lte(F_out, params...); 
    return sum((L_out[661:901] - ltemodel[661:901]).^2)  + sum((L_out[1321:1579] - ltemodel[1321:1579]).^2)
end

initial_lte = [ 0.03543453182,    0.25199904432,  -0.20547232043,    
0.05712064935,    0.15191982489,
0.05380426894,    1.29335625628,
0.05750569478,    1.18127713312,
0.05715109507,    1.81600815878,
0.01788368605,    0.51965320022,
0.25299931155,   -1.69184555805,
0.20624320577,    0.80434802004,
0.11047742120,    1.60023730551,
0.14906865746,   -2.26057785793,
0.15001724486,    1.38472072080,
0.27429256062,    0.05374279754,
0.12185868172,    0.79440025450,
0.10114022823599941, -24.23959170160515, 
-0.06781561164694189, 35.04435367593589,
   1.16981187834,    0.22975540684,   -0.85748371909,
   0.28329025627,    0.09805382039,   -2.43050881486
]

include("initial_lte.jl");

res = optimize(lte_error_criterion, initial_lte, NelderMead(), options1)
initial_lte = Optim.minimizer(res)

L_out = lte(F_out, initial_lte... );
println("initial_lte = ", initial_lte)


display(plot(timeline, ylims=(-8,8), ltemodel)); display(plot!(timeline, L_out))
println("CC=", cor(ltemodel,L_out))
open("initial_lte.jl", "w") do f
    println(f, "initial_lte = ", initial_lte)
end