% RUN_JA_CROSSMODEL_FIT_PAR_local.m
% Parallelisierte Parameter-Transformation (parpool 'local')

%% 1) SETTINGS
close all; clear all;

% Reference model
modelA = @JAModel_Thieltges25;  

% Model to optimize
modelB_list = { ...
    @JAModel_Bergqvist96
};
modelB_names = {'Bergqvist96'};

% Reference parameter
JA_A.Ms    = 800e3;
JA_A.a     = 12e3;
JA_A.k     = 3e3;
JA_A.c     = 0.2;
JA_A.alpha = -0.01;

% Excitation
H_list   = [20e3, 20e3, 20e3, 20e3, 20e3, 20e3, 20e3, 20e3, 20e3, 20e3];  % [A/m]
f        = 1;         %[Hz]
nCycles  = 2;        
nSteps   = 512;        
normalizeB = true;     

% DE-setup
DE.NP     = 60;
DE.F      = 0.8;
DE.CR     = 0.9;
DE.MaxGen = 200;

% DE search space
boundsB = [ ...
    200e3, 1900e3;  % Ms
    5e3,     20e3;  % a
    1e3,      5e3;  % k
    0,        0.5;  % c
   -0.1,      0.1   % alpha
];

% Output/Loging
timestamp   = datestr(now,'yyyymmdd_HHMMSS');
out_matfile = sprintf('JA_CrossModelFit_%s.mat', timestamp);
out_csv     = sprintf('JA_CrossModelFit_summary_%s.csv', timestamp);

doPostPlots = false;  % no plots for parfor

%% using Matlab parallel computing
if isempty(gcp('nocreate'))
    parpool('local');  % start local worker process
end

%% prepare model combinations
nModels = numel(modelB_list);
nH      = numel(H_list);
nComb   = nModels * nH;

% pre-allocate 
results(nComb) = struct( ...
    'modelB_name',      '', ...
    'H_amp',            [], ...
    'JA_A',             [], ...
    'JA_B_best',        [], ...
    'bestErr_RMSE',     [], ...
    'fitness_log',      [], ...
    'pop_log',          [], ...
    'H_A_last',         [], ...
    'B_A_last',         [], ...
    'H_B_last',         [], ...
    'B_B_last',         [], ...
    'nCycles',          nCycles, ...
    'nSteps',           nSteps, ...
    'normalizeB',       normalizeB, ...
    'boundsB',          boundsB, ...
    'DE',               DE ...
);

baseSeed = 1;

%% start loop for optimization
parfor idx = 1:nComb
    rng(baseSeed + idx, 'twister');

    [iModel, iH] = ind2sub([nModels, nH], idx);
    modelB       = modelB_list{iModel};
    modelB_name  = modelB_names{iModel};
    H_amp        = H_list(iH);

    % simulate refecerence B(H) curve
    [H_A, ~, B_A, ~] = modelA(JA_A, H_amp, f, nCycles, nSteps);
    [H_A_last, B_A_last] = func_get_last_cycle(H_A, B_A, nSteps);

    % Objective
    obj = @(p) func_objective_Bcurve_match( ...
                  p, modelB, H_amp, f, nCycles, nSteps, ...
                  H_A_last, B_A_last, normalizeB);

    % DE-optimization
    [bestP, bestErr, fitness_log, pop_log] = func_runDE_5param(obj, boundsB, DE);

    JA_Bbest = struct('Ms',    bestP(1), ...
                      'a',     bestP(2), ...
                      'k',     bestP(3), ...
                      'c',     bestP(4), ...
                      'alpha', bestP(5));

    [H_B, ~, B_B, ~] = modelB(JA_Bbest, H_amp, f, nCycles, nSteps);
    [H_B_last, B_B_last] = func_get_last_cycle(H_B, B_B, nSteps);

    % Preallocate results
    r = struct();
    r.modelB_name    = modelB_name;
    r.H_amp          = H_amp;
    r.JA_A           = JA_A;
    r.JA_B_best      = JA_Bbest;
    r.bestErr_RMSE   = bestErr;
    r.fitness_log    = fitness_log;
    r.pop_log        = pop_log;
    r.H_A_last       = H_A_last;
    r.B_A_last       = B_A_last;
    r.H_B_last       = H_B_last;
    r.B_B_last       = B_B_last;
    r.nCycles        = nCycles;
    r.nSteps         = nSteps;
    r.normalizeB     = normalizeB;
    r.boundsB        = boundsB;
    r.DE             = DE;

    results(idx) = r;
end

%% Summary
sum_model   = strings(nComb,1);
sum_H       = zeros(nComb,1);
sum_RMSE    = zeros(nComb,1);
sum_Ms      = zeros(nComb,1);
sum_a       = zeros(nComb,1);
sum_k       = zeros(nComb,1);
sum_c       = zeros(nComb,1);
sum_alpha   = zeros(nComb,1);

for i = 1:nComb
    sum_model(i) = string(results(i).modelB_name);
    sum_H(i)     = results(i).H_amp;
    sum_RMSE(i)  = results(i).bestErr_RMSE;
    sum_Ms(i)    = results(i).JA_B_best.Ms;
    sum_a(i)     = results(i).JA_B_best.a;
    sum_k(i)     = results(i).JA_B_best.k;
    sum_c(i)     = results(i).JA_B_best.c;
    sum_alpha(i) = results(i).JA_B_best.alpha;
end

%% Dump outputs
Tsummary = table(sum_model, sum_H, sum_RMSE, ...
                 sum_Ms, sum_a, sum_k, sum_c, sum_alpha, ...
                 'VariableNames', {'ModelB','H_amp_Aperm','RMSE', ...
                                   'Ms','a','k','c','alpha'});

save(out_matfile, 'results', 'Tsummary', 'modelA', 'modelB_names', 'H_list', 'JA_A', 'DE', 'boundsB', '-v7.3');
writetable(Tsummary, out_csv);

fprintf('\n== Fertig. Ergebnisse gespeichert in ==\n  %s\n  %s\n', out_matfile, out_csv);

%%Plot results
if doPostPlots
    pick = find(strcmp(Tsummary.ModelB,'Jiles92') & Tsummary.H_amp_Aperm==1e4, 1); % For example for H=10 kA/m
    if ~isempty(pick)
        r = results(pick);
        figure('Name','B(H) letzter Zyklus: Referenz vs. gefittet'); hold on; grid on; box on;
        plot(r.H_A_last, r.B_A_last, 'LineWidth', 1.6);
        plot(r.H_B_last, r.B_B_last, 'LineWidth', 1.6);
        xlabel('H [A/m]'); ylabel('B [T]');
        title(sprintf('Ref=A vs. B=%s @ H_0=%.0f A/m', r.modelB_name, r.H_amp));
        legend('Referenz A','Gefittet B','Location','best');
    end
end
