text:
//
//  Copyright (C) 2003-2017 Greg Landrum and Rational Discovery LLC
//
//   @@ All Rights Reserved @@
//  This file is part of the RDKit.
//  The contents are covered by the terms of the BSD license
//  which is included in the file license.txt, found at the root
//  of the RDKit source tree.
//
#include <GraphMol/RDKitBase.h>
#include <GraphMol/QueryOps.h>
#include <GraphMol/Rings.h>
#include <RDGeneral/types.h>
#include <boost/dynamic_bitset.hpp>
#include <set>

// introduced for the sake of efficiency
// this is the maximum ring size that will be considered
// as a candiate for fused-ring aromaticity. This is picked to
// be a bit bigger than the outer ring in a porphyrin
// This came up while fixing sf.net issue249
const unsigned int maxFusedAromaticRingSize = 24;

/****************************************************************
Here are some molecule that have trouble us aromaticity wise

1. We get 2 aromatic rings for this molecule, but daylight says none. Also
   if we replace [O+] with the N daylight says we have two aromatic rings
   [O-][N+]1=CC2=CC=C[O+]2[Cu]13[O+]4C=CC=C4C=[N+]3[O-]

2. We the fused ring system with all the rings in it is considered aroamtic by
our code.
   This is because we count electrons from burried atoms as well when
   we are dealing with fused rings

3. Here's a weird fused ring case, Pattern NAN, A overall
   O=C3C2=CC1=CC=COC1=CC2=CC=C3
 ************************************************************/
namespace RingUtils {
using namespace RDKit;

void pickFusedRings(int curr, const INT_INT_VECT_MAP &neighMap, INT_VECT &res,
                    boost::dynamic_bitset<> &done, int depth) {
  auto pos = neighMap.find(curr);
  PRECONDITION(pos != neighMap.end(), "bad argument");
  done[curr] = 1;
  res.push_back(curr);

  const INT_VECT &neighs = pos->second;
#if 0
    std::cerr<<"depth: "<<depth<<" ring: "<<curr<<" size: "<<res.size()<<" neighs: "<<neighs.size()<<std::endl;
    std::cerr<<"   ";
    std::copy(neighs.begin(),neighs.end(),std::ostream_iterator<int>(std::cerr," "));
    std::cerr<<"\n";
#endif
  for (int neigh : neighs) {
    if (!done[neigh]) {
      pickFusedRings(neigh, neighMap, res, done, depth + 1);
    }
  }
}

bool checkFused(const INT_VECT &rids, INT_INT_VECT_MAP &ringNeighs) {
  INT_INT_VECT_MAP_CI nci;
  int nrings = rdcast<int>(ringNeighs.size());
  boost::dynamic_bitset<> done(nrings);
  int rid;
  INT_VECT fused;

  // mark all rings in the system other than those in rids as done
  for (nci = ringNeighs.begin(); nci != ringNeighs.end(); nci++) {
    rid = (*nci).first;
    if (std::find(rids.begin(), rids.end(), rid) == rids.end()) {
      done[rid] = 1;
    }
  }

  // then pick a fused system from the remaining (i.e. rids)
  // If the rings in rids are fused we should get back all of them
  // in fused
  // if we get a smaller number in fused then rids are not fused
  pickFusedRings(rids.front(), ringNeighs, fused, done);

  CHECK_INVARIANT(fused.size() <= rids.size(), "");
  return (fused.size() == rids.size());
}

void makeRingNeighborMap(const VECT_INT_VECT &brings,
                         INT_INT_VECT_MAP &neighMap, unsigned int maxSize,
                         unsigned int maxOverlapSize) {
  int nrings = rdcast<int>(brings.size());
  int i, j;
  INT_VECT ring1;
  for (i = 0; i < nrings; i++) {
    INT_VECT neighs;
    neighMap[i] = neighs;
  }

  for (i = 0; i < nrings; i++) {
    if (maxSize && brings[i].size() > maxSize) continue;
    ring1 = brings[i];
    for (j = i + 1; j < nrings; j++) {
      if (maxSize && brings[j].size() > maxSize) continue;
      INT_VECT inter;
      Intersect(ring1, brings[j], inter);
      if (inter.size() > 0 &&
          (!maxOverlapSize || inter.size() <= maxOverlapSize)) {
        neighMap[i].push_back(j);
        neighMap[j].push_back(i);
      }
    }
  }
#if 0
    for (i = 0; i < nrings; i++) {
      std::cerr<<"**************\n    "<<i<<"\n*************\n";
      std::copy(neighMap[i].begin(),neighMap[i].end(),std::ostream_iterator<int>(std::cerr," "));
      std::cerr<<"\n";
    }
#endif
}

}  // end of namespace RingUtils

// local utility namespace
namespace {
using namespace RDKit;

typedef enum {
  VacantElectronDonorType,
  OneElectronDonorType,
  TwoElectronDonorType,
  OneOrTwoElectronDonorType,
  AnyElectronDonorType,
  NoElectronDonorType
} ElectronDonorType;  // used in setting aromaticity
typedef std::vector<ElectronDonorType> VECT_EDON_TYPE;
typedef VECT_EDON_TYPE::iterator VECT_EDON_TYPE_I;
typedef VECT_EDON_TYPE::const_iterator VECT_EDON_TYPE_CI;

/******************************************************************************
 * SUMMARY:
 *  Apply Huckel rule to the specified ring and mark the bonds and atoms in the
 *  ring accordingly.
 *
 * ARGUMENTS:
 *  mol - molecule of interest
 *  ring - list of atoms that form the simple ring
 *  edon - list of electron donor type of the atoms in the molecule
 *
 * RETURN :
 *  none
 *
 * ASSUMPTIONS:
 *  the ring has to a simple ring and the electron donor type of each atom
 *  in the ring have been computed
 ******************************************************************************/
static bool applyHuckel(ROMol &mol, const INT_VECT &ring,
                        const VECT_EDON_TYPE &edon);

static void applyHuckelToFused(
    ROMol &mol,                   // molecule of interets
    const VECT_INT_VECT &srings,  // list of all ring as atom IDS
    const VECT_INT_VECT &brings,  // list of all rings as bond ids
    const INT_VECT &fused,       // list of ring ids in the current fused system
    const VECT_EDON_TYPE &edon,  // eletron donar state for each atom
    INT_INT_VECT_MAP &ringNeighs,
    int &narom,  // number of aromatic ring so far
    unsigned int maxNumFusedRings, unsigned int minRingSize);

void markAtomsBondsArom(ROMol &mol, const VECT_INT_VECT &srings,
                        const VECT_INT_VECT &brings, const INT_VECT &ringIds,
                        std::set<unsigned int> &doneAtms) {
  INT_VECT aring, bring;
  INT_VECT_CI ri, ai, bi;

  // first mark the atoms in the rings
  for (ri = ringIds.begin(); ri != ringIds.end(); ri++) {
    aring = srings[*ri];

    // first mark the atoms in the ring
    // std::cerr << "aring:";
    for (ai = aring.begin(); ai != aring.end(); ai++) {
      // std::cerr << " " << *ai;
      mol.getAtomWithIdx(*ai)->setIsAromatic(true);
      doneAtms.insert(*ai);
    }
    // std::cerr << std::endl;
  }

  // mark the bonds
  // here we want to be more careful. We don't want to mark the fusing bonds
  // as aromatic - only the outside bonds in a fused system are marked aromatic.
  // - loop through the rings and count the number of times each bond appears in
  //   all the fused rings.
  // - bonds that appeard only once are marked aromatic
  INT_MAP_INT bndCntr;
  INT_MAP_INT_I bci;

  for (ri = ringIds.begin(); ri != ringIds.end(); ri++) {
    bring = brings[*ri];
    for (bi = bring.begin(); bi != bring.end(); bi++) {
      if (bndCntr.find(*bi) == bndCntr.end()) {
        bndCntr[*bi] = 1;
      } else {
        bndCntr[*bi] += 1;
      }
    }
  }
  // now mark bonds that have a count of 1 to be aromatic;
  // std::cerr << "bring:";
  for (bci = bndCntr.begin(); bci != bndCntr.end(); bci++) {
    // std::cerr << " " << bci->first << "(" << bci->second << ")";
    if ((*bci).second == 1) {
      Bond *bond = mol.getBondWithIdx(bci->first);
      bond->setIsAromatic(true);
      switch (bond->getBondType()) {
        case Bond::SINGLE:
        case Bond::DOUBLE:
          bond->setBondType(Bond::AROMATIC);
          break;
        default:
          break;
      }
    }
  }
  // std::cerr << std::endl;
}

void getMinMaxAtomElecs(ElectronDonorType dtype, int &atlw, int &atup) {
  switch (dtype) {
    case AnyElectronDonorType:
      atlw = 0;
      atup = 2;
      break;
    case OneOrTwoElectronDonorType:
      atlw = 1;
      atup = 2;
      break;
    case OneElectronDonorType:
      atlw = atup = 1;
      break;
    case TwoElectronDonorType:
      atlw = atup = 2;
      break;
    case NoElectronDonorType:
    case VacantElectronDonorType:
    default:
      atlw = atup = 0;
      break;
  }
}

bool incidentNonCyclicMultipleBond(const Atom *at, int &who) {
  PRECONDITION(at, "bad atom");
  // check if "at" has an non-cyclic multiple bond on it
  // if yes check which atom this bond goes to
  // and record the atomID in who
  const ROMol &mol = at->getOwningMol();
  ROMol::OEDGE_ITER beg, end;
  boost::tie(beg, end) = mol.getAtomBonds(at);
  while (beg != end) {
    if (!mol.getRingInfo()->numBondRings(mol[*beg]->getIdx())) {
      if (mol[*beg]->getValenceContrib(at) >= 2.0) {
        who = mol[*beg]->getOtherAtomIdx(at->getIdx());
        return true;
      }
    }
    ++beg;
  }
  return false;
}

bool incidentCyclicMultipleBond(const Atom *at) {
  PRECONDITION(at, "bad atom");
  const ROMol &mol = at->getOwningMol();
  ROMol::OEDGE_ITER beg, end;
  boost::tie(beg, end) = mol.getAtomBonds(at);
  while (beg != end) {
    if (mol.getRingInfo()->numBondRings(mol[*beg]->getIdx())) {
      if (mol[*beg]->getValenceContrib(at) >= 2.0) {
        return true;
      }
    }
    beg++;
  }
  return false;
}

bool incidentMultipleBond(const Atom *at) {
  PRECONDITION(at, "bad atom");
  int deg = at->getDegree() + at->getNumExplicitHs();
  ROMol::OEDGE_ITER beg, end;
  boost::tie(beg, end) = at->getOwningMol().getAtomBonds(at);
  while (beg != end) {
    Bond *bond = at->getOwningMol()[*beg];
    if (!std::lround(bond->getValenceContrib(at))) --deg;
    ++beg;
  }
  return at->getExplicitValence() != static_cast<int>(deg);
}

bool applyHuckel(ROMol &mol, const INT_VECT &ring, const VECT_EDON_TYPE &edon,
                 unsigned int minRingSize) {
  RDUNUSED_PARAM(mol);
  if (ring.size() < minRingSize) return false;
  int atlw, atup, rlw, rup, rie;
  bool aromatic = false;
  rlw = 0;
  rup = 0;
  for (auto idx : ring) {
    ElectronDonorType edonType = edon[idx];
    getMinMaxAtomElecs(edonType, atlw, atup);
    rlw += atlw;
    rup += atup;
  }

  if (rup >= 6) {
    for (rie = rlw; rie <= rup; rie++) {
      if ((rie - 2) % 4 == 0) {
        aromatic = true;
        break;
      }
    }
  } else if (rup == 2) {
    aromatic = true;
  }
#if 0
    std::cerr <<" ring: ";
    std::copy(ring.begin(),ring.end(),std::ostream_iterator<int>(std::cerr," "));
    std::cerr <<" rlw: "<<rlw<<" rup: "<<rup<<" aromatic? "<<aromatic<<std::endl;
#endif
  return aromatic;
}

void applyHuckelToFused(
    ROMol &mol,                   // molecule of interets
    const VECT_INT_VECT &srings,  // list of all ring as atom IDS
    const VECT_INT_VECT &brings,  // list of all rings as bond ids
    const INT_VECT &fused,       // list of ring ids in the current fused system
    const VECT_EDON_TYPE &edon,  // eletron donor state for each atom
    INT_INT_VECT_MAP &ringNeighs,  // list of neighbors for eac candidate ring
    int &narom,                    // number of aromatic ring so far
    unsigned int maxNumFusedRings, unsigned int minRingSize = 0) {
  // this function check huckel rule on a fused system it starts
  // with the individual rings in the system and then proceeds to
  // check for larger system i.e. if we have a 3 ring fused system,
  // huckel rule checked first on all teh 1 ring subsystems then 2
  // rung subsystems etc.

  INT_VECT aromRings;
  aromRings.resize(0);
  unsigned int nrings = rdcast<unsigned int>(fused.size());
  INT_VECT curRs;
  INT_VECT_CI mri;
  curRs.push_back(fused.front());
  int pos;
  unsigned int i, curSize = 0;
  INT_VECT comb;
  pos = -1;
  INT_VECT unionAtms;
  Union(srings, unionAtms);
  unsigned int nAtms = rdcast<unsigned int>(unionAtms.size());
  std::set<unsigned int> doneAtoms;
  while (1) {
    if (pos == -1) {
      curSize++;
      // check is we are done with all the atoms in the fused
      // system, if so quit. This is a fix for Issue252 REVIEW: is
      // this check sufficient or should we add an additional
      // constraint on the number of combinations of rings in a
      // fused system that we will try. The number of combinations
      // can obviously be quite large when the number of rings in
      // the fused system is large
      if ((doneAtoms.size() == nAtms) ||
          (curSize > std::min(nrings, maxNumFusedRings))) {
        break;
      }
      comb.resize(curSize);
      pos = 0;
      for (i = 0; i < curSize; i++) {
        comb[i] = i;
      }
    } else {
      pos = nextCombination(comb, nrings);
    }

    if (pos == -1) {
      continue;
    }

    curRs.resize(0);
    for (i = 0; i < comb.size(); i++) {
      curRs.push_back(fused[comb[i]]);
    }

    // check if the picked subsystem is fused
    if (ringNeighs.size() && !RingUtils::checkFused(curRs, ringNeighs)) {
      continue;
    }

    // check aromaticity on the current fused system
    INT_VECT exclude;
    for (i = 0; i < srings.size(); i++) {
      if (std::find(curRs.begin(), curRs.end(), static_cast<int>(i)) ==
          curRs.end()) {
        exclude.push_back(i);
      }
    }
    INT_VECT unon;
    Union(srings, unon, &exclude);

    if (applyHuckel(mol, unon, edon, minRingSize)) {
      // mark the atoms and bonds in these rings to be aromatic
      markAtomsBondsArom(mol, srings, brings, curRs, doneAtoms);

      // add the ring IDs to the aromatic rings found so far
      // avoid duplicates
      for (mri = curRs.begin(); mri != curRs.end(); mri++) {
        if (std::find(aromRings.begin(), aromRings.end(), (*mri)) ==
            aromRings.end()) {
          aromRings.push_back(*mri);
        }
      }
    }  // end check huckel rule
  }    // end while(1)
  narom += rdcast<int>(aromRings.size());
}

bool isAtomCandForArom(const Atom *at, const ElectronDonorType edon,
                       bool allowThirdRow = true, bool allowTripleBonds = true,
                       bool allowHigherExceptions = true, bool onlyCorN = false,
                       bool allowExocyclicMultipleBonds = true) {
  PRECONDITION(at, "bad atom");
  if (onlyCorN && at->getAtomicNum() != 6 && at->getAtomicNum() != 7)
    return false;
  if (!allowThirdRow && at->getAtomicNum() > 10) return false;

  // limit aromaticity to:
  //   - the first two rows of the periodic table
  //   - Se and Te
  if (at->getAtomicNum() > 18 &&
      (!allowHigherExceptions ||
       (at->getAtomicNum() != 34 && at->getAtomicNum() != 52))) {
    return false;
  }
  switch (edon) {
    case VacantElectronDonorType:
    case OneElectronDonorType:
    case TwoElectronDonorType:
    case OneOrTwoElectronDonorType:
    case AnyElectronDonorType:
      break;
    default:
      return (false);
  }

  // atoms that aren't in their default valence state also get shut out
  int defVal = PeriodicTable::getTable()->getDefaultValence(at->getAtomicNum());
  if (defVal > 0 && rdcast<int>(at->getTotalValence()) >
                        (PeriodicTable::getTable()->getDefaultValence(
                            at->getAtomicNum() - at->getFormalCharge()))) {
    return false;
  }

  // heteroatoms or charged carbons with radicals also disqualify us from being
  // considered. This was github issue 432 (heteroatoms) and 1936 (charged
  // carbons)
  if (at->getNumRadicalElectrons() &&
      (at->getAtomicNum() != 6 || at->getFormalCharge())) {
    return false;
  }

  // We are going to explicitly disallow atoms that have more
  // than one double or triple bond. This is to handle
  // the situation:
  //   C1=C=NC=N1 (sf.net bug 1934360)
  int nUnsaturations = at->getExplicitValence() - at->getDegree();
  if (nUnsaturations > 1) {
    unsigned int nMult = 0;
    const ROMol &mol = at->getOwningMol();
    ROMol::OEDGE_ITER beg, end;
    boost::tie(beg, end) = mol.getAtomBonds(at);
    while (beg != end) {
      switch (mol[*beg]->getBondType()) {
        case Bond::SINGLE:
        case Bond::AROMATIC:
          break;
        case Bond::DOUBLE:
          ++nMult;
          break;
        case Bond::TRIPLE:
          if (!allowTripleBonds) return false;
          ++nMult;
          break;
        default:
          // hopefully we had enough sense that we don't even
          // get here with these bonds... If we do land here,
          // just bail... I have no good answer for them.
          break;
      }
      if (nMult > 1) break;
      ++beg;
    }
    if (nMult > 1) return (false);
  }

  if (!allowExocyclicMultipleBonds) {
    const ROMol &mol = at->getOwningMol();
    ROMol::OEDGE_ITER beg, end;
    boost::tie(beg, end) = mol.getAtomBonds(at);
    while (beg != end) {
      const Bond *bnd = mol[*beg];
      if ((bnd->getBondType() == Bond::DOUBLE ||
           bnd->getBondType() == Bond::TRIPLE) &&
          !queryIsBondInRing(bnd))
        return false;
      ++beg;
    }
  }

  return (true);
}

ElectronDonorType getAtomDonorTypeArom(
    const Atom *at, bool exocyclicBondsStealElectrons = true) {
  PRECONDITION(at, "bad atom");
  if (at->getAtomicNum() == 0) {
    // dummies can be anything:
    return AnyElectronDonorType;
  }

  ElectronDonorType res = NoElectronDonorType;
  int nelec = MolOps::countAtomElec(at);
  int who = -1;
  const ROMol &mol = at->getOwningMol();
  if (nelec < 0) {
    res = NoElectronDonorType;
  } else if (nelec == 0) {
    if (incidentNonCyclicMultipleBond(at, who)) {
      // This is borderline:  no electron to spare but may have an empty
      // p-orbital
      // Not sure if this should return vacantElectronDonorType
      // FIX: explicitly doing this as a note for potential problems
      //
      res = VacantElectronDonorType;
    } else if (incidentCyclicMultipleBond(at)) {
      // no electron but has one in a in cycle multiple bond
      res = OneElectronDonorType;
    } else {
      // no multiple bonds no electrons
      res = NoElectronDonorType;
    }
  } else if (nelec == 1) {
    if (incidentNonCyclicMultipleBond(at, who)) {
      // the only available electron is going to be from the
      // external multiple bond this electron will not be available
      // for aromaticity if this atom is bonded to a more electro
      // negative atom
      const Atom *at2 = mol.getAtomWithIdx(who);
      if (exocyclicBondsStealElectrons &&
          PeriodicTable::getTable()->moreElectroNegative(at2->getAtomicNum(),
                                                         at->getAtomicNum())) {
        res = VacantElectronDonorType;
      } else {
        res = OneElectronDonorType;
      }
    } else {
      // require that the atom have at least one multiple bond
      if (incidentMultipleBond(at)) {
        res = OneElectronDonorType;
      }
      // account for the tropylium and cyclopropenyl cation cases
      else if (at->getFormalCharge() == 1) {
        res = VacantElectronDonorType;
      }
    }
  } else {
    if (incidentNonCyclicMultipleBond(at, who)) {
      // for cases with more than one electron :
      // if there is an incident multiple bond with an element that
      // is more electronegative than the this atom, count one less
      // electron
      const Atom *at2 = mol.getAtomWithIdx(who);
      if (exocyclicBondsStealElectrons &&
          PeriodicTable::getTable()->moreElectroNegative(at2->getAtomicNum(),
                                                         at->getAtomicNum())) {
        nelec--;
      }
    }
    if (nelec % 2 == 1) {
      res = OneElectronDonorType;
    } else {
      res = TwoElectronDonorType;
    }
  }
  return (res);
}
}  // namespace

namespace RDKit {
namespace MolOps {
int countAtomElec(const Atom *at) {
  PRECONDITION(at, "bad atom");

  // default valence :
  int dv = PeriodicTable::getTable()->getDefaultValence(at->getAtomicNum());
  if (dv <= 1) {
    // univalent elements can't be either aromatic or conjugated
    return -1;
  }

  // total atom degree:
  int degree = at->getDegree() + at->getTotalNumHs();

  ROMol::OEDGE_ITER beg, end;
  boost::tie(beg, end) = at->getOwningMol().getAtomBonds(at);
  while (beg != end) {
    Bond *bond = at->getOwningMol()[*beg];
    // don't count bonds that aren't actually contributing to the valence here:
    if (!std::lround(bond->getValenceContrib(at))) {
      --degree;
    }
    ++beg;
  }

  // if we are more than 3 coordinated we should not be aromatic
  if (degree > 3) {
    return -1;
  }

  // number of lone pair electrons = (outer shell elecs) - (default valence)
  int nlp;
  nlp = PeriodicTable::getTable()->getNouterElecs(at->getAtomicNum()) - dv;

  // subtract the charge to get the true number of lone pair electrons:
  nlp = std::max(nlp - at->getFormalCharge(), 0);

  int nRadicals = at->getNumRadicalElectrons();

  // num electrons available for donation into the pi system:
  int res = (dv - degree) + nlp - nRadicals;

  if (res > 1) {
    // if we have an incident bond with order higher than 2,
    // (e.g. triple or higher), we only want to return 1 electron
    // we detect this using the total unsaturation, because we
    // know that there aren't multiple unsaturations (detected
    // above in isAtomCandForArom())
    int nUnsaturations = at->getExplicitValence() - at->getDegree();
    if (nUnsaturations > 1) {
      res = 1;
    }
  }

  return res;
}

namespace {
int mdlAromaticityHelper(RWMol &mol, const VECT_INT_VECT &srings) {
  int narom = 0;
  // loop over all the atoms in the rings that can be candidates
  // for aromaticity
  // Atoms are candidates if
  //   - it is part of ring
  //   - has one or more electron to donate or has empty p-orbitals
  int natoms = mol.getNumAtoms();
  boost::dynamic_bitset<> acands(natoms);
  boost::dynamic_bitset<> aseen(natoms);
  VECT_EDON_TYPE edon(natoms);

  VECT_INT_VECT cRings;  // holder for rings that are candidates for aromaticity
  for (auto &sring : srings) {
    bool allAromatic = true;
    bool allDummy = true;
    for (auto firstIdx : sring) {
      Atom *at = mol.getAtomWithIdx(firstIdx);

      if (allDummy && at->getAtomicNum() != 0) {
        allDummy = false;
      }

      if (aseen[firstIdx]) {
        if (!acands[firstIdx]) allAromatic = false;
        continue;
      }
      aseen[firstIdx] = 1;

      // now that the atom is part of ring check if it can donate
      // electron or has empty orbitals. Record the donor type
      // information in 'edon' - we will need it when we get to
      // the Huckel rule later
      edon[firstIdx] = getAtomDonorTypeArom(at, false);
      // we only accept one electron donors?
      if (edon[firstIdx] != OneElectronDonorType) {
        allAromatic = false;
        continue;
      }
      const bool allowThirdRow = false;
      const bool allowTripleBonds = false;
      const bool allowHigherExceptions = false;
      const bool onlyCorN = true;
      const bool allowExocyclicMultipleBonds = false;
      acands[firstIdx] = isAtomCandForArom(
          at, edon[firstIdx], allowThirdRow, allowTripleBonds,
          allowHigherExceptions, onlyCorN, allowExocyclicMultipleBonds);
      if (!acands[firstIdx]) allAromatic = false;
    }
    if (allAromatic && !allDummy) {
      cRings.push_back(sring);
    }
  }

  // first convert all rings to bonds ids
  VECT_INT_VECT brings;
  RingUtils::convertToBonds(cRings, brings, mol);

  // make the neighbor map for the rings
  // i.e. a ring is a neighbor a another candidate ring if
  // shares at least one bond
  // useful to figure out fused systems
  INT_INT_VECT_MAP neighMap;
  RingUtils::makeRingNeighborMap(brings, neighMap, maxFusedAromaticRingSize, 1);

  // now loop over all the candidate rings and check the
  // huckel rule - of course paying attention to fused systems.
  INT_VECT doneRs;
  int curr = 0;
  int cnrs = rdcast<int>(cRings.size());
  boost::dynamic_bitset<> fusDone(cnrs);
  INT_VECT fused;
  while (curr < cnrs) {
    fused.resize(0);
    RingUtils::pickFusedRings(curr, neighMap, fused, fusDone);
    const unsigned int maxFused = 6;
    const unsigned int minRingSize = 6;
    applyHuckelToFused(mol, cRings, brings, fused, edon, neighMap, narom,
                       maxFused, minRingSize);

    int rix;
    for (rix = 0; rix < cnrs; rix++) {
      if (!fusDone[rix]) {
        curr = rix;
        break;
      }
    }
    if (rix == cnrs) {
      break;
    }
  }

  mol.setProp(common_properties::numArom, narom, true);

  return narom;
}

// use minRingSize=0 or maxRingSize=0 to ignore these constraints
int aromaticityHelper(RWMol &mol, const VECT_INT_VECT &srings,
                      unsigned int minRingSize, unsigned int maxRingSize,
                      bool includeFused) {
  int narom = 0;
  // loop over all the atoms in the rings that can be candidates
  // for aromaticity
  // Atoms are candidates if
  //   - it is part of ring
  //   - has one or more electron to donate or has empty p-orbitals
  int natoms = mol.getNumAtoms();
  boost::dynamic_bitset<> acands(natoms);
  boost::dynamic_bitset<> aseen(natoms);
  VECT_EDON_TYPE edon(natoms);

  VECT_INT_VECT cRings;  // holder for rings that are candidates for aromaticity
  for (auto &sring : srings) {
    size_t ringSz = sring.size();
    // test ring size:
    if ((minRingSize && ringSz < minRingSize) ||
        (maxRingSize && ringSz > maxRingSize))
      continue;

    bool allAromatic = true;
    bool allDummy = true;
    for (auto firstIdx : sring) {
      Atom *at = mol.getAtomWithIdx(firstIdx);

      if (allDummy && at->getAtomicNum() != 0) {
        allDummy = false;
      }

      if (aseen[firstIdx]) {
        if (!acands[firstIdx]) allAromatic = false;
        continue;
      }
      aseen[firstIdx] = 1;

      // now that the atom is part of ring check if it can donate
      // electron or has empty orbitals. Record the donor type
      // information in 'edon' - we will need it when we get to
      // the Huckel rule later
      edon[firstIdx] = getAtomDonorTypeArom(at);
      acands[firstIdx] = isAtomCandForArom(at, edon[firstIdx]);
      if (!acands[firstIdx]) allAromatic = false;
    }
    if (allAromatic && !allDummy) {
      cRings.push_back(sring);
    }
  }

  // first convert all rings to bonds ids
  VECT_INT_VECT brings;
  RingUtils::convertToBonds(cRings, brings, mol);

  if (!includeFused) {
    // now loop over all the candidate rings and check the
    // huckel rule - skipping fused systems
    INT_INT_VECT_MAP neighMap;
    for (size_t ri = 0; ri < cRings.size(); ++ri) {
      INT_VECT fused;
      fused.push_back(ri);
      const unsigned int maxFused = 6;
      const unsigned int minRingSize = 0;
      applyHuckelToFused(mol, cRings, brings, fused, edon, neighMap, narom,
                         maxFused, minRingSize);
    }
  } else {
    // make the neighbor map for the rings
    // i.e. a ring is a neighbor a another candidate ring if
    // shares at least one bond
    // useful to figure out fused systems
    INT_INT_VECT_MAP neighMap;
    RingUtils::makeRingNeighborMap(brings, neighMap, maxFusedAromaticRingSize, 1);

    // now loop over all the candidate rings and check the
    // huckel rule - of course paying attention to fused systems.
    INT_VECT doneRs;
    int curr = 0;
    int cnrs = rdcast<int>(cRings.size());
    boost::dynamic_bitset<> fusDone(cnrs);
    INT_VECT fused;
    while (curr < cnrs) {
      fused.resize(0);
      RingUtils::pickFusedRings(curr, neighMap, fused, fusDone);
      applyHuckelToFused(mol, cRings, brings, fused, edon, neighMap, narom, 6);

      int rix;
      for (rix = 0; rix < cnrs; rix++) {
        if (!fusDone[rix]) {
          curr = rix;
          break;
        }
      }
      if (rix == cnrs) {
        break;
      }
    }
  }

  mol.setProp(common_properties::numArom, narom, true);

  return narom;
}

}  // end of anonymous namespace

int setAromaticity(RWMol &mol, AromaticityModel model, int (*func)(RWMol &)) {
  // This function used to check if the input molecule came
  // with aromaticity information, assumed it is correct and
  // did not touch it. Now it ignores that information entirely.

  // first find the all the simple rings in the molecule
  VECT_INT_VECT srings;
  if (mol.getRingInfo()->isInitialized()) {
    srings = mol.getRingInfo()->atomRings();
  } else {
    MolOps::symmetrizeSSSR(mol, srings);
  }

  int res;
  switch (model) {
    case AROMATICITY_DEFAULT:
    case AROMATICITY_RDKIT:
      res = aromaticityHelper(mol, srings, 0, 0, true);
      break;
    case AROMATICITY_SIMPLE:
      res = aromaticityHelper(mol, srings, 5, 6, false);
      break;
    case AROMATICITY_MDL:
      res = mdlAromaticityHelper(mol, srings);
      break;
    case AROMATICITY_CUSTOM:
      PRECONDITION(
          func,
          "function must be set when aromaticity model is AROMATICITY_CUSTOM");
      res = func(mol);
      break;
    default:
      throw ValueErrorException("Bad AromaticityModel");
  }
  return res;
}

};  // end of namespace MolOps
};  // end of namespace RDKit

max_stars_repo_name: j-wags/rdkit

max_issues_repo_name: j-wags/rdkit

max_forks_repo_name: j-wags/rdkit


text: #include <boost/random/poisson_distribution.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
//==============================================================================
//         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
// cover for functor ellik in scalar mode
#include <nt2/elliptic/include/functions/ellipke.hpp>

#include <boost/simd/sdk/simd/io.hpp>
#include <cmath>
#include <iostream>
#include <nt2/include/constants/eps.hpp>
#include <nt2/include/functions/ellpk.hpp>
#include <nt2/include/functions/ellpe.hpp>
#include <nt2/include/constants/zero.hpp>
#include <nt2/include/functions/oneminus.hpp>
#include <nt2/sdk/meta/as_integer.hpp>
#include <nt2/sdk/unit/args.hpp>
#include <nt2/sdk/unit/module.hpp>
#include <nt2/sdk/unit/tests/cover.hpp>
#include <vector>

template < class T > struct ellipk1
{
  T operator()(const T & x)
  {
    T a1 = nt2::Zero<T>();
    T a0 = nt2::ellipke(x, nt2::Eps<T>(),a1);
    return a0;
  }
};

 template < class T > struct ellipk2
{
  T operator()(const T & x)
  {
    T a1 = nt2::Zero<T>();
    T a0 = nt2::ellipke(x, nt2::Eps<T>(), a1);
    return a1;
  }
};

NT2_TEST_CASE_TPL(ellpk_0,  NT2_SIMD_REAL_TYPES)
{
  using nt2::unit::args;
  const std::size_t NR = args("samples", NT2_NB_RANDOM_TEST);
  const double ulpd = args("ulpd", 2);

  typedef typename nt2::meta::as_integer<T>::type iT;
  const T min0 = args("min0", T(0));
  const T max0 = args("max0", T(1));
  std::cout << "Argument samples #0 chosen in range: [" << min0 << ",  " << max0 << "]" << std::endl;
  NT2_CREATE_BUF(a0,T, NR, min0, max0);

  std::vector<T> ref1(NR), ref2(NR);
  for(std::size_t i=0; i!=NR; ++i)
  {
    ref1[i] = nt2::ellpk(nt2::oneminus(a0[i]));
    ref2[i] = nt2::ellpe(nt2::oneminus(a0[i]));
  }

  NT2_COVER_FN_ULP_EQUAL(ellipk1<T>() , ((T, a0)), ref1, ulpd);
  NT2_COVER_FN_ULP_EQUAL(ellipk2<T>() , ((T, a0)), ref2, ulpd);

}

ext: cpp


text:
#ifndef TRIUMF_BNMR_SLR_GAUSS_DIST_EXP_HPP
#define TRIUMF_BNMR_SLR_GAUSS_DIST_EXP_HPP

#include <boost/math/constants/constants.hpp>
#include <cmath>
#include <triumf/bnmr/slr/common.hpp>
#include <triumf/math/faddeeva.hpp>
#include <iostream>

// TRIUMF: Canada's particle accelerator centre
namespace triumf {

// β-detected nuclear magnetic resonance (β-NMR)
namespace bnmr {

// spin-lattice relaxation (SLR)
namespace slr {

/// pulsed Gaussian distribution of exponentials integral
/// (from 0 to time_p <= time)
template <typename T = double>
T pulsed_gauss_dist_exp_integral(T time, T time_p, T nuclear_lifetime,
                                 T slr_rate, T sigma) {
  // imaginary number i
  std::complex<T> i(0.0, -1.0);
  // common exponential factor
  T expo =
      std::exp(-slr_rate / nuclear_lifetime / sigma / sigma -
               1.0 / 2.0 / sigma / sigma / nuclear_lifetime / nuclear_lifetime -
               slr_rate * slr_rate / 2.0 / sigma / sigma);
  std::cout << "expo = " << expo << "\n";
  // argument for the first error function
  std::complex<T> erf_arg1 =
      (((boost::math::constants::root_two<T>() * i * time -
         boost::math::constants::root_two<T>() * i * time_p) *
            sigma * sigma -
        boost::math::constants::root_two<T>() * i * slr_rate) *
           nuclear_lifetime -
       boost::math::constants::root_two<T>() * i) /
      (2.0 * sigma * nuclear_lifetime);
  std::cout << "erf_arg1 = " << erf_arg1 << "\n";
  // argument for the second error function
  std::complex<T> erf_arg2 =
      ((i * time * sigma * sigma - i * slr_rate) * nuclear_lifetime - i) /
      (boost::math::constants::root_two<T>() * sigma * nuclear_lifetime);
  std::cout << "erf_arg2 = " << erf_arg2 << "\n";
  // evaluate the result
  std::complex<T> result = boost::math::constants::root_pi<T>() * i * expo *
                           (triumf::math::faddeeva::erf(erf_arg1) -
                            triumf::math::faddeeva::erf(erf_arg2)) /
                           (boost::math::constants::root_two<T>() * sigma);
  std::cout << "result = " << result << "\n";
  // return only the real part
  return result.real();
}

/// pulsed Gaussian distribution of exponentials
template <typename T = double>
T pulsed_gauss_dist_exp(T time, T nuclear_lifetime, T pulse_length, T asymmetry,
                        T slr_rate, T sigma) {
  if (time == 0.0) {
    return asymmetry;
  } else if (time > 0.0 and time <= pulse_length) {
    return asymmetry *
           pulsed_gauss_dist_exp_integral(time, time, nuclear_lifetime,
                                          slr_rate, sigma) /
           normalization(time, nuclear_lifetime);
  } else if (time > pulse_length) {
    return (asymmetry *
            pulsed_gauss_dist_exp_integral(pulse_length, pulse_length,
                                           nuclear_lifetime, slr_rate, sigma) /
            normalization(pulse_length, nuclear_lifetime)) /
           std::exp(-(time - pulse_length) / nuclear_lifetime);
  } else {
    return 0.0;
  }
}

/// pulsed Gaussian distribution of exponentials (ROOT)
template <typename T = double>
T pulsed_gauss_dist_exp(const T *x, const T *par) {
  return pulsed_gauss_dist_exp<T>(*x, par[0], par[1], par[2], par[3], par[4]);
}

} // namespace slr

} // namespace bnmr

} // namespace triumf

#endif // TRIUMF_BNMR_SLR_GAUSS_DIST_EXP_HPP

max_stars_repo_name: rmlmcfadden/triumfpp

max_issues_repo_name: rmlmcfadden/triumfpp

max_forks_repo_name: rmlmcfadden/triumfpp


text:
// Copyright (c) 2009-2012 Bitcoin Developers
// Copyright (c) 2011-2012 Litecoin Developers
// Copyright (c) 2013-2014 Phoenixcoin Developers
// Distributed under the MIT/X11 software license, see the accompanying
// file LICENCE or http://www.opensource.org/licenses/mit-license.php

#include "init.h" // for pwalletMain
#include "bitcoinrpc.h"
#include "ui_interface.h"
#include "base58.h"

#include <boost/lexical_cast.hpp>

#define printf OutputDebugStringF

using namespace json_spirit;
using namespace std;

class CTxDump
{
public:
    CBlockIndex *pindex;
    int64 nValue;
    bool fSpent;
    CWalletTx* ptx;
    int nOut;
    CTxDump(CWalletTx* ptx = NULL, int nOut = -1)
    {
        pindex = NULL;
        nValue = 0;
        fSpent = false;
        this->ptx = ptx;
        this->nOut = nOut;
    }
};

Value importprivkey(const Array& params, bool fHelp)
{
    if (fHelp || params.size() < 1 || params.size() > 2)
        throw runtime_error(
            "importprivkey <phoenixcoinprivkey> [label]\n"
            "Adds a private key (as returned by dumpprivkey) to your wallet.");

    string strSecret = params[0].get_str();
    string strLabel = "";
    if (params.size() > 1)
        strLabel = params[1].get_str();
    CBitcoinSecret vchSecret;
    bool fGood = vchSecret.SetString(strSecret);

    if (!fGood) throw JSONRPCError(-5,"Invalid private key");

    CKey key;
    bool fCompressed;
    CSecret secret = vchSecret.GetSecret(fCompressed);
    key.SetSecret(secret, fCompressed);
    CKeyID vchAddress = key.GetPubKey().GetID();
    {
        LOCK2(cs_main, pwalletMain->cs_wallet);

        pwalletMain->MarkDirty();
        pwalletMain->SetAddressBookName(vchAddress, strLabel);

        if (!pwalletMain->AddKey(key))
            throw JSONRPCError(-4,"Error adding key to wallet");

        pwalletMain->ScanForWalletTransactions(pindexGenesisBlock, true);
        pwalletMain->ReacceptWalletTransactions();
    }

    return Value::null;
}

Value dumpprivkey(const Array& params, bool fHelp)
{
    if (fHelp || params.size() != 1)
        throw runtime_error(
	    "dumpprivkey <phoenixcoinaddress>\n"
	    "Reveals the private key corresponding to <phoenixcoinaddress>.");

    string strAddress = params[0].get_str();
    CBitcoinAddress address;
    if (!address.SetString(strAddress))
	throw JSONRPCError(-5, "Invalid Phoenixcoin address");
    CKeyID keyID;
    if (!address.GetKeyID(keyID))
        throw JSONRPCError(-3, "Address does not refer to a key");
    CSecret vchSecret;
    bool fCompressed;
    if (!pwalletMain->GetSecret(keyID, vchSecret, fCompressed))
        throw JSONRPCError(-4,"Private key for address " + strAddress + " is not known");
    return CBitcoinSecret(vchSecret, fCompressed).ToString();
}

/* Wallet integrity check */
Value checkwallet(const Array& params, bool fHelp) {

    if(fHelp || (params.size() > 0)) throw runtime_error(
      "checkwallet\n"
      "Wallet integrity check.\n");

    int nMismatchSpent;
    int nOrphansFound;
    int64 nBalanceInQuestion;
    pwalletMain->FixSpentCoins(nMismatchSpent, nOrphansFound, nBalanceInQuestion, true);
    Object result;
    if(!nMismatchSpent && !nOrphansFound) {
        result.push_back(Pair("wallet check passed", true));
    } else {
        if(nMismatchSpent) {
            result.push_back(Pair("mismatched spent coins", nMismatchSpent));
            result.push_back(Pair("amount in question", ValueFromAmount(nBalanceInQuestion)));
        }
        if(nOrphansFound) {
            result.push_back(Pair("orphans found", nOrphansFound));
        }
    }

    return(result);
}

/* Wallet repair (removal of failed and orphaned transactions) */
Value repairwallet(const Array& params, bool fHelp) {

    if(fHelp || (params.size() > 0)) throw runtime_error(
      "repairwallet\n"
      "Wallet repair if any mismatches found.\n");

    int nMismatchSpent;
    int nOrphansFound;
    int64 nBalanceInQuestion;
    pwalletMain->FixSpentCoins(nMismatchSpent, nOrphansFound, nBalanceInQuestion, false);
    Object result;
    if(!nMismatchSpent && !nOrphansFound) {
        result.push_back(Pair("wallet check passed", true));
    } else {
        if(nMismatchSpent) {
            result.push_back(Pair("mismatched spent coins", nMismatchSpent));
            result.push_back(Pair("amount affected by repair", ValueFromAmount(nBalanceInQuestion)));
        }
        if(nOrphansFound) {
            result.push_back(Pair("orphans removed", nOrphansFound));
        }
    }

    return(result);
}

/* Resend unconfirmed wallet transactions immediately */
Value resendtx(const Array& params, bool fHelp) {

    if(fHelp || (params.size() > 1)) throw runtime_error(
      "resendtx\n"
      "Resend unconfirmed transactions immediately.\n");

    ResendWalletTransactions(true);

    return(Value::null);
}

max_stars_repo_name: fossabot/Phoenixcoin

max_issues_repo_name: fossabot/Phoenixcoin

max_forks_repo_name: fossabot/Phoenixcoin


text:
//==============================================================================
//         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_DISPATCH_META_UNPROXY_HPP_INCLUDED
#define BOOST_DISPATCH_META_UNPROXY_HPP_INCLUDED

/*!
  @file
  @brief Defines and implements the unproxify function
**/

#include <boost/fusion/include/copy.hpp>
#include <boost/dispatch/meta/proxy.hpp>

namespace boost { namespace dispatch
{
  namespace meta
  {
    /*!
      @brief Metafunction computing the concrete type associated to a potential
      proxy type.
    **/
    template<class T> struct unproxy
    {
      template<class U> struct value_ { typedef typename U::value_type type; };

      typedef typename mpl::eval_if < is_proxy<T>
                                    , value_<T>
                                    , mpl::identity<T const&>
                                    >::type         type;
    };
  }

  /*!
    @brief Unproxify a value

    For any given value @c t of type @T, return either @c t if @T is not a
    proxy type or a value of the concrete type of @c T if @c T is a proxy type.

    @param t Value to unproxify
  **/
  template<class T>
  BOOST_FORCEINLINE typename meta::unproxy<T>::type unproxy(T const& t)
  {
    return t;
  }
} }

#endif


text:
#pragma once
#include <boost/program_options.hpp>
#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <math.h>
#include "image.pb.h"

/*
README
To avoid confusion in indexing:
  x goes along heigth
  y goes along width
*/

namespace utils {

  struct Options {
    // Command line arguments (with default values)
    unsigned num_ipus;
    unsigned num_iterations;
    float alpha;
    std::size_t height;
    std::size_t width;
    std::string vertex;
    std::string in_file;
    std::string out_file;
    bool cpu;
    // Other arguments which are a consequence of the code or environment
    std::size_t side;
    std::size_t tiles_per_ipu = 0;
    std::size_t num_tiles_available = 0;
    std::size_t halo_volume = 0;
    std::vector<std::size_t> splits = {0,0};
    std::vector<std::size_t> smallest_slice = {std::numeric_limits<size_t>::max(),1};
    std::vector<std::size_t> largest_slice = {0,0};
    denoise::Image image; // copy of image (if input jpg image is used)
  };

  inline
  Options parseOptions(int argc, char** argv) {
    Options options;
    namespace po = boost::program_options;
    po::options_description desc("Flags");
    // Construction of exactly cubic sub-grids
    desc.add_options()
    ("help", "Show command help.")
    (
      "num-ipus",
      po::value<unsigned>(&options.num_ipus)->default_value(1),
      "Number of IPUs to use."
    )
    (
      "num-iterations",
      po::value<unsigned>(&options.num_iterations)->default_value(1000),
      "PDE: number of iterations to execute on grid."
    )
    (
      "height",
      po::value<std::size_t>(&options.height)->default_value(8000),
      "Heigth of a custom 2D grid. Will be overwritten if --in-file is used."
    )
    (
      "width",
      po::value<std::size_t>(&options.width)->default_value(8000),
      "Width of a custom 2D grid. Will be overwritten if --in-file is used."
    )
    (
      "alpha",
      po::value<float>(&options.alpha)->default_value(0.1),
      "PDE: float. kappa * delta t / h^2"
    )
    (
      "cpu",
      po::bool_switch(&options.cpu)->default_value(false),
      "Also perform CPU execution to control results from IPU."
    )
    (
      "vertex",
      po::value<std::string>(&options.vertex)->default_value("HeatEquationOptimized"),
      "Name of vertex (from codelets.cpp) to use for the computation."
    )
    (
      "in-file",
      po::value<std::string>(&options.in_file)->default_value(""),
      "Filename of image (.bin): if empty string, then a random grid will be generated instead."
    )
    (
      "out-file",
      po::value<std::string>(&options.out_file)->default_value("denoised.jpg.bin"),
      "Output jpg file."
    ); // NOTE: remember to remove this semicolon if more options are added in future
    po::variables_map vm;
    po::store(po::parse_command_line(argc, argv, desc), vm);
    if (vm.count("help")) {
      std::cout << desc << "\n";
      throw std::runtime_error("Show help");
    }
    po::notify(vm);

    return options;
  }

} // End of namespace Utils

poplar::Device getDevice(unsigned numIpus) {
  /* return a Poplar device with the desired number of IPUs */
  auto manager = poplar::DeviceManager::createDeviceManager();
  auto devices = manager.getDevices(poplar::TargetType::IPU, numIpus);
  // Use the first available device
  for (auto &device : devices)
    if (device.attach())
      return std::move(device);

  throw std::runtime_error("No hardware device available.");
}

inline float randomFloat() {
  return static_cast <float> (rand() / static_cast <float> (RAND_MAX));
}

inline static unsigned index(unsigned x, unsigned y, unsigned width) {
  return y + (x*width);
}

inline static unsigned block_low(unsigned id, unsigned p, unsigned n) {
  return (id*n)/p;
}

inline static unsigned block_high(unsigned id, unsigned p, unsigned n) {
  return block_low(id+1, p, n);
}

inline static unsigned block_size(unsigned id, unsigned p, unsigned n) {
  return block_high(id, p, n) - block_low(id, p, n);
}

std::size_t area(std::vector<std::size_t> shape) {
  // return area of shape vector (2D)
  return shape[0]*shape[1];
}

void workDivision(utils::Options &options) {
  /* Function UPDATES options.splits
   * 1) all tiles will be used, hence options.num_tiles_available must
   *    previously be updated by using the target object
   * 2) choose the splits that minimizes the halo regions
   */
  std::size_t tile_count = options.tiles_per_ipu;
  std::size_t height = (options.height - 2); // Actual height (boundaries wont be computed)
  std::size_t width = (options.width - 2) / options.num_ipus; // Actual width
  float smallest_halo_region = std::numeric_limits<float>::infinity();

  // Try all unique combinations where i*j = tile_count
  for (std::size_t i = 1; i <= tile_count; ++i) {
    if ((tile_count % i) == 0) { // then i is a factor
      std::size_t j = tile_count / i; // and j must be the other factor
      std::size_t slice_height = height / i;
      std::size_t slice_width = width / j;
      std::size_t halo_region = 2.0*(slice_height + slice_width) - 4;
      if (halo_region <= smallest_halo_region) {
        smallest_halo_region = halo_region;
        options.splits[0] = i;
        options.splits[1] = j;
      }
    }
  }
}

std::vector<float> imageToVector(const denoise::Image& image) {
  /* Convert data from image into a vector<float> */
  std::vector<float> values(image.values().size());
  auto bytes = reinterpret_cast<const unsigned char *>(image.values().c_str());
  for (int i = 0; i < values.size(); ++i)
    values[i] = bytes[i]/255.0f;

  return values;
}

void vectorToImage(utils::Options &options, std::vector<float> &values) {
  /* Convert data from a vector<float> into image*/
  auto output = reinterpret_cast<unsigned char *>(&(*options.image.mutable_values())[0]);
  for (int i = 0; i < values.size(); ++i) {
    output[i] = 255.0f*values[i];
  }
}

std::vector<float> initializeValues(utils::Options &options) {
  if (options.in_file.length() != 0) {
    // Load jpg image, convert to vector, and update height and width
    std::fstream image_pb(options.in_file, std::ios::in | std::ios::binary);
    if (!options.image.ParseFromIstream(&image_pb))
      std::cerr << "Failed to read: " << options.in_file << std::endl;
    options.height = options.image.shape(1);
    options.width = options.image.shape(2);
    return imageToVector(options.image);
  } else {
    // Initialize random values
    std::vector<float> initial_values(options.height*options.width);
    for (std::size_t i = 0; i < options.height*options.width; ++i)
      initial_values[i] = randomFloat();
    return initial_values;
  }
}

void saveImage(utils::Options &options, std::vector<float> results) {
  vectorToImage(options, results);
  std::fstream denoised_pb(options.out_file, std::ios::out | std::ios::trunc | std::ios::binary);
  if (!options.image.SerializeToOstream(&denoised_pb)) {
    std::cerr << "Failed to write: " << options.out_file << "\n";
  }
}

std::vector<float> heatEquationCpu(
  const std::vector<float> initial_values,
  const utils::Options &options) {
  /*
   * Heat Equation Vertex, computed on the CPU, unparallelized.
   * The purpose of this function is to serve as the "true"
   * solution.
   * NOTE: can be very slow for large grids/large number of iterations
   * Approach:
   * for iter = 0, maxIter:
   *    - Heat Eqaution on inner elements of initalValues -> Results
   *    - Copy inner elements of results -> initial_values
   */
  const float gamma = 1.0 - 4.0*options.alpha;
  unsigned h = options.height;
  unsigned w = options.width;
  unsigned iter = options.num_iterations;
  std::vector<float> a(initial_values.size());
  std::vector<float> b(initial_values.size());

  // initial copy to include edges
  for (std::size_t i = 0; i < initial_values.size(); ++i)  {
    a[i] = initial_values[i];
    b[i] = initial_values[i];
  }

  // Heat Equation iterations
  for (std::size_t t = 0; t < iter; ++t) {
    for (std::size_t x = 1; x < h - 1; ++x) {
      for (std::size_t y = 1; y < w - 1; ++y) {
        a[index(x,y,w)] = gamma*b[index(x,y,w)] +
          options.alpha*(
            b[index(x+1,y,w)] +
            b[index(x-1,y,w)] +
            b[index(x,y+1,w)] +
            b[index(x,y-1,w)]
          );
      }
    }
    for (std::size_t x = 1; x < h - 1; ++x) {
      for (std::size_t y = 1; y < w - 1; ++y) {
        b[index(x,y,w)] = a[index(x,y,w)];
      }
    }
  }
  return a;
}

void printMeanSquaredError(
  std::vector<float> a,
  std::vector<float> b,
  utils::Options &options) {
  /*
   * Compute the MSE on only the inner elements of two 2D grids
   */
  double squared_error = 0, diff;
  std::size_t h = options.height;
  std::size_t w = options.width;
  for (std::size_t x = 1; x < h - 1; ++x) {
    for (std::size_t y = 1; y < w - 1; ++y) {
      diff = double(a[index(x,y,w)] - b[index(x,y,w)]);
      squared_error += diff*diff;
    }
  }
  double mean_squared_error = squared_error / (double) ((h-2)*(w-2));

  std::cout << "\nMean Squared Error (IPU vs. CPU) = " << mean_squared_error;
  if (mean_squared_error == double(0.0))
    std::cout << " (exactly)";
  std::cout << "\n";
}

void printResults(utils::Options &options, double wall_time) {

  // Calculate metrics
  double inner_area = (double)(options.height - 2)*(double)(options.width - 2);
  double flops_per_element = 6.0;
  double flops = inner_area*options.num_iterations*flops_per_element / wall_time;
  double internal_communication_ops = 2.0*(double)options.halo_volume*(double)options.num_ipus;
  double external_communication_ops = 4.0*(double)options.height*(options.num_ipus - 1.0); // 2 load and 2 stores of a height (partition along width)
  double bandwidth = (5.0*inner_area + internal_communication_ops + external_communication_ops)*(double)options.num_iterations*sizeof(float)/wall_time;
  double tflops = flops*1e-12;
  double bandwidth_TB_s = bandwidth*1e-12;

  std::cout << "2D Isotropic Diffusion"
    << "\n----------------------"
    << "\nVertex             = " << options.vertex
    << "\nNo. IPUs           = " << options.num_ipus
    << "\nNo. Tiles          = " << options.num_tiles_available
    << "\nTotal Grid         = " << options.height << "*" << options.width << " = "
                                 << options.height*options.width*1e-6 << " million elements"
    << "\nSmallest tile grid = " << options.smallest_slice[0] << "*" << options.smallest_slice[1]
    << "\nLargest tile grid  = " << options.largest_slice[0] << "*" << options.largest_slice[1]
    << "\nalpha              = " << options.alpha
    << "\nNo. Iterations     = " << options.num_iterations
    << "\n"
    << "\nLaTeX Tabular Row"
    << "\n-----------------"
    << "\nNo. IPUs & Grid & No. Iterations & Time [s] & Throughput [TFLOPS] & Minimum Bandwidth [TB/s] \\\\\n"
    << options.num_ipus << " & "
    << "$" << options.height << "\\times " << options.width << "$ & "
    << options.num_iterations << " & " << std::fixed
    << std::setprecision(2) << wall_time << " & "
    << std::setprecision(2) << tflops << " & "
    << std::setprecision(2) << bandwidth_TB_s << " \\\\"
    << "\n";
}

max_stars_repo_name: simehaa/IPU

max_issues_repo_name: simehaa/IPU

max_forks_repo_name: simehaa/IPU


text:
#pragma once

#include <Tanker/Crypto/Hash.hpp>
#include <Tanker/Crypto/PrivateSignatureKey.hpp>
#include <Tanker/Crypto/PublicEncryptionKey.hpp>
#include <Tanker/Crypto/PublicSignatureKey.hpp>
#include <Tanker/Crypto/SealedPrivateEncryptionKey.hpp>
#include <Tanker/Crypto/Signature.hpp>
#include <Tanker/Serialization/SerializedSource.hpp>
#include <Tanker/Trustchain/Actions/DeviceCreation/v1.hpp>
#include <Tanker/Trustchain/Actions/DeviceCreation/v3.hpp>
#include <Tanker/Trustchain/Actions/Nature.hpp>
#include <Tanker/Trustchain/Preprocessor/Actions/VariantImplementation.hpp>
#include <Tanker/Trustchain/UserId.hpp>

#include <boost/variant2/variant.hpp>
#include <nlohmann/json_fwd.hpp>

#include <cstddef>
#include <cstdint>
#include <utility>
#include <vector>

namespace Tanker
{
namespace Trustchain
{
namespace Actions
{
class DeviceCreation
{
  TANKER_TRUSTCHAIN_ACTION_VARIANT_IMPLEMENTATION(
      DeviceCreation,
      (DeviceCreation1, DeviceCreation3),
      (trustchainId, TrustchainId),
      (ephemeralPublicSignatureKey, Crypto::PublicSignatureKey),
      (userId, UserId),
      (delegationSignature, Crypto::Signature),
      (publicSignatureKey, Crypto::PublicSignatureKey),
      (publicEncryptionKey, Crypto::PublicEncryptionKey),
      (author, Crypto::Hash),
      (hash, Crypto::Hash),
      (signature, Crypto::Signature))

public:
  // v2 is missing, it's on purpose. We removed the "reset" feature, and v2 can
  // be converted into a v1 if the lastReset field is zero-filled.
  using v1 = DeviceCreation1;
  using v3 = DeviceCreation3;

  using DeviceType = v3::DeviceType;

  Nature nature() const;
  bool isGhostDevice() const;

  std::vector<std::uint8_t> delegationSignatureData() const;

private:
  friend std::uint8_t* to_serialized(std::uint8_t*, DeviceCreation const&);
  friend std::size_t serialized_size(DeviceCreation const&);
  friend void to_json(nlohmann::json&, DeviceCreation const&);
};

std::uint8_t* to_serialized(std::uint8_t* it, DeviceCreation const& dc);
std::size_t serialized_size(DeviceCreation const& dc);

void to_json(nlohmann::json& j, DeviceCreation const& dc);
}
}
}

max_stars_repo_name: TankerHQ/sdk-native

max_issues_repo_name: TankerHQ/sdk-native

max_forks_repo_name: TankerHQ/sdk-native


text:
//=======================================================================
// Copyright (c)
//
// Distributed under 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)
//=======================================================================
/**
 * @files zelikovsky_11_per_6_long_test.cpp
 * @brief
 * @author Piotr Wygocki
 * @version 1.0
 * @date 2013-02-04
 */

#include <boost/test/unit_test.hpp>

// this include must be here! //hack for clang
#include "paal/steiner_tree/zelikovsky_11_per_6.hpp"

#include "test_utils/logger.hpp"
#include "test_utils/read_steinlib.hpp"
#include "test_utils/test_result_check.hpp"

#include "paal/utils/irange.hpp"
#include "paal/data_structures/bimap.hpp"

#include <vector>

BOOST_AUTO_TEST_CASE(zelikovsky_11_per_6_test) {
    std::vector<paal::steiner_tree_test_with_metric> data;
    LOGLN("READING INPUT...");
    read_steinlib_tests(data);
    for (auto const &test : data) {
        LOGLN("TEST " << test.test_name);
        LOGLN("OPT " << test.optimal);

        using Metric = decltype(test.metric);
        using voronoiT = paal::data_structures::voronoi<Metric>;
        using FSet = typename voronoiT::GeneratorsSet;
        voronoiT voronoi(
            FSet(test.terminals.begin(), test.terminals.end()),
            FSet(test.steiner_points.begin(), test.steiner_points.end()),
            test.metric);
        std::vector<int> selected_steiner_points;
        paal::steiner_tree_zelikovsky11per6approximation(
            test.metric, voronoi, std::back_inserter(selected_steiner_points));

        auto res_range = boost::join(test.terminals, selected_steiner_points);
        LOG_COPY_RANGE_DEL(res_range, ",");
        paal::data_structures::bimap<int> idx;
        auto g = paal::data_structures::metric_to_bgl_with_index(
            test.metric, res_range, idx);
        std::vector<int> pm(res_range.size());
        boost::prim_minimum_spanning_tree(g, &pm[0]);
        auto idx_m = paal::data_structures::make_metric_on_idx(test.metric, idx);
        int res(0);
        for (int i : paal::irange(pm.size())) {
            if (pm[i] != i) {
                res += idx_m(i, pm[i]);
            }
        }
        check_result(res, test.optimal, 11. / 6);
    }
}

max_stars_repo_name: Kommeren/AA

max_issues_repo_name: Kommeren/AA

max_forks_repo_name: Kommeren/AA


text:
//==================================================================================================
/*!
  @file

  Defines the as_unsigned meta-function

  @copyright 2016 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_DETAIL_DISPATCH_META_AS_UNSIGNED_HPP_INCLUDED
#define BOOST_SIMD_DETAIL_DISPATCH_META_AS_UNSIGNED_HPP_INCLUDED

#include <boost/simd/detail/dispatch/meta/factory_of.hpp>
#include <boost/simd/detail/dispatch/meta/primitive_of.hpp>
#include <boost/simd/detail/dispatch/meta/is_natural.hpp>
#include <boost/simd/detail/nsm.hpp>

namespace boost { namespace dispatch
{
  namespace detail
  {
    template<typename T, typename Enable = void>
    struct as_unsigned
    {
      using f_t  = boost::dispatch::factory_of<T>;
      using p_t  = boost::dispatch::primitive_of_t<T>;
      using type = typename f_t::template apply<typename detail::as_unsigned<p_t>::type>::type;
    };

    template<typename T>
    struct as_unsigned<T, typename std::enable_if<std::is_fundamental<T>::value>::type>
    {
      template<typename U> struct id { using type = T;};
      using result  = std::conditional<is_natural<T>::value,std::make_unsigned<T>,id<T>>;
      using type    = typename result::type::type;
    };
  }

  /*!
    @ingroup group-generation
    @brief Turn a type into an equivalent unsigned type

    Converts any type which Primitive is Fundamental into a compatible unsigned type.

    @tparam T     Type to convert.
  **/
  template<typename T>
  struct as_unsigned : detail::as_unsigned<typename std::decay<T>::type>
  {};

  template<typename T> using as_unsigned_t= typename as_unsigned<T>::type;
} }

#endif

max_stars_repo_name: SylvainCorlay/pythran

max_issues_repo_name: SylvainCorlay/pythran

max_forks_repo_name: SylvainCorlay/pythran


text:
#include <ros/ros.h>

#include <vector>
#include <boost/algorithm/string.hpp>

#include <pulson_ros/P4xx/RCM/rcmIf.h>
#include <pulson_ros/P4xx/RCM/rcm.h>
#include <pulson_ros/P4xx/RCM/hostInterfaceCommon.h>
#include <pulson_ros/RangeArray.h>
#include <pulson_ros/RangeMeasurement.h>

void print_status(rcmMsg_GetStatusInfoConfirm *statusInfo);
void print_configuration(rcmConfiguration *config);
void error_check(int r, const char *msg);

int main(int argc, char **argv)
{

    ros::init(argc, argv, "reset");
    ros::NodeHandle nh("~");

    // parameters
    int node_id;
    std::string dev_path;

    nh.param("node_id", node_id, 100);
    nh.param("dev_path", dev_path, std::string("/dev/ttyACM0"));

    int r;

    // initialize interface
    r = rcmIfInit(rcmIfUsb, (char*) dev_path.c_str());
    error_check(r, "Initialization Failed");

    // put in idle mode during configuration
    r = rcmSleepModeSet(RCM_SLEEP_MODE_IDLE);
    error_check(r, "Time out waiting for sleep mode set");

    // make sure correct internal opmode is set
    r = rcmOpModeSet(RCM_OPMODE_RCM);
    error_check(r, "Timeout waiting for opmode set");

    // execute Built-In Test
    int status;
    r = rcmBit(&status);
    error_check(r, "Timeout waiting for test");
    error_check(status, "Built-In Test Failed");

    // get status information
    rcmMsg_GetStatusInfoConfirm statusInfo;
    r = rcmStatusInfoGet(&statusInfo);
    error_check(r, "Timeout waiting for status information");
    print_status(&statusInfo);

    // get configuration from RCM
    rcmConfiguration config;
    r = rcmConfigGet(&config);
    error_check(r, "Timeout waiting for configuration");
    print_configuration(&config);

    // default configuration
    config.nodeId = node_id;
    config.integrationIndex = 7;
    config.codeChannel = 0;
    config.electricalDelayPsA = 0;
    config.electricalDelayPsB = 0;
    config.flags = 0;
    config.txGain = RCM_TXGAIN_MAX;
    config.persistFlag = RCM_PERSIST_ALL;

    r = rcmConfigSet(&config);
    error_check(r, "Time out waiting for rcmConfig confirm");
    print_configuration(&config);

    // cleanup
    rcmIfFlush();
    rcmIfClose();

}

void print_status(rcmMsg_GetStatusInfoConfirm *statusInfo)
{
    char buf[1024];
    sprintf(buf, "\n\tStatus Information:\n");
    sprintf(buf, "%s\t\tPackage Version: %s\n", buf, statusInfo->packageVersionStr);
    sprintf(buf, "%s\t\tRCM Version: %d.%d Build %d\n", buf,
            statusInfo->appVersionMajor,
            statusInfo->appVersionMinor,
            statusInfo->appVersionBuild);
    sprintf(buf, "%s\t\tUWB Kernel Version: %d.%d Build %d\n", buf,
            statusInfo->uwbKernelVersionMajor,
            statusInfo->uwbKernelVersionMinor,
            statusInfo->uwbKernelVersionBuild);
    sprintf(buf, "%s\t\tFirmware Version: %x/%x/%x Ver. %X\n", buf,
            statusInfo->firmwareMonth,
            statusInfo->firmwareDay,
            statusInfo->firmwareYear,
            statusInfo->firmwareVersion);
    sprintf(buf, "%s\t\tSerial Number: %08X\n", buf, statusInfo->serialNum);
    sprintf(buf, "%s\t\tBoard Revision: %c\n", buf, statusInfo->boardRev);
    sprintf(buf, "%s\t\tTemperature: %.2f degC\n", buf, statusInfo->temperature/4.0);

    ROS_INFO("[rcm_node] %s", buf);

}

void print_configuration(rcmConfiguration *config)
{
    char buf[1024];
    sprintf(buf, "\n\tConfiguration:\n");
    sprintf(buf, "%s\t\tNode ID: %d\n", buf, config->nodeId);
    sprintf(buf, "%s\t\tIntegration Index: %d\n", buf, config->integrationIndex);
    sprintf(buf, "%s\t\tAntenna Mode: %d\n", buf, config->antennaMode);
    sprintf(buf, "%s\t\tCode Channel: %d\n", buf, config->codeChannel);
    sprintf(buf, "%s\t\tElectrical Delay PsA: %d\n", buf, config->electricalDelayPsA);
    sprintf(buf, "%s\t\tElectrical Delay PsB: %d\n", buf, config->electricalDelayPsB);
    sprintf(buf, "%s\t\tFlags: 0x%X\n", buf, config->flags);
    sprintf(buf, "%s\t\ttxGain: %d\n", buf, config->txGain);

    ROS_INFO("[rcm_node] %s", buf);
}

void error_check(int r, char const *msg)
{
    if (r != 0)
    {
        ROS_ERROR("[rcm_node] %s\n", msg);
        ros::requestShutdown();
    }

}


text:
//  Copyright 2010-2011 Vicente J. Botet Escriba
//  Distributed under the Boost Software License, Version 1.0.
//  See http://www.boost.org/LICENSE_1_0.txt

#define BOOST_CHRONO_VERSION 2

#include <iostream>
#include <boost/type_traits/is_same.hpp>
#include <boost/chrono/stopwatches/strict_stopwatch.hpp>
#include "../cycle_count.hpp"
#include <boost/chrono/stopwatches/reporters/stopwatch_reporter.hpp>
#include <boost/chrono/stopwatches/reporters/system_default_formatter.hpp>

#include <boost/chrono/chrono_io.hpp>
#include <boost/system/system_error.hpp>
#include <boost/detail/lightweight_test.hpp>

#if !defined(BOOST_NO_STATIC_ASSERT)
#define NOTHING ""
#endif

using namespace boost::chrono;

template <typename Reporter>
void check_invariants()
{
  typedef Reporter Stopwatch;
    BOOST_CHRONO_STATIC_ASSERT((boost::is_same<typename Stopwatch::rep, typename Stopwatch::clock::duration::rep>::value), NOTHING, ());
    BOOST_CHRONO_STATIC_ASSERT((boost::is_same<typename Stopwatch::period, typename Stopwatch::clock::duration::period>::value), NOTHING, ());
    BOOST_CHRONO_STATIC_ASSERT((boost::is_same<typename Stopwatch::duration, typename Stopwatch::clock::time_point::duration>::value), NOTHING, ());
    BOOST_CHRONO_STATIC_ASSERT(Stopwatch::is_steady == Stopwatch::clock::is_steady, NOTHING, ());
    BOOST_CHRONO_STATIC_ASSERT((boost::is_same<typename Stopwatch::stopwatch_type, strict_stopwatch<typename Stopwatch::clock> >::value), NOTHING, ());
}

template <typename Reporter>
void check_default_constructor()
{
  Reporter _;
}

struct file_line {
  file_line(const char* file, std::size_t line)
  : fmt("%1%[%2%] Elapsed time:")
  {
    fmt % file % line;
  }
  ~file_line()
  {
    std::cout << fmt;
  }
  boost::format fmt;

};

template <typename Reporter>
void check_file_line2()
{
  Reporter _("%1%\n");
  file_line fl(__FILE__, __LINE__);
  ex::sleep_for<typename Reporter::clock>(milliseconds(100));

}
template <typename Reporter>
void check_file_line()
{
  Reporter rp("%1%[%2%] Elapsed time: %3%\n");
  rp.format() % __FILE__ % __LINE__;

  ex::sleep_for<typename Reporter::clock>(milliseconds(100));

}

#if !defined BOOST_CHRONO_DONT_PROVIDE_HYBRID_ERROR_HANDLING
template <typename Reporter>
void check_constructor_ec()
{
  boost::system::error_code ec;
  Reporter _(ec);
  BOOST_TEST(ec.value()==0);
}

template <typename Reporter>
void check_constructor_throws()
{
  Reporter _(boost::throws());
}
#endif

template <typename Reporter>
void check_elapsed(bool check=true)
{
  Reporter sw;
  ex::sleep_for<typename Reporter::clock>(milliseconds(100));
  typename Reporter::duration d=sw.elapsed();
  std::cout << d << std::endl;
  if (check)
    BOOST_TEST(d >= milliseconds(100));
}

template <typename Reporter>
void check_report()
{
  Reporter sw;
  ex::sleep_for<typename Reporter::clock>(milliseconds(100));
  sw.report();
}

template <typename Clock>
void check_all(bool check=true)
{
  typedef stopwatch_reporter<strict_stopwatch<Clock> > Reporter;
  typedef stopwatch_reporter<strict_stopwatch<Clock>, elapsed_formatter > ReporterE;

  check_invariants<Reporter>();
  check_default_constructor<Reporter>();
#if !defined BOOST_CHRONO_DONT_PROVIDE_HYBRID_ERROR_HANDLING
  check_constructor_ec<Reporter>();
  check_constructor_throws<Reporter>();
#endif
  check_elapsed<Reporter>(check);
  check_report<Reporter>();
  check_file_line<ReporterE>();
}

int main()
{
  typedef strict_stopwatch<high_resolution_clock > Stopwatch;
  typedef basic_stopwatch_reporter_default_formatter<char, Stopwatch>::type Formatter;
  typedef stopwatch_reporter<Stopwatch> Reporter;
  static Formatter fmtr;

  Reporter _(fmtr);

  check_all<ex::cycle_count<1500> >(true);

#if 0
  std::cout << "high_resolution_clock=\n";

  check_all<high_resolution_clock>();
#ifdef BOOST_CHRONO_HAS_CLOCK_STEADY
  std::cout << "steady_clock=\n";
  check_all<steady_clock>(false);
#endif
  std::cout << "system_clock=\n";
  check_all<system_clock>(false);

#if defined(BOOST_CHRONO_HAS_THREAD_CLOCK)
  std::cout << "thread_clock=\n";
  check_all<thread_clock>(false);
#endif

#if defined(BOOST_CHRONO_HAS_PROCESS_CLOCKS)
  std::cout << "process_real_cpu_clock=\n";
  check_all<process_real_cpu_clock>(false);
  std::cout << "process_user_cpu_clock=\n";
  check_all<process_user_cpu_clock>(false);
  std::cout << "process_system_cpu_clock=\n";
  check_all<process_system_cpu_clock>(false);
  std::cout << "process_cpu_clock=\n";
  check_all<process_cpu_clock>(false);
#endif
#endif

  return boost::report_errors();
}

max_stars_repo_name: juslee/boost-svn

max_issues_repo_name: juslee/boost-svn

max_forks_repo_name: juslee/boost-svn


text:
/*

Copyright (c) 2005-2017, University of Oxford.
All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

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 the name of the University of Oxford nor the names of its
   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 COPYRIGHT HOLDER 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.

*/

#ifndef ABSTRACTSIMPLECELLCYCLEMODEL_HPP_
#define ABSTRACTSIMPLECELLCYCLEMODEL_HPP_

#include "ChasteSerialization.hpp"
#include "ClassIsAbstract.hpp"

#include <boost/serialization/base_object.hpp>

#include <vector>

#include "AbstractCellCycleModel.hpp"
#include "CellCyclePhases.hpp"
#include "SimulationTime.hpp"

/**
 * This class contains basic information to all cell-cycle models
 * that do NOT explicitly include distinct phases (G1, S, G2, M).
 */
class AbstractSimpleCellCycleModel : public AbstractCellCycleModel
{
private:

    /** Needed for serialization. */
    friend class boost::serialization::access;
    /**
     * Archive the object and its member variables.
     *
     * @param archive the archive
     * @param version the current version of this class
     */
    template<class Archive>
    void serialize(Archive & archive, const unsigned int version)
    {
        archive & boost::serialization::base_object<AbstractCellCycleModel>(*this);
        archive & mCellCycleDuration;
    }

protected:

    /**
     * Duration of cell cycle.
     * May be used as a mean duration for stochastic cell-cycle models.
     */
    double mCellCycleDuration;

    /**
     * Protected copy-constructor for use by CreateCellCycleModel.
     * The only way for external code to create a copy of a cell cycle model
     * is by calling that method, to ensure that a model of the correct subclass is created.
     * This copy-constructor helps subclasses to ensure that all member variables are correctly copied when this happens.
     *
     * This method is called by child classes to set member variables for a daughter cell upon cell division.
     * Note that the parent cell cycle model will have had ResetForDivision() called just before CreateCellCycleModel() is called,
     * so performing an exact copy of the parent is suitable behaviour. Any daughter-cell-specific initialisation
     * can be done in InitialiseDaughterCell().
     *
     * @param rModel the cell cycle model to copy.
     */
    AbstractSimpleCellCycleModel(const AbstractSimpleCellCycleModel& rModel);

public:

    /**
     * Default constructor - creates an AbstractSimpleCellCycleModel.
     */
    AbstractSimpleCellCycleModel();

    /**
     * Destructor.
     */
    virtual ~AbstractSimpleCellCycleModel();

    /**
     * See AbstractCellCycleModel::ResetForDivision()
     *
     * @return whether the cell is ready to divide (enter M phase).
     */
    virtual bool ReadyToDivide();

    /** See AbstractCellCycleModel::ResetForDivision() */
    virtual void ResetForDivision();

    /**
     * Overridden InitialiseDaughterCell() method.
     *
     * Set the new cell's cell cycle duration once it has been created after division.
     * This is by calling SetCellCycleDuration() defined in child classes.
     */
    void InitialiseDaughterCell();

    /** See AbstractPhaseBasedCellCycleModel::Initialise()
     *
     * Calls SetCellCycleDuration() defined in child classes.
     */
    virtual void Initialise();

    /**
     * This method is implemented in Subclasses to set the cell cycle duration of the cell.
     * It is called on the beginning of the simulation (intialisation) when initialising
     * daughter cells and when reseting parent cells
     *
     * It must be implemented in child classes.
     */
    virtual void SetCellCycleDuration() = 0;

    /**
     * @return mCellCycleDuration
     */
    double GetCellCycleDuration() const;

    /**
     * Outputs cell cycle model parameters to file.
     *
     * As this method is pure virtual, it must be overridden
     * in subclasses.
     *
     * @param rParamsFile the file stream to which the parameters are output
     */
    virtual void OutputCellCycleModelParameters(out_stream& rParamsFile)=0;
};

CLASS_IS_ABSTRACT(AbstractSimpleCellCycleModel)

#endif /*ABSTRACTSIMPLECELLCYCLEMODEL_HPP_*/

max_stars_repo_name: gonayl/Chaste

max_issues_repo_name: gonayl/Chaste

max_forks_repo_name: gonayl/Chaste