вторник, 14 марта 2017 г.

Гамма-регрессия онлайн


Расчет реализован на основе написанного на Object Pascal математического модуля DMath (автор: Dr Jean DEBORD, Laboratoire de Pharmacologie, Faculte de Medecine 2 Rue du Docteur Marcland, 87025 Limoges, France).
Очень полезный вид регрессионного анализа. Хорошо подходит для оцифровки «крутых» функций, которые для наглядности принято отображать с логарифмической шкалой по оси ординат.
Это был один из самых нудных рефакторингов, однако в итоге получен вполне работоспособный PHP-вариант расчета.
Искомая зависимость:


Ссылка на страницу с онлайн-расчетом:



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


Исходный код (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'];// "рудимент" для совместимости с остальными расчетами. Передается 1.
      //
      // уборщик "мусора" в 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); // удаляем файл
                  }//if
              } //if
            }//foreach
      }
      //
      $rcount = count($dimdata)/2;
      $Lb = 0; $Ub = $rcount-1;
      $ErrCode = '0';
      $filepng = 'xx';
      $B1 = 0;
      $B2 = 0;
      $B3 = 0;
      $B4 = 0;
      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];
      }
      $xmin = $xx[0];
      $xmax = $xx[0];
      $ymin = $yy[0];
      for ($K=1; $K<=$Ub; $K++){
            if($xmin > $xx[$K]){$xmin = $xx[$K];}
            if($xmax < $xx[$K]){$xmax = $xx[$K];}
            if($ymin > $yy[$K]){$ymin = $yy[$K];}
      }
      $delta = ($xmax - $xmin)/99;
      $WW = $xmin;
      $B[2] = 0.999 * $WW;
      $P = -1;         
      for ($K=0; $K<=$Ub; $K++){
            if($Y[$K] > 0){
                  $P++;
                  $WW = $X[$K] - $B[2];
                  $U[$P][1] = Log($WW);
                  $U[$P][2] = $WW;
                  $Y1[$P] = Log($Y[$K]);
                  $S1[$P] = 1.0 / $Y[$K];
            }
      }
      for ($K=0; $K<=$P; $K++){
            if($S1[$K] <= 0){
                  $ErrCode = 'Квази-сингулярная матрица';
                  Exit;
            }
      }
      for ($K=0; $K<=$P; $K++){
            $W[$K] = 1.0 / ($S1[$K]*$S1[$K]);
      }
      $Nvar = 2;
      for ($I=0; $I<=$Nvar; $I++){
            for ($J=0; $J<=$Nvar; $J++){
                  $V[$I][$J] = 0.0;
            }
            $A[$I] = 0.0;
      }
      for ($K=0; $K<=$P; $K++){
            $V[0][0] = $V[0][0] + $W[$K];
            for ($J=1; $J<=$Nvar; $J++){
                  $V[0][$J] = $V[0][$J] + $W[$K] * $U[$K][$J];
            }
            $A[0] = $A[0] + $W[$K] * $Y1[$K];
      }
      for ($J=1; $J<=$Nvar; $J++){
            $V[$J][0] = $V[0][$J];
      }
      for ($K=0; $K<=$P; $K++){
            for ($I=1; $I<=$Nvar; $I++){
                  $WX = $W[$K] * $U[$K][$I];
                  for ($J=$I; $J<=$Nvar; $J++){
                        $V[$I][$J] = $V[$I][$J] + $WX * $U[$K][$J];
                  }
                  $A[$I] = $A[$I] + $WX * $Y1[$K];
            }
      }
      for ($I=2; $I<=$Nvar; $I++){
            for ($J=1; $J<=($I-1); $J++){
                  $V[$I][$J] = $V[$J][$I];
            }
      }
      for ($I=0; $I<=$Nvar; $I++){
            for ($J=0; $J<=$Nvar; $J++){
                  $AA[$I][$J] = $V[$I][$J];
            }
      }          
      $Det = 1.0;
            $K = 0;
      $ErrCode = '0';        
      while ($K <= $Nvar){
            $Pvt = $AA[$K][$K];         
            $Ik = $K;
            $Jk = $K;
            for ($I=$K; $I<=$Nvar; $I++){
                  for ($J=$K; $J<=$Nvar; $J++){
                        if(abs($AA[$I][$J]) > abs($Pvt)){
                             $Pvt = $AA[$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=0; $J<=$Nvar; $J++){
                        $Temp = $AA[$Ik][$J];
                        $AA[$Ik][$J] = $AA[$K][$J];
                        $AA[$K][$J] = $Temp;                          
                  }
                  $Temp = $A[$Ik];
                  $A[$Ik] = $A[$K];
                  $A[$K] = $Temp;
            }
            if($Jk != $K){
                  for ($I=0; $I<=$Nvar; $I++){
                        $Temp = $AA[$I][$Jk];
                        $AA[$I][$Jk] = $AA[$I][$K];
                        $AA[$I][$K] = $Temp;
                  }
            }
            for ($I=0; $I<=$Nvar; $I++){
                  if($I != $K){
                        $MCol[$I] = $AA[$I][$K];
                  $AA[$I][$K] = 0;
                  } else {
                        $MCol[$I] = 0;
                  $AA[$I][$K] = 1;
                  }
            }
            $T = 1 / $Pvt;
            for ($J=0; $J<=$Nvar; $J++){
                  $AA[$K][$J] = $T * $AA[$K][$J];
            }
            $A[$K] = $T * $A[$K];
            for ($I=0; $I<=$Nvar; $I++){
                  if($I != $K){
                        $T = $MCol[$I];
                        for ($J=0; $J<=$Nvar; $J++){
                             $AA[$I][$J] = $AA[$I][$J] - $T * $AA[$K][$J];
                        }
                        $A[$I] = $A[$I] - $T * $A[$K];
                  }
            }
            $K++;
      }
      for ($I=$Nvar; $I>=0; $I--){
            $Ik = $PCol[$I];
            $Temp = $A[$I];
            $A[$I] = $A[$Ik];
            $A[$Ik] = $Temp;
      }                                 
      $B[1] = exp($A[0]);
      $B[3] = $A[1];
      $B[4] = - 1 / $A[2];
      //         
      if($ErrCode == '0'){
          // Коэффициенты:                                                 
            $B1 = $B[1];
            $B2 = $B[2];
            $B3 = $B[3];
            $B4 = $B[4];
            $ErrCode = '1';             
            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;
                  $DataSet->AddPoint($xdot,"abscise2");
                  $xval = $B[1]*pow(($xdot-$B[2]), $B[3])*exp(-($xdot-$B[2])/$B[4]);
                  $DataSet->AddPoint($xval,"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);
            }
      $B1 = str_replace(".", ",", $B1);
      $B2 = str_replace(".", ",", $B2);
      $B3 = str_replace(".", ",", $B3);
      $B4 = str_replace(".", ",", $B4);
      //
      echo json_encode(array($ErrCode, $B1, $B2, $B3, $B4, $filepng));
}
}
?>