# file KLEIDER/web/src/pinw/Ausgleichskurve.pm
# Ausgleichskurve stückweise aus Polynomen 3. Grades
# 2020-07-30 Herbert Schiemann <h.schiemann@herbaer.de>

package Herbaer::Ausgleichskurve ;

use utf8;
use POSIX qw(floor);

# sect: [a(0), ..., a(n-1)]    aufsteigende Folge der Zerlegungspunkte
sub new {
   my ($class, $sect) = @_;
   my $n = @$sect - 1;
   my $self = {
      "sect"         => $sect,
      "numintervals" => $n,
      
      # 2n + 2 linear unabhängige Funktionen
      # "Zielvektor"
      "vect"         => [ (0) x ($n + $n + 2) ],
      
      # Einträge {"points" => [[x,y], ...], "weights" => [w, ...]}
      "pointdata"  => [],

      # Punkte mit vorgegebenen Werten [ [x, val], ... ]
      "fixvalues"    => [],

      # Punkte mit vorgegebenen Ableitungen [ [x, derivation], ... ]
      "fixderivs"    => [],
      
      # symmetrische Matrix
      # Zahl der Matrixelemente
      # (2n + 2) (2n + 3) / 2
      # (n + 1) (2n + 3)
      # 2n**2 + 5n + 3
      "symmatrix"    => [ (0) x (($n + $n + 5) * $n + 3)],

      # Koeffizienten 0, ..., 2n + 1
      # [ y(x0), y'(x0), y(x1), y'(x1), ... y(xn), y'(xn) ]
      "coeff" => [],
      
      # 3n + 1 Stützpunkte für kubische Bezier-Kurven
      # [ [x0, y(x0)], ., ., [x1, y(x1)], ... , [xn, y(xn)] ]
      "cbezier" => undef,
      
      # 0: neu angelegt
      # 1: Gewichte sind schon addiert, Zerlegung kann nicht mehr geändert werden
      # 2: Fixwerte gesetzt, Gewichte können nicht mehr hinzugefügt werden
      # 3: Matrix ist invertiert, Ausgleichskurve berechnet
      "status"  => 0,
      "verbose" => 0,
   };
   return bless ($self, $class);
} # new

sub verbose {
   my ($self, $verb) = @_;
   $r = $self -> {"verbose"};
   $self -> {"verbose"} = $verb if defined ($verb);
   return $r;
} # verbose

=for comment
Zweite Ableitung ("Krümmung") gewichten.
$w kann eine positive Zahl oder ein Array positiver Zahlen sein.
=cut
sub add_curvweight {
   my ($self, $w) = @_;
   my $verb = $self -> {"verbose"};
   if ($self -> {"status"} > 1) {
      print STDERR "Ausgleichskurve::add_curvweight: ",
         "Gewichtungen können nicht mehr hinzugefügt werden\n"
         if $verb;
      return $self;
   }
   $self -> {"status"} = 1;
   if (! defined ($w)) {
      $w = [1];
   }
   elsif ( ref ($w) ne "ARRAY" ) {
      $w = [$w];
   }
   my $sect = $self -> {"sect"};
   my $n    = $self -> {"numintervals"};
   my $a    = $self -> {"symmatrix"};
=for comment
      Intervall a(k) .. a(k+1)
      Funktionen
      f(2k)     : (x - a(k+1)) ** 2 * (1 + 2 * (x - a(k))/d) / d ** 2
                  (x - a(k) - d) ** 2 * (1 + 2 * (x - a(k))/d) / d ** 2
                  ((x - a(k)) ** 2 - 2 * d * (x - a(k)) + d ** 2 ) * (1 + 2 * (x - a(k)/d) / d ** 2
                  (z ** 2 - 2 * z + 1) * (1 + 2 * z)  mit z := (x - a(k))/d
                  2 * z ** 3 - 4 * z ** 2 + 2 * z + z ** 2 - 2 * z + 1
                  2 * z ** 3 - 3 * z ** 2 + 1
      f'(2k)    : ( 6 * z ** 2 - 6 * z ) / d
      f''(2k)   : ( 12 * z - 6 ) / d ** 2
      f(2k+1)   : (x - a(k+1)) ** 2 * (x - a(k)) / d ** 2
                  (x - a(k) - d)) ** 2 * (x - a(k)/ d ** 2
                  (z - 1) ** 2 * z * d
                  (z ** 2 - 2 * z + 1) * z * d
                  (z ** 3 - 2 * z ** 2 + z) * d
      f'(2k+1)  : 3 * z ** 2 - 4 * z + 1
      f''(2k+1) : (6 * z - 4) / d
      f(2k+2)   : (x - a(k)) ** 2 * (1 - 2 * (x - a(k+1))/d) / d ** 2
                  z ** 2 * (1 - 2 * (x - a(k) - d)/d)
                  z ** 2 * (1 - 2 * (z - 1))
                  z ** 2 * (1 - 2 * z + 2)
                  -2 * z ** 3 + 3 * z ** 2
      f'(2k+2)  : (-6 * z ** 2 + 6 * z) / d
      f''(2k+2) : (-12 * z + 6) / d ** 2
      f(2k+3)   : (x - a(k)) ** 2 * (x - a(k+1)) / d ** 2
                  z ** 2 * (x - a(k) - d)
                  z ** 2 * (z - 1) * d
                  (z ** 3 - z ** 2) * d
      f'(2k+3)  : 3 * z**2 - 2 * z
      f''(2k+3) : (6 * z - 2) / d
=cut
   my $nw = @$w;
   my $k;
   my $cw;
   my $d;          # Länge eines Intervalls
   my ($d2, $d3);  # Potenzen
   my $r;          # Matrix - Index
   for ($k = 0; $k < $n; ++$k) {
      $cw = $w -> [$k % $nw];
      $d = $sect -> [$k + 1] - $sect -> [$k];
      $d2 = $d * $d;
      $d3 = $d2 * $d;
      
      # A(i,j) = .. + cw * Integral (f''(i)*f''(j))

=for comment      
      i = 2k, j = 2k
      r = 2k * (2k + 1) / 2 + 2k
          2 * k**2 + k + 2k
          2 * k**2 + 3 * k
      f''(i) * f''(j) = ( 12 * z - 6 ) ** 2 / d ** 4
                      = (144 * z**2 - 144 * z + 36) / d ** 4
      Stammfkt          ( 144/3 * z**3 - 72 * z**2 + 36z ) / d ** 3
      Integral          ( 144/3 - 72 + 36) / d ** 3
                        (144/3 - 36) / d ** 3
                        ( (144 - 108)/3 ) / d ** 3
                        12 / d ** 3
=cut
      $r = (2 * $k + 3) * $k;
      $a -> [$r] += $cw * 12 / $d3;
=for comment
      i = 2k, j = 2k+1
      r = j * (j+1) / 2 + i 
          (2k + 1) (2k + 2) / 2 + 2k
          2k**2 + 3k + 1 + 2k
          (2k**2 + 3k) + 2k + 1
          2k**2 + 5k + 1
      f''(i) * f''(j) = ( 12 * z - 6 ) * (6 * z - 4) / d ** 3
                      = (72 * z**2 - 48 * z - 36 * z + 24) / d ** 3
                      = (72 * z**2 - 84 * z + 24) / d ** 3
      Stammfkt          ( 144/3 * z**3 - 42 * z**2 + 24z ) / d ** 2
      Integral          ( 144/3 - 42 + 24 ) / d ** 2
                        ( 144/3 - 18 ) / d ** 2
                        ( (144 - 54) / 3 / d ** 2
                        30 / d ** 2
=cut
      $r += 2 * $k + 1;
      $a -> [$r] += $cw * 30 / $d2;
=for comment
      i = 2k, j = 2k+2
      r = j * (j+1) / 2 + i 
          (2k + 2) (2k + 3) / 2 + 2k
          (k + 1) (2k + 3) + 2k
          2k**2 + 5k + 3 + 2k
          2k**2 + 7k + 3
          2k**2 + 5k + 1 + 2k + 2
      f''(i) * f''(j) = ( 12 * z - 6 ) * (-12 * z + 6) / d ** 4
                      = (-144 * z**2 + 144 * z - 36) / d ** 4
      Stammfkt          ( -144/3 * z**3 + 72 * z**2 - 36 * z) / d ** 3
      Integral          ( -144/3 + 72 - 36 ) / d ** 3
                        ( -48 + 72 - 36 ) / d ** 3
                        - 12 / d ** 3
=cut
      $r += 2 * $k + 2;
      $a -> [$r] -= $cw * 12 / $d3;
=for comment
      i = 2k, j = 2k+3
      r = j * (j+1) / 2 + i 
          (2k + 3) (2k + 4) / 2 + 2k
          (2k + 3) (k + 2) + 2k
          2k**2 + 7k + 6 + 2k
          2k**2 + 9k + 6
          2k**2 + 7k + 3 + 2k + 3
      f''(i) * f''(j) = ( 12 * z - 6 ) (6 * z - 2) / d ** 3
                      = ( 72 * z**2 - 36 * z - 24 * z + 12 ) / d ** 3
                      = ( 72 * z**2 - 60 * z + 12 ) / d ** 3
      Stammfkt          ( 24 * z**3 - 30 * z**2 + 12 * z ) / d ** 2
      Integral          ( 24 - 30 + 12 ) / d ** 2
                        6 / d ** 2
=cut
      $r += 2 * $k + 3;
      $a -> [$r] += $cw * 6 / $d2 ;
=for comment
      i = 2k+1, j = 2k+1
      r = j * (j+1) / 2 + i 
          (2k + 1) (2k + 2) / 2 + 2k + 1
          2k**2 + 3k + 1 + 2k + 1
          2k**2 + 5k + 2
          2k**2 + 9k + 6 - 4k - 4
      f''(i) * f''(j) = (6 * z - 4) ** 2 / d ** 2
                      = (36 * z**2 - 48 * z + 16) / d ** 2
      Stammfkt          (12 * z**3 - 24 * z**2 + 16 * z) / d
      Integral          (12 - 24 + 16) / d
                        4 / d
=cut
      $r -= 4 * ($k + 1);
      $a -> [$r] += $cw * 4 / $d;
=for comment
      i = 2k+1, j = 2k+2
      r = j * (j+1) / 2 + i 
          (2k + 2) (2k + 3) / 2 + 2k + 1
          (k + 1) (2k + 3) + 2k + 1
          2k**2 + 5k + 3 + 2k + 1
          2k**2 + 7k + 4
          2k**2 + 5k + 2 + 2k + 2
      f''(i) * f''(j) = (6 * z - 4) (-12 * z + 6) / d ** 3
                      = (-72 * z**2 + 48 * z + 36 * z - 24) / d ** 3
                      = (-72 * z**2 + 84 * z - 24) / d ** 3
      Stammfkt          (-24 * z**3 + 42 * z**2 - 24 * z) / d ** 2
      Integral          (-24 + 42 - 24) / d ** 2
                        -6 / d ** 2
=cut
      $r += 2 * $k + 2;
      $a -> [$r] -= $cw * 6 / $d2 ;
=for comment
      i = 2k+1, j = 2k+3
      r = j * (j+1) / 2 + i 
          (2k + 3) (2k + 4) / 2 + 2k + 1
          (2k + 3) (k + 2) + 2k + 1
          2k**2 + 7k + 6 + 2k + 1
          2k**2 + 9k + 7
          2k**2 + 7k + 4 + 2k + 3
      f''(i) * f''(j) = (6 * z - 4) (6 * z - 2) / d ** 2
                      = (36 * z**2 - 12 * z - 24 * z + 8) / d ** 2
                      = (36 * z**2 - 36 * z + 8) / d ** 2
      Stammfkt          (12 * z**3 - 18 * z**2 + 8 * z) / d
      Integral          (12 - 18 + 8) / d
                        2 / d
=cut
      $r += 2 * $k + 3;
      $a -> [$r] += $cw * 2 / $d;
=for comment
      i = 2k+2, j = 2k+2
      r = j * (j+1) / 2 + i 
          (2k + 2) (2k + 3) / 2 + 2k + 2
          (k + 1) (2k + 3) + 2k + 2
          2k**2 + 5k + 3 + 2k + 2
          2k**2 + 7k + 5
          2k**2 + 9k + 7 - 2k - 2
      f''(i) * f''(j) = (-12 * z + 6) (-12 * z + 6) / d ** 4
                      = 36 * (2z - 1) (2z - 1) / d ** 4
                      = 36 * (4 * z**2 - 4 * z + 1) / d ** 4
                      = (144 * z**2 - 144 * z + 36) / d ** 4
      Stammfkt          (48 * z ** 3 - 72 * z**2 + 36 * z) / d ** 3
      Integral          (48 - 72 + 36) / d ** 3
                        12 / d ** 3
=cut
      $r -= 2 * $k + 2;
      $a -> [$r] += $cw * 12 / $d3;
=for comment
      i = 2k+2, j = 2k+3
      r = j * (j+1) / 2 + i 
          (2k + 3) (2k + 4) / 2 + 2k + 2
          (2k + 3) (k + 2) +  2k + 2
          2k**2 + 7k + 6 + 2k + 2
          2k**2 + 9k + 8
          2k**2 + 7k + 5 + 2k + 3
      f''(i) * f''(j) = (-12 * z + 6) (6 * z - 2) / d ** 3
                      = (-72 * z**2 - 24 * z + 36 * z - 12) / d ** 3
                      = (-72 * z**2 + 12 * z - 12) / d ** 3
      Stammfkt          (-24 * z**3 +  6 * z**2 - 12 * z) / d ** 2
      Integral          (-24 + 6 - 12) / d ** 2
                        -30 / d ** 2
=cut
      $r += 2 * $k + 3;
      $a -> [$r] -= $cw * 30 / $d2;
=for comment
      i = 2k+3, j = 2k+3
      r = j * (j+1) / 2 + i
          (2k + 3) (2k + 4) / 2 + 2k + 3
      f''(i) * f''(j) = (6 * z - 2) (6 * z - 2) / d ** 2
                      = (36 * z**2 - 24 * z + 4) / d ** 2
      Stammfkt          (12 * z**3 - 12 * z**4 + 4 * z) / d
      Integral          4 / d
=cut
      $r += 1;
      $a -> [$r] += $cw * 4 / $d;
   }
   return $self;
} # add_curvweight

sub get_cbezier {
   my $self = shift;
   my $verb = $self -> {"verbose"};   
   $self -> solve() if $self -> {"status"} < 3;
   if ($self -> {"status"} < 3) {
      print STDERR "Ausgleichskurve::get_cbezier: ",
         "Problem bei der Berechnung\n"
         if $verb;
      return undef;
   }
   if (! $self -> {"cbezier"}) {
      my $c = $self -> {"coeff"};
      my $s = $self -> {"sect"};
      my $n = $self -> {"numintervals"};
      my $b = [ [0, 0] x (3 * $n + 1) ];
      my $i3 = 0;
      my $i2 = 0;
      for ($i = 0; $i < $n;) {
         $x0 = $s->[$i];
         $x1 = $s->[++$i];
         $d = ($x1 - $x0) / 3;
         $b -> [$i3++] = [$x0, $c -> [$i2]];
         $b -> [$i3++] = [$x0 + $d, $c -> [$i2]   + $d * $c -> [$i2 + 1]];
         $b -> [$i3++] = [$x1 - $d, $c -> [$i2+2] - $d * $c -> [$i2 + 3]];
         $i2 += 2;
      }
      $b -> [$i3] = [$s -> [$n], $c -> [2 * $n]];
      $self -> {"cbezier"} = $b;
   }
   return $self -> {"cbezier"};
} # get_cbezier

sub add_points {
   my ($self, $points, $w) = @_;
   my $verb = $self -> {"verbose"};
   print STDERR "Ausgleichskurve::add_points\n" if $verb;
   if ($self -> {"status"} > 1) {
      print STDERR "Ausgleichskurve::add_points: ",
         "Gewichtungen können nicht mehr hinzugefügt werden\n"
         if $verb;
      return $self;
   }
   $self -> {"status"} = 1;
   if (! defined ($w)) {
      $w = [1];
   }
   elsif ( ref ($w) ne "ARRAY" ) {
      $w = [$w];
   }
   elsif ( ! @$w ) {
      $w = [1];
   }
   push (@{$self -> {"pointdata"}},  {"points" => $points, "weights" => $w});
   my $sect = $self -> {"sect"};
   my $n    = $self -> {"numintervals"};
   my $a    = $self -> {"symmatrix"};
   my $b    = $self -> {"vect"};
   my $nw   = @$w;
   my $cw;         # Gewicht
   my ($x0, $x1);  # Intervallgrenzen
   my $xe   = $sect -> [$n];
   my ($x, $y);
   my $d;          # Intervalllänge
   my $z;
   my ($f0, $f1, $f2, $f3);
   
=for comment
Ausgleichskurve Summe (0 <= j < 2*n + 1) c(j) * f(j)

zu minimieren:
        
        w(i)(c(j) * f(j)(x(i)) - y(i)) ** 2 + crvweight * Integral (c(j) * f''(j)) ** 2
      = w(i) * ( (c(j)*c(k) * f(j)(x(i))*f(k)(x(i)) - 2 * c(j)*f(j)(x(i))*y(i) + y(i)**2 )
        + crvweight * c(j)*c(k) * Integral (f''(j)*f''(k))
      = c(j)c(k) * ( w(i)*f(j)(x(i))*f(k)(x(i)) + crvweight * Integral (f''(j)*f''(k)) )
        - 2 * c(j) * ( w(i)*f(j)(x(i))*y(i) )
        + w(i) * y(i) ** 2
      = c(j)c(k) * A(j,k) - 2 * c(j) * B(j) + C

mit
        A(j,k) = w(i) * f(j)(x(i)) * f(k)(x(i)) + crvweight * Integral (f''(j)*f''(k))
        B(j)   = w(i) * f(j)(x(i))*y(i)
        C      = w(i) * y(i) ** 2
        (d(j) * A(j,k) * c(k) + c(j) * A(j,k) * d(k)) - 2 * d(j) * B(j) = 0
        A(j,k) * c(k) - B(j) = 0
=cut
   my $m;             # Punkt-Index
   my $p = @$points;  # Zahl der Punkte
   
=for comment
   links vom Anfangspunkt
   f(0)    = 1
   f(1)    = z = x - a[0]
=cut
   $x1 = $sect -> [0];
   for ($m = 0; $m < $p; ++$m) {
      ($x, $y) = @{$points -> [$m]};
      last if $x > $x1;
      $cw = $w -> [ $m % $nw ];
      $z = $x - $x1;
=for comment
      i = 0, j = 0
      r = j * (j+1) / 2 + i 
          0
      f0 = 1
=cut
      $a -> [0] += $cw ;           # $f0 * $f0 * $cw;
=for comment
      i = 0, j = 1
      r = j * (j+1) / 2 + i 
          1
      f(1) = z
=cut
      $a -> [1] += $z * $cw ;      # $f0 * $f1 * $cw;
=for comment
      i = 1, j = 1
      r = j * (j+1) / 2 + i 
          2
=cut
      $a -> [2] += $z * $z * $cw ;  # $f1 * $f1 * $cw;
      $b -> [0] += $y;
      $b -> [1] += $z * $y;
   }
#  mittlerer Teil
   $k = 0;
   for (; $m < $p; ++$m) {
      ($x, $y) = @{$points -> [$m]};
      last if $x > $xe;
      while ($x > $x1) {
         $x0 = $x1;
         $x1 = $sect -> [++$k];
         $d  = $x1 - $x0;
      }
      $z = ($x - $x0) / $d;
      $z2 = $z * $z;
      $cw = $w -> [ $m % $nw ];
=for comment
      f(2k-2) = 2 * z ** 3 - 3 * z ** 2 + 1
      f(2k-1) = (z ** 3 - 2 * z ** 2 + z) * d
      f(2k)   = -2 * z ** 3 + 3 * z ** 2
      f(2k+1) = (z ** 3 - z ** 2) * d
=cut
      $f0 = ((2 * $z) - 3) * $z2 + 1;
      $f1 = (($z - 2) * $z + 1) * $z * $d;
      $f2 = (-2 * $z + 3) * $z2;
      $f3 = ($z - 1) * $z2 * $d;
=for comment
      A(i,j) = f(i)(x(m)) * f(j)(x(m)) * w(m)
      B(j)   = f(j)(x(m))*y(m)         * w(m)
=cut
=for comment      
      i = 2k-2, j = 2k-2
      r = j * (j+1) / 2 + i 
          (2k - 2) (2k - 1) / 2 + 2k - 2
          (k - 1) (2k - 1) + 2k -2
          2k**2 - 3k + 1 + 2k - 2
          2k**2 - k - 1
=cut
      $r = (2 * $k - 1) * $k - 1;
      $a -> [$r] += $f0 * $f0 * $cw;
=for comment
      i = 2k-2, j = 2k-1
      r = j * (j+1) / 2 + i 
          (2k - 1) * k  + 2k - 2
          2k**2 - k + 2k - 2
          2k**2 + k - 2
          2k**2 - k - 1 + 2k - 1
=cut
      $r += 2 * $k - 1;
      $a -> [$r] += $f0 * $f1 * $cw;
=for comment
      i = 2k-2, j = 2k
      r = j * (j+1) / 2 + i 
          (2k) (2k + 1) / 2 + 2k - 2
          2k**2 + k + 2k - 2
          2k**2 + 3k - 2
          2k**2 + k - 2 + 2k
=cut
      $r += 2 * $k;
      $a -> [$r] += $f0 * $f2 * $cw;
=for comment
      i = 2k-2, j = 2k+1
      r = j * (j+1) / 2 + i 
          (2k + 1) (2k + 2) / 2 + 2k - 2
          (2k + 1) (k + 1) + 2k - 2
          2k**2 + 3k + 1 + 2k - 2
          2k**2 + 5k - 1
          2k**2 + 3k - 2 + 2k + 1
=cut
      $r += 2 * $k + 1;
      $a -> [$r] += $f0 * $f3 * $cw;
=for comment
      i = 2k-1, j = 2k-1
      r = j * (j+1) / 2 + i 
          (2k - 1) (2k) / 2 + 2k - 1
          2k**2 - k + 2k - 1
          2k**2 + k - 1
          2k**2 + 5k - 1 - 4k
=cut
      $r -= 4 * $k;
      $a -> [$r] += $f1 * $f1 * $cw;
=for comment
      i = 2k-1, j = 2k
      r = j * (j+1) / 2 + i 
          (2k) (2k + 1) / 2 + 2k - 1
          2k**2 + k + 2k - 1
          2k**2 + 3k - 1
          2k**2 + k - 1 + 2k
=cut
      $r += 2 * $k;
      $a -> [$r] += $f1 * $f2 * $cw;
=for comment
      i = 2k-1, j = 2k+1
      r = j * (j+1) / 2 + i 
          (2k + 1) (2k + 2) / 2 + 2k - 1
          (2k + 1) (k + 1) + 2k - 1
          2k**2 + 3k + 1 + 2k - 1
          2k**2 + 5k
          2k**2 + 3k - 1 + 2k + 1
=cut
      $r += 2 * $k + 1;
      $a -> [$r] += $f1 * $f3 * $cw;
=for comment
      i = 2k. j = 2k
      r = j * (j+1) / 2 + i 
          (2k) (2k + 1) / 2 + 2k
          2k**2 + k + 2k
          2k**2 + 3k
          2k**2 + 5k - 2k
=cut
      $r -= 2 * $k;
      $a -> [$r] += $f2 * $f2 * $cw;
=for comment
      i = 2k, j = 2k+1
      r = j * (j+1) / 2 + i 
          (2k + 1) (2k + 2) / 2 + 2k
          (2k + 1) (k + 1) + 2k
          2k**2 + 3k + 1 + 2k
          2k**2 + 5k + 1
          2k**2 + 3k + 2k + 1
=cut
      $r += 2 * $k + 1;
      $a -> [$r] += $f2 * $f3 * $cw;
=for comment
      i = 2k+1, j = 2k+1
      r = j * (j+1) / 2 + i
          (2k + 1) (2k + 2) / 2 + 2k + 1
          (2k + 1) (k + 1) + 2k + 1
          2k**2 + 3k + 1 + 2k + 1
          2k**2 + 5k + 2
          2k**2 + 5k + 1 + 1
=cut
      $r += 1;
      $a -> [$r] += $f3 * $f3 * $cw;
      $y *= $cw;
      $r = $k + $k - 2;
      $b -> [$r++] += $f0 * $y;
      $b -> [$r++] += $f1 * $y;
      $b -> [$r++] += $f2 * $y;
      $b -> [$r++] += $f3 * $y;
   } 
=for comment
   rechts vom Endpunkt
   f(2n)   = 1
   f(2n+1) = z = x - xend
=cut
   for (; $m < $p; ++$m) {
      ($x, $y) = @{$points -> [$m]};
      $cw = $w -> [ $m % $nw ];
      $z = $x - $xe;
=for comment
      i = 2n, j = 2n
      r = j * (j+1) / 2 + i 
          n * (2n + 1) + 2n
          2n**2 + 3n
=cut
      $r = ($n + $n + 3) * $n;
      $a -> [$r] += $cw ;           # $f0 * $f0 * $cw;
=for comment
      i = 2n, j = 2n + 1
      r = j * (j+1) / 2 + i 
          (2n + 1) * (n + 1) + 2n
          2n**2 + 3n + 1 + 2n
          2n**2 + 5n + 1
          2n**2 + 3n  + 2n + 1
=cut
      $r += $n + $n + 1;
      $a -> [$r] += $z * $cw ;      # $f0 * $f1 * $cw;
=for comment
      i = 2n + 1, j = 2n + 1
      r = j * (j+1) / 2 + i 
          (2n + 1) * (2n + 2) / 2 + 2n + 1
          2n**2 + 3n + 1 + 2n + 1
          2n**2 + 5n + 1 + 1
=cut
      ++$r;
      $a -> [$r] += $z *$z * $cw ;      # $f1 * $f1 * $cw;
      $y *= $cw;
      $b -> [2 * $n]     += $y;
      $b -> [2 * $n + 1] += $z * $y;
   }
   return $self;
};

sub get_pointdata {
   my $self = shift;
   return $self -> {"pointdata"};
}; # get_pointdata

# vorgegebene feste Werte
sub fix_coeff {
   my ($self, $i, $v) = @_;
   my $verb = $self -> {"verbose"};
   print STDERR "Ausgleichskurve::fix_coeff ($i, $v)\n" if $verb;
   # 0: neu angelegt
   # 1: Gewichte sind schon addiert, Zerlegung kann nicht mehr geändert werden
   # 2: Fixwerte gesetzt, Gewichte können nicht mehr hinzugefügt werden
   # 3: Matrix ist invertiert, Ausgleichskurve berechnet
   my $stat = $self -> {"status"};
   if ($stat == 0) {
      print STDERR "Ausgleichskurve::fix_coeff keine Gewichtungen\n"
         if $verb;
      return $self;
   }
   elsif ($stat > 2) {
      print STDERR
         "Ausgleichskurve::fix_coeff Matrix invertiert\n"
         if $verb;
      return $self;
   }
   $self -> {"status"} = 2;
   my $a = $self -> {"symmatrix"};
   my $b = $self -> {"vect"};
=for comment
       Einschränkung

       c(i) = v  vorgegeben für 0 <= i < 2n 

       c(j)c(k) * A(j,k) - c(j) * B(j)       // 0 <= j,k < 2 * n
       c(j)c(k) * A(j,k) + c(i)c(k) * A(i,k) + c(j)c(i) * A(j,i) + c(i)c(i) * A(i,i) - c(j) * B(j) - c(i) * B(i) 

       c(j)c(k) * A(j,k) + c(k) * c(i)A(k,i) + c(j)c(i) * A(j,i) + c(i)c(i) * A(i,i) - c(j) * B(j) - c(i) * B(i)

       c(j)c(k) * A(j,k) - c(j) * (B(j) - 2 * v * A(j,i)) + ....
=cut
   my $j;
   my $r;
   my $v2 = $v + $v;
   my $n  = @$b;
   my $o = $i * ($i+1) / 2;
   for ($j = 0; $j < $i; ++$j) {
      $r = $o + $j;
      $b -> [$j] -= $v2 * $a -> [$r];
      $a -> [$r] = 0;
   }
   for ($j = $i + 1; $j < $n; ++$j) {
      $r = $j * ($j + 1) / 2 + $i;
      $b -> [$j] -= $v2 * $a -> [$r];
      $a -> [$r] = 0;
   }
   $b -> [$i] = $v;
   $a -> [$i * ($i + 1) / 2 + $i] = 1;
} # fix_coeff

sub fix_value_at_point {
   my ($self, $x, $v) = @_;
   my $s = $self -> {"sect"};
   my $n = $self -> {"numintervals"} + 1;
   push (@{$self -> {"fixvalues"}}, [$x, $v]);
   my $i = 0;
   for ($i = 0; $i < $n; ++$i) {
      if ($s -> [$i] == $x) {
         $self -> fix_coeff ($i + $i, $v);
         last;
      }
   }
   if ($i == $n) {
      print STDERR "Ausgleichskurve::fix_value_at_point ",
         "$x ist kein Zerlegungspunkt\n"
         if $self -> {"verbose"};
   }
   return $self;
} # fix_value_at_point

sub fix_derivation_at_point {
   my ($self, $x, $d) = @_;
   my $s = $self -> {"sect"};
   my $n = $self -> {"numintervals"} + 1;
   push (@{$self -> {"fixderivs"}}, [$x, $d]);
   my $i = 0;
   for ($i = 0; $i < $n; ++$i) {
      if ($s -> [$i] == $x) {
         $self -> fix_coeff ($i + $i + 1, $d);
         last;
      }
   }
   if ($i == $n) {
      print STDERR "Ausgleichskurve::fix_derivative_at_point ",
         "$x ist kein Zerlegungspunkt\n"
         if $self -> {"verbose"};
   }
   return $self;
} # fix_derivation_at_point

sub solve {
   my $self = shift;
   my $verb = $self -> {"verbose"};
   # 0: neu angelegt
   # 1: Gewichte sind schon addiert, Zerlegung kann nicht mehr geändert werden
   # 2: Fixwerte gesetzt, Gewichte können nicht mehr hinzugefügt werden
   # 3: Matrix ist invertiert, Ausgleichskurve berechnet
   my $stat = $self -> {"status"};
   if ($stat == 0) {
      print STDERR "Ausgleichskurve::solve keine Gewichtungen\n"
         if $verb;
      return $self;
   }
   elsif ($stat > 2) {
      print STDERR
         "Ausgleichskurve::solve Matrix invertiert\n"
         if $verb;
      return $self;
   }
   $self -> {"status"} = 3;
   my $a = $self -> {"symmatrix"};
   my $b = $self -> {"vect"};
   my $s = $self -> {"sect"};
   my $n = $self -> {"numintervals"};
   my ($i, $i2, $i3);
   my ($x0, $x1);  # Intervall-Endpunkte
   my $d;          # Intervalllänge / 3
   symmatrix_invert ($a);
   my $c = symmatrix_mult_vector ($a, $b);
   $self -> {"coeff"} = $c;
   return $self;
} # solve

# invertiert eine symmetrische Matrix "inplace"
sub symmatrix_invert {
   my $a = shift;
   my $n2 = @$a;
=for comment
   n2 = n * (n + 1) / 2
   2 * n2 = n**2 + n
   8 * n2 + 1 = (4 * n**2 + 4 * n + 1) = (2 * n + 1) ** 2
   2 * n + 1 = sqrt (8 * n2 + 1)
   n = (sqrt (8 * n2 + 1) - 1) / 2
=cut
   my $n = floor (sqrt (8 * $n2 - 1) / 2 + 0.5);
=for comment
   a(i,j) * x(j) = b(i)
   tausche x(k) und b(k):
       x'(k) = b(k)
       b'(k) = - x(k)
       x'(i) = x(i) für i != k
       b'(i) = b(i)
   a(i,j) * x(j) + a(i,k) * x(k) = b(i)         # Summe über j != k
   a(k,j) * x(j) + a(k,k) * x(k) = b(k)
   x(k) = b(k)/a(k,k) - (a(k,j)/a(k,k)) * x(j)
   a(i,j) * x(j) + a(i,k)/a(k,k) * b(k) - a(i,k) * a(k,j) / a(k,k) * x(j) = b(i)
   (a(i,j) - a(i,k) * a(k,j) / a(k,k)) * x(j)  + a(i,k)/a(k,k) * b(k) = b(i)
            a(k,j)/a(k,k)              * x(j)  -      1/a(k,k) * b(k) = - x(k)
 
   a'(i,j) * x'(j) = b'(i)
   a'(i,j) = a(i,j) - a(i,k) * a(k,j) / a(k,k)
   a'(i,k) = a(i,k) / a(k,k)
   a'(k,j) = a(k,j) / a(k,k)
   a'(k,k) = - 1/a(k,k)
=cut
   my $inv = [ (0) x $n ]; # markiert die getauschten Indizes
   my $c = 0;              # Anzahl der "getauschten" Indizex
   my $d;                  # maximaler Wert auf der "Diagonalen"
   my $i;
   my $j;
   my $k;
   my ($r, $ri, $rj);      # Index in der Matrix
   for ($c = 0; $c < $n; ++ $c) {
      $d = 0;
      $k = 0;
=for comment
   r(i-1) = (i-1) * i / 2 + i - 1
          = 1/2 * i**2 - 1/2 * i + i - 1
          = 1/2 * i**2 + 1/2 * i - 1
   r(i)   = i * (i+1) / 2 + i
          = 1/2 * i**2 + 1/2 * i + i
          = r(i-1) + i + 1
   r(-1)  = (-2) * (-1) / 2 + (-1) - 1
          = 1 - 1 - 1
          = -1
=cut

      $r = -1;
      for ($i = 0; $i < $n; ++$i) {
         $r += $i + 1;
         next if $inv -> [$i];
         if (abs ($d) < abs ($a -> [$r])) {
            $d = $a -> [$r];
            $k = $i;
         }
      }
      $inv -> [$k] = 1;
      if ($d == 0) {
         print STDERR "ERROR: Matrix nicht invertierbar\n";
         last;
      }

# =for comment    
      $rj = $k * ($k + 1) / 2;
      $ri = $rj;
      $r  = 0;
      for ($j = 0; $j < $k; ++$j) {
         for ($i = 0; $i <= $j; ++$i) {
            $a -> [$r] -= $a -> [$ri] * $a -> [$rj] / $d;
            ++$ri;
            ++$r;
         }
         $ri -= $j + 1;
         ++$rj;
      }
      $rj += $k + 1;
      $r  += $k + 1;
      for ($j = $k + 1; $j < $n; ++$j) {
         $ri = $k * ($k + 1) / 2;
         for ($i = 0; $i < $k; ++$i) {
            $a -> [$r] -= $a -> [$ri] * $a -> [$rj] / $d;
            ++$ri;
            ++$r;
         }
         $ri += $k + 1;
         ++$r;
         for ($i = $k + 1; $i <= $j; ++$i) {
            $a -> [$r] -= $a -> [$ri] * $a -> [$rj] / $d;
            $ri += $i + 1;
            ++$r;
         }
         $rj += $j + 1;
      }
      $r = $k * ($k + 1) / 2;
      for ($i = 0; $i < $k; ++$i) {
         $a -> [$r] /= $d;
         ++$r;
      }
      $a -> [$r] = - 1 / $d;
      $r += $k + 1;
      for ($i = $k + 1; $i < $n; ++$i) {
         $a -> [$r] /= $d;
         $r += $i + 1;
      }
# =cut

=for comment
      # zum Verständnis derselbe Inhalt nicht optimiert codiert
      for ($j = 0; $j < $n; ++$j) {
         next if $j == $k;
         $rj = $j < $k ? $k * ($k + 1) / 2 + $j : $j * ($j + 1) / 2 + $k;
         for ($i = 0; $i <= $j; ++$i) {
            next if $i == $k;
            $ri = $i < $k ? $k * ($k + 1) / 2 + $i : $i * ($i + 1) / 2 + $k;
            $r  = $i < $j ? $j * ($j + 1) / 2 + $i : $i * ($i + 1) / 2 + $j;
            $a -> [$r] -= $a -> [$ri] * $a -> [$rj] / $d;
         }
      }
      for ($i = 0; $i < $n; ++$i) {
         $r = $i < $k ? $k * ($k + 1) / 2 + $i : $i * ($i + 1) / 2 + $k;
         if ($i == $k) {
            $a -> [$r] = - 1/$d;
         }
         else {
            $a -> [$r] /= $d;
         }
      }
=cut

   }
   for ($r = 0; $r < $n2; ++$r) {
      $a -> [$r] = - $a -> [$r];
   }
} # symmatrix_invert

# ergibt einen neuen Vektor
sub symmatrix_mult_vector {
   my ($m, $v) = @_;
   my $n = @$v;
   my $p = [ 0 x $n ];
   my $e = $n * ($n + 1) / 2;
   my $i = 0;
   my $j = 0;
   for ($r = 0; $r < $e; ++$r) {
      $p -> [$i] += $m -> [$r] * $v -> [$j];
      if ($i == $j) {
         ++$i;
         $j = 0;
      }
      else {
         $p -> [$j] += $m -> [$r] * $v -> [$i];
         ++$j;
      }
   }
   return $p;
} # symmatrix_mult_vector

# ergibt einen neuen Vektor
sub symmatrix_show {
   my ($a, $h) = @_;
   my $n2 = @$a;
   my $n = floor (sqrt (8 * $n2 - 1) / 2 + 0.5);
   my $i;
   my $j;
   my $r = 0;
   for ($i = 0; $i < $n; ++$i) {
      for ($j = 0; $j <= $i; ++$j) {
         print $h ", " if ($j > 0);
         print $a -> [$r++];
      }
      print "\n";
   }
   return $a;
} # symmatrix_mult_vector

sub show_matrix {
   my ($self, $h) = @_;
   symmatrix_show ($self -> {"symmatrix"}, $h);
   return $self;
} # show_matrix

# mode: 0 y-Werte, 1 [x, y] - Paare
sub _get_list {
   my ($self, $start, $step, $end, $mode) = @_;
   my $verb = $self -> {"verbose"};
   print STDERR "Ausgleichskurve::_get_list\n" if $verb;
   $self -> solve() if $self -> {"status"} < 3;
   if ($self -> {"status"} < 3) {
      print STDERR "Ausgleichskurve::_get_list: ",
         "Problem bei der Berechnung\n"
         if $verb;
      return undef;
   }
   if ($step <= 0) {
      print STDERR "Ausgleichskurve::_get_list: ",
         "Schrittweite $step nicht positiv\n"
         if $verb;
      return undef;
   }
   my $sect = $self -> {"sect"};
   my $n    = $self -> {"numintervals"};
   my $c    = $self -> {"coeff"};
   my $r    = [];  # Ergebnis
   my $k    = -2;  # Koeffizienten-Index
   my $i    = 0;   # Intervall-Index
   my $xe   = $sect -> [$n];
   my ($x, $y);
   my ($x0, $x1);  # Intervallgrenzen
   my $d;          # Intervallänge
   my ($z, $z2);
   $x1 = $sect -> [0];
   for ($x = $start; $x <= $end && $x <= $x1; $x += $step) {
      $z = $x - $x1;
      $y = $c -> [0] + $c -> [1] * $z;
      push (@$r, ($mode ? [$x, $y] : $y));
   }
   for (; $x <= $end && $x <= $xe; $x += $step) {
      while ($x > $x1) {
         $x0 = $x1;
         $x1 = $sect -> [++$i];
         $d  = $x1 - $x0;
         $k += 2;
      }
      $z  = ($x - $x0) / $d;
      $z2 = $z * $z;
=for comment
      f(2k)   = 2 * z ** 3 - 3 * z ** 2 + 1
      f(2k+1) = (z ** 3 - 2 * z ** 2 + z) * d
      f(2k+2) = -2 * z ** 3 + 3 * z ** 2
      f(2k+3) = (z ** 3 - z ** 2) * d
=cut
      $y = $c -> [$k]     * ((2 * $z - 3) * $z2 + 1)
         + $c -> [$k + 1] * (($z - 2) * $z + 1) * $z * $d
         + $c -> [$k + 2] * (3 - 2 * $z) * $z2
         + $c -> [$k + 3] * ($z - 1) * $z2 * $d;
      push (@$r, ($mode ? [$x, $y] : $y));
   }
   $k = 2 * $n;
   for (; $x <= $end; $x += $step) {
      $z = $x - $xe;
      $y = $c -> [$k] + $c -> [$k + 1] * $z;
      push (@$r, ($mode ? [$x, $y] : $y));
   }
   return $r;
} # _get_list

sub get_values {
   my ($self, $start, $step, $end) = @_;
   my $verb = $self -> {"verbose"};
   print STDERR "Ausgleichskurve::get_values\n" if $verb;
   $self -> _get_list ($start, $step, $end, 0);
} # get_values

sub get_polygon {
   my ($self, $start, $step, $end) = @_;
   my $verb = $self -> {"verbose"};
   print STDERR "Ausgleichskurve::get_polygon\n" if $verb;
   $self -> _get_list ($start, $step, $end, 1);
} # get_polygon

1;

# end of file KLEIDER/web/src/pinw/Ausgleichskurve.pm