Browse code

2.15.1: Added MAGELLANs sources and functionality from MAGELLAN

ramon diaz-uriarte (at Phelsuma) authored on 02/07/2019 14:55:40
Showing76 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.15.0
5
-Date: 2019-03-18
4
+Version: 2.15.1
5
+Date: 2019-06-06
6 6
 Authors@R: c(person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),
7 7
 		     email = "rdiaz02@gmail.com"),
8 8
 	      person("Mark", "Taylor", role = "ctb", email = "ningkiling@gmail.com"))
... ...
@@ -10,6 +10,7 @@ export("oncoSimulPop", "oncoSimulIndiv", "samplePop",
10 10
        "rfitness",
11 11
        "plotFitnessLandscape",
12 12
        "to_Magellan",
13
+       "Magellan_stats",
13 14
        "sampledGenotypes"
14 15
        , "POM", "LOD"
15 16
        , "diversityPOM", "diversityLOD"
... ...
@@ -78,7 +78,7 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
78 78
 
79 79
         ## Might not be needed with the proper gfm object (so gmf <- x)
80 80
         ## but is needed if arbitrary matrices.
81
-        gfm <- allGenotypes_to_matrix(afe) 
81
+        gfm <- allGenotypes_to_matrix(afe)
82 82
     } else if(inherits(x, "fitnessEffects")) {
83 83
         if(!is.null(x$orderEffects) )
84 84
             stop("We cannot yet deal with order effects")
... ...
@@ -105,16 +105,16 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
105 105
         ## Assume a two-column data frame of genotypes as character
106 106
         ## vectors and fitness
107 107
         if(colnames(x)[2] != "Fitness")
108
-            stop("We cannot guess what you are passing here") 
108
+            stop("We cannot guess what you are passing here")
109 109
         afe <- evalAllGenotypes(allFitnessEffects(genotFitness = x),
110 110
                                 order = FALSE, addwt = TRUE,
111 111
                                 max = max_num_genotypes)
112 112
         gfm <- allGenotypes_to_matrix(afe)
113 113
     } else {
114
-        stop("We cannot guess what you are passing here") 
114
+        stop("We cannot guess what you are passing here")
115 115
     }
116 116
     return(list(gfm = gfm, afe = afe))
117
-}   
117
+}
118 118
 
119 119
 ## Based on from_genotype_fitness
120 120
 ## but if we are passed a fitness landscapes as produced by
... ...
@@ -125,17 +125,17 @@ to_genotFitness_std <- function(x, simplify = TRUE,
125 125
                                 sort_gene_names = TRUE) {
126 126
     ## Would break with output from allFitnessEffects and
127 127
     ## output from allGenotypeAndMut
128
-    
128
+
129 129
     ## For the very special and weird case of
130 130
     ## a matrix but only a single gene so with a 0 and 1
131 131
     ## No, this is a silly and meaningless case.
132 132
     ## if( ( ncol(x) == 2 ) && (nrow(x) == 1) && (x[1, 1] == 1) ) {
133
-    
134
-    ## } else  blabla: 
135
-    
133
+
134
+    ## } else  blabla:
135
+
136 136
     if(! (inherits(x, "matrix") || inherits(x, "data.frame")) )
137 137
         stop("Input must inherit from matrix or data.frame.")
138
-    
138
+
139 139
     ## if((ncol(x) > 2) && !(inherits(x, "matrix"))
140 140
     ##     stop(paste0("Genotype fitness input either two-column data frame",
141 141
     ##          " or a numeric matrix with > 2 columns."))
... ...
@@ -143,7 +143,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
143 143
     ##     stop(paste0("It looks like you have a matrix for a single genotype",
144 144
     ##                 " of a single gene. For this degenerate cases use",
145 145
     ##                 " a data frame specification."))
146
-    
146
+
147 147
     if(ncol(x) > 2) {
148 148
         if(inherits(x, "matrix")) {
149 149
             if(!is.numeric(x))
... ...
@@ -152,14 +152,14 @@ to_genotFitness_std <- function(x, simplify = TRUE,
152 152
             if(!all(unlist(lapply(x, is.numeric))))
153 153
                 stop("A genotype fitness matrix/data.frame must be numeric.")
154 154
         }
155
-        
155
+
156 156
         ## We are expecting here a matrix of 0/1 where columns are genes
157 157
         ## except for the last column, that is Fitness
158 158
         ## Of course, can ONLY work with epistastis, NOT order
159 159
         ## return(genot_fitness_to_epistasis(x))
160 160
         if(any(duplicated(colnames(x))))
161 161
             stop("duplicated column names")
162
-        
162
+
163 163
         cnfl <- which(colnames(x)[-ncol(x)] == "")
164 164
         if(length(cnfl)) {
165 165
             freeletter <- setdiff(LETTERS, colnames(x))[1]
... ...
@@ -176,7 +176,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
176 176
                 x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
177 177
             }
178 178
         }
179
-        
179
+
180 180
         if(is.null(colnames(x))) {
181 181
             ncx <- (ncol(x) - 1)
182 182
             message("No column names: assigning gene names from LETTERS")
... ...
@@ -185,7 +185,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
185 185
                      " as you see fit.")
186 186
             colnames(x) <- c(LETTERS[1:ncx], "Fitness")
187 187
         }
188
-        
188
+
189 189
         if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
190 190
             stop("First ncol - 1 entries not in {0, 1}.")
191 191
     } else {
... ...
@@ -202,7 +202,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
202 202
             stop(paste0("genotFitness: first column of data frame is numeric.",
203 203
                         " Ambiguous and suggests possible error. If sure,",
204 204
                         " enter that column as character"))
205
-        
205
+
206 206
         omarker <- any(grepl(">", x[, 1], fixed = TRUE))
207 207
         emarker <- any(grepl(",", x[, 1], fixed = TRUE))
208 208
         nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))
... ...
@@ -262,17 +262,17 @@ to_genotFitness_std <- function(x, simplify = TRUE,
262 262
 ## from_genotype_fitness <- function(x) {
263 263
 ##     ## Would break with output from allFitnessEffects and
264 264
 ##     ## output from allGenotypeAndMut
265
-    
265
+
266 266
 ##     ## For the very special and weird case of
267 267
 ##     ## a matrix but only a single gene so with a 0 and 1
268 268
 ##     ## No, this is a silly and meaningless case.
269 269
 ##     ## if( ( ncol(x) == 2 ) && (nrow(x) == 1) && (x[1, 1] == 1) ) {
270
-    
271
-##     ## } else  blabla: 
272
-    
270
+
271
+##     ## } else  blabla:
272
+
273 273
 ##     if(! (inherits(x, "matrix") || inherits(x, "data.frame")) )
274 274
 ##         stop("Input must inherit from matrix or data.frame.")
275
-    
275
+
276 276
 ##     ## if((ncol(x) > 2) && !(inherits(x, "matrix"))
277 277
 ##     ##     stop(paste0("Genotype fitness input either two-column data frame",
278 278
 ##     ##          " or a numeric matrix with > 2 columns."))
... ...
@@ -280,7 +280,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
280 280
 ##     ##     stop(paste0("It looks like you have a matrix for a single genotype",
281 281
 ##     ##                 " of a single gene. For this degenerate cases use",
282 282
 ##     ##                 " a data frame specification."))
283
-    
283
+
284 284
 ##     if(ncol(x) > 2) {
285 285
 ##         if(inherits(x, "matrix")) {
286 286
 ##             if(!is.numeric(x))
... ...
@@ -289,7 +289,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
289 289
 ##             if(!all(unlist(lapply(x, is.numeric))))
290 290
 ##                 stop("A genotype fitness matrix/data.frame must be numeric.")
291 291
 ##         }
292
-        
292
+
293 293
 ##         ## We are expecting here a matrix of 0/1 where columns are genes
294 294
 ##         ## except for the last column, that is Fitness
295 295
 ##         ## Of course, can ONLY work with epistastis, NOT order
... ...
@@ -304,7 +304,7 @@ to_genotFitness_std <- function(x, simplify = TRUE,
304 304
 ##             stop(paste0("genotFitness: first column of data frame is numeric.",
305 305
 ##                         " Ambiguous and suggests possible error. If sure,",
306 306
 ##                         " enter that column as character"))
307
-        
307
+
308 308
 ##         omarker <- any(grepl(">", x[, 1], fixed = TRUE))
309 309
 ##         emarker <- any(grepl(",", x[, 1], fixed = TRUE))
310 310
 ##         nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))
... ...
@@ -356,7 +356,7 @@ genot_fitness_to_epistasis <- function(x) {
356 356
     ## you use this is because you say "this is the mapping genotype ->
357 357
     ## fitness. Period." so we should not combine other terms (or other
358 358
     ## terms that involve these genes)
359
-    
359
+
360 360
     nr <- nrow(x)
361 361
     if(nr < (2^(ncol(x) - 1)))
362 362
         message("Number of genotypes less than 2^L.",
... ...
@@ -379,7 +379,7 @@ genot_fitness_to_epistasis <- function(x) {
379 379
                 "Dividing all fitnesses by fitness of wildtype.")
380 380
         f <- f/fwt
381 381
     }
382
-    
382
+
383 383
     if(is.null(colnames(x)) || any(grepl("^$", colnames(x))) ) {
384 384
         message("Setting/resetting gene names because one or more are missing.",
385 385
                 " If this is not what you want, pass a matrix",
... ...
@@ -403,10 +403,10 @@ genot_fitness_to_epistasis <- function(x) {
403 403
 
404 404
 
405 405
 
406
-allGenotypes_to_matrix <- function(x) { 
406
+allGenotypes_to_matrix <- function(x) {
407 407
     ## Makes no sense to allow passing order: the matrix would have
408 408
     ## repeated rows. A > B and B > A both have exactly A and B
409
-    
409
+
410 410
     ## Take output of evalAllGenotypes or identical data frame and return
411 411
     ## a matrix with 0/1 in a column for each gene and a final column of
412 412
     ## Fitness
... ...
@@ -452,15 +452,88 @@ allGenotypes_to_matrix <- function(x) {
452 452
 }
453 453
 
454 454
 
455
+
456
+Magellan_stats <- function(x, max_num_genotypes = 2000,
457
+                           verbose = FALSE,
458
+                           use_log = TRUE,
459
+                           short = TRUE,
460
+                           replace_missing = FALSE) {
461
+    ## I always use
462
+    ## if(!is.null(x) && is.null(file))
463
+    ##     stop("one of object or file name")
464
+    ## if(is.null(file))
465
+    fn <- tempfile()
466
+    fnret <- tempfile()
467
+    if(verbose)
468
+        cat("\n Using input file", fn, " and output file ", fnret, "\n")
469
+
470
+    if(use_log) {
471
+        logarg <- "-l"
472
+    } else {
473
+        logarg <- NULL
474
+    }
475
+    if(short) {
476
+        shortarg <- "-s"
477
+    } else {
478
+        shortarg <- NULL
479
+    }
480
+
481
+    if(replace_missing) {
482
+        zarg <- "-z"
483
+    } else {
484
+        zarg <- NULL
485
+    }
486
+
487
+    to_Magellan(x, fn, max_num_genotypes = max_num_genotypes)
488
+    ## newer versions of Magellan provide some extra values to standard output
489
+    call_M <- system2(fl_statistics_binary(),
490
+                      args = paste(shortarg, logarg, zarg, "-o", fnret, fn),
491
+                      stdout = NULL)
492
+    if(short) {
493
+        ## tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[-1])
494
+
495
+        tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[c(-1)])
496
+        ## ## Make names more explicit, but check we have what we think we have
497
+        ## ## New versions of Magellan produce different output apparently of variable length
498
+        ## stopifnot(length(tmp) >= 23) ## 23) ## variable length
499
+        ## stopifnot(identical(names(tmp)[1:13], ## only some
500
+        ##                     c("ngeno", "npeaks", "nsinks", "gamma", "gamma.", "r.s",
501
+        ##                       "nchains", "nsteps", "nori", "depth", "magn", "sign",
502
+        ##                       "rsign"))) ## , "w.1.", "w.2.", "w.3..", "mode_w", "s.1.",
503
+        ## ## "s.2.", "s.3..", "mode_s", "pairs_s", "outD_v")))
504
+        ## if(length(tmp) >= 24) ## the new version
505
+        ##     stopifnot(identical(names(tmp)[c(20, 24)],
506
+        ##                         c("steps_m", "mProbOpt_0")))
507
+        ## ## steps_m: the mean number of steps over the entire landscape to
508
+        ## ## reach the global optimum
509
+        ## ## mProbOpt_0: The mean probability to
510
+        ## ## reach that optimum (again averaged over the entire landscape).
511
+        ## ## but if there are multiple optima, there are many of these
512
+        ## names(tmp)[1:13] <- c("n_genotypes", "n_peaks", "n_sinks", "gamma", "gamma_star",
513
+        ##                 "r/s","n_chains", "n_steps", "n_origins", "max_depth",
514
+        ##                 "epist_magn", "epist_sign", "epist_rsign")## ,
515
+        ##                 ## "walsh_coef_1", "walsh_coef_2", "walsh_coef_3", "mode_walsh",
516
+        ##                 ## "synerg_coef_1", "synerg_coef_2", "synerg_coef_3", "mode_synerg",
517
+        ## ## "std_dev_pairs", "var_outdegree")
518
+    } else {
519
+        message("Output file: ", fnret)
520
+        tmp <- readLines(fnret)
521
+    }
522
+    return(tmp)
523
+}
524
+
525
+
526
+## Former version, that always tries to give a vector
527
+## and that uses an external executable.
455 528
 ## Magellan_stats and Magellan_draw cannot be tested
456 529
 ## routinely, as they depend on external software
457
-Magellan_stats <- function(x, max_num_genotypes = 2000,
530
+Magellan_stats_former <- function(x, max_num_genotypes = 2000,
458 531
                            verbose = FALSE,
459 532
                            use_log = TRUE,
460 533
                            short = TRUE,
461 534
                            fl_statistics = "fl_statistics",
462 535
                            replace_missing = FALSE) { # nocov start
463
-    ## I always use 
536
+    ## I always use
464 537
     ## if(!is.null(x) && is.null(file))
465 538
     ##     stop("one of object or file name")
466 539
     ## if(is.null(file))
... ...
@@ -468,7 +541,7 @@ Magellan_stats <- function(x, max_num_genotypes = 2000,
468 541
     fnret <- tempfile()
469 542
     if(verbose)
470 543
         cat("\n Using input file", fn, " and output file ", fnret, "\n")
471
-    
544
+
472 545
     if(use_log) {
473 546
         logarg <- "-l"
474 547
     } else {
... ...
@@ -479,7 +552,7 @@ Magellan_stats <- function(x, max_num_genotypes = 2000,
479 552
     } else {
480 553
         shortarg <- NULL
481 554
     }
482
-    
555
+
483 556
     if(replace_missing) {
484 557
         zarg <- "-z"
485 558
     } else {
... ...
@@ -493,7 +566,7 @@ Magellan_stats <- function(x, max_num_genotypes = 2000,
493 566
                       stdout = NULL)
494 567
     if(short) {
495 568
         ## tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[-1])
496
-        
569
+
497 570
         tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[c(-1)])
498 571
         ## Make names more explicit, but check we have what we think we have
499 572
         ## New versions of Magellan produce different output apparently of variable length
... ...
@@ -538,14 +611,14 @@ Magellan_draw <- function(x, max_num_genotypes = 2000,
538 611
     fn_out <- paste0(fn, ".svg")
539 612
     if(verbose)
540 613
         cat("\n Using input file", fn, " and output file ", fn_out, "\n")
541
-       
614
+
542 615
     to_Magellan(x, fn, max_num_genotypes = max_num_genotypes)
543 616
     call_M <- system2(fl_draw, args = paste(fn, args), wait = FALSE)
544 617
     call_view <- system2(svg_open, args = fn_out, wait = FALSE,
545 618
                          stdout = ifelse(verbose, "", FALSE),
546 619
                          stderr = ifelse(verbose, "", FALSE))
547
-    
548
-    invisible() 
620
+
621
+    invisible()
549 622
 } # nocov end
550 623
 
551 624
 
... ...
@@ -14,6 +14,48 @@
14 14
 ## along with this program.  If not, see <http://www.gnu.org/licenses/>.
15 15
 
16 16
 
17
+create_eq_ref <- function(g) {
18
+    ## "random" gives more prob. to genotypes with
19
+    ## number of mutated genes close to g/2.
20
+    ## This gives equal prob to having the reference
21
+    ## be of any of the possible number of mutated genes.
22
+    nm <- sample(g, 1)
23
+    ref <- c(rep(1, nm), rep(0, g - nm))
24
+    sample(ref)
25
+}
26
+
27
+
28
+
29
+get_magellan_binaries <-
30
+    function(names_binaries = c("fl_statistics", "fl_generate", "fl_genchains"))
31
+{
32
+    rarch <- Sys.getenv('R_ARCH')
33
+    nn_names_binaries <- names_binaries
34
+    if(.Platform$OS.type == "windows")
35
+        names_binaries <- paste0(names_binaries, ".exe")
36
+    if(nzchar(rarch)) {
37
+        rarch <- sub("^/", "", rarch)
38
+        magellan_binaries <-  system.file(package = "OncoSimulR", "exec",
39
+                                          rarch, names_binaries)
40
+    } else {
41
+        magellan_binaries <-  system.file(package = "OncoSimulR", "exec",
42
+                                          names_binaries)
43
+    }
44
+    names(magellan_binaries) <- nn_names_binaries
45
+    return(magellan_binaries)
46
+}
47
+
48
+## pkg.env <- new.env()
49
+## ## The next will not install with error
50
+## ## ERROR: hard-coded installation path:
51
+## please report to the package maintainer and use ‘--no-staged-install’
52
+## pkg.env <- c(pkg.env, get_magellan_binaries())
53
+
54
+
55
+
56
+fl_statistics_binary <- function() get_magellan_binaries("fl_statistics")
57
+fl_generate_binary <- function() get_magellan_binaries("fl_generate")
58
+
17 59
 
18 60
 rfitness <- function(g, c= 0.5,
19 61
                      sd = 1,
... ...
@@ -21,7 +63,7 @@ rfitness <- function(g, c= 0.5,
21 63
                      reference = "random", ## "random", "max", or the vector,
22 64
                                      ## e.g., rep(1, g). If random, a
23 65
                                      ## random genotype is chosen as
24
-                                     ## reference. If "max" this is rep(1, g)
66
+                     ## reference. If "max" this is rep(1, g)
25 67
                      scale = NULL, ## a two-element vector: min and max
26 68
                      wt_is_1 = c("subtract", "divide", "force", "no"),
27 69
                      ## wt_is_1 = TRUE, ## wt has fitness 1
... ...
@@ -30,7 +72,10 @@ rfitness <- function(g, c= 0.5,
30 72
                                  ## wt_is_1 = divide
31 73
                      min_accessible_genotypes = NULL,
32 74
                      accessible_th = 0,
33
-                     truncate_at_0 = TRUE) {
75
+                     truncate_at_0 = TRUE,
76
+                     K = 1,
77
+                     r = TRUE,
78
+                     model = c("RMF", "NK")) {
34 79
     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene
35 80
     ## and Crona, 2014. And this allows moving from HoC to purely additive
36 81
     ## changing c and sd.
... ...
@@ -38,29 +83,77 @@ rfitness <- function(g, c= 0.5,
38 83
     ## FIXME future: do this with order too?
39 84
     ##    - do not generate a matrix of genotypes but a matrix of ordered genot.
40 85
     wt_is_1 = match.arg(wt_is_1)
86
+    model = match.arg(model)
41 87
     if(is_null_na(g)) stop("Number of genes argument (g) is null or NA")
42 88
     m <- generate_matrix_genotypes(g)
43 89
     done <- FALSE
44 90
     ## attempts <- 0 ## for debugging/tracking purposes
45 91
     while(!done) {
46
-        ## attempts <- attempts + 1
47
-        f_r <- rnorm(nrow(m), mean = mu, sd = sd)
48
-        if(inherits(reference, "character") && length(reference) == 1) {
49
-            if(reference == "random") {
50
-                referenceI <- m[sample(nrow(m), 1), ]
51
-            } else if(reference == "max") {
52
-                referenceI <- rep(1, g)
53
-            } else if(reference == "random2") {
54
-                referenceI <- create_eq_ref(g)
55
-            }
56
-        } else {
57
-            referenceI <- reference
92
+        if(model == "RMF") {
93
+            ## attempts <- attempts + 1
94
+            f_r <- rnorm(nrow(m), mean = mu, sd = sd)
95
+            if(inherits(reference, "character") && length(reference) == 1) {
96
+                if(reference == "random") {
97
+                    referenceI <- m[sample(nrow(m), 1), ]
98
+                } else if(reference == "max") {
99
+                    referenceI <- rep(1, g)
100
+                } else if(reference == "random2") {
101
+                    referenceI <- create_eq_ref(g)
102
+                }
103
+            } else {
104
+                referenceI <- reference
58 105
             }
59
-        d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))
60
-        f_det <- -c * d_reference
61
-        ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona
62
-        fi <- f_r + f_det
63
-        
106
+            d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))
107
+            f_det <- -c * d_reference
108
+            ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona
109
+            fi <- f_r + f_det
110
+        } else if(model == "NK") {
111
+            if(K >= g) stop("It makes no sense to have K >= g")
112
+            argsnk <- paste0("-K ", K,
113
+                             ifelse(r, " -r ", " "),
114
+                             g, " 2")
115
+            fl1 <- system2(fl_generate_binary(), args = argsnk, stdout = TRUE)[-1]
116
+            fl1 <- matrix(
117
+                as.numeric(unlist(strsplit(paste(fl1, collapse = " "), " "))),
118
+                ncol = g + 1, byrow = TRUE)
119
+            m1 <- fl1[, 1:g]
120
+            fi <- fl1[, g + 1]
121
+
122
+            ## For scaling, etc, all that matters, if anything, is the wildtype
123
+
124
+            ## We could order by doing this
125
+            ## But I am not 100% sure this will always be the same as
126
+            ## generate_matrix_genotypes
127
+            ## oo <- do.call(order,
128
+            ##               c(list(muts),
129
+            ##                 as.list(data.frame(m1[, rev(1:ncol(m1))]))))
130
+            ## m2 <- m1[oo, ]
131
+            ## Or create an id and order by it?
132
+
133
+            ## When we get to 20 genes, this is slow, about 18 seconds each
134
+            ## apply. Matching is fast (< 0.5 seconds)
135
+            gtstring <- apply(m, 1, function(x) paste0(x, collapse = ""))
136
+            gtstring2 <- apply(m1, 1, function(x) paste0(x, collapse = ""))
137
+            oo <- match(gtstring, gtstring2)
138
+            fi <- fi[oo]
139
+            ## make sure no left overs
140
+            rm(gtstring, gtstring2, oo, fl1, m1)
141
+
142
+            ## Had we not ordered, do this!!!
143
+            ## Which one is WT?
144
+            ## muts <- rowSums(m1)
145
+            ## w_wt <- which(muts == 0)
146
+            ## if(w_wt != 1) {
147
+            ##     f_a <- fi[1]
148
+            ##     fi[1] <- fi[w_wt]
149
+            ##     fi[w_wt] <- f_a
150
+            ##     rm(f_a)
151
+            ## }
152
+            ## m[] <- m1
153
+            ## rm(m1)
154
+            ## rm(fl1)
155
+        }
156
+
64 157
         if(!is.null(scale)) {
65 158
             fi <- (fi - min(fi))/(max(fi) - min(fi))
66 159
             fi <- scale[1] + fi * (scale[2] - scale[1])
... ...
@@ -119,16 +212,114 @@ rfitness <- function(g, c= 0.5,
119 212
 }
120 213
 
121 214
 
122
-create_eq_ref <- function(g) {
123
-    ## "random" gives more prob. to genotypes with
124
-    ## number of mutated genes close to g/2.
125
-    ## This gives equal prob to having the reference
126
-    ## be of any of the possible number of mutated genes.
127
-    nm <- sample(g, 1)
128
-    ref <- c(rep(1, nm), rep(0, g - nm))
129
-    sample(ref)
130
-}
131 215
 
132 216
 
133 217
 
134 218
 
219
+
220
+
221
+
222
+
223
+## rfitness <- function(g, c= 0.5,
224
+##                      sd = 1,
225
+##                      mu = 1,
226
+##                      reference = "random", ## "random", "max", or the vector,
227
+##                                      ## e.g., rep(1, g). If random, a
228
+##                                      ## random genotype is chosen as
229
+##                                      ## reference. If "max" this is rep(1, g)
230
+##                      scale = NULL, ## a two-element vector: min and max
231
+##                      wt_is_1 = c("subtract", "divide", "force", "no"),
232
+##                      ## wt_is_1 = TRUE, ## wt has fitness 1
233
+##                      log = FALSE, ## log only makes sense if all values >
234
+##                                  ## 0. scale with min > 0, and/or set
235
+##                                  ## wt_is_1 = divide
236
+##                      min_accessible_genotypes = NULL,
237
+##                      accessible_th = 0,
238
+##                      truncate_at_0 = TRUE) {
239
+##     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene
240
+##     ## and Crona, 2014. And this allows moving from HoC to purely additive
241
+##     ## changing c and sd.
242
+
243
+##     ## FIXME future: do this with order too?
244
+##     ##    - do not generate a matrix of genotypes but a matrix of ordered genot.
245
+##     wt_is_1 = match.arg(wt_is_1)
246
+##     if(is_null_na(g)) stop("Number of genes argument (g) is null or NA")
247
+##     m <- generate_matrix_genotypes(g)
248
+##     done <- FALSE
249
+##     ## attempts <- 0 ## for debugging/tracking purposes
250
+##     while(!done) {
251
+##         ## attempts <- attempts + 1
252
+##         f_r <- rnorm(nrow(m), mean = mu, sd = sd)
253
+##         if(inherits(reference, "character") && length(reference) == 1) {
254
+##             if(reference == "random") {
255
+##                 referenceI <- m[sample(nrow(m), 1), ]
256
+##             } else if(reference == "max") {
257
+##                 referenceI <- rep(1, g)
258
+##             } else if(reference == "random2") {
259
+##                 referenceI <- create_eq_ref(g)
260
+##             }
261
+##         } else {
262
+##             referenceI <- reference
263
+##             }
264
+##         d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))
265
+##         f_det <- -c * d_reference
266
+##         ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona
267
+##         fi <- f_r + f_det
268
+
269
+##         if(!is.null(scale)) {
270
+##             fi <- (fi - min(fi))/(max(fi) - min(fi))
271
+##             fi <- scale[1] + fi * (scale[2] - scale[1])
272
+##         }
273
+##         if(wt_is_1 == "divide") {
274
+##             ## We need to shift to avoid ratios involving negative numbers and
275
+##             ## we need to avoid having any fitness at 0, specially the wt.  If
276
+##             ## any negative value, add the min, and shift by the magnitude of
277
+##             ## the min to avoid any at 0.
278
+
279
+##             ## If you use scale and wt_is_1, this will move the scale. It is
280
+##             ## not possible to obtain a linear transformation that will keep
281
+##             ## the min and max of the scale, and wt at 1.
282
+##             min_fi <- min(fi)
283
+##             if(min_fi < 0)
284
+##                 fi <- fi + 2 * abs(min(fi))
285
+##             fi <- fi/fi[1]
286
+##         } else if (wt_is_1 == "subtract") {
287
+##             fi <- fi - fi[1] + 1.0
288
+##         } else if(wt_is_1 == "force") {
289
+##             fi[1] <- 1.0
290
+##             if(!is.null(scale)) {
291
+##                 if( (1 > scale[2]) || (1 < scale[1]))
292
+##                     warning("Using wt_is_1 = force and scale, but scale does ",
293
+##                             "not include 1")
294
+##             }
295
+##         }
296
+##         if(truncate_at_0) {
297
+##             fi[fi <= 0] <- 1e-9
298
+##         }
299
+##         if(log) {
300
+##             fi <- log(fi/fi[1]) + 1
301
+##         }
302
+##         m <- cbind(m, Fitness = fi)
303
+##         if(!is_null_na(min_accessible_genotypes)) {
304
+##             ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)
305
+##             ## Recall accessibleGenotypes includes the wt, if accessible.
306
+##             num_accessible_genotypes <- length(wrap_accessibleGenotypes(m, accessible_th)) - 1
307
+##             ## message("\n     num_accessible_genotypes = ", num_accessible_genotypes, "\n")
308
+##             if(num_accessible_genotypes >= min_accessible_genotypes) {
309
+##                 done <- TRUE
310
+##                 attributes(m) <- c(attributes(m),
311
+##                                    accessible_genotypes = num_accessible_genotypes,
312
+##                                    accessible_th = accessible_th)
313
+##             } else {
314
+##                 ## Cannot start again with a fitness column
315
+##                 m <- m[, -ncol(m), drop = FALSE]
316
+##             }
317
+##         } else {
318
+##             done <- TRUE
319
+##         }
320
+##     }
321
+##     ## message("\n number of attempts = ", attempts, "\n")
322
+##     class(m) <- c(class(m), "genotype_fitness_matrix")
323
+##     return(m)
324
+## }
325
+
... ...
@@ -1,3 +1,9 @@
1
+Changes in version 2.15.1 (2019-06-06):
2
+	- Added MAGELLAN's sources and functionality from MAGELLAN.
3
+	
4
+Changes in version 2.15.0 (2019-06-06):
5
+	- Bumped version to match current Biocdevel.
6
+	
1 7
 Changes in version 2.13.2 (2019-03-18):
2 8
 	- changes in behavior of sample
3 9
 	  (see NEWS and https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17494) were
... ...
@@ -5,7 +5,7 @@
5 5
 \title{Generate random fitness.}
6 6
 
7 7
 \description{ Generate random fitness landscapes under a House of Cards,
8
-  Rough Mount Fuji, or additive model.  }
8
+  Rough Mount Fuji, additive model, and Kauffman's NK model.  }
9 9
 
10 10
 
11 11
 \usage{
... ...
@@ -13,7 +13,8 @@
13 13
 rfitness(g, c = 0.5, sd = 1, mu = 1, reference = "random", scale = NULL,
14 14
          wt_is_1 = c("subtract", "divide", "force", "no"),
15 15
          log = FALSE, min_accessible_genotypes = NULL,
16
-         accessible_th = 0, truncate_at_0 = TRUE)
16
+         accessible_th = 0, truncate_at_0 = TRUE,
17
+         K = 1, r = TRUE, model = c("RMF", "NK"))
17 18
 }
18 19
 
19 20
 
... ...
@@ -51,7 +52,7 @@ mean \code{mu} and standard deviation \code{sd}).}
51 52
   two-element vector, fitness is re-scaled between \code{scale[1]} (the
52 53
   minimum) and \code{scale[2]} (the maximum).}
53 54
 
54
-\item{wt_is_1}{If "divide" (the default) the fitness of all genotypes is
55
+\item{wt_is_1}{If "divide" the fitness of all genotypes is
55 56
   divided by the fitness of the wildtype (after possibly adding a value
56 57
   to ensure no negative fitness) so that the wildtype (the genotype with
57 58
   no mutations) has fitness 1. This is a case of scaling, and it is
... ...
@@ -60,7 +61,7 @@ mean \code{mu} and standard deviation \code{sd}).}
60 61
   likely that the final fitness will not respect the limits in
61 62
   \code{scale}.
62 63
 
63
-  If "subtract" we shift all the fitness values (subtracting fitness of
64
+  If "subtract" (the default) we shift all the fitness values (subtracting fitness of
64 65
   the wildtype and adding 1) so that the wildtype ends up with a fitness
65 66
   of 1. This is also applied after \code{scale}, so if you specify both
66 67
   "wt_is_1 = 'subtract'" and use an argument for \code{scale} it is most
... ...
@@ -114,13 +115,23 @@ mean \code{mu} and standard deviation \code{sd}).}
114 115
   log) with values <=0. Or we might have trouble if we want to log the
115 116
   fitness.}
116 117
 
118
+\item{K}{K for NK model; K is the number of loci with which each locus
119
+  interacts, and the larger the K the larger the ruggedness of the
120
+  landscape.}
121
+
122
+\item{r}{For the NK model, whether interacting loci are chosen at random
123
+  (\code{r = TRUE}) or are neighbors (\code{r = FALSE}).}
124
+
125
+\item{model}{One of "RMF" (default), for Rough Mount Fuji, or "NK", for
126
+  Kauffman's NK model.}
117 127
 } 
118 128
 
119 129
 
120 130
 \details{
121 131
 
122
-  The model used here follows the Rough Mount Fuji model in Szendro et
123
-  al., 2013 or Franke et al., 2011. Fitness is given as
132
+  When using \code{model = "RMF"}, the model used here follows
133
+  the Rough Mount Fuji model in Szendro et al., 2013 or Franke et al.,
134
+  2011. Fitness is given as
124 135
 
125 136
   \deqn{f(i) = -c d(i, reference) + x_i}
126 137
 
... ...
@@ -144,6 +155,15 @@ mean \code{mu} and standard deviation \code{sd}).}
144 155
   is different from zero. In this case, with \code{c} large, the range
145 156
   of the data can be large, specially if \code{g} (the number of genes)
146 157
   is large.
158
+
159
+
160
+  When using \code{model = "NK"}, the model used is Kauffman's NK model
161
+  (see details in Ferretti et al., or Brouillet et al., below), as
162
+  implemented in MAGELLAN
163
+  (\url{http://wwwabi.snv.jussieu.fr/public/Magellan/}). This fitness
164
+  landscape is generated by directly calling the \code{fl_generate}
165
+  function of MAGELLAN. Fitness is drawn from a uniform (0, 1)
166
+  distribution.
147 167
   
148 168
 } 
149 169
 
... ...
@@ -159,7 +179,12 @@ mean \code{mu} and standard deviation \code{sd}).}
159 179
   \code{accessible_th} that show the number of accessible
160 180
   genotypes under the specified  threshold.
161 181
 }
162
-  
182
+
183
+
184
+\note{MAGELLAN uses its own random number generating functions; using
185
+  \code{set.seed} does not allow to obtain the same fitness landscape
186
+  repeatedly.}
187
+
163 188
 \references{
164 189
 
165 190
   Szendro I.~G. et al. (2013). Quantitative analyses of empirical
... ...
@@ -169,9 +194,23 @@ fitness landscapes. \emph{Journal of Statistical Mehcanics: Theory and
169 194
 Franke, J. et al. (2011). Evolutionary accessibility of mutational
170 195
 pathways. \emph{PLoS Computational Biology\/}, \bold{7}(8), 1--9.
171 196
 
197
+Brouillet, S. et al. (2015). MAGELLAN: a tool to explore small fitness
198
+landscapes. \emph{bioRxiv},
199
+\bold{31583}. \url{http://doi.org/10.1101/031583}
200
+
201
+Ferretti, L., Schmiegelt, B., Weinreich, D., Yamauchi, A., Kobayashi,
202
+Y., Tajima, F., & Achaz, G. (2016). Measuring epistasis in fitness
203
+landscapes: The correlation of fitness effects of mutations. \emph{Journal of
204
+Theoretical Biology\/}, \bold{396}, 132--143. \url{https://doi.org/10.1016/j.jtbi.2016.01.037}
205
+
206
+MAGELLAN web site: \url{http://wwwabi.snv.jussieu.fr/public/Magellan/}
207
+
172 208
 }
173 209
 
174
-\author{ Ramon Diaz-Uriarte
210
+\author{ Ramon Diaz-Uriarte for the RMF and general wrapping
211
+code. S. Brouillet, G. Achaz, S. Matuszewski, H. Annoni, and L. Ferreti
212
+for the MAGELLAN code.
213
+
175 214
 }
176 215
 
177 216
 \seealso{
... ...
@@ -192,6 +231,11 @@ r1 <- rfitness(4)
192 231
 plot(r1)
193 232
 oncoSimulIndiv(allFitnessEffects(genotFitness = r1))
194 233
 
234
+
235
+## NK model
236
+rnk <- rfitness(5, K = 3, model = "NK")
237
+plot(rnk)
238
+oncoSimulIndiv(allFitnessEffects(genotFitness = rnk))
195 239
 }
196 240
 
197 241
 \keyword{ datagen }
... ...
@@ -1,16 +1,24 @@
1 1
 \name{to_Magellan}
2 2
 \alias{to_Magellan}
3
+\alias{Magellan_stats}
3 4
 
4
-\title{ Create output for MAGELLAN.  }
5
+\title{ Create output for MAGELLAN and obtain MAGELLAN statistics.  }
5 6
 
6
-\description{ Create a fitness landscape in a format that is understood
7
-  by MAGELLAN \url{http://wwwabi.snv.jussieu.fr/public/Magellan/}. 
7
+\description{ Export a fitness landscape in a format that is understood
8
+  by MAGELLAN \url{http://wwwabi.snv.jussieu.fr/public/Magellan/} and
9
+  obtain fitness landscape statistics from MAGELLAN. 
8 10
 
9 11
 }
10 12
 
11 13
 \usage{
12 14
 to_Magellan(x, file,
13 15
             max_num_genotypes = 2000)
16
+
17
+Magellan_stats(x, max_num_genotypes = 2000,
18
+               verbose = FALSE,
19
+               use_log = TRUE,
20
+               short = TRUE,
21
+               replace_missing = FALSE)
14 22
 }
15 23
 
16 24
 \arguments{
... ...
@@ -49,15 +57,31 @@ to_Magellan(x, file,
49 57
   created using \code{\link{tempfile}}.}
50 58
 
51 59
   \item{max_num_genotypes}{Maximum allowed number of genotypes. For some
52
-  types of input, we make a call to \code{\link{evalAllGenotypes}}, and
53
-  use this as the maximum.}
54
-
60
+    types of input, we make a call to \code{\link{evalAllGenotypes}}, and
61
+    use this as the maximum.}
62
+  
63
+  \item{verbose}{If TRUE provide additional information about names of
64
+    intermediate files.}
65
+  \item{use_log}{Use log fitness when computing statistics.}
66
+  \item{short}{Give short output when computing statistics.}
67
+
68
+  \item{replace_missing}{From MAGELLAN's \code{fl_statistics}: replace
69
+  missing fitness values with 0 (otherwise check that all values are
70
+  specified).}
55 71
 }
56 72
 
57 73
 \value{
58 74
 
59
-  A file is written to disk. You can then plot and/or show summary
75
+  \code{to_Magellan}: A file is written to disk. You can then plot and/or show summary
60 76
   statistics using MAGELLAN.
77
+
78
+  \code{Magellan_stats}: MAGELLAN's statistics for fitness
79
+  landscapes. If you use \code{short = TRUE} a vector of statistics is
80
+  returned. If \code{short = FALSE}, MAGELLAN returns a file with
81
+  detailed statistics that cannot be turned into a simple vector of
82
+  statistics. The returned object uses \code{readLines} and, as a
83
+  message, you are also shown the path of the file, in case you want to
84
+  process it yourself.
61 85
 }
62 86
 
63 87
 \note{
... ...
@@ -101,6 +125,21 @@ cs <-  data.frame(parent = c(rep("Root", 3), "a", "d", "c"),
101 125
 
102 126
 to_Magellan(allFitnessEffects(cs), NULL)
103 127
 
128
+## Default, short output
129
+Magellan_stats(allFitnessEffects(cs))
130
+
131
+## Long output; since it is a > 200 lines file,
132
+## place in an object. Name of output file is given as message
133
+statslong <- Magellan_stats(allFitnessEffects(cs), short = FALSE)
134
+
135
+
136
+## Default, short output of two NK fitness landscapes
137
+rnk1 <- rfitness(6, K = 1, model = "NK")
138
+Magellan_stats(rnk1)
139
+
140
+rnk2 <- rfitness(6, K = 4, model = "NK")
141
+Magellan_stats(rnk2)
142
+
104 143
 }
105 144
 
106 145
 \keyword{ manip }
107 146
new file mode 100644
... ...
@@ -0,0 +1,956 @@
1
+/*
2
+	All functions needed for Linear Algrebra
3
+	This is especially used for Least Square Estimates
4
+*/
5
+
6
+/*
7
+	All functions needed for Linear Algrebra
8
+	This is especially used for Least Square Estimates
9
+	
10
+	
11
+	date: Nov 2103
12
+	author: gachaz
13
+*/
14
+
15
+
16
+#include <stdio.h>
17
+#include <stdlib.h>
18
+#include <math.h>
19
+
20
+#include "LinearAlgebra.h"
21
+
22
+
23
+/*
24
+	Memory functions
25
+*/
26
+
27
+struct matrix MemMat( int m, int n ){
28
+
29
+	int i;
30
+	struct matrix M;
31
+	
32
+	M.r=m;
33
+	M.c=n;
34
+	
35
+	M.val = (double **) malloc( (size_t) M.r * sizeof(double *) );
36
+	if( ! M.val ) fprintf(stderr, "MemMat: cannot allocate matrix, bye\n"), exit(1);
37
+	
38
+	for(i=0; i< M.r ; i++){
39
+		M.val[i] = (double *) malloc( (size_t) M.c * sizeof(double) );
40
+		if(!M.val[i]) fprintf(stderr, "MemMat: cannot allocate M[%d], bye\n",i), exit(1);
41
+	}
42
+
43
+	return M;
44
+}
45
+
46
+void IdentMat(struct matrix *M)
47
+{
48
+	
49
+	int i, j;
50
+	
51
+	for (i = 0; i < M->r; i++)
52
+	{
53
+		for (j = 0; j < M->c; j++)
54
+		{
55
+			if (i==j)
56
+			{
57
+				M->val[i][i] = 1.;
58
+			}
59
+			else
60
+			{
61
+				M->val[i][j] = 0.;
62
+			}
63
+		}
64
+	}
65
+}
66
+
67
+
68
+void FreeMat( struct matrix *M ){
69
+
70
+	int i;
71
+	
72
+	if(M->val == NULL)
73
+		return;
74
+	
75
+	for(i=0;i<M->r; i++)
76
+		free( M->val[i] );
77
+
78
+	free(M->val);
79
+	
80
+	M->val =NULL;
81
+}
82
+
83
+
84
+/*
85
+	Output a matrix
86
+*/
87
+void PrintMat( struct matrix *M, char *name ){
88
+
89
+	int i,j;
90
+
91
+	printf("%s [%d,%d]\n", name, M->r, M->c  );
92
+
93
+	for(i=0;i<M->r; i++, printf("\n") )
94
+		for(j=0;j<M->c;j++)
95
+			printf("%8.3f", M->val[i][j]);
96
+
97
+	printf("//\n");
98
+}
99
+
100
+
101
+
102
+struct matrix MatrixTranspose( struct matrix *M, char opt_free ){
103
+
104
+	int i,j;
105
+	struct matrix T = MemMat( M->c, M->r );	
106
+	
107
+	for(i=0 ; i<M->r ; i++)
108
+		for(j=0 ; j<M->c ; j++)
109
+			T.val[j][i] = M->val[i][j];
110
+	
111
+	
112
+	if(opt_free)
113
+		FreeMat(M);
114
+	
115
+	return T;
116
+	
117
+}
118
+
119
+/*
120
+	From a vector, create a Diag Matrix
121
+*/
122
+struct matrix CreateMatrixDiag( struct matrix *V, char opt_free )
123
+{
124
+	int i,j;
125
+	
126
+	struct matrix D = MemMat( V->r,V->r );
127
+	
128
+	for(i=0; i<D.r; i++ )
129
+	{
130
+		for(j=0; j<D.c; j++ )
131
+		{
132
+			D.val[i][j]=0;
133
+			
134
+			if( i == j )
135
+			{
136
+				D.val[i][j] = V->val[i][0];
137
+			}
138
+		}
139
+		
140
+	}
141
+	
142
+	if(opt_free)
143
+	{
144
+		FreeMat(V);
145
+	}
146
+	
147
+	return D;
148
+}
149
+
150
+
151
+struct matrix CreateMatrixDiagMatrix( struct matrix *M1, char opt_free )
152
+{
153
+	int i,j;
154
+	
155
+	struct matrix D = MemMat( M1->r,M1->c );
156
+	
157
+	for(i=0; i<D.r; i++ )
158
+	{
159
+		for(j=0; j<D.c; j++ )
160
+		{
161
+			D.val[i][j] = 0;
162
+			if( i == j )
163
+			{
164
+				D.val[i][j] = M1->val[i][i];
165
+			}
166
+			
167
+		}
168
+	}
169
+	
170
+	if(opt_free)
171
+	{
172
+		FreeMat(M1);
173
+	}
174
+	
175
+		return (D);
176
+}
177
+
178
+
179
+
180
+struct matrix Invert_MatrixDiag(struct matrix *D, char opt_free )
181
+{
182
+	int i,j;
183
+	struct matrix M;
184
+	
185
+	M = MemMat( D->r, D->c );
186
+		
187
+	for(i=0; i<M.r; i++ )
188
+		for(j=0; j<M.c; j++ )
189
+		{
190
+			M.val[i][j]=0;
191
+			
192
+			if( i == j )
193
+				M.val[i][j] = 1.0/D->val[i][j];
194
+			
195
+		}
196
+	
197
+	if(opt_free)
198
+		FreeMat(D);
199
+	
200
+
201
+	return M;
202
+		
203
+}
204
+
205
+
206
+struct matrix MatrixProduct(struct matrix *M_left, struct matrix *M_right, char opt_free)
207
+{
208
+
209
+	int i,j,k;
210
+	struct matrix P;
211
+	
212
+	if( M_left->c != M_right->r )
213
+		fprintf(stderr, "MatrixProduct: The matrices cannot be multiplied together. Incompatible dimensions, bye\n"),exit(1);
214
+	
215
+	
216
+	P=MemMat( M_left->r, M_right->c );
217
+	
218
+	for(i=0; i<P.r; i++ )
219
+	{
220
+		for(j=0; j<P.c; j++ )
221
+		{
222
+			P.val[i][j]=0;
223
+			
224
+			for(k=0; k<M_left->c; k++ )
225
+			{
226
+				P.val[i][j] += M_left->val[i][k] * M_right->val[k][j];
227
+			}
228
+		}
229
+	}
230
+	
231
+	if(opt_free)
232
+	{
233
+		FreeMat(M_right);
234
+		FreeMat(M_left);
235
+	}
236
+	
237
+	return P;
238
+			
239
+}
240
+
241
+
242
+void MatrixProductRef(struct matrix *M_left, struct matrix *M_right, struct matrix * Prod, char opt_free)
243
+{
244
+	
245
+	int i,j,k;
246
+	
247
+	if( M_left->c != M_right->r )
248
+		fprintf(stderr, "MatrixProduct: The matrices cannot be multiplied together. Incompatible dimensions, bye\n"),exit(1);
249
+	
250
+		
251
+		for(i=0; i<Prod->r; i++ )
252
+		{
253
+			for(j=0; j<Prod->c; j++ )
254
+			{
255
+				Prod->val[i][j]=0;
256
+				
257
+				for(k=0; k<M_left->c; k++ )
258
+				{
259
+					Prod->val[i][j] += M_left->val[i][k] * M_right->val[k][j];
260
+				}
261
+			}
262
+		}
263
+	
264
+	if(opt_free)
265
+	{
266
+		FreeMat(M_right);
267
+		FreeMat(M_left);
268
+	}
269
+	
270
+	
271
+}
272
+
273
+struct matrix MatrixSum(struct matrix *M1, struct matrix *M2, char opt_free)
274
+{
275
+
276
+	int i,j;
277
+	struct matrix S;
278
+	
279
+	
280
+	if( M1->r != M2->r || M1->c != M2->c )
281
+		fprintf(stderr, "MatrixSum: The matrices cannot be summed together. Incompatible dimensions, bye\n"),exit(1);
282
+	
283
+
284
+	S = MemMat( M1->r, M1->c );
285
+		
286
+	for(i=0; i<S.r; i++ )
287
+	{
288
+		for(j=0; j<S.c; j++ )
289
+		{
290
+			S.val[i][j] = M1->val[i][j] + M2->val[i][j];
291
+		}
292
+	}
293
+	
294
+	if(opt_free)
295
+	{
296
+		FreeMat(M1);
297
+		FreeMat(M2);
298
+	}
299
+	
300
+	return S;
301
+}
302
+
303
+void MatrixSumRef(struct matrix *M1, struct matrix *M2, struct matrix *Sum, char opt_free)
304
+{
305
+	
306
+	int i,j;
307
+	
308
+	if( M1->r != M2->r || M1->c != M2->c )
309
+		fprintf(stderr, "MatrixSum: The matrices cannot be summed together. Incompatible dimensions, bye\n"),exit(1);
310
+		
311
+	
312
+		
313
+		for(i=0; i<Sum->r; i++ )
314
+		{
315
+			for(j=0; j<Sum->c; j++ )
316
+			{
317
+				Sum->val[i][j] = M1->val[i][j] + M2->val[i][j];
318
+			}
319
+		}
320
+				
321
+				if(opt_free)
322
+				{
323
+					FreeMat(M1);
324
+					FreeMat(M2);
325
+				}
326
+}
327
+
328
+
329
+
330
+
331
+struct matrix MatrixDiff(struct matrix *M1, struct matrix *M2, char opt_free)
332
+{
333
+	
334
+	int i,j;
335
+	struct matrix S;
336
+	
337
+	
338
+	if( M1->r != M2->r || M1->c != M2->c )
339
+		fprintf(stderr, "MatrixSum: The matrices cannot be substracted. Incompatible dimensions, bye\n"),exit(1);
340
+		
341
+		
342
+		S = MemMat( M1->r, M1->c );
343
+		
344
+		for(i=0; i<S.r; i++ )
345
+		{
346
+			for(j=0; j<S.c; j++ )
347
+			{
348
+				S.val[i][j] = M1->val[i][j] - M2->val[i][j];
349
+			}
350
+		}
351
+				
352
+				if(opt_free)
353
+				{
354
+					FreeMat(M1);
355
+					FreeMat(M2);
356
+				}
357
+	
358
+	return S;
359
+}
360
+
361
+void MatrixDiffRef(struct matrix *M1, struct matrix *M2, struct matrix * Diff, char opt_free)
362
+{
363
+	
364
+	int i,j;
365
+
366
+	if( M1->r != M2->r || M1->c != M2->c )
367
+		fprintf(stderr, "MatrixSum: The matrices cannot be substracted. Incompatible dimensions, bye\n"),exit(1);
368
+	
369
+	
370
+	for(i=0; i<Diff->r; i++ )
371
+	{
372
+		for(j=0; j<Diff->c; j++ )
373
+		{
374
+			Diff->val[i][j] = M1->val[i][j] - M2->val[i][j];
375
+		}
376
+	}
377
+	
378
+	
379
+	if(opt_free)
380
+	{
381
+		FreeMat(M1);
382
+		FreeMat(M2);
383
+	}
384
+	
385
+}
386
+
387
+
388
+struct matrix MatrixHadarmardProduct(struct matrix *M1, struct matrix *M2, char opt_free)
389
+{
390
+	struct matrix HadarmardProduct;
391
+	int i, j;
392
+
393
+	if( M1->r != M2->r || M1->c != M2->c )
394
+	{
395
+		fprintf(stderr, "MatrixHadarmardProduct: The matrices cannot be multiplied together. Incompatible dimensions, bye\n"),exit(1);
396
+	}
397
+
398
+	HadarmardProduct = MemMat(M1->r, M1->c);
399
+	
400
+	for(i = 0; i < M1->r; i++)
401
+	{
402
+		for(j = 0; j < M1->c; j++)
403
+		{
404
+			HadarmardProduct.val[i][j] = M1->val[i][j] * M2->val[i][j];
405
+		}
406
+	}
407
+	
408
+	if(opt_free)
409
+	{
410
+		FreeMat(M1);
411
+		FreeMat(M2);
412
+	}
413
+	
414
+	return (HadarmardProduct);
415
+}
416
+
417
+
418
+void MatrixHadarmardProductRef(struct matrix *M1, struct matrix *M2, struct matrix * HadarmardProduct, char opt_free)
419
+{
420
+	int i, j;
421
+	
422
+	if( M1->r != M2->r || M1->c != M2->c )
423
+	{
424
+		fprintf(stderr, "MatrixHadarmardProduct: The matrices cannot be multiplied together. Incompatible dimensions, bye\n"),exit(1);
425
+	}
426
+	
427
+	for(i = 0; i < M1->r; i++)
428
+	{
429
+		for(j = 0; j < M1->c; j++)
430
+		{
431
+			HadarmardProduct->val[i][j] = M1->val[i][j] * M2->val[i][j];
432
+		}
433
+	}
434
+	
435
+	if(opt_free)
436
+	{
437
+		FreeMat(M1);
438
+		FreeMat(M2);
439
+	}
440
+	
441
+}
442
+
443
+
444
+struct matrix scaleMatrix(struct matrix *M1, float scale, char opt_free)
445
+{
446
+	struct matrix M2;
447
+	int i, j;
448
+	
449
+	M2 = MemMat(M1->r, M1->c);
450
+	for(i = 0; i < M1->r; i++)
451
+	{
452
+		for(j = 0; j < M1->c; j++)
453
+		{
454
+			M2.val[i][j] = scale * M1->val[i][j];
455
+		}
456
+	}
457
+	
458
+	if(opt_free)
459
+	{
460
+		FreeMat(M1);
461
+	}
462
+	
463
+	return (M2);
464
+}
465
+
466
+void scaleMatrixRef(struct matrix *M1, 	struct matrix * M2, float scale, char opt_free)
467
+{
468
+	int i, j;
469
+	
470
+	for(i = 0; i < M1->r; i++)
471
+	{
472
+		for(j = 0; j < M1->c; j++)
473
+		{
474
+			M2->val[i][j] = scale * M1->val[i][j];
475
+		}
476
+	}
477
+	
478
+	if(opt_free)
479
+	{
480
+		FreeMat(M1);
481
+	}
482
+	
483
+}
484
+
485
+/*
486
+	Do the SVD decomposition of any A matrix --dimension (m,n)--
487
+	A = U.w.Vt
488
+	U is orthogonal and dimension (m,n)
489
+	w is a diagonal matrix of size (n,n)
490
+	V is a square orthogonal matrix (n,n)
491
+	
492
+	The code has been adapted from numrec 92 -- reset mat from 0 to n-1 as in std C --
493
+	
494
+ */
495
+
496
+#define CONVERGENCE 100
497
+
498
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
499
+
500
+static int iminarg1,iminarg2;
501
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
502
+(iminarg1) : (iminarg2))
503
+
504
+
505
+static double dsqrarg;
506
+#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
507
+
508
+static double dmaxarg1,dmaxarg2;
509
+#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
510
+(dmaxarg1) : (dmaxarg2))
511
+
512
+
513
+/*
514
+	Apparently clever way to do (a^2 + b^2)^1/2
515
+ */
516
+static float pythag(float a, float b)
517
+{
518
+	float absa=fabs(a),
519
+	absb=fabs(b);
520
+	
521
+	if (absa > absb)
522
+		return absa*sqrt(1.0+DSQR(absb/absa));
523
+	else
524
+		return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb)));
525
+}
526
+
527
+static void svdcmp(double **A, int m, int n, double **u, double *w, double **v)
528
+{
529
+	int flag,i,its,j,jj,k,l,nm=0;
530
+	
531
+	double anorm,c,f,g,h,s,scale,x,y,z,
532
+	*rv1;
533
+	
534
+	if(u==NULL || w==NULL || v==NULL)
535
+		fprintf(stderr, "you should get the memory before using this function, bye\n"),exit(3);
536
+	
537
+	
538
+	for(i=0;i<m;i++)
539
+		for(j=0;j<n;j++)
540
+			u[i][j]= A[i][j];
541
+	
542
+	rv1=(double *)malloc( (size_t) n*sizeof(double) );
543
+	if(!rv1)fprintf(stderr, "cannot allocate rv1, bye");
544
+	
545
+	g=scale=anorm=0.0;
546
+	
547
+	for (i=0;i<n;i++) {
548
+		
549
+		l=i+1;
550
+		rv1[i]=scale*g;
551
+		g=s=scale=0.0;
552
+		
553
+		if (i < m)
554
+		{
555
+			
556
+			for (k=i;k<m;k++)
557
+				scale += fabs(u[k][i]);
558
+			
559
+			if (scale)
560
+			{
561
+				
562
+				for (k=i;k<m;k++)
563
+				{
564
+					u[k][i] /= scale;
565
+					s += u[k][i]*u[k][i];
566
+				}
567
+				
568
+				f=u[i][i];
569
+				g = -SIGN(sqrt(s),f);
570
+				h=f*g-s;
571
+				u[i][i]=f-g;
572
+				
573
+				for (j=l;j<n;j++)
574
+				{
575
+					for (s=0.0,k=i;k<m;k++)
576
+						s += u[k][i]*u[k][j];
577
+					f=s/h;
578
+					
579
+					for (k=i;k<m;k++)
580
+						u[k][j] += f*u[k][i];
581
+				}
582
+				
583
+				for (k=i;k<m;k++)
584
+					u[k][i] *= scale;
585
+			}
586
+		}
587
+		
588
+		
589
+		w[i]=scale *g;
590
+		g=s=scale=0.0;
591
+		
592
+		if (i < m && i != n-1)
593
+		{
594
+			for (k=l;k<n;k++)
595
+				scale += fabs(u[i][k]);
596
+			
597
+			if (scale)
598
+			{
599
+				
600
+				for (k=l;k<n;k++)
601
+				{
602
+					u[i][k] /= scale;
603
+					s += u[i][k]*u[i][k];
604
+				}
605
+				
606
+				f=u[i][l];
607
+				g = -SIGN(sqrt(s),f);
608
+				h=f*g-s;
609
+				u[i][l]=f-g;
610
+				
611
+				for (k=l;k<n;k++)
612
+					rv1[k]=u[i][k]/h;
613
+				
614
+				for (j=l;j<m;j++)
615
+				{
616
+					for (s=0.0,k=l;k<n;k++)
617
+						s += u[j][k]*u[i][k];
618
+					
619
+					for (k=l;k<n;k++)
620
+						u[j][k] += s*rv1[k];
621
+				}
622
+				
623
+				for (k=l;k<n;k++)
624
+					u[i][k] *= scale;
625
+			}
626
+		}
627
+		
628
+		anorm = DMAX( anorm , (fabs(w[i])+fabs(rv1[i])) );
629
+		
630
+	}
631
+	
632
+	
633
+	
634
+	for (i=n-1;i>=0;i--)
635
+	{
636
+		if (i < n-1){
637
+			
638
+			if (g) {
639
+				for (j=l;j<n;j++)
640
+					v[j][i]=(u[i][j]/u[i][l])/g;
641
+				
642
+				for (j=l;j<n;j++)
643
+				{
644
+					for (s=0.0,k=l;k<n;k++)
645
+						s += u[i][k]*v[k][j];
646
+					
647
+					for (k=l;k<n;k++)
648
+						v[k][j] += s*v[k][i];
649
+				}
650
+			}
651
+			
652
+			for (j=l;j<n;j++)
653
+				v[i][j]=v[j][i]=0.0;
654
+		}
655
+		v[i][i]=1.0;
656
+		g=rv1[i];
657
+		l=i;
658
+	}
659
+	
660
+	
661
+	for (i=IMIN(m-1,n-1);i>=0;i--) {
662
+		
663
+		l=i+1;
664
+		g=w[i];
665
+		
666
+		for (j=l;j<n;j++)
667
+			u[i][j]=0.0;
668
+		
669
+		if (g)
670
+		{
671
+			g=1.0/g;
672
+			
673
+			for (j=l;j<n;j++)
674
+			{
675
+				for (s=0.0,k=l;k<m;k++)
676
+					s += u[k][i]*u[k][j];
677
+				
678
+				f=(s/u[i][i])*g;
679
+				
680
+				for (k=i;k<m;k++)
681
+					u[k][j] += f*u[k][i];
682
+			}
683
+			
684
+			for (j=i;j<m;j++)
685
+				u[j][i] *= g;
686
+		}
687
+		else
688
+			for (j=i;j<m;j++)
689
+				u[j][i]=0.0;
690
+		++u[i][i];
691
+	}
692
+	
693
+	
694
+	
695
+	for (k=n-1;k>=0;k--) {
696
+		
697
+		for (its=1;its<=CONVERGENCE;its++) {
698
+			
699
+			flag=1;
700
+			
701
+			for (l=k;l>0;l--) {
702
+				
703
+				nm=l-1;
704
+				if ( (fabs(rv1[l])+anorm) == anorm ){
705
+					flag=0;
706
+					break;
707
+				}
708
+				
709
+				if ( (fabs(w[nm])+anorm) == anorm )
710
+					break;
711
+			}
712
+			
713
+			if (flag) {
714
+				
715
+				c=0.0;
716
+				s=1.0;
717
+				
718
+				for (i=l;i<k;i++) {
719
+					
720
+					f=s*rv1[i];
721
+					rv1[i]=c*rv1[i];
722
+					if ( (fabs(f)+anorm) == anorm ) break;
723
+					g=w[i];
724
+					h=pythag(f,g);
725
+					w[i]=h;
726
+					h=1.0/h;
727
+					c=g*h;
728
+					s = -f*h;
729
+					
730
+					for (j=0;j<m;j++) {
731
+						
732
+						y=u[j][nm];
733
+						z=u[j][i];
734
+						u[j][nm]=y*c+z*s;
735
+						u[j][i]=z*c-y*s;
736
+					}
737
+				}
738
+			}
739
+			
740
+			z=w[k];
741
+			
742
+			if (l == k) {
743
+				
744
+				if (z < 0.0)
745
+				{
746
+					w[k] = -z;
747
+					for (j=0;j<n;j++)
748
+						v[j][k] = -v[j][k];
749
+				}
750
+				break;
751
+			}
752
+			
753
+			if (its == CONVERGENCE)
754
+				fprintf(stderr, "no convergence in %d svdcmp iterations, bye\n", CONVERGENCE), exit(5);
755
+			
756
+			x=w[l];
757
+			nm=k-1;
758
+			y=w[nm];
759
+			g=rv1[nm];
760
+			h=rv1[k];
761
+			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
762
+			g=pythag(f,1.0);
763
+			f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
764
+			c=s=1.0;
765
+			
766
+			for (j=l;j<=nm;j++) {
767
+				
768
+				i=j+1;
769
+				g=rv1[i];
770
+				y=w[i];
771
+				h=s*g;
772
+				g=c*g;
773
+				z=pythag(f,h);
774
+				rv1[j]=z;
775
+				c=f/z;
776
+				s=h/z;
777
+				f=x*c+g*s;
778
+				g = g*c-x*s;
779
+				h=y*s;
780
+				y *= c;
781
+				
782
+				for (jj=0;jj<n;jj++) {
783
+					
784
+					x=v[jj][j];
785
+					z=v[jj][i];
786
+					v[jj][j]=x*c+z*s;
787
+					v[jj][i]=z*c-x*s;
788
+				}
789
+				
790
+				z=pythag(f,h);
791
+				w[j]=z;
792
+				
793
+				if (z)
794
+				{
795
+					z=1.0/z;
796
+					c=f*z;
797
+					s=h*z;
798
+				}
799
+				
800
+				f=c*g+s*y;
801
+				x=c*y-s*g;
802
+				
803
+				for (jj=0;jj<m;jj++)
804
+				{
805
+					y=u[jj][j];
806
+					z=u[jj][i];
807
+					u[jj][j]=y*c+z*s;
808
+					u[jj][i]=z*c-y*s;
809
+				}
810
+			}
811
+			rv1[l]=0.0;
812
+			rv1[k]=f;
813
+			w[k]=x;
814
+		}
815
+	}
816
+	
817
+	free(rv1);
818
+}
819
+
820
+
821
+
822
+/*
823
+	Calculating the inverse of a matrix M1 by SVD
824
+*/
825
+
826
+// :: FOR LARGE MATRICES THIS BECOMES NUMERICALLY UNSTABLE :: ANY BETTER OPTION AVAILABLE? :: ALTERNATIVELY ONE COULD JUST DO BRUTE FORCE MATRIX MULTIPLICATION :: WOULD BE OF COMPLEXITY O(x n^2) FOR A MATRIX OF DIMENSION n THAT IS RAISED TO THE POWER x ::
827
+
828
+void inverseMatrix(struct matrix * M1, struct matrix * inverseM)
829
+{
830
+	
831
+	struct matrix U,w,V;    /* misc , SVD decomposition (U,w,V) and matrix product (Tmp) */
832
+	
833
+	int i,j,k;                  /* counters */
834
+	
835
+	/*
836
+		get memory for the matrices
837
+	 */
838
+	U = MemMat( M1->r, M1->c);
839
+	w = MemMat( 1, M1->c );
840
+	V = MemMat( M1->c, M1->c );
841
+
842
+	
843
+	/*
844
+	 svd decomposition --adapted from numrec--
845
+	 svdcmp(double **A, int m, int n, double **u, double *w, double **v)
846
+	 */
847
+	
848
+	svdcmp(M1->val, M1->r, M1->c, U.val, w.val[0], V.val);
849
+
850
+	/*
851
+		Check the w's for 0 values
852
+	 */
853
+	for(k = 0; k < inverseM->r; k++)
854
+	{
855
+		if( w.val[0][k] < 1e-6 )
856
+		{
857
+			fprintf(stderr, "warning w[%d]=%e, need extra check\n", k, w.val[0][k]);
858
+		}
859
+	}
860
+			
861
+			
862
+	/*
863
+	 This is equivalent to
864
+	 V . Diag(w)^-1 . Ut
865
+	 but more efficient in computation time O(m n^2)
866
+	 */
867
+	for(i = 0; i < inverseM->r; i++ )
868
+	{
869
+		for(j = 0; j < inverseM->c; j++ )
870
+		{
871
+			inverseM->val[i][j]=0;
872
+			
873
+			for(k = 0; k < inverseM->r; k++ )
874
+			{
875
+				inverseM->val[i][j] += V.val[i][k]*U.val[j][k]/w.val[0][k];
876
+			}
877
+		}
878
+	}
879
+	
880
+	
881
+	FreeMat(&U);
882
+	FreeMat(&V);
883
+	FreeMat(&w);
884
+}
885
+
886
+
887
+
888
+struct matrix LeastSquare( struct matrix *X, struct matrix *Y ){
889
+
890
+
891
+	struct matrix beta;          /* this will contains the estimates you aim at */
892
+	
893
+	struct matrix Tmp, U,w,V;    /* misc , SVD decomposition (U,w,V) and matrix product (Tmp) */
894
+	
895
+	int i,j,k;                  /* counters */
896
+
897
+
898
+	/*
899
+		get memory for the matrices
900
+	*/
901
+	beta = MemMat( X->c, 1);
902
+	Tmp = MemMat( X->c, X->r);
903
+	U = MemMat( X->r, X->c);
904
+	w = MemMat( 1, X->c );
905
+	V = MemMat( X->c, X->c );
906
+	
907
+	
908
+	/*
909
+	 	svd decomposition --adapted from numrec--
910
+	 	svdcmp(double **A, int m, int n, double **u, double *w, double **v)
911
+	*/
912
+
913
+	svdcmp(X->val,  X->r, X->c, U.val, w.val[0], V.val);
914
+	
915
+	
916
+	/*
917
+		Check the w's for 0 values
918
+	*/
919
+	for(k=0; k<Tmp.r; k++)
920
+		if( w.val[0][k] < 1e-6 )
921
+			fprintf(stderr, "warning w[%d]=%e, need extra check\n", k, w.val[0][k]);
922
+	
923
+	
924
+	
925
+	/*
926
+		This is equivalent to 
927
+		V . Diag(w)^-1 . Vt
928
+		but more efficient in computation time O(m n^2)
929
+	*/
930
+	for(i=0; i<Tmp.r; i++ )
931
+		for(j=0; j<Tmp.c; j++ )
932
+		{
933
+			Tmp.val[i][j]=0;
934
+			
935
+			for(k=0; k<Tmp.r; k++ )
936
+				Tmp.val[i][j] += V.val[i][k]*U.val[j][k]/w.val[0][k];
937
+
938
+		}
939
+	
940
+	
941
+	
942
+	/*
943
+		From there, estimate the beta's
944
+	*/
945
+	beta = MatrixProduct( &Tmp , Y, 0 );
946
+	
947
+	
948
+	FreeMat(&U);
949
+	FreeMat(&V);
950
+	FreeMat(&w);
951
+	FreeMat(&Tmp);
952
+
953
+	return beta;
954
+
955
+}
956
+
0 957
new file mode 100644
... ...
@@ -0,0 +1,71 @@
1
+/*
2
+	All functions needed for Linear Algrebra
3
+	This is especially used for Least Square Estimates
4
+*/
5
+
6
+
7
+#ifndef __LINEAR_ALGEBRA__
8
+#define __LINEAR_ALGEBRA__
9
+
10
+struct matrix {
11
+	int r; 
12
+	int c;
13
+	double **val;
14
+};
15
+
16
+/*
17
+	Memory functions
18
+*/
19
+
20
+struct matrix MemMat( int m, int n );
21
+void IdentMat(struct matrix *M);
22
+void FreeMat( struct matrix *M );
23
+
24
+/*
25
+	Output a matrix
26
+*/
27
+void PrintMat( struct matrix *M, char *name );
28
+
29
+/*
30
+	Basic matrix operators
31
+	when opt_free is set to 1, it also free the entry matrix
32
+*/
33
+
34
+// :: WHAT ABOUT MEMORY HERE?:: WOULD IT BE BETTER TO MAKE ALL THESE FUNCTIONS VOID AND PASS THE MATRIX AS REFERENCE?:: //
35
+
36
+struct matrix MatrixProduct(struct matrix *M_left, struct matrix *M_right, char opt_free);
37
+void MatrixProductRef(struct matrix *M_left, struct matrix *M_right, struct matrix *Prod, char opt_free);
38
+
39
+struct matrix MatrixSum(struct matrix *M1, struct matrix *M2, char opt_free);
40
+void MatrixSumRef(struct matrix *M1, struct matrix *M2, struct matrix *Sum, char opt_free);
41
+
42
+struct matrix MatrixDiff(struct matrix *M1, struct matrix *M2, char opt_free);
43
+void MatrixDiffRef(struct matrix *M1, struct matrix *M2, struct matrix *Diff, char opt_free);
44
+
45
+
46
+struct matrix MatrixTranspose( struct matrix *M, char opt_free );
47
+struct matrix Invert_MatrixDiag(struct matrix *D, char opt_free );
48
+struct matrix MatrixHadarmardProduct(struct matrix *M1, struct matrix *M2, char opt_free);
49
+void MatrixHadarmardProductRef(struct matrix *M1, struct matrix *M2, struct matrix *HadarmardProduct, char opt_free);
50
+struct matrix scaleMatrix(struct matrix *M1, float scale, char opt_free);
51
+void scaleMatrixRef(struct matrix *M1, 	struct matrix *M2, float scale, char opt_free);
52
+void inverseMatrix(struct matrix *M1, struct matrix *inverseM);
53
+
54
+/*
55
+	From a vector, create a Diag Matrix
56
+*/
57
+struct matrix CreateMatrixDiag( struct matrix *Vector, char opt_free );
58
+
59
+/*
60
+	From a Matrix, create a Diag Matrix
61
+ */
62
+struct matrix CreateMatrixDiagMatrix( struct matrix *M1, char opt_free );
63
+
64
+/*
65
+	Return a vector [n,1] from a matrix X [m,n] and a vector Y[m,1]
66
+	m i the number of observation
67
+	n is the number of parameter to estimate
68
+*/
69
+struct matrix LeastSquare( struct matrix *X, struct matrix *Y );
70
+
71
+#endif
0 72
new file mode 100644
... ...
@@ -0,0 +1,169 @@
1
+
2
+# ATTENTION pour pouvoir sauvegarder en pdf installer ca.. sudo apt-get install librsvg2-bin
3
+
4
+
5
+BIN= fl_statistics fl_generate fl_draw  fl_genchains  fl_webdraw.cgi fl_webmodel.cgi fl_drawfile.cgi #test_gamma fl_gamma
6
+
7
+
8
+
9
+HOME_BIN=~/bin/
10
+
11
+# SGI:    RANLIB=touch  
12
+# LINUX:  RANLIB=ranlib 
13
+# MACOSX: RANLIB=ranlib 
14
+RANLIB=ranlib
15
+
16
+CC = gcc
17
+ALLCFLAGS= -O -Wall 
18
+
19
+SHELL=/bin/sh
20
+LIBS = liblandscape.a
21
+LIBS_DIR = -L.
22
+LFLAGS = -llandscape
23
+
24
+INCLUDE_DIR = -I.
25
+
26
+AR = ar
27
+
28
+SRC = gamma.c\
29
+      landscape.c \
30
+      random.c\
31
+      sort.c\
32
+      chain.c\
33
+      input.c\
34
+      common_drawings.c\
35
+      genotypes.c\
36
+      ordered_pairs.c\
37
+      calculus.c\
38
+      vector.c\
39
+      generalized_chain.c\
40
+      LinearAlgebra.c\
41
+      summary_statistics.c\
42
+      epistasis_type.c\
43
+      decomposition.c\
44
+      web/web_common.c\
45
+      models.c
46
+
47
+OBJECTS = $(SRC:.c=.o)
48
+
49
+
50
+all: liblandscape.a $(BIN)
51
+
52
+%.o: %.c
53
+	$(CC) $(ALLCFLAGS) $(INCLUDE_DIR)  -c -o $@ $<;
54
+
55
+#
56
+# Check these two and move them into fl_real
57
+# 
58
+
59
+#test_gamma: test_gamma.c $(LIBS)
60
+#	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
61
+
62
+fl_genchains: main_generalizedchains_test.c $(LIBS)
63
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm #-g
64
+
65
+#
66
+# Have a look and update if nescessary - very likelily deprecated
67
+#
68
+
69
+#fl_gamma: test_gamma.c $(LIBS)
70
+#	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm
71
+
72
+
73
+#
74
+# Command line exe
75
+# 
76
+
77
+fl_statistics: main_fl_statistics.c $(LIBS)
78
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
79
+
80
+fl_generate: main_fl_generate.c  $(LIBS)
81
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
82
+
83
+fl_draw: main_draw.c $(LIBS)
84
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
85
+
86
+
87
+#
88
+# Web exe
89
+# 
90
+
91
+fl_webdraw.cgi: main_webdraw.c  $(LIBS)
92
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
93
+
94
+fl_webmodel.cgi: main_webmodels.c   $(LIBS)
95
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
96
+
97
+fl_drawfile.cgi:	web/draw_custom.c  $(LIBS)
98
+	$(CC) $(ALLCFLAGS) -o $@ $(INCLUDE_DIR) $(LIBS_DIR) $< $(LFLAGS) -lm;
99
+
100
+#
101
+# Librarys
102
+#
103
+
104
+liblandscape.a: $(OBJECTS)