Files
jens c1c1d60c1e - fixed substance enable
- re-enable banned substance on successful calculation
- refactored names
2023-06-15 08:30:00 +02:00

95 lines
2.2 KiB
Objective-C

## 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} =} vecfit (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn
## Author: Jens <jens@orion>
## Created: 2023-06-10
function [P1, na_1, D] = vecfit (Pp, element_mask, Pt, P0, na_0, step)
[N_subst, N_elem] = size(Pp);
# Create default element_mask
if isempty(element_mask),
element_mask = ones(1,N_elem);
end
# A-Priory state with mask
P1 = P0 + sum(Pp.*na_0.*element_mask);
# Calc distance to target
d = (Pt-P1);
D = norm(d);
dn = d/D;
# Enable all substances
subst_mask = ones(1, N_subst);
while 1,
# Find best prototype
best_i = 0;
best_k = 10000;
best_dir = 1;
for i=1:N_subst,
if subst_mask(i) == 1
p = Pp(i,:);
pn = p/norm(p);
# compare d and p
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
end
# Calc next Pa
na_1 = na_0;
# Calc P1
P1 = P0 + sum(Pp.*na_0);
if best_i,
# Get value from best prototype
delta = best_dir*step;
dps = delta*Pp(best_i,:);
# Check for negative n
accu = na_1(best_i) + delta;
if accu < 0,
subst_mask(best_i) = 0;
continue;
end
# Update point
P1 += dps;
# Update Point accu
na_1(best_i) = accu;
# Re-enable substance
subst_mask = ones(1, N_subst);
break;
end
end
endfunction