## Copyright (C) 2023 Jens
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <https://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {} {@var{retval} =} water_optimizer (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn

## Author: Jens <jens@orion>
## Created: 2023-06-08

function water_optimizer (TEST_SUITE, k_ve, learning_rate)

# Load constants
water_constants()

# Profile from tap water (Bad Toelz)
k_Ca_0_mg   = 94.0;	% mg/l
k_Mg_0_mg   = 18.1;	% mg/l
k_Na_0_mg   =  3.6;	% mg/l
k_Cl_0_mg   =  4.6;	% mg/l
k_SO4_0_mg    = 44.1;	% mg/l
k_HCO3_0_mmol =  5.5;	% mmol/l

if 0
# Profile from tap water (Debug)
k_Ca_0_mg = 100; % mg/l
k_Mg_0_mg = 20; % mg/l
k_Na_0_mg = 5; % mg/l
k_Cl_0_mg = 5; % mg/l
k_SO4_0_mg = 50; % mg/l
end

# Target Profile
k_Ca_1_mg = 100; % mg/l
k_Mg_1_mg = 20; % mg/l
k_Na_1_mg = 40; % mg/l
k_Cl_1_mg = 100; % mg/l
k_SO4_1_mg = 50; % mg/l

# Tap water: Calc Stoffmenge/l [mmol/l]
n_Ca_0 = k_Ca_0_mg/sum(M_EL_Ca);
n_Mg_0 = k_Mg_0_mg/sum(M_EL_Mg);
n_Na_0 = k_Na_0_mg/sum(M_EL_Na);
n_Cl_0 = k_Cl_0_mg/sum(M_EL_Cl);
n_SO4_0 = k_SO4_0_mg/sum(M_EL_SO4);

# Target: Calc Stoffmenge/l [mmol/l]
n_Ca_1 = k_Ca_1_mg/sum(M_EL_Ca);
n_Mg_1 = k_Mg_1_mg/sum(M_EL_Mg);
n_Na_1 = k_Na_1_mg/sum(M_EL_Na);
n_Cl_1 = k_Cl_1_mg/sum(M_EL_Cl);
n_SO4_1 = k_SO4_1_mg/sum(M_EL_SO4);

# Initial profile
m_A0 = k_ve*create_mol([EL_Ca EL_Mg EL_Na EL_Cl EL_SO4], [n_Ca_0 n_Mg_0 n_Na_0 n_Cl_0 n_SO4_0], mol_masses)

# Target profile
m_T =       create_mol([EL_Ca EL_Mg EL_Na EL_Cl EL_SO4], [n_Ca_1 n_Mg_1 n_Na_1 n_Cl_1 n_SO4_1], mol_masses);

if TEST_SUITE == 1,
# Create molar masses of allowed substances, easy
Mp = [M_EL_Ca; M_EL_Mg; M_EL_Na; M_EL_Cl; M_EL_SO4];
end

if TEST_SUITE == 2,
# Create molar masses of allowed substances, difficult
Mp = [M_CaSO4+2*M_H2O; M_CaCl2+2*M_H2O; M_MgSO4+7*M_H2O; M_NaCl; M_HCl];
end

# Init vecfit vars
[N_subst, N_elem] = size(Mp);
# create stoff mengen 'n' [mol]
na = zeros(N_subst, 1);

# ignore oxygen
ignore_mask = m_T > 0;

_dn = [];
CONV_COUNTER_RELOAD = 20;
Dn0 = 1;
conv_counter = CONV_COUNTER_RELOAD;
while true,
  [m_A1, na, Dn1] = vecfit(Mp, ignore_mask, m_T, m_A0, na, learning_rate.*(1-exp(-Dn0)));
  Dc = round_n(Dn1,2);
  if Dn0 == Dc,
    if conv_counter == 0,
      break;
    else
      conv_counter -= 1;
    end
  else
    conv_counter = CONV_COUNTER_RELOAD;
  end
  Dn0 = Dc;
  _dn = [_dn Dn1];
end
plot(1:length(_dn), _dn); grid();
title("Distance"); xlabel("Iteration"); ylabel("Dn")
Dn1
m_elem=sum(Mp.*na);
m_subst_add_g_per_hl = na.*sum(Mp,2)/1000*100
m_T
m_A1

endfunction

function zr = round_n(z, n)
   zr = round(z*10^n)/10^n;
endfunction