diff --git a/nbs/marl.ipynb b/nbs/marl.ipynb index 000869c..099ec3b 100644 --- a/nbs/marl.ipynb +++ b/nbs/marl.ipynb @@ -1,154 +1,30 @@ { "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "5a947995", + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, "id": "8d6b4579", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "(2, 3, 77)\n" - ] - }, - { - "data": { - "text/plain": [ - "array([[[-0. , -0.16740674, -0.06527205, 0.07791231,\n", - " -0.02334672, 0.10651681, -0.04032374, -0.09000908,\n", - " -0. , -0. , -0.14130992, -0.10716073,\n", - " -0. , -0. , 0. , 0.07216024,\n", - " 0. , 0.01171841, -0.10167756, -0.0860365 ,\n", - " 0. , -0.00251561, 0.16547332, 0.05970075,\n", - " 0. , 0. , 0. , 0.13273609,\n", - " 0.01546233, 0. , 0. , -0.07790096,\n", - " -0. , -0. , 0. , -0. ,\n", - " -0. , -0. , 0.08016346, 0. ,\n", - " -0. , 0.04080275, 0. , 0. ,\n", - " -0. , 0.05893229, 0. , 0.03338421,\n", - " -0. , -0.0037284 , 0. , 0.07843556,\n", - " -0.19221128, 0.11882027, 0.17016669, -0.19694884,\n", - " 0.13400745, 0. , -0. , -0. ,\n", - " 0. , -0. , 0. , -0. ,\n", - " -0. , 0. , -0.09685217, 0. ,\n", - " 0. , -0. , 0. , -0.11594658,\n", - " -0. , -0.08525832, 0.1080862 , -0.05032869,\n", - " 0.11070766],\n", - " [ 0. , -0.13369633, 0.03570481, -0. ,\n", - " 0. , -0. , -0. , 0.07353768,\n", - " 0.06817306, 0.06399359, 0.10936508, 0.02181026,\n", - " 0.07215239, -0.05659802, 0.0698835 , -0.17156897,\n", - " -0.10751673, -0.07604196, 0.0771776 , -0.02900459,\n", - " -0. , 0. , 0.00680535, 0. ,\n", - " 0.01449199, 0.04864531, 0.02191286, 0.0063002 ,\n", - " 0. , -0.01836898, 0. , -0. ,\n", - " 0. , 0.18066614, 0.06187952, -0.07484735,\n", - " -0. , 0.01840071, 0. , 0.12390175,\n", - " 0. , -0. , -0.02858164, 0.10949008,\n", - " -0.03000907, 0.04536685, -0. , -0.08964098,\n", - " -0. , -0.06219959, 0. , -0. ,\n", - " 0. , -0. , -0.11418024, -0. ,\n", - " -0. , 0.10855307, -0. , 0. ,\n", - " 0.07895877, 0.0450811 , 0. , 0. ,\n", - " -0.01556281, -0. , -0.04141774, 0. ,\n", - " -0.13117491, 0.14870693, -0. , -0. ,\n", - " -0.00230331, 0.10070572, 0. , -0.03584967,\n", - " -0.0358523 ],\n", - " [ 0.00346828, -0. , -0. , -0.06835983,\n", - " 0.15181963, -0.12193283, -0. , -0. ,\n", - " -0. , 0. , 0. , 0.04357489,\n", - " -0.06102486, 0.01911994, 0.01863411, 0.02873777,\n", - " -0. , -0. , 0.03432964, -0.07652427,\n", - " 0.1450196 , -0. , 0.01890378, 0. ,\n", - " 0. , 0.00632307, -0.02393692, 0. ,\n", - " 0.03479244, -0.20124773, -0. , -0.09872656,\n", - " 0.01674399, 0. , 0.05757198, -0. ,\n", - " 0.07212994, 0. , 0.04404994, 0. ,\n", - " -0. , 0. , -0.02517718, 0.03949109,\n", - " 0. , -0.10251922, 0.1057609 , 0. ,\n", - " 0. , -0.00070022, 0.01298889, -0. ,\n", - " -0.12222556, -0. , -0.03469444, 0. ,\n", - " 0. , 0. , -0. , 0.01280573,\n", - " -0.0581795 , 0. , -0. , 0. ,\n", - " 0.02232887, -0. , 0.13227224, -0.07421724,\n", - " 0.10009447, 0. , 0. , 0. ,\n", - " -0. , -0.06519596, -0. , 0.00896831,\n", - " 0.07257131]],\n", - "\n", - " [[-0. , -0.16740674, -0.06527205, 0.07791231,\n", - " -0. , 0. , -0. , -0.09000908,\n", - " -0.12261757, -0.04929118, -0. , -0.10716073,\n", - " -0.26610878, -0.07733363, 0. , 0.07216024,\n", - " 0.02077423, 0. , -0. , -0.0860365 ,\n", - " 0.06077938, -0.00251561, 0.16547332, 0.05970075,\n", - " 0. , 0.03002638, 0. , 0. ,\n", - " 0.01546233, 0. , 0. , -0.07790096,\n", - " -0.03088952, -0.07117025, 0.04056206, -0. ,\n", - " -0.00716205, -0. , 0. , 0.0040351 ,\n", - " -0. , 0.04080275, 0.01734744, 0.0983687 ,\n", - " -0.0151504 , 0. , 0.03784797, 0.03338421,\n", - " -0. , -0. , 0. , 0.07843556,\n", - " -0.19221128, 0. , 0. , -0.19694884,\n", - " 0. , 0.15261368, -0. , -0.07315045,\n", - " 0.17382503, -0.17040975, 0. , -0. ,\n", - " -0.10972709, 0.0094038 , -0.09685217, 0. ,\n", - " 0. , -0. , 0. , -0. ,\n", - " -0. , -0.08525832, 0. , -0. ,\n", - " 0.11070766],\n", - " [ 0.11801033, -0.13369633, 0. , -0.0687048 ,\n", - " 0.11730684, -0.0916871 , -0.04051587, 0.07353768,\n", - " 0. , 0.06399359, 0.10936508, 0. ,\n", - " 0.07215239, -0.05659802, 0. , -0. ,\n", - " -0. , -0. , 0.0771776 , -0. ,\n", - " -0.11570772, 0. , 0.00680535, 0.0482616 ,\n", - " 0. , 0. , 0.02191286, 0.0063002 ,\n", - " 0. , -0. , 0.05764654, -0.147537 ,\n", - " 0.0093513 , 0. , 0.06187952, -0. ,\n", - " -0. , 0.01840071, 0.19468755, 0.12390175,\n", - " 0. , -0. , -0.02858164, 0. ,\n", - " -0.03000907, 0.04536685, -0. , -0. ,\n", - " -0. , -0. , 0. , -0.0706455 ,\n", - " 0.02757949, -0. , -0. , -0. ,\n", - " -0. , 0.10855307, -0.07525277, 0. ,\n", - " 0. , 0.0450811 , 0.06722807, 0. ,\n", - " -0. , -0.05833585, -0. , 0.16977923,\n", - " -0.13117491, 0.14870693, -0. , -0.07954684,\n", - " -0. , 0.10070572, 0. , -0. ,\n", - " -0.0358523 ],\n", - " [ 0. , -0. , -0. , -0.06835983,\n", - " 0.15181963, -0.12193283, -0.09164272, -0. ,\n", - " -0. , 0.05030431, 0.10992882, 0. ,\n", - " -0.06102486, 0. , 0. , 0.02873777,\n", - " -0.10900734, -0.0848316 , 0. , -0. ,\n", - " 0. , -0.02109335, 0. , 0. ,\n", - " 0.1189612 , 0. , -0.02393692, 0. ,\n", - " 0.03479244, -0. , -0. , -0.09872656,\n", - " 0. , 0.03525269, 0.05757198, -0.00105693,\n", - " 0. , 0. , 0.04404994, 0. ,\n", - " -0. , 0.07952407, -0.02517718, 0.03949109,\n", - " 0. , -0. , 0.1057609 , 0.05089853,\n", - " 0. , -0.00070022, 0.01298889, -0. ,\n", - " -0.12222556, -0. , -0.03469444, 0.10075818,\n", - " 0. , 0.11131384, -0.03292956, 0.01280573,\n", - " -0. , 0. , -0. , 0. ,\n", - " 0. , -0. , 0. , -0. ,\n", - " 0.10009447, 0.01522033, 0. , 0.01894175,\n", - " -0.03899706, -0. , -0. , 0. ,\n", - " 0.07257131]]])" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", - " \n", - " \n", + "from fastcore.basics import patch\n", + "import uuid\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", " \n", "class ChewC:\n", " def __init___(self):\n", @@ -189,19 +65,28 @@ " return np.sum(np.array([self @ ind + self.intercept for ind in other.individuals]), axis=(1,2))\n", " \n", "class Population:\n", - " def __init__(self, genome, size):\n", + " def __init__(self, genome, size=None):\n", " self.genome = genome\n", - " self.size = size\n", + "# self.size = size\n", " self.ploidy = 2\n", - " self.individuals = self._create_initial_population()\n", - " self.chewc = None\n", + " if size:\n", + " self.size = size\n", + " self.individuals = self._create_initial_population()\n", + " else:\n", + " pass\n", "\n", " def _create_initial_population(self):\n", " \"\"\"Create an initial population of founder individuals.\"\"\"\n", - " return [Individual(self.genome, chewc=self) for _ in range(self.size)]\n", + " return [Individual(self.genome) for _ in range(self.size)]\n", " \n", " def get_haplo(self):\n", " return np.array([x.haplotype for x in self.individuals])\n", + " \n", + " def trial(self, h2):\n", + " [x.phenotype(h2) for x in self.individuals]\n", + "\n", + " def __getitem__(self, index):\n", + " return self.individuals[index]\n", "\n", " def __repr__(self):\n", " return f'Population of size: {self.size}'\n", @@ -222,7 +107,6 @@ " scaling_factors = np.sqrt(self.target_variance / founder_var)\n", " random_effects *= scaling_factors\n", " self.effects = random_effects\n", - " \n", " self.intercept = target_mean - founder_mean\n", " \n", " def __matmul__(self,other):\n", @@ -238,8 +122,8 @@ " self.source = source # 'founder', 'cross', 'self', 'dh'\n", " self.descendents = []\n", " self.fitness = 0\n", - "\n", " \n", + "\n", " # Logic for setting haplotype,mother,father\n", " if self.source == 'founder':\n", " self.haplotype= self._generate_random_haplotype()\n", @@ -258,116 +142,61 @@ " \n", " def __repr__(self):\n", " return f'Individual with haplotype shape: {self.haplotype.shape}'\n", + " \n", + " def make_id(self):self.id = uuid.uuid4().hex\n", + "\n", "\n", " def _generate_random_haplotype(self):\n", " \"\"\"Generate a random haplotype for the individual.\"\"\"\n", " return np.random.choice([0, 1], size=(self.genome.ploidy, self.genome.n_chr, self.genome.n_loci))\n", " \n", - " def __add__(self, h2):\n", - " self.fitness = self.phenotype + h2\n", - " \n", - " def gametes(self):\n", - " haplotypes = self.haplotype\n", - " def shuffle_chr(chromosome_pair):\n", - " \"\"\"\n", - " Perform crossover on a pair of chromosomes.\n", - "\n", - " Parameters:\n", - " chromosome_pair (list): A list of two equal length lists representing chromosomes.\n", - "\n", - " Returns:\n", - " list: A new chromosome formed by shuffling the given chromosome pair.\n", - " \"\"\"\n", - " # Ensure the chromosome pair contains two chromosomes of equal length\n", - " assert len(chromosome_pair) == 2\n", - " assert len(chromosome_pair[0]) == len(chromosome_pair[1])\n", - "\n", - " # Number of crossover points, sampled from a Poisson distribution with λ=1.3\n", - " n_crossover = np.random.poisson(1.3)\n", - "\n", - " # Determine crossover locations, sampled without replacement from chromosome length\n", - " chromosome_length = len(chromosome_pair[0])\n", - " crossover_locs = np.sort(np.random.choice(chromosome_length, n_crossover, replace=False))\n", - "\n", - " # Initialize the new chromosome and set the current chromosome to the first one\n", - " new_chromosome = []\n", - " current_chr = 0\n", - "\n", - " # Perform crossover by alternating segments between the two chromosomes\n", - " last_loc = 0\n", - " for loc in crossover_locs:\n", - " new_chromosome.extend(chromosome_pair[current_chr][last_loc:loc])\n", - " current_chr = 1 - current_chr # Switch to the other chromosome\n", - " last_loc = loc\n", - "\n", - " # Append the remaining segment\n", - " new_chromosome.extend(chromosome_pair[current_chr][last_loc:])\n", - " return np.array(new_chromosome)\n", - "\n", - "\n", - " # Initialize an empty array to store the shuffled chromosomes\n", - " shuffled_haplotypes = np.zeros_like(haplotypes)\n", - " ploidy, n_chr, n_loci = chewc.genome.shape\n", - " # Iterate over each chromosome and apply the shuffle_chr function\n", - " for i in range(haplotypes.shape[1]): # Iterate over the chromosomes\n", - " chromosome_pair = haplotypes[:, i, :] # Extract the chromosome pair (2, 77)\n", - " shuffled_chromosome = shuffle_chr(chromosome_pair) # Shuffle the chromosome pair \n", - " shuffled_haplotypes[:, i, :] = np.array(shuffled_chromosome).reshape(1, n_loci) # Store the shuffled chromosome\n", - " return shuffled_haplotypes[0,:,:]\n", + " def __mul__(self, partner):\n", + " if isinstance(partner,Individual):\n", + " source = 'cross'\n", + " mother, father = self.gamete(), partner.gamete()\n", "\n", " def phenotype(self, h2):\n", " breeding_value = chewc.trait @ self\n", - " print(breeding_value.shape)\n", " genetic_variance = np.var(breeding_value, ddof=1)\n", - " return breeding_value\n", - " \n", + " environmental_variance = (1 - h2) / h2 * genetic_variance\n", + " phenotype_value = breeding_value + np.random.normal(0, np.sqrt(environmental_variance))\n", + " self.fitness = phenotype_value\n", + "\n", + "\n", " \n", "\n", - " # Example usage\n", + "#Define the Simulation Parameters\n", "g = Genome(3, 77)\n", - "#make population with 100 founder individuals\n", "population = Population(g, size=100)\n", - "#make a trait\n", "trait = Trait(g, population,0,1)\n", - "# how to do trait @ population to get the trait values for a population?\n", "\n", + "#Plug them into ChewC\n", "chewc = ChewC()\n", "chewc.trait = trait\n", "chewc.population = population\n", - "chewc.genome = g\n", - "trait @ chewc.population\n", - "chewc.population.individuals[0].phenotype(1)" + "chewc.genome = g" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02e6f280", + "metadata": {}, + "outputs": [], + "source": [ + "chewc.population.trial(h2=1)" ] }, { "cell_type": "code", "execution_count": null, - "id": "e4cb0c2a", + "id": "e47e36d7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([ 0.90058078, 1.30441654, 1.55869168, 1.61792792, 0.52265153,\n", - " 1.54847419, 0.6922385 , 4.24468214, 0.3032351 , 1.2531959 ,\n", - " -0.37246337, 1.12530289, 0.07714408, 0.96496298, 0.65905819,\n", - " -1.15774067, -0.23176033, 0.40473708, -0.58547203, 0.91608648,\n", - " 1.74541433, 1.68537778, 0.02966238, 1.66886463, 0.12316054,\n", - " 1.04415549, 1.47553108, 1.75153837, 0.44607381, 0.79847286,\n", - " 2.85142217, 0.70347733, 1.04798263, 1.08622491, 2.73747453,\n", - " 0.05264775, -0.16099917, 1.26926665, 2.64583469, -0.43661853,\n", - " 0.05883183, 0.39777427, 1.12004998, -0.69285545, 0.23794893,\n", - " -0.48797881, -0.19745621, 1.88694249, 0.89300059, 0.41219806,\n", - " 1.5067358 , 0.8994578 , 1.86402182, 1.77564769, -0.08073889,\n", - " 0.19697015, 2.08561036, 0.3340788 , -1.00671384, -0.39244753,\n", - " 0.7540654 , 1.8896223 , -1.01793872, 0.83940053, 0.56879004,\n", - " -0.98796068, 0.99136101, 0.94027401, 0.33533645, 0.49182625,\n", - " 2.99311411, 2.02599197, -0.02020336, 2.63750165, 1.90660723,\n", - " 1.51705176, 1.12961731, 1.90010729, 0.7215806 , 0.27591076,\n", - " -0.5051664 , 1.52952333, 1.442353 , -0.09284849, 1.97394406,\n", - " 1.68704844, 1.53499888, -0.95910299, 0.1874767 , 1.01349149,\n", - " 1.65511075, 0.70985102, -0.16009287, 1.41368641, 1.98022853,\n", - " 0.57819343, 1.20019563, -0.62919673, -0.0437496 , 2.02325201])" + "-0.6592775099412544" ] }, "execution_count": null, @@ -376,13 +205,119 @@ } ], "source": [ - "chewc.trait @ chewc.population" + "chewc.population[9].fitness.sum()" ] }, { "cell_type": "code", "execution_count": null, - "id": "a0c7b773", + "id": "113ac472", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e833ae1", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from fastcore.basics import patch\n", + " \n", + "@patch\n", + "def gamete(self:Individual):\n", + " haplotypes = self.haplotype\n", + " def shuffle_chr(chromosome_pair):\n", + " \"\"\"\n", + " Perform crossover on a pair of chromosomes.\n", + "\n", + " Parameters:\n", + " chromosome_pair (list): A list of two equal length lists representing chromosomes.\n", + "\n", + " Returns:\n", + " list: A new chromosome formed by shuffling the given chromosome pair.\n", + " \"\"\"\n", + " # Ensure the chromosome pair contains two chromosomes of equal length\n", + " assert len(chromosome_pair) == 2\n", + " assert len(chromosome_pair[0]) == len(chromosome_pair[1])\n", + "\n", + " # Number of crossover points, sampled from a Poisson distribution with λ=1.3\n", + " n_crossover = np.random.poisson(1.3)\n", + "\n", + " # Determine crossover locations, sampled without replacement from chromosome length\n", + " chromosome_length = len(chromosome_pair[0])\n", + " crossover_locs = np.sort(np.random.choice(chromosome_length, n_crossover, replace=False))\n", + "\n", + " # Initialize the new chromosome and set the current chromosome to the first one\n", + " new_chromosome = []\n", + " current_chr = 0\n", + "\n", + " # Perform crossover by alternating segments between the two chromosomes\n", + " last_loc = 0\n", + " for loc in crossover_locs:\n", + " new_chromosome.extend(chromosome_pair[current_chr][last_loc:loc])\n", + " current_chr = 1 - current_chr # Switch to the other chromosome\n", + " last_loc = loc\n", + "\n", + " # Append the remaining segment\n", + " new_chromosome.extend(chromosome_pair[current_chr][last_loc:])\n", + " return np.array(new_chromosome)\n", + "\n", + "\n", + " # Initialize an empty array to store the shuffled chromosomes\n", + " shuffled_haplotypes = np.zeros_like(haplotypes)\n", + " ploidy, n_chr, n_loci = chewc.genome.shape\n", + " # Iterate over each chromosome and apply the shuffle_chr function\n", + " for i in range(haplotypes.shape[1]): # Iterate over the chromosomes\n", + " chromosome_pair = haplotypes[:, i, :] # Extract the chromosome pair (2, 77)\n", + " shuffled_chromosome = shuffle_chr(chromosome_pair) # Shuffle the chromosome pair \n", + " shuffled_haplotypes[:, i, :] = np.array(shuffled_chromosome).reshape(1, n_loci) # Store the shuffled chromosome\n", + " return shuffled_haplotypes[0,:,:]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c2eec44c", + "metadata": {}, + "outputs": [], + "source": [ + "chewc.population[0] * chewc.population[6]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0070cfa6", + "metadata": {}, + "outputs": [], + "source": [ + "new_pop = Population(chewc)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74cf1536", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "107594fa", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7702d741", "metadata": {}, "outputs": [], "source": []