text:
#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE MyTest
#include<boost/test/unit_test.hpp>

#include <climits>
#include <boost/test/unit_test.hpp>

#include "SeqLib/BWAWrapper.h"
#include "SeqLib/BamReader.h"
#include "SeqLib/BamWriter.h"
#include "SeqLib/ReadFilter.h"
#include "SeqLib/FermiAssembler.h"
#include "SeqLib/SeqPlot.h"
#include "SeqLib/RefGenome.h"

#define GZBED "test_data/test.bed.gz"
#define GZVCF "test_data/test.vcf.gz"
#define SBAM "test_data/small.bam"
#define OBAM "test_data/small_out.bam"
#define OCRAM "test_data/small_out.cram"
#define HGREF "/seq/references/Homo_sapiens_assembly19/v1/Homo_sapiens_assembly19.fasta"
#define TREF "test_data/test_ref.fa"
#define OREF "tmp_output.fa"
#define BEDFILE "test_data/test.bed"
#define VCFFILE "test_data/test.vcf"
#define JSON1 "test_data/example4.json"

using namespace SeqLib::Filter;
using namespace SeqLib;

#include <fstream>
#include "SeqLib/BFC.h"

BOOST_AUTO_TEST_CASE( read_gzbed ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  SeqLib::GRC g(GZBED, br.Header());

  BOOST_CHECK_EQUAL(g.size(), 3);

  BOOST_CHECK_EQUAL(g[2].chr, 21);

  SeqLib::GRC v(GZVCF, br.Header());
  BOOST_CHECK_EQUAL(v.size(), 57);

  BOOST_CHECK_EQUAL(v[29].chr, 0);
}

BOOST_AUTO_TEST_CASE ( bfc ) {

  BFC b;

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  SeqLib::BamRecord rec;
  BamRecordVector brv, brv2;
  size_t count = 0;
  while(br.GetNextRecord(rec) && count++ < 10000)
    brv.push_back(rec);

  for (int i = 5000; i < brv.size(); ++i)
    brv2.push_back(brv[i]);

  count = 0;
  while(br.GetNextRecord(rec) && count++ < 10000)
    brv2.push_back(rec);

  std::ofstream orig("orig.fa");
  std::ofstream corr("corr.fa");

  for (auto& i : brv)
    orig << ">" << i.Qname() << std::endl << i.Sequence() << std::endl;

  // add the seqs
  for (auto& r : brv)
    b.AddSequence(r.Sequence().c_str(), r.Qualities().c_str(), r.Qname().c_str());

  b.Train();
  b.clear();

  //b.ErrorCorrectToTag(brv2, "KC");

  //UnalignedSequenceVector v;
  //b.GetSequences(v);

  // write to corrected
  //for (auto& i : v) {
  //  corr << ">" << i.Name << std::endl << i.Seq << std::endl;
  //}
  //orig.close();
  //corr.close();

  //
  //v.clear();
  //b.FilterUnique();
  //b.GetSequences(v);

  // do everything at once
  //b.TrainAndCorrect(brv2);

  // do everything in place
  //b.TrainCorrection(brv2);
  //b.ErrorCorrectInPlace(brv2);
}

BOOST_AUTO_TEST_CASE( correct_and_assemble ) {

  BFC b;

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  SeqLib::BamRecord rec;
  BamRecordVector brv, brv2;
  size_t count = 0;
  while(br.GetNextRecord(rec) && count++ < 10000)
    b.AddSequence(rec.Sequence().c_str(), rec.Qualities().c_str(), rec.Qname().c_str());

  b.Train();
  b.ErrorCorrect();

  float kcov = b.GetKCov();
  int   kmer = b.GetKMer();

  UnalignedSequenceVector v;

  v.clear();
  std::string seq, name;
  while (b.GetSequence(seq, name))
    v.push_back({name, seq});

  //std::ofstream filt("filt.fa");
  //for (auto& i : v) {
  //  filt << ">" << i.Name << std::endl << i.Seq << std::endl;
  //}
  //filt.close();

  FermiAssembler f;
  f.AddReads(v);
  f.DirectAssemble(kcov);

  // retrieve the contigs
  std::vector<std::string> contigs = f.GetContigs();

  std::ofstream cont("contigs.fa");
  size_t cc = 0;
  for (auto& i : f.GetContigs()) {
    ++cc;
    cont << ">" << cc << std::endl << i << std::endl;
  }
  cont.close();

}

BOOST_AUTO_TEST_CASE( header_check ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  // a BAM header check
  BOOST_CHECK_EQUAL(h.GetSequenceLength(0), 249250621);
  BOOST_CHECK_EQUAL(h.GetSequenceLength(3), 191154276);
  BOOST_CHECK_EQUAL(h.GetSequenceLength("1"), 249250621);
  BOOST_CHECK_EQUAL(h.GetSequenceLength("4"), 191154276);
  BOOST_CHECK_EQUAL(h.GetSequenceLength("d4"), -1);
  BOOST_CHECK_EQUAL(h.GetSequenceLength(10000), -1);

  BOOST_CHECK_EQUAL(h.GetHeaderSequenceVector().size(), h.NumSequences());
  BOOST_CHECK_EQUAL(h.GetHeaderSequenceVector().begin()->Length, 249250621);

  SeqLib::BamRecord rec;
  size_t count = 0;

  while(br.GetNextRecord(rec) && count++ < 10000) {

  }

}

BOOST_AUTO_TEST_CASE( merge ) {

  SeqLib::GRC grc;
  // add two more that we know of
  grc.add(SeqLib::GenomicRegion(23, 10,100));
  grc.add(SeqLib::GenomicRegion(23, 20,110));

  grc.add(SeqLib::GenomicRegion(2, 10,100));
  grc.add(SeqLib::GenomicRegion(2, 20,110));
  grc.add(SeqLib::GenomicRegion(2, 200,310));

  grc.MergeOverlappingIntervals();
  BOOST_CHECK_EQUAL(grc.size(), 3);
  BOOST_CHECK_EQUAL(grc[0].chr, 2);
  BOOST_CHECK_EQUAL(grc[1].chr, 2);
  BOOST_CHECK_EQUAL(grc[2].chr, 23);
  BOOST_CHECK_EQUAL(grc[2].pos2, 110);
  BOOST_CHECK_EQUAL(grc[2].pos1, 10);
}

BOOST_AUTO_TEST_CASE ( interval_queries ) {

  SeqLib::GRC grc;

  // create a large GRC
  for (int i = 0; i < 10; ++i) {
    int chr = rand() % 23;
    int pos = rand() % 10000;
    grc.add(SeqLib::GenomicRegion(chr, pos, pos + 100));
  }
  grc.MergeOverlappingIntervals();

  // add two more that we know of
  grc.add(SeqLib::GenomicRegion(23, 10,100));
  grc.add(SeqLib::GenomicRegion(23, 20,110));

  // create the interval tree
  grc.CreateTreeMap();

  SeqLib::GRC results = grc.FindOverlaps(SeqLib::GenomicRegion(23, 10, 100), true);

  for (auto& i : results)
    std::cerr << " GRC overlaps results " << i << std::endl;

  BOOST_CHECK_EQUAL(results.size(), 2);
  BOOST_CHECK_EQUAL(results[1].pos2, 100);

  grc.MergeOverlappingIntervals();
  grc.CreateTreeMap();

  for(auto& r : grc)
    std::cerr << r << std::endl;

  std::vector<int32_t> q, s;
  results = grc.FindOverlaps(grc, q, s, true);

  std::cerr << " results.size " << results.size() << " Input size " << grc.size() << std::endl;
  BOOST_CHECK_EQUAL(results.size(), grc.size());
  BOOST_CHECK_EQUAL(results.TotalWidth(), grc.TotalWidth());

}

BOOST_AUTO_TEST_CASE( json_parse_from_file ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  std::string rules = "{\"global\" : {\"!anyflag\" : 1536}, \"\" : { \"rules\" : [{\"ic\" : true}, {\"clip\" : 5}, {\"ins\" : true}, {\"del\" : true}, {\"mapped\": true , \"mate_mapped\" : false}, {\"mate_mapped\" : true, \"mapped\" : false}]}}";
  ReadFilterCollection rfc(rules, br.Header());

  std::cerr << rfc << std::endl;

  SeqLib::BamRecord rec;
  size_t count = 0;

  while(br.GetNextRecord(rec) && count++ < 10000) {
    if (!rfc.isValid(rec))
      continue;
    // test global flag rule
    if ( (rec.QCFailFlag() || rec.DuplicateFlag())) {
      std::cerr << rec << std::endl;
      assert(false);
    }
  }

  /// direct from string
  ReadFilterCollection rfc2(rules, br.Header());

  while(br.GetNextRecord(rec) && count++ < 10000) {
    if (!rfc.isValid(rec))
      continue;
    // test global flag rule
    if ( (rec.QCFailFlag() || rec.DuplicateFlag())) {
      std::cerr << rec << std::endl;
      assert(false);
    }
  }

  // check that a bad key throws error
  //rules = "{\"global\" : {\"!anyflagf\" : 1536}, \"\" : { \"rules\" : [{\"ic\" : true}, {\"clip\" : 5}, {\"ins\" : true}, {\"del\" : true}, {\"mapped\": true , \"mate_mapped\" : false}, {\"mate_mapped\" : true, \"mapped\" : false}]}}";
  //BOOST_CHECK_THROW(ReadFilterCollection rfc2(rules, br.Header()), std::invalid_argument);
  // bad value, expected int
  //rules = "{\"global\" : {\"!anyflag\" : \"BAD\"}, \"\" : { \"rules\" : [{\"ic\" : true}, {\"clip\" : 5}, {\"ins\" : true}, {\"del\" : true}, {\"mapped\": true , \"mate_mapped\" : false}, {\"mate_mapped\" : true, \"mapped\" : false}]}}";
  //BOOST_CHECK_THROW(ReadFilterCollection rfc3(rules, br.Header()), std::invalid_argument);
  // bad JSON itself
  rules = "{\"global\" : \"!anyflag\" : 1536}, \"\" : { \"rules\" : [{\"ic\" : true}, {\"clip\" : 5}, {\"ins\" : true}, {\"del\" : true}, {\"mapped\": true , \"mate_mapped\" : false}, {\"mate_mapped\" : true, \"mapped\" : false}]}}";
  BOOST_CHECK_THROW(ReadFilterCollection rfc4(rules, br.Header()), std::invalid_argument);
  // bad value, expected range
  rules = "{\"global\" : {\"!anyflag\" : 1536}, \"\" : { \"rules\" : [{\"isize\" : \"BAD\"}]}}";
  BOOST_CHECK_THROW(ReadFilterCollection rfc4(rules, br.Header()), std::invalid_argument);

}

BOOST_AUTO_TEST_CASE( sw_alignment ) {

  const std::string ref = "ACTGCGAGCGACTAGCTCGTAGCTAGCTAGCTAGCTAGTGACTGCGGGCGATCATCGATCTTTTATTATCGCGATCGCTACGAC";
  const std::string seq = "ACTGCGAGCGACTAGCTCGTAGCTAGCTAGCTAGCTAGTGACTGCGGGCGATCATCGATCTTTTATTATCGCGATCGCTACGAC";
  //const std::string seq =                "CTCGTAGCTAGCTGCTAGCTAGTGACTGCGGGCGATCATCGATCTTTTATTATCGCG";
  const SeqLib::GenomicRegion gr(0,0,0);
  SeqLib::BamRecord b("test_name", seq, ref, &gr);

  std::cerr << " SMITH WATERMAN " << std::endl;
  std::cerr << b << std::endl;
}

BOOST_AUTO_TEST_CASE( read_filter_1 ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  SeqLib::GRC g;
  g.add(SeqLib::GenomicRegion(h.Name2ID("X"), 1100000, 1800000));
  g.CreateTreeMap();

  // make a new rule set
  ReadFilterCollection rfc;

  // make a new filter region
  ReadFilter rf;

  // add an isize rule on whole-genome
  AbstractRule ar;
  ar.isize = Range(200, 600, false); // 200 to 600, not inverted
  ar.mapq  = Range(10, 50, false); // 200 to 600, not inverted
  ar.nm    = Range(1, 1, false); // 200 to 600, not inverted
  rf.AddRule(ar);

  rf.setRegions(g);

  // add to the filter collection
  rfc.AddReadFilter(rf);

  SeqLib::GRC gback = rfc.getAllRegions();
  BOOST_CHECK_EQUAL(gback.size(), g.size());
  for (size_t i = 0; i < gback.size(); ++i)
    assert(g[i] == gback[i]);

  // display
  std::cerr << br.PrintRegions() << std::endl;

  // read / filter the reads
  SeqLib::BamRecord rec;
  size_t count = 0;

  while(br.GetNextRecord(rec) && count++ < 10000) {

    if (!rfc.isValid(rec))
      continue;

    // test isize rule
    if (!(rec.FullInsertSize() >= 200 || rec.FullInsertSize() <= 600)) {
      std::cerr << rec.FullInsertSize() << std::endl;
      assert(false);
    }
    // test mapq rule
    if (!(rec.MapQuality() >= 10 || rec.MapQuality() <= 50)) {
      std::cerr << rec.MapQuality() << std::endl;
      assert(false);
    }
    // test nm rule
    int32_t nm;
    rec.GetIntTag("NM", nm);
    if (nm == 1)
      assert(false);
  }
}

BOOST_AUTO_TEST_CASE ( fermi_add_reads ) {

  SeqLib::FermiAssembler f;

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamRecord r;
  size_t count = 0;
  while (br.GetNextRecord(r) && count++ < 1000)
    f.AddRead(r);
  while (br.GetNextRecord(r) && count++ < 2000)
    f.AddRead(SeqLib::UnalignedSequence(r.Qname(), r.Sequence(), r.Qualities()));

  f.CorrectReads();
  f.PerformAssembly();
  std::vector<std::string> out = f.GetContigs();

}

BOOST_AUTO_TEST_CASE ( seq_utils ) {

  // add commas
  BOOST_CHECK_EQUAL(SeqLib::AddCommas(1),"1");
  BOOST_CHECK_EQUAL(SeqLib::AddCommas(1000000),"1,000,000");

  // percent calc
  BOOST_CHECK_EQUAL(SeqLib::percentCalc(10,100), 10);
  BOOST_CHECK_EQUAL(SeqLib::percentCalc(7,8), 87);
  BOOST_CHECK_EQUAL(SeqLib::percentCalc(9,10), 90);
  BOOST_CHECK_EQUAL(SeqLib::percentCalc(2,3), 66);

  // scrub string
  BOOST_CHECK_EQUAL(SeqLib::scrubString("chr1", "chr"), "1");
  BOOST_CHECK_EQUAL(SeqLib::scrubString("chr1", ""), "chr1");
  BOOST_CHECK_EQUAL(SeqLib::scrubString("chr1", "dd"), "chr1");
  BOOST_CHECK_EQUAL(SeqLib::scrubString("chr1", "1"), "chr");

}

BOOST_AUTO_TEST_CASE( bam_record ) {

  // get a record
  SeqLib::BamReader br, br2;

  // try multiple open
  std::vector<std::string> bs = {"test_data/small.bam", "test_data/small.bam"};
  BOOST_CHECK(!br2.Open(bs)); // should be false, no dups
  bs = {"test_data/small.bam", "test_data/small.cram"};
  BOOST_CHECK(br.Open(bs)); // should be true
  SeqLib::BamRecord r;

  SeqLib::BamRecordVector brv;

  size_t count = 0;
  br.GetNextRecord(r);

  BOOST_CHECK_EQUAL(r.AsGenomicRegion().chr, 0);
  BOOST_CHECK_EQUAL(r.AsGenomicRegion().pos1,9995);
  BOOST_CHECK_EQUAL(r.AsGenomicRegion().pos2,10075);
  BOOST_CHECK_EQUAL(r.AsGenomicRegion().strand,'-');

  BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().chr, 15);
  BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().pos1,67300983);
  BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().pos2,67301134);
  BOOST_CHECK_EQUAL(r.AsGenomicRegionMate().strand,'-');

  BOOST_CHECK_EQUAL(std::floor(r.MeanPhred()), 15);

  BOOST_CHECK_EQUAL(r.CountNBases(), 0);

  r.SetQname("testq");
  BOOST_CHECK_EQUAL(r.Qname(), "testq");

  const std::string s = "ACTGCTAGCTAGCTACTCTGCTACTATATTAGCGCGCATTCGC";
  r.SetSequence(s);
  BOOST_CHECK_EQUAL(r.Sequence(), s);

  r.SmartAddTag("ST", "1");
  r.SmartAddTag("ST", "3");
  r.SmartAddTag("ST", "5");
  r.SmartAddTag("S2", "5");

  BOOST_CHECK_EQUAL(r.GetSmartIntTag("ST").size(), 3);
  BOOST_CHECK_EQUAL(r.GetSmartIntTag("ST").at(2), 5);
  BOOST_CHECK_EQUAL(r.GetSmartDoubleTag("ST").size(), 3);
  BOOST_CHECK_EQUAL(r.GetSmartDoubleTag("ST").at(2), 5.0);

  BOOST_CHECK_EQUAL(r.GetSmartStringTag("ST").size(), 3);
  BOOST_CHECK_EQUAL(r.GetSmartStringTag("ST")[1], "3");

  BOOST_CHECK_EQUAL(r.GetSmartDoubleTag("S2").at(0), 5.0);
  BOOST_CHECK_EQUAL(r.GetSmartIntTag("S2").at(0), 5);
}

BOOST_AUTO_TEST_CASE( fermi_assemble ) {

  SeqLib::FermiAssembler f;

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamRecord r;

  SeqLib::BamRecordVector brv;

  size_t count = 0;
  while(br.GetNextRecord(r) && count++ < 1000) {
    brv.push_back(r);
  }

  f.AddReads(brv);

  f.CorrectReads();

  SeqLib::UnalignedSequenceVector reads = f.GetSequences();
  BOOST_CHECK_EQUAL(reads.size(), brv.size());

  for (int i = 0; i < reads.size(); ++i) {
    if (brv[i].Sequence() != reads[i].Seq) {
      std::cerr << "************" << std::endl;
      std::cerr << brv[i].Sequence() << std::endl;
      std::cerr << reads[i].Seq << std::endl;
    }
  }

  // peform the assembly
  std::cerr << "...performing assembly" << std::endl;
  f.PerformAssembly();

  // retrieve the contigs
  std::vector<std::string> contigs = f.GetContigs();

}

BOOST_AUTO_TEST_CASE( bam_header_stdout ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  std::cout << h.AsString() << std::endl;
}

BOOST_AUTO_TEST_CASE( bam_header_name2id ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  BOOST_CHECK_EQUAL(h.Name2ID("2"), 1);
  BOOST_CHECK_EQUAL(h.Name2ID("23"), -1);

}

BOOST_AUTO_TEST_CASE( bam_header_id2name ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  BOOST_CHECK_EQUAL(h.IDtoName(2), "3");
  BOOST_CHECK_THROW(h.IDtoName(100), std::out_of_range);
  BOOST_CHECK_THROW(h.IDtoName(-1), std::invalid_argument);
  BOOST_CHECK_THROW(SeqLib::BamHeader().IDtoName(1), std::out_of_range);
}

BOOST_AUTO_TEST_CASE( genomic_ranges_string_constructor) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  const std::string in = "2:1,000,000-2,000,000";
  SeqLib::GenomicRegion gr(in, h);
  BOOST_CHECK_EQUAL(gr.chr, 1);
  BOOST_CHECK_EQUAL(gr.pos1, 1000000);
  BOOST_CHECK_EQUAL(gr.pos2, 2000000);

  BOOST_CHECK_THROW(SeqLib::GenomicRegion(in, SeqLib::BamHeader()), std::invalid_argument);

  BOOST_CHECK_EQUAL(gr.ChrName(h), "2");
  BOOST_CHECK_EQUAL(gr.ChrName(SeqLib::BamHeader()), "2");
  gr.chr = 1000;
  BOOST_CHECK_THROW(gr.ChrName(h), std::invalid_argument);

}

BOOST_AUTO_TEST_CASE( genomic_region_less_than ) {

  SeqLib::GenomicRegion gr1(0, 1, 2);
  SeqLib::GenomicRegion gr2(1, 1, 2);
  SeqLib::GenomicRegion gr3(1, 2, 2);
  SeqLib::GenomicRegion gr4(1, 6, 6);

  BOOST_CHECK(gr1 < gr2);
  BOOST_CHECK(gr2 > gr1);
  BOOST_CHECK(!(gr1 > gr2));

  BOOST_CHECK(gr2 < gr3);
  BOOST_CHECK(gr3 > gr2);
  BOOST_CHECK(!(gr2 > gr3));

  BOOST_CHECK(gr3 < gr4);
  BOOST_CHECK(!(gr4 == gr3));
  BOOST_CHECK(!(gr3 > gr4));
  BOOST_CHECK(gr4 > gr3);

  BOOST_CHECK(!(gr1 < gr1));
  BOOST_CHECK(!(gr1 > gr1));

  BOOST_CHECK(!(gr1 != gr1));
  BOOST_CHECK(gr2 != gr1);
  BOOST_CHECK(gr3 != gr1);
  BOOST_CHECK(gr4 != gr3);

  BOOST_CHECK(gr1 >= gr1);
  BOOST_CHECK(gr2 >= gr2);
  BOOST_CHECK(gr3 >= gr3);
  BOOST_CHECK(gr4 >= gr4);

  BOOST_CHECK(gr1 <= gr1);
  BOOST_CHECK(gr2 <= gr2);
  BOOST_CHECK(gr3 <= gr3);
  BOOST_CHECK(gr4 <= gr4);

  BOOST_CHECK(gr1 <= gr2);
  BOOST_CHECK(gr2 >= gr1);

  BOOST_CHECK(gr2 <= gr3);
  BOOST_CHECK(gr3 >= gr2);

}

BOOST_AUTO_TEST_CASE( genomic_region_distance ) {

  SeqLib::GenomicRegion gr1(0, 10, 100);
  SeqLib::GenomicRegion gr2(0, 10, 200);
  SeqLib::GenomicRegion gr3(1, 10, 100);
  SeqLib::GenomicRegion gr4(0, 100, 100);

  BOOST_CHECK_EQUAL(gr1.DistanceBetweenEnds(gr3), -1);
  BOOST_CHECK_EQUAL(gr1.DistanceBetweenEnds(gr1), 0);
  BOOST_CHECK_EQUAL(gr1.DistanceBetweenEnds(gr2), 100);
  BOOST_CHECK_EQUAL(gr1.DistanceBetweenEnds(gr4), 0);

  BOOST_CHECK_EQUAL(gr1.DistanceBetweenStarts(gr3), -1);
  BOOST_CHECK_EQUAL(gr1.DistanceBetweenStarts(gr1), 0);
  BOOST_CHECK_EQUAL(gr1.DistanceBetweenStarts(gr2), 0);
  BOOST_CHECK_EQUAL(gr1.DistanceBetweenStarts(gr4), 90);

}

BOOST_AUTO_TEST_CASE( small_trie_from_file) {

  AbstractRule ar;
  const bool inverted = false;
  ar.addMotifRule("test_data/motif.txt", inverted);

  ReadFilterCollection rfc;
  ReadFilter rf;
  rf.AddRule(ar);
  rfc.AddReadFilter(rf);

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  SeqLib::BamRecord rec;
  bool rule;
  size_t count = 0;
  while (br.GetNextRecord(rec) && count++ < 1000){
  }

}

BOOST_AUTO_TEST_CASE( large_trie ) {

  const std::string dictionary = "ACTG";

  const int string_size = 20;
  const int string_count = 10000;

  SeqLib::Filter::AhoCorasick aho;

  std::vector<std::string> k;

  std::cerr << "...generating key" << std::endl;
  for (int i = 0; i < string_count; ++i) {
    char* c = (char*) malloc(string_size + 1);
    for (int j = 0; j < string_size; ++j)
      c[j] = dictionary.at(rand() % 4);
    c[string_size] = '\0';
    k.push_back(std::string(c));
    free(c);
  }
  std::cerr << "...done with key" << std::endl;

  std::cerr << "...generating trie" << std::endl;
  for (auto& i : k)
    aho.AddMotif(i);
  std::cerr << "...done generating trie" << std::endl;

  std::cerr << "...querying trie" << std::endl;
  auto result = aho.aho_trie->parse_text(k[0]);
  std::cerr << "...querying trie fast" << std::endl;
  for (int i = 0; i < string_count; ++i) {
    //if (i % 20000 == 0)
    //  std::cerr << "... " << i << std::endl;
    auto result = aho.aho_trie->parse_text(k[i]);
  }

}

BOOST_AUTO_TEST_CASE( genomic_region_constructors ) {

  // GenomicRegion Constructors
  SeqLib::GenomicRegion gr(0, 0, 10, '+');
  BOOST_CHECK_EQUAL(gr.Width(), 11);

  SeqLib::GenomicRegion gr_empty;
  BOOST_TEST(gr_empty.IsEmpty());

  SeqLib::GenomicRegion gr2("chrX", "0", "10", SeqLib::BamHeader());
  BOOST_CHECK_EQUAL(gr2.Width(), 11);
  BOOST_CHECK_EQUAL(gr2.chr, 22);

  SeqLib::GenomicRegion gr3("X", "0", "10", SeqLib::BamHeader());
  BOOST_TEST(gr2 == gr3);

  BOOST_CHECK_THROW(SeqLib::GenomicRegion gr3("X", "a", "10", SeqLib::BamHeader()), std::invalid_argument);
  BOOST_CHECK_THROW(SeqLib::GenomicRegion gr3("X", "1000000000000000000000000000000000000000000000000000000000000000000000000000000", "10", SeqLib::BamHeader()), std::out_of_range);

  BOOST_CHECK_EQUAL(gr.DistanceBetweenStarts(gr2), -1);
  BOOST_CHECK_EQUAL(gr2.DistanceBetweenStarts(gr), -1);

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  BOOST_CHECK_EQUAL(SeqLib::GenomicRegion("X","1","100", br.Header()).chr, 22);

  // check negative inputs
  SeqLib::GenomicRegion grn(-1,-11,-10);
  BOOST_CHECK_EQUAL(grn.chr, -1);
  BOOST_CHECK_EQUAL(grn.pos1, -11);
  BOOST_CHECK_EQUAL(grn.pos2, -10);

  // check strand constructions
  SeqLib::GenomicRegion gra(0,0,0);
  SeqLib::GenomicRegion grb(0,10000,10001, '+');
  SeqLib::GenomicRegion grc(0,0,3, '-');
  BOOST_CHECK_EQUAL(gra.strand, '*');
  BOOST_CHECK_EQUAL(grb.strand, '+');
  BOOST_CHECK_EQUAL(grc.strand, '-');

  // check point string
  BOOST_CHECK_EQUAL(grb.PointString(), "1:10,000(+)");

  // check pretty string
  std::stringstream ss;
  ss << grb;
  BOOST_CHECK_EQUAL(ss.str(), "1:10,000-10,001(+)");

}

BOOST_AUTO_TEST_CASE( genomic_region_bad_inputs ) {

  BOOST_CHECK_THROW(SeqLib::GenomicRegion(0, 10, 9), std::invalid_argument);

  BOOST_CHECK_THROW(SeqLib::GenomicRegion(0,0,0,'P'), std::invalid_argument);

}

// BOOST_AUTO_TEST_CASE( genomic_region_random ) {

//   SeqLib::GenomicRegion gr;
//   std::srand(42);
//   gr.random();
//   BOOST_CHECK_EQUAL(gr.pointString(), "9:69,477,830(*)");

// }

BOOST_AUTO_TEST_CASE( genomic_region_range_operations ) {

  SeqLib::GenomicRegion gr(0,1,10);
  SeqLib::GenomicRegion gr2(0,1,11);
  gr.Pad(3);
  gr2.Pad(-3);
  BOOST_CHECK_EQUAL(gr.pos1,-2);
  BOOST_CHECK_EQUAL(gr.pos2,13);
  BOOST_CHECK_EQUAL(gr2.pos1,4);
  BOOST_CHECK_EQUAL(gr2.pos2,8);

  BOOST_CHECK_THROW(gr.Pad(-10), std::out_of_range);

}

BOOST_AUTO_TEST_CASE( genomic_check_overlaps ) {

  SeqLib::GenomicRegion gr1(0, 0, 10, '+');
  SeqLib::GenomicRegion gr2(1, 0, 10, '+');

  SeqLib::GenomicRegion gr3(0, 10, 20, '+');
  SeqLib::GenomicRegion gr4(1, 4, 10, '+');

  SeqLib::GenomicRegion gr5(1, 11, 12, '+');

  // partial overlaps should be one
  BOOST_CHECK_EQUAL(gr1.GetOverlap(gr3), 1);

  // argument contained gets 2
  BOOST_CHECK_EQUAL(gr2.GetOverlap(gr4), 2);

  // object contained gets 3
  BOOST_CHECK_EQUAL(gr4.GetOverlap(gr2), 3);

  // same chr, no overlap
  BOOST_CHECK_EQUAL(gr4.GetOverlap(gr5), 0);
  BOOST_CHECK_EQUAL(gr5.GetOverlap(gr4), 0);

}

BOOST_AUTO_TEST_CASE( bwa_wrapper ) {

  SeqLib::BWAWrapper bwa;

  // Set some options
  bwa.SetGapOpen(32);
  bwa.SetGapExtension(1);
  bwa.SetMismatchPenalty(18);
  bwa.SetAScore(2);
  bwa.SetZDropoff(100);
  bwa.Set3primeClippingPenalty(5);
  bwa.Set5primeClippingPenalty(5);
  bwa.SetBandwidth(1000);
  bwa.SetReseedTrigger(1.5);

  BOOST_CHECK_THROW(bwa.SetGapOpen(-1), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.SetGapExtension(-1), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.SetMismatchPenalty(-18), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.SetAScore(-2), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.SetZDropoff(-100), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.Set3primeClippingPenalty(-5), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.Set5primeClippingPenalty(-5), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.SetBandwidth(-1000), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.SetReseedTrigger(-1.5), std::invalid_argument);

  // no index loaded yet
  BOOST_CHECK_THROW(bwa.ChrIDToName(1), std::runtime_error);

  // load a test index
  BOOST_TEST(SeqLib::read_access_test(TREF));
  bwa.LoadIndex(TREF);

  BOOST_CHECK_EQUAL(bwa.NumSequences(), 2);

  BOOST_CHECK_EQUAL(bwa.ChrIDToName(0), "ref1");
  BOOST_CHECK_EQUAL(bwa.ChrIDToName(1), "ref2");
  BOOST_CHECK_THROW(bwa.ChrIDToName(2), std::out_of_range);

  BOOST_CHECK(!bwa.LoadIndex("test_data/small.bam"));

  SeqLib::BamHeader hh = bwa.HeaderFromIndex();
  BOOST_CHECK_EQUAL(hh.NumSequences(), 2);

  // error check the index construction
  SeqLib::UnalignedSequenceVector usv_bad1, usv_bad2, usv;;
  usv_bad1.push_back(SeqLib::UnalignedSequence("ref1","ACATGGCGAGCACTTCTAGCATCAGCTAGCTACGATCGATCGATCGATCGTAGC", std::string()));
  usv_bad1.push_back(SeqLib::UnalignedSequence("ref4", std::string(), std::string()));
  usv_bad1.push_back(SeqLib::UnalignedSequence("ref5","CGATCGTAGCTAGCTGATGCTAGAAGTGCTCGCCATGT", std::string()));
  usv_bad2.push_back(SeqLib::UnalignedSequence(std::string(), "ACATGGCGAGCACTTCTAGCATCAGCTAGCTACGATCGATCGATCGATCGTAGC", std::string()));
  usv_bad2.push_back(SeqLib::UnalignedSequence("ref4","ACCATCGCAGCAGCTATCTATTATATCGGCAGCATCTAGC", std::string()));
  usv_bad2.push_back(SeqLib::UnalignedSequence("ref5","CGATCGTAGCTAGCTGATGCTAGAAGTGCTCGCCATGT", std::string()));
  BOOST_CHECK_THROW(bwa.ConstructIndex(usv_bad1), std::invalid_argument);
  BOOST_CHECK_THROW(bwa.ConstructIndex(usv_bad2), std::invalid_argument);

  // construct a normal index
  usv.push_back(SeqLib::UnalignedSequence("ref3","ACATGGCGAGCACTTCTAGCATCAGCTAGCTACGATCGATCGATCGATCGTAGC", std::string()));
  usv.push_back(SeqLib::UnalignedSequence("ref4","CTACTTTATCATCTACACACTGCCTGACTGCGGCGACGAGCGAGCAGCTACTATCGACT", std::string()));
  usv.push_back(SeqLib::UnalignedSequence("ref5","CGATCGTAGCTAGCTGATGCTAGAAGTGCTCGCCATGT", std::string()));
  usv.push_back(SeqLib::UnalignedSequence("ref6","TATCTACTGCGCGCGATCATCTAGCGCAGGACGAGCATC" + std::string(100,'N') + "CGATCGTTATTATCGAGCGACGATCTACTACGT", std::string()));

  bwa.ConstructIndex(usv);

  BOOST_CHECK_EQUAL(bwa.NumSequences(), 4);
  bwa.ChrIDToName(1);

  BOOST_CHECK_THROW(bwa.ChrIDToName(-1), std::out_of_range);
  BOOST_CHECK_THROW(bwa.ChrIDToName(10000), std::out_of_range);

  std::cout << bwa.ChrIDToName(2) << std::endl;

  BOOST_CHECK_EQUAL(bwa.ChrIDToName(0), "ref3");
  BOOST_CHECK_EQUAL(bwa.ChrIDToName(1), "ref4");
  BOOST_CHECK_EQUAL(bwa.ChrIDToName(2), "ref5");
  BOOST_CHECK_EQUAL(bwa.ChrIDToName(3), "ref6");
  BOOST_CHECK_THROW(bwa.ChrIDToName(4), std::out_of_range);

  // write the index
  BOOST_CHECK(bwa.WriteIndex(OREF));

  // write the fasta
  std::ofstream os;
  os.open(OREF);
  os          << ">" << usv[0].Name << std::endl << usv[0].Seq <<
    std::endl << ">" << usv[1].Name << std::endl << usv[1].Seq <<
    std::endl << ">" << usv[2].Name << std::endl << usv[2].Seq <<
    std::endl << ">" << usv[3].Name << std::endl << usv[3].Seq <<
    std::endl;
  os.close();

  // read it back
  BOOST_CHECK(bwa.LoadIndex(OREF));

  // check that its good
  BOOST_CHECK_EQUAL(bwa.ChrIDToName(0), "ref3");
  BOOST_CHECK_EQUAL(bwa.ChrIDToName(1), "ref4");

  // try aligning a sequence
  std::cerr << "...aligning sequences" << std::endl;
  SeqLib::BamRecordVector brv, brv2;
  bool hardclip = false;
  bwa.AlignSequence("ACATGGCGAGCACTTCTAGCATCAGCTAGCTACGATCG", "name", brv, 0.9, hardclip, 1);
  // reverse complement
  bwa.AlignSequence("CGATCGTAGCTAGCTGATGCTAGAAGTGCTCGC", "name", brv2, 0.9, hardclip, 2);

  BOOST_CHECK_EQUAL(brv[0].Qname(), "name");
  BOOST_CHECK_EQUAL(brv[0].ChrID(), 2);
  BOOST_CHECK_EQUAL(brv[0].Sequence(), "CGATCGTAGCTAGCTGATGCTAGAAGTGCTCGCCATGT");
  std::cerr << " brv[0].GetCigar() " << brv[0].GetCigar() << std::endl;
  BOOST_CHECK_EQUAL(brv[0].GetCigar()[0].Type(), 'M');
  BOOST_CHECK_EQUAL(brv[0].GetCigar()[0].Length(), 38);

  // check from iterator
  SeqLib::Cigar ccc = brv[0].GetCigar();

  //SeqLib::Cigar::const_iterator f = brv[0].GetCigar().begin();
  BOOST_CHECK_EQUAL(ccc.begin()->Length(), 38);

  // check that it got both alignments
  BOOST_CHECK_EQUAL(brv2.size(), 2);

  // print info
  std::cerr << bwa << std::endl;
}

BOOST_AUTO_TEST_CASE( bam_reader ) {

  // print empty read
  std::cerr << SeqLib::BamRecord() << std::endl;

  SeqLib::BamReader bw;
  bw.Open(SBAM);

  // open index
  bw.SetRegion(SeqLib::GenomicRegion(22, 1000000, 1001000));

  // make a set of locations
  SeqLib::GRC grc;
  for (size_t i = 0; i < 24; ++i)
    grc.add(SeqLib::GenomicRegion(i, 1, 100));

  // set regions
  bw.SetMultipleRegions(grc);

  // write index of new bam
  // should print a warning since no write bam is specified
  //bw.BuildIndex();

  // open an output BAM
  //bw.OpenWriteBam(OBAM);

  // set tags to strip
  //bw.setStripTags("OQ,BI");

  // loop through and grab some reads
  SeqLib::BamRecord r;
  size_t count = 0;
  while (bw.GetNextRecord(r)) {
    //if (++count % 10 == 0)
    //  bw.WriteRecord(r);
  }

  // display info about BAM
  std::cerr << bw << std::endl;

  // write index of new bam
  //bw.BuildIndex();

  // reset the walker
  bw.Reset();

  // set a smaller region
  bw.SetRegion(SeqLib::GenomicRegion(0, 1, 100));
  std::cerr << bw << std::endl;

  // write as a cram
  //bw.OpenWriteBam(OCRAM);

  //
  //bw.setCram(OCRAM, HGREF);

  // print cram writer
  //std::cerr << bw << std::endl;
  // write the CRAM
  //while (bw.GetNextRecord(r, rule)) {
  //  if (++count % 10 == 0) {
  //    std::cerr << count << std::endl;
  //    bw.WriteRecord(r);
  //  }
  //}

}

BOOST_AUTO_TEST_CASE( set_qualities ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  SeqLib::BamRecord r;
  while (br.GetNextRecord(r)) {
    r.SetQualities("", 0);
    BOOST_CHECK_EQUAL(r.Qualities().at(0), '!');
    r.SetQualities(std::string(r.Length(), '#'), 33);
    BOOST_CHECK_EQUAL(r.Qualities(), std::string(r.Length(), '#'));
    BOOST_CHECK_THROW(r.SetQualities(std::string(8, '#'), 0), std::invalid_argument);
    break;
  }
}

BOOST_AUTO_TEST_CASE( header_constructor ) {

  SeqLib::HeaderSequenceVector hsv;
  hsv.push_back(SeqLib::HeaderSequence("1", 1000));
  hsv.push_back(SeqLib::HeaderSequence("chr2", 1200));
  SeqLib::BamHeader hdr(hsv);

}

BOOST_AUTO_TEST_CASE( overlapping_coverage ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamRecordVector brv;
  size_t count = 0;
  SeqLib::BamRecord r;
  while(br.GetNextRecord(r) && ++count < 4) {
    std::cout << " r " << r << std::endl;
    brv.push_back(r);
  }
  BOOST_CHECK_EQUAL(brv[0].OverlappingCoverage(brv[2]), 78);

}

BOOST_AUTO_TEST_CASE( gr_chr_region_set) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");

  SeqLib::GenomicRegion gr("1", br.Header());
  BOOST_CHECK_EQUAL(gr.chr, 0);
  BOOST_CHECK_EQUAL(gr.pos2, 249250621);
  BOOST_CHECK_EQUAL(gr.pos1, 1);

  BOOST_CHECK_THROW(SeqLib::GenomicRegion gr2("-1", br.Header()), std::invalid_argument);

}

BOOST_AUTO_TEST_CASE( sequtils ) {

  std::string seq = "actgACGTnTCN";

  SeqLib::rcomplement(seq);

  BOOST_CHECK_EQUAL(seq, "NGAnACGTcagt");

}

BOOST_AUTO_TEST_CASE( bam_write ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  SeqLib::BamRecord rec;

  // empty constructor
  SeqLib::BamWriter w;

  BOOST_CHECK(!w.WriteHeader());
  BOOST_CHECK(!w.Close());
  BOOST_CHECK(!w.BuildIndex());
  BOOST_CHECK(!w.WriteRecord(rec));

  w.Open("tmp_out.bam");

  // check that set CRAM fails
  BOOST_CHECK(!w.SetCramReference("dummy"));
  BOOST_CHECK(!w.WriteHeader());

  w.SetHeader(h);

  w.WriteHeader();

  size_t count = 0;

  while(br.GetNextRecord(rec) && count++ < 10000)
    w.WriteRecord(rec);

  BOOST_CHECK(!w.BuildIndex());
  w.Close();

  w.BuildIndex();

  // check that write header now fails
  BOOST_CHECK(!w.WriteHeader());

  // check that set CRAM fails
  BOOST_CHECK(!w.SetCramReference("badref"));

  // print some info
  std::cerr << w << std::endl;

}

BOOST_AUTO_TEST_CASE( bam_record_more ) {

  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamHeader h = br.Header();

  SeqLib::BamRecord rec;
  size_t count = 0;

  while(br.GetNextRecord(rec) && count++ < 100) {
    rec.ClearSeqQualAndTags();
    assert(rec.Sequence().empty());
    assert(rec.Qualities().empty());
    //int32_t nm;
    //rec.GetIntTag("NM", nm);
    //assert(!nm);
    std::string xa;
    rec.GetZTag("XA", xa);
    BOOST_CHECK_EQUAL(xa.empty(), rec.CountBWASecondaryAlignments()==0);
    rec.CountBWAChimericAlignments();
  }

  br.Reset();

  SeqLib::Filter::ReadFilterCollection rf;

}

BOOST_AUTO_TEST_CASE( bam_record_manipulation ) {

  SeqLib::Cigar cig;

  // manually construct a cigar
  cig.add(SeqLib::CigarField('M', 10));
  cig.add(SeqLib::CigarField('I', 1));
  cig.add(SeqLib::CigarField('M', 10));
  cig.add(SeqLib::CigarField('D', 1));
  cig.add(SeqLib::CigarField('M', 10));
  cig.add(SeqLib::CigarField('S', 10));

  // check that coversion to cigar data strutur (uint32_t) worked
  SeqLib::CigarField cm('M', 1);
  SeqLib::CigarField ci('I', 1);
  SeqLib::CigarField cd('D', 1);
  SeqLib::CigarField cn('N', 1);
  SeqLib::CigarField cs('S', 1);
  SeqLib::CigarField ch('H', 1);
  SeqLib::CigarField cp('P', 1);
  SeqLib::CigarField ce('=', 1);
  SeqLib::CigarField cx('X', 1);
  SeqLib::CigarField cb('B', 1);

  BOOST_CHECK_EQUAL(cm.Type(), 'M');
  BOOST_CHECK_EQUAL(ci.Type(), 'I');
  BOOST_CHECK_EQUAL(cd.Type(), 'D');
  BOOST_CHECK_EQUAL(cn.Type(), 'N');
  BOOST_CHECK_EQUAL(cs.Type(), 'S');
  BOOST_CHECK_EQUAL(ch.Type(), 'H');
  BOOST_CHECK_EQUAL(cp.Type(), 'P');
  BOOST_CHECK_EQUAL(ce.Type(), '=');
  BOOST_CHECK_EQUAL(cx.Type(), 'X');

  // check invalid constructions
  BOOST_CHECK_THROW(SeqLib::CigarField('L', 1), std::invalid_argument);

  // make a sequence
  const std::string seq = std::string(10, 'A') + std::string(1, 'T') + std::string(10, 'C') + std::string(10, 'G') + std::string(10, 'A');

  // check
  BOOST_CHECK_EQUAL(cig.NumQueryConsumed(), 41);
  BOOST_CHECK_EQUAL(cig.NumReferenceConsumed(), 31);

  std::stringstream ss;
  ss << cig;

  // cigar from string
  SeqLib::Cigar cig2(ss.str());

  // check that the string from / to are consistent
  assert(cig == cig2);
  assert(!(cig != cig2));
  for (int i = 0; i < cig.size(); ++i)
    assert(cig[i] == cig2[i]);
  for (int i = 0; i < cig.size(); ++i)
    assert(!(cig[i] != cig2[i]));

  // manually make a read
  SeqLib::GenomicRegion gr_wrong(0, 100, 131);
  SeqLib::GenomicRegion gr(0, 100, 130);

  BOOST_CHECK_THROW(SeqLib::BamRecord("dumname", seq, &gr_wrong, cig), std::invalid_argument);
  BOOST_CHECK_THROW(SeqLib::BamRecord("dumname", seq + "A", &gr, cig), std::invalid_argument);

  SeqLib::BamRecord br("dumname", seq, &gr, cig);

  BOOST_CHECK_EQUAL(br.Sequence(), seq);
  BOOST_CHECK_EQUAL(br.GetCigar(), cig);
  BOOST_CHECK_EQUAL(br.Qname(), "dumname");
  BOOST_CHECK_EQUAL(br.Position(), 100);
  BOOST_CHECK_EQUAL(br.Length(), 41);
  BOOST_CHECK_EQUAL(br.ChrID(), 0);

}

BOOST_AUTO_TEST_CASE( change_bam_record ) {

  // get a record
  SeqLib::BamReader br;
  br.Open("test_data/small.bam");
  SeqLib::BamRecord r;

  SeqLib::BamRecordVector brv;

  size_t count = 0;
  br.GetNextRecord(r);

  SeqLib::Cigar c = r.GetCigar();
  std::cerr << c << std::endl;

  // try replace with cigar of same size
  SeqLib::Cigar c2;
  c2.add(SeqLib::CigarField('S', 101));
  r.SetCigar(c2);
  std::cerr << r << std::endl;

  // try replace with new cigar
  SeqLib::Cigar c3;
  c3.add(SeqLib::CigarField('S', 10));
  c3.add(SeqLib::CigarField('M', 91));
  r.SetCigar(c3);
  std::cerr << r << std::endl;

  const std::string new_seq = "ACTGGACTACAC";

  r.SetSequence(new_seq);
  std::cerr << r << std::endl;

  r.SetQname("dummy_qname");
  std::cerr << r << std::endl;

}

BOOST_AUTO_TEST_CASE( stdinput ) {

#ifdef RUN_STDIN
  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("-");

  // write it back out
  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    std::cerr << " STDIN " << r << std::endl;
  }
#endif
}

BOOST_AUTO_TEST_CASE( cramin ) {

  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.cram");

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    std::cerr << "CRAM " << r << std::endl;
  }
}

BOOST_AUTO_TEST_CASE( cramin_new_ref ) {

  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.cram");

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 10) {
    std::cerr << "CRAM " << r << std::endl;
  }

  b.Reset();

  // should fail

}

BOOST_AUTO_TEST_CASE( bamin ) {

  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.bam");

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    std::cerr << "BAM " << r << std::endl;
  }
}

BOOST_AUTO_TEST_CASE( samin ) {

  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.sam");

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    std::cerr << "SAM " << r << std::endl;
  }
}

BOOST_AUTO_TEST_CASE( bamout ) {

  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.sam");

  SeqLib::BamWriter w(SeqLib::BAM);
  //SeqLib::BamWriter w;
  w.Open("tmp_out.bam");
  w.SetHeader(b.Header());
  w.WriteHeader();

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    w.WriteRecord(r);
  }
  w.Close();

}

BOOST_AUTO_TEST_CASE( samout ) {

  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.sam");

  SeqLib::BamWriter w(SeqLib::SAM);
  w.Open("tmp_out.sam");
  w.SetHeader(b.Header());
  w.WriteHeader();

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    w.WriteRecord(r);
  }
  w.Close();
  b.Close();
}

BOOST_AUTO_TEST_CASE( cramout ) {

  SeqLib::BamReader b;
  b.Open("test_data/small.cram");

  SeqLib::BamWriter w(SeqLib::CRAM);
  w.Open("tmp_out.cram");
  w.SetHeader(b.Header());
  w.WriteHeader();

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    w.WriteRecord(r);
  }
  w.Close();

}

BOOST_AUTO_TEST_CASE( samout_to_stdout ) {

#ifdef RUN_SAM_STDOUT
  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.sam");

  SeqLib::BamWriter w(SeqLib::SAM);
  w.Open("-");
  w.SetHeader(b.Header());
  w.WriteHeader();

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    w.WriteRecord(r);
  }
  w.Close();
#endif
}

BOOST_AUTO_TEST_CASE( bamout_to_stdout ) {

  //
  // dont actually run every time
  // too much stdout-ing
  //

#ifdef RUN_BAM_STDOUT
  // read a BAM from stdin
  SeqLib::BamReader b;
  b.Open("test_data/small.sam");

  SeqLib::BamWriter w(SeqLib::BAM);
  w.Open("-");
  w.SetHeader(b.Header());
  w.WriteHeader();

  SeqLib::BamRecord r;
  size_t count = 0;
  while(b.GetNextRecord(r) && count++ < 1) {
    w.WriteRecord(r);
  }
  w.Close();
#endif

}

BOOST_AUTO_TEST_CASE( bam_poly ) {

  SeqLib::BamReader r;

  BOOST_CHECK(r.Open("test_data/small.bam"));
  BOOST_CHECK(r.Open("test_data/small.cram"));

  BOOST_CHECK(r.SetRegion(SeqLib::GenomicRegion(r.Header().Name2ID("X"),1001000, 1001100)));
  BOOST_CHECK(!r.SetRegion(SeqLib::GenomicRegion(1000, 1001000, 1001100))); // should fail

  SeqLib::BamWriter w(SeqLib::BAM);
  w.Open("tmp_out_poly.bam");
  w.SetHeader(r.Header());
  w.WriteHeader();

  SeqLib::BamRecord rec;
  while(r.GetNextRecord(rec)) {
    w.WriteRecord(rec);
  }

  BOOST_CHECK(r.Reset("test_data/small.bam"));
  BOOST_CHECK(!r.Reset("dum"));

  BOOST_CHECK(r.Close("test_data/small.bam"));
  BOOST_CHECK(r.Close("test_data/small.cram"));

  // problematic here FIXME
  //SeqLib::BamReader r2;
  //BOOST_CHECK(r2.Open("tmp_out_poly.bam"));
  // should return false, no index
  //BOOST_CHECK(!r2.SetRegion(SeqLib::GenomicRegion(r.Header().Name2ID("X"),1001000, 1001100)));

}

BOOST_AUTO_TEST_CASE( plot_test ) {

  SeqLib::BamReader r;
  r.Open("test_data/small.bam");

  // should return false on empty region
  BOOST_CHECK(!r.SetMultipleRegions(SeqLib::GRC()));

  SeqLib::GenomicRegion gr("X:1,002,942-1,003,294", r.Header());
  r.SetRegion(gr);

  SeqLib::SeqPlot s;

  s.SetView(gr);

  SeqLib::BamRecord rec;
  SeqLib::BamRecordVector brv;
  while(r.GetNextRecord(rec))
    if (!rec.CountNBases() && rec.MappedFlag())
      brv.push_back(rec);

  s.SetPadding(20);

  std::cout << s.PlotAlignmentRecords(brv);

}

// CURRENTLY DOES NOT WORK
// need to find how to do reset
// BOOST_AUTO_TEST_CASE ( reset_works ) {

//   SeqLib::BamReader r;
//   r.Open("test_data/small.bam");
//   //r.Open("test_data/small.cram");

//   SeqLib::BamRecord rec1, rec2;
//   r.GetNextRecord(rec1);
//   r.Reset();
//   std::cerr << " AFTER RESET " << std::endl;
//   std::cerr << r.GetNextRecord(rec2) << std::endl;

//   BOOST_CHECK_EQUAL(rec1.Qname(), rec2.Qname());
//   }

BOOST_AUTO_TEST_CASE (json_parse) {

  SeqLib::BamReader r;
  r.Open("test_data/small.bam");
  ReadFilterCollection rfc(JSON1, r.Header());

  ReadFilter rf;
  SeqLib::GRC g(VCFFILE, r.Header());
  rf.addRegions(g);
  AbstractRule ar;
  ar.isize = Range(10,100, false);
  rf.SetMateLinked(true);
  rf.AddRule(ar);
  rfc.AddReadFilter(rf);

  std::cout << rfc << std::endl;

  SeqLib::BamRecord rec;
  size_t count = 0;
  int start, end;
  while(r.GetNextRecord(rec) && ++count < 10) {
    rec.QualityTrimmedSequence(4, start, end); // phred trim first
    rfc.isValid(rec);
  }

  // empty
  ReadFilterCollection rfc2("", r.Header());

}

BOOST_AUTO_TEST_CASE ( ref_genome ) {

  //SeqLib::RefGenome r("test_data/test_ref.fa");
  SeqLib::RefGenome r;
  r.LoadIndex("test_data/test_ref.fa");

  BOOST_CHECK(!r.IsEmpty());

  std::string out = r.QueryRegion("ref1", 0, 5);
  BOOST_CHECK_EQUAL(out, "ATCGAC");

  BOOST_CHECK_THROW(r.QueryRegion("ref1", 5,4), std::invalid_argument);
  BOOST_CHECK_THROW(r.QueryRegion("ref1", -1,4), std::invalid_argument);

  SeqLib::RefGenome r2;
  BOOST_CHECK_THROW(r2.QueryRegion("ref1",1,2), std::invalid_argument);

  // reload
  r2.LoadIndex("test_data/test_ref.fa");
}

BOOST_AUTO_TEST_CASE ( set_cigar ) {

  SeqLib::BamReader rr;
  rr.Open(SBAM);
  SeqLib::BamRecord rec;
  size_t count = 0;
  while (rr.GetNextRecord(rec) && ++count < 10) {
      SeqLib::Cigar c;
      c.add(SeqLib::CigarField('M', 70));
      c.add(SeqLib::CigarField('I', 80));
      c.add(SeqLib::CigarField('M',1));
      rec.SetCigar(c);
      std::cerr << rec << std::endl;
  }

}

max_stars_repo_name: agraubert/SeqLib

max_issues_repo_name: GenePlus/ncsv

max_forks_repo_name: GenePlus/ncsv


text:
#pragma once

#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <array>
#include <vector>

namespace ipc::rigid {

class MeshSelector {
public:
    MeshSelector() {}
    MeshSelector(
        size_t num_vertices,
        const Eigen::MatrixXi& E,
        const Eigen::MatrixXi& F);

    size_t vertex_to_edge(size_t vi) const { return m_vertex_to_edge[vi]; }
    size_t vertex_to_face(size_t vi) const { return m_vertex_to_face[vi]; }

    size_t face_to_edge(size_t fi, size_t fi_ei) const
    {
        return m_faces_to_edges(fi, fi_ei);
    }

    const Eigen::MatrixXi& face_to_edges() const { return m_faces_to_edges; }

    size_t edge_to_face(size_t ei) const { return m_edge_to_face[ei]; }

    size_t codim_vertices_to_vertices(size_t vi) const
    {
        return m_codim_vertices_to_vertices[vi];
    }

    size_t codim_edges_to_edges(size_t ei) const
    {
        return m_codim_edges_to_edges[ei];
    }

    const std::vector<size_t>& codim_edges_to_edges() const
    {
        return m_codim_edges_to_edges;
    }

    size_t num_codim_vertices() const
    {
        return m_codim_vertices_to_vertices.size();
    }

    size_t num_codim_edges() const { return m_codim_edges_to_edges.size(); }

protected:
    /// Map from vertices to the minimum index edge containing the vertex.
    void init_vertex_to_edge(size_t num_vertices, const Eigen::MatrixXi& E);

    /// Map from vertices to the minimum index face containing the vertex.
    void init_vertex_to_face(size_t num_vertices, const Eigen::MatrixXi& F);

    /// Each row is a face with the three indices of the edges.
    void init_face_to_edge(
        const Eigen::MatrixXi& F,
        const Eigen::SparseMatrix<size_t>& vertices_to_edge);

    void init_edge_to_face(const size_t num_edges);

    void init_codim_vertices_to_vertices(
        size_t num_vertices, const Eigen::MatrixXi& F);

    void init_codim_edges_to_edges(size_t num_edges);

    /// A mapping from a vertex id to the smallest adjacent face id.
    std::vector<size_t> m_vertex_to_edge;
    /// A mapping from a vertex id to the smallest adjacent face id.
    std::vector<size_t> m_vertex_to_face;
    /// A mapping from a edge id to the smallest adjacent face id.
    std::vector<size_t> m_edge_to_face;
    /// A mapping from a face to the indices of the face's edges.
    Eigen::MatrixXi m_faces_to_edges;
    /// A map from codimensional vertices ids to full vertices ids.
    std::vector<size_t> m_codim_vertices_to_vertices;
    /// A map from codimensional edges ids to full edges ids.
    std::vector<size_t> m_codim_edges_to_edges;
};

} // namespace ipc::rigid

ext: hpp

max_stars_repo_name: ipc-sim/rigid-ipc

max_issues_repo_name: ipc-sim/rigid-ipc

max_forks_repo_name: ipc-sim/rigid-ipc


text:
/*
Copyright (c) 2015, Wator Vapor
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 the name of wator 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.
*/
#pragma once
#include <memory>
#include <vector>
#include <map>
#include <cinttypes>
using namespace std;

#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/json_parser.hpp>
using namespace boost::property_tree;

namespace Wator {
    class LayerBase;
    typedef shared_ptr<LayerBase> LayerBasePtr;
    /**
     * BaseLayer Base of all type Net layers.
     **/
    class LayerBase {
        friend class ImageLayer;
    public:
        /**
         * round
         * @return None.
         **/
        virtual void round(void);

        /**
         * forward
         * @return None.
         **/
        virtual void forward(void);

        /**
         * bottom
         * @param [in] lyr layer.
         * @return None.
         **/
        void bottom(LayerBasePtr lyr);

        /**
         * top
         * @param [in] lyr layer.
         * @return None.
         **/
        void top(LayerBasePtr lyr);

    protected:
        /**
         * Constructor
        **/
        explicit LayerBase();

        /**
         * Constructor
         * @param [in] info layer parameter.
        **/
        explicit LayerBase(const ptree& info);

    protected:
        vector<LayerBasePtr> bottom_;
        vector<LayerBasePtr> top_;
    };

    /**
     * LayerInput Base of Input layers.
     **/
    class LayerInput :public LayerBase{
        friend class CoulombLayer;
        friend class NewtonLayer;
    public:

        /**
         * load
         * @return None.
         **/
        virtual void load(bool train);
    protected:
        /**
         * Constructor
         **/
        explicit LayerInput();

        /**
         * Constructor
         * @param [in] info layer parameter.
         **/
        explicit LayerInput(const ptree& info);

        virtual ~LayerInput();

        /**
         * get ptr
         * @return None.
         **/
        shared_ptr<Blob<uint8_t>> getBlob(const LayerBase *who);
    protected:
        int threshold_;
        vector<shared_ptr<Blob<uint8_t>>> blobs_;
    };

    class ImplicitMemory;
    typedef unique_ptr<ImplicitMemory> ImplicitMemoryPtr;
    /**
     * LayerOutput Base of Hidden layers.
     **/
    class LayerHidden :public LayerBase{
    public:
    protected:
        /**
         * Constructor
         **/
        explicit LayerHidden();

        /**
         * Constructor
         * @param [in] info layer parameter.
         **/
        explicit LayerHidden(const ptree& info);

        /**
         * Desconstructor
         **/
        virtual ~LayerHidden();

        /**
         * get ptr
         * @return None.
         **/
        shared_ptr<Blob<bool>> getBlob(const LayerBase *who);
    protected:
        vector<shared_ptr<Blob<bool>>> pinchs_;
        vector<shared_ptr<Blob<int>>> raws_;
        vector<shared_ptr<Blob<bool>>> blobs_;
        ImplicitMemoryPtr memory_;
    };

    /**
     * LayerOutput Base of Output layers.
     **/
    class LayerOutput :public LayerBase{
    public:
    protected:
        /**
         * Constructor
         **/
        explicit LayerOutput();

        /**
         * Constructor
         * @param [in] info layer parameter.
         **/
        explicit LayerOutput(const ptree& info);

    protected:
    };

}


text:
/*

Copyright (c) 2005-2016, 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 TESTDELTANOTCHODESYSTEM_HPP_
#define TESTDELTANOTCHODESYSTEM_HPP_

#include <cxxtest/TestSuite.h>

#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include <vector>
#include <iostream>

#include "OutputFileHandler.hpp"
#include "DeltaNotchOdeSystem.hpp"
#include "RungeKutta4IvpOdeSolver.hpp"
#include "RungeKuttaFehlbergIvpOdeSolver.hpp"
#include "BackwardEulerIvpOdeSolver.hpp"
#include "CvodeAdaptor.hpp"
#include "Timer.hpp"

//This test is always run sequentially (never in parallel)
#include "FakePetscSetup.hpp"

class TestDeltaNotchOdeSystem : public CxxTest::TestSuite
{
public:

    void TestDeltaNotchOdeSystemSetup() throw(Exception)
    {
#ifdef CHASTE_CVODE
        DeltaNotchOdeSystem ode_system;

        double h_value = 1.0; // Was 0.0001 for some utterly bizarre reason
        CvodeAdaptor cvode_solver;
        OdeSolution solutions;

        std::vector<double> initial_conditions = ode_system.GetInitialConditions();

        TS_ASSERT_DELTA(initial_conditions[0], 1.0, 1e-6);
        TS_ASSERT_DELTA(initial_conditions[1], 1.0, 1e-6);

        Timer::Reset();
        solutions = cvode_solver.Solve(&ode_system, initial_conditions, 0.0, 100.0, h_value, h_value);
        Timer::Print("1. Cvode for 100 hours");

        // Test solutions are OK for a small time increase...
        int end = solutions.rGetSolutions().size() - 1;
        // Tests the simulation is ending at the right time...(going into S phase at 7.8 hours)
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 100, 1e-2);

        // Decent results
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 0.9615, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 0.0107, 1e-4);
#else
        std::cout << "CVODE is not enabled. " << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_CVODE
    }

    void TestArchiving()
    {
#ifdef CHASTE_CVODE
        OutputFileHandler handler("archive", false);
        std::string archive_filename = handler.GetOutputDirectoryFullPath() + "delta_notch_ode.arch";

        {
            std::vector<double> state_variables;
            state_variables.push_back(3.0);
            state_variables.push_back(4.0);

            DeltaNotchOdeSystem ode_system(state_variables);

            ode_system.SetDefaultInitialCondition(1, 3.25);

            std::vector<double> initial_conditions = ode_system.GetInitialConditions();

            // These are the initial conditions hard-coded in the constructor.
            TS_ASSERT_EQUALS(initial_conditions.size(), 2u);
            TS_ASSERT_DELTA(initial_conditions[0], 1.0, 1e-6);
            TS_ASSERT_DELTA(initial_conditions[1], 3.25, 1e-6);

            ode_system.SetParameter("Mean Delta", 10.0);

            // Create an output archive
            std::ofstream ofs(archive_filename.c_str());
            boost::archive::text_oarchive output_arch(ofs);

            // Archive ODE system
            AbstractOdeSystem* const p_const_ode_system = &ode_system;
            output_arch << p_const_ode_system;
        }

        {
            AbstractOdeSystem* p_ode_system;

            // Create an input archive
            std::ifstream ifs(archive_filename.c_str(), std::ios::binary);
            boost::archive::text_iarchive input_arch(ifs);

            // Restore from the archive
            input_arch >> p_ode_system;

            // Check that archiving worked correctly
            std::vector<double> initial_conditions = p_ode_system->GetInitialConditions();

            TS_ASSERT_EQUALS(initial_conditions.size(), 2u);
            TS_ASSERT_DELTA(initial_conditions[0], 1.0, 1e-6);
            TS_ASSERT_DELTA(initial_conditions[1], 1.0, 1e-6);

            double mean_delta = p_ode_system->GetParameter("Mean Delta");

            TS_ASSERT_DELTA(mean_delta, 10.0, 1e-3);

            // Get state variables
            double var1 = p_ode_system->GetStateVariable(0);
            double var2 = p_ode_system->GetStateVariable(1);

            TS_ASSERT_DELTA(var1, 3.0, 1e-3);
            TS_ASSERT_DELTA(var2, 4.0, 1e-3);

            // Tidy up
            delete p_ode_system;
        }
#endif //CHASTE_CVODE
    }

    void TestSetStateVariables()
    {
#ifdef CHASTE_CVODE

        std::vector<double> state_vars;
        state_vars.push_back(0.0);
        state_vars.push_back(1.0);
        DeltaNotchOdeSystem ode_system(state_vars);

        TS_ASSERT_EQUALS(ode_system.GetStateVariable(0),0.0);
        TS_ASSERT_EQUALS(ode_system.GetStateVariable(1),1.0);

#endif //CHASTE_CVODE
   }
};

#endif /*TESTDELTANOTCHODESYSTEM_HPP_*/

max_stars_repo_name: uofs-simlab/ChasteOS

max_issues_repo_name: uofs-simlab/ChasteOS

max_forks_repo_name: uofs-simlab/ChasteOS


text:
#include <iostream>
#include <boost/lexical_cast.hpp>
#include <string>

using namespace std;

int main(){
    int x       = boost::lexical_cast<int>("100");
    long y      = boost::lexical_cast<long>("200000");
    float pai   = boost::lexical_cast<float>("3.1415926");
    double e    = boost::lexical_cast<double>("2.71828");

    double r    = boost::lexical_cast<double>("1.414,x",5);
    cout << r << endl;

    string str = boost::lexical_cast<string>(2323);
    cout << boost::lexical_cast<string>(0.168)  << endl;
    cout << boost::lexical_cast<string>(0x10)   << endl;
    return 0;
}

ext: cpp


text:
#include "db.h"
#include <asio_ipfs.h>
#include "republisher.h"
#include "btree.h"
#include "../or_throw.h"

#include <boost/asio/io_service.hpp>

#include <algorithm>
#include <fstream>

using namespace std;
using namespace ouinet;

static const unsigned int BTREE_NODE_SIZE=64;

static BTree::CatOp make_cat_operation(asio_ipfs::node& ipfs_node)
{
    return [&ipfs_node] (const BTree::Hash& hash, asio::yield_context yield) {
        return ipfs_node.cat(hash, yield);
    };
}

static BTree::AddOp make_add_operation(asio_ipfs::node& ipfs_node)
{
    return [&ipfs_node] (const BTree::Value& value, asio::yield_context yield) {
        sys::error_code ec;

        auto ret = ipfs_node.add(value, yield[ec]);
        if (ec) return or_throw(yield, ec, move(ret));

        ipfs_node.pin(ret, yield[ec]);
        if (ec) return or_throw(yield, ec, move(ret));

        return ret;
    };
}

static BTree::RemoveOp make_remove_operation(asio_ipfs::node& ipfs_node)
{
    return [&ipfs_node] (const BTree::Value& hash, asio::yield_context yield) {
        ipfs_node.unpin(hash, yield);
    };
}

static string path_to_db(const string& path_to_repo, const string& ipns)
{
    return path_to_repo + "/ipfs_cache_db." + ipns;
}

static void load_db( BTree& db_map
                   , const string& path_to_repo
                   , const string& ipns
                   , asio::yield_context yield)
{
    string path = path_to_db(path_to_repo, ipns);

    ifstream file(path);

    if (!file.is_open()) {
        cerr << "Warning: Couldn't open " << path << endl;
        return;
    }

    try {
        string ipfs;
        file >> ipfs;

        if (ipfs.substr(0, 2) != "Qm") {
            throw runtime_error("Content doesn't start with 'Qm'");
        }

        if (ipfs.size() != asio_ipfs::node::CID_SIZE) {
            throw runtime_error("Content doesn't appear to be a CID hash");
        }

        sys::error_code ec;
        db_map.load(ipfs, yield[ec]);
    }
    catch (const std::exception& e) {
        cerr << "ERROR: parsing " << path << ": " << e.what() << endl;
    }
}

static void save_db( const string& path_to_repo
                   , const string& ipns
                   , const string& ipfs)
{
    string path = path_to_db(path_to_repo, ipns);

    ofstream file(path, std::ofstream::trunc);

    if (!file.is_open()) {
        cerr << "ERROR: Saving " << path << endl;
        return;
    }

    file << ipfs;
    file.close();
}

ClientDb::ClientDb(asio_ipfs::node& ipfs_node, string path_to_repo, string ipns)
    : _path_to_repo(move(path_to_repo))
    , _ipns(move(ipns))
    , _ipfs_node(ipfs_node)
    , _was_destroyed(make_shared<bool>(false))
    , _download_timer(_ipfs_node.get_io_service())
    , _db_map(make_unique<BTree>( make_cat_operation(ipfs_node)
                                , nullptr
                                , nullptr
                                , BTREE_NODE_SIZE))
{
    auto d = _was_destroyed;

    asio::spawn(get_io_service(), [=](asio::yield_context yield) {
            if (*d) return;
            load_db(*_db_map, _path_to_repo, _ipns, yield);
            if (*d) return;
            continuously_download_db(yield);
        });
}

InjectorDb::InjectorDb(asio_ipfs::node& ipfs_node, string path_to_repo)
    : _path_to_repo(move(path_to_repo))
    , _ipns(ipfs_node.id())
    , _ipfs_node(ipfs_node)
    , _republisher(new Republisher(_ipfs_node))
    , _has_callbacks(_ipfs_node.get_io_service())
    , _was_destroyed(make_shared<bool>(false))
    , _db_map(make_unique<BTree>( make_cat_operation(ipfs_node)
                                , make_add_operation(ipfs_node)
                                , make_remove_operation(ipfs_node)
                                , BTREE_NODE_SIZE))
{
    auto d = _was_destroyed;

    asio::spawn(get_io_service(), [=](asio::yield_context yield) {
            if (*d) return;
            load_db(*_db_map, _path_to_repo, _ipns, yield);
        });
}

const string ipfs_uri_prefix = "ipfs:/ipfs/";

void InjectorDb::update(string key, string value, asio::yield_context yield)
{
    auto wd = _was_destroyed;
    sys::error_code ec;

    _db_map->insert(move(key), move(value), yield[ec]);

    if (!ec && *wd) ec = asio::error::operation_aborted;
    if (ec) return or_throw(yield, ec);

    upload_database(yield[ec]);

    if (!ec && *wd) ec = asio::error::operation_aborted;
    return or_throw(yield, ec);
}

void InjectorDb::upload_database(asio::yield_context yield)
{
    string db_ipfs_id = _db_map->root_hash();

    if (db_ipfs_id.empty()) {
        return;
    }

    save_db(_path_to_repo, _ipns, db_ipfs_id);

    _republisher->publish(move(db_ipfs_id), yield);
}

static string query_(string key, BTree& db, asio::yield_context yield)
{
    sys::error_code ec;

    auto val = db.find(key, yield[ec]);

    if (ec) return or_throw<string>(yield, ec);

    return val;
}

string InjectorDb::query(string key, asio::yield_context yield)
{
    return query_(move(key), *_db_map, yield);
}

string ClientDb::query(string key, asio::yield_context yield)
{
    return query_(move(key), *_db_map, yield);
}

void ClientDb::continuously_download_db(asio::yield_context yield)
{
    auto d = _was_destroyed;

    while(true) {
        sys::error_code ec;

        auto ipfs_id = _ipfs_node.resolve(_ipns, yield[ec]);
        if (*d) return;

        if (!ec) {
            _ipfs = ipfs_id;

            _db_map->load(ipfs_id, yield[ec]);

            if (*d) return;
        }

        save_db(_path_to_repo, _ipns, ipfs_id);

        if (ec) {
            _download_timer.expires_from_now(chrono::seconds(5));
            _download_timer.async_wait(yield[ec]);
            if (*d) return;
            continue;
        }

        flush_db_update_callbacks(sys::error_code());

        _download_timer.expires_from_now(chrono::seconds(5));
        _download_timer.async_wait(yield[ec]);
        if (*d) return;
    }
}

void ClientDb::wait_for_db_update(asio::yield_context yield)
{
    using Handler = asio::handler_type<asio::yield_context,
          void(sys::error_code)>::type;

    Handler h(yield);
    asio::async_result<Handler> result(h);
    _on_db_update_callbacks.push([ h = move(h)
                                 , w = asio::io_service::work(get_io_service())
                                 ] (auto ec) mutable { h(ec); });
    result.get();
}

void ClientDb::flush_db_update_callbacks(const sys::error_code& ec)
{
    auto& q = _on_db_update_callbacks;

    while (!q.empty()) {
        auto c = move(q.front());
        q.pop();
        get_io_service().post([c = move(c), ec] () mutable { c(ec); });
    }
}

asio::io_service& ClientDb::get_io_service() {
    return _ipfs_node.get_io_service();
}

asio::io_service& InjectorDb::get_io_service() {
    return _ipfs_node.get_io_service();
}

ClientDb::~ClientDb() {
    *_was_destroyed = true;
    flush_db_update_callbacks(asio::error::operation_aborted);
}

InjectorDb::~InjectorDb() {
    *_was_destroyed = true;

    for (auto& cb : _upload_callbacks) {
        get_io_service().post([cb = move(cb)] {
                cb(asio::error::operation_aborted);
            });
    }
}


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

    @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_FUNCTION_IF_ELSE_NAN_HPP_INCLUDED
#define BOOST_SIMD_FUNCTION_IF_ELSE_NAN_HPP_INCLUDED

#if defined(DOXYGEN_ONLY)
namespace boost { namespace simd
{

  /*!
    @ingroup group-boolean
    This function object conditionally returns @c x (respectively @ref Nan)
    if @c c is @ref True (respectively  @ref False)

    This is an alias of @ref if_else_allbits

    @par Header <boost/simd/function/if_else_nan.hpp>

    @par Example:

      @snippet if_else_nan.cpp if_else_nan

    @par Possible output:

      @snippet if_else_nan.txt if_else_nan

  **/
  Value1 if_else_nan(Value0 const& c, Value1 const& x);
} }
#endif

#include <boost/simd/function/scalar/if_else_allbits.hpp>
#include <boost/simd/function/simd/if_else_nan.hpp>

#endif

max_stars_repo_name: SylvainCorlay/pythran

max_issues_repo_name: SylvainCorlay/pythran

max_forks_repo_name: SylvainCorlay/pythran


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

  @copyright 2015 NumScale SAS
  @copyright 2015 J.T. Lapreste

  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_ARCH_COMMON_SCALAR_FUNCTION_COPYSIGN_HPP_INCLUDED
#define BOOST_SIMD_ARCH_COMMON_SCALAR_FUNCTION_COPYSIGN_HPP_INCLUDED

#include <boost/simd/function/scalar/abss.hpp>
#include <boost/simd/function/scalar/bitwise_notand.hpp>
#include <boost/simd/function/scalar/bitofsign.hpp>
#include <boost/simd/function/scalar/signnz.hpp>
#include <boost/simd/constant/signmask.hpp>
#include <boost/simd/detail/math.hpp>
#include <boost/dispatch/function/overload.hpp>
#include <boost/config.hpp>

namespace boost { namespace simd { namespace ext
{
  namespace bd = boost::dispatch;
  BOOST_DISPATCH_OVERLOAD ( copysign_
                          , (typename A0)
                          , bd::cpu_
                          , bd::scalar_< bd::arithmetic_<A0> >
                          , bd::scalar_< bd::arithmetic_<A0> >
                          )
  {
    BOOST_FORCEINLINE A0 operator() ( A0 a0, A0 a1) const BOOST_NOEXCEPT
    {
      return abss(a0)*signnz(a1);
    }
  };

  BOOST_DISPATCH_OVERLOAD ( copysign_
                          , (typename A0)
                          , bd::cpu_
                          , bd::scalar_< bd::double_<A0> >
                          , bd::scalar_< bd::double_<A0> >
                          )
  {
    BOOST_FORCEINLINE A0 operator() ( A0 a0, A0 a1) const BOOST_NOEXCEPT
    {
    #ifdef BOOST_SIMD_HAS_COPYSIGN
      return ::copysign(a0, a1);
      // _copysign appears to be bogus for a1 = -0 in old MSVCRT
    #elif defined(BOOST_SIMD_HAS__COPYSIGN) && !defined(__MSVCRT__)
      return ::_copysign(a0, a1);
    #else
      return bitwise_or(bitofsign(a1), bitwise_notand(Signmask<A0>(), a0));
    #endif
    }
  };

  BOOST_DISPATCH_OVERLOAD ( copysign_
                          , (typename A0)
                          , bd::cpu_
                          , bd::scalar_< bd::single_<A0> >
                          , bd::scalar_< bd::single_<A0> >
                                    )
  {
    BOOST_FORCEINLINE A0 operator() ( A0 a0, A0 a1) const BOOST_NOEXCEPT
    {
    #ifdef BOOST_SIMD_HAS_COPYSIGNF
      return ::copysignf(a0, a1);
    #elif defined(BOOST_SIMD_HAS__COPYSIGNF)
      return ::_copysignf(a0, a1);
    #else
       return bitwise_or(bitofsign(a1), bitwise_notand(Signmask<A0>(), a0));
    #endif
    }
  };
} } }

#endif


text:
///////////////////////////////////////////////////////////////
//  Copyright 2012 John Maddock. 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_
//
// Comparison operators for cpp_int_backend:
//
#ifndef BOOST_MP_CPP_INT_MISC_HPP
#define BOOST_MP_CPP_INT_MISC_HPP

#include <boost/multiprecision/detail/bitscan.hpp> // lsb etc

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4702)
#endif

namespace boost{ namespace multiprecision{ namespace backends{

template <class R, class CppInt>
void check_in_range(const CppInt& val, const mpl::int_<checked>&)
{
   typedef typename boost::multiprecision::detail::canonical<R, CppInt>::type cast_type;
   if(val.sign())
   {
      if(val.compare(static_cast<cast_type>((std::numeric_limits<R>::min)())) < 0)
         BOOST_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
   }
   else
   {
      if(val.compare(static_cast<cast_type>((std::numeric_limits<R>::max)())) > 0)
         BOOST_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
   }
}
template <class R, class CppInt>
inline void check_in_range(const CppInt& /*val*/, const mpl::int_<unchecked>&) BOOST_NOEXCEPT {}

inline void check_is_negative(const mpl::true_&) BOOST_NOEXCEPT {}
inline void check_is_negative(const mpl::false_&)
{
   BOOST_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
}

template <class Integer>
inline Integer negate_integer(Integer i, const mpl::true_&) BOOST_NOEXCEPT
{
   return -i;
}
template <class Integer>
inline Integer negate_integer(Integer i, const mpl::false_&) BOOST_NOEXCEPT
{
   return ~(i-1);
}

template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<is_integral<R>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, void>::type
   eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) BOOST_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
{
   typedef mpl::int_<Checked1> checked_type;
   check_in_range<R>(backend, checked_type());

   *result = static_cast<R>(backend.limbs()[0]);
   unsigned shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   for(unsigned i = 1; (i < backend.size()) && (shift < static_cast<unsigned>(std::numeric_limits<R>::digits)); ++i)
   {
      *result += static_cast<R>(backend.limbs()[i]) << shift;
      shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   }
   if(backend.sign())
   {
      check_is_negative(boost::is_signed<R>());
      *result = negate_integer(*result, boost::is_signed<R>());
   }
}

template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<is_floating_point<R>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, void>::type
   eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) BOOST_NOEXCEPT_IF(is_arithmetic<R>::value)
{
   typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs();
   unsigned shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   *result = static_cast<R>(*p);
   for(unsigned i = 1; i < backend.size(); ++i)
   {
      *result += static_cast<R>(std::ldexp(static_cast<long double>(p[i]), shift));
      shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   }
   if(backend.sign())
      *result = -*result;
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type
   eval_is_zero(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) BOOST_NOEXCEPT
{
   return (val.size() == 1) && (val.limbs()[0] == 0);
}
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, int>::type
   eval_get_sign(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) BOOST_NOEXCEPT
{
   return eval_is_zero(val) ? 0 : val.sign() ? -1 : 1;
}
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_abs(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) BOOST_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
{
   result = val;
   result.sign(false);
}

//
// Get the location of the least-significant-bit:
//
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
   eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
{
   using default_ops::eval_get_sign;
   if(eval_get_sign(a) == 0)
   {
      BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
   }
   if(a.sign())
   {
      BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
   }

   //
   // Find the index of the least significant limb that is non-zero:
   //
   unsigned index = 0;
   while(!a.limbs()[index] && (index < a.size()))
      ++index;
   //
   // Find the index of the least significant bit within that limb:
   //
   unsigned result = boost::multiprecision::detail::find_lsb(a.limbs()[index]);

   return result + index * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
}

//
// Get the location of the most-significant-bit:
//
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
   eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
{
   using default_ops::eval_get_sign;
   if(eval_get_sign(a) == 0)
   {
      BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
   }
   if(a.sign())
   {
      BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
   }

   //
   // Find the index of the most significant bit that is non-zero:
   //
   return (a.size() - 1) * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits + boost::multiprecision::detail::find_msb(a.limbs()[a.size() - 1]);
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type
   eval_bit_test(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index) BOOST_NOEXCEPT
{
   unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
   if(offset >= val.size())
      return false;
   return val.limbs()[offset] & mask ? true : false;
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_bit_set(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index)
{
   unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
   if(offset >= val.size())
   {
      unsigned os = val.size();
      val.resize(offset + 1, offset + 1);
      if(offset >= val.size())
         return;  // fixed precision overflow
      for(unsigned i = os; i <= offset; ++i)
         val.limbs()[i] = 0;
   }
   val.limbs()[offset] |= mask;
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_bit_unset(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index) BOOST_NOEXCEPT
{
   unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
   if(offset >= val.size())
      return;
   val.limbs()[offset] &= ~mask;
   val.normalize();
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_bit_flip(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index)
{
   unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
   limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
   if(offset >= val.size())
   {
      unsigned os = val.size();
      val.resize(offset + 1, offset + 1);
      if(offset >= val.size())
         return;  // fixed precision overflow
      for(unsigned i = os; i <= offset; ++i)
         val.limbs()[i] = 0;
   }
   val.limbs()[offset] ^= mask;
   val.normalize();
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_qr(
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& y,
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) BOOST_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
{
   divide_unsigned_helper(&q, x, y, r);
   q.sign(x.sign() != y.sign());
   r.sign(x.sign());
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_qr(
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
      limb_type y,
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) BOOST_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
{
   divide_unsigned_helper(&q, x, y, r);
   q.sign(x.sign());
   r.sign(x.sign());
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class U>
inline typename enable_if_c<is_integral<U>::value>::type eval_qr(
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
      U y,
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) BOOST_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
{
   using default_ops::eval_qr;
   cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(y);
   eval_qr(x, t, q, r);
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
inline typename enable_if_c<is_unsigned<Integer>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, Integer>::type
   eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, Integer val)
{
   if((sizeof(Integer) <= sizeof(limb_type)) || (val <= (std::numeric_limits<limb_type>::max)()))
   {
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> d;
      divide_unsigned_helper(static_cast<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>*>(0), x, static_cast<limb_type>(val), d);
      return d.limbs()[0];
   }
   else
   {
      return default_ops::eval_integer_modulus(x, val);
   }
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
BOOST_MP_FORCEINLINE typename enable_if_c<is_signed<Integer>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, Integer>::type
   eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, Integer val)
{
   return eval_integer_modulus(x, boost::multiprecision::detail::unsigned_abs(val));
}

inline limb_type integer_gcd_reduce(limb_type u, limb_type v)
{
   do
   {
      if(u > v)
         std::swap(u, v);
      if(u == v)
         break;
      v -= u;
      v >>= boost::multiprecision::detail::find_lsb(v);
   } while(true);
   return u;
}

inline double_limb_type integer_gcd_reduce(double_limb_type u, double_limb_type v)
{
   do
   {
      if(u > v)
         std::swap(u, v);
      if(u == v)
         break;
      if(v <= ~static_cast<limb_type>(0))
      {
         u = integer_gcd_reduce(static_cast<limb_type>(v), static_cast<limb_type>(u));
         break;
      }
      v -= u;
      while((static_cast<unsigned>(v) & 1u) == 0)
         v >>= 1;
   } while(true);
   return u;
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_gcd(
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
      limb_type v)
{
   using default_ops::eval_lsb;
   using default_ops::eval_is_zero;
   using default_ops::eval_get_sign;

   int shift;

   cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> u(a);

   int s = eval_get_sign(u);

   /* GCD(0,x) := x */
   if(s < 0)
   {
      u.negate();
   }
   else if(s == 0)
   {
      result = v;
      return;
   }
   if(v == 0)
   {
      result = u;
      return;
   }

   /* Let shift := lg K, where K is the greatest power of 2
   dividing both u and v. */

   unsigned us = eval_lsb(u);
   unsigned vs = boost::multiprecision::detail::find_lsb(v);
   shift = (std::min)(us, vs);
   eval_right_shift(u, us);
   if(vs)
      v >>= vs;

   do
   {
      /* Now u and v are both odd, so diff(u, v) is even.
      Let u = min(u, v), v = diff(u, v)/2. */
      if(u.size() <= 2)
      {
         if(u.size() == 1)
            v = integer_gcd_reduce(*u.limbs(), v);
         else
         {
            double_limb_type i;
            i = u.limbs()[0] | (static_cast<double_limb_type>(u.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
            v = static_cast<limb_type>(integer_gcd_reduce(i, static_cast<double_limb_type>(v)));
         }
         break;
      }
      eval_subtract(u, v);
      us = eval_lsb(u);
      eval_right_shift(u, us);
   }
   while(true);

   result = v;
   eval_left_shift(result, shift);
}
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
inline typename enable_if_c<is_unsigned<Integer>::value && (sizeof(Integer) <= sizeof(limb_type)) && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_gcd(
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
      const Integer& v)
{
   eval_gcd(result, a, static_cast<limb_type>(v));
}
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
inline typename enable_if_c<is_signed<Integer>::value && (sizeof(Integer) <= sizeof(limb_type)) && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_gcd(
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
      const Integer& v)
{
   eval_gcd(result, a, static_cast<limb_type>(v < 0 ? -v : v));
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_gcd(
      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
{
   using default_ops::eval_lsb;
   using default_ops::eval_is_zero;
   using default_ops::eval_get_sign;

   if(a.size() == 1)
   {
      eval_gcd(result, b, *a.limbs());
      return;
   }
   if(b.size() == 1)
   {
      eval_gcd(result, a, *b.limbs());
      return;
   }

   int shift;

   cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> u(a), v(b);

   int s = eval_get_sign(u);

   /* GCD(0,x) := x */
   if(s < 0)
   {
      u.negate();
   }
   else if(s == 0)
   {
      result = v;
      return;
   }
   s = eval_get_sign(v);
   if(s < 0)
   {
      v.negate();
   }
   else if(s == 0)
   {
      result = u;
      return;
   }

   /* Let shift := lg K, where K is the greatest power of 2
   dividing both u and v. */

   unsigned us = eval_lsb(u);
   unsigned vs = eval_lsb(v);
   shift = (std::min)(us, vs);
   eval_right_shift(u, us);
   eval_right_shift(v, vs);

   do
   {
      /* Now u and v are both odd, so diff(u, v) is even.
      Let u = min(u, v), v = diff(u, v)/2. */
      s = u.compare(v);
      if(s > 0)
         u.swap(v);
      if(s == 0)
         break;
      if(v.size() <= 2)
      {
         if(v.size() == 1)
            u = integer_gcd_reduce(*v.limbs(), *u.limbs());
         else
         {
            double_limb_type i, j;
            i = v.limbs()[0] | (static_cast<double_limb_type>(v.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
            j = (u.size() == 1) ? *u.limbs() : u.limbs()[0] | (static_cast<double_limb_type>(u.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
            u = integer_gcd_reduce(i, j);
         }
         break;
      }
      eval_subtract(v, u);
      vs = eval_lsb(v);
      eval_right_shift(v, vs);
   }
   while(true);

   result = u;
   eval_left_shift(result, shift);
}
//
// Now again for trivial backends:
//
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
BOOST_MP_FORCEINLINE typename enable_if_c<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
   eval_gcd(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) BOOST_NOEXCEPT
{
   *result.limbs() = boost::math::gcd(*a.limbs(), *b.limbs());
}
// This one is only enabled for unchecked cpp_int's, for checked int's we need the checking in the default version:
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
BOOST_MP_FORCEINLINE typename enable_if_c<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && (Checked1 == unchecked)>::type
   eval_lcm(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) BOOST_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
{
   *result.limbs() = boost::math::lcm(*a.limbs(), *b.limbs());
   result.normalize(); // result may overflow the specified number of bits
}

inline void conversion_overflow(const mpl::int_<checked>&)
{
   BOOST_THROW_EXCEPTION(std::overflow_error("Overflow in conversion to narrower type"));
}
inline void conversion_overflow(const mpl::int_<unchecked>&){}

template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<
            is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
            && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
         >::type
   eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
{
   typedef typename common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type common_type;
   if(std::numeric_limits<R>::is_specialized && (static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)())))
   {
      conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
      *result = (std::numeric_limits<R>::max)();
   }
   else
      *result = static_cast<R>(*val.limbs());
   if(val.isneg())
   {
      check_is_negative(mpl::bool_<boost::is_signed<R>::value || boost::is_floating_point<R>::value>());
      *result = negate_integer(*result, mpl::bool_<boost::is_signed<R>::value || boost::is_floating_point<R>::value>());
   }
}

template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<
            is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
            && is_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
         >::type
   eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
{
   typedef typename common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type common_type;
   if(std::numeric_limits<R>::is_specialized && (static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)())))
   {
      conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
      *result = (std::numeric_limits<R>::max)();
   }
   else
      *result = static_cast<R>(*val.limbs());
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
   eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
{
   using default_ops::eval_get_sign;
   if(eval_get_sign(a) == 0)
   {
      BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
   }
   if(a.sign())
   {
      BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
   }
   //
   // Find the index of the least significant bit within that limb:
   //
   return boost::multiprecision::detail::find_lsb(*a.limbs());
}

template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
inline typename enable_if_c<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
   eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
{
   using default_ops::eval_get_sign;
   if(eval_get_sign(a) == 0)
   {
      BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
   }
   if(a.sign())
   {
      BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
   }
   //
   // Find the index of the least significant bit within that limb:
   //
   return boost::multiprecision::detail::find_msb(*a.limbs());
}

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

}}} // namespaces

#endif

max_stars_repo_name: bowlofstew/omim

max_issues_repo_name: bowlofstew/omim

max_forks_repo_name: bowlofstew/omim