आर

Question

FYI में तेज़ कोड: मैंने अपने पहले संस्करण के बाद इसे महत्वपूर्ण रूप से संपादित किया है। यह अनुकरण 14 घंटे से 14 मिनट लेने में कम कर दिया गया है।

मैं प्रोग्रामिंग के लिए नया हूं लेकिन मैंने एक सिमुलेशन बनाया है जो जीव में असमान प्रतिकृति का पालन करने की कोशिश करता है और माता-पिता और बेटी जीवों के बीच गुणसूत्र संख्या में मतभेदों को मापता है। सिमुलेशन बहुत धीमी गति से चलता है। इसे पूरा करने में लगभग 6 घंटे लगते हैं। मैं जानना चाहता था कि सिमुलेशन को तेजी से चलाने का सबसे अच्छा तरीका क्या होगा।

इन डिजिटल जीवों में गुणसूत्रों की संख्या x है। अधिकांश जीवों के विपरीत गुणसूत्र एक-दूसरे से स्वतंत्र होते हैं, इसलिए उनके पास बेटी जीव में स्थानांतरित होने का बराबर मौका होता है।

इस मामले में, एक बेटी सेल में गुणसूत्रों का वितरण 0.5 की संभावना के साथ एक द्विपदीय वितरण का पालन करता है।

फ़ंक्शन sim_repo गुणसूत्रों की एक ज्ञात संख्या के साथ डिजिटल जीवों का एक मैट्रिक्स लेता है और उन्हें प्रतिकृति के 12 पीढ़ियों के माध्यम से रखता है। यह इन गुणसूत्रों को डुप्लिकेट करता है और फिर rbinom फ़ंक्शन का उपयोग यादृच्छिक रूप से एक संख्या उत्पन्न करने के लिए करता है। यह संख्या तब बेटी सेल को सौंपी जाती है। चूंकि असामान्य प्रजनन के दौरान कोई गुणसूत्र खो नहीं जाता है, अन्य बेटी कोशिका शेष गुणसूत्र प्राप्त करती है। इसके बाद जी पीढ़ियों की संख्या के लिए दोहराया जाता है। फिर मैट्रिक्स में प्रत्येक पंक्ति से एक एकल मान नमूना लिया जाता है।

sim_repo = function(x1, G=12, k=1, t=25, h=1000) { 

      # x1 is the list of copy numbers for a somatic chromosome 
      # G is the number of generations, default is 12 
      # k is the transfer size, default is 1 
      # t is the number of transfers, default is 25 
      # h is the number of times to replicate, default is 1000 

      dup <- x1 * 2 # duplicate the initial somatic chromosome copy number for replication 
      pop <- 1 # set generation time 
      set.seed(11) 
      z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) # amount of somatic chromosome is distributed to one of the daughter cells 
      z1 <- dup - z # as no somatic chromosomes are lost, the other daughter cells receives the remainder somatic chromosomes 
      x1 <- cbind(z, z1) # put both in a matrix 

      for (pop in 1:G) { # this loop does the replication for each cell in each generation 
       pop <- 1 + pop # number of generations. This is a count for the for loop 
       dup <- x1 * 2 # double the somatic chromosomes for replication 
       set.seed(11) 
       z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) # amount of somatic c hromosomes distributed to one of the daughter cells 
       z1 <- dup - z # as no somatic chromosomes are lost, the other daughter cells receives the remainder somatic chromosomes 
       x1 <- cbind(z, z1) # put both in a matrix 
       } 

      # the following for loop randomly selects one cell in the population that was created 
      # the output is a matrix of 1 column 
      x1 <- matrix(apply(x1, 1, sample, size=k), ncol=1) 
      x1 
    } 

मैं अपने अनुसंधान मैं प्रारंभिक पैतृक जीवों और इस अनुकरण में अंतिम समय बिंदु के गुणसूत्रों में विचरण में परिवर्तन में दिलचस्पी है। निम्न कार्य एक सेल को एक नए जीवित वातावरण में स्थानांतरित करने का प्रतिनिधित्व करता है। यह sim_re पी फ़ंक्शन से आउटपुट लेता है और अधिक पीढ़ियों को उत्पन्न करने के लिए इसका उपयोग करता है। इसके बाद पहले और अंतिम मैट्रिक्स कॉलम में पंक्तियों के बीच भिन्नता पाई जाती है और उनके बीच अंतर मिलता है।

# The following function is mostly the same as I talked about in the description. 
    # The only difference is I changed some aspects to take into account I am using 
    # matrices and not lists. 
    # The function outputs the difference between the intial variance component between 
    # 'cell lines' with the final variance after t number of transfers 

sim_exp = function(x1, G=12, k=1, t=25, h=1000) { 

    xn <- matrix(NA, nrow(x1), t) 
    x <- x1 
    xn[,1] <- x1 
    for (l in 2:t) { 
     x <- sim_repo(x, G, k, t, h) 
     xn[, l] <- x 
    } 

    colvar <- matrix(apply(xn,2,var),ncol=ncol(xn)) 
    ivar <- colvar[,1] 
    fvar <- colvar[,ncol(xn)] 
    deltavar <- fvar - ivar 
    deltavar 
} 

मैं इस अनुकरण ज समय की संख्या को दोहराने की जरूरत है। इस प्रकार मैंने निम्न कार्य किया जो फंक्शन sim_expएच कई बार कॉल करेगा।

sim_1000 = function(x1, G=12, k=1, t=25, h=1000) { 
    xn <- vector(length=h) 
    for (l in 2:h) { 
     x <- sim_exp(x1, G, k, t, h) 
     xn[l] <- x 
    } 
     xn 
} 

जब मैं sim_exp फ़ंक्शन को 6 मानों के मान के साथ कॉल करता हूं तो इसे पूरा करने में लगभग 52 सेकंड लगते हैं।

x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) 
system.time(sim_1000(x1,h=1)) 
    user system elapsed 
    1.280 0.105 1.369 

यदि मैं इसे तेज़ी से प्राप्त कर सकता हूं तो मैं इन सिमुलेशन में से अधिक पूरा कर सकता हूं और सिमुलेशन पर एक चयन मॉडल लागू कर सकता हूं।

मेरे इनपुट निम्नलिखित x1, उसकी पंक्ति में प्रत्येक पैतृक जीव के साथ एक मैट्रिक्स तरह दिखेगा:

x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) # a matrix of 6 organisms 

जब मैं चलाएँ:

a <- sim_repo(x1, G=12, k=1) 

मेरे उम्मीद उत्पादन होगा:

a 
    [,1] 
[1,] 137 
[2,] 82 
[3,] 89 
[4,] 135 
[5,] 89 
[6,] 109 

system.time(sim_repo(x1)) 
    user system elapsed 
    1.969 0.059 2.010 

जब मैं sim_exp फ़ंक्शन को कॉल करता हूं,

ख < - sim_exp (x1, जी = 12 k = 1, टी = 25)

यह कहता है sim_repo समारोह जी बार और आउटपुट:

b 
[1] 18805.47 

जब मैं sim_1000 फ़ंक्शन को कॉल करें, मैं आमतौर पर एच से 1000 सेट कर दूंगा, लेकिन यहां मैं इसे 2 पर सेट कर दूंगा। इसलिए यहां सिम_1000 sim_exp पर कॉल करेगा और इसे दो बार दोहराएगा।

c <- sim_1000(x1, G=12, k=1, t=25, h=2) 
c 
[1] 18805.47 18805.47 

Answer 1

, टिप्पणियों में अन्य लोगों ने उल्लेख किया है कि अगर हम समारोह sim_repo पर केवल देखने के लिए, और रेखा की जगह: लाइनों

dup <- x1 * 2 

साथ

dup <- apply(x1, c(1,2),"*",2) 

z <- apply(dup,c(1,2),rbinom,n=1,prob=0.5) 

z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) 

और

x1 <- matrix(apply(x1,1,sample,size = 1), ncol=1) 

साथ पाश के लिए आंतरिक के साथ मैं एक, अच्छी तरह से, बड़े गति वृद्धि मिलती है:

system.time(sim_exp(x1)) 
    user system elapsed 
    0.655 0.017 0.686 
> system.time(sim_expOld(x1)) 
    user system elapsed 
21.445 0.128 21.530 

और मैं सत्यापित है कि यह एक ही बात कर रहा है:

set.seed(123) 
out1 <- sim_exp(x1) 

set.seed(123) 
out2 <- sim_expOld(x1) 

all.equal(out1,out2) 
> TRUE 

और यह भी delving नहीं है प्री-आवंटन में, जो वास्तव में आपके कोड को संरचित करने के तरीके के अनुसार चीजों को पूरी तरह से फिर से डिजाइन किए बिना कठिन हो सकता है।

और वह भी नहीं पर आप वास्तव में सभी तीन कार्यों की जरूरत है कि क्या देखने के लिए ...