"Fossies" - the Fresh Open Source Software Archive

Member "gama-2.05/tests/gama-local/scripts/compare_xml_adjustment.cpp" (10 May 2019, 19064 Bytes) of package /linux/privat/gama-2.05.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "compare_xml_adjustment.cpp": 2.01_vs_2.02.

    1 /* GNU Gama -- testing adjustment results from different algorithms
    2    Copyright (C) 2012, 2014, 2018  Ales Cepek <cepek@gnu.org>
    3 
    4    This file is part of the GNU Gama C++ library.
    5 
    6    This library is free software; you can redistribute it and/or modify
    7    it under the terms of the GNU General Public License as published by
    8    the Free Software Foundation; either version 3, or (at your option)
    9    any later version.
   10 
   11    This program is distributed in the hope that it will be useful,
   12    but WITHOUT ANY WARRANTY; without even the implied warranty of
   13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   14    GNU General Public License for more details.
   15 
   16    You should have received a copy of the GNU General Public License
   17    along with this program; if not, write to the Free Software Foundation,
   18    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
   19 
   20 #include <iostream>
   21 #include <iomanip>
   22 #include <fstream>
   23 #include <gnu_gama/local/float.h>
   24 #include "compare_xml_adjustment.h"
   25 
   26 using GNU_gama::LocalNetworkAdjustmentResults;
   27 
   28 int compare_xml_adjustment(GNU_gama::LocalNetworkAdjustmentResults* html,
   29                            GNU_gama::LocalNetworkAdjustmentResults* xml,
   30                            double covmat_tol)
   31 {
   32   int result_gp = 0, rcoord = 0, robs = 0;
   33 
   34   // needs to fix trailing white spaces and HTML tags
   35   //if (html->description != xml->description) result_gp = 1;
   36 
   37   // general parameters
   38   {
   39     bool test =
   40       html->coordinates_summary.adjusted.xyz    == xml->coordinates_summary.adjusted.xyz   &&
   41       html->coordinates_summary.adjusted.xy     == xml->coordinates_summary.adjusted.xy    &&
   42       html->coordinates_summary.adjusted.z      == xml->coordinates_summary.adjusted.z     &&
   43       html->coordinates_summary.constrained.xyz == xml->coordinates_summary.constrained.xyz&&
   44       html->coordinates_summary.constrained.xy  == xml->coordinates_summary.constrained.xy &&
   45       html->coordinates_summary.constrained.z   == xml->coordinates_summary.constrained.z  &&
   46       html->coordinates_summary.fixed.xyz       == xml->coordinates_summary.fixed.xyz      &&
   47       html->coordinates_summary.fixed.xy        == xml->coordinates_summary.fixed.xy       &&
   48       html->coordinates_summary.fixed.z         == xml->coordinates_summary.fixed.z;
   49     if (!test) {
   50       std::cout << "         coordinate summary failed\n";
   51       result_gp = 1;
   52     }
   53   }
   54 
   55   {
   56     bool test =
   57       html->observations_summary.distances  ==  xml->observations_summary.distances  &&
   58       html->observations_summary.directions ==  xml->observations_summary.directions &&
   59       html->observations_summary.angles     ==  xml->observations_summary.angles     &&
   60       html->observations_summary.xyz_coords ==  xml->observations_summary.xyz_coords &&
   61       html->observations_summary.h_diffs    ==  xml->observations_summary.h_diffs    &&
   62       html->observations_summary.z_angles   ==  xml->observations_summary.z_angles   &&
   63       html->observations_summary.s_dists    ==  xml->observations_summary.s_dists    &&
   64       html->observations_summary.vectors    ==  xml->observations_summary.vectors;
   65     if (!test) {
   66       std::cout << "         observation summary failed\n";
   67       result_gp = 1;
   68     }
   69   }
   70 
   71   {
   72     bool test =
   73       html->project_equations.equations          ==  xml->project_equations.equations           &&
   74       html->project_equations.unknowns           ==  xml->project_equations.unknowns            &&
   75       html->project_equations.degrees_of_freedom ==  xml->project_equations.degrees_of_freedom  &&
   76       html->project_equations.defect             ==  xml->project_equations.defect; 
   77     if (!test) {
   78       std::cout << "         project equations failed\n";
   79       result_gp = 1;
   80     }
   81 
   82     {
   83       double pvv = std::abs(html->project_equations.sum_of_squares
   84                             - xml->project_equations.sum_of_squares);
   85       double qvv = (html->project_equations.sum_of_squares
   86                     + xml->project_equations.sum_of_squares)/2;
   87       if (qvv > 1e-5) pvv /= qvv;
   88       if (pvv > 1e-5)
   89         {
   90           std::cout << "         sum of squares failed\n";
   91           result_gp = 1;
   92         }
   93     }
   94     {
   95       double mapr = std::abs(html->standard_deviation.apriori
   96                              - xml->standard_deviation.apriori);
   97       double qapr = (html->standard_deviation.apriori
   98                      + xml->standard_deviation.apriori)/2;
   99       if (qapr) mapr /= qapr;
  100       if (mapr > 1e-3)
  101         {
  102           std::cout << "         apriori standard deviation failed\n";
  103           result_gp = 1;
  104         }
  105     }
  106     {
  107       double mapo = std::abs(html->standard_deviation.aposteriori
  108                              - xml->standard_deviation.aposteriori);
  109       double qapo = (html->standard_deviation.aposteriori
  110                      + xml->standard_deviation.aposteriori)/2;
  111       if (qapo) mapo /= qapo;
  112       if (mapo > 1e-3)
  113         {
  114           std::cout << "         aposteriori standard deviation failed\n";
  115           result_gp = 1;
  116         }
  117     }
  118 
  119     if (html->project_equations.connected_network
  120         != xml->project_equations.connected_network)
  121       {
  122         std::cout << "         network connectivity test failed\n";
  123         result_gp = 1;
  124       }
  125 
  126     if (html->standard_deviation.using_aposteriori
  127         != xml->standard_deviation.using_aposteriori)
  128       {
  129         std::cout << "         a priori/a posteriori test failed\n";
  130         result_gp = 1;
  131       }
  132 
  133     {
  134       double mprob = std::abs(html->standard_deviation.probability
  135                              - xml->standard_deviation.probability);
  136       double qprob = (html->standard_deviation.probability
  137                      + xml->standard_deviation.probability)/2;
  138       if (qprob) mprob /= qprob;
  139       if (mprob > 1e-5)
  140         {
  141           std::cout << "         probability failed\n";
  142           result_gp = 1;
  143         }
  144     }
  145 
  146     if (html->standard_deviation.using_aposteriori)
  147       {
  148         if (html->standard_deviation.passed
  149             != xml->standard_deviation.passed)
  150           {
  151             std::cout << "         m0 in (lower, upper) test failed\n";
  152             result_gp = 1;
  153           }
  154         {
  155           double mratio = std::abs(html->standard_deviation.probability
  156                                   - xml->standard_deviation.probability);
  157           double qratio = (html->standard_deviation.probability
  158                           + xml->standard_deviation.probability)/2;
  159           if (qratio) mratio /= qratio;
  160           if (mratio > 1e-5)
  161             {
  162               std::cout << "         ratio m0'/m0 failed\n";
  163               result_gp = 1;
  164             }
  165         }
  166         {
  167           double mlower = std::abs(html->standard_deviation.lower
  168                                   - xml->standard_deviation.lower);
  169           double qlower = (html->standard_deviation.lower
  170                           + xml->standard_deviation.lower)/2;
  171           if (qlower) mlower /= qlower;
  172           if (mlower > 1e-5)
  173             {
  174               std::cout << "         lower limit failed\n";
  175               result_gp = 1;
  176             }
  177         }
  178         {
  179           double mupper = std::abs(html->standard_deviation.upper
  180                                   - xml->standard_deviation.upper);
  181           double qupper = (html->standard_deviation.upper
  182                           + xml->standard_deviation.upper)/2;
  183           if (qupper) mupper /= qupper;
  184           if (mupper > 1e-5)
  185             {
  186               std::cout << "         upper limit failed\n";
  187               result_gp = 1;
  188             }
  189         }
  190         {
  191           double mcscale = std::abs(html->standard_deviation.confidence_scale
  192                                   - xml->standard_deviation.confidence_scale);
  193           double qcscale = (html->standard_deviation.confidence_scale
  194                           + xml->standard_deviation.confidence_scale)/2;
  195           if (qcscale) mcscale /= qcscale;
  196           if (mcscale > 1e-5)
  197             {
  198               std::cout << "         confidence scale failed\n";
  199               result_gp = 1;
  200             }
  201         }
  202       } // if (html->standard_deviation.using_aposteriori)
  203 
  204     if (result_gp == 0) std::cout << "         general parameters"
  205                                   << "                 passed\n";
  206   } // general parameters
  207 
  208 
  209   // fixed coordinates
  210   {
  211     double dfix = 0;
  212     LocalNetworkAdjustmentResults::Point p, q;
  213 
  214     for (size_t i=0; i<html->fixed_points.size(); i++)
  215       {
  216         p = html->fixed_points[i];
  217         q = html ->fixed_points[i];
  218 
  219         if (p.hxy != q.hxy || p.hz != q.hz)
  220           {
  221             rcoord = 1;
  222           }
  223         if (p.hxy)
  224           {
  225             double d;
  226             d = p.x - q.x;
  227             if (std::abs(d) > std::abs(dfix)) dfix = d;
  228             d = p.y - q.y;
  229             if (std::abs(d) > std::abs(dfix)) dfix = d;
  230           }
  231         if (p.hz)
  232           {
  233             double d;
  234             d = p.z - q.z;
  235             if (std::abs(d) > std::abs(dfix)) dfix = d;
  236           }
  237       }
  238 
  239     std::cout << "         fixed points       "
  240               << std::scientific << std::setprecision(3) << std::setw(11)
  241               << dfix << " [m] ";
  242     if (std::abs(dfix) < 1e-5)
  243       std::cout << "passed\n";
  244     else
  245       {
  246         std::cout << "failed\n";
  247         rcoord = 1;
  248       }
  249   } // fixed coordinates
  250 
  251 
  252   // adjusted coordinated
  253   {
  254     double aprdif = 0;
  255     double adjdif = 0;
  256 
  257     if (html->approximate_points.size() != xml->approximate_points.size() )
  258       {
  259         std::cout << "         approximate coordinates dimensions "
  260                   << html->approximate_points.size() << " "
  261                   << xml->approximate_points.size()
  262                   << " failed\n";
  263         rcoord = 1;
  264         return rcoord;
  265       }
  266     if (html->adjusted_points.size() != xml->adjusted_points.size() )
  267       {
  268         std::cout << "         adjusted coordinates dimensions "
  269                   << html->approximate_points.size() << " "
  270                   << xml->approximate_points.size()
  271                   << " failed\n";
  272         rcoord = 1;
  273         return rcoord;
  274       }
  275 
  276 
  277     for (size_t n=0; n<html->adjusted_points.size(); n++)
  278       {
  279         LocalNetworkAdjustmentResults::Point& P=html->adjusted_points[n];
  280         LocalNetworkAdjustmentResults::Point& Q=xml ->adjusted_points[n];
  281 
  282         LocalNetworkAdjustmentResults::Point& p=html->approximate_points[n];
  283         LocalNetworkAdjustmentResults::Point& q=xml ->approximate_points[n];
  284 
  285         if (P.id != Q.id)
  286           {
  287             std::cout << "         matching point ids (index " << n << ") "
  288                       << P.id << " " << Q.id << " failed\n";
  289             rcoord = 1;
  290             return rcoord;
  291           }
  292         // std::cout << "ID   " << P.id   << " " <<  Q.id   << "\n";
  293         // std::cout << "id   " << p.id   << " " <<  q.id   << "\n";
  294         // std::cout << "HXY  " << P.hxy  << " " <<  Q.hxy  << "\n";
  295         // std::cout << "hxy  " << p.hxy  << " " <<  q.hxy  << "\n";
  296         // std::cout << "Z    " << P.hz   << " " <<  Q.hz   << "\n";
  297         // std::cout << "z    " << p.hz   << " " <<  q.hz   << "\n";
  298         // std::cout << "CXY  " << P.cxy  << " " <<  Q.cxy  << "\n";
  299         // std::cout << "CZ   " << P.cz   << " " <<  Q.cz   << "\n";
  300         // std::cout << "INDX " << P.indx << " " <<  Q.indx << "\n";
  301         // std::cout << "INDY " << P.indy << " " <<  Q.indy << "\n";
  302         // std::cout << "INDZ " << P.indz << " " <<  Q.indz << "\n";
  303         if (P.id   != Q.id   || p.id   != q.id   ||
  304             P.hxy  != Q.hxy  || p.hxy  != q.hxy  ||
  305             P.hz   != Q.hz   || p.hz   != q.hz   ||
  306             P.cxy  != Q.cxy  ||
  307             P.cz   != Q.cz   /*||
  308             P.indx != Q.indx ||
  309             P.indy != Q.indy || ?xml err?
  310             P.indz != Q.indz */ )
  311           {
  312             std::cout << "         matching point atttributes failed\n";
  313             rcoord = 1;
  314             return rcoord;
  315           }
  316 
  317         double d;
  318 
  319         if (p.hxy)
  320           {
  321             d = p.x - q.x;
  322             if (std::abs(d) > std::abs(aprdif)) aprdif = d;
  323             d = p.y - q.y;
  324             if (std::abs(d) > std::abs(aprdif)) aprdif = d;
  325           }
  326         if (p.hz)
  327           {
  328             d = p.z - q.z;
  329             if (std::abs(d) > std::abs(aprdif)) aprdif = d;
  330             }
  331         if (P.hxy)
  332           {
  333             d = P.x - Q.x;
  334             if (std::abs(d) > std::abs(adjdif)) adjdif = d;
  335             d = P.y - Q.y;
  336             if (std::abs(d) > std::abs(adjdif)) adjdif = d;
  337           }
  338         if (P.hz)
  339           {
  340             d = P.z - Q.z;
  341             if (std::abs(d) > std::abs(adjdif)) adjdif = d;
  342             }
  343 
  344       }
  345     /* test on approximate coordinates is irelevant
  346     std::cout << "         approx.coordinates "
  347               << std::scientific << std::setprecision(3) << std::setw(11)
  348               << aprdif << " [m] ";
  349     if (std::abs(aprdif) < 1e-5)
  350       std::cout << "passed\n";
  351     else
  352       {
  353         std::cout << "failed\n";
  354         rcoord = 1;
  355       }
  356     */
  357     std::cout << "         adjusted coords.   "
  358               << std::scientific << std::setprecision(3) << std::setw(11)
  359               << adjdif << " [m] ";
  360     if (std::abs(adjdif) < 1e-5)
  361       std::cout << "passed\n";
  362     else
  363       {
  364         std::cout << "failed\n";
  365         rcoord = 1;
  366       }
  367   }// adjusted coordinated
  368 
  369   // original index list
  370   {
  371     int tori = 0;
  372     if (html->original_index.size() != xml->original_index.size()) {
  373       tori = 1;
  374     }
  375     else
  376       for (size_t i=0; i<html->original_index.size(); i++)
  377         if (html->original_index[i] != xml->original_index[i])
  378           {
  379             tori = 1;
  380             break;
  381           }
  382 
  383     if (tori) {
  384       std::cout << "         original index list failed\n";
  385       rcoord = 1;
  386     }
  387   }  // original index list
  388 
  389 
  390   // adjusted orientations
  391   {
  392     double oridif = 0;
  393 
  394     if (html->orientations.size() != xml->orientations.size())
  395       {
  396         std::cout << "         adj. orientations dimension test failed\n";
  397         return 1;
  398       }
  399 
  400     for (size_t i=0; i<html->orientations.size(); i++)
  401       {
  402         double d = html->orientations[i].adj - xml->orientations[i].adj;
  403         if (std::abs(d) > std::abs(oridif)) oridif = d;
  404       }
  405 
  406     std::cout << "         adj. orientations  "
  407               << std::scientific << std::setprecision(3) << std::setw(11)
  408               << oridif << " [g] ";
  409     if (std::abs(oridif) < 1e-5)
  410       std::cout << "passed\n";
  411     else
  412       {
  413         std::cout << "failed\n";
  414         rcoord = 1;
  415       }
  416   } // adjusted orientations
  417 
  418   // covariance band
  419   if (1) {
  420     double dcov = 0;
  421     double dmax = 1/covmat_tol;
  422 
  423     const int dim  = html->cov.dim();
  424     const int band = std::min(html->cov.bandWidth(), xml->cov.bandWidth());
  425 
  426     for (int i=1; i<=dim; i++)
  427       for (int j=0; j<=band; j++)
  428         if (i+j <= dim) {
  429           double d = html->cov(i,i+j);
  430           if (d > dmax) dmax = d;
  431 
  432           d -= xml->cov(i,i+j);
  433           if (std::abs(d) > std::abs(dcov)) dcov = d;
  434         }
  435 
  436     std::cout << "         cov. matrix band   "
  437               << std::scientific << std::setprecision(3) << std::setw(11)
  438               << dcov/dmax << " [%] ";
  439     if (std::abs(dcov)/dmax < covmat_tol)
  440       std::cout << "passed\n";
  441     else
  442       {
  443         std::cout << "failed\n";
  444         rcoord = 1;
  445       }
  446    } // covariance band
  447 
  448 
  449   // observations
  450   {
  451     double dang = 0;
  452     double dlin = 0;
  453     double fpar = 0;
  454     double sres = 0;
  455 
  456     if (html->obslist.size() != xml->obslist.size())
  457       {
  458         std::cout << "            ###  failed obslist dimension "
  459                   << html->obslist.size() << " "
  460                   << xml ->obslist.size() << "\n";
  461         return 1;
  462       }
  463 
  464     for (size_t i=0; i<html->obslist.size(); i++)
  465       {
  466         LocalNetworkAdjustmentResults::Observation H = html->obslist[i];
  467         LocalNetworkAdjustmentResults::Observation X = xml ->obslist[i];
  468         if (H.xml_tag != X.xml_tag)
  469           {
  470             std::cout << "            ###  failed xml_tag obs #"
  471                       << i << " " << H.xml_tag << " " << X.xml_tag << "\n";
  472             return 1;
  473           }
  474         if (H.from != X.from)
  475           {
  476             std::cout << "            ###  failed from obs #"
  477                       << i << " " << H.from << " " << X.from << "\n";
  478             return 1;
  479           }
  480         if (H.to != X.to)
  481           {
  482             std::cout << "            ###  failed to obs #"
  483                       << i << " " << H.to << " " << X.to << "\n";
  484             return 1;
  485           }
  486         if (H.left != X.left)
  487           {
  488             std::cout << "            ###  failed left obs #"
  489                       << i << " " << H.left << " " << X.left << "\n";
  490             return 1;
  491           }
  492         if (H.right != X.right)
  493           {
  494             std::cout << "            ###  failed right obs #"
  495                       << i << " " << H.right << " " << X.right << "\n";
  496             return 1;
  497           }
  498 
  499         double dobs = H.obs - X.obs;
  500         double dadj = H.adj - X.adj;
  501         if (H.xml_tag == "angle" ||
  502             H.xml_tag == "direction" ||
  503             H.xml_tag == "zenith-angle" )
  504           {
  505             dobs = std::asin(std::sin(dobs*G2R))*R2G;
  506             dadj = std::asin(std::sin(dadj*G2R))*R2G;
  507 
  508             if (std::abs(dobs) > std::abs(dang)) dang = dobs;
  509             if (std::abs(dadj) > std::abs(dang)) dang = dadj;
  510           }
  511         else
  512           {
  513             if (std::abs(dobs) > std::abs(dlin)) dlin = dobs;
  514             if (std::abs(dadj) > std::abs(dlin)) dlin = dadj;
  515           }
  516 
  517         double df = H.f - X.f;
  518         if (std::abs(df) > std::abs(fpar)) fpar = df;
  519 
  520         double ds = H.std_residual - X.std_residual;
  521         if (std::abs(ds) > std::abs(sres)) sres = ds;
  522       }
  523 
  524     std::cout << "         angular obs.       "
  525               << std::scientific << std::setprecision(3) << std::setw(11)
  526               << dang << " [g] ";
  527     if (std::abs(dang) < 1e-6)
  528       std::cout << "passed\n";
  529     else
  530       {
  531         std::cout << "failed\n";
  532         rcoord = 1;
  533       }
  534 
  535     std::cout << "         linear observations"
  536               << std::scientific << std::setprecision(3) << std::setw(11)
  537               << dlin << " [m] ";
  538     if (std::abs(dlin) < 1e-5)
  539       std::cout << "passed\n";
  540     else
  541       {
  542         std::cout << "failed\n";
  543         rcoord = 1;
  544       }
  545 
  546     std::cout << "         f% obs. checked    "
  547               << std::scientific << std::setprecision(3) << std::setw(11)
  548               << fpar << " [%] ";
  549     if (std::abs(fpar) < 1e-1)
  550       std::cout << "passed\n";
  551     else
  552       {
  553         std::cout << "failed\n";
  554         rcoord = 1;
  555       }
  556 
  557     std::cout << "         standardized resid."
  558               << std::scientific << std::setprecision(3) << std::setw(11)
  559               << sres << "     ";
  560     if (std::abs(sres) < 1e-1)
  561       std::cout << "passed\n";
  562     else
  563       {
  564         std::cout << "failed\n";
  565         rcoord = 1;
  566       }
  567 
  568   } // observations
  569 
  570   return result_gp + rcoord + robs;
  571 }