Browse code

v. 2.99.6: no v.1

ramon diaz-uriarte (at Phelsuma) authored on 30/12/2020 00:29:19
Showing26 changed files

... ...
@@ -1,8 +1,8 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 
4
-Version: 2.99.5
5
-Date: 2020-12-18
4
+Version: 2.99.6
5
+Date: 2020-12-13
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -499,6 +499,8 @@ oncoSimulIndiv <- function(fp,
499 499
     if(inherits(fp, "fitnessEffects")) {
500 500
         s <- sh <- NULL ## force it soon!
501 501
     }
502
+    if(!inherits(fp, "fitnessEffects")) 
503
+        stop("v.1 functionality has been removed. Please use v.2")
502 504
 
503 505
     ## legacies from poor name choices
504 506
     typeFitness <- switch(model,
... ...
@@ -549,8 +551,8 @@ oncoSimulIndiv <- function(fp,
549 551
             stop("Unknown model")
550 552
     }
551 553
 
552
-    if( (length(mu) > 1) && !inherits(fp, "fitnessEffects"))
553
-        stop("Per-gene mutation rates cannot be used with the old poset format")
554
+    ## if( (length(mu) > 1) && !inherits(fp, "fitnessEffects"))
555
+    ##     stop("Per-gene mutation rates cannot be used with the old poset format")
554 556
 
555 557
     if(any(mu < 0)) {
556 558
         stop("(at least one) mutation rate (mu) is negative")
... ...
@@ -584,83 +586,83 @@ oncoSimulIndiv <- function(fp,
584 586
 
585 587
     if(is_null_na(sampleEvery)) stop("sampleEvery cannot be NULL or NA")
586 588
     
587
-    if(!inherits(fp, "fitnessEffects")) {
588
-        if(any(unlist(lapply(list(fp, 
589
-                                  numPassengers,
590
-                                  s, sh), is.null)))) {
591
-            m <- paste("You are using the old poset format.",
592
-                       "You must specify all of poset, numPassengers",
593
-                       "s, and sh.")
594
-            stop(m)
589
+    ## if(!inherits(fp, "fitnessEffects")) {
590
+        ## if(any(unlist(lapply(list(fp, 
591
+        ##                           numPassengers,
592
+        ##                           s, sh), is.null)))) {
593
+        ##     m <- paste("You are using the old poset format.",
594
+        ##                "You must specify all of poset, numPassengers",
595
+        ##                "s, and sh.")
596
+        ##     stop(m)
595 597
            
596
-        }
597
-        if(AND_DrvProbExit) {
598
-            stop("The AND_DrvProbExit = TRUE setting is invalid",
599
-                 " with the old poset format.")
600
-        }
601
-        if(!is.null(muEF))
602
-            stop("Mutator effects cannot be specified with the old poset format.")
603
-        if( length(initMutant) > 0)  
604
-            warning("With the old poset format you can no longer use initMutant.",
605
-                    " The initMutant you passed will be ignored.")
606
-        ## stop("With the old poset, initMutant can only take a single value.")
607
-        if(!is_null_na(fixation))
608
-            stop("'fixation' cannot be specified with the old poset format.")
609
-        ## Seeding C++ is now much better in new version
610
-        if(is.null(seed) || (seed == "auto")) {## passing a null creates a random seed
611
-            ## name is a legacy. This is really the seed for the C++ generator.
612
-            ## Nope, we cannot use 2^32, because as.integer will fail.
613
-            seed <- as.integer(round(runif(1, min = 0, max = 2^16)))
614
-        }
615
-        if(verbosity >= 2)
616
-            cat(paste("\n Using ", seed, " as seed for C++ generator\n"))
617
-
618
-        if(!is_null_na(detectionProb)) stop("detectionProb cannot be used in v.1 objects")
619
-        ## if(message.v1)
620
-        ##     message("You are using the old poset format. Consider using the new one.")
598
+        ## }
599
+        ## if(AND_DrvProbExit) {
600
+        ##     stop("The AND_DrvProbExit = TRUE setting is invalid",
601
+        ##          " with the old poset format.")
602
+        ## }
603
+        ## if(!is.null(muEF))
604
+        ##     stop("Mutator effects cannot be specified with the old poset format.")
605
+        ## if( length(initMutant) > 0)  
606
+        ##     warning("With the old poset format you can no longer use initMutant.",
607
+        ##             " The initMutant you passed will be ignored.")
608
+        ## ## stop("With the old poset, initMutant can only take a single value.")
609
+        ## if(!is_null_na(fixation))
610
+        ##     stop("'fixation' cannot be specified with the old poset format.")
611
+        ## ## Seeding C++ is now much better in new version
612
+        ## if(is.null(seed) || (seed == "auto")) {## passing a null creates a random seed
613
+        ##     ## name is a legacy. This is really the seed for the C++ generator.
614
+        ##     ## Nope, we cannot use 2^32, because as.integer will fail.
615
+        ##     seed <- as.integer(round(runif(1, min = 0, max = 2^16)))
616
+        ## }
617
+        ## if(verbosity >= 2)
618
+        ##     cat(paste("\n Using ", seed, " as seed for C++ generator\n"))
619
+
620
+        ## if(!is_null_na(detectionProb)) stop("detectionProb cannot be used in v.1 objects")
621
+        ## ## if(message.v1)
622
+        ## ##     message("You are using the old poset format. Consider using the new one.")
621 623
    
622 624
     
623
-        ## A simulation stops if cancer or finalTime appear, the first
624
-        ## one. But if we set onlyCnacer = FALSE, we also accept simuls
625
-        ## without cancer (or without anything)
625
+        ## ## A simulation stops if cancer or finalTime appear, the first
626
+        ## ## one. But if we set onlyCnacer = FALSE, we also accept simuls
627
+        ## ## without cancer (or without anything)
626 628
         
627
-        op <- try(oncoSimul.internal(poset = fp, ## restrict.table = rt,
628
-                                     ## numGenes = numGenes,
629
-                                     numPassengers = numPassengers,
630
-                                     typeCBN = "CBN",
631
-                                     birth = birth,
632
-                                     s = s,
633
-                                     death = death,  
634
-                                     mu =  mu,  
635
-                                     initSize =  initSize, 
636
-                                     sampleEvery =  sampleEvery,  
637
-                                     detectionSize =  detectionSize, 
638
-                                     finalTime = finalTime, 
639
-                                     initSize_species = 2000, 
640
-                                     initSize_iter = 500, 
641
-                                     seed = seed, 
642
-                                     verbosity = verbosity, 
643
-                                     speciesFS = 10000,  
644
-                                     ratioForce = 2,
645
-                                     typeFitness = typeFitness,
646
-                                     max.memory = max.memory,
647
-                                     mutationPropGrowth = mutationPropGrowth,                                   
648
-                                     initMutant = -1, 
649
-                                     max.wall.time = max.wall.time,
650
-                                     max.num.tries = max.num.tries,
651
-                                     keepEvery = keepEvery,  
652
-                                     ## alpha = 0.0015,  
653
-                                     sh = sh,
654
-                                     K = K, 
655
-                                     minDetectDrvCloneSz = minDetectDrvCloneSz,
656
-                                     extraTime = extraTime,
657
-                                     detectionDrivers = detectionDrivers,
658
-                                     onlyCancer = onlyCancer,
659
-                                     errorHitWallTime = errorHitWallTime,
660
-                                     errorHitMaxTries = errorHitMaxTries),
661
-                  silent = !verbosity)
662
-        objClass <- "oncosimul"
663
-    } else {
629
+        ## op <- try(oncoSimul.internal(poset = fp, ## restrict.table = rt,
630
+        ##                              ## numGenes = numGenes,
631
+        ##                              numPassengers = numPassengers,
632
+        ##                              typeCBN = "CBN",
633
+        ##                              birth = birth,
634
+        ##                              s = s,
635
+        ##                              death = death,  
636
+        ##                              mu =  mu,  
637
+        ##                              initSize =  initSize, 
638
+        ##                              sampleEvery =  sampleEvery,  
639
+        ##                              detectionSize =  detectionSize, 
640
+        ##                              finalTime = finalTime, 
641
+        ##                              initSize_species = 2000, 
642
+        ##                              initSize_iter = 500, 
643
+        ##                              seed = seed, 
644
+        ##                              verbosity = verbosity, 
645
+        ##                              speciesFS = 10000,  
646
+        ##                              ratioForce = 2,
647
+        ##                              typeFitness = typeFitness,
648
+        ##                              max.memory = max.memory,
649
+        ##                              mutationPropGrowth = mutationPropGrowth,                                   
650
+        ##                              initMutant = -1, 
651
+        ##                              max.wall.time = max.wall.time,
652
+        ##                              max.num.tries = max.num.tries,
653
+        ##                              keepEvery = keepEvery,  
654
+        ##                              ## alpha = 0.0015,  
655
+        ##                              sh = sh,
656
+        ##                              K = K, 
657
+        ##                              minDetectDrvCloneSz = minDetectDrvCloneSz,
658
+        ##                              extraTime = extraTime,
659
+        ##                              detectionDrivers = detectionDrivers,
660
+        ##                              onlyCancer = onlyCancer,
661
+        ##                              errorHitWallTime = errorHitWallTime,
662
+        ##                              errorHitMaxTries = errorHitMaxTries),
663
+        ##           silent = !verbosity)
664
+        ## objClass <- "oncosimul"
665
+    ## } else {
664 666
         s <- sh <- NULL ## force it.
665 667
         if(numPassengers != 0)
666 668
             warning(paste("Specifying numPassengers has no effect",
... ...
@@ -728,7 +730,7 @@ oncoSimulIndiv <- function(fp,
728 730
                                         fixation = fixation),
729 731
                   silent = !verbosity)
730 732
         objClass <- c("oncosimul", "oncosimul2")
731
-    }
733
+    ## }
732 734
     if(inherits(op, "try-error")) {
733 735
         ##         if(length(grep("BAIL OUT NOW", op)))
734 736
         stop(paste("Unrecoverable error:", op ))
... ...
@@ -1763,156 +1765,156 @@ get.mut.vector <- function(x, timeSample, typeSample,
1763 1765
 
1764 1766
 
1765 1767
 
1766
-oncoSimul.internal <- function(poset, ## restrict.table,
1767
-                               numPassengers, 
1768
-                               ## numGenes,
1769
-                               typeCBN,
1770
-                               birth, 
1771
-                               s,
1772
-                               death,
1773
-                               mu,
1774
-                               initSize,
1775
-                               sampleEvery,
1776
-                               detectionSize,
1777
-                               finalTime,
1778
-                               initSize_species,
1779
-                               initSize_iter,
1780
-                               seed,
1781
-                               verbosity,
1782
-                               speciesFS,
1783
-                               ratioForce,
1784
-                               typeFitness,
1785
-                               max.memory,
1786
-                               mutationPropGrowth, ## make it explicit
1787
-                               initMutant,
1788
-                               max.wall.time,
1789
-                               keepEvery,
1790
-                               alpha,
1791
-                               sh,                               
1792
-                               K,
1793
-                               ## endTimeEvery,
1794
-                               detectionDrivers,
1795
-                               onlyCancer,
1796
-                               errorHitWallTime,
1797
-                               max.num.tries,
1798
-                               errorHitMaxTries,
1799
-                               minDetectDrvCloneSz,
1800
-                               extraTime) {
1768
+## oncoSimul.internal <- function(poset, ## restrict.table,
1769
+##                                numPassengers, 
1770
+##                                ## numGenes,
1771
+##                                typeCBN,
1772
+##                                birth, 
1773
+##                                s,
1774
+##                                death,
1775
+##                                mu,
1776
+##                                initSize,
1777
+##                                sampleEvery,
1778
+##                                detectionSize,
1779
+##                                finalTime,
1780
+##                                initSize_species,
1781
+##                                initSize_iter,
1782
+##                                seed,
1783
+##                                verbosity,
1784
+##                                speciesFS,
1785
+##                                ratioForce,
1786
+##                                typeFitness,
1787
+##                                max.memory,
1788
+##                                mutationPropGrowth, ## make it explicit
1789
+##                                initMutant,
1790
+##                                max.wall.time,
1791
+##                                keepEvery,
1792
+##                                alpha,
1793
+##                                sh,                               
1794
+##                                K,
1795
+##                                ## endTimeEvery,
1796
+##                                detectionDrivers,
1797
+##                                onlyCancer,
1798
+##                                errorHitWallTime,
1799
+##                                max.num.tries,
1800
+##                                errorHitMaxTries,
1801
+##                                minDetectDrvCloneSz,
1802
+##                                extraTime) {
1801 1803
 
1802
-    ## the value of 20000, in megabytes, for max.memory sets a limit of ~ 20 GB
1804
+##     ## the value of 20000, in megabytes, for max.memory sets a limit of ~ 20 GB
1803 1805
   
1804 1806
 
1805
-    ## if(keepEvery < sampleEvery)
1806
-    ##     warning("setting keepEvery to sampleEvery")
1807
-
1808
-    ## a backdoor to allow passing restrictionTables directly
1809
-    if(inherits(poset, "restrictionTable"))
1810
-        restrict.table <- poset
1811
-    else
1812
-        restrict.table <- poset.to.restrictTable(poset)
1813
-    numDrivers <- nrow(restrict.table)
1814
-    numGenes <- (numDrivers + numPassengers)
1815
-    if(numGenes < 2)
1816
-        stop("There must be at least two genes (loci) in the fitness effects.",
1817
-             "If you only care about a degenerate case with just one,",
1818
-             "you can enter a second gene",
1819
-             "with fitness effect of zero.")
1820
-    if(numGenes > 64)
1821
-        stop("Largest possible number of genes (loci) is 64 for version 1.",
1822
-             "You are strongly encouraged to use the new specification",
1823
-             "as in version 2.")
1824
-
1825
-    ## These can never be set by the user
1826
-    ## if(initSize_species < 10) {
1827
-    ##     warning("initSize_species too small?")
1828
-    ## }
1829
-    ## if(initSize_iter < 100) {
1830
-    ##     warning("initSize_iter too small?")
1831
-    ## }
1807
+##     ## if(keepEvery < sampleEvery)
1808
+##     ##     warning("setting keepEvery to sampleEvery")
1809
+
1810
+##     ## a backdoor to allow passing restrictionTables directly
1811
+##     if(inherits(poset, "restrictionTable"))
1812
+##         restrict.table <- poset
1813
+##     else
1814
+##         restrict.table <- poset.to.restrictTable(poset)
1815
+##     numDrivers <- nrow(restrict.table)
1816
+##     numGenes <- (numDrivers + numPassengers)
1817
+##     if(numGenes < 2)
1818
+##         stop("There must be at least two genes (loci) in the fitness effects.",
1819
+##              "If you only care about a degenerate case with just one,",
1820
+##              "you can enter a second gene",
1821
+##              "with fitness effect of zero.")
1822
+##     if(numGenes > 64)
1823
+##         stop("Largest possible number of genes (loci) is 64 for version 1.",
1824
+##              "You are strongly encouraged to use the new specification",
1825
+##              "as in version 2.")
1826
+
1827
+##     ## These can never be set by the user
1828
+##     ## if(initSize_species < 10) {
1829
+##     ##     warning("initSize_species too small?")
1830
+##     ## }
1831
+##     ## if(initSize_iter < 100) {
1832
+##     ##     warning("initSize_iter too small?")
1833
+##     ## }
1832 1834
 
1833
-    ## numDrivers <- nrow(restrict.table)
1834
-    if(length(unique(restrict.table[, 1])) != numDrivers)
1835
-        stop("BAIL OUT NOW: length(unique(restrict.table[, 1])) != numDrivers)")
1836
-    ddr <- restrict.table[, 1]
1837
-    if(any(diff(ddr) != 1))
1838
-        stop("BAIL OUT NOW:  any(diff(ddr) != 1")
1839
-    ## sanity checks
1840
-    if(max(restrict.table[, 1]) != numDrivers)
1841
-        stop("BAIL OUT NOW: max(restrict.table[, 1]) != numDrivers")
1842
-    if(numDrivers > numGenes)
1843
-        stop("BAIL OUT NOW: numDrivers > numGenes")
1835
+##     ## numDrivers <- nrow(restrict.table)
1836
+##     if(length(unique(restrict.table[, 1])) != numDrivers)
1837
+##         stop("BAIL OUT NOW: length(unique(restrict.table[, 1])) != numDrivers)")
1838
+##     ddr <- restrict.table[, 1]
1839
+##     if(any(diff(ddr) != 1))
1840
+##         stop("BAIL OUT NOW:  any(diff(ddr) != 1")
1841
+##     ## sanity checks
1842
+##     if(max(restrict.table[, 1]) != numDrivers)
1843
+##         stop("BAIL OUT NOW: max(restrict.table[, 1]) != numDrivers")
1844
+##     if(numDrivers > numGenes)
1845
+##         stop("BAIL OUT NOW: numDrivers > numGenes")
1844 1846
     
1845
-    non.dep.drivers <- restrict.table[which(restrict.table[, 2] == 0), 1]
1847
+##     non.dep.drivers <- restrict.table[which(restrict.table[, 2] == 0), 1]
1846 1848
 
1847 1849
 
1848 1850
 
1849 1851
 
1850
-    ## if( (is.null(endTimeEvery) || (endTimeEvery > 0)) &&
1851
-    ##    (typeFitness %in% c("bozic1", "exp") )) {
1852
-    ##     warning(paste("endTimeEvery will take a positive value. ",
1853
-    ##                   "This will make simulations not stop until the next ",
1854
-    ##                   "endTimeEvery has been reached. Thus, in simulations ",
1855
-    ##                   "with very fast growth, simulations can take a long ",
1856
-    ##                   "time to finish, or can hit the wall time limit. "))
1857
-    ## }
1858
-    ## if(is.null(endTimeEvery))
1859
-    ##     endTimeEvery <- keepEvery
1860
-    ## if( (endTimeEvery > 0) && (endTimeEvery %% keepEvery) )
1861
-    ##     warning("!(endTimeEvery %% keepEvery)")
1862
-    ## a sanity check in restricTable, so no neg. indices for the positive deps
1863
-    neg.deps <- function(x) {
1864
-        ## checks a row of restrict.table
1865
-        numdeps <- x[2]
1866
-        if(numdeps) {
1867
-            deps <- x[3:(3+numdeps - 1)]
1868
-            return(any(deps < 0))
1869
-        } else FALSE
1870
-    }
1871
-    if(any(apply(restrict.table, 1, neg.deps)))
1872
-        stop("BAIL OUT NOW: Negative dependencies in restriction table")
1852
+##     ## if( (is.null(endTimeEvery) || (endTimeEvery > 0)) &&
1853
+##     ##    (typeFitness %in% c("bozic1", "exp") )) {
1854
+##     ##     warning(paste("endTimeEvery will take a positive value. ",
1855
+##     ##                   "This will make simulations not stop until the next ",
1856
+##     ##                   "endTimeEvery has been reached. Thus, in simulations ",
1857
+##     ##                   "with very fast growth, simulations can take a long ",
1858
+##     ##                   "time to finish, or can hit the wall time limit. "))
1859
+##     ## }
1860
+##     ## if(is.null(endTimeEvery))
1861
+##     ##     endTimeEvery <- keepEvery
1862
+##     ## if( (endTimeEvery > 0) && (endTimeEvery %% keepEvery) )
1863
+##     ##     warning("!(endTimeEvery %% keepEvery)")
1864
+##     ## a sanity check in restricTable, so no neg. indices for the positive deps
1865
+##     neg.deps <- function(x) {
1866
+##         ## checks a row of restrict.table
1867
+##         numdeps <- x[2]
1868
+##         if(numdeps) {
1869
+##             deps <- x[3:(3+numdeps - 1)]
1870
+##             return(any(deps < 0))
1871
+##         } else FALSE
1872
+##     }
1873
+##     if(any(apply(restrict.table, 1, neg.deps)))
1874
+##         stop("BAIL OUT NOW: Negative dependencies in restriction table")
1873 1875
 
1874
-    ## transpose the table
1875
-    rtC <- convertRestrictTable(restrict.table)
1876
+##     ## transpose the table
1877
+##     rtC <- convertRestrictTable(restrict.table)
1876 1878
 
1877 1879
     
1878
-    return(c(
1879
-        BNB_Algo5(restrictTable = rtC,
1880
-        numDrivers = numDrivers,
1881
-        numGenes = numGenes,
1882
-        typeCBN_= typeCBN,
1883
-        s = s, 
1884
-        death = death,
1885
-        mu = mu,
1886
-        initSize = initSize,
1887
-        sampleEvery = sampleEvery,
1888
-        detectionSize = detectionSize,
1889
-        finalTime = finalTime,
1890
-        initSp = initSize_species,
1891
-        initIt = initSize_iter,
1892
-        seed = seed,
1893
-        verbosity = verbosity,
1894
-        speciesFS = speciesFS,
1895
-        ratioForce = ratioForce,
1896
-        typeFitness_ = typeFitness,
1897
-        maxram = max.memory,
1898
-        mutationPropGrowth = as.integer(mutationPropGrowth),
1899
-        initMutant = initMutant,
1900
-        maxWallTime = max.wall.time,
1901
-        keepEvery = keepEvery,
1902
-        sh = sh,
1903
-        K = K,
1904
-        detectionDrivers = detectionDrivers,
1905
-        onlyCancer = onlyCancer,
1906
-        errorHitWallTime = errorHitWallTime,
1907
-        maxNumTries = max.num.tries,
1908
-        errorHitMaxTries = errorHitMaxTries,
1909
-        minDetectDrvCloneSz = minDetectDrvCloneSz,
1910
-        extraTime = extraTime
1911
-    ),
1912
-    NumDrivers = numDrivers
1913
-             ))
1880
+##     return(c(
1881
+##         BNB_Algo5(restrictTable = rtC,
1882
+##         numDrivers = numDrivers,
1883
+##         numGenes = numGenes,
1884
+##         typeCBN_= typeCBN,
1885
+##         s = s, 
1886
+##         death = death,
1887
+##         mu = mu,
1888
+##         initSize = initSize,
1889
+##         sampleEvery = sampleEvery,
1890
+##         detectionSize = detectionSize,
1891
+##         finalTime = finalTime,
1892
+##         initSp = initSize_species,
1893
+##         initIt = initSize_iter,
1894
+##         seed = seed,
1895
+##         verbosity = verbosity,
1896
+##         speciesFS = speciesFS,
1897
+##         ratioForce = ratioForce,
1898
+##         typeFitness_ = typeFitness,
1899
+##         maxram = max.memory,
1900
+##         mutationPropGrowth = as.integer(mutationPropGrowth),
1901
+##         initMutant = initMutant,
1902
+##         maxWallTime = max.wall.time,
1903
+##         keepEvery = keepEvery,
1904
+##         sh = sh,
1905
+##         K = K,
1906
+##         detectionDrivers = detectionDrivers,
1907
+##         onlyCancer = onlyCancer,
1908
+##         errorHitWallTime = errorHitWallTime,
1909
+##         maxNumTries = max.num.tries,
1910
+##         errorHitMaxTries = errorHitMaxTries,
1911
+##         minDetectDrvCloneSz = minDetectDrvCloneSz,
1912
+##         extraTime = extraTime
1913
+##     ),
1914
+##     NumDrivers = numDrivers
1915
+##              ))
1914 1916
 
1915
-}
1917
+## }
1916 1918
 
1917 1919
 OncoSimulWide2Long <- function(x) {
1918 1920
     ## Put data in long format, for ggplot et al
... ...
@@ -5,10 +5,6 @@ nr_BNB_Algo5 <- function(rFE, mu_, death, initSize_, sampleEvery, detectionSize,
5 5
     .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu_, death, initSize_, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list)
6 6
 }
7 7
 
8
-BNB_Algo5 <- function(restrictTable, numDrivers, numGenes, typeCBN_, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime) {
9
-    .Call('OncoSimulR_BNB_Algo5', PACKAGE = 'OncoSimulR', restrictTable, numDrivers, numGenes, typeCBN_, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime)
10
-}
11
-
12 8
 evalRGenotype <- function(rG, rFE, spPop, verbose, prodNeg, calledBy_, currentTime) {
13 9
     .Call('OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, spPop, verbose, prodNeg, calledBy_, currentTime)
14 10
 }
... ...
@@ -1,3 +1,6 @@
1
+Changes in version 2.99.6 (2020-12-13):
2
+	- No longer v.1 functionality.
3
+
1 4
 Changes in version 2.99.5 (2020-12-18):
2 5
 	- Random timeouts when building in tokay2 (Windows, BioC); fixing
3 6
 	seed in vignette.
... ...
@@ -8,7 +8,7 @@
8 8
 }
9 9
 \description{
10 10
   
11
-  Conver the \code{pops.by.time} component from its "wide" format (with
11
+  Convert the \code{pops.by.time} component from its "wide" format (with
12 12
   one column for time, and as many columns as clones/genotypes) into
13 13
   "long" format, so that it can be used with other functions, for
14 14
   instance for plots.
... ...
@@ -43,16 +43,8 @@ OncoSimulWide2Long(x)
43 43
 }
44 44
 
45 45
 \examples{
46
-data(examplePosets)
47
-## An object of class oncosimul
48
-p705 <- examplePosets[["p705"]]
49
-p1 <- oncoSimulIndiv(p705)
50
-class(p1)
51
-lp1 <- OncoSimulWide2Long(p1)
52
-head(lp1)
53
-summary(lp1)
54
-
55
-## An object of class oncosimul2
46
+
47
+
56 48
 data(examplesFitnessEffects)
57 49
 
58 50
 sm <-  oncoSimulIndiv(examplesFitnessEffects$cbn1,
... ...
@@ -867,74 +867,13 @@ of the individual done, and the number of attempts and time used.}
867 867
 
868 868
 
869 869
 \seealso{
870
-  \code{\link{plot.oncosimul}}, \code{\link{examplePosets}},
870
+  \code{\link{plot.oncosimul}}, 
871 871
   \code{\link{samplePop}}, \code{\link{allFitnessEffects}}
872 872
   
873 873
 
874 874
 }
875 875
 \examples{
876 876
 
877
-#################################
878
-#####
879
-##### Examples using v.1
880
-#####
881
-#################################
882
-
883
-
884
-## use poset p701
885
-data(examplePosets)
886
-p701 <- examplePosets[["p701"]]
887
-
888
-## Exp Model
889
-
890
-b1 <- oncoSimulIndiv(p701)
891
-summary(b1)
892
-
893
-plot(b1, addtot = TRUE)
894
-
895
-## McFarland; use a small sampleEvery, but also a reasonable
896
-##   keepEvery.
897
-## We also modify mutation rate to values similar to those in the
898
-##   original paper.
899
-## Note that detectionSize will play no role
900
-## finalTime is large, since this is a slower process
901
-## initSize is set to 4000 so the default K is larger and we are likely
902
-## to reach cancer. Alternatively, set K = 2000.
903
-
904
-m1 <- oncoSimulIndiv(p701,
905
-                     model = "McFL",
906
-                     mu = 5e-7,
907
-                     initSize = 4000,
908
-                     sampleEvery = 0.025,
909
-                     finalTime = 15000,
910
-                     keepEvery = 10,
911
-                     onlyCancer = FALSE)
912
-plot(m1, addtot = TRUE, log = "")
913
-
914
-
915
-
916
-
917
-## Simulating 4 individual trajectories
918
-## (I set mc.cores = 2 to comply with --as-cran checks, but you
919
-##  should either use a reasonable number for your hardware or
920
-##  leave it at its default value).
921
-
922
-
923
-p1 <- oncoSimulPop(4, p701,
924
-                   keepEvery = 10,
925
-                   mc.cores = 2)
926
-summary(p1)
927
-samplePop(p1)
928
-
929
-p2 <- oncoSimulSample(4, p701)
930
-
931
-
932
-#########################################
933
-######
934
-###### Examples using v.2:
935
-######
936
-#########################################
937
-
938 877
 
939 878
 
940 879
 #### A model similar to the one in McFarland. We use 270 genes.
... ...
@@ -1001,7 +940,7 @@ object.size(oiP1)
1001 940
 
1002 941
 
1003 942
 
1004
-######## Using a poset for pancreatic cancer from Gerstung et al.
943
+######## Using an extended poset for pancreatic cancer from Gerstung et al.
1005 944
 ###      (s and sh are made up for the example; only the structure
1006 945
 ###       and names come from Gerstung et al.)
1007 946
 
... ...
@@ -1042,7 +981,7 @@ femuv <- allFitnessEffects(noIntGenes = ni)
1042 981
 oncoSimulIndiv(femuv, mu = muv)
1043 982
 }
1044 983
 
1045
-#########Frequency dependent fitness examples
984
+######### Frequency dependent fitness examples
1046 985
 
1047 986
 ## An example with cooperation. Presence of WT favours all clones
1048 987
 ## and all clones have a positive effect on themselves
... ...
@@ -324,34 +324,6 @@
324 324
   \code{\link{oncoSimulIndiv}}
325 325
 }
326 326
 \examples{
327
-\dontrun{
328
-data(examplePosets)
329
-p701 <- examplePosets[["p701"]]
330
-
331
-## Simulate and plot a single individual, including showing
332
-## Shannon's diversity index
333
-b1 <- oncoSimulIndiv(p701)
334
-plot(b1, addtot = TRUE, plotDiversity = TRUE)
335
-
336
-## A stacked area plot
337
-plot(b1, type = "stacked", plotDiversity = TRUE)
338
-
339
-## And what if I show a stream plot?
340
-plot(b1, type = "stream", plotDiversity = TRUE)
341
-
342
-## Simulate and plot 2 individuals
343
-## (I set mc.cores = 2 to comply with --as-cran checks, but you
344
-##  should either use a reasonable number for your hardware or
345
-##  leave it at its default value).
346
-
347
-p1 <- oncoSimulPop(2, p701, mc.cores = 2)
348
-
349
-par(mfrow = c(1, 2))
350
-plot(p1, ask = FALSE)
351
-
352
-## Stacked; we cannot log here, and harder to see patterns
353
-plot(p1, ask = FALSE, type = "stacked")
354
-}
355 327
 
356 328
 ## Show individual genotypes and drivers for an
357 329
 ## epistasis case with at most eight genotypes
... ...
@@ -41,6 +41,10 @@
41 41
   This specification of restrictions is for version 1. See
42 42
   \code{\link{allFitnessEffects}} for a much more flexible one for
43 43
   version 2. Both can be used with  \code{\link{oncoSimulIndiv}}.
44
+
45
+
46
+  Note that simulating using posets directly is no longer
47
+  supported. This function is left here only for historical purposes.
44 48
 }
45 49
 
46 50
 \references{
... ...
@@ -168,14 +168,26 @@ sampledGenotypes(y, genes = NULL)
168 168
 }
169 169
 
170 170
 \examples{
171
-data(examplePosets)
172
-p705 <- examplePosets[["p705"]]
171
+######## Using an extended poset for pancreatic cancer from Gerstung et al.
172
+###      (s and sh are made up for the example; only the structure
173
+###       and names come from Gerstung et al.)
174
+
175
+
176
+pancr <- allFitnessEffects(data.frame(parent = c("Root", rep("KRAS", 4), "SMAD4", "CDNK2A",
177
+                                          "TP53", "TP53", "MLL3"),
178
+                                      child = c("KRAS","SMAD4", "CDNK2A",
179
+                                          "TP53", "MLL3",
180
+                                          rep("PXDN", 3), rep("TGFBR2", 2)),
181
+                                      s = 0.15,
182
+                                      sh = -0.3,
183
+                                      typeDep = "MN"))
184
+
173 185
 
174 186
 ## (I set mc.cores = 2 to comply with --as-cran checks, but you
175 187
 ##  should either use a reasonable number for your hardware or
176 188
 ##  leave it at its default value).
177 189
 
178
-p1 <- oncoSimulPop(4, p705, mc.cores = 2)
190
+p1 <- oncoSimulPop(4, pancr, mc.cores = 2)
179 191
 (sp1 <- samplePop(p1))
180 192
 sampledGenotypes(sp1)
181 193
 
... ...
@@ -190,7 +202,7 @@ sampledGenotypes(sp2)
190 202
 
191 203
 ## Now single cell sampling
192 204
 
193
-r1 <- oncoSimulIndiv(p705)
205
+r1 <- oncoSimulIndiv(pancr)
194 206
 samplePop(r1, typeSample = "single")
195 207
 
196 208
 sampledGenotypes(samplePop(r1, typeSample = "single"))
197 209
deleted file mode 100644
... ...
@@ -1,2314 +0,0 @@
1
-//     Copyright 2013-2021 Ramon Diaz-Uriarte
2
-
3
-//     This program is free software: you can redistribute it and/or modify
4
-//     it under the terms of the GNU General Public License as published by
5
-//     the Free Software Foundation, either version 3 of the License, or
6
-//     (at your option) any later version.
7
-
8
-//     This program is distributed in the hope that it will be useful,
9
-//     but WITHOUT ANY WARRANTY; without even the implied warranty of
10
-//     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11
-//     GNU General Public License for more details.
12
-
13
-//     You should have received a copy of the GNU General Public License
14
-//     along with this program.  If not, see <http://www.gnu.org/licenses/>.
15
-
16
-
17
-#include "debug_common.h"
18
-#include "common_classes.h"
19
-#include "bnb_common.h"
20
-#include <cfloat> 
21
-#include <limits>
22
-#include <iostream>
23
-#include <random>
24
-#include <bitset>
25
-#include <set>
26
-#include <iterator>
27
-#include <map>
28
-#include <sstream>
29
-#include <string>
30
-#include <ctime>
31
-#include <sys/time.h> 
32
-
33
-#include <stdexcept>
34
-
35
-using namespace Rcpp;
36
-using std::vector;
37
-
38
-
39
-// To track if mutation is really much smaller than birth/death
40
-#define MIN_RATIO_MUTS
41
-#ifdef MIN_RATIO_MUTS
42
-// There is really no need for these to be globals?
43
-// Unless I wanted to use them inside some function. So leave as globals.
44
-double g_min_birth_mut_ratio = DBL_MAX;
45
-double g_min_death_mut_ratio = DBL_MAX;
46
-double g_tmp1 = DBL_MAX;
47
-#endif
48
-
49
-
50
-typedef std::bitset<64> Genotype64;
51
-
52
-
53
-// Format of restrictTable
54
-// - mutations in columns
55
-// - first row, the number
56
-// - second row, the number of dependencies
57
-//   - rest of rows, the id of the dependency
58
-//   - past number of dependencies: a -9
59
-// In fact, the first row is redundant. Leave it, just in case.
60
-
61
-
62
-// Genotypes: the first (or column 0) genotype is the all ceros.
63
-// would not be needed, but makes Algo5 a lot simpler.
64
-
65
-// But mutatedPos start at 0. 
66
-// Will need to add 1 when plotting and analyzing with R.
67
-
68
-// Ojo: typeCBN is going to be an int.
69
-// But from R we pass a string, and that determined the integer.
70
-
71
-
72
-static void fitness(spParamsP& tmpP,
73
-		    const spParamsP& parentP,
74
-		    const int& mutatedPos, 
75
-		    Rcpp::IntegerMatrix restrictTable,
76
-		    const std::string& typeCBN,
77
-		    const Genotype64& newGenotype,
78
-		    // const double& birthRate, 
79
-		    const double& s,
80
-		    // const double& death,
81
-		    const int& numDrivers,
82
-		    const std::string& typeFitness,
83
-		    // const double& genTime,
84
-		    // const double& adjust_fitness_B,
85
-		    const double& sh){
86
-  //const double& adjust_fitness_MF) {
87
-
88
-  using namespace Rcpp;
89
-  // Two pieces: split into two functions??
90
-  //    - checking restrictions
91
-  //    - returning actual fitness according
92
-
93
-
94
-  int numDependencies;
95
-  int sumDriversMet = 0;
96
-  int sumDriversNoMet = 0;
97
-  int sumDependenciesMet = 0;
98
-
99
-  // set appropriate defaults. Change only needed stuff.
100
-  tmpP.birth = parentP.birth;
101
-  tmpP.death = parentP.death;
102
-  tmpP.absfitness = parentP.absfitness;
103
-
104
-
105
-
106
-
107
-
108
-  //      **** Are driver constraints met? ***
109
-
110
-
111
-  // Two cases: same s, sh, sp or different ones. If same, return three
112
-  // integers: sumDriversMet, sumDriversNoMet, sumPassengers.  If
113
-  // different, return three vectors, filled with the non-zero
114
-  // entries. These vectors then are combined as dictated by the fintness
115
-  // functions.
116
-
117
-  // If same single s, sh, sp: function takes three integers. O.w. it
118
-  // takes three integer vectors.
119
-
120
-  
121
-
122
-  if(mutatedPos >= numDrivers) { //the new mutation is a passenger
123
-    return;
124
-  } else {
125
-    for(int m = 0; m < numDrivers; ++m) {
126
-      if( newGenotype[m] ) { // this m is mutated
127
-	const Rcpp::IntegerMatrix::Column thisRestrict = 
128
-	  restrictTable(_, m);
129
-	numDependencies = thisRestrict[1];
130
-	if(!numDependencies) { // this driver has no dependencies
131
-	  sumDriversMet++;
132
-#ifdef DEBUGZ
133
-	  Rcpp::Rcout << "\n No dependencies:  ";
134
-	  DP2(sumDriversMet);
135
-#endif
136
-
137
-	}
138
-	else {
139
-	  sumDependenciesMet = 0;
140
-	  for(int i = 2; i < (2 + numDependencies); i++) {
141
-	    sumDependenciesMet += newGenotype[ thisRestrict[i] ];
142
-	  }
143
-	  if( ( (typeCBN == "Multiple") && (sumDependenciesMet) ) ||
144
-	      ( (typeCBN == "CBN") && (sumDependenciesMet == numDependencies) )) {
145
-	    sumDriversMet++;   
146
-	  } else {
147
-	    sumDriversNoMet++;
148
-	  }
149
-	}
150
-      }
151
-    }
152
-  }
153
-
154
-#ifdef DEBUGZ
155
-  DP2(sumDriversMet);
156
-  DP2(sumDriversNoMet);
157
-  DP2(sh);
158
-  DP2(typeFitness);
159
-#endif
160
-
161
-  // if sh < 0 : we do not allow any unment dependencies.  
162
-  // if sh = 0: no penalty for unmet dependencies
163
-
164
-
165
-  // FIXME: why not just pass the birth and death rates, and combine them
166
-  // in arbitrary ways? Might even allow to pass on death and birth rates
167
-  // from R. Only need care when any are density dependent.
168
-
169
-
170
-  // Beware: doing it this way with Bozic1 is kind of questionable because
171
-  // if birth = 1, there is no immediate extinction. In fact, there never
172
-  // is.
173
-  if((sh < 0) && sumDriversNoMet) {
174
-    tmpP.absfitness = 0.0;
175
-    tmpP.death = 1.0;
176
-    tmpP.birth = 0.0; // this is what really matters so that
177
-    // the pop does not get added.
178
-    // Line with comment "fitness is 0"
179
-  } else {
180
-    if(typeFitness == "bozic1") {
181
-      tmpP.death = pow( 1.0 - s, sumDriversMet) * 
182
-	pow( 1.0 + sh, sumDriversNoMet);
183
-      tmpP.birth = 1.0;
184
-    // } else if (typeFitness == "bozic2") {
185
-    //   double pp = pow( 1.0 - s, sumDriversMet) * 
186
-    // 	pow( 1.0 + sh, sumDriversNoMet);
187
-    //   tmpP.birth = (1.0/genTime) * (1.0 - 0.5 * pp );
188
-    //   tmpP.death = (0.5/genTime) * pp;
189
-    // } else if(typeFitness == "beerenwinkel") {
190
-    //   // like Datta et al., 2013
191
-    //   tmpP.absfitness = pow(1.0 + s, sumDriversMet) * 
192
-    // 	pow( 1.0 - sh, sumDriversNoMet);
193
-    //   tmpP.birth = adjust_fitness_B * tmpP.absfitness;
194
-    // } else if(typeFitness == "mcfarland0") {
195
-    //   tmpP.absfitness = pow(1.0 + s, sumDriversMet) / 
196
-    // 	pow( 1.0 + sh, sumDriversNoMet);
197
-    //   tmpP.birth = adjust_fitness_MF * tmpP.absfitness;
198
-    // } else if(typeFitness == "mcfarland") {
199
-    //   tmpP.birth = pow(1.0 + s, sumDriversMet) / 
200
-    // 	pow( 1.0 + sh, sumDriversNoMet);
201
-    } else if(typeFitness == "mcfarlandlog") {
202
-      tmpP.birth = pow(1.0 + s, sumDriversMet) / 
203
-	pow( 1.0 + sh, sumDriversNoMet);
204
-    } else if (typeFitness == "exp") { 
205
-      // Also like Datta et al., 2013 An additional driver gene mutation
206
-      // increases a cell’s fitness by a factor of (1+sd), whereas an
207
-      // additional housekeeper gene mutation decreases fitness by a
208
-      // factor of (1-sh) and the effect of multiple mutations is
209
-      // multiplicative
210
-      tmpP.birth = pow(1.0 + s, sumDriversMet) * 
211
-	pow( 1.0 - sh, sumDriversNoMet);
212
-
213
-#ifdef DEBUGZ
214
-      double posi = pow(1.0 + s, sumDriversMet);
215
-      double negi = pow( 1.0 - sh, sumDriversNoMet);
216
-      DP2(posi);
217
-      DP2(negi);
218
-#endif
219
-
220
-    } // else if (typeFitness == "log") {
221
-    //   tmpP.birth = birthRate+ s * log1p(sumDriversMet) - 
222
-    // 	sh * log(1 + sumDriversNoMet);
223
-    // } else { // linear
224
-    //   tmpP.birth = birthRate + s * static_cast<double>(sumDriversMet) - 
225
-    // 	sh * static_cast<double>(sumDriversNoMet);
226
-    // } 
227
-  }
228
-}
229
-// Notice: if restriction is 3 -> 4 -> 5
230
-// and one has 5 and 4, only 4 is unmet. Beware of that.
231
-// So we talk about the immediate dependency or restriction.
232
-// Not the whole transitive closure.
233
-
234
-// When birth == 0, popSize should become 0 immediately. 
235
-// No evaluation through random numbers, etc.
236
-// This is how we do it.
237
-
238
-// How small can they get?
239
-// d1 <- function(s, mut) { (1 - s)^mut}
240
-// d2 <- function(s, mut) {  (0.5/4) * ((1 - s)^mut) }
241
-
242
-
243
-// limited benchmarks suggest the following is slower
244
-// static inline void new_sp_bitset2(unsigned int& sp, const Genotype64& newGenotype,
245
-// 			   const std::vector<Genotype64>& Genotypes) {
246
-//   sp = std::distance(Genotypes.begin(),
247
-// 		     std::find(Genotypes.begin(), 
248
-// 			       Genotypes.end(), newGenotype));
249
-// }
250
-
251
-
252
-static inline void new_sp_bitset(unsigned int& sp, const Genotype64& newGenotype,
253
-			  const std::vector<Genotype64>& Genotypes) {
254
-  sp = 0;
255
-
256
-  for(sp = 0; sp < Genotypes.size(); ++sp) {
257
-    if( newGenotype == Genotypes[sp] )
258
-      break;
259
-  }
260
-}
261
-
262
-
263
-
264
-static void getMutatedPos_bitset(int& mutatedPos, int& numMutablePosParent,
265
-				 //gsl_rng *r,
266
-				 std::mt19937& ran_generator, 
267
-				 std::vector<int>& mutablePos,
268
-				 const Genotype64& nextMutantGenotype,
269
-				 // const int& nextMutant,
270
-				 // const std::vector<Genotype64>& Genotypes,
271
-				 const int& numGenes) {
272
-  // We want mutatedPos and numMutablePosParent
273
-
274
-  // Note: impossible to have a second recorded mutation in
275
-  // the same gene.  
276
-
277
-  // Remember numMutablePosParent is the number of mutable positions in
278
-  // the parent!  so after mutation is one less, but we do not decrease it
279
-  // here.
280
-  
281
-  numMutablePosParent = 0;
282
-  for(int i = 0; i < numGenes; ++i) {
283
-    if( !nextMutantGenotype.test(i) ) { 
284
-      mutablePos[numMutablePosParent] = i;
285
-      ++numMutablePosParent;
286
-    }
287
-  }
288
-  
289
-  if(numMutablePosParent > 1) {
290
-    std::uniform_int_distribution<int> unif(0, numMutablePosParent - 1);
291
-    mutatedPos = mutablePos[unif(ran_generator)];
292
-  } else {
293
-    mutatedPos = mutablePos[0];
294
-  } 
295
-
296
-  // if(numMutablePosParent > 1) {
297
-  //   mutatedPos = mutablePos[gsl_rng_uniform_int(r, numMutablePosParent)];
298
-  // } else {
299
-  //   mutatedPos = mutablePos[0];
300
-  // } 
301
-
302
-
303
-#ifdef DEBUGV
304
-      Rcpp::Rcout << "\n numMutablePosParent = " << numMutablePosParent;
305
-      Rcpp::Rcout << "\n mutatedPos = " << mutatedPos  << "\n";
306
-      
307
-#endif
308
-
309
-  // if(numMutablePos > 1) {
310
-  //   mutatedPos = mutablePos[gsl_rng_uniform_int(r, numMutablePos)];
311
-  // } else if (numMutablePos == 1) {
312
-  //   mutatedPos = mutablePos[0];
313
-  // } else {
314
-  //   // Should never happen, as mutation = 0 if no mutable positions.
315
-  //   throw std::out_of_range("Algo5: run out of mutable places!!??");
316
-  // }
317
-
318
-}
319
-
320
-
321
-static void remove_zero_sp_v7(std::vector<int>& sp_to_remove,
322
-			      std::vector<Genotype64>& Genotypes,
323
-			      std::vector<spParamsP>& popParams,
324
-			      std::multimap<double, int>& mapTimes) {
325
-  //  here("entering remove_zero_sp_v7");
326
-  std::vector<spParamsP>::iterator popParams_begin = popParams.begin();
327
-  std::vector<Genotype64>::iterator Genotypes_begin = Genotypes.begin();
328
-  std::vector<int>::reverse_iterator r = sp_to_remove.rbegin();
329
-  int remove_this;
330
-  // for(r = sp_to_remove.rbegin(); r != sp_to_remove.rend(); ++r) {
331
-  while(r != sp_to_remove.rend() ) {
332
-    remove_this = *r;
333
-    mapTimes.erase(popParams[remove_this].pv);
334
-    popParams.erase(popParams_begin + remove_this);
335
-    Genotypes.erase(Genotypes_begin + remove_this);
336
-    ++r;
337
-  }
338
-  // here("exiting remove_zero_sp_v7");
339
-
340
-}
341
-
342
-static inline int count_NDrivers(const Genotype64& Genotype,
343
-				 const int& NumDrivers) {
344
-  int totalDr = 0;
345
-  for(int i = 0; i < NumDrivers; ++i)
346
-    totalDr += Genotype[i];
347
-  return totalDr;
348
-}
349
-
350
-static void totPopSize_and_fill_out_crude_P(int& outNS_i,
351
-					    double& totPopSize, 
352
-					    double& lastStoredSample,
353
-					    std::vector<Genotype64>& genot_out,
354
-					    //std::vector<unsigned long long>& sp_id_out,
355
-					    std::vector<double>& popSizes_out,
356
-					    std::vector<int>& index_out,
357
-					    std::vector<double>& time_out,
358
-					    std::vector<double>& sampleTotPopSize,
359
-					    std::vector<double>& sampleLargestPopSize,
360
-					    std::vector<int>& sampleMaxNDr,
361
-					    std::vector<int>& sampleNDrLargestPop,
362
-					    bool& simulsDone,
363
-					    bool& reachDetection,
364
-					    int& lastMaxDr,
365
-					    double& done_at,
366
-					    const std::vector<Genotype64>& Genotypes,
367
-					    const std::vector<spParamsP>& popParams, 
368
-					    const double& currentTime,
369
-					    const int& NumDrivers,
370
-					    const double& keepEvery,
371
-					    const double& detectionSize,
372
-					    const double& finalTime,
373
-					    // const double& endTimeEvery,
374
-					    const int& detectionDrivers,
375
-					    const int& verbosity,
376
-					    const double& minDetectDrvCloneSz,
377
-					    const double& extraTime,
378
-					    const double& fatalPopSize = 1e15) {
379
-  // Fill out, but also compute totPopSize
380
-  // and return sample summaries for popsize, drivers.
381
-  
382
-  // This determines if we are done or not by checking popSize, number of
383
-  // drivers, etc
384
-  
385
-  // static int lastMaxDr = 0; // preserves value across calls to Algo5 from R.
386
-  // so can not use it.
387
-  bool storeThis = false;
388
-  totPopSize = 0.0;
389
-  
390
-   // DP2(lastMaxDr);
391
-  // DP2(detectionDrivers);
392
-  // DP2(currentTime);
393
-  // DP2((lastStoredSample + endTimeEvery));
394
-  // DP2(detectionSize);
395
-
396
-  // this could all be part of popSize_over_m_dr, with a better name
397
-  int tmp_ndr = 0;
398
-  int max_ndr = 0;
399
-  double popSizeOverDDr = 0.0;
400
-
401
-  for(size_t i = 0; i < popParams.size(); ++i) {
402
-    totPopSize += popParams[i].popSize;
403
-    tmp_ndr = count_NDrivers(Genotypes[i], NumDrivers);
404
-    if(tmp_ndr > max_ndr) max_ndr = tmp_ndr;
405
-    if(tmp_ndr >= detectionDrivers) popSizeOverDDr += popParams[i].popSize;
406
-  }
407
-  lastMaxDr = max_ndr;
408
-
409
-  
410
-  if (keepEvery < 0) {
411
-    storeThis = false;
412
-  } else if( currentTime >= (lastStoredSample + keepEvery) ) {
413
-    storeThis = true;
414
-  }
415
-
416
-  if( (totPopSize <= 0.0) || (currentTime >= finalTime)  ) {
417
-    simulsDone = true;
418
-  }
419
-
420
-
421
-  // if( (totPopSize >= detectionSize) ||
422
-  //     ( (lastMaxDr >= detectionDrivers) &&
423
-  //       (popSizeOverDDr >= minDetectDrvCloneSz) ) ) {
424
-  //   simulsDone = true;
425
-  //   reachDetection = true;
426
-  // }
427
-
428
-  if(extraTime > 0) {
429
-    if(done_at <  0) {
430
-      if( (totPopSize >= detectionSize) ||
431
-	  ( (lastMaxDr >= detectionDrivers) &&
432
-	    (popSizeOverDDr >= minDetectDrvCloneSz) ) ) {
433
-	done_at = currentTime + extraTime;
434
-      }
435
-    } else if (currentTime >= done_at) {
436
-	simulsDone = true;
437
-	reachDetection = true; 
438
-      }
439
-  } else if( (totPopSize >= detectionSize) ||
440
-	     ( (lastMaxDr >= detectionDrivers) &&
441
-	       (popSizeOverDDr >= minDetectDrvCloneSz) ) ) {
442
-    simulsDone = true;
443
-    reachDetection = true; 
444
-  }
445
-  
446
- 
447
-
448
-    
449
-  // This is no longer used.
450
-  // // Beware: this can lead to never stopping if
451
-  // // decreases in popSize or drivers
452
-
453
-  // // Logic: if a period k you meet any condition, recheck again at k +
454
-  // // endTimeEvery, and if conditions met exit. Prevents exiting if you
455
-  // // reach the cancer state almost by chance. But this is way too
456
-  // // paranoid. The idea is of application mainly for McF and Beeren
457
-  // // models, so we do not bail out as soon as just a single cell with one
458
-  // // new driver. But this makes things very slow.
459
-
460
-  // // Thus, never pass an endTimeEvery > 0, but use detectionDrivers = 1 +
461
-  // // intended final Drivers.
462
-
463
-  // // FIXME
464
-  // // Ideally, we would check, for McFL, that popsize of the pop with
465
-  // // required number of drivers is at least, say, > initSize.
466
-  // // But that is not trivial, as requires more accounting. Do later.
467
-
468
-  
469
-  // if(endTimeEvery > 0) {
470
-  //   if(done_at <= 0 ) {
471
-  //     if( (totPopSize >= detectionSize) ||
472
-  // 	   (lastMaxDr >= detectionDrivers)  )
473
-  // 	done_at = currentTime + endTimeEvery;
474
-  //   } else if (currentTime >= done_at) {
475
-  //     if( (totPopSize >= detectionSize) ||
476
-  // 	  (lastMaxDr >= detectionDrivers)  ) {
477
-  // 	simulsDone = true;
478
-  // 	reachDetection = true;
479
-  //     }
480
-  //     else
481
-  // 	done_at = -9;
482
-  //   }
483
-  // } else if( (totPopSize >= detectionSize) ||
484
-  // 	     (lastMaxDr >= detectionDrivers) )  {
485
-    
486
-  //     simulsDone = true;
487
-  //     reachDetection = true;
488
-  // }
489
-
490
-  
491
-  if(totPopSize >= fatalPopSize) {
492
-    Rcpp::Rcout << "\n\totPopSize > " << fatalPopSize
493
-		<<". You are likely to loose precision and run into numerical issues\n";
494
-       }
495
-  
496
-  if(simulsDone)
497
-    storeThis = true;
498
-
499
-
500
-  if( storeThis ) {
501
-    lastStoredSample = currentTime;
502
-    outNS_i++;
503
-    int ndr_lp = 0;
504
-    double l_pop_s = 0.0;
505
-    
506
-    time_out.push_back(currentTime);
507
-    
508
-    for(size_t i = 0; i < popParams.size(); ++i) {
509
-      genot_out.push_back(Genotypes[i]);
510
-      popSizes_out.push_back(popParams[i].popSize);
511
-      index_out.push_back(outNS_i);
512
-      
513
-      if(popParams[i].popSize > l_pop_s) {
514
-	l_pop_s = popParams[i].popSize;
515
-	ndr_lp = count_NDrivers(Genotypes[i], NumDrivers);
516
-      }
517
-    }
518
-    sampleTotPopSize.push_back(totPopSize);
519
-    sampleLargestPopSize.push_back(l_pop_s);
520
-    sampleMaxNDr.push_back(max_ndr);
521
-    sampleNDrLargestPop.push_back(ndr_lp);
522
-  } 
523
-
524
-
525
-  
526
-  // if( storeThis ) {
527
-  //   lastStoredSample = currentTime;
528
-  //   outNS_i++;
529
-  //   int tmp_ndr = 0;
530
-  //   int max_ndr = 0;
531
-  //   int ndr_lp = 0;
532
-  //   double l_pop_s = 0.0;
533
-    
534
-  //   time_out.push_back(currentTime);
535
-    
536
-  //   for(size_t i = 0; i < popParams.size(); ++i) {
537
-  //     genot_out.push_back(Genotypes[i]);
538
-  //     popSizes_out.push_back(popParams[i].popSize);
539
-  //     index_out.push_back(outNS_i);
540
-  //     // I have to repeat the counting of drivers here.
541
-  //     tmp_ndr = count_NDrivers(Genotypes[i], NumDrivers); 
542
-  //     if(tmp_ndr > max_ndr) max_ndr = tmp_ndr;
543
-  //     if(popParams[i].popSize > l_pop_s) {
544
-  // 	l_pop_s = popParams[i].popSize;
545
-  // 	ndr_lp = tmp_ndr;
546
-  // 	// ndr_lp = count_NDrivers(Genotypes[i], NumDrivers); 
547
-  //     }
548
-  //     // lastMaxDr = max_ndr; // and this should have been out of the
549
-  //     // popParams.size() loop
550
-  //   }
551
-  //   // lastMaxDr = max_ndr;
552
-  //   sampleTotPopSize.push_back(totPopSize);
553
-  //   sampleLargestPopSize.push_back(l_pop_s);
554
-  //   sampleMaxNDr.push_back(max_ndr);
555
-  //   sampleNDrLargestPop.push_back(ndr_lp);
556
-  // }//  else if (keepEvery < 0) {
557
-  //   // FIXME keepEvery
558
-  //   // must keep track of results to bail out
559
-
560
-  //   // FIXME counting max drivers should be done always, like counting
561
-  //   // totPopSize.
562
-    
563
-  //   int tmp_ndr = 0;
564
-  //   int max_ndr = 0;
565
-   
566
-  //   for(size_t i = 0; i < popParams.size(); ++i) {
567
-  //     tmp_ndr = count_NDrivers(Genotypes[i], NumDrivers);
568
-  //     if(tmp_ndr > max_ndr) max_ndr = tmp_ndr;
569
-  //     // lastMaxDr = max_ndr;
570
-  //   }
571
-  //   lastMaxDr = max_ndr;
572
-  // }
573
-
574
-  
575
-    
576
-  
577
-  if( !std::isfinite(totPopSize) ) {
578
-    throw std::range_error("totPopSize not finite");
579
-  }
580
-  if( std::isnan(totPopSize) ) {
581
-    throw std::range_error("totPopSize is NaN");
582
-  }
583
-  
584
-  if(totPopSize > (4.0 * 1e15)) {
585
-    if(verbosity > 0)
586
-      Rcpp::Rcout << "\nWARNING: popSize > 4e15. Likely loss of precission\n";
587
-  }
588
-}
589
-
590
-// FIXME: I might want to return the actual drivers in each period
591
-// and the actual drivers in the population with largest popsize
592
-// Something like what we do now with whichDrivers
593
-// and count_NumDrivers
594
-
595
-
596
-// static inline void fill_SStats(Rcpp::NumericMatrix& perSampleStats,
597
-// 			       const std::vector<double>& sampleTotPopSize,
598
-// 			       const std::vector<double>& sampleLargestPopSize,
599
-// 			       const std::vector<double>& sampleLargestPopProp,
600
-// 			       const std::vector<int>& sampleMaxNDr,
601
-// 			       const std::vector<int>& sampleNDrLargestPop){
602
-
603
-//   for(size_t i = 0; i < sampleTotPopSize.size(); ++i) {
604
-//     perSampleStats(i, 0) = sampleTotPopSize[i];
605
-//     perSampleStats(i, 1) = sampleLargestPopSize[i]; // Never used in R FIXME: remove!!
606
-//     perSampleStats(i, 2) = sampleLargestPopProp[i]; // Never used in R
607
-//     perSampleStats(i, 3) = static_cast<double>(sampleMaxNDr[i]);
608
-//     perSampleStats(i, 4) = static_cast<double>(sampleNDrLargestPop[i]);
609
-//   }
610
-// }
611
-
612
-inline void reshape_to_outNS(Rcpp::NumericMatrix& outNS,
613
-			     const std::vector<unsigned long long>& uniqueGenotV,
614
-			     const std::vector<unsigned long long>& genot_out_ul,
615
-			     const std::vector<double>& popSizes_out,
616
-			     const std::vector<int>& index_out,
617
-			     const std::vector<double>& time_out){
618
-  
619
-  std::vector<unsigned long long>::const_iterator fbeg = uniqueGenotV.begin();
620
-  std::vector<unsigned long long>::const_iterator fend = uniqueGenotV.end();
621
-
622
-  int column;
623
-
624
-  for(size_t i = 0; i < genot_out_ul.size(); ++i) {
625
-    column = std::distance(fbeg, lower_bound(fbeg, fend, genot_out_ul[i]) );
626
-    // here("   looping over i ");
627
-    outNS(index_out[i], column + 1) =  popSizes_out[i];
628
-  }
629
-
630
-  for(size_t j = 0; j < time_out.size(); ++j)
631
-    outNS(j, 0) = time_out[j];
632
-}
633
-
634
-static inline void find_unique_genotypes(std::set<unsigned long long>& uniqueGenotypes,
635
-				  const std::vector<unsigned long long>& genot_out_l) {
636
-  for(size_t i = 0; i < genot_out_l.size(); ++i) 
637
-    uniqueGenotypes.insert( genot_out_l[i] );
638
-}
639
-
640
-static inline void genot_out_to_ullong(std::vector<unsigned long long>& go_l,
641
-			       const std::vector<Genotype64>& go) {
642
-  for(size_t i = 0; i < go.size(); ++i)
643
-    go_l[i] = go[i].to_ullong();
644
-}
645
-
646
-
647
-static inline void uniqueGenotypes_to_vector(std::vector<unsigned long long>& ugV,
648
-				      const std::set<unsigned long long>& uniqueGenotypes) {
649
-  ugV.assign(uniqueGenotypes.begin(), uniqueGenotypes.end() );
650
-}
651
-
652
-
653
-static inline void create_returnGenotypes(Rcpp::IntegerMatrix& returnGenotypes,
654
-					  const int& numGenes,
655
-					  const std::vector<unsigned long long>& uniqueGenotypesV){
656
-  // In C++, as the original were bitsets, pos 0 is at the right
657
-  // In R, pos 0 is at the left
658
-
659
-  for(size_t i = 0; i < uniqueGenotypesV.size(); ++i) {
660
-    Genotype64 tmpbs(uniqueGenotypesV[i]);
661
-    for(int j = 0; j < numGenes; ++j) {
662
-      returnGenotypes(j, i) = tmpbs[j];
663
-    }
664
-  }
665
-}
666
-
667
-// FIXME: change this, now that we keep a count of drivers?
668
-static inline void count_NumDrivers(int& maxNumDrivers, 
669
-				    std::vector<int>& countByDriver,
670
-				    Rcpp::IntegerMatrix& returnGenotypes,
671
-				    const int& numDrivers){
672
-  //				    Rcpp::IntegerVector& totDrivers){
673
-
674
-  // At the end, over all genotypes
675
-  // Using a returnGenotypes object as input
676
-  // Redundant? 
677
-  maxNumDrivers = 0;
678
-  int tmpdr = 0;
679
-  
680
-  for(int j = 0; j < returnGenotypes.ncol(); ++j) {
681
-    tmpdr = 0;
682
-    for(int i = 0; i < numDrivers; ++i) {
683
-      tmpdr += returnGenotypes(i, j);
684
-      countByDriver[i] += returnGenotypes(i, j);
685
-    }
686
-    // totDrivers(j) = tmpdr;
687
-    if(tmpdr > maxNumDrivers) maxNumDrivers = tmpdr;
688
-  }
689
-}
690
-
691
-static inline void whichDrivers(int& totalPresentDrivers,
692
-				std::string& strDrivers,
693
-				const std::vector<int>& countByDriver){
694
-  std::string comma = "";
695
-  for(size_t i = 0; i < countByDriver.size(); ++i) {
696
-    if(countByDriver[i] > 0) {
697
-      strDrivers += (comma + std::to_string(i + 1)); //SSTR(i + 1));
698
-// #ifdef _WIN32  
699
-//       strDrivers += (comma + SSTR(i + 1));
700
-// #endif
701
-// #ifndef _WIN32
702
-//       strDrivers += (comma + std::to_string(i + 1)); //SSTR(i + 1));
703
-// #endif      
704
-      comma = ", ";
705
-      ++totalPresentDrivers;
706
-    }
707
-  }
708
-  if(totalPresentDrivers == 0) strDrivers = "NA";
709
-}
710
-
711
-
712
-static void sample_all_pop_P(std::vector<int>& sp_to_remove,
713
-			     std::vector<spParamsP>& popParams,
714
-			     // next only used with DEBUGV
715
-			     const std::vector<Genotype64>& Genotypes,
716
-			     const double& tSample,
717
-			     const int& mutationPropGrowth){
718
-
719
-  // here("entering sample_all_pop_P");
720
-  // currentTime = tSample;
721
-  sp_to_remove.clear();
722
-  // sp_to_remove.push_back(0);
723
-  // sp_to_remove[0] = 0;
724
-
725
-  for(size_t i = 0; i < popParams.size(); i++) {
726
-    //STOPASSERT(popParams[i].Flag == false);
727
-    STOPASSERT(popParams[i].timeLastUpdate >= 0.0);
728
-    STOPASSERT(tSample - popParams[i].timeLastUpdate >= 0.0);
729
-#ifdef DEBUGV
730
-    Rcpp::Rcout << "\n\n     ********* 5.9 ******\n " 
731
-	      << "     Species  = " << i 
732
-	      << "\n      Genotype = " << Genotypes[i]
733
-	      << "\n      sp_id = " << Genotypes[i].to_ullong() // sp_id[i]  
734
-	      << "\n      pre-update popSize = " 
735
-	      << popParams[i].popSize 
736
-	      << "\n      time of sample = " << tSample 
737
-	      << "\n      popParams[i].timeLastUpdate = " 
738
-	      << popParams[i].timeLastUpdate 
739
-	      << ";\n     t for Algo2 = " 
740
-	      << tSample - popParams[i].timeLastUpdate 
741
-	      << " \n     species R " << popParams[i].R
742
-	      << " \n     species W " << popParams[i].W
743
-	      << " \n     species death " << popParams[i].death
744
-	      << " \n     species birth " << popParams[i].birth;
745
-    // << " \n     species nextMutationTime " 
746
-    // << popParams[i].nextMutationTime;
747
-#endif
748
-
749
-    // Account for forceSampling. When 
750
-    // forceSampling, popSize for at least one species
751
-    // was updated in previous loop, so we skip that one
752
-    if(tSample > popParams[i].timeLastUpdate) {
753
-      popParams[i].popSize = 
754
-	Algo2_st(popParams[i], tSample, mutationPropGrowth);
755
-    }
756
-    if( popParams[i].popSize <=  0.0 ) {
757
-      // this i has never been non-zero in any sampling time
758
-      // eh??
759
-
760
-      // If it is 0 here, remove from _current_ population. Anything that
761
-      // has had a non-zero size at sampling time is preserved (if it
762
-      // needs to be preserved, because it is keepEvery time).
763
-
764
-      // sp_to_remove[0]++;
765
-      sp_to_remove.push_back(i);
766
-      // sp_to_remove[sp_to_remove[0]] = i;
767
-#ifdef DEBUGV
768
-      Rcpp::Rcout << "\n\n     Removing species i = " << i 
769
-		<< " with sp_id = " << Genotypes[i].to_ullong(); //sp_id[i];
770
-#endif
771
-    } // else {
772
-    //   popParams[i].Flag = true;
773
-    // }
774
-#ifdef DEBUGV
775
-    Rcpp::Rcout << "\n\n   post-update popSize = " 
776
-	      << popParams[i].popSize << "\n";
777
-#endif
778
-  }
779
-  // here("exiting sample_all_pop_P");
780
-}
781
-
782
-
783
-static void innerBNB(const int& numGenes,
784
-		     const double& initSize,
785
-		     const double& K,
786
-		     // const double& alpha,
787
-		     const std::string& typeCBN,
788
-		     //		     const double& genTime,
789
-		     const std::string& typeFitness,
790
-		     const int& mutationPropGrowth,
791
-		     const double& mu,
792
-		     const double& sh,
793
-		     const double& s,
794
-		     const double& death,
795
-		     // const double& birthRate,
796
-		     const double& keepEvery,
797
-		     const double& sampleEvery,		     
798
-		     const int& numDrivers,
799
-		     const int& initMutant,
800
-		     const time_t& start_time,
801
-		     const double& maxWallTime,
802
-		     const double& finalTime,
803
-		     const double& detectionSize,
804
-		     // const double& endTimeEvery,
805
-		     const int& detectionDrivers,
806
-		     const double& minDetectDrvCloneSz,
807
-		     const double& extraTime,
808
-		     const int& verbosity,
809
-		     double& totPopSize,
810
-		     double& em1,
811
-		     double& em1sc,
812
-		     // double& n_1,
813
-		     // double& en1,
814
-		     double& ratioForce,
815
-		     double& currentTime,
816
-		     int& speciesFS,
817
-		     int& outNS_i,
818
-		     int& iter,
819
-		     std::vector<Genotype64>& genot_out,
820
-		     std::vector<double>& popSizes_out,
821
-		     std::vector<int>& index_out,
822
-		     std::vector<double>& time_out,
823
-		     std::vector<double>& sampleTotPopSize,
824
-		     std::vector<double>& sampleLargestPopSize,
825
-		     std::vector<int>& sampleMaxNDr,
826
-		     std::vector<int>& sampleNDrLargestPop,
827
-		     bool& reachDetection,
828
-		     std::mt19937& ran_generator,
829
-		     double& runningWallTime,
830
-		     bool& hittedWallTime,
831
-		     Rcpp::IntegerMatrix restrictTable) {
832
-		     //bool& anyForceRerunIssues
833
-  //  if(numRuns > 0) {
834
-
835
-  double dummyMutationRate = std::max(mu/1000, 1e-13);
836
-  // double dummyMutationRate = 1e-10;
837
-  // ALWAYS initialize this here, or reinit or rezero
838
-  genot_out.clear();
839
-  popSizes_out.clear();
840
-  index_out.clear();
841
-  time_out.clear();
842
-  totPopSize = 0.0;
843
-  sampleTotPopSize.clear();
844
-  currentTime = 0.0;
845
-  iter = 0;
846
-
847
-  outNS_i = -1;
848
-
849
-  sampleTotPopSize.clear();
850
-  sampleLargestPopSize.clear();
851
-  sampleMaxNDr.clear();
852
-  sampleNDrLargestPop.clear();
853
-  // end of rezeroing.
854
-
855
-  
856
-  // }
857
-  // anyForceRerunIssues = false;
858
-  
859
-  bool forceSample = false;
860
-  bool simulsDone = false;
861
-  double lastStoredSample = 0.0;
862
-
863
-
864
-  double minNextMutationTime;
865
-  double mutantTimeSinceLastUpdate;
866
-  double timeNextPopSample;
867
-  double tSample;
868
-
869
-  int nextMutant;
870
-  unsigned int numSpecies = 0;
871
-  int numMutablePosParent = 0;
872
-  int mutatedPos = 0;
873
-  //int indexMutatedPos = 0;
874
-
875
-  unsigned int sp = 0;
876
-  //int type_resize = 0;
877
-
878
-  int iterL = 1000;
879
-  int speciesL = 1000; 
880
-  //int timeL = 1000;
881
-  
882
-  double tmpdouble1 = 0.0;
883
-  double tmpdouble2 = 0.0;
884
-  
885
-  std::vector<int>sp_to_remove(1);
886
-  sp_to_remove.reserve(10000);
887
-
888
-  // those to update
889
-  int to_update = 1; //1: one species; 2: 2 species; 3: all.
890
-  int u_1 = -99;
891
-  int u_2 = -99;
892
-
893
-  Genotype64 newGenotype;
894
-  std::vector<Genotype64> Genotypes(1);
895
-  Genotypes[0].reset();
896
-  
897
-  std::vector<spParamsP> popParams(1);
898
-
899
-      // // FIXME debug
900
-      // Rcpp::Rcout << "\n popSize[0]  at 1 ";
901
-      // print_spP(popParams[0]);
902
-      // // end debug
903
-  
904
-  
905
-  const int sp_per_period = 5000;
906
-  popParams.reserve(sp_per_period);
907
-  Genotypes.reserve(sp_per_period);
908
-
909
-
910
-      // // FIXME debug
911
-      // Rcpp::Rcout << "\n popSize[0]  at 01 ";
912
-      // print_spP(popParams[0]);
913
-      // // end debug
914
-
915
-  
916
-  spParamsP tmpParam; 
917
-  init_tmpP(tmpParam);
918
-  init_tmpP(popParams[0]);
919
-  popParams[0].popSize = initSize;
920
-  totPopSize = initSize;
921
-
922
-
923
-      // // FIXME debug
924
-      // Rcpp::Rcout << "\n popSize[0]  at 10000 ";
925
-      // print_spP(popParams[0]);
926
-      // // end debug
927
-
928
-
929
-  
930
-  std::vector<int>mutablePos(numGenes); // could be inside getMuatedPos_bitset
931
-
932
-    // multimap to hold nextMutationTime
933
-  std::multimap<double, int> mapTimes;
934
-  //std::multimap<double, int>::iterator m1pos;
935
-
936
-  int ti_dbl_min = 0;
937
-  int ti_e3 = 0;
938
-
939
-
940
-      // // FIXME debug
941
-      // Rcpp::Rcout << "\n popSize[0]  at 10002 ";
942
-      // print_spP(popParams[0]);
943
-      // // end debug
944
-
945
-
946
-  
947
-  // // Beerenwinkel
948
-  // double adjust_fitness_B = -std::numeric_limits<double>::infinity();
949
-  //McFarland
950
-  double adjust_fitness_MF = -std::numeric_limits<double>::infinity();
951
-
952
-  // for McFarland error
953
-  em1 = 0.0;
954
-  em1sc = 0.0;
955
-  // n_0 = 0.0;
956
-  // n_1 = 0.0;
957
-  // double tps_0; //, tps_1; 
958
-  // tps_0 = totPopSize;
959
-  // tps_1 = totPopSize;
960
-
961
-  
962
-  double totPopSize_previous = totPopSize;
963
-  double DA_previous = log1p(totPopSize_previous/K);
964
-  
965
-      // // FIXME debug
966
-      // Rcpp::Rcout << "\n popSize[0]  at 10004 ";
967
-      // print_spP(popParams[0]);
968
-      // // end debug
969
-
970
-  
971
-  
972
-  int lastMaxDr = 0;
973
-  double done_at = -9;
974
-
975
-#ifdef MIN_RATIO_MUTS
976
-  g_min_birth_mut_ratio = DBL_MAX;
977
-  g_min_death_mut_ratio = DBL_MAX;
978
-  g_tmp1 = DBL_MAX;
979
-#endif
980
-
981
-      // // FIXME debug
982
-      // Rcpp::Rcout << " popSize[0]  at 1b ";
983
-      // print_spP(popParams[0]);
984
-      // // end debug
985
-
986
-    // This long block, from here to X1, is ugly and a mess!
987
-  // This is what takes longer to figure out whenever I change
988
-  // anything. FIXME!!
989
-
990
-  if(initMutant >= 0)
991
-    throw std::invalid_argument("initMutant no longer allowed. Bug in R code.");
992
-  // if(initMutant >= 0) {
993
-  //   popParams[0].numMutablePos = numGenes - 1;
994
-  //   Genotypes[0].set(initMutant);
995
-  //   // if(typeFitness == "beerenwinkel") {
996
-  //   //   popParams[0].death = 1.0; //note same is in McFarland.
997
-  //   //   // But makes sense here; adjustment in beerenwinkel is via fitness
998
-      
999
-  //   //   // initialize to prevent birth/mutation warning with Beerenwinkel
1000
-  //   //   // when no mutator. O.w., the defaults
1001
-  //   //   if(!mutationPropGrowth)
1002
-  //   // 	popParams[0].mutation = mu * popParams[0].numMutablePos;
1003
-  //   //   popParams[0].absfitness = 1.0 + s;
1004
-  //   //   updateRatesBeeren(popParams, adjust_fitness_B, initSize,
1005
-  //   // 			currentTime, alpha, initSize, 
1006
-  //   // 			mutationPropGrowth, mu);
1007
-  //   // } else if(typeFitness == "mcfarland0") {
1008
-  //   //   // death equal to birth of a non-mutant.
1009
-  //   //   popParams[0].death = log1p(totPopSize/K); // log(2.0), except rare cases
1010
-  //   //   if(!mutationPropGrowth)
1011
-  //   // 	popParams[0].mutation = mu * popParams[0].numMutablePos;
1012
-  //   //   popParams[0].absfitness = 1.0 + s;
1013
-  //   //   updateRatesMcFarland0(popParams, adjust_fitness_MF, K, 
1014
-  //   // 			    totPopSize,
1015
-  //   // 			    mutationPropGrowth, mu);
1016
-  //   // } else if(typeFitness == "mcfarland") {
1017
-  //   //   popParams[0].death = totPopSize/K;
1018
-  //   //   popParams[0].birth = 1.0 + s;
1019
-  //   // } else if(typeFitness == "mcfarlandlog") {
1020
-  //   if(typeFitness == "mcfarlandlog") {      
1021
-  //     popParams[0].death = log1p(totPopSize/K);
1022
-  //     popParams[0].birth = 1.0 + s;
1023
-  //   } else if(typeFitness == "bozic1") {
1024
-  //     tmpParam.birth =  1.0;
1025
-  //     tmpParam.death = -99.9;
1026
-  //   // } else if (typeFitness == "bozic2") {
1027
-  //   //   tmpParam.birth =  -99;
1028
-  //   //   tmpParam.death = -99;
1029
-  //   } else if (typeFitness == "exp") {
1030
-  //     tmpParam.birth =  -99;
1031
-  //     tmpParam.death = death;
1032
-  //   } // else { // linear or log
1033
-  //   //   tmpParam.birth =  -99;
1034
-  //   //   tmpParam.death = death;
1035
-  //   // } 
1036
-  //   // if( (typeFitness != "beerenwinkel") && (typeFitness != "mcfarland0") 
1037
-  //   // 	&& (typeFitness != "mcfarland") && (typeFitness != "mcfarlandlog")) // wouldn't matter
1038
-  //   //   fitness(popParams[0], tmpParam, initMutant, restrictTable,
1039
-  //   // 	      typeCBN, Genotypes[0], birthRate, s, numDrivers, 
1040
-  //   // 	      typeFitness, genTime, adjust_fitness_B, sh,
1041
-  //   // 	      adjust_fitness_MF);
1042
-
1043
-  //   if( typeFitness != "mcfarlandlog") // wouldn't matter
1044
-  //     fitness(popParams[0], tmpParam, initMutant, restrictTable,
1045
-  // 	      typeCBN, Genotypes[0], // birthRate,
1046
-  // 	      s, numDrivers, 
1047
-  // 	      typeFitness, // genTime, adjust_fitness_B,
1048
-  // 	      sh);
1049
-  //   // adjust_fitness_MF);
1050
-  //   // we pass as the parent the tmpParam; it better initialize
1051
-  //   // everything right, or that will blow. Reset to init
1052
-  //   init_tmpP(tmpParam);
1053
-  // } else {
1054
-    popParams[0].numMutablePos = numGenes;
1055
-    // if(typeFitness == "beerenwinkel") {
1056
-    //   popParams[0].death = 1.0;
1057
-    //   // initialize to prevent birth/mutation warning with Beerenwinkel
1058
-    //   // when no mutator. O.w., the defaults
1059
-    //   if(!mutationPropGrowth)
1060
-    // 	popParams[0].mutation = mu * popParams[0].numMutablePos;
1061
-    //   popParams[0].absfitness = 1.0;
1062
-    //   updateRatesBeeren(popParams, adjust_fitness_B, initSize,
1063
-    // 			currentTime, alpha, initSize, 
1064
-    // 			mutationPropGrowth, mu);
1065
-    // } else if(typeFitness == "mcfarland0") {
1066
-    //   popParams[0].death = log1p(totPopSize/K);
1067
-    //   if(!mutationPropGrowth)
1068
-    // 	popParams[0].mutation = mu * popParams[0].numMutablePos;
1069
-    //   popParams[0].absfitness = 1.0;
1070
-    //   updateRatesMcFarland0(popParams, adjust_fitness_MF, K, 
1071
-    // 			    totPopSize,
1072
-    // 			    mutationPropGrowth, mu);
1073
-    // } else if(typeFitness == "mcfarland") {
1074
-    //   popParams[0].birth = 1.0;
1075
-    //   popParams[0].death = totPopSize/K;
1076
-    //   // no need to call updateRates
1077
-    // } else if(typeFitness == "mcfarlandlog") {
1078
-    if(typeFitness == "mcfarlandlog") {      
1079
-      popParams[0].birth = 1.0;
1080
-      popParams[0].death = log1p(totPopSize/K);
1081
-      // no need to call updateRates
1082
-    } else if(typeFitness == "bozic1") {
1083
-      popParams[0].birth = 1.0;
1084
-      popParams[0].death = 1.0;
1085
-    // } else if (typeFitness == "bozic2") {
1086
-    //   popParams[0].birth = 0.5/genTime;
1087
-    //   popParams[0].death = 0.5/genTime;
1088
-    } else if (typeFitness == "exp") {
1089
-      popParams[0].birth = 1.0;
1090
-      popParams[0].death = death;
1091
-    } // else { // linear or log
1092
-    //   popParams[0].birth = birthRate;
1093
-    //   popParams[0].death = death;
1094
-    // }
1095
-    //  }
1096
-
1097
-
1098
-  // // FIXME debug
1099
-  //     Rcpp::Rcout << " popSize[0]  at 2 ";
1100
-  //     print_spP(popParams[0]);
1101
-  //     // end debug
1102
-  
1103
-
1104
-  
1105
-  // these lines (up to, and including, R_F_st)
1106
-  // not needed with mcfarland0 or beerenwinkel
1107
-  if(mutationPropGrowth)
1108
-    popParams[0].mutation = mu * popParams[0].birth * popParams[0].numMutablePos;
1109
-  else
1110
-    popParams[0].mutation = mu * popParams[0].numMutablePos;
1111
-
1112
-  W_f_st(popParams[0]);
1113
-  R_f_st(popParams[0]);
1114
-
1115
-  // // FIXME debug
1116
-  //     Rcpp::Rcout << " popSize[0]  at 3 ";
1117
-  //     print_spP(popParams[0]);
1118
-  //     // end debug
1119
-  
1120
-
1121
-  // X1: end of mess of initialization block
1122
-
1123
-  popParams[0].pv = mapTimes.insert(std::make_pair(-999, 0));
1124
-
1125
-  if( keepEvery > 0 ) {
1126
-    // We keep the first ONLY if we are storing more than one.
1127
-    outNS_i++;
1128
-    time_out.push_back(currentTime);
1129
-    
1130
-    genot_out.push_back(Genotypes[0]);
1131
-    popSizes_out.push_back(popParams[0].popSize);
1132
-    index_out.push_back(outNS_i);
1133
-
1134
-    sampleTotPopSize.push_back(popParams[0].popSize);
1135
-    sampleLargestPopSize.push_back(popParams[0].popSize);
1136
-    sampleMaxNDr.push_back(count_NDrivers(Genotypes[0], numDrivers));
1137
-    sampleNDrLargestPop.push_back(sampleMaxNDr[0]);
1138
-  }
1139
-  // FIXME: why next line and not just genot_out.push_back(Genotypes[i]);
1140
-  // if keepEvery > 0? We do that already.
1141
-  // It is just ugly to get a 0 in that first genotype when keepEvery < 0
1142
-  // uniqueGenotypes.insert(Genotypes[0].to_ullong());
1143
-  timeNextPopSample = currentTime + sampleEvery;
1144
-  numSpecies = 1;
1145
-
1146
-  
1147
-#ifdef DEBUGV
1148
-  Rcpp::Rcout << "\n the initial species\n";
1149
-  print_spP(popParams[0]);
1150
-#endif
1151
-
1152
-
1153
-  // // FIXME debug
1154
-  //     Rcpp::Rcout << " popSize[0]  at 4 ";
1155
-  //     print_spP(popParams[0]);
1156
-  //     // end debug
1157
-  
1158
-
1159
-  
1160
-  while(!simulsDone) {
1161
-    // Check how we are doing with time as first thing.
1162
-    runningWallTime = difftime(time(NULL), start_time);
1163
-    if( runningWallTime > maxWallTime ) {
1164
-      hittedWallTime = true;
1165
-      forceSample = true;
1166
-      simulsDone = true;
1167
-    }
1168
-    
1169
-    iter++;
1170
-    if(verbosity > 1) {
1171
-      if(! (iter % iterL) ) {
1172
-	Rcpp::Rcout << "\n\n    ... iteration " << iter;
1173
-	Rcpp::Rcout << "\n    ... currentTime " << currentTime <<"\n";
1174
-      }
1175
-      if(!(numSpecies % speciesL )) {
1176
-	Rcpp::Rcout << "\n\n    ... iteration " << iter;
1177
-	Rcpp::Rcout << "\n\n    ... numSpecies " << numSpecies << "\n";
1178
-      }
1179
-    }
1180
-
1181
-    //  ************   5.2   ***************
1182
-    if(verbosity >= 2)
1183
-      Rcpp::Rcout <<"\n\n\n*** Looping through 5.2. Iter = " << iter << " \n";
1184
-
1185
-    tSample = std::min(timeNextPopSample, finalTime);
1186
-
1187
-#ifdef DEBUGV
1188
-    Rcpp::Rcout << " DEBUGV\n";
1189
-    Rcpp::Rcout << "\n ForceSample? " << forceSample 
1190
-	      << "  tSample " << tSample 
1191
-	      << "  currentTime " << currentTime;
1192
-#endif
1193
-
1194
-    if(iter == 1) { // handle special case of first iter
1195
-      // // FIXME debug
1196
-      // Rcpp::Rcout << " popSize[0] ";
1197
-      // print_spP(popParams[0]);
1198
-      // // end debug
1199
-      tmpdouble1 = ti_nextTime_tmax_2_st(popParams[0],
1200
-					 currentTime,
1201
-					 tSample, 
1202
-					 ti_dbl_min, ti_e3);
1203
-      mapTimes_updateP(mapTimes, popParams, 0, tmpdouble1);
1204
-      //popParams[0].Flag = false;
1205
-      popParams[0].timeLastUpdate = currentTime;
1206
-    } else { // any other iter
1207
-      if(to_update == 1) {
1208
-	// we did not sample in previous period.
1209
-	tmpdouble1 = ti_nextTime_tmax_2_st(popParams[u_1],
1210
-					   currentTime,
1211
-					   tSample, 
1212
-					   ti_dbl_min, ti_e3);
1213
-	mapTimes_updateP(mapTimes, popParams, u_1, tmpdouble1);
1214
-	popParams[u_1].timeLastUpdate = currentTime;
1215
-
1216
-#ifdef DEBUGV
1217
-	Rcpp::Rcout << "\n\n     ********* 5.2: call to ti_nextTime ******\n For to_update = \n " 
1218
-		  << "     tSample  = " << tSample
1219
-	    
1220
-		  << "\n\n**   Species  = " << u_1 
1221
-		  << "\n       genotype =  " << Genotypes[u_1] 
1222
-		  << "\n       popSize = " << popParams[u_1].popSize 
1223
-		  << "\n       currentTime = " << currentTime 
1224
-		  << "\n       popParams[i].nextMutationTime = " 
1225
-		  << tmpdouble1
1226
-		  << " \n     species R " << popParams[u_1].R
1227
-		  << " \n     species W " << popParams[u_1].W
1228
-		  << " \n     species death " << popParams[u_1].death
1229
-		  << " \n     species birth " << popParams[u_1].birth;
1230
-#endif
1231
-
1232
-      } else if(to_update == 2) {
1233
-	// we did not sample in previous period.
1234
-	tmpdouble1 = ti_nextTime_tmax_2_st(popParams[u_1],
1235
-					   currentTime,
1236
-					   tSample, ti_dbl_min, ti_e3);
1237
-	mapTimes_updateP(mapTimes, popParams, u_1, tmpdouble1);
1238