Merge pull request #326 from VibekeSkytt/AIC

sbriseid · web-flow · commit bbeebb4e28c1 · 2022-04-08T13:31:34.000+02:00
Aic
diff --git a/lrsplines2D/app/PointCloud2LR.C b/lrsplines2D/app/PointCloud2LR.C
@@ -102,6 +102,7 @@ void print_help_text()
   std::cout << "-toldoc: Documentation on file format for tolerance domains. \n";
   std::cout << "-featuredoc: Show feature documentation \n";
   std::cout << "-refstrat: Print parameter information to define refinement strategy. \n";
+  std::cout << "-AIC <filename>: Compute AIC at each iteration level and output to the given file. Time consuming. Default: No AIC computations \n";
   std::cout << "-h or --help : Write this text\n";
 }
 
@@ -231,6 +232,8 @@ int main(int argc, char *argv[])
   double minsize = -1.0;
   int feature_out = -1;
   int verbose = 0;
+  int AIC = 0;
+  char *AIC_file = 0;
   
   int initncoef = 10;
   int distribute_ncoef = 0;
@@ -514,7 +517,14 @@ int main(int argc, char *argv[])
 	  if (stat < 0)
 	    return 1;
 	}
-	
+      else if (arg == "-AIC")
+	{
+	  AIC = 1;
+	  int stat = fetchCharParameter(argc, argv, ki, AIC_file, 
+					nmb_par, par_read);
+	  if (stat < 0)
+	    return 1;
+	}
     }
 
   // Read remaining parameters
@@ -848,6 +858,10 @@ int main(int argc, char *argv[])
       approx->addUpperConstraint(extent[5] + zfac);
       approx->setLocalConstraint(zfac);
     }
+
+  if (AIC)
+    approx->setAIC(true);
+  
 #ifdef DEBUG
   std::cout << "Range: " << extent[1]-extent[0] << ", " << extent[3]-extent[2];
   std::cout << ", " << extent[5]-extent[4] << std::endl;
@@ -930,7 +944,25 @@ int main(int argc, char *argv[])
       regular_cloud.write(ofr);
     }
       
-
+  if (AIC)
+    {
+      // Fetch Mahalanobis distance and log likelyhood
+      vector<double> AICinfo;
+      vector<int> ncoef;
+      approx->getAICInfo(AICinfo, ncoef);
+      std::ofstream ofAIC(AIC_file);
+      ofAIC.precision(10);
+      ofAIC << "For each iteration level the following is reported: " << std::endl;
+      ofAIC << "level, number of coefficients, Mahalanobix distance with students t-distribution, \n";
+      ofAIC << "AIC student t-distribution, Mahalanobix distance normal distribution, AIC normal distribution" << std::endl;
+      for (size_t ki=0; ki<ncoef.size(); ++ki)
+	{
+	  ofAIC << ki << "\t" << ncoef[ki];
+	  for (size_t kj=0; kj<4; ++kj)
+	    ofAIC << "\t" << AICinfo[4*ki+kj];
+	  ofAIC << std::endl;
+	}
+    }
  
   std::ofstream sfout(surffile);     // Surface output stream
 
diff --git a/lrsplines2D/app/loglikelyhood.C b/lrsplines2D/app/loglikelyhood.C
@@ -0,0 +1,77 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no                      
+ * SINTEF ICT, Department of Applied Mathematics,                         
+ * P.O. Box 124 Blindern,                                                 
+ * 0314 Oslo, Norway.                                                     
+ *
+ * This file is part of GoTools.
+ *
+ * GoTools is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version. 
+ *
+ * GoTools 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with GoTools. If not, see
+ * <http://www.gnu.org/licenses/>.
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using GoTools.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the GoTools library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT. 
+ */
+
+#include "GoTools/geometry/GoTools.h"
+#include "GoTools/geometry/FileUtils.h"
+#include "GoTools/lrsplines2D/LogLikelyhood.h"
+#include <iostream>
+#include <fstream>
+
+using std::vector;
+using namespace Go;
+
+int main(int argc, char *argv[])
+{
+
+  if (argc != 2)
+    {
+      std::cout << "Parameters: residuals (x,y,z,r)" << std::endl;
+      return 1;
+    }
+  
+  std::ifstream input(argv[1]);
+  double Tny = 10.0;  // Guess value
+
+  // Read point cloud with residuals
+  vector<double> data;
+  vector<double> extent(8);   // Limits for points in all coordinates
+  int nmb_pts = 0;
+  FileUtils::readTxtPointFile(input, 4, data, nmb_pts, extent);
+
+  std::cout << "Number of points: " << nmb_pts << std::endl;
+  // Extract residuals
+  vector<double> residuals(nmb_pts);
+  for (size_t ki=0; ki<nmb_pts; ++ki)
+    residuals[ki] = data[4*ki+3];
+  std::cout << "Min: " << extent[6] << ", max: " << extent[7] << std::endl;
+
+  double loglh2 = 0.0;
+  double loglh = LogLikelyhood::compute(residuals, Tny, true, loglh2);
+  printf("Loglikelyhood: %7.13f, %7.13f \n",loglh, loglh2);
+}
diff --git a/lrsplines2D/app/loglikelyhoodFixed.C b/lrsplines2D/app/loglikelyhoodFixed.C
@@ -0,0 +1,77 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no                      
+ * SINTEF ICT, Department of Applied Mathematics,                         
+ * P.O. Box 124 Blindern,                                                 
+ * 0314 Oslo, Norway.                                                     
+ *
+ * This file is part of GoTools.
+ *
+ * GoTools is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version. 
+ *
+ * GoTools 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with GoTools. If not, see
+ * <http://www.gnu.org/licenses/>.
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using GoTools.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the GoTools library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT. 
+ */
+
+#include "GoTools/geometry/GoTools.h"
+#include "GoTools/geometry/FileUtils.h"
+#include "GoTools/lrsplines2D/LogLikelyhood.h"
+#include <iostream>
+#include <fstream>
+
+using std::vector;
+using namespace Go;
+
+int main(int argc, char *argv[])
+{
+
+  if (argc != 3)
+    {
+      std::cout << "Parameters: residuals (x,y,z,r), degrees of freedom in T-distribution" << std::endl;
+      return 1;
+    }
+  
+  std::ifstream input(argv[1]);
+  double Tny = atof(argv[2]);
+
+  // Read point cloud with residuals
+  vector<double> data;
+  vector<double> extent(8);   // Limits for points in all coordinates
+  int nmb_pts = 0;
+  FileUtils::readTxtPointFile(input, 4, data, nmb_pts, extent);
+
+  std::cout << "Number of points: " << nmb_pts << std::endl;
+  // Extract residuals
+  vector<double> residuals(nmb_pts);
+  for (size_t ki=0; ki<nmb_pts; ++ki)
+    residuals[ki] = data[4*ki+3];
+  std::cout << "Min: " << extent[6] << ", max: " << extent[7] << std::endl;
+
+  double loglh2 = 0.0;
+  double loglh = LogLikelyhood::compute(residuals, Tny, false, loglh2);
+  printf("Loglikelyhood: %7.13f, %7.13f \n",loglh, loglh2);
+}
diff --git a/lrsplines2D/include/GoTools/lrsplines2D/LRSurfApprox.h b/lrsplines2D/include/GoTools/lrsplines2D/LRSurfApprox.h
@@ -507,6 +507,19 @@ class LRSurfApprox
       write_feature_ = false;
     }
 
+  void setAIC(bool AIC)
+  {
+    compute_AIC_ = AIC;
+  }
+
+  void getAICInfo(std::vector<double>& AIC_info, std::vector<int>& ncoef)
+  {
+    AIC_info.clear();
+    ncoef.clear();
+    AIC_info.insert(AIC_info.end(), AIC_.begin(), AIC_.end());
+    ncoef.insert(ncoef.end(), ncoef_.begin(), ncoef_.end());
+  }
+
 private:
     shared_ptr<LRSplineSurface> srf_;
     shared_ptr<Eval1D3DSurf> evalsrf_;
@@ -593,6 +606,11 @@ class LRSurfApprox
     bool write_feature_;
     int ncell_;
 
+  // AIC
+  bool compute_AIC_;
+  std::vector<double> AIC_;
+  std::vector<int> ncoef_;
+
     void initDefaultParams();
 
     /// Define free and fixed coefficients
diff --git a/lrsplines2D/include/GoTools/lrsplines2D/LogLikelyhood.h b/lrsplines2D/include/GoTools/lrsplines2D/LogLikelyhood.h
@@ -0,0 +1,64 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no                      
+ * SINTEF ICT, Department of Applied Mathematics,                         
+ * P.O. Box 124 Blindern,                                                 
+ * 0314 Oslo, Norway.                                                     
+ *
+ * This file is part of GoTools.
+ *
+ * GoTools is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version. 
+ *
+ * GoTools 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with GoTools. If not, see
+ * <http://www.gnu.org/licenses/>.
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using GoTools.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the GoTools library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT. 
+ */
+
+#ifndef LOGLIKELYHOOD_H
+#define LOGLIKELYHOOD_H
+
+#include <vector>
+#include <iostream>
+
+namespace LogLikelyhood
+{
+  double compute(const std::vector<double>& residual, double indegT,
+		 bool iterateT, double& llh2);
+
+  void getMeanAndVariance(const std::vector<double>& residual,
+			  double& mean, double& variance);
+    
+  void getMeanAndVarianceIterated(const std::vector<double>& residual,
+				  double& degreeT,
+				  double& mean, double& variance,
+				  bool iterateT);
+    
+  void iterateDegreeT(const std::vector<double>& residual,
+		      double mean, double variance,
+		      double& degreeT);
+};
+
+#endif
diff --git a/lrsplines2D/src/LRSurfApprox.C b/lrsplines2D/src/LRSurfApprox.C
diff --git a/lrsplines2D/src/LogLikelyhood.C b/lrsplines2D/src/LogLikelyhood.C
