воскресенье, 12 марта 2017 г.

Расчет коэффициентов полинома онлайн



UPD 13.05.2017. Приношу свои извинения. В листингах была ошибка при инвертировании матрицы выходных коэффициентов (заключительный этап расчета). Исправлено.
Расчет реализован на основе написанного на Object Pascal математического модуля DMath (автор: Dr Jean DEBORD, Laboratoire de Pharmacologie, Faculte de Medecine 2 Rue du Docteur Marcland, 87025 Limoges, France).
Искомая зависимость:
Ссылка на страницу с онлайн-расчетом:


В реализованном на сайте расчете максимально возможная степень полинома – пятая (в оригинальном модуле DMath допустимы полиномы двенадцатой степени).
Что имеем в результате расчета (по выбору пользователя можно просто вывести в небольшом окне, а можно сохранить вместе с исходными данными в RTF-файл):


Процедура на Object Pascal:
procedure Tprocmod.PolFit;
// Polynomial regression.
// Base: mathematical module Dmath. Author: Dr Jean DEBORD, Laboratoire de Pharmacologie,
// Faculte de Medecine 2 Rue du Docteur Marcland, 87025 Limoges (France).
// Искомая зависимость:
// Y(x) = A0+A1*x+A2*x^2+..+A6*x^6
// Примечание. Программные ограничения:
// - индексы исходных массивов - от 1 до 100;
// - полином шестой степени
// (в принципе допустим расчет коэффициентов полинома до 12 степени).
// Исходные данные:
// main.deg - степень полинома;
// main.xpol[], main.ypol[] - статические массивы с исходными данными.
// Результат расчета:
// main.b0 .. main.b6 - коэффициенты полинома (A0..A6).
var
rcount, I, J, K, Ik, Jk:Integer;
X, Y:Array[0..100] of Single;
BA, BB, PRow, PCol, MCol:Array[0..6] of Single;
V, A:Array[0..6, 0..6] of Single;
XI, Det, Pvt, MachEp, Temp, T:Single;
begin
if (main.deg>1) and (main.deg<7) then
   begin
   for K:=1 to 100 do
       begin
       X[K]:=0; Y[K]:=0;
       end;
   for I:=0 to 6 do
       begin
       for J:=0 to 6 do
           begin
           V[I, J]:=0;
           A[I, J]:=0;
           end;
       BA[I]:=0;
       BB[I]:=0;
       PRow[I]:=0;
       PCol[I]:=0;
       MCol[I]:=0;
       end;
   //
   rcount:=-1;
   for K:=1 to 100 do
       begin
       if (main.xpol[K]<>0) or (main.ypol[K]<>0) then
          begin
          Inc(rcount);
          X[rcount]:=main.xpol[rcount];
          Y[rcount]:=main.ypol[rcount];
          end;
       end;
   //
   if main.deg<=rcount then
      begin
      main.b0:=0;
        main.b1:=0;
        main.b2:=0;
        main.b3:=0;
        main.b4:=0;
        main.b5:=0;
      main.b6:=0;
      //
      Lb:=0;
        Ub:=rcount;
      //
      V[0, 0]:=Ub-Lb+1;
      for K:=Lb to Ub do
          begin
          XI:=X[K];
          BA[0]:=BA[0]+Y[K];
          V[0, 1]:=V[0, 1]+XI;
          BA[1]:=BA[1]+XI*Y[K];
          //
          for I:=2 to main.deg do
              begin
              XI:=XI*X[K];
              V[0, I]:=V[0, I]+XI;
              BA[I]:=BA[I]+XI*Y[K];
              end;
          //
          for I:=1 to main.deg do
              begin
              XI:=XI*X[K];
              V[I, main.deg]:=V[I, main.deg]+XI;
              end;
          end;
      //
      for I:=1 to main.deg do
          begin
          for J:=0 to (main.deg-1) do
              begin
              V[I, J]:=V[I-1, J+1];
              end;
          end;
      for I:=0 to main.deg do
          begin
          for J:=0 to main.deg do
              begin
              A[I, J]:=V[I, J];
              end;
          end;
     //
     Det:=1;
     K:=Lb;
     main.ErrCode:='0';
     while K<=main.deg do
         begin
         Pvt:=A[K, K];
         Ik:=K; Jk:=K;
         for I:=K to main.deg do
             begin
             for J:=K to main.deg do
                 begin
                 if Abs(A[I, J])>Abs(Pvt) then
                    begin
                    Pvt:=A[I, J];
                    Ik:=I; Jk:=J;
                    end;
                 end;
             end;
         //
         PRow[K]:=Ik;
         PCol[K]:=Jk;
         Det:=Det*Pvt;
         if Ik<>K then Det:=-Det;
         if Jk<>K then Det:=-Det;
         MachEp:=2.5E-30;
         if Abs(Pvt)<MachEp then
            begin
            main.ErrCode:='Квази-сингулярная матрица';
            Exit;
            end;
         //
         if Ik<>K then
            begin
            for J:=Lb to main.deg do
                begin
                Temp:=A[Ik, J];
                A[Ik, J]:=A[K, J];
                A[K, J]:=Temp;
                end;
            Temp:=BA[Ik];
            BA[Ik]:=BA[K];
            BA[K]:=Temp;
            end;
        //
        if Jk<>K then
           begin
           for I:=Lb to main.deg do
               begin
               Temp:=A[I, Jk];
               A[I, Jk]:=A[I, K];
               A[I, K]:=Temp;
               end;
           end;
        //
        for I:=Lb to main.deg do
            begin
            if I<>K then
               begin
               MCol[I]:=A[I, K];
               A[I, K]:=0;
               end
               else
               begin
               MCol[I]:=0;
               A[I, K]:=1;
               end;
            end;
         //
         T:=1/Pvt;
         for J:=Lb to main.deg do
             begin
             A[K, J]:=T*A[K, J];
             end;
         //
         BA[K]:=T*BA[K];
         for I:=Lb to main.deg do
             begin
             if I<>K then
                begin
                T:=MCol[I];
                for J:=Lb to main.deg do
                    begin
                    A[I, J]:=A[I, J]-T*A[K, J];
                    end;
                BA[I]:=BA[I]-T*BA[K];
                end;
             end;
         //
         Inc(K);
         end;
     for I := main.deg downto Lb do
          begin
          Ik := PCol[I];
          Temp:=B[I];
          B[I]:=B[Ik];
          B[Ik]:=Temp;
                   end;
     //
     if main.ErrCode='0' then
        begin
        main.b0:=BB[0];
        main.b1:=BB[1];
        if main.deg>1 then main.b2:=BB[2];
        if main.deg>2 then main.b3:=BB[3];
        if main.deg>3 then main.b4:=BB[4];
        if main.deg>4 then main.b5:=BB[5];
        if main.deg>5 then main.b6:=BB[6];
        main.ErrCode:='1';
        end
        else
        ShowMessage('Упс! Неустановленная ошибка расчета!');
      end
      else
      ShowMessage('Упс! Недостаточно данных для расчета!');
   end
   else
   ShowMessage('Упс! Некорректная степень полинома!');
end;


Исходный код PHP:
<?php
// Полиномиальная регрессия
header('Content-Type^ text/html; charset=utf-8');
if ($_SERVER['HTTP_X_REQUESTED_WITH']=='XMLHttpRequest'){
      if ($_POST){
            $dimdata = $_POST['dimdata'];// массив с данными (абцисса и ордината) из таблицы
            $Deg = $_POST['deg'];// степень полинома (от 2 до 5)
            // уборщик "мусора" в temp
            // источник:
            // http://itheap.info/blog/web-master-blog/154.html
            function clearTEMP() {
            $dir = 'temp';
                  $files = scandir($dir); // сканирование файлов в папке и создание списка
            unset($files[0], $files[1]); // удаление из списка (массива) элементов ".",".."
            foreach ($files as $fname) { // цикл по списку файлов
                $filename = $dir .'/'. $fname; // создаем полный путь к файлу
                if (file_exists($filename)) { // проверяем наличии файла, за это время он уже мог удалиться
                    $ftime = filemtime($filename); // определяем время создания/изменения файла
                    $d_time = time() - $ftime; //определяем разницу между текущим временем и временем создания файла
                    if ($d_time > 600) { //условие - если время больше чем [разница в сек] (10мин)
                        @unlink($filename); // удаляем файл
                    }
                }
            }
      }
            $B0 = 0;
            $B1 = 0;
            $B2 = 0;
            $B3 = 0;
            $B4 = 0;
            $B5 = 0;
            $filepng = 'xx'; 
            $rcount = count($dimdata)/2;// количество пар значений абциссы и ординаты
            $Lb = 0;
            $Ub = $rcount-1;
            for ($K=0; $K<=$Ub; $K++){
                  $X[$K] = $dimdata[$K*2];
                  $Y[$K] = $dimdata[$K*2+1];
                  $xx[$K] = $dimdata[$K*2];
                  $yy[$K] = $dimdata[$K*2+1];
            }
            for ($I=0; $I<=$Deg; $I++){
                  for ($J=0; $J<=$Deg; $J++){
                        $V[$I][$J] = 0.0;
                  }
                  $B[$I] = 0.0;
            }
            $V[0][0] = $Ub - $Lb + 1;
            for ($K=$Lb; $K<=$Ub; $K++){
                  $XI = $X[$K];                
            $B[0] = $B[0] + $Y[$K];
            $V[0][1] = $V[0][1] + $XI;
            $B[1] = $B[1] + $XI * $Y[$K];
                  //
                  for ($I=2; $I<=$Deg; $I++){
                        $XI = $XI * $X[$K];
                        $V[0][$I] = $V[0][$I] + $XI;
                        $B[$I] = $B[$I] + $XI * $Y[$K];
                  }
                  for ($I=1; $I<=$Deg; $I++){
                        $XI = $XI * $X[$K];
                        $V[$I][$Deg] = $V[$I][$Deg] + $XI;
                  }
            }
            $D1 = $Deg - 1;
            for ($I=1; $I<=$Deg; $I++){
                  $I1 = $I - 1;
                  for ($J=0; $J<=$D1; $J++){
                        $V[$I][$J] = $V[$I1][$J+1];
                  }
            }    
            for ($I=0; $I<=$Deg; $I++){
                  for ($J=0; $J<=$Deg; $J++){
                        $A[$I][$J] = $V[$I][$J];
                  }
            }          
            $Det = 1.0;
            $K = $Lb;
            $ErrCode = '0';        
            while ($K <= $Deg){
                  $Pvt = $A[$K][$K];           
                  $Ik = $K;
                  $Jk = $K;
                  for ($I=$K; $I<=$Deg; $I++){
                        for ($J=$K; $J<=$Deg; $J++){
                              if(abs($A[$I][$J]) > abs($Pvt)){
                                    $Pvt = $A[$I][$J];
                              $Ik = $I;
                              $Jk = $J;
                              }
                        }
                  }                
                  $PRow[$K] = $Ik;
                  $PCol[$K] = $Jk;
                  $Det = $Det * $Pvt;
                  if($Ik != $K){
                        $Det = - $Det;
                  }
                  if($Jk != $K){
                        $Det = - $Det;
                  }          
                  $MachEp = 2.220446049250313E-16;
                  if(abs($Pvt) < $MachEp){
                        $ErrCode = 'Квази-сингулярная матрица';
                        Exit;
                  }
                  if($Ik != $K){
                        for ($J=$Lb; $J<=$Deg; $J++){
                              $Temp = $A[$Ik][$J];
                              $A[$Ik][$J] = $A[$K][$J];
                              $A[$K][$J] = $Temp;                            
                        }
                        $Temp = $B[$Ik];
                        $B[$Ik] = $B[$K];
                        $B[$K] = $Temp;
                  }
                  if($Jk != $K){
                        for ($I=$Lb; $I<=$Deg; $I++){
                              $Temp = $A[$I][$Jk];
                              $A[$I][$Jk] = $A[$I][$K];
                              $A[$I][$K] = $Temp;
                        }
                  }
                  for ($I=$Lb; $I<=$Deg; $I++){
                        if($I != $K){
                              $MCol[$I] = $A[$I][$K];
                        $A[$I][$K] = 0;
                        } else {
                              $MCol[$I] = 0;
                        $A[$I][$K] = 1;
                        }
                  }
                  $T = 1.0 / $Pvt;
                  for ($J=$Lb; $J<=$Deg; $J++){
                        $A[$K][$J] = $T * $A[$K][$J];
                  }
                  $B[$K] = $T * $B[$K];
                  for ($I=$Lb; $I<=$Deg; $I++){
                        if($I != $K){
                              $T = $MCol[$I];
                              for ($J=$Lb; $J<=$Deg; $J++){
                                    $A[$I][$J] = $A[$I][$J] - $T * $A[$K][$J];
                              }
                              $B[$I] = $B[$I] - $T * $B[$K];
                        }
                  }
                  $K++;
            }
            for ($I=$Deg; $I>=$Lb; $I--){
                  $Ik = $PCol[$I];
                  $Temp = $B[$I];
                  $B[$I] = $B[$Ik];
                  $B[$Ik] = $Temp;
            }           
            if($ErrCode == '0'){                                  
                  $ErrCode = '1';
                  // коэффициенты полинома:
                  $B0 = $BB[0];
                  $B1 = $BB[1];
                  if($Deg > 1){$B2 = $BB[2];}
                  if($Deg > 2){$B3 = $BB[3];}
                  if($Deg > 3){$B4 = $BB[4];}
                  if($Deg > 4){$B5 = $BB[5];}
                  //
                  $xmin = $xx[0];
                  $xmax = $xx[0];
                  for ($K=0; $K<=$Ub; $K++){
                        if($xmin > $xx[$K]){$xmin = $xx[$K];}
                        if($xmax < $xx[$K]){$xmax = $xx[$K];}
                  }
                  $delta = ($xmax - $xmin)/99;
                  //
                  clearTEMP();
            // Построение графика
                  // Утилита PChart. Для использования в php.ini
                  // необходимо подключить php_gd2.dll
                  include("pChart/pData.class");
                  include("pChart/pChart.class");
                  // Dataset definition
                  $DataSet = new pData;
                  for ($K=0; $K<=$Ub; $K++){
                        $DataSet->AddPoint($xx[$K],"abscise");
                        $DataSet->AddPoint($yy[$K],"indata");                      
                  }
                  for ($K=0; $K<=99; $K++){
                        $xdot = $xmin + $delta*$K;
                        $zz0 = $B0;
                        $zz1 = $B1*$xdot;
                        if($Deg > 1){$zz2 = $xdot*$xdot*$B2;} else {$zz2 = 0;}     
                        if($Deg > 2){$zz3 = $xdot*$xdot*$xdot*$B3;} else {$zz3 = 0;}
                        if($Deg > 3){$zz4 = $xdot*$xdot*$xdot*$xdot*$B4;} else {$zz4 = 0;}
                        if($Deg > 4){$zz5 = $xdot*$xdot*$xdot*$xdot*$xdot*$B5;} else {$zz5 = 0;}
                        $zall = $zz0 + $zz1 + $zz2 + $zz3 + $zz4 + $zz5;
                        $DataSet->AddPoint($xdot,"abscise2");
                        $DataSet->AddPoint($zall,"outdata");
                  }
                  $DataSet->AddAllSeries();    
                  $DataSet->SetAbsciseLabelSerie();    
                  $DataSet->SetSerieName("Исходные данные","indata");    
                  $DataSet->SetSerieName("Результат расчета","outdata");
                  $DataSet->SetXAxisName("абцисса");      
                  $DataSet->SetYAxisName("ордината"); 
                  // Initialise the graph
                  $Test = new pChart(600,233);
                  $Test->setColorPalette(0,0,0,200);
                  $Test->setColorPalette(1,200,0,0);    
                  $Test->setFontProperties("Fonts/tahoma.ttf",8);    
                  $Test->setGraphArea(70,30,580,190);    
                  $Test->drawFilledRoundedRectangle(7,7,593,228,1,240,240,240);    
                  $Test->drawRoundedRectangle(5,5,595,230,3,128,57,67);    
                  $Test->drawGraphArea(250,250,250,FALSE); 
                  $Test->drawXYScale($DataSet->GetData(),$DataSet->GetDataDescription(),"indata","abscise",100,100,100,TRUE,0);                    
                  $Test->drawGrid(3,FALSE,128,57,67,0);                      
                  // Draw the 0 line    
                  $Test->setFontProperties("Fonts/tahoma.ttf",6);    
                  $Test->drawTreshold(0,50,50,50,TRUE,TRUE);       
                  // Draw the line graph 
                  $Test->drawXYPlotGraph($DataSet->GetData(),$DataSet->GetDataDescription(),"indata","abscise",0,2,1,230,230,230);        
                  $Test->drawXYGraph($DataSet->GetData(),$DataSet->GetDataDescription(),"outdata","abscise2",1);
                  // Finish the graph    
                  $Test->setFontProperties("Fonts/tahoma.ttf",8);
                  $DataSet->RemoveSerie("abscise");
                  $DataSet->RemoveSerie("abscise2");
                  $Test->drawLegend(75,35,$DataSet->GetDataDescription(),240,240,240);    
                  $Test->setFontProperties("Fonts/tahoma.ttf",10);    
                  $Test->drawTitle(60,22,"Полиномиальная регрессия",10,10,10,585);
                  //
                  $file = substr(md5(microtime() . rand(0, 9999)), 0, 12);
                  //
                  $filepng = 'temp/chart_'.$file.'.'.'png';    
                  $Test->Render($filepng);
                  }
            $B0 = str_replace(".", ",", $B0);
            $B1 = str_replace(".", ",", $B1);
            $B2 = str_replace(".", ",", $B2);
            $B3 = str_replace(".", ",", $B3);
            $B4 = str_replace(".", ",", $B4);
            $B5 = str_replace(".", ",", $B5);
            echo json_encode(array($ErrCode, $B0, $B1, $B2, $B3, $B4, $B5, $filepng));         
      }
}
?>