def fail_simulation(n,s):
    t = 0
    r = 0
    trep = 2**64-1 # officialy INF
    t_list = []
    for _ in xrange(n):
        t_list.append(random.expovariate(1))

    t_list.sort()

    while True:
        # t1 < tmax
        t1 = t_list[0]
        
        if t1 < trep:
            t = t1
            r += 1
            if r == s+1:
                return t
            elif r < s+1:
                x = random.expovariate(1)
                t_list.pop(0)
                t_list.append(t+x)
                t_list.sort()
            if r == 1:
                y = random.expovariate(8)
                trep = t + y
        else:
            t = trep
            r -= 1
            if r > 0:
                y = random.expovariate(8)
                trep = t + y
            elif r == 0:
                trep = 2**64-1

def bankQueue(timeMax = 600): # t in minutes
    ta = 0       # time to next arrival in minutes
    ts = 0       # time to end of next service in minutes
    q = 1        # number of people in the queue
    c = 0        # time on the clock in minutes
    nserved = 0  # total number of people served
    cs = []
    ss = []
    qs = []
    while c < timeMax:
        if q == 0:
            c += ta
            q += 1
            ta = random.expovariate(1./avgArrive)
        elif ta<ts:
            ts -= ta
            c += ta
            q += 1
            ta = random.expovariate(1./avgArrive)
        else:
            ta -= ts
            c += ts
            q -= 1
            ts = random.expovariate(1./avgService)
            nserved += 1
        cs.append(c/60.)
        ss.append(nserved)
        qs.append(q)
    # end while
    p = biggles.FramedPlot()
    c1 = biggles.Curve(cs,ss,color='red')
    c2 = biggles.Curve(cs,qs,color='blue')
    p.add(c1,c2)
    p.show()

 def test_all(self):
     for test in ['', 'a', 'b', 'abc', 'abc'*50, 'hello world']:
         #print test
         #print sha256.sha256(test).hexdigest()
         #print hashlib.sha256(test).hexdigest()
         #print
         assert sha256.sha256(test).hexdigest() == hashlib.sha256(test).hexdigest()
     def random_str(l):
         return ''.join(chr(random.randrange(256)) for i in xrange(l))
     for length in xrange(150):
         test = random_str(length)
         a = sha256.sha256(test).hexdigest()
         b = hashlib.sha256(test).hexdigest()
         #print length, a, b
         if a != b:
             print 'ERROR!'
             raise ValueError()
     for i in xrange(100):
         test = random_str(int(random.expovariate(1/100)))
         test2 = random_str(int(random.expovariate(1/100)))
         
         a = sha256.sha256(test)
         a = a.copy()
         a.update(test2)
         a = a.hexdigest()
         
         b = hashlib.sha256(test)
         b = b.copy()
         b.update(test2)
         b = b.hexdigest()
         #print a, b
         if a != b:
             print 'ERROR!'
             raise ValueError()

def getFalseValue(provided_false_values, bin, lambda_bin_index):
    '''to pick a false value among the set
    The probability to select a false value, already provided, is higher than selecting a new one.
    provided_false_values is used for this reason
    '''
    value = 'null'
    if (len(set(provided_false_values)) < max_false_domain_cardinality) and (random.random() < 0.4): #selection not in VALUES ALREADY PROVIDED
        bin_index = int(random.expovariate(lambda_bin_index))
        while bin_index < 0 or bin_index >= len(bin):
            bin_index = int(random.expovariate(lambda_bin_index))
        if len(bin[bin_index])==0:
            return None
        if (len(bin[bin_index]) == 1):
            value_index = 0
        else:
            value_index = int(random.randint(0, len(bin[bin_index])-1))#Return a random integer N such that a <= N <= b.
        value = (bin[bin_index])[value_index]
    else: #selection among values already provided
        if (len(provided_false_values) > 0):
            if (len(provided_false_values) == 1):
                value_index = 0
            else:
                value_index = random.randint(0, len(provided_false_values)-1)

            value = provided_false_values[value_index]

    provided_false_values.append(value)
    return [value, provided_false_values]

 def get_interstorm_event_duration(self):
     """ This method is the interstorm duration generator
 
 This method takes one argument, the mean_interstorm parameter.
 (In Eagleson (1978), this parameter was called Tb.)
 
 This method is modeled identically to get_precipitation_event_duration()
 
 This method finds a random value for interstorm_duration
 based on the poisson distribution about the mean.
 This is accomplished using the expovariate function
 from the "random" standard library.
 Additionally, it is rounded to contain 4 significant figures, for neatness. 
 
 The if-else statement is very important here. Values of 0
 can exist in the Poission distribution, but it does not make
 sense to have 0 hour interstorm durations.
 To avoid 0 hour interstorm durations, we return a 
 interstorm duration IF it is greater than 0,
 otherwise, recursion is employed to re-call the interstorm
 duration generator function and get a new value.
 
 :returns: interstorm_duration as a float"""
 
     interstorm = round(random.expovariate(1/self.mean_interstorm),2)
     while interstorm == 0:
         interstorm = round(random.expovariate(1/self.mean_interstorm),2)
     self.interstorm_duration = interstorm
     return self.interstorm_duration

def generate_site_mutations(tree, position, mu, site_table, mutation_table,
                            multiple_per_node=True):
    """
    Generates mutations for the site at the specified position on the specified
    tree. Mutations happen at rate mu along each branch. The site and mutation
    information are recorded in the specified tables.  Note that this records
    more than one mutation per edge.
    """
    assert tree.interval[0] <= position < tree.interval[1]
    states = {"A", "C", "G", "T"}
    state = random.choice(list(states))
    site_table.add_row(position, state)
    site = site_table.num_rows - 1
    stack = [(tree.root, state, msprime.NULL_MUTATION)]
    while len(stack) != 0:
        u, state, parent = stack.pop()
        if u != tree.root:
            branch_length = tree.branch_length(u)
            x = random.expovariate(mu)
            new_state = state
            while x < branch_length:
                new_state = random.choice(list(states - set(state)))
                if multiple_per_node and (state != new_state):
                    mutation_table.add_row(site, u, new_state, parent)
                    parent = mutation_table.num_rows - 1
                    state = new_state
                x += random.expovariate(mu)
            else:
                if (not multiple_per_node) and (state != new_state):
                    mutation_table.add_row(site, u, new_state, parent)
                    parent = mutation_table.num_rows - 1
                    state = new_state
        stack.extend(reversed([(v, state, parent) for v in tree.children(u)]))

 def zombie_interval(self):
     if (datetime.now() - self.start_time).total_seconds() > 200:
         interval = random.expovariate(1/250)
     else:
         interval = 1000 - 50*self.kills
         interval = random.expovariate(1/max(interval, 500))
     self.timer.setInterval(interval)

    def generate(self, service):
        index = 0
        global buffer_size 
        global dropped_packets
        buffer_size = 10
        print ('Buffer Size is %d' %(buffer_size))
        while index < 6 :
                for x in range(1, total_packets):
			        arrivalTime = random.expovariate(lambda_list[index])
			        packet = Packet(env, arrivalTime, x, service)
			        yield env.timeout(arrivalTime)
                print ('Lambda: %f : Packet Loss Probability: %f'  %(lambda_list[index],dropped_packets *1.0/total_packets))
                dropped_packets = 0
                index = index + 1
        index = 0
        buffer_size = 50
        print ('Buffer Size is %d' %(buffer_size))
        while index < 6 :
                for x in range(1, total_packets):
                        arrivalTime = random.expovariate(lambda_list[index])
                        packet = Packet(env, arrivalTime, x, service)
                        yield env.timeout(arrivalTime)
                print ('Lambda: %f : Packet Loss Probability: %f'  %(lambda_list[index],dropped_packets *1.0/total_packets))
                dropped_packets = 0
                index = index + 1

    def simulate_once(self):
        """Simulate once.

        Return
        ------
            This will return the amount of operational loss simulated for the
            given time period.

        Example:
        r=OpRiskModel(stor4)
        lower, mu, upper = r.simulate_many()
        print "Bootstrap: ",lower,mu,upper

        Output:
        0.68% between 127760271.155 and 162467836.895
        0.8% between 122874286.419 and 167353821.63
        0.9% between 116569621.33 and 173658486.72
        0.95% between 111101264.604 and 179126843.445
        0.98% between 104743118.138 and 185484989.911
        0.99% between 100413671.581 and 189814436.469
        0.995% between 96401399.4999 and 193826708.549
        0.998% between 91486833.5654 and 198741274.484
        0.999% between 88010967.5982 and 202217140.451
        0.9999% between 77597567.1919 and 212630540.857
        0.99999% between 68459385.7079 and 221768722.341
        Bootstrap:  138714608.714 145114054.025 150873917.501 
        """
        t=random.expovariate(self.params.lamb)
        loss=0.0
        while t<self.params.days:
            t+=random.expovariate(self.params.lamb)
            amount=self.params.xm*random.paretovariate(self.params.alpha)
            loss+=amount
        return loss

def main():
    # initialize defaults
    meanoff = 1.0 # default mean off is 1 second
    meanon = 20.0 # default mean on is 20 seconds
    
    opts, args = getopt.getopt(sys.argv[1:], 'h', ['mean-off=', 'mean-on='])
    name = None
    for o, a in opts:
        if o == '-h': 
            print usage()
            sys.exit(0)
        elif o == '--mean-off': 
            meanoff = float(a)
        elif o == '--mean-on': 
            meanon = float(a)
        else: 
            print "Error parsing command-line option: unknown/unrecognized option %s" % o 
            #print >> sys.stderr, usage()
            sys.exit(1)
    if len(args):
        print "Error parsing command-line arguments: no arguments allowed" 
        #print >> sys.stderr, usage()
        sys.exit(2)
        
    seed()
    
    while 1:
        print time()
        system("cnistnet -a 10.2.34.150 192.168.5.10 --drop 0%; cnistnet -a 192.168.5.10 10.2.34.150 --drop 0%")
        sleep(expovariate(1.0/meanon))
        print time()
        system("cnistnet -a 10.2.34.150 192.168.5.10 --drop 100%; cnistnet -a 192.168.5.10 10.2.34.150 --drop 100%")
        sleep(expovariate(1.0/meanoff))

def req_generator(topology, n_contents, alpha, seed, rate=12.0, duration_warmup=9000, duration_real=36000):
    """This function generates events on the fly, i.e. instead of creating an 
    event schedule to be kept in memory, returns an iterator that generates
    events when needed.
    
    This is useful for running large schedules of events where RAM is limited
    as its memory impact is considerably lower.
    """
    random.seed(seed)
    
    rate = 4.0
    warmup = 0
    duration = 25000000
    
    receivers = [v for v in topology.nodes_iter() if topology.node[v]['stack'][0] == 'receiver']
    zipf = ZipfDistribution(alpha, n_contents)
    
    t_event = (expovariate(rate))
    while t_event < warmup + duration:
        recv = choice(receivers)
        content = int(zipf.rand_val())
        log = (t_event > warmup)
        event = {'receiver': recv, 'content': content, 'log': log}
        yield (t_event, event)
        t_event += (expovariate(rate))
    raise StopIteration()

    def random_gen(self, filepath="random.replay", avgTransmit=30.0, avgWait=0, totalTime=1800.0):
        f = open(filepath, "w")
        replay = list()
        count = 0
        port = 1234

        for cl in self.clients:
            curTime = 0
            port = 1234
            while True:
                if avgWait == 0:
                    startTime = 0
                else:
                    startTime = random.expovariate(1 / avgWait)

                txTime = random.expovariate(1 / avgTransmit)

                # Generate traffic until total time
                if curTime + startTime + txTime > totalTime:
                    break
                else:
                    replay.append((curTime + startTime, cl.IP() + ":" + str(port), curTime + startTime + txTime))
                    port += 1
                    curTime += startTime + txTime

        # Write replay to file sorted by timestamp
        replay.sort()
        count = 0
        port = 1234
        for startTime, dstIP, endTime in replay:
            srcIP = self.servers[count % len(self.servers)].IP() + ":" + str(port + count)
            f.write("%f %s %s %f\n" % (startTime, srcIP, dstIP, endTime))
            count += 1
        f.close()

def sample_rate_matrix(fs, Q_mut):
    nstates = len(Q_mut)
    # sample the selection parameters
    if fs.really_low_var:
        v = 0.04
    elif fs.low_var:
        v = 0.2
    elif fs.medium_var:
        v = 1
    elif fs.high_var:
        v = 5.0
    elif fs.really_high_var:
        v = 25.0
    s = math.sqrt(v)
    if fs.neg_skew:
        sels = [-random.expovariate(1/s) for i in range(nstates)]
    elif fs.no_skew:
        sels = [random.gauss(0, s) for i in range(nstates)]
    elif fs.pos_skew:
        sels = [random.expovariate(1/s) for i in range(nstates)]
    # define the mutation-selection rate matrix using Halpern-Bruno
    Q = np.zeros_like(Q_mut)
    for i in range(nstates):
        for j in range(nstates):
            if i != j:
                tau = math.exp(-(sels[j] - sels[i]))
                coeff = math.log(tau) / (1 - 1/tau)
                Q[i, j] = Q_mut[i, j] * coeff
    for i in range(nstates):
        Q[i, i] = -np.sum(Q[i])
    return Q

def poisson_train(start, end, rate):
    newspike = start+rnd.expovariate(rate)
    strain = []
    while newspike < end:
        strain.append(newspike)
        newspike += rnd.expovariate(rate)
    return strain

def similar_items(num_resources, min_fill, base_item, dev=.05, minr=MIN_RES, maxr=MAX_RES):
    """
    Generates a bin. Bin capacities are not correlated.
    All items are similar to the given base item
    base*dev is the standard deviation
    """
    mf = min(min_fill, 1-1e-15) # helps getting rid of numerical instabilities
    
    base = random.randint(minr,maxr)
    mean = base*dev
    lbd = 1.0/mean
    bin_cap = base_item.requirements[:]
    for i,v in enumerate(bin_cap):
        base = 5*v
        mean = base*dev
        lbd = 1.0/mean        
        bin_cap[i] = max(0,int(round(base+(random.expovariate(lbd)-mean))))
    
    b = Bin(bin_cap)
    item_vol = 0.0
    items = []
    vl = 0
    while item_vol < mf:
        for tr in xrange(MAX_TRY):
            it_res = base_item.requirements[:]
            for i,base in enumerate(it_res):
                mean = base*dev
                lbd = 1.0/mean
                it_res[i] = max(0,int(round(base+(random.expovariate(lbd)-mean))))
            vl = update(items, it_res, b, item_vol)
            if vl : break
        item_vol = vl
        if not item_vol: break

    return items, b

def gen_restock_list(cat):
    np = NeoPage(base_url='https://items.jellyneo.net')
    path = '/search/'
    args = []
    args.append(f'cat[]={cat}')

    args.append('min_rarity=1')
    args.append('max_rarity=100')
    args.append('status=1')

    args.append('sort=5')
    args.append('sort_dir=desc')

    start = 0
    while True:
        np.get(path, *args, f'start={start}')
        time.sleep(min(60, random.expovariate(30)) + 60)
        last_page = not re.search(r'<li class="arrow"><a href=".*?">&raquo;</a>', np.content)
        referer = np.referer
        results = re_jn_item.findall(np.content)
        for item_id, name, price_updated, price in results:
            if price == 'Inflation Notice':
                np.set_referer(referer)
                np.get(f'/item/{item_id}')
                time.sleep(min(60, random.expovariate(30)) + 60)
                price = re.search('<div class="price-row">(.*?) NP', np.content)[1]
            try:
                price = lib.amt(price)
            except (ValueError, TypeError, IndexError):
                # safe assumption that unpriceable items are at least 10M NP?
                price = 10000001
            print(f'\'{name}\': {price},')
        start += 50
        if last_page: break

def correlated_capacities(num_resources, min_fill, dev=.05,
                          rem_cons=1.0, correlated_items = False, minr=MIN_RES, maxr=MAX_RES):
    """
    Generates a bin. Bin capacities are correlated,
    and if correlated_items is True item requirements are correlated
    Either Volume(items) > min_fill or a 0 weighted item was generated
    base*dev is the standard deviation
    Item capacities are generated in [0,rem_cons*b.remaining]
    """
    mf = min(min_fill, 1-1e-15) # helps getting rid of numerical instabilities
    base = random.randint(minr,maxr)
    mean = base*dev
    lbd = 1.0/mean
    bin_cap = [max(0,int(round(base+(random.expovariate(lbd)-mean))))
               for x in xrange(num_resources)]
    b = Bin(bin_cap)
    item_vol = 0.0
    items = []
    while item_vol < mf:
        if correlated_items:
            for tr in xrange(MAX_TRY):
                base = random.randint(1,int(rem_cons*min(bin_cap)))
                mean = base*dev
                lbd = 1.0/mean
                it_res = [max(0,int(round(base+(random.expovariate(lbd)-mean))))
                      for x in xrange(num_resources)]
                vl = update(items, it_res, b, item_vol)
                if vl: break
            item_vol = vl
        else:
            it_res = [random.randint(0,int(rem_cons*x)) for x in b.remaining]
            item_vol = update(items, it_res, b, item_vol)
        if not item_vol: break

    return items, b

 def helper(size, mononomial):
     if size <= 1:
         if random.randint(0, 1):
             var = random.randrange(0, num_vars)
             return (OpTag.variable, var)
         k = random.expovariate(mean_constant)
         if random.randint(0, 1):
             k = int(k)
         return k
     tag = random_tag(weights, mononomial)
     if tag == OpTag.addition or tag == OpTag.product:
         a, b = random_split(size - 1)
         l = helper(a, mononomial)
         r = helper(b, mononomial)
         return (tag, l, r)
     if tag == OpTag.power:
         while True:
             k = random.expovariate(mean_exponent)
             if (not generalized) or random.randint(0, 1):
                 k = int(k + 2)
             if k < 40:
                 break
         x = helper(size - 1, mononomial)
         return (OpTag.power, x, k)
     if tag == OpTag.ratio:
         a, b = random_split(size - 1)
         l = helper(a, mononomial)
         r = helper(b, True)
         return (OpTag.ratio, l, r)
     assert(False)

	def infect(self,t,mu_input,sigma_e,dur_p,Dmax,D50,Dk,dur_d,rho,beta_p,betamax,beta50,betak,beta_d):
		self.Infected=True
		self.mu=mu_input
		self.timeInfection=t
		print("timeInfection: "+str(t))
		print("mu: "+str(self.mu))
# 	def setViralLoad(self,sigma_e):
		self.ViralLoad= random.normalvariate(self.mu, sigma_e)		
		print("ViralLoad: "+str(self.ViralLoad))
# 	def setTimes(self,dur_p,Dmax,D50,Dk,dur_d,rho):
		self.timeAsymptomatic=nextStage(self.timeInfection,self.timeDeath,random.expovariate(1.0/dur_p))
		print("timeAsymptomatic: "+str(self.timeAsymptomatic))
		
		expectedDurationAsymptomatic=Dmax*pow(D50,Dk)/(pow(10,self.ViralLoad*Dk) + pow(D50,Dk))
		print("expectedDurationAsymptomatic: "+str(expectedDurationAsymptomatic))
		actualDurationAsymptomatic=expectedDurationAsymptomatic*pow(-math.log(random.uniform(0,1)),1.0/rho)/math.gamma(1+1.0/rho)
		
		self.timeDisease=nextStage(self.timeAsymptomatic,self.timeDeath,actualDurationAsymptomatic)
		print("timeDisease: "+str(self.timeDisease))	
			
		self.timeDeathAIDS=nextStage(self.timeDisease,self.timeDeath,random.expovariate(1.0/dur_d))
		if not(self.timeDeathAIDS is None) and self.timeDeathAIDS < self.timeDeath:
			self.timeDeath = self.timeDeathAIDS
		print("timeDeath: "+ str(self.timeDeath))
# 	def setInfectivity(self,beta_p,betamax,beta50,betak,beta_d):
		self.infectivityPrimary= beta_p
		self.infectivityAsymptomatic= betamax*pow(10,self.ViralLoad*betak)/(pow(10,self.ViralLoad*betak)+pow(beta50,betak))
		self.infectivityDisease= beta_d		

    def __init__(self, space_width, space_height, animal_count, pred_count, genomes):
        self.pool = Pool()
        self.manager = Manager()

        self.width = space_width
        self.height = space_height
        self.genomes = genomes
        self.predators = [
            Predator(
                random.randrange(0, self.width, 1),
                random.randrange(0, self.height, 1)
            ) for _ in range(pred_count)
        ]
        self.animals = [ Animal(random.choice(genomes), PHEROMONES) for _ in range(animal_count) ]
        self.objects = [
            Food(
                random.randrange(0, self.width, 1),
                random.randrange(0, self.height, 1),
                10
            ) for _ in range(animal_count)
        ]
        for a in self.animals:
            a.teleport(
                random.randrange(0, self.width, 1),
                random.randrange(0, self.height, 1),
                random.uniform(0, 6.28)
            )
        self.next_food = random.expovariate(1.0/FOOD_PERIOD)
        self.next_poison = random.expovariate(1.0/POISON_PERIOD)
        self.pheromones = self.manager.list([])
        self.next_breed = random.expovariate(1.0/BREEDING_PERIOD)
        self.new_genomes = []