- added mask

brew/eval_vecfit.m

@@ -22,21 +22,20 @@
 ## Author: Jens <jens@orion>
 ## Created: 2023-06-10
 
-function eval_vecfit ()
-P0 = [-1 -1];
+P0 = [1 1];
 P1 = [2 4];
 P2 = [4 2];
 
 Pp = [P0; P1; P2];
 
 T = [15 12];
-A = [0 0];
+A = [15 20];
 
 Dn = 1;
 Pa = [0 0 0]';
 
 for i=1:50,
-  [A1, Pa, Dn] = vecfit(Pp,T,A,Pa,0.1.*(1-exp(-Dn)))
+  [A1, Pa, Dn] = vecfit(Pp,[],T,A,Pa,0.1.*(1-exp(-Dn)))
 end
 
-endfunction
+

brew/vecfit.m

@@ -22,10 +22,17 @@
 ## Author: Jens <jens@orion>
 ## Created: 2023-06-10
 
-function [P1, Pa1, D] = vecfit (Pp, Pt, P0, Pa0, step)
+function [P1, Pa1, D] = vecfit (Pp, mask, Pt, P0, Pa0, step)
 
-  # A-Priory state
-  P1 = P0 + sum(Pp.*Pa0);
+  [N_subst, N_elem] = size(Pp);
+
+  # Create default mask
+  if isempty(mask),
+    mask = ones(1,N_elem);
+  end
+
+  # A-Priory state with mask
+  P1 = P0 + sum(Pp.*Pa0.*mask);
 
   # Calc distance to target
   d = (Pt-P1);
@@ -35,30 +42,39 @@ function [P1, Pa1, D] = vecfit (Pp, Pt, P0, Pa0, step)
   # Find best prototype
   best_i = 0;
   best_k = 10000;
+  best_dir = 1;
 
-  for i=1:size(Pp)(1),
+  for i=1:N_subst,
     p = Pp(i,:);
     pn = p/norm(p);
+
     # compare d and p
-    k = norm(pn - dn);
-    if k < best_k,
-      best_k = k;
-      best_i = i;
+    for vecdir=[1,-1],
+      k = norm(vecdir*pn - dn);
+      if k < best_k,
+        best_k = k;
+        best_i = i;
+        best_dir = vecdir;
+      end
     end
   end
 
   # Calc next Pa
   Pa1 = Pa0;
 
+  # Calc P1
+  P1 = P0 + sum(Pp.*Pa0);
+
   if best_i,
     # Get value from best prototype
-    dps = step*Pp(best_i,:);
+    delta = best_dir*step;
+    dps = delta*Pp(best_i,:);
 
     # Update point
     P1 += dps;
 
     # Update Point accu
-    Pa1(best_i) += step;
+    Pa1(best_i) += delta;
   end
 
 endfunction

brew/water_optimizer.m

@@ -45,7 +45,7 @@ M_H20 = create_mol([i_H i_O], [2 1], mol_masses)
 
 m_CaSO4 = sum(M_CaSO4)
 
-mask = create_mol([i_Ca i_Mg i_S i_Cl i_Na], [1 1 1 1 1], ones(1, length(mol_masses)))
+mask = create_mol([i_Ca i_Mg i_S i_Cl i_Na], [1 1 1 1 1], mol_masses) > 0
 
 m_empty = create_mol([], [], mol_masses)
 
@@ -54,21 +54,18 @@ molare_masse_H2O = sum(M_H20)
 
 T = create_mol([i_Ca i_Mg i_Na i_S i_Cl], [0.05 0.01 0.02 0.05 0.10], mol_masses)
 A = create_mol([i_Ca i_Mg i_Na i_S i_Cl], [0.094 0.018 0.004 0.044 0.005], mol_masses)
-A = 4*T
+A = 0*T
 
-M_ve = create_mol([i_Ca i_Mg i_Na i_S i_Cl], [1 1 1 1 1]/1000, mol_masses)
 Pa = [0 0 0 0]';
-Pp = [M_CaSO4; M_CaCl2; M_MgSO4; -M_ve];
-
-Pa = [0]';
-Pp = [-M_ve]
+Pp = [M_CaSO4; M_CaCl2; M_MgSO4; M_NaCl];
 
 Dn = 1;
-while Dn > 0.05,
-  T = T
-  [A1, Pa, Dn] = vecfit(T,A,Pa,Pp, 0.1*Dn)
+while Dn > 0.5,
+  [A1, Pa, Dn] = vecfit(Pp,mask,T,A,Pa, 0.0001.*(1-exp(-Dn)))
 end
-
+T=T
+A1=A1
+Pa=Pa
 endfunction