1import functools
  2import numpy as np
  3
  4
  5from aegis_sim.parameterization import parametermanager
  6
  7
  8def stan_age(age):
  9    return age / parametermanager.parameters.AGE_LIMIT
 10
 11
 12class GPM_decoder:
 13    """Converter of simple user input into a genotype-phenotype map (as a list or a matrix).
 14
 15    --- Example input ---
 16
 17    PHENOMAP:
 18        "AP1, 7": # name of block, number of sites
 19            - ["surv", "agespec", 0.1] # trait, age-dependency, scale
 20            - ["repr", "agespec", -0.1]
 21        "MA1, 13":
 22            - ["surv", "agespec", -0.1]
 23
 24    # Make sure you set:
 25    GENARCH_TYPE: modifying
 26    BITS_PER_LOCUS: 1
 27
 28    # Make sure your genome size is big enough for all sites set up
 29    MODIF_GENOME_SIZE: 2000
 30
 31    ---
 32    """
 33
 34    def __init__(self, config_PHENOMAP):
 35        self.blocks = {}
 36        self.n = 0
 37
 38        for namen, quartets in config_PHENOMAP.items():
 39            # YML dict cannot take tuple as key, so it is fused as f'{name}, {n}'
 40            name, n = namen.split(",")
 41            n = int(n)
 42            self.__add_genblock(name, n)
 43            for trait, agefunc, funcparam in quartets:
 44                self.__add_encoding(name, trait, agefunc, funcparam)
 45
 46    def get_total_phenolist(self):
 47        phenolist = []
 48        for block in self.blocks.values():
 49            phenolist.extend(block.get_phenolist())
 50        return phenolist
 51
 52    def __len__(self):
 53        return self.n
 54
 55    # PRIVATE
 56    def __add_genblock(self, name, n):
 57        assert name not in self.blocks.keys(), f"Block with name {name} already exists."
 58
 59        genblock = Genblock(name=name, n=n, position=self.n)
 60        self.blocks[name] = genblock
 61        self.n += genblock.n
 62
 63    def __add_encoding(self, name, trait, agefunc, funcparam):
 64        self.blocks[name].add_encoding(trait, agefunc, funcparam)
 65
 66
 67class Genblock:
 68    def __init__(self, name, n, position):
 69        self.name = name
 70        self.n = n
 71
 72        self.position = position
 73
 74        self.encodings = []
 75
 76        self.phenolist = []
 77
 78    def add_encoding(self, trait, agefunc, funcparam):
 79        encoding = {
 80            "trait": trait,
 81            "agefunc": agefunc,
 82            "funcparam": funcparam,
 83        }
 84        self.encodings.append(encoding)
 85        self.__decode(encoding=encoding)
 86
 87    def get_phenolist(self):
 88        for index, trait, age, magnitude in self.phenolist:
 89            yield index + self.position, trait, age, magnitude
 90
 91    def __decode(self, encoding):
 92
 93        def add_to_phenolist(index, trait, age, magnitude):
 94            self.phenolist.append([index, trait, age, magnitude])
 95
 96        for i in range(self.n):
 97            func = Genblock.__resolve_function(encoding["agefunc"], encoding["funcparam"])
 98            if encoding["agefunc"] == "agespec":
 99                age = np.random.randint(0, parametermanager.parameters.AGE_LIMIT)
100                magnitude = func(age)
101                add_to_phenolist(i, encoding["trait"], age, magnitude)
102            else:
103                for age in range(parametermanager.parameters.AGE_LIMIT):
104                    magnitude = func(age)
105                    add_to_phenolist(i, encoding["trait"], age, magnitude)
106
107    @staticmethod
108    def __resolve_function(func, params):
109
110        scale = params
111
112        # PARAMETER DISTRIBUTIONS
113        def _beta_symmetrical(a=3, b=3):
114            """
115            a=3, b=3 -> bell-shape
116            a=1, b=3 -> decay (highest around 0, decreases with distance from 0)
117            a=1, b=1 -> flat
118            """
119            return np.random.beta(a, b) * 2 - 1
120
121        def _beta_onesided(a, b):
122            return np.random.beta(a, b)
123
124        # PHENOTYPE(AGE) functions
125        # These functions are used to compute phenotypic effects for every age
126        # They are going to be partial-ed with beta-shaped parameters (above)
127        # And then they take in age as a variable
128
129        def acc(intercept, slope, acc, age):
130            x = stan_age(age)
131            y = intercept + slope * x + acc * x**3
132            return scale * y
133
134        def exp(exp, age):
135            x = stan_age(age)
136            sign = [-1, 1][exp > 0]
137            y = sign * ((x + 1) ** (abs(exp) * 5) - 1)
138            return scale * y
139
140        def hump(amplitude, x0, width, age):
141            x = stan_age(age)
142            midpoint = 0.5
143            y = amplitude * np.exp(-((x - midpoint - x0) ** 2) / (0.02 + (width * 0.2) ** 2))
144            return scale * y
145
146        def sigm(slope, shift, age):
147            x = stan_age(age)
148            k = slope * 50
149            midpoint = 0.5
150            y = 1 / (1 + np.exp((-1) * k * (x - midpoint - shift / 2)))
151            return scale * y - scale / 2
152
153        # Assembling it together
154        if func == "agespec":
155            return functools.partial(acc, _beta_onesided(a=1, b=3), 0, 0)
156        elif func == "lin":
157            # linearly changing with age
158            return functools.partial(acc, _beta_symmetrical() / 2, _beta_symmetrical(), 0)
159        elif func == "acc":
160            # exponentially changing with age
161            return functools.partial(acc, _beta_symmetrical(), 0, _beta_symmetrical())
162        elif func == "exp":
163            # exponentially changing with age
164            return functools.partial(exp, _beta_symmetrical())
165        elif func == "hump":
166            # not age-specific but with an age-hump
167            return functools.partial(hump, _beta_symmetrical(), _beta_symmetrical(1, 1), _beta_symmetrical())
168        elif func == "sigm":
169            # sigmoidal across age
170            return functools.partial(sigm, _beta_symmetrical(), _beta_symmetrical(1, 1))
171        else:
172            raise Exception(f"{func} not allowed.")