text:
/*
* LEGAL NOTICE
* This computer software was prepared by Battelle Memorial Institute,
* hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830
* with the Department of Energy ( DOE ). NEITHER THE GOVERNMENT NOR THE
* CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
* LIABILITY FOR THE USE OF THIS SOFTWARE. This notice including this
* sentence must appear on any copies of this computer software.
*
* EXPORT CONTROL
* User agrees that the Software will not be shipped, transferred or
* exported into any country or used in any manner prohibited by the
* United States Export Administration Act or any other applicable
* export laws, restrictions or regulations (collectively the "Export Laws").
* Export of the Software may require some form of license or other
* authority from the U.S. Government, and failure to obtain such
* export control license may result in criminal liability under
* U.S. laws. In addition, if the Software is identified as export controlled
* items under the Export Laws, User represents and warrants that User
* is not a citizen, or otherwise located within, an embargoed nation
* (including without limitation Iran, Syria, Sudan, Cuba, and North Korea)
*     and that User is not otherwise prohibited
* under the Export Laws from receiving the Software.
*
* Copyright 2011 Battelle Memorial Institute.  All Rights Reserved.
* Distributed as open-source under the terms of the Educational Community
* License version 2.0 (ECL 2.0). http://www.opensource.org/licenses/ecl2.php
*
* For further details, see: http://www.globalchange.umd.edu/models/gcam/
*
*/

/*!
* \file logbroyden.cpp
* \ingroup objects
* \brief LogBroyden class (Broyden's method solver) source file
* \author Robert Link
*/

#include "util/base/include/definitions.h"
#include <string>
#include <algorithm>
#include <iomanip>
#include <math.h>
#include <xercesc/dom/DOMNode.hpp>
#include <xercesc/dom/DOMNodeList.hpp>

#include "solution/solvers/include/solver_component.h"
#include "solution/solvers/include/logbroyden.hpp"
#include "solution/util/include/calc_counter.h"
#include "marketplace/include/marketplace.h"
#include "containers/include/world.h"
#include "solution/util/include/solution_info_set.h"
#include "solution/util/include/solution_info.h"
#include "solution/util/include/solver_library.h"
#include "util/base/include/util.h"
#include "util/base/include/configuration.h"
#include "util/logger/include/ilogger.h"
#include "util/base/include/xml_helper.h"
#include "solution/util/include/solution_info_filter_factory.h"
#include "solution/util/include/solvable_nr_solution_info_filter.h"

#include "solution/util/include/functor-subs.hpp"
#include "solution/util/include/linesearch.hpp"
#include "solution/util/include/fdjac.hpp"
#include "solution/util/include/edfun.hpp"
#include "solution/util/include/ublas-helpers.hpp"
#include "util/base/include/fltcmp.hpp"
#include "solution/util/include/jacobian-precondition.hpp"

#if USE_LAPACK
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/bindings/lapack/gesvd.hpp>
#include "solution/util/include/svd_invert_solve.hpp"
#else
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/lu.hpp>
#endif

#include "util/base/include/timer.h"

using namespace xercesc;

std::string LogBroyden::SOLVER_NAME = "broyden-solver-component";

#if USE_LAPACK
#define UBMATRIX boost::numeric::ublas::matrix<double,boost::numeric::ublas::column_major>
#else
#define UBMATRIX boost::numeric::ublas::matrix<double>
#endif
#define UBVECTOR boost::numeric::ublas::vector<double>

namespace {
  // helper functions for the std::transform algorithm
  inline double SI2lgprice (const SolutionInfo &si) {
    double p = std::max(si.getPrice(), util::getTinyNumber());
    return log( p );
  }
  inline double SI2price (const SolutionInfo &si) {return si.getPrice();}

  // read-only accessor for solutionInfoSet (used to prepare log outputs)
  const SolutionInfoSet *cSolInfo=0;

  // utility function for finding the minimimum and maximum absolute
  // value entries in a vector.
  void locate_vector_minmax(const UBVECTOR &v, double &vmax, double &vmin, int &imax, int &imin)
  {
    vmax = vmin = fabs(v[0]);
    imax = imin = 0;
    for(int i=1; i<v.size(); ++i) {
      double vabs = fabs(v[i]);
      if(vabs < vmin) {
        vmin = vabs;
        imin = i;
      }
      if(vabs > vmax) {
        vmax = vabs;
        imax = i;
      }
    }
  }
}

int LogBroyden::mLastPer = 0;
int LogBroyden::mPerIter = 0;

bool LogBroyden::XMLParse( const DOMNode* aNode ) {
    // assume we were passed a valid node.
    assert( aNode );

    // get the children of the node.
    DOMNodeList* nodeList = aNode->getChildNodes();

    // loop through the children
    for ( unsigned int i = 0; i < nodeList->getLength(); ++i ){
        DOMNode* curr = nodeList->item( i );
        std::string nodeName = XMLHelper<std::string>::safeTranscode( curr->getNodeName() );

        if( nodeName == "#text" ) {
            continue;
        }
        else if( nodeName == "max-iterations" ) {
            mMaxIter = XMLHelper<unsigned int>::getValue( curr );
        }
        else if( nodeName == "ftol" ) {
            mFTOL = XMLHelper<double>::getValue( curr );
        }
        else if( nodeName == "solution-info-filter" ) {
            mSolutionInfoFilter.reset(
                                      SolutionInfoFilterFactory::createSolutionInfoFilterFromString( XMLHelper<std::string>::getValue( curr ) ) );
        }
        else if(nodeName == "linear-price") {
          mLogPricep = false;
        }
        else if(nodeName == "log-price") {
          mLogPricep = true;    // not strictly necessary, as this is the default.
        }
        else if( SolutionInfoFilterFactory::hasSolutionInfoFilter( nodeName ) ) {
            mSolutionInfoFilter.reset( SolutionInfoFilterFactory::createAndParseSolutionInfoFilter( nodeName, curr ) );
        }
        else {
            ILogger& mainLog = ILogger::getLogger( "main_log" );
            mainLog.setLevel( ILogger::WARNING );
            mainLog << "Unrecognized text string: " << nodeName << " found while parsing "
                    << getXMLName() << "." << std::endl;
        }
    }
    return true;
}

/*! \brief Broyden's method solver.
 * \details Attempts to solve the selected markets using Broyden's
 * method (see Numerical Recipes sectn. 9.7).  Broyden's method is
 * broadly similar to the Newton-Raphson Method.  At each iteration we
 * solve B . dx = -F to get the solution step for that iteration, we
 * backtrack as necessary to ensure that F . F decreases at each step,
 * and we iterate to convergence.  The difference is that in place of
 * the exact Jacobian, J, we use an approximate Jacobian, B.  At each
 * iteration we update B using Broyden's secant condition: B(i+1) =
 * B( i ) + ((dF( i ) - B( i ) . dx( i )) X dx( i )) / (dx . dx).  Since this
 * update does not require any model evaluations, it is *much* faster
 * than computing finite-difference Jacobians.
 *
 * We still need an initial approximation to the Jacobian.  Whenever
 * possible we'll get that from estimates of derivatives that we've
 * previously squirreled away in the SolutionInfo objects.  When we
 * add a new market to the solution set, we may not have good
 * derivatives for the newcomer.  In that case we'll do finite
 * difference approximations for the new column, and we'll zero the
 * off-diagonal terms of the new row.
 *
 * The solver can run in either log-log mode or linear-linear mode.
 *
 * \author Robert Link
 * \param solnset An initial set of SolutionInfo objects representing all of the markets we will attempt to solve
 * \param period Model time period
 * \return Status code indicating whether the algorithm was successful or not.
 */
SolverComponent::ReturnCode LogBroyden::solve(SolutionInfoSet &solnset, int period) {
  ReturnCode code = SolverComponent::ORIGINAL_STATE;

    // If all markets are solved, then return with success code.
    if( solnset.isAllSolved() ){
        return code = SolverComponent::SUCCESS;
    }

    startMethod();
    if(period != mLastPer) {
        // reset our internal counters
        mPerIter = 0;
        mLastPer = period;
    }

    // Update the solution vector for the correct markets to solve.
    // Need to update solvable status before starting solution (Ignore return code)
    solnset.updateSolvable( mSolutionInfoFilter.get() );

    ILogger& solverLog = ILogger::getLogger( "solver_log" );
    solverLog.setLevel( ILogger::NOTICE );
    solverLog << "Beginning Broyden solution for period " << period
              << "Solving " << solnset.getNumSolvable() << "markets.\n";
    if( mLogPricep ) {
      solverLog << "Log price in effect\n";
    }
    else {
      solverLog << "Linear price in effect\n";
    }

    ILogger& worstMarketLog = ILogger::getLogger( "worst_market_log" );
    worstMarketLog.setLevel( ILogger::DEBUG );
    ILogger& singleLog = ILogger::getLogger( "single_market_log" );
    singleLog.setLevel( ILogger::DEBUG );

    size_t nsolv = solnset.getNumSolvable();
    if( nsolv == 0 ){
      solverLog << "No markets were assigned to this solver.  Exiting." << std::endl;
        return SUCCESS;
    }

    solverLog << "Initial market state:\nmkt    \tprice   \tsupply  \tdemand\n";
    std::vector<SolutionInfo> solvables = solnset.getSolvableSet();
    for(size_t i=0; i<solvables.size(); ++i) {
        solverLog << std::setw( 8 ) << i << "\t"
                  << std::setw( 8 ) << solvables[i].getPrice() << "\t"
                  << std::setw( 8 ) << solvables[i].getSupply() << "\t"
                  << std::setw( 8 ) << solvables[i].getDemand()
                  << "\t\t" << solvables[i].getName() << "\n";
    }

    Timer& solverTimer = TimerRegistry::getInstance().getTimer( TimerRegistry::SOLVER );
    solverTimer.start();

    UBVECTOR x( nsolv ), fx( nsolv );
    int neval = 0;

    // set our initial x from the solutionInfoSet
    std::vector<SolutionInfo> smkts(solnset.getSolvableSet());
    if( mLogPricep ) {
      std::transform(smkts.begin(), smkts.end(), x.begin(), SI2lgprice);
    }
    else {
      std::transform(smkts.begin(), smkts.end(), x.begin(), SI2price);
    }

    // This is the closure that will evaluate the ED function
    LogEDFun F(solnset, world, marketplace, period, mLogPricep);
    // check the assumptions:  narg==nrtn==nsolv
    if(F.narg() != nsolv || F.nrtn() != nsolv) {
      solverLog.setLevel(ILogger::SEVERE);
      solverLog << "size mismatch in logbroyden:  nsolv= " << F.narg()
                << "  nrtn= " << F.nrtn()
                << "  nsolv= " << nsolv
                << std::endl;
      abort();
    }

    // scale the initial guess for use in the solver algorithm
    F.scaleInitInputs( x );

    // Call F( x ), store the result in fx
    F(x,fx);

    solverLog.setLevel(ILogger::DEBUG);
    solverLog << "Initial guess:\n" << x << "\nInitial F( x ):\n" << fx << "\n";
    solnset.printMarketInfo("Broyden-initial", calcCounter->getPeriodCount(), singleLog);

    // Precondition the x values to avoid singular columns in the Jacobian
    solverLog.setLevel(ILogger::DEBUG);
    UBMATRIX J(F.narg(), F.nrtn());
    fdjac(F, x, fx, J, true);

    solverLog << ">>>> Main loop jacobian called.\n";
    int pcfail = jacobian_precondition(x, fx, J, F, &solverLog, mLogPricep);

    if( pcfail ) {
      solverLog.setLevel(ILogger::WARNING);
      solverLog << "Unable to find nonsingular initial guess for one or more markets.  bsolve() will probably fail.\n";
      solverLog.setLevel(ILogger::DEBUG);
    }
    else {
      solverLog << "Revised guess:\n" << x << "\nRevised F( x ):\n" << fx << "\n";
    }
    solnset.printMarketInfo("Broyden-preconditioned", calcCounter->getPeriodCount(), singleLog);
    cSolInfo = &solnset;        // make available for log outputs

    // call the solver
    int bstatus = bsolve(F, x, fx, J, neval);
    mPerIter++;                 // increment the iteration count.  This should produce a visible gap in the trace plots.

    solverTimer.stop();

    solverLog.setLevel(ILogger::NOTICE);
    solverLog << "Broyden solver:  neval= " << neval << "\nResult:  ";
    if(bstatus == 0) {
        solverLog << "Broyden solution success.\n";
        code = SUCCESS;
    }
    else if(bstatus == -1) {
        code = FAILURE_ITER_MAX_REACHED;
        solverLog << "Broyden solution failed: Iteration max reached.\n";
    }
    else if(bstatus == -3) {
        code = FAILURE_ZERO_GRADIENT;
        solverLog << "Broyden solution failed:  Encountered zero gradient in F*F.\n";
    }
    else if(bstatus == -4) {
        code = FAILURE_POOR_PROGRESS;
        solverLog << "Broyden solution failed:  repeated poor progress.\n";
    }
    else if(bstatus > 0) {
        code = FAILURE_SINGULAR_MATRIX;
        int singrow = bstatus-1; // L-U decomp returns row number as 1..N numbering
        solverLog << "Broyden solution failed:  Encountered singular matrix (row " << singrow << ").\n";
        solverLog << smkts[singrow] << std::endl;
    }
    else {
        code = FAILURE_UNKNOWN;
        solverLog << "Broyden solution failed for unknown reason.\n";
    }
    if(!solnset.isAllSolved()) {
        solverLog << "The following markets were not solved:\n";
        solnset.printUnsolved( solverLog );
    }

    solverLog << std::endl;

    // log some final debugging info
    const SolutionInfo* maxred = solnset.getWorstSolutionInfo();
    addIteration(maxred->getName(), maxred->getRelativeED());
    if( mLogPricep ) {
      worstMarketLog << "###Broyden-end-logPrice:  " << *maxred << std::endl;
    }
    else {
      worstMarketLog << "###Broyden-end-linearPrice:  " << *maxred << std::endl;
    }

    solnset.printMarketInfo("Broyden-end ", calcCounter->getPeriodCount(), singleLog);
    singleLog << std::endl;

    return code;
}

int LogBroyden::bsolve(VecFVec<double,double> &F, UBVECTOR &x, UBVECTOR &fx,
                       UBMATRIX & B, int &neval)
{
#if !USE_LAPACK
  using boost::numeric::ublas::permutation_matrix;
  using boost::numeric::ublas::lu_factorize;
  using boost::numeric::ublas::lu_substitute;
#endif
  using boost::numeric::ublas::axpy_prod;
  using boost::numeric::ublas::inner_prod;
  int nrow = B.size1(), ncol = B.size2();
  int ageB = 0;   // number of iterations since the last reset on B
  // svd decomposition elements (note nrow == ncol)
#if USE_LAPACK
  UBMATRIX Usv(nrow,ncol),VTsv(ncol,ncol);
  UBVECTOR Ssv( ncol );
#else
  permutation_matrix<int> p(F.narg()); // permutation vector for pivoting in L-U decomposition
#endif

  UBMATRIX Btmp(nrow, ncol);
  ILogger &solverLog = ILogger::getLogger("solver_log");
  ILogger& worstMarketLog = ILogger::getLogger( "worst_market_log" );
  worstMarketLog.setLevel( ILogger::DEBUG );
  ILogger& singleLog = ILogger::getLogger( "single_market_log" );
  singleLog.setLevel( ILogger::DEBUG );

  // market ids and solvable flag for just the solvable markets
  std::vector<int> mktids_solv;
  cSolInfo->getMarketIDs(mktids_solv,true);
  std::vector<bool> issolvable_solv(mktids_solv.size(), true);
  // market ids and solvable flag for all markets (solvable and unsolvable)
  std::vector<int> mktids_all;
  cSolInfo->getMarketIDs(mktids_all, false);
  std::vector<bool> issolvable_all(mktids_all.size());
  for(unsigned i=0; i<issolvable_solv.size(); ++i) { // solvable markets are at the beginning; set their flag to true
      issolvable_all[i] = true;
  }
  for(unsigned i=issolvable_solv.size(); i<issolvable_all.size(); ++i) {
      // unsolvable markets at the end of the array; set their flag to false
      issolvable_all[i] = false;
  }
  // working space for variables that will be printed for solvable and unsolvable markets
  UBVECTOR rptvec_all(mktids_all.size());

  F(x,fx);

  solverLog.setLevel(ILogger::DEBUG);

  neval += 1 + x.size();        // initial function evaluation + jacobian calculations

  const double FTINY = mFTOL*mFTOL;

  UBVECTOR dx(F.narg());
  UBVECTOR xnew(F.narg());
  UBVECTOR gx(F.narg());
  assert(F.nrtn() == F.narg());
  assert(x.size() == F.narg());
  assert(fx.size() == F.nrtn());
  UBVECTOR jdiag(F.narg());
  // do the asserts manually, since we can't turn on normal assertions
  // (GCAM is riddled with asserts that fail under normal
  // circumstances).
  if(F.nrtn() != F.narg() || x.size() != F.narg() || fx.size() != F.nrtn()) {
    solverLog.setLevel(ILogger::SEVERE);
    solverLog << "size mismatch:  nrtn= " << F.nrtn()
              << "  narg= " << F.narg()
              << "  x.size= " << x.size()
              << "  fx.size= " << fx.size()
              << std::endl;
    abort();
  }

  // We create a functor that computes f( x ) = F( x )*F( x ).  It also
  // stores the value of F that it produces as an intermediate.
  FdotF<double,double> fnorm( F );
  double f0 = inner_prod(fx,fx); // already have a value of F on input, so no need to call fnorm yet
  if(f0 < FTINY) {
    // Guard against F=0 since it can cause a NaN in our solver.  This
    // is a more stringent test than our regular convergence test
    return 0;
  }

  bool lsfail = false;        // flag indicating whether we have had a line-search failure
  for(int iter=0; iter<mMaxIter; ++iter) {
    // log some debug info

    solverLog << "Broyden iter= " << iter << "\tneval= " << neval << "\n";
    solverLog << "Internal iteration count ( mPerIter )= " << mPerIter << "\n";
    cSolInfo->printMarketInfo("Broyden ", calcCounter->getPeriodCount(), singleLog);
    for(int j=0;j<F.narg();++j) {
      // double bjj= B(j,j);
      // jdiag[j] = bjj;
      jdiag[j] = B(j,j);
    }
    static_cast<LogEDFun&>(F).setSlope(jdiag);
    double jdmax=0.0, jdmin=0.0;
    int jdjmax=0, jdjmin=0;
    locate_vector_minmax(jdiag, jdmax, jdmin, jdjmax, jdjmin);
    solverLog << "diag( B ):\n" << jdiag << "\n";
    solverLog << "maxval= " << jdmax << " jmax= " << jdjmax << "  "
              << "minval= " << jdmin << "  jmin= " << jdjmin << "\n";

    axpy_prod(fx,B,gx);         // compute the gradient of F*F (= fx^T * B == B^T * fx)
                                // axpy_prod clears gx on entry, so we don't have to do it.
                                // NB: the order of fx and B in that last call is significant!

    // Check for zero gradient.  This indicates a local minimum in f,
    // from which we are unlikely to escape.  We will need to try
    // again with a different initial guess.  (This should be very
    // uncommon)
    double gmag2 = inner_prod(gx,gx);
    if(gmag2 / (f0+FTINY) < mFTOL) {
      solverLog << "**** ||gx|| = 0.  Returning.\n";
      return -3;
    }

    Btmp = B;                   // save the jacobian approximant
#if USE_LAPACK /* Solve using SVD */
    int ierr = boost::numeric::bindings::lapack::gesvd('O','A','A', // control parameters
                                                       B,           // input matrix
                                                       Ssv,Usv,VTsv); // outputs
    if(ierr>0) {
      // svd failed.  It's not even clear under what circumstances
      // this can happen
      solverLog.setLevel(ILogger::SEVERE);
      solverLog << "****************SVD failed.  This shouldn't happen.  It can't mean anything good.\n";
      return ierr;
    }

    // At this point, U, S, and VT contain the SVD of the original Jacobian
    solverLog.setLevel(ILogger::DEBUG);
    dx = -1.0*fx;
    int nsing = svdInvertSolve(Usv,Ssv,VTsv,dx, solverLog);

    solverLog << "\nIteration " << iter << "\nf0= " << f0
              << "\tnsing= " << nsing
              << "\nx: " << x << "\nF( x ): " << fx << "\ndx: " << dx << "\n";

#else /* No USE_LAPACK.  Solve using L-U decomposition */
    int itrial = 0;
    /* If the L-U decomposition fails the first time around, we will
       invoke the jacobian preconditioner and try again.  If it fails
       a second time, we bail out */
    do {
      for(size_t i=0; i<p.size(); ++i) {
        p[i] = i;
      }
      int sing = lu_factorize(B,p);
      if(sing>0) {
        int fail=1;
        B = Btmp;           // restore Jacobian
        if(itrial == 0) {
            solverLog << "Salvaging Jacobian.\n";
            fail = jacobian_precondition(x, fx, B, F, &solverLog, mLogPricep);
            f0 = inner_prod(fx,fx);

            // log the diagonal of the new jacobian
            for(int j=0; j<F.narg(); ++j) {
                jdiag[j] = B(j,j);
            }
            solverLog << "After jacobian salvage.  diag( B )=\n" << jdiag << "\n";

        }

        if( fail ) {
            solverLog.setLevel(ILogger::WARNING);
            solverLog << "Singular Jacobian:\n" << B << "\n";
            return sing;
        }
      }
      else {
        // L-U decomp was successful.  Continue with the next phase of the algorithm.
        break;
      }
    } while(++itrial < 2);

    // J now holds the L-U decomposition of the Jacobian.  Attempt backsubstitution
    dx = -1.0*fx;
    try {
      lu_substitute(B,p,dx);    // solve dx = J^-1 F
    }
    catch (const boost::numeric::ublas::internal_logic &err) {
      // This error seems to be thrown when the Jacobian is
      // ill-conditioned.  We let it go because often the solver will
      // muddle through to a solution.  If not, then it will
      // eventually stop with a genuinely singular matrix.
    }
    solverLog << "dx: " << dx << "\n";
#endif /* USE_LAPACK */

    // log the proposal step
    solverLog << "Proposal step magnitude dxmag= " << sqrt(inner_prod(dx,dx)) << "\n\n";
    reportVec("dxprop", dx, mktids_solv, issolvable_solv);

    // dx now holds the newton step.  Execute the line search along
    // that direction.
    double fnew;
    int lserr = linesearch(fnorm,x,f0,gx,dx, xnew,fnew, neval, &solverLog);

    if(lserr != 0) {
      // line search failed.  There are a couple of things that could
      // be happening here.

      // 1) B is not a descent direction.  If this is the first
      // failure starting from this x value, try a finite difference
      // jacobian
      if(!lsfail) {
        solverLog << "**Failed line search. Evaluating fdjac\n";
        lsfail = true;
        // call fdjac such that it re-calculates the model at x as linesearch will
        // have left off on some other price vector thus we could have bad state
        // data from which we calculate derivatives
        fdjac(F,x,B);
        neval += x.size();
        ageB = 0;  // reset the age on B

        // Log the diagonal of the new jacobian after the failed line search
        for(int j=0; j<F.narg(); ++j) {
            jdiag[j] = B(j,j);
        }
        solverLog << "New Jacobian: diag( B )=\n" << jdiag << "\n";
        static_cast<LogEDFun&>(F).setSlope(jdiag);

        // start the next iteration *without* updating x
        continue;
      }

      // 2) The descent only continues for a very short distance
      // (roughly TOL*x0).  Maybe we're really close to a solution and
      // generating very tiny dx values.  Make a relaxed convergence
      // test and return if we have a "close enough" solution.
      double msf = f0/fx.size();
      if(msf < mFTOL) {
        // basically, we're letting ourselves converge to the sqrt of
        // our intended tolerance.
        return 0;
      }

      // 3) Neither of the above.  The most likely way for this to
      // happen is if we're close to a discontinuity in (probably
      // several components of) the Jacobian.  Using smoother supply
      // and/or demand functions might help here.

      // We're not close enough to the solution, and we don't have a
      // good descent direction.  There are no good options at this point
      // so kick out and hope that the preconditioner can set us straight.
      solverLog << "linesearch failure\n";
      return -4;
    }
    else {
      // reset the line search fail flag
      lsfail = false;
    }

    UBVECTOR xstep(xnew-x);    // step in x eventually taken
    double lambda = fabs(dx[0]) > 0.0 ? xstep[0] / dx[0] : 0.0;
    solverLog << "################Return from linesearch\nfold= " << f0 << "\tfnew= " << fnew
              << "\tlambda= " << lambda << "\n";

    UBVECTOR fxnew(fx.size());
    fnorm.lastF( fxnew );            // get the last value of big-F
    solverLog << "\nxnew: " << xnew << "\nfxnew: " << fxnew << "\n";
    UBVECTOR fxstep(fxnew -fx); // change in F( x ).  We will need this for the secant update

    // log the worst market info
    const SolutionInfo* maxred = cSolInfo->getWorstSolutionInfo();
    addIteration(maxred->getName(), maxred->getRelativeED());
    if( mLogPricep ) {
      worstMarketLog << "Broyden-logPrice:  " << *maxred << "\n";
    }
    else {
      worstMarketLog << "Broyden-linearPrice:  " << *maxred << "\n";
    }

    // test for convergence
    double maxval = fabs(fxnew[0]);
    double imaxval = 0;
    for(size_t i=1; i<fxnew.size(); ++i) {
      double val = fabs(fxnew[i]);
      if(val > maxval) {
        maxval = val;
        imaxval = i;
      }
    }

    solverLog << "Convergence test maxval: " << maxval << "  imaxval= " << imaxval << "\n";
    solverLog << "\tx[i]= " << xnew[imaxval] << "  dx[i]= " << dx[imaxval] << "  xstep[i]= "
              << xstep[imaxval] << "\n";
    if(maxval <= mFTOL) {
      solverLog << "Solution successful.\n";
      x = xnew;
      fx = fxnew;
      return 0;                 // SUCCESS
    }

#if 0
    // secondary convergence test based on an estimated change in the
    // price vector.  Only test this if we have a "fresh" jacobian
    if(ageB == 0)  {
        maxval = fabs(fxnew[0]) / (util::getSmallNumber() + fabs(Btmp(0,0)));
        imaxval = 0;
        for(size_t i=1; i<fxnew.size(); ++i) {
            double val = fabs(fxnew[i]) / (util::getSmallNumber() + fabs(Btmp(i,i)));
            if(val > maxval) {
                maxval = val;
                imaxval = i;
            }
        }
        solverLog << "Secondary convergence test maxval:  " << maxval << "  imaxval= " << imaxval << "\n";
        if(maxval <= mFTOL) {   // XXX should put an x-tolerance in here
            solverLog << "Solved using secondary criterion.\n";
            x = xnew;
            fx = fxnew;
            return 0;
        }
    }
#endif

    // update B for next iteration
    double fratio_cutoff = 1.0 - 1.0/nrow;
    if(fnew/f0 < fratio_cutoff) { // making adequate progress with the Broyden formula
      double dx2 = inner_prod(xstep,xstep);
      UBVECTOR Bdx(F.nrtn());
      B = Btmp;
      fxstep -= axpy_prod(B, xstep, Bdx);
      fxstep /= dx2;
      B += outer_prod(fxstep, xstep);
      ageB++;                // increment the age of B
    }
    else {
      // Progress using the Broyden formula is anemic.  This usually
      // happens near discontinuities in the Jacobian matrix.  If B is
      // old, try a finite-difference jacobian to get us back on track.
      if(ageB > 0) {
        solverLog << "Insufficient progress with Broyden formula.  Resetting the Jacobian.\n(f0= " << f0 << ", fnew= " << fnew << ")\n";
        // just in case call fdjac such that it re-calculates the model at xnew
        // otherwise we could have bad state data from which we calculate derivatives
        fdjac(F,xnew,B);
        neval += x.size();
        ageB = 0;

        // Log the results of the Jacobian reset
        for(int j=0; j<F.narg(); ++j) {
            jdiag[j] = B(j,j);
        }

        solverLog << "New Jacobian:  diag( B )=\n" << jdiag << "\n";
        static_cast<LogEDFun&>(F).setSlope(jdiag);

      }
      else {
        // just did a reset, and it didn't help us.  Probably we've
        // got a very ill-behaved value in one of the variables.  Kick
        // it out and see if the bracketing routine can fix it.
        solverLog << "Repeated poor progress in Broyden solver.  Returning.\n";
        return -4;
      }
    }

    // log the data trace before we do the update
    reportVec("x", xnew, mktids_solv, issolvable_solv);
    reportVec("fx", fxnew, mktids_solv, issolvable_solv);
    reportVec("deltax", xnew-x, mktids_solv, issolvable_solv);    // xstep may have been modified above
    reportVec("deltafx", fxnew-fx, mktids_solv, issolvable_solv); // fxstep definitely modified above
    reportVec("diagB", jdiag, mktids_solv, issolvable_solv);
    reportPSD(rptvec_all, mktids_all, issolvable_all);                // report price, supply, and demand.
    mPerIter++;

    // update x, fx, f0 for next iteration
    f0 = fnew;
    x  = xnew;
    fx = fxnew;

  }

  // if we get here, then we didn't converge in the number of
  // iterations allowed us.  Return an error code
  solverLog << "\n****************Maximum solver iterations exceeded.\nlastx: " << x
            << "\nlastF: " << fx << "\n";
  return -1;
}

/*! \brief Write a vector into the solver data log
 *
 *  \details We write the solver data log in "long" format; i.e., with
 *           one completely-specified data item per line.  This allows
 *           us to simply skip variables that we don't have data for
 *           (vs. filling in NaN values), and we don't have to commit
 *           to a fixed number of variables; e.g., we could have more
 *           variables for price, supply, and demand.
 *
 *           The output columns are:
 *           period, iteration, variable name, market id, solvable (T/F), value
 *
 */
void LogBroyden::reportVec(const std::string &aname, const UBVECTOR &av, const std::vector<int> &amktids,
                           const std::vector<bool> &aissolvable)
{
    ILogger &datalog = ILogger::getLogger("solver-data-log");
    datalog.setLevel(ILogger::DEBUG);
    // Skip preparing the output if it won't be printed
    if(!datalog.wouldPrint(ILogger::DEBUG)) {
        return;
    }

    for(unsigned int i=0; i<av.size(); ++i) {
        datalog << mLastPer <<  " , " << mPerIter
                << ",\"" << aname << "\""
                << ", " << amktids[i]
                << (aissolvable[i] ? ", T" : ", F")
                << ", " << av[i] << "\n";
    }
}

void LogBroyden::reportPSD(UBVECTOR &arptvec, const std::vector<int> &amktids, const std::vector<bool> &aissolvable)
{
    // Skip preparing the log output if we're not going to print it anyhow
    if(!ILogger::getLogger("solver-data-log").wouldPrint(ILogger::DEBUG)) {
        return;
    }

    if(!cSolInfo) {
        ILogger &solverlog = ILogger::getLogger("solver_log");
        ILogger::WarningLevel olvl = solverlog.setLevel(ILogger::ERROR);
        solverlog << "cSolInfo is not set.  Data log will not include price, supply, or demand." << std::endl;
        solverlog.setLevel( olvl );
        return;
    }

    std::vector<SolutionInfo> solvable = cSolInfo->getSolvableSet();
    std::vector<SolutionInfo> unsolvable = cSolInfo->getUnsolvableSet();

    unsigned int i;
    unsigned int j;
    // log prices
    for(i=0; i<solvable.size(); ++i) {
        arptvec[i] = solvable[i].getPrice();
    }
    for(i=0,j=solvable.size(); i<unsolvable.size(); ++i,++j) {
        arptvec[j] = unsolvable[i].getPrice();
    }

    reportVec("price", arptvec, amktids, aissolvable);

    // log supply
    for(i=0; i<solvable.size(); ++i) {
        arptvec[i] = solvable[i].getSupply();
    }
    for(i=0,j=solvable.size(); i<unsolvable.size(); ++i,++j) {
        arptvec[j] = unsolvable[i].getSupply();
    }
    reportVec("supply", arptvec, amktids, aissolvable);

    // log demand
    for(i=0; i<solvable.size(); ++i) {
        arptvec[i] = solvable[i].getDemand();
    }
    for(i=0,j=solvable.size(); i<unsolvable.size(); ++i,++j) {
        arptvec[j] = unsolvable[i].getDemand();
    }
    reportVec("demand", arptvec, amktids, aissolvable);

}


text:
#pragma once

#include <boost/filesystem.hpp>
#include <rai/common/errors.hpp>
#include <rai/common/util.hpp>
#include <rai/node/node.hpp>
#include <rai/node/rpc.hpp>

namespace rai
{
class Daemon
{
public:
    rai::ErrorCode Run(const boost::filesystem::path&,
                       const boost::filesystem::path&);
};

class DaemonConfig
{
public:
    DaemonConfig(const boost::filesystem::path&);
    rai::ErrorCode DeserializeJson(bool&, rai::Ptree&);
    void SerializeJson(rai::Ptree&) const;
    rai::ErrorCode UpgradeJson(bool&, uint32_t, rai::Ptree&) const;

    rai::NodeConfig node_;
    rai::NodeRpcConfig rpc_;
};
}

ext: hpp


text: #include <boost/typeof/vector200.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Copyright 2012 The MITRE Corporation                                      *
 *                                                                           *
 * Licensed under the Apache License, Version 2.0 (the "License");           *
 * you may not use this file except in compliance with the License.          *
 * You may obtain a copy of the License at                                   *
 *                                                                           *
 *     http://www.apache.org/licenses/LICENSE-2.0                            *
 *                                                                           *
 * Unless required by applicable law or agreed to in writing, software       *
 * distributed under the License is distributed on an "AS IS" BASIS,         *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  *
 * See the License for the specific language governing permissions and       *
 * limitations under the License.                                            *
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#include <Eigen/Dense>

#include <openbr/plugins/openbr_internal.h>

using namespace cv;

namespace br
{

/*!
 * \ingroup transforms
 * \brief Sliding window feature extraction from a multi-channel integral image.
 * \author Josh Klontz \cite jklontz
 */
class IntegralSamplerTransform : public UntrainableTransform
{
    Q_OBJECT
    Q_PROPERTY(int scales READ get_scales WRITE set_scales RESET reset_scales STORED false)
    Q_PROPERTY(float scaleFactor READ get_scaleFactor WRITE set_scaleFactor RESET reset_scaleFactor STORED false)
    Q_PROPERTY(float stepFactor READ get_stepFactor WRITE set_stepFactor RESET reset_stepFactor STORED false)
    Q_PROPERTY(int minSize READ get_minSize WRITE set_minSize RESET reset_minSize STORED false)
    Q_PROPERTY(bool secondOrder READ get_secondOrder WRITE set_secondOrder RESET reset_secondOrder STORED false)
    BR_PROPERTY(int, scales, 6)
    BR_PROPERTY(float, scaleFactor, 2)
    BR_PROPERTY(float, stepFactor, 0.75)
    BR_PROPERTY(int, minSize, 8)
    BR_PROPERTY(bool, secondOrder, false)

    void project(const Template &src, Template &dst) const
    {
        typedef Eigen::Map< const Eigen::Matrix<qint32,Eigen::Dynamic,1> > InputDescriptor;
        typedef Eigen::Map< const Eigen::Matrix<float,Eigen::Dynamic,1> > SecondOrderInputDescriptor;
        typedef Eigen::Map< Eigen::Matrix<float,Eigen::Dynamic,1> > OutputDescriptor;

        const Mat &m = src.m();
        if (m.depth() != CV_32S) qFatal("Expected CV_32S matrix depth.");
        const int channels = m.channels();
        const int rowStep = channels * m.cols;

        int descriptors = 0;
        float idealSize = min(m.rows, m.cols)-1;
        for (int scale=0; scale<scales; scale++) {
            const int currentSize(idealSize);
            const int numDown = 1+(m.rows-currentSize-1)/int(idealSize*stepFactor);
            const int numAcross = 1+(m.cols-currentSize-1)/int(idealSize*stepFactor);
            descriptors += numDown*numAcross;
            if (secondOrder) descriptors += numDown*(numAcross-1) + (numDown-1)*numAcross;
            idealSize /= scaleFactor;
            if (idealSize < minSize) break;
        }
        Mat n(descriptors, channels, CV_32FC1);

        const qint32 *dataIn = (qint32*)m.data;
        float *dataOut = (float*)n.data;
        idealSize = min(m.rows, m.cols)-1;
        int index = 0;
        for (int scale=0; scale<scales; scale++) {
            const int currentSize(idealSize);
            const int currentStep(idealSize*stepFactor);
            for (int i=currentSize; i<m.rows; i+=currentStep) {
                for (int j=currentSize; j<m.cols; j+=currentStep) {
                    InputDescriptor a(dataIn+((i-currentSize)*rowStep+(j-currentSize)*channels), channels, 1);
                    InputDescriptor b(dataIn+((i-currentSize)*rowStep+ j             *channels), channels, 1);
                    InputDescriptor c(dataIn+(i              *rowStep+(j-currentSize)*channels), channels, 1);
                    InputDescriptor d(dataIn+(i              *rowStep+ j             *channels), channels, 1);
                    OutputDescriptor y(dataOut+(index*channels), channels, 1);
                    y = (d-b-c+a).cast<float>()/(currentSize*currentSize);
                    index++;
                }
            }
            if (secondOrder) {
                const int numDown = 1+(m.rows-currentSize-1)/currentStep;
                const int numAcross = 1+(m.cols-currentSize-1)/currentStep;
                const float *dataIn = n.ptr<float>(index - numDown*numAcross);
                for (int i=0; i<numDown; i++) {
                    for (int j=0; j<numAcross; j++) {
                        SecondOrderInputDescriptor a(dataIn + (i*numAcross+j)*channels, channels, 1);
                        if (j < numAcross-1) {
                            OutputDescriptor y(dataOut+(index*channels), channels, 1);
                            y = a - SecondOrderInputDescriptor(dataIn + (i*numAcross+j+1)*channels, channels, 1);
                            index++;
                        }
                        if (i < numDown-1) {
                            OutputDescriptor y(dataOut+(index*channels), channels, 1);
                            y = a - SecondOrderInputDescriptor(dataIn + ((i+1)*numAcross+j)*channels, channels, 1);
                            index++;
                        }
                    }
                }
            }
            idealSize /= scaleFactor;
            if (idealSize < minSize) break;
        }

        if (descriptors != index)
            qFatal("Allocated %d descriptors but computed %d.", descriptors, index);

        dst.m() = n;
    }
};

BR_REGISTER(Transform, IntegralSamplerTransform)

} // namespace br

#include "imgproc/integralsampler.moc"

ext: cpp

max_stars_repo_name: kassemitani/openbr

max_issues_repo_name: kassemitani/openbr

max_forks_repo_name: kassemitani/openbr


text:
//==================================================================================================
/**
  Copyright 2017 NumScale SAS

  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)
**/
//==================================================================================================

#ifndef BOOST_SIMD_FUNCTION_SCALAR_COMPARE_EQUAL_HPP_INCLUDED
#define BOOST_SIMD_FUNCTION_SCALAR_COMPARE_EQUAL_HPP_INCLUDED

#include <boost/simd/function/definition/compare_equal.hpp>
#include <boost/simd/arch/common/scalar/function/compare_equal.hpp>

#endif

max_stars_repo_name: SylvainCorlay/pythran

max_issues_repo_name: SylvainCorlay/pythran

max_forks_repo_name: SylvainCorlay/pythran


text:
//
// Created by Marco Bassaletti on 07-03-21.
//
#include <sstream>
#include <fstream>
#include <boost/json.hpp>
#include "JsonWriter.h"

namespace NSPass::Storage {
	using namespace std;

	JsonWriter::JsonWriter(boost::json::value& root) : root{ root }
	{
	}

	std::string JsonWriter::save()
	{
		ostringstream os;
		os << root;
		return os.str();
	}
}

max_stars_repo_name: mbassale/ownpass

max_issues_repo_name: mbassale/ownpass

max_forks_repo_name: mbassale/ownpass


text:
/*
 *
 * Copyright (c) 2002
 * John Maddock
 *
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0. (See accompanying file
 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 *
 */

#include "./regex_comparison.hpp"
#include <boost/timer.hpp>
#include <boost/regex.hpp>

namespace b{

double time_match(const std::string& re, const std::string& text)
{
    boost::regex e(re, boost::regex_constants::ECMAScript | boost::regex_constants::optimize);
    boost::smatch what;
    boost::timer tim;
    int iter = 1;
    int counter, repeats;
    double result = 0;
    double run;
    do
    {
        tim.restart();
        for(counter = 0; counter < iter; ++counter)
        {
            boost::regex_match(text, what, e);
        }
        result = tim.elapsed();
        iter *= 2;
    }while(result < 0.5);
    iter /= 2;

    // repeat test and report least value for consistency:
    for(repeats = 0; repeats < REPEAT_COUNT; ++repeats)
    {
        tim.restart();
        for(counter = 0; counter < iter; ++counter)
        {
            boost::regex_match(text, what, e);
        }
        run = tim.elapsed();
        result = (std::min)(run, result);
    }
    return result / iter;
}

//bool dummy_grep_proc(const boost::smatch&)
//{ return true; }

struct noop
{
    void operator()( boost::smatch const & ) const
    {
    }
};

double time_find_all(const std::string& re, const std::string& text)
{
    boost::regex e(re, boost::regex_constants::ECMAScript | boost::regex_constants::optimize);
    boost::smatch what;
    boost::timer tim;
    int iter = 1;
    int counter, repeats;
    double result = 0;
    double run;
    do
    {
        tim.restart();
        for(counter = 0; counter < iter; ++counter)
        {
            boost::sregex_iterator begin( text.begin(), text.end(), e ), end;
            std::for_each( begin, end, noop() );
            //boost::regex_grep(&dummy_grep_proc, text, e);
        }
        result = tim.elapsed();
        iter *= 2;
    }while(result < 0.5);
    iter /= 2;

    if(result >10)
        return result / iter;

    // repeat test and report least value for consistency:
    for(repeats = 0; repeats < REPEAT_COUNT; ++repeats)
    {
        tim.restart();
        for(counter = 0; counter < iter; ++counter)
        {
            boost::sregex_iterator begin( text.begin(), text.end(), e ), end;
            std::for_each( begin, end, noop() );
            //boost::regex_grep(&dummy_grep_proc, text, e);
        }
        run = tim.elapsed();
        result = (std::min)(run, result);
    }
    return result / iter;
}

}


text:
//
// Copyright © 2017 Arm Ltd. All rights reserved.
// SPDX-License-Identifier: MIT
//
#include "armnn/Tensor.hpp"
#include "armnn/Utils.hpp"
#include "armnn/Exceptions.hpp"
#include "armnn/TypesUtils.hpp"

#include <boost/assert.hpp>
#include <boost/log/trivial.hpp>
#include <boost/numeric/conversion/cast.hpp>

#include <sstream>

namespace armnn
{

// ---
// --- TensorShape
// ---

TensorShape::TensorShape()
 : m_NumDimensions(0)
{
}

TensorShape::TensorShape(unsigned int numDimensions)
 : m_NumDimensions(numDimensions)
{
    if (numDimensions < 1)
    {
        throw InvalidArgumentException("Tensor numDimensions must be greater than 0");
    }

    if (numDimensions > MaxNumOfTensorDimensions)
    {
        throw InvalidArgumentException("Tensor numDimensions must be less than or equal to MaxNumOfTensorDimensions");
    }

    std::fill(m_Dimensions.begin(), m_Dimensions.begin() + m_NumDimensions, 0);
}

TensorShape::TensorShape(const unsigned int numDimensions, const unsigned int* const dimensionSizes)
 : m_NumDimensions(numDimensions)
{
    if (numDimensions < 1)
    {
        throw InvalidArgumentException("Tensor numDimensions must be greater than 0");
    }

    if (numDimensions > MaxNumOfTensorDimensions)
    {
        throw InvalidArgumentException("Tensor numDimensions must be less than or equal to MaxNumOfTensorDimensions");
    }

    if (dimensionSizes == nullptr)
    {
        throw InvalidArgumentException("Tensor dimensionSizes must not be NULL");
    }

    std::copy(dimensionSizes, dimensionSizes + numDimensions, m_Dimensions.begin());
}

TensorShape::TensorShape(std::initializer_list<unsigned int> dimensionSizeList)
 : TensorShape(boost::numeric_cast<unsigned int>(dimensionSizeList.size()), dimensionSizeList.begin())
{
}

TensorShape::TensorShape(const TensorShape& other)
 : m_NumDimensions(other.m_NumDimensions)
{
    std::copy(other.m_Dimensions.cbegin(), other.m_Dimensions.cbegin() + other.m_NumDimensions, m_Dimensions.begin());
}

TensorShape& TensorShape::operator =(const TensorShape& other)
{
    m_NumDimensions = other.m_NumDimensions;
    std::copy(other.m_Dimensions.cbegin(), other.m_Dimensions.cbegin() + other.m_NumDimensions, m_Dimensions.begin());
    return *this;
}

unsigned int TensorShape::operator[](unsigned int i) const
{
    CheckDimensionIndex(i);
    return m_Dimensions.at(i);
}

unsigned int& TensorShape::operator[](unsigned int i)
{
    CheckDimensionIndex(i);
    return m_Dimensions.at(i);
}

bool TensorShape::operator==(const TensorShape& other) const
{
    return ((m_NumDimensions == other.m_NumDimensions) &&
        std::equal(m_Dimensions.cbegin(), m_Dimensions.cbegin() + m_NumDimensions, other.m_Dimensions.cbegin()));
}

bool TensorShape::operator!=(const TensorShape& other) const
{
    return !(*this == other);
}

unsigned int TensorShape::GetNumElements() const
{
    if (m_NumDimensions == 0)
    {
        return 0;
    }

    unsigned int count = 1;
    for (unsigned int i = 0; i < m_NumDimensions; i++)
    {
        count *= m_Dimensions[i];
    }

    return count;
}

void TensorShape::CheckDimensionIndex(unsigned int i) const
{
    if (i >= m_NumDimensions)
    {
        std::stringstream errorMessage;
        errorMessage << "Invalid dimension index: " << i << " (number of dimensions is " << m_NumDimensions << ")";
        throw InvalidArgumentException(errorMessage.str(), CHECK_LOCATION());
    }
}

// ---
// --- TensorInfo
// ---

TensorInfo::TensorInfo()
: m_DataType(DataType::Float32)
{
}

TensorInfo::TensorInfo(const TensorShape& shape, DataType dataType,
    float quantizationScale, int32_t quantizationOffset)
 : m_Shape(shape)
 , m_DataType(dataType)
{
    m_Quantization.m_Scale = quantizationScale;
    m_Quantization.m_Offset = quantizationOffset;
}

TensorInfo::TensorInfo(unsigned int numDimensions, const unsigned int* dimensionSizes, DataType dataType,
    float quantizationScale, int32_t quantizationOffset)
 : m_Shape(numDimensions, dimensionSizes)
 , m_DataType(dataType)
{
    m_Quantization.m_Scale = quantizationScale;
    m_Quantization.m_Offset = quantizationOffset;
}

TensorInfo::TensorInfo(const TensorInfo& other)
: m_Shape(other.m_Shape)
, m_DataType(other.m_DataType)
, m_Quantization(other.m_Quantization)
{
}

TensorInfo& TensorInfo::operator=(const TensorInfo& other)
{
    m_Shape = other.m_Shape;
    m_DataType = other.m_DataType;
    m_Quantization = other.m_Quantization;
    return *this;
}

bool TensorInfo::operator==(const TensorInfo& other) const
{
    return ((m_Shape == other.m_Shape) &&
            (m_DataType == other.m_DataType) &&
            (m_Quantization == other.m_Quantization));
}

bool TensorInfo::operator!=(const TensorInfo& other) const
{
    return !(*this == other);
}

unsigned int TensorInfo::GetNumBytes() const
{
    return GetDataTypeSize(m_DataType) * GetNumElements();
}

bool TensorInfo::IsTypeSpaceMatch(const TensorInfo& other) const
{
    bool match = true;

    match &= m_DataType == other.m_DataType;

    if (IsQuantized())
    {
        match &= GetQuantizationScale() == other.GetQuantizationScale() &&
                 GetQuantizationOffset() == other.GetQuantizationOffset();
    }
    return match;
}

// ---
// --- BaseTensor
// ---

template<typename MemoryType>
BaseTensor<MemoryType>::BaseTensor()
 : m_MemoryArea(nullptr)
{
}

template<typename MemoryType>
BaseTensor<MemoryType>::BaseTensor(const TensorInfo& info, MemoryType memoryArea)
 : m_MemoryArea(memoryArea)
 , m_Info(info)
{
}

template<typename MemoryType>
BaseTensor<MemoryType>::BaseTensor(const BaseTensor<MemoryType>& other)
 : m_MemoryArea(other.m_MemoryArea)
 , m_Info(other.GetInfo())
{
}

template<typename MemoryType>
BaseTensor<MemoryType>& BaseTensor<MemoryType>::operator =(const BaseTensor<MemoryType>& other)
{
    m_Info = other.m_Info;
    m_MemoryArea = other.m_MemoryArea;
    return *this;
}

// Explicit instantiations.
template class BaseTensor<const void*>;
template class BaseTensor<void*>;

} // namespace armnn

max_stars_repo_name: VinayKarnam/armnn

max_issues_repo_name: VinayKarnam/armnn

max_forks_repo_name: VinayKarnam/armnn


text:
/* Author: Wolfgang Bangerth, Texas A&M University, 2006, 2007 */

/*    $Id: step-27.cc 27657 2012-11-21 13:19:08Z bangerth $       */
/*                                                                */
/*    Copyright (C) 2006-2012 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

// @sect3{Include files}

// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>

// These are the new files we need. The first one provides an alternative to
// the usual SparsityPattern class and the CompressedSparsityPattern class
// already discussed in step-11 and step-18. The last two provide hp
// versions of the DoFHandler and FEValues classes as described in the
// introduction of this program.
#include <deal.II/lac/compressed_set_sparsity_pattern.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_values.h>

// The last set of include files are standard C++ headers. We need support for
// complex numbers when we compute the Fourier transform.
#include <fstream>
#include <iostream>
#include <complex>

// Finally, this is as in previous programs:
namespace Step27
{
  using namespace dealii;

  // @sect3{The main class}

  // The main class of this program looks very much like the one already used
  // in the first few tutorial programs, for example the one in step-6. The
  // main difference is that we have merged the refine_grid and output_results
  // functions into one since we will also want to output some of the
  // quantities used in deciding how to refine the mesh (in particular the
  // estimated smoothness of the solution). There is also a function that
  // computes this estimated smoothness, as discussed in the introduction.
  //
  // As far as member variables are concerned, we use the same structure as
  // already used in step-6, but instead of a regular DoFHandler we use an
  // object of type hp::DoFHandler, and we need collections instead of
  // individual finite element, quadrature, and face quadrature objects. We
  // will fill these collections in the constructor of the class. The last
  // variable, max_degree, indicates the maximal polynomial
  // degree of shape functions used.
  template <int dim>
  class LaplaceProblem
  {
  public:
    LaplaceProblem ();
    ~LaplaceProblem ();

    void run ();

  private:
    void setup_system ();
    void assemble_system ();
    void solve ();
    void create_coarse_grid ();
    void estimate_smoothness (Vector<float> &smoothness_indicators) const;
    void postprocess (const unsigned int cycle);

    Triangulation<dim>   triangulation;

    hp::DoFHandler<dim>      dof_handler;
    hp::FECollection<dim>    fe_collection;
    hp::QCollection<dim>     quadrature_collection;
    hp::QCollection<dim-1>   face_quadrature_collection;

    ConstraintMatrix     constraints;

    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;

    Vector<double>       solution;
    Vector<double>       system_rhs;

    const unsigned int max_degree;
  };

  // @sect3{Equation data}
  //
  // Next, let us define the right hand side function for this problem. It is
  // $x+1$ in 1d, $(x+1)(y+1)$ in 2d, and so on.
  template <int dim>
  class RightHandSide : public Function<dim>
  {
  public:
    RightHandSide () : Function<dim> () {}

    virtual double value (const Point<dim>   &p,
                          const unsigned int  component) const;
  };

  template <int dim>
  double
  RightHandSide<dim>::value (const Point<dim>   &p,
                             const unsigned int  /*component*/) const
  {
    double product = 1;
    for (unsigned int d=0; d<dim; ++d)
      product *= (p[d]+1);
    return product;
  }

  // @sect3{Implementation of the main class}

  // @sect4{LaplaceProblem::LaplaceProblem}

  // The constructor of this class is fairly straightforward. It associates
  // the hp::DoFHandler object with the triangulation, and then sets the
  // maximal polynomial degree to 7 (in 1d and 2d) or 5 (in 3d and higher). We
  // do so because using higher order polynomial degrees becomes prohibitively
  // expensive, especially in higher space dimensions.
  //
  // Following this, we fill the collections of finite element, and cell and
  // face quadrature objects. We start with quadratic elements, and each
  // quadrature formula is chosen so that it is appropriate for the matching
  // finite element in the hp::FECollection object.
  template <int dim>
  LaplaceProblem<dim>::LaplaceProblem ()
    :
    dof_handler (triangulation),
    max_degree (dim <= 2 ? 7 : 5)
  {
    for (unsigned int degree=2; degree<=max_degree; ++degree)
      {
        fe_collection.push_back (FE_Q<dim>(degree));
        quadrature_collection.push_back (QGauss<dim>(degree+1));
        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
      }
  }

  // @sect4{LaplaceProblem::~LaplaceProblem}

  // The destructor is unchanged from what we already did in step-6:
  template <int dim>
  LaplaceProblem<dim>::~LaplaceProblem ()
  {
    dof_handler.clear ();
  }

  // @sect4{LaplaceProblem::setup_system}
  //
  // This function is again an almost verbatim copy of what we already did in
  // step-6. The first change is that we append the Dirichlet boundary
  // conditions to the ConstraintMatrix object, which we consequently call
  // just constraints instead of
  // hanging_node_constraints. The second difference is that we
  // don't directly build the sparsity pattern, but first create an
  // intermediate object that we later copy into the usual SparsityPattern
  // data structure, since this is more efficient for the problem with many
  // entries per row (and different number of entries in different rows). In
  // another slight deviation, we do not first build the sparsity pattern and
  // then condense away constrained degrees of freedom, but pass the
  // constraint matrix object directly to the function that builds the
  // sparsity pattern. We disable the insertion of constrained entries with
  // <tt>false</tt> as fourth argument in the DoFTools::make_sparsity_pattern
  // function. All of these changes are explained in the introduction of this
  // program.
  //
  // The last change, maybe hidden in plain sight, is that the dof_handler
  // variable here is an hp object -- nevertheless all the function calls we
  // had before still work in exactly the same way as they always did.
  template <int dim>
  void LaplaceProblem<dim>::setup_system ()
  {
    dof_handler.distribute_dofs (fe_collection);

    solution.reinit (dof_handler.n_dofs());
    system_rhs.reinit (dof_handler.n_dofs());

    constraints.clear ();
    DoFTools::make_hanging_node_constraints (dof_handler,
                                             constraints);
    VectorTools::interpolate_boundary_values (dof_handler,
                                              0,
                                              ZeroFunction<dim>(),
                                              constraints);
    constraints.close ();

    CompressedSetSparsityPattern csp (dof_handler.n_dofs(),
                                      dof_handler.n_dofs());
    DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
    sparsity_pattern.copy_from (csp);

    system_matrix.reinit (sparsity_pattern);
  }

  // @sect4{LaplaceProblem::assemble_system}

  // This is the function that assembles the global matrix and right hand side
  // vector from the local contributions of each cell. Its main working is as
  // has been described in many of the tutorial programs before. The
  // significant deviations are the ones necessary for hp finite
  // element methods. In particular, that we need to use a collection of
  // FEValues object (implemented through the hp::FEValues class), and that we
  // have to eliminate constrained degrees of freedom already when copying
  // local contributions into global objects. Both of these are explained in
  // detail in the introduction of this program.
  //
  // One other slight complication is the fact that because we use different
  // polynomial degrees on different cells, the matrices and vectors holding
  // local contributions do not have the same size on all cells. At the
  // beginning of the loop over all cells, we therefore each time have to
  // resize them to the correct size (given by
  // dofs_per_cell). Because these classes are implement in such
  // a way that reducing the size of a matrix or vector does not release the
  // currently allocated memory (unless the new size is zero), the process of
  // resizing at the beginning of the loop will only require re-allocation of
  // memory during the first few iterations. Once we have found in a cell with
  // the maximal finite element degree, no more re-allocations will happen
  // because all subsequent reinit calls will only set the size
  // to something that fits the currently allocated memory. This is important
  // since allocating memory is expensive, and doing so every time we visit a
  // new cell would take significant compute time.
  template <int dim>
  void LaplaceProblem<dim>::assemble_system ()
  {
    hp::FEValues<dim> hp_fe_values (fe_collection,
                                    quadrature_collection,
                                    update_values    |  update_gradients |
                                    update_quadrature_points  |  update_JxW_values);

    const RightHandSide<dim> rhs_function;

    FullMatrix<double>   cell_matrix;
    Vector<double>       cell_rhs;

    std::vector<unsigned int> local_dof_indices;

    typename hp::DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
    for (; cell!=endc; ++cell)
      {
        const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;

        cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
        cell_matrix = 0;

        cell_rhs.reinit (dofs_per_cell);
        cell_rhs = 0;

        hp_fe_values.reinit (cell);

        const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();

        std::vector<double>  rhs_values (fe_values.n_quadrature_points);
        rhs_function.value_list (fe_values.get_quadrature_points(),
                                 rhs_values);

        for (unsigned int q_point=0;
             q_point<fe_values.n_quadrature_points;
             ++q_point)
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
              for (unsigned int j=0; j<dofs_per_cell; ++j)
                cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
                                     fe_values.shape_grad(j,q_point) *
                                     fe_values.JxW(q_point));

              cell_rhs(i) += (fe_values.shape_value(i,q_point) *
                              rhs_values[q_point] *
                              fe_values.JxW(q_point));
            }

        local_dof_indices.resize (dofs_per_cell);
        cell->get_dof_indices (local_dof_indices);

        constraints.distribute_local_to_global (cell_matrix, cell_rhs,
                                                local_dof_indices,
                                                system_matrix, system_rhs);
      }

    // Now with the loop over all cells finished, we are done for this
    // function. The steps we still had to do at this point in earlier
    // tutorial programs, namely condensing hanging node constraints and
    // applying Dirichlet boundary conditions, have been taken care of by the
    // ConstraintMatrix object constraints on the fly.
  }

  // @sect4{LaplaceProblem::solve}

  // The function solving the linear system is entirely unchanged from
  // previous examples. We simply try to reduce the initial residual (which
  // equals the $l_2$ norm of the right hand side) by a certain factor:
  template <int dim>
  void LaplaceProblem<dim>::solve ()
  {
    SolverControl           solver_control (system_rhs.size(),
                                            1e-8*system_rhs.l2_norm());
    SolverCG<>              cg (solver_control);

    PreconditionSSOR<> preconditioner;
    preconditioner.initialize(system_matrix, 1.2);

    cg.solve (system_matrix, solution, system_rhs,
              preconditioner);

    constraints.distribute (solution);
  }

  // @sect4{LaplaceProblem::postprocess}

  // After solving the linear system, we will want to postprocess the
  // solution. Here, all we do is to estimate the error, estimate the local
  // smoothness of the solution as described in the introduction, then write
  // graphical output, and finally refine the mesh in both $h$ and $p$
  // according to the indicators computed before. We do all this in the same
  // function because we want the estimated error and smoothness indicators
  // not only for refinement, but also include them in the graphical output.
  template <int dim>
  void LaplaceProblem<dim>::postprocess (const unsigned int cycle)
  {
    // Let us start with computing estimated error and smoothness indicators,
    // which each are one number for each active cell of our
    // triangulation. For the error indicator, we use the KellyErrorEstimator
    // class as always. Estimating the smoothness is done in the respective
    // function of this class; that function is discussed further down below:
    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
    KellyErrorEstimator<dim>::estimate (dof_handler,
                                        face_quadrature_collection,
                                        typename FunctionMap<dim>::type(),
                                        solution,
                                        estimated_error_per_cell);

    Vector<float> smoothness_indicators (triangulation.n_active_cells());
    estimate_smoothness (smoothness_indicators);

    // Next we want to generate graphical output. In addition to the two
    // estimated quantities derived above, we would also like to output the
    // polynomial degree of the finite elements used on each of the elements
    // on the mesh.
    //
    // The way to do that requires that we loop over all cells and poll the
    // active finite element index of them using
    // cell-@>active_fe_index(). We then use the result of this
    // operation and query the finite element collection for the finite
    // element with that index, and finally determine the polynomial degree of
    // that element. The result we put into a vector with one element per
    // cell. The DataOut class requires this to be a vector of
    // float or double, even though our values are
    // all integers, so that it what we use:
    {
      Vector<float> fe_degrees (triangulation.n_active_cells());
      {
        typename hp::DoFHandler<dim>::active_cell_iterator
        cell = dof_handler.begin_active(),
        endc = dof_handler.end();
        for (unsigned int index=0; cell!=endc; ++cell, ++index)
          fe_degrees(index)
            = fe_collection[cell->active_fe_index()].degree;
      }

      // With now all data vectors available -- solution, estimated errors and
      // smoothness indicators, and finite element degrees --, we create a
      // DataOut object for graphical output and attach all data. Note that
      // the DataOut class has a second template argument (which defaults to
      // DoFHandler@<dim@>, which is why we have never seen it in previous
      // tutorial programs) that indicates the type of DoF handler to be
      // used. Here, we have to use the hp::DoFHandler class:
      DataOut<dim,hp::DoFHandler<dim> > data_out;

      data_out.attach_dof_handler (dof_handler);
      data_out.add_data_vector (solution, "solution");
      data_out.add_data_vector (estimated_error_per_cell, "error");
      data_out.add_data_vector (smoothness_indicators, "smoothness");
      data_out.add_data_vector (fe_degrees, "fe_degree");
      data_out.build_patches ();

      // The final step in generating output is to determine a file name, open
      // the file, and write the data into it (here, we use VTK format):
      const std::string filename = "solution-" +
                                   Utilities::int_to_string (cycle, 2) +
                                   ".vtk";
      std::ofstream output (filename.c_str());
      data_out.write_vtk (output);
    }

    // After this, we would like to actually refine the mesh, in both $h$ and
    // $p$. The way we are going to do this is as follows: first, we use the
    // estimated error to flag those cells for refinement that have the
    // largest error. This is what we have always done:
    {
      GridRefinement::refine_and_coarsen_fixed_number (triangulation,
                                                       estimated_error_per_cell,
                                                       0.3, 0.03);

      // Next we would like to figure out which of the cells that have been
      // flagged for refinement should actually have $p$ increased instead of
      // $h$ decreased. The strategy we choose here is that we look at the
      // smoothness indicators of those cells that are flagged for refinement,
      // and increase $p$ for those with a smoothness larger than a certain
      // threshold. For this, we first have to determine the maximal and
      // minimal values of the smoothness indicators of all flagged cells,
      // which we do using a loop over all cells and comparing current minimal
      // and maximal values. (We start with the minimal and maximal values of
      // all cells, a range within which the minimal and maximal values
      // on cells flagged for refinement must surely lie.) Absent any better
      // strategies, we will then set the threshold above which will increase
      // $p$ instead of reducing $h$ as the mean value between minimal and
      // maximal smoothness indicators on cells flagged for refinement:
      float max_smoothness = *std::min_element (smoothness_indicators.begin(),
                                                smoothness_indicators.end()),
                             min_smoothness = *std::max_element (smoothness_indicators.begin(),
                                                                 smoothness_indicators.end());
      {
        typename hp::DoFHandler<dim>::active_cell_iterator
        cell = dof_handler.begin_active(),
        endc = dof_handler.end();
        for (unsigned int index=0; cell!=endc; ++cell, ++index)
          if (cell->refine_flag_set())
            {
              max_smoothness = std::max (max_smoothness,
                                         smoothness_indicators(index));
              min_smoothness = std::min (min_smoothness,
                                         smoothness_indicators(index));
            }
      }
      const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;

      // With this, we can go back, loop over all cells again, and for those
      // cells for which (i) the refinement flag is set, (ii) the smoothness
      // indicator is larger than the threshold, and (iii) we still have a
      // finite element with a polynomial degree higher than the current one
      // in the finite element collection, we then increase the polynomial
      // degree and in return remove the flag indicating that the cell should
      // undergo bisection. For all other cells, the refinement flags remain
      // untouched:
      {
        typename hp::DoFHandler<dim>::active_cell_iterator
        cell = dof_handler.begin_active(),
        endc = dof_handler.end();
        for (unsigned int index=0; cell!=endc; ++cell, ++index)
          if (cell->refine_flag_set()
              &&
              (smoothness_indicators(index) > threshold_smoothness)
              &&
              (cell->active_fe_index()+1 < fe_collection.size()))
            {
              cell->clear_refine_flag();
              cell->set_active_fe_index (cell->active_fe_index() + 1);
            }
      }

      // At the end of this procedure, we then refine the mesh. During this
      // process, children of cells undergoing bisection inherit their mother
      // cell's finite element index:
      triangulation.execute_coarsening_and_refinement ();
    }
  }

  // @sect4{LaplaceProblem::create_coarse_grid}

  // The following function is used when creating the initial grid. It is a
  // specialization for the 2d case, i.e. a corresponding function needs to be
  // implemented if the program is run in anything other then 2d. The function
  // is actually stolen from step-14 and generates the same mesh used already
  // there, i.e. the square domain with the square hole in the middle. The
  // meaning of the different parts of this function are explained in the
  // documentation of step-14:
  template <>
  void LaplaceProblem<2>::create_coarse_grid ()
  {
    const unsigned int dim = 2;

    static const Point<2> vertices_1[]
      = {  Point<2> (-1.,   -1.),
           Point<2> (-1./2, -1.),
           Point<2> (0.,    -1.),
           Point<2> (+1./2, -1.),
           Point<2> (+1,    -1.),

           Point<2> (-1.,   -1./2.),
           Point<2> (-1./2, -1./2.),
           Point<2> (0.,    -1./2.),
           Point<2> (+1./2, -1./2.),
           Point<2> (+1,    -1./2.),

           Point<2> (-1.,   0.),
           Point<2> (-1./2, 0.),
           Point<2> (+1./2, 0.),
           Point<2> (+1,    0.),

           Point<2> (-1.,   1./2.),
           Point<2> (-1./2, 1./2.),
           Point<2> (0.,    1./2.),
           Point<2> (+1./2, 1./2.),
           Point<2> (+1,    1./2.),

           Point<2> (-1.,   1.),
           Point<2> (-1./2, 1.),
           Point<2> (0.,    1.),
           Point<2> (+1./2, 1.),
           Point<2> (+1,    1.)
        };
    const unsigned int
    n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
    const std::vector<Point<dim> > vertices (&vertices_1[0],
                                             &vertices_1[n_vertices]);
    static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
    = {{0, 1, 5, 6},
      {1, 2, 6, 7},
      {2, 3, 7, 8},
      {3, 4, 8, 9},
      {5, 6, 10, 11},
      {8, 9, 12, 13},
      {10, 11, 14, 15},
      {12, 13, 17, 18},
      {14, 15, 19, 20},
      {15, 16, 20, 21},
      {16, 17, 21, 22},
      {17, 18, 22, 23}
    };
    const unsigned int
    n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);

    std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
    for (unsigned int i=0; i<n_cells; ++i)
      {
        for (unsigned int j=0;
             j<GeometryInfo<dim>::vertices_per_cell;
             ++j)
          cells[i].vertices[j] = cell_vertices[i][j];
        cells[i].material_id = 0;
      }

    triangulation.create_triangulation (vertices,
                                        cells,
                                        SubCellData());
    triangulation.refine_global (3);
  }

  // @sect4{LaplaceProblem::run}

  // This function implements the logic of the program, as did the respective
  // function in most of the previous programs already, see for example
  // step-6.
  //
  // Basically, it contains the adaptive loop: in the first iteration create a
  // coarse grid, and then set up the linear system, assemble it, solve, and
  // postprocess the solution including mesh refinement. Then start over
  // again. In the meantime, also output some information for those staring at
  // the screen trying to figure out what the program does:
  template <int dim>
  void LaplaceProblem<dim>::run ()
  {
    for (unsigned int cycle=0; cycle<6; ++cycle)
      {
        std::cout << "Cycle " << cycle << ':' << std::endl;

        if (cycle == 0)
          create_coarse_grid ();

        setup_system ();

        std::cout << "   Number of active cells:       "
                  << triangulation.n_active_cells()
                  << std::endl
                  << "   Number of degrees of freedom: "
                  << dof_handler.n_dofs()
                  << std::endl
                  << "   Number of constraints       : "
                  << constraints.n_constraints()
                  << std::endl;

        assemble_system ();
        solve ();
        postprocess (cycle);
      }
  }

  // @sect4{LaplaceProblem::estimate_smoothness}

  // This last function of significance implements the algorithm to estimate
  // the smoothness exponent using the algorithms explained in detail in the
  // introduction. We will therefore only comment on those points that are of
  // implementational importance.
  template <int dim>
  void
  LaplaceProblem<dim>::
  estimate_smoothness (Vector<float> &smoothness_indicators) const
  {
    // The first thing we need to do is to define the Fourier vectors ${\bf
    // k}$ for which we want to compute Fourier coefficients of the solution
    // on each cell. In 2d, we pick those vectors ${\bf k}=(\pi i, \pi j)^T$
    // for which $\sqrt{i^2+j^2}\le N$, with $i,j$ integers and $N$ being the
    // maximal polynomial degree we use for the finite elements in this
    // program. The 3d case is handled analogously. 1d and dimensions higher
    // than 3 are not implemented, and we guard our implementation by making
    // sure that we receive an exception in case someone tries to compile the
    // program for any of these dimensions.
    //
    // We exclude ${\bf k}=0$ to avoid problems computing $|{\bf k}|^{-mu}$
    // and $\ln |{\bf k}|$. The other vectors are stored in the field
    // k_vectors. In addition, we store the square of the
    // magnitude of each of these vectors (up to a factor $\pi^2$) in the
    // k_vectors_magnitude array -- we will need that when we
    // attempt to find out which of those Fourier coefficients corresponding
    // to Fourier vectors of the same magnitude is the largest:
    const unsigned int N = max_degree;

    std::vector<Tensor<1,dim> > k_vectors;
    std::vector<unsigned int>   k_vectors_magnitude;
    switch (dim)
      {
      case 2:
      {
        for (unsigned int i=0; i<N; ++i)
          for (unsigned int j=0; j<N; ++j)
            if (!((i==0) && (j==0))
                &&
                (i*i + j*j < N*N))
              {
                k_vectors.push_back (Point<dim>(numbers::PI * i,
                                                numbers::PI * j));
                k_vectors_magnitude.push_back (i*i+j*j);
              }

        break;
      }

      case 3:
      {
        for (unsigned int i=0; i<N; ++i)
          for (unsigned int j=0; j<N; ++j)
            for (unsigned int k=0; k<N; ++k)
              if (!((i==0) && (j==0) && (k==0))
                  &&
                  (i*i + j*j + k*k < N*N))
                {
                  k_vectors.push_back (Point<dim>(numbers::PI * i,
                                                  numbers::PI * j,
                                                  numbers::PI * k));
                  k_vectors_magnitude.push_back (i*i+j*j+k*k);
                }

        break;
      }

      default:
        Assert (false, ExcNotImplemented());
      }

    // After we have set up the Fourier vectors, we also store their total
    // number for simplicity, and compute the logarithm of the magnitude of
    // each of these vectors since we will need it many times over further
    // down below:
    const unsigned n_fourier_modes = k_vectors.size();
    std::vector<double> ln_k (n_fourier_modes);
    for (unsigned int i=0; i<n_fourier_modes; ++i)
      ln_k[i] = std::log (k_vectors[i].norm());

    // Next, we need to assemble the matrices that do the Fourier transforms
    // for each of the finite elements we deal with, i.e. the matrices ${\cal
    // F}_{{\bf k},j}$ defined in the introduction. We have to do that for
    // each of the finite elements in use. Note that these matrices are
    // complex-valued, so we can't use the FullMatrix class. Instead, we use
    // the Table class template.
    std::vector<Table<2,std::complex<double> > >
    fourier_transform_matrices (fe_collection.size());

    // In order to compute them, we of course can't perform the Fourier
    // transform analytically, but have to approximate it using quadrature. To
    // this end, we use a quadrature formula that is obtained by iterating a
    // 2-point Gauss formula as many times as the maximal exponent we use for
    // the term $e^{i{\bf k}\cdot{\bf x}}$:
    QGauss<1>      base_quadrature (2);
    QIterated<dim> quadrature (base_quadrature, N);

    // With this, we then loop over all finite elements in use, reinitialize
    // the respective matrix ${\cal F}$ to the right size, and integrate each
    // entry of the matrix numerically as ${\cal F}_{{\bf k},j}=\sum_q
    // e^{i{\bf k}\cdot {\bf x}}\varphi_j({\bf x}_q) w_q$, where $x_q$ are the
    // quadrature points and $w_q$ are the quadrature weights. Note that the
    // imaginary unit $i=\sqrt{-1}$ is obtained from the standard C++ classes
    // using std::complex@<double@>(0,1).

    // Because we work on the unit cell, we can do all this work without a
    // mapping from reference to real cell and consequently do not need the
    // FEValues class.
    for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
      {
        fourier_transform_matrices[fe].reinit (n_fourier_modes,
                                               fe_collection[fe].dofs_per_cell);

        for (unsigned int k=0; k<n_fourier_modes; ++k)
          for (unsigned int j=0; j<fe_collection[fe].dofs_per_cell; ++j)
            {
              std::complex<double> sum = 0;
              for (unsigned int q=0; q<quadrature.size(); ++q)
                {
                  const Point<dim> x_q = quadrature.point(q);
                  sum += std::exp(std::complex<double>(0,1) *
                                  (k_vectors[k] * x_q)) *
                         fe_collection[fe].shape_value(j,x_q) *
                         quadrature.weight(q);
                }
              fourier_transform_matrices[fe](k,j)
                = sum / std::pow(2*numbers::PI, 1.*dim/2);
            }
      }

    // The next thing is to loop over all cells and do our work there, i.e. to
    // locally do the Fourier transform and estimate the decay coefficient. We
    // will use the following two arrays as scratch arrays in the loop and
    // allocate them here to avoid repeated memory allocations:
    std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
    Vector<double>                     local_dof_values;

    // Then here is the loop:
    typename hp::DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
    for (unsigned int index=0; cell!=endc; ++cell, ++index)
      {
        // Inside the loop, we first need to get the values of the local
        // degrees of freedom (which we put into the
        // local_dof_values array after setting it to the right
        // size) and then need to compute the Fourier transform by multiplying
        // this vector with the matrix ${\cal F}$ corresponding to this finite
        // element. We need to write out the multiplication by hand because
        // the objects holding the data do not have vmult-like
        // functions declared:
        local_dof_values.reinit (cell->get_fe().dofs_per_cell);
        cell->get_dof_values (solution, local_dof_values);

        for (unsigned int f=0; f<n_fourier_modes; ++f)
          {
            fourier_coefficients[f] = 0;

            for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
              fourier_coefficients[f] +=
                fourier_transform_matrices[cell->active_fe_index()](f,i)
                *
                local_dof_values(i);
          }

        // The next thing, as explained in the introduction, is that we wanted
        // to only fit our exponential decay of Fourier coefficients to the
        // largest coefficients for each possible value of $|{\bf k}|$. To
        // this end, we create a map that for each magnitude $|{\bf k}|$
        // stores the largest $|\hat U_{{\bf k}}|$ found so far, i.e. we
        // overwrite the existing value (or add it to the map) if no value for
        // the current $|{\bf k}|$ exists yet, or if the current value is
        // larger than the previously stored one:
        std::map<unsigned int, double> k_to_max_U_map;
        for (unsigned int f=0; f<n_fourier_modes; ++f)
          if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
               k_to_max_U_map.end())
              ||
              (k_to_max_U_map[k_vectors_magnitude[f]] <
               std::abs (fourier_coefficients[f])))
            k_to_max_U_map[k_vectors_magnitude[f]]
              = std::abs (fourier_coefficients[f]);
        // Note that it comes in handy here that we have stored the magnitudes
        // of vectors as integers, since this way we do not have to deal with
        // round-off-sized differences between different values of $|{\bf
        // k}|$.

        // As the final task, we have to calculate the various contributions
        // to the formula for $\mu$. We'll only take those Fourier
        // coefficients with the largest magnitude for a given value of $|{\bf
        // k}|$ as explained above:
        double  sum_1           = 0,
                sum_ln_k        = 0,
                sum_ln_k_square = 0,
                sum_ln_U        = 0,
                sum_ln_U_ln_k   = 0;
        for (unsigned int f=0; f<n_fourier_modes; ++f)
          if (k_to_max_U_map[k_vectors_magnitude[f]] ==
              std::abs (fourier_coefficients[f]))
            {
              sum_1 += 1;
              sum_ln_k += ln_k[f];
              sum_ln_k_square += ln_k[f]*ln_k[f];
              sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
              sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) *
                               ln_k[f];
            }

        // With these so-computed sums, we can now evaluate the formula for
        // $\mu$ derived in the introduction:
        const double mu
          = (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
             *
             (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));

        // The final step is to compute the Sobolev index $s=\mu-\frac d2$ and
        // store it in the vector of estimated values for each cell:
        smoothness_indicators(index) = mu - 1.*dim/2;
      }
  }
}

// @sect3{The main function}

// The main function is again verbatim what we had before: wrap creating and
// running an object of the main class into a try block and catch
// whatever exceptions are thrown, thereby producing meaningful output if
// anything should go wrong:
int main ()
{
  try
    {
      using namespace dealii;
      using namespace Step27;

      deallog.depth_console (0);

      LaplaceProblem<2> laplace_problem;
      laplace_problem.run ();
    }
  catch (std::exception &exc)
    {
      std::cerr << std::endl << std::endl
                << "----------------------------------------------------"
                << std::endl;
      std::cerr << "Exception on processing: " << std::endl
                << exc.what() << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;

      return 1;
    }
  catch (...)
    {
      std::cerr << std::endl << std::endl
                << "----------------------------------------------------"
                << std::endl;
      std::cerr << "Unknown exception!" << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;
      return 1;
    }

  return 0;
}

ext: cc

max_stars_repo_name: wathen/PhD

max_issues_repo_name: wathen/PhD

max_forks_repo_name: wathen/PhD


text:
/*ckwg +29
 * Copyright 2019 by Kitware, Inc.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 *  * Neither name of Kitware, Inc. nor the names of any contributors may be used
 *    to endorse or promote products derived from this software without specific
 *    prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <stdexcept>

#include <vital/vital_types.h>
#include <vital/types/homography.h>
#include <vital/types/homography_f2f.h>

#include <Eigen/Core>

#include <pybind11/pybind11.h>

namespace kv = kwiver::vital;
namespace py = pybind11;

using f2f_homography = kv::f2f_homography;

PYBIND11_MODULE(homography_f2f, m)
{
  // This should wrap all of f2f_homography except for the (templated)
  // constructor directly from an Eigen::Matrix and the copy
  // constructor
  py::class_<f2f_homography, std::shared_ptr<f2f_homography>>(m, "F2FHomography")
    .def(py::init<kv::homography_sptr const&, kv::frame_id_t, kv::frame_id_t>())
    .def(py::init<kv::frame_id_t>())
    .def_property_readonly("homography", &f2f_homography::homography)
    .def_property_readonly("from_id", &f2f_homography::from_id)
    .def_property_readonly("to_id", &f2f_homography::to_id)
    .def("inverse", &f2f_homography::inverse)
    .def("__mul__", &f2f_homography::operator*)
    .def("get",
	 [] (f2f_homography const& self, int r, int c)
	 {
	   auto m = self.homography()->matrix();
	   if(0 <= r && r < m.rows() && 0 <= c && c < m.cols())
	   {
	     return m(r, c);
	   }
	   throw std::out_of_range("Tried to perform get() out of bounds");
	 },
	 "Convenience method that returns the underlying coefficient"
	 " at the given row and column")
    ;
}

ext: cxx

max_stars_repo_name: neal-siekierski/kwiver

max_issues_repo_name: neal-siekierski/kwiver

max_forks_repo_name: neal-siekierski/kwiver


text:
#include <assert.h>
#include <boost/algorithm/string/predicate.hpp>
#include <boost/program_options.hpp>
#include <fmt/core.h>
#include <iostream>
#include <string>
#include <time.h>

#include "run_dedoppler.h"

using namespace std;

namespace po = boost::program_options;

const string VERSION = "0.0.6";

// This method just handles command line parsing, and the real work is done
// via the dedoppler function.
int main(int argc, char* argv[]) {
  po::options_description desc("seticore options");
  desc.add_options()
    ("help,h", "produce help message")
    ("input", po::value<string>(), "alternate way of setting the input .h5 file")
    ("output", po::value<string>(), "the output .dat file. if not provided, uses the input filename but replaces its suffix with .dat")
    ("max_drift,M", po::value<double>()->default_value(10.0), "maximum drift in Hz/sec")
    ("min_drift,m", po::value<double>()->default_value(0.0001), "minimum drift in Hz/sec")
    ("snr,s", po::value<double>()->default_value(25.0), "minimum SNR to report a hit")
    ;

  po::positional_options_description p;
  p.add("input", -1);

  po::variables_map vm;
  po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm);
  po::notify(vm);

  if (!vm.count("input") || vm.count("help")) {
    cerr << "usage: seticore [input]\n";
    cerr << "seticore version: " << VERSION << endl;
    cerr << desc << "\n";
    return 1;
  }

  string input = vm["input"].as<string>();
  string output;
  if (!vm.count("output")) {
    // By default, output to a .dat file in the same location as the input file
    auto index = input.find_last_of(".");
    if (index >= input.size()) {
      cerr << "unrecognized input filename: " << input << endl;
      return 1;
    }
    output = input.substr(0, index) + ".dat";
  } else {
    output = vm["output"].as<string>();
  }

  double max_drift = vm["max_drift"].as<double>();
  double snr = vm["snr"].as<double>();
  double min_drift = vm["min_drift"].as<double>();

  cout << "welcome to seticore, version " << VERSION << endl;
  cout << "loading input from " << input << endl;
  cout << fmt::format("dedoppler parameters: max_drift={:.2f} min_drift={:.4f} snr={:.2f}\n",
                      max_drift, min_drift, snr);
  cout << "writing output to " << output << endl;
  int tstart = time(NULL);
  runDedoppler(input, output, max_drift, min_drift, snr);
  int tstop = time(NULL);
  cerr << fmt::format("dedoppler elapsed time {:d}s\n", tstop - tstart);
}

max_stars_repo_name: lacker/seticore

max_issues_repo_name: lacker/seticore

max_forks_repo_name: lacker/seticore