https://github.com/cran/GenAlgo
Tip revision: 35a421aca28768b31853a0b29cdbab60e40e7569 authored by Kevin R. Coombes on 18 May 2018, 14:29:40 UTC
version 2.1.4
version 2.1.4
Tip revision: 35a421a
featureSelection.html
<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<title>Simulating a Data Set</title>
<script type="text/javascript">
window.onload = function() {
var imgs = document.getElementsByTagName('img'), i, img;
for (i = 0; i < imgs.length; i++) {
img = imgs[i];
// center an image if it is the only element of its parent
if (img.parentElement.childElementCount === 1)
img.parentElement.style.textAlign = 'center';
}
};
</script>
<!-- Styles for R syntax highlighter -->
<style type="text/css">
pre .operator,
pre .paren {
color: rgb(104, 118, 135)
}
pre .literal {
color: #990073
}
pre .number {
color: #099;
}
pre .comment {
color: #998;
font-style: italic
}
pre .keyword {
color: #900;
font-weight: bold
}
pre .identifier {
color: rgb(0, 0, 0);
}
pre .string {
color: #d14;
}
</style>
<!-- R syntax highlighter -->
<script type="text/javascript">
var hljs=new function(){function m(p){return p.replace(/&/gm,"&").replace(/</gm,"<")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
hljs.initHighlightingOnLoad();
</script>
<style type="text/css">
body, td {
font-family: sans-serif;
background-color: white;
font-size: 13px;
}
body {
max-width: 800px;
margin: auto;
padding: 1em;
line-height: 20px;
}
tt, code, pre {
font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
}
h1 {
font-size:2.2em;
}
h2 {
font-size:1.8em;
}
h3 {
font-size:1.4em;
}
h4 {
font-size:1.0em;
}
h5 {
font-size:0.9em;
}
h6 {
font-size:0.8em;
}
a:visited {
color: rgb(50%, 0%, 50%);
}
pre, img {
max-width: 100%;
}
pre {
overflow-x: auto;
}
pre code {
display: block; padding: 0.5em;
}
code {
font-size: 92%;
border: 1px solid #ccc;
}
code[class] {
background-color: #F8F8F8;
}
table, td, th {
border: none;
}
blockquote {
color:#666666;
margin:0;
padding-left: 1em;
border-left: 0.5em #EEE solid;
}
hr {
height: 0px;
border-bottom: none;
border-top-width: thin;
border-top-style: dotted;
border-top-color: #999999;
}
@media print {
* {
background: transparent !important;
color: black !important;
filter:none !important;
-ms-filter: none !important;
}
body {
font-size:12pt;
max-width:100%;
}
a, a:visited {
text-decoration: underline;
}
hr {
visibility: hidden;
page-break-before: always;
}
pre, blockquote {
padding-right: 1em;
page-break-inside: avoid;
}
tr, img {
page-break-inside: avoid;
}
img {
max-width: 100% !important;
}
@page :left {
margin: 15mm 20mm 15mm 10mm;
}
@page :right {
margin: 15mm 10mm 15mm 20mm;
}
p, h2, h3 {
orphans: 3; widows: 3;
}
h2, h3 {
page-break-after: avoid;
}
}
</style>
</head>
<body>
<pre><code>## Loading required package: Umpire
</code></pre>
<p>In this vignette, we ilustrate how to apply the <code>GenAlgo</code> package
to the problem of feature selection in an “omics-scale” data set. We
start by loading the packages that we will need.</p>
<pre><code class="r">library(GenAlgo)
library(Umpire)
library(oompaBase)
</code></pre>
<h1>Simulating a Data Set</h1>
<p>We will use the <code>Umpire</code> package to simulate a comnplex enough data set
to stress our feature selection algorithm. We begin by setting up a
time-to-event model, built on an exponential baseline.</p>
<pre><code class="r">set.seed(391629)
sm <- SurvivalModel(baseHazard=1/5, accrual=5, followUp=1)
</code></pre>
<p>Next, we build a “cancer model” with six subtypes.</p>
<pre><code class="r">nBlocks <- 20 # number of possible hits
cm <- CancerModel(name="cansim",
nPossible=nBlocks,
nPattern=6,
OUT = function(n) rnorm(n, 0, 1),
SURV= function(n) rnorm(n, 0, 1),
survivalModel=sm)
### Include 100 blocks/pathways that are not hit by cancer
nTotalBlocks <- nBlocks + 100
</code></pre>
<p>Now define the hyperparameters for the models.</p>
<pre><code class="r">### block size
blockSize <- round(rnorm(nTotalBlocks, 100, 30))
### log normal mean hypers
mu0 <- 6
sigma0 <- 1.5
### log normal sigma hypers
rate <- 28.11
shape <- 44.25
### block corr
p <- 0.6
w <- 5
</code></pre>
<p>Now set up the baseline “Engine”.</p>
<pre><code class="r">rho <- rbeta(nTotalBlocks, p*w, (1-p)*w)
base <- lapply(1:nTotalBlocks,
function(i) {
bs <- blockSize[i]
co <- matrix(rho[i], nrow=bs, ncol=bs)
diag(co) <- 1
mu <- rnorm(bs, mu0, sigma0)
sigma <- matrix(1/rgamma(bs, rate=rate, shape=shape), nrow=1)
covo <- co *(t(sigma) %*% sigma)
MVN(mu, covo)
})
eng <- Engine(base)
</code></pre>
<p>We alter the means if there is a hit, or else build it using the original
engine components.</p>
<pre><code class="r">altered <- alterMean(eng, normalOffset, delta=0, sigma=1)
object <- CancerEngine(cm, eng, altered)
summary(object)
</code></pre>
<pre><code>## A 'CancerEngine' using the cancer model:
## --------------
## cansim , a CancerModel object constructed via the function call:
## CancerModel(name = "cansim", nPossible = nBlocks, nPattern = 6, SURV = function(n) rnorm(n, 0, 1), OUT = function(n) rnorm(n, 0, 1), survivalModel = sm)
##
## Pattern prevalences:
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
##
## Survival effects:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.13455 -0.06092 0.30563 0.27837 0.73876 2.51446
##
## Outcome effects:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.335728 -0.590101 0.009406 -0.104634 0.143112 1.454560
## --------------
##
## Base expression given by:
## An Engine with 120 components.
##
## Altered expression given by:
## An Engine with 120 components.
</code></pre>
<pre><code class="r">rm(altered, base, blockSize, cm, eng, mu0, nBlocks, nTotalBlocks,
p, rate, rho, shape, sigma0, sm, w)
</code></pre>
<p>Now we can use this elaborate setup to generate the simulated data.</p>
<pre><code class="r">train <- rand(object, 198)
tdata <- train$data
pid <- paste("PID", sample(1001:9999, 198+93), sep='')
rownames(train$clinical) <- colnames(tdata) <- pid[1:198]
</code></pre>
<p>Of course, to make things harder, we will add noise to the simulated measurements.</p>
<pre><code class="r">noise <- NoiseModel(3, 1, 1e-16)
train$data <- log2(blur(noise, 2^(tdata)))
sum(is.na(train$data))
</code></pre>
<pre><code>## [1] 0
</code></pre>
<pre><code class="r">rm(tdata)
summary(train$clinical)
</code></pre>
<pre><code>## CancerSubType Outcome LFU Event
## Min. :1.000 Bad :106 Min. : 0.00 Mode :logical
## 1st Qu.:2.000 Good: 92 1st Qu.: 3.00 FALSE:88
## Median :4.000 Median :19.00 TRUE :110
## Mean :3.662 Mean :22.32
## 3rd Qu.:5.000 3rd Qu.:36.75
## Max. :6.000 Max. :67.00
</code></pre>
<pre><code class="r">summary(train$data[, 1:3])
</code></pre>
<pre><code>## PID7842 PID6153 PID2085
## Min. : 1.410 Min. : 1.337 Min. : 1.724
## 1st Qu.: 5.021 1st Qu.: 5.008 1st Qu.: 5.071
## Median : 6.084 Median : 6.048 Median : 6.117
## Mean : 6.127 Mean : 6.099 Mean : 6.163
## 3rd Qu.: 7.172 3rd Qu.: 7.122 3rd Qu.: 7.192
## Max. :11.705 Max. :12.687 Max. :12.170
</code></pre>
<p>Now we can also simualte a validation data set.</p>
<pre><code class="r">valid <- rand(object, 93)
vdata <- valid$data
vdata <- log2(blur(noise, 2^(vdata))) # add noise
sum(is.na(vdata))
</code></pre>
<pre><code>## [1] 0
</code></pre>
<pre><code class="r">vdata[is.na(vdata)] <- 0.26347
valid$data <- vdata
colnames(valid$data) <- rownames(valid$clinical) <- pid[199:291]
rm(vdata, noise, object, pid)
summary(valid$clinical)
</code></pre>
<pre><code>## CancerSubType Outcome LFU Event
## Min. :1.000 Bad :54 Min. : 0.00 Mode :logical
## 1st Qu.:2.000 Good:39 1st Qu.: 3.00 FALSE:28
## Median :3.000 Median :15.00 TRUE :65
## Mean :3.172 Mean :22.32
## 3rd Qu.:5.000 3rd Qu.:35.00
## Max. :6.000 Max. :71.00
</code></pre>
<pre><code class="r">summary(valid$data[, 1:3])
</code></pre>
<pre><code>## PID6584 PID5256 PID2944
## Min. : 1.287 Min. : 1.317 Min. : 1.760
## 1st Qu.: 4.957 1st Qu.: 4.997 1st Qu.: 5.022
## Median : 5.995 Median : 6.030 Median : 6.064
## Mean : 6.067 Mean : 6.083 Mean : 6.121
## 3rd Qu.: 7.109 3rd Qu.: 7.116 3rd Qu.: 7.151
## Max. :12.205 Max. :12.223 Max. :12.341
</code></pre>
<h1>Setting up the Genetic Algorithm</h1>
<p>Now we can start using the <code>GenAlgo</code> package. The key step is to define
sensible functions that can measure the “fitness” of a solution and to
introduce “mutations”. When these functions are called, they are passed a
<code>context</code> argument that can be used to access extra information about
how to proceed. In this case, that context will be the <code>train</code> object,
which includes the clinical information about the samples.</p>
<h2>Fitness</h2>
<p>Now we can define the fitness function. The idea is to compute the Mahalanobis
distance between the two groups (of “Good” or “Bad” outcome samples) in the
space defined by the selected features.</p>
<pre><code class="r">measureFitness <- function(arow, context) {
predictors <- t(context$data[arow, ]) # space defined by features
groups <- context$clinical$Outcome # good or bad outcome
maha(predictors, groups, method='var')
}
</code></pre>
<h2>Mutations</h2>
<p>The mutation function randomly chooses any other feature/row to swap out
possible predictors of the outcome.</p>
<pre><code class="r">mutator <- function(allele, context) {
sample(1:nrow(context$data),1)
}
</code></pre>
<h2>Initialization</h2>
<p>We need to decide how many features to include in a potential predictor
(here we use ten). We also need to decide how big a population of feature-sets
(here we use 200) should be used in each generation of the genetic algorithm.</p>
<pre><code class="r">set.seed(821831)
n.individuals <- 200
n.features <- 10
y <- matrix(0, n.individuals, n.features)
for (i in 1:n.individuals) {
y[i,] <- sample(1:nrow(train$data), n.features)
}
</code></pre>
<p>Having chosen the staring population, we can run the first step of the
genetic algorithm.</p>
<pre><code class="r">my.ga <- GenAlg(y, measureFitness, mutator, context=train) # initialize
summary(my.ga)
</code></pre>
<pre><code>## An object representing generation 1 in a genetic algorithm.
## Population size: 200
## Mutation probability: 0.001
## Crossover probability: 0.5
## Fitness distribution:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02137 0.14909 0.22210 0.24071 0.31384 0.58730
</code></pre>
<p>To be able to evaluate things later, we save the starting generation</p>
<pre><code class="r">recurse <- my.ga
pop0 <- sort(table(as.vector(my.ga@data)))
</code></pre>
<h2>Multiple Generations</h2>
<p>Realistically, we probably want to run a couple of hundred or even a
couple thousand iterations of the algorithm. But, in interests of making
the vignette complete ins a reasonable amount of time, we are only going
to terate through 20 generations.</p>
<pre><code class="r">NGEN <- 20
diversity <- meanfit <- fitter <- rep(NA, NGEN)
for (i in 1:NGEN) {
recurse <- newGeneration(recurse)
fitter[i] <- recurse@best.fit
meanfit[i] <- mean(recurse@fitness)
diversity[i] <- popDiversity(recurse)
}
</code></pre>
<p>Plot max and mean fitness by generation. This figure shows that both the mean
and the maximum fitness are increasing.</p>
<pre><code class="r">plot(fitter, type='l', ylim=c(0, 1.5), xlab="Generation", ylab="Fitness")
abline(h=max(fitter), col='gray', lty=2)
lines(fitter)
lines(meanfit, col='gray')
points(meanfit, pch=16, col=jetColors(NGEN))
legend("bottomleft", c("Maximum", "Mean"), col=c("black", "blue"), lwd=2)
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAA2FBMVEUAAAAAADoAAGYAAJkAAMwAAP8AM/8AOmYAOpAAZrYAZv8Amf8AzP8A//8z/8w6AAA6ADo6AGY6OmY6OpA6ZpA6ZrY6kJA6kNtmAABmADpmAGZmOgBmOjpmZjpmZmZmZrZmtttmtv9m/5mQOgCQOjqQOmaQZgCQZpCQkGaQtpCQ29uQ2/+Z/2a2ZgC2Zjq2tma225C2/7a2/9u2//++vr7MAADM/zPbkDrbkGbb25Db/7bb/9vb////AAD/MwD/ZgD/mQD/tmb/zAD/25D//wD//7b//9v///8x3bEkAAAACXBIWXMAAAsSAAALEgHS3X78AAAODUlEQVR4nO3diVrb2AGGYUGWQjZCaQph2mkJSTsdhzZTnNBOySSOY+v+76g+krxv0lmkI/3f98wEjC1LzpsjyRa2kpQkS5peAGom4EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRQNeNOBFA1404EUDXjTgRXOBTyjmAsI7TEuhA1404EUDXjTgRXOFH13epqOLJDm6rz4tNZgHeGOfDl9Vn5YazAP88Ow+H/lVp6UGc4a/OPj5nRnxZ2vreuBjzn3nbnydHKeDw7UBD3zUsVcvGvCi+YJf2LkreRiAGo0RLxrwogEvGvCiub+AU+zJrT+RBz7mnEf8+PrcelpqLvdV/eh1z3paaiy28aIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBL1pA+L6Jr5F8TZKVy7vxGPFdKJm09rM9k7jMzmFa8tYG9OzHe6ZymaPDtOSnLeop8B1u0wp+4do9E7vM2GFacmwnenYDp6uDTUsu7VVPge9gJdRT4DtXOXbgu1WZdfz0pk5XB5uWLCqvngLfnSqxA9+RKqzjp1M4XR1sWqpSZfXUA/zwxLxClBzeWkxLPrJhd4cfX19lXwdH95WnJbeSPMuJna5O09Hl7dLXKtOSTckst7txulpwxDv+pZedNNma7ZxXZ+B09aTRhcg2fvnv3Q4hm2L/dP54d8zD6epg00bVduRqg3B+w13T1KGeAr+vUq6l9FdusfVfUtklc8wX/MLOne+tUZNVeAy7H/CGK9cnqPOvjBG/q4oPYav99p8nG7+vIeB3ZPOC2IZp9qwKyu3w+Q747dk9gFXDMlv/Bv6qgN+a/fIni/vvXpbFf86v3F0Ue3LrT+Rjfcglc3xhrFiD+1mWADmP+PH1ufW0Mee+9HE/r3Ff1Y9e96ynjbdWL3yZ2MZvKuqx6ifgN9TeJS8f8Ou1dsGrBPxabV3uagG/WksXu2rAr9TOpa4e8Mu1cqFtAn6pNi6zXcAv1sJFtg34hdq3xPYBPy3ul9a9B7ypI78oViXgBdFN2vCi6CZdeGF0U4vht7/LqFSBly72WgsvL+dYW+Fhd6yd8Ax351oJD7t7LYRnuPuodfCw+6lt8LB7ql3wDHdvtQoedn+1CJ7h7rP2wMPutbbAM9w9Fx88h1RqKTp4iOspNnjca6oc/N3R/V2SXHm9a48TUfVKwY9e9yb/DU/XP7bU4a69TUM2lYO/vJ2M+Rrgca+tkqv65KA3CL+qx72+otq5A76+YoLHPWAfP35cuhzRXj3uAfv4cUU+nr163ENmCV/HXj3w+/ry5YvNZH2THXwde/W47+vLl7Lynz9/zr/JzPNv7bbxVlWaFve9ZfBTx2m//fbb2g0/mxbMNxYJPO776veLEd+flWbuc/npzzP4ffdXDn58nRz9uu0zay3v2vKmkuXKq2v6yU8z+IV/CFn+4MfX58Oz+w2nlnO5a8ub6rVrlb004meVcC+9Vz+B33AySZe7trqlXjs31Ju38aWqMOLvQo143Nf7+vVrunuwO1Z6G58kFd1Le+K+3ldTOPU0ir164NfL4IPOoXl43JfLdtBjgR9sOd+Qy11XulWn+/btW/7NwpOy0O4l9+ovKr5aW+auK9yo230z7XuhzXtln875v+sKN+p2GXztcy23qr/Zdooxh7sufZsuZ8Z5xPDFWQXdtvG8QWal2cq9CfcI9uq1+v79u/lS9xZ9vQrb+HAv2er03dQ4uqkE/Oz0seEO0siUwTe9EFlN79XrlD1haxl8kLvufJ8+fcq/WTxaHot7qVX95S8+9urV+mSq/XWZ0jHiQ5XBN70Q2wM+TP1+++Htdu2U4fPVe9TuZeGHZxV/C2P/XXe1WDfqKwHvt3aop8B76v3792lrBntWlVfueDq3rfemFqmn7NVb9ObNm6XLk3GewTe0OJYBX7U3bzL5pXcyAe9r2ojL4FdX661zB75q/WLEtz3gq5St2DvhDnyF2rXbvifgS9am5+hlAr5UHVNPgd/TixcvujfYs4Df1QtTB9VT4Bd6/vz50uXJQM/gG1qcwAE/7fnzTH7pFTngA9x1ZPUz+NWteWfdReGfPXs2+346vosRr5Ik/DNTOl+j50m5C8N3c2e9dLrwTS9EwwnCT9buuOvBd/FVOJvE4GGfJgUP+zwleNgX0oFnuC+lAg/7St2Hf/r0KezrdR7+qanphYgw4EXrNny/3wd+cx2GL4684b6xrsJ38hckfdYt+CdPnqQxfF5oC+oU/BMT6KVyhTcfl2E+OWHDx502BF/7XNuZB/jsU7GGr6pP67s+8OXzAJ99Ps6Gz0SrGd6s4nEvnTP8xcHP78yIX/90pFrh2bJXzH3nbnydHKeDDR+MVCM87JXrwl497Ba1Hx52q3zBL+zc1XqmIdgta++If/z4Mez2tRb+sSnoHLqdM/zwZNvHnQIfc67w4+v8tLOD9ddsQ8L3+8C75eO1+sWvVaa1LT/4hrtTrRvxHHL1k/M2fvt5Z0PAo+6rduzVP3r0KGWwe60V8I9MqHutPfDe7o1MbYDvA++/+OHNKh5370UOz/5cqKKGRz1cEcPDHrJY4VnHBy46+IcPHzLYayg2+Icm2MMXJbzDTKlkwIsWGXw/xb2eooJnn66+YoKHvcbigWe411os8LDXXCTwsNddFPAM9/qLAR72BmoenuHeSI3CP3jwgOHeUE3CPzA5zIAcAl404EVrDn6yU4d7czUFz758wzUDD3vjNQEPewQ1AA97DNUOz3CPo5rhYY+l2uAPDw9hj6i64A9NDndGngNeNOBFq3UbT/HU/C9iUCMBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwooWEryWHBZSu7SMeeMuAFw140YAXDXjRgBcNeNGAFw140YAXLQ744cnx5M+bw9ulH57ebrm53UxosUjgX/7hPh39qYy0/UxosXrg9x1fGZ7+uZcOf5zAD0+S5Cq9OU/vjicjfvz2pyQ5H0z+z8b/8k+cF1C6WuD3Hlgbnv7rKv3PP05vR697hnd0+cvlrWG+Ps42A5Nvp/DznzgvoHSRrOpPP7wav/uQY44ub9O7fIyP3/ZS8//k38NsxM9+4ryA0sUC/++///qXbBTfJMlkH2/4sgd80GKBv/3nT+dmHX9xla3Fb348Bj5o0cAPDnrTLfnL3vDsf297y8xG+u4QeF9FA1+o3yXJ73746/VVOjj67/L4nlzzx0vgfRUHfOQz6WLAi+YMb15xSbI98erTegh4y1zhx5PNsWlwdF95Wh8Bb5krvHm1ZfFrlWl9BLxljHjRnLfxowsP2/iNh2VLBbxlcezVc1i29uqBf7Sx+fXzw7KTFYgZ9/nR2eHZ38yXUAsonS/4hZ279cOum92X4IvDsvmR+LQ4Ojs8Od+082C1gLRcJKv66WFZI178GxpdTo/Bh1pA6WKBLw7LZruKB73i6Czw4Yrjlbv5YdlitOdHZ4EPVxzP42eHZc02fnJX06OzwAcrjlfu5odlJ+t6s6bPjs5eAR+uOEa8fcBbFscrd/YBb1kce/WRz6SLAS8a8KIBLxrwogEvWkh4irlw8KHrzholwpkALzoT4EVnArzoTIAXnQnwojMBXnQmMcNTwIAXDXjRgBcNeNGAFw140YAXDXjRgBctWvi7ZOObOLyWvXtrdJHsfoe+j5mEfTT5p0xUeyTRwt/s/qAMHw2MhXmL2N1x4JmEfTTZp0y87FV7JLHCm4+zDdzNgfmEDvNm0FInx3GZSdhHM8g+WOqq2iOJFT57x17oQZ+9WffsfvaBycFmEv7RmM99r/RIYoU3JzIIPuqNiXkXcHD44I9mfH1e8ZHECp8Vejtf24jPCvhoRhfnacVHIg8feBtfC/zwxNx1N7bxZr01flfD0zmzkgy5Vz/bnoR7NLl7xUcSK7x55nsQese+1ufx4R7NXfbGmauOPI+nsAEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvGvCiAS8a8KIBLxrwogEvmhb8IEmSEr+HWpwir9NJwQ8O81+n3lfn1VMt+PzcidkbTrKT2E/PXl9c+v0Ph7fZG47N5Q+zU92XOMd9G1OCnw/k/CT207PXTy9dzU5rn6/qi9Pg7j/HfRuTgj8r+IqT2Bdntl24lF15WcCbtySZN6HuP/9tG5OCN4QnyUGvOIn9FH5+aX5a+2LQj98C3/7ybfxkEBdnxp7Czy/NT2vPiO9SA/P+NfPHwknsl09pPzut/frPO5YUfPYpQYf5GyUPetMRv3CpOK39+Hpxrx546lDAiwa8aMCLBrxowIsGvGjAiwa8aMCLBrxowIsGvGjAiwa8aMCLBrxowIv2f4BomHT3Dc0hAAAAAElFTkSuQmCC" alt="Fitness by generation."/></p>
<p>Plot the diversity of the population, to see that it is deceasing.</p>
<pre><code class="r">plot(diversity, col='gray', type='l', ylim=c(0,10), xlab='', ylab='', yaxt='n')
points(diversity, pch=16, col=jetColors(NGEN))
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAAkFBMVEUAAAAAADoAAGYAAJkAAMwAAP8AM/8AOpAAZrYAZv8Amf8AzP8A//8z/8w6AAA6ADo6AGY6OpA6kNtmAABmADpmAGZmtv9m/5mQOgCQOjqQOmaQtpCQ2/+Z/2a2ZgC2//++vr7MAADM/zPbkDrb////AAD/MwD/ZgD/mQD/tmb/zAD/25D//wD//7b//9v////PBLrnAAAACXBIWXMAAAsSAAALEgHS3X78AAAIU0lEQVR4nO3d60ITRwBA4cjNcrOttmBb0QUlILd9/7drdoIRFCGb7Oxczjk/0ezs+DmZWUUzaQ3ZJPUNWJqEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIe2DvzEci4i/BqvtdgJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMe2mjwGxsba1zLhm4s+I0N5bNqVPhmjavZsI274pvQ4oubm5trXN/WKcUe39z7b24qn6yEp/pG+ISlfJwL8I0bf5KSPsfPF/yjbd9GKpM/wNF+7DKBbxf2W1tba4xpy5YPfBvst7aUH6Ws4Nuw4IUfoyzhPe7FLzf4xR4vftyyg3/YAn97e3vti9mjsobv6vC3t5Ufuuzhu4QfvnLg3fIHrQj4+z3e896AlQH/Pe0HqjT4dmG/s7MTaQBEBcK3wX5nR/l1KhO+DQte+DUqG94Nf+WKhb/f4z3pr1i58N/TfoVqgG8X9ru7uyMOWnSVwLfBfndX+WWrB74NC174JasP3g1/qaqCv9/jPewtUV3wi6R/qUrhXfYvVS18+81+b28v9Y3kWM3wbWe/t6f8U1UOHxa88E8kPLTq4cMe70Hvp+qHnyf9D1Hgpf8hDrz0jyLBS/8gFrz0i2jwc/r9/f3Ut5E6HvyMfn9feSJ8G+B/+lucg4ODNLeTJC58u/gvNue/BQ4OUPJI+J/3+Jm98D1+ONprExTgG8z36Qu/aLHgH+gfHh6mu6GoCf+Lgv7hYbXywj+X8BEuXULCR7h0EXXudZ72hF+iGumFX6r66IVfstrkhV+2yha98MtXFb3wfaqIXvh+zeiPj49T38UACd+34+Mq5IXvm/BrXrrUAnz5e73wvQsLvvhjnvCrVji98KtXNL3w61Tw92kJv2al0gu/dh39yclJ6tvomfAD1JycFCcv/BAJP9RrC0v4oV5bWp17Wcc84QesJHnhh6ygRS/8sBVDL/zQFSIv/OCVseiFj1AJ9MJHKX954eM0W/RnZ2ep7+KZhI/V2VnW8sLHSnhoAT7fvV74aIUFn+336AgfvTzthR+jDOmFH6fslr3woxXsb25uUt/HPOHHrLm5yUVe+FETHlqAz2K7F37c5gs+g6Oe8IlKTS98stIue+FTNre/uroaf2jhE9c0V1cp5IVPXoAf/V1f+OTNV3zz+FNRLi8v444qfPoevNN/+1Csy8vY8sLnWCM8tAAf9XFP+DybL/iIH4YmfO5Fshe+gO7tLy4uhrum8GU0s7+4eEJ+Op2udkHhiynAN4+f96fTVeWFL6aHK/7bp2ALT+jpd/rpdKXDn/BlFxb8Kgd/4euot73w1dSPXviK6rPsha+rYH9+fv7iTxS+uprz8yXkha8v4aEF+Jd2e+ErLCz4Fw56wtfbs/TC19wz9MLX3S8f7YWvvqfphQfU0Z+enj76mvCImtPTH+SFZyQ8NOGpucdbSHhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NBiwlvOxYOP3Sj3Rh1EeOggwkMHER46iPDQQYSHDiI8dJCc4S1iwkMTHprw0ISHJjw04aEJD014aMJDyxb+82Qy2fgUd4zrP2YD3L6bvP4Se5C4s7l+M5kc9ZxJtvAfj6IP8bWzuHt/1H7+LfIgcWdz+9eH9vr3D/1mkiv83T8fYg/x8dV/s8V4+/en+aKMOUjc2XzttD8e9ZtJrvCzt63w9hW17pfp+s8vYclEHST+bGZT6DeTXOFnb13xV31n8vX1CPDRZ3P3/m3PmeQKH4q9z4+24kMRZ3P77m3bcyZ4+Mh7/Cjw12+6S9exx3fvW3f/jvA4171JxjzVL/aTeLOZu/ecSa7w3ZPvq9gH+1Gf4+PN5nP4hzNHlTzHW9yEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJDEx6a8NCEhyY8NOGhCQ9NeGjCQxMemvDQhIcmPDThoQkPTXhowkMTHprw0ISHJjw04aEJD014aMJD+x+2uL+O/nbV2QAAAABJRU5ErkJggg==" alt="Diversity."/></p>
<p>See which predictors get selected most frequently in the latest
generation.</p>
<pre><code class="r">sort(table(as.vector(recurse@data)))
</code></pre>
<pre><code>##
## 1137 2839 3318 5599 9337 2852 3439 4248 5851 7642 7720 9039
## 1 1 1 1 1 2 2 2 2 2 2 2
## 9771 11427 49 918 1599 2390 5344 6130 6333 8322 8725 8807
## 2 2 3 3 3 3 3 3 3 3 3 3
## 9150 9512 10186 10884 11139 1958 4632 4991 5094 6601 9728 11091
## 3 3 3 3 3 4 4 4 4 4 4 4
## 12304 93 456 3801 3857 6386 6657 7680 8741 10825 11512 4100
## 4 5 5 5 5 5 5 5 5 5 5 6
## 9044 12161 2845 6022 9923 7846 5152 5473 9574 81 10220 10965
## 6 6 7 7 7 8 9 9 9 11 11 11
## 2335 3055 9565 10435 1176 6440 8749 449 10765 1534 6959 7306
## 12 12 12 12 14 15 16 17 17 18 18 19
## 11164 7000 8651 3838 8973 9524 5086 8832 11464 8683 1451 4976
## 20 21 21 22 23 24 27 27 28 29 32 34
## 4953 6568 8926 7348 1394 9874 10247 171 6117 8600 1855 1935
## 38 44 49 61 62 65 72 84 88 101 120 152
## 936 6252
## 158 159
</code></pre>
</body>
</html>
