text:
########################################################################################
# Author: Sandhya Saisubramanian
# Description: Implements environment shaping and actor-designer coordination
# 			   for driving domain (multiple actor setting)
########################################################################################
import numpy as np
import sys
import os
import random
import time

from driving import Driving
from designer import Designer
from domain_helper import *
from sklearn.ensemble import RandomForestClassifier

current_file_path = os.path.dirname(os.path.realpath(__file__))
sys.path.append(os.path.join(current_file_path, '..'))
DRIVING_MAP_PATH = os.path.join(current_file_path, '..', 'maps','driving')

random.seed(100)

def Predict(x_train,y_train,x_test):
	final_model = RandomForestClassifier(n_estimators=100)
	x=np.array(x_train)
	y=np.array(y_train)
	final_model.fit(x_train,y_train)

	test_label = final_model.predict(x_test)

	return test_label

def actor_final_NSE_driving(map, start, goal,policy):
	dr = Driving(map,start,goal)
	return dr.calculate_final_cost_NSE(policy)

def actor_driving(map, start, goal,trajectory_budget=-1):
	dr = Driving(map,start,goal)
	if trajectory_budget > -1:
		return dr.solve_simulate(trajectory_budget)
	return dr.solve()

def solve_disapproved_actions(filename,start,goal,disapproved_actions={}):
	dr = Driving(filename,start,goal)
	return dr.solve_feedback(disapproved_actions)

def getE0():
	E0 = os.path.join(DRIVING_MAP_PATH,'grid-1.dr')
	return E0

def setup_driving():
	start_locations = []
	goal_locations = []
	maps = [os.path.join(DRIVING_MAP_PATH,f) for f in os.listdir(DRIVING_MAP_PATH) if os.path.isfile(os.path.join(DRIVING_MAP_PATH,f))]
	E0 = os.path.join(DRIVING_MAP_PATH,'grid-1.dr')
	modifications = {
			'null':  E0,
			'reduce_speed_all': os.path.join(DRIVING_MAP_PATH,'grid-1_reduced_speed.dr'),
			'fill_potholes': os.path.join(DRIVING_MAP_PATH,'grid-1_fill.dr'),
			'fill_deep_potholes': os.path.join(DRIVING_MAP_PATH,'grid-1_fill_deep.dr'),
			'reduce_speed_fill_deep': os.path.join(DRIVING_MAP_PATH,'grid-1_fill_reduce.dr'), #reduces speed at all shallow potholes and fills deep potholes
			'reduce_speed_zone1': os.path.join(DRIVING_MAP_PATH,'grid-1_reduced_speed1.dr'),
			'reduce_speed_zone2': os.path.join(DRIVING_MAP_PATH,'grid-1_reduced_speed2.dr'),
			'reduce_speed_zone3': os.path.join(DRIVING_MAP_PATH,'grid-1_reduced_speed3.dr'),
			'reduce_speed_zone4': os.path.join(DRIVING_MAP_PATH,'grid-1_reduced_speed4.dr')
			}

	# Cost of each modification proportional to pothole area in E0
	# Alternatively, this can be defined as per unit cost and
	# the exact value can be extracted for each configuration.
	cost_modifications = {
		'null': 0,
		'reduce_speed_zone1': 10,
		'reduce_speed_zone2': 16,
		'reduce_speed_zone3': 30,
		'reduce_speed_zone4': 28,
		'reduce_speed_all': 84, #cost/unit=2
		'fill_potholes': 168, #cost/unit=4
		'fill_deep_potholes': 76,
		'reduce_speed_fill_deep': 122
		}

	dummy_start = (1,1)
	dummy_goal = (13,22)
	dr = Driving(E0,dummy_start, dummy_goal)
	states = dr.getStates()

	for a in range(number_actors):
		s,g = get_start_goal_driving(states)
		start_locations.append(s)
		goal_locations.append(g)

	return maps, modifications, cost_modifications, E0, states, start_locations,goal_locations

def solve_driving(trajectory_budget_list, start_locations=[], goal_locations=[],cluster=False):
	k = 10000 # NSE penalty when the actor is unable to reach the goal
	trajectory_budget_modifications = []
	trajectory_budget_NSE = []
	trajectory_budget_costs = []

	maps, modifications, cost_modifications, E0, states,start_loc,goal_loc = setup_driving()
	if start_locations == []:
		start_locations = start_loc
		goal_locations = goal_loc

	for trajectory_budget in trajectory_budget_list:
		visited = []
		delta_actor_arr = []
		curr_policies = []
		current_NSE = []

		best_total_nse = 0
		violation = False
		best_modification = 'null'
		expected_costs = []
		nse_val = 0
		num_tested = 0

		designer = Designer(maps, modifications, cost_modifications, E0, delta_designer_percentage, domain_name, states)

		start = time.time()
		for actor in range(number_actors):
			policy,visitation_freq, expected_cost = actor_driving(E0, start_locations[actor], goal_locations[actor],trajectory_budget)
			curr_policies.append(policy)
			delta_actor = delta_actor_percentage * expected_cost
			delta_actor_arr.append(delta_actor)
			expected_costs.append(expected_cost)

			designer.populateParameters(policy,visitation_freq)
			nse_val += designer.getNSE(E0, policy)
			current_NSE.append(nse_val)

		# Computes modifications that do not require testing since they are too
		# similar to other modifications with better utility.
		if(cluster):
			print("Shaping budget = ",shaping_budget)
			visited = designer.cluster_best_design(curr_policies,shaping_budget,number_actors)

		start1 = time.time()

		delta_designer = delta_designer_percentage * nse_val
		if nse_val > delta_designer:
			violation = True
		best_total_nse = nse_val

		while violation:
			violation = False
			current_NSE = []
			utility, modification, updated_NSE_actors = designer.best_design_multiple_actors(visited, curr_policies)
			visited.append(modification)
			num_tested += 1
			if modification == 'null':
				break

			for actor in range(number_actors):
				updated_policy, visitation_freq, updated_cost = actor_driving(modifications[modification], start_locations[actor],\
																			 goal_locations[actor],trajectory_budget)

				if  updated_cost - expected_costs[actor] > delta_actor_arr[actor]:
					updated_NSE_actors[actor] = k
				else:
					policy = updated_policy

			if sum(updated_NSE_actors) > delta_designer:
				violation = True

			current_NSE = updated_NSE_actors

			if sum(current_NSE) < best_total_nse:
				best_modification = modification
				best_config = modifications[modification]
				best_NSE = current_NSE

		print("***************************************************************************")
		print("Time taken (s) =%s time taken (wo similarity calculation) = %s num_tested=%s\n"%((time.time() -start),(time.time()-start1),num_tested))
		print("Best modification = %s\n"%best_modification)
		expected_costs = []
		final_nse = []
		trajectory_budget_modifications.append(best_modification)
		for actor in range(number_actors):
			policy, expected_cost = actor_driving(modifications[best_modification], start_locations[actor], goal_locations[actor])
			expected_costs.append(expected_cost)
			avg_cost, std_cost, avg_nse, std_nse = actor_final_NSE_driving(modifications[best_modification],\
														start_locations[actor], goal_locations[actor],policy)
			final_nse.append(avg_nse)

		trajectory_budget_NSE.append(final_nse)
		trajectory_budget_costs.append(expected_costs)

	return trajectory_budget_NSE, trajectory_budget_costs, trajectory_budget_modifications, start_locations, goal_locations

def get_start_goal_driving(states):
	E0 = os.path.join(DRIVING_MAP_PATH,'grid-1.dr')
	grid = readmap(E0)
	walls = WallLoc_driving(grid)
	valid_state = False
	while(valid_state == False):
		start_state = random.choice(states)
		goal_state = random.choice(states)
		if start_state and goal_state not in walls:
			valid_state = True
	return start_state, goal_state

def noShaping_driving(start_locations, goal_locations):
	baseline_nse = []
	baseline_costs = []

	E0 = os.path.join(DRIVING_MAP_PATH,'grid-1.dr')
	for actor in range(number_actors):
		policy,expected_cost = actor_driving(E0, start_locations[actor], goal_locations[actor])
		avg_cost, std_cost, avg_nse, std_nse = actor_final_NSE_driving(E0, start_locations[actor], \
																		goal_locations[actor],policy)
		baseline_nse.append(avg_nse)
		baseline_costs.append(expected_cost)

	baseline_avg_nse = sum(baseline_nse)/number_actors
	baseline_std_nse = np.std(np.array(baseline_nse))
	baseline_avg_cost = sum(baseline_costs)/number_actors

	return baseline_avg_nse, baseline_std_nse, baseline_avg_cost

def DisapprovedActions(approval,all_states,sa_list,all_actions,generalize_feedback=False):
	x_train = []
	y_train = []
	x_test = []
	testing_sa = []
	disapproved_actions = {s:[] for s,state in enumerate(all_states)}
	# Based on gathered feedback
	for info in approval:
		temp=[]
		s = info[0]
		temp.append(s[0])
		temp.append(s[1])
		temp.append(all_actions.index(info[1]))
		x_train.append(temp)
		y_train.append(int(info[2]))

		state_index = all_states.index(info[0])
		if int(info[2]) == 0:
			disapproved_actions[state_index].append(info[1])

	# Generalize the gathered data to unseen states:
	if generalize_feedback == True:
		for sa in sa_list:
			sa_val = sa[1:]
			temp=[]
			if sa_val not in x_train:
				s = info[0]
				temp.append(s[0])
				temp.append(s[1])
				temp.append(all_actions.index(sa_val[1]))
				x_test.append(temp)
				testing_sa.append(sa_val)

		y_label = Predict(x_train, y_train, x_test)

		for i in range(len(y_label)):
			if y_label[i] == 0:
				sa_val = testing_sa[i]
				state_index = all_states.index(sa_val[0])
				disapproved_actions[state_index].append(sa_val[1])

	return disapproved_actions

def feedback(trajectory_budget_list,start_locations, goal_locations,generalize_feedback=False):
	feedback_budget = 500
	E0 = getE0()

	grid = readmap(E0)
	dr = Driving(E0,start_locations[0], goal_locations[0])
	all_states = dr.getStates()
	all_actions = dr.getActions()
	sa_list = dr.generate_state_actions()
	NSE_locations = PotholeLoc_driving(grid)
	NSE_values = []
	std_NSE = 0

	for trajectory_budget in trajectory_budget_list:
		total_nse = []
		for actor in range(number_actors):
			policy,visitation_freq, expected_cost = actor_driving(E0, start_locations[actor], goal_locations[actor],trajectory_budget)
			nse_penalty = NSE_penalty_driving(all_states, policy,NSE_locations,visitation_freq)
			delta_designer = delta_designer_percentage * nse_penalty
			feedback_count = 0

			approval = []
			if nse_penalty > delta_designer:
				for s,state in enumerate(all_states):
					if feedback_count >= feedback_budget:
						break
					if s in policy:
						temp = []
						action = policy[s]
						if mild_NSE_driving(state,action, NSE_locations) or severe_NSE_driving(state,action, NSE_locations):
							temp.append(state)
							temp.append(action)
							temp.append(0)
						else:
							temp.append(state)
							temp.append(action)
							temp.append(1)
						approval.append(temp)
						feedback_count += 1
				disapproved_actions = DisapprovedActions(approval,all_states,sa_list,all_actions,generalize_feedback)
				updated_policy,updated_expected_cost = solve_disapproved_actions(E0,start_locations[actor], \
																			goal_locations[actor],disapproved_actions)
				if updated_policy != None and updated_expected_cost - expected_cost <= delta_actor_percentage * expected_cost:
					avg_cost, std_cost, avg_nse, std_nse = actor_final_NSE_driving(E0, start_locations[actor], \
																			goal_locations[actor],updated_policy)

					total_nse.append(avg_nse)
				else:
					policy,expected_cost = actor_driving(E0, start_locations[actor], goal_locations[actor])
					avg_cost, std_cost, avg_nse, std_nse = actor_final_NSE_driving(E0, start_locations[actor], \
																			goal_locations[actor],policy)
					total_nse.append(avg_nse)
			else:
				policy,expected_cost = actor_driving(E0, start_locations[actor], goal_locations[actor])
				avg_cost, std_cost, avg_nse, std_nse = actor_final_NSE_driving(E0, start_locations[actor], \
																			goal_locations[actor],policy)

				total_nse.append(avg_nse)
			std_NSE = np.std(np.array(total_nse))
		NSE_values.append(sum(total_nse)/number_actors)

	return NSE_values, std_NSE

def main():
	avg_nse = []
	std_nse = []
	avg_cost = []

	shaping_budget_avg_nse = []
	shaping_budget_std_nse = []
	shaping_budget_avg_cost = []

	trajectory_budget_list = [100]

	op_file = "../results/multiple_actors_driving_trials.txt"
	print("Shaping with exhaustive search..\n")
	arr_NSE, arr_cost,best_modification, start_locations, goal_locations = solve_driving(trajectory_budget_list)

	print("Shaping with budget..\n")
	cluster_arr_NSE, cluster_arr_cost,cluster_best_modification, \
	start_locations, goal_locations = solve_driving(trajectory_budget_list,start_locations, goal_locations,True)

	baseline_avg_nse, baseline_std_nse, baseline_avg_cost = noShaping_driving(start_locations, goal_locations)

	feedback_nse, feedback_std_nse = feedback(trajectory_budget_list,start_locations, goal_locations)

	feedback_gen_nse, feedback_gen_std_nse = feedback(trajectory_budget_list,start_locations, goal_locations,True)

	for t in range(len(trajectory_budget_list)):
		avg_nse.append(sum(arr_NSE[t])/number_actors)
		std_nse.append(np.std(np.array(arr_NSE[t])))
		avg_cost.append(sum(arr_cost[t])/number_actors)

		shaping_budget_avg_nse.append(sum(cluster_arr_NSE[t])/number_actors)
		shaping_budget_std_nse.append(np.std(np.array(cluster_arr_NSE[t])))
		shaping_budget_avg_cost.append(sum(cluster_arr_cost[t])/number_actors)

	file = open(op_file,"a+")
	file.write("#Actors=%s\n"%number_actors)
	file.write("Designer slack percentage=%s\n"%delta_designer_percentage)
	file.write("Actor slack percentage=25\n")
	file.write("Baseline NSE =%s\n"%baseline_avg_nse)
	file.write("Baseline_std_NSE =%s\n"%baseline_std_nse)
	file.write("Baseline costs =%s\n"%baseline_avg_cost)

	file.write("Shaping_NSE=%s\n"%avg_nse)
	file.write("Shaping_std_NSE=%s\n"%std_nse)
	file.write("Shaping_actor_costs=%s\n"%avg_cost)

	file.write("Shaping_budget_NSE=%s\n"%shaping_budget_avg_nse)
	file.write("Shaping_budget_std_NSE=%s\n"%shaping_budget_std_nse)
	file.write("Shaping_budget_actor_costs=%s\n"%shaping_budget_avg_cost)

	file.write("Feedback_NSE=%s\n"%feedback_nse)
	file.write("Feedback_std_NSE=%s\n"%feedback_std_nse)
	file.write("Feedback_gen_NSE=%s\n"%feedback_gen_nse)
	file.write("Feedback_gen_std_NSE=%s\n"%feedback_gen_std_nse)
	file.write("***********************************\n")
	file.close()

if __name__ == '__main__':
	domain_name = "driving"
	shaping_budget = 4
	delta_designer_percentage = 0
	delta_actor_percentage = 0.25
	number_actors = 10
	if len(sys.argv) > 1:
		number_actors = int(sys.argv[1])
	main()


text:
import numpy as np
import scipy.signal

neighbors = np.full((3, 3), 1)
neighbors[1, 1] = 0

o = np.array(list(map(int, open("input-11.txt").read().replace("\n", "")))).reshape(10, 10)

total_flashes = 0
for step in range(1, 10000):
    o += 1

    flashed = np.full(o.shape, False)
    while True:
        flashed_now = (o > 9) & np.logical_not(flashed)
        if not flashed_now.any():
            break

        flashed |= flashed_now
        o += scipy.signal.convolve2d(flashed_now, neighbors, mode="same")

    total_flashes += flashed.sum()
    o[flashed] = 0

    if flashed.all():
        print("all flashed", step)
        break

    if step == 100:
        print(total_flashes)

ext: py

max_stars_repo_name: lkesteloot/advent-of-code

max_issues_repo_name: lkesteloot/advent-of-code

max_forks_repo_name: lkesteloot/advent-of-code

reason: import numpy,import scipy


text:
import math
import sys
from itertools import chain
from typing import Optional

import numpy as np
import torch
from torch import nn

from ifaces import DistributionSampler

class UnitGaussianSampler(nn.Module, DistributionSampler):
    def __init__(self, device: str = 'cpu'):
        nn.Module.__init__(self)
        self.device = device

    def forward(self, x):
        return self.sample(x)

    def sample(self, shape_or_x: tuple or torch.Tensor) -> torch.Tensor:
        if type(shape_or_x) == tuple:
            return torch.randn(shape_or_x).to(self.device)
        return torch.randn_like(shape_or_x)

    def log_likelihood(self, samples: torch.Tensor, x: torch.Tensor = None) -> torch.Tensor:
        """
        Given samples as rows of a matrix, returns their log-likelihood under the zero mean unit covariance
        Gaussian as a COLUMN vector.
        """
        return -0.5 * samples.shape[1] * math.log(2 * math.pi) - 0.5 * torch.sum(samples.square(), dim=1)

# noinspection PyMethodMayBeStatic
class Exponential(nn.Module):
    def __init__(self):
        super().__init__()

    def forward(self, x):
        return torch.exp(x)

class GaussianSampler(nn.Module, DistributionSampler):
    def __init__(self, h_network: nn.Module, mean_network: nn.Module, sigma_network: nn.Module):
        nn.Module.__init__(self)
        self.h_network = h_network
        self.mean_network = mean_network
        self.sigma_network = sigma_network
        self.prior = UnitGaussianSampler()

    def forward(self, x):
        return self.sample(x)

    def get_mean_sigma(self, x: torch.Tensor):
        """
        Returns the mean and the square root of the covariance of the Gaussian.
        :param torch.Tensor x: input to the input layer
        """
        h = self.h_network(x)
        mean = self.mean_network(h)
        sigma = self.sigma_network(h)
        return mean, sigma

    def get_mean(self, x: torch.Tensor):
        """
        :param torch.Tensor x: input to the input layer
        """
        h = self.h_network(x)
        return self.mean_network(h)

    def sample(self, shape_or_x: tuple or torch.Tensor):
        assert type(shape_or_x) == torch.Tensor
        mean, sigma = self.get_mean_sigma(shape_or_x)
        unit_gaussian_samples = torch.randn_like(mean)
        return sigma * unit_gaussian_samples + mean

    @staticmethod
    def log_likelihood_for_mean_sigma(samples: torch.Tensor, mean: torch.Tensor, sigma: torch.Tensor):
        return -0.5 * samples.shape[1] * math.log(2 * np.pi) \
               - 0.5 * torch.sum(((samples - mean) / (sigma + sys.float_info.epsilon)).square() +
                                 2.0 * torch.log(sigma + sys.float_info.epsilon), dim=1)

    def log_likelihood(self, samples: torch.Tensor, x: torch.Tensor):
        mean, sigma = self.get_mean_sigma(x)
        return GaussianSampler.log_likelihood_for_mean_sigma(samples, mean, sigma)

    def first_linear_layer_weights_np(self) -> np.ndarray:
        assert type(self.h_network[0]) == nn.Linear
        return self.h_network[0].weight.data.clone().detach().cpu().numpy().T

    @staticmethod
    @torch.no_grad()
    def init_weights(m):
        if type(m) == nn.Linear:
            torch.nn.init.xavier_uniform_(m.weight)
            m.bias.data.fill_(0.0)

    @staticmethod
    def random(n_units: list, mean: Optional[torch.Tensor or float] = None):
        h_network = nn.Sequential(*list(chain.from_iterable(
            [nn.Linear(in_features=h_in, out_features=h_out), nn.Tanh()]
            for h_in, h_out in zip(n_units[:-2], n_units[1:-1])
        )))
        h_network.apply(GaussianSampler.init_weights)
        mean_network = nn.Linear(n_units[-2], n_units[-1])
        mean_network.apply(GaussianSampler.init_weights)
        if mean is not None:
            with torch.no_grad():
                mean_network.bias.copy_(mean)
        sigma_network = nn.Sequential(
            nn.Linear(n_units[-2], n_units[-1]),
            Exponential(),
        )
        sigma_network.apply(GaussianSampler.init_weights)
        return GaussianSampler(h_network, mean_network, sigma_network)

max_stars_repo_name: NickSmyr/iwae-pytorch

max_issues_repo_name: NickSmyr/iwae-pytorch

max_forks_repo_name: NickSmyr/iwae-pytorch

reason: import numpy


text:
from numpy import uint64
from pyteal import *
from algosdk.future import transaction
from algosdk.v2client import algod
import os
import pathlib
from util import deploy, get_private_key_from_mnemonic

#if Resource: application_args = (user_address_bytes, asset_total, asset_unit_name, asset_name, asset_url, asset_metadata_hash, asset_manager_address)
# NOT TRUE Resource must include the address of the apat field of their no_op txn
#if User: application_args = ()
# User must include the ASA ID in the apas field of their no_op txn

# ASA Byteslice scheme: user_address[32B], asset_id[8B], round_created[8B]
ADDRESS_LEN = Int(32)
ASSET_ID_LEN = Int(8)
ROUND_LEN = Int(8)
FREE_AFTER_ROUNDS = Int(72)
class LocalState:
    SCHEMA: transaction.StateSchema = transaction.StateSchema(num_uints=0, num_byte_slices=0)

class GlobalState:
    SCHEMA: transaction.StateSchema = transaction.StateSchema(num_uints=0, num_byte_slices=64)
    class Variables:
        ASSET_PREFIX: TealType.bytes = Bytes("asset_")

def approval_program():

    asa_slot = ScratchVar(TealType.uint64)
    i = ScratchVar(TealType.uint64)
    curr_slot = ScratchVar(TealType.bytes)
    curr_i_round = ScratchVar(TealType.uint64)
    oldest_created = ScratchVar(TealType.uint64)

    create_asa = Seq([
        asa_slot.store(Int(64)),
        oldest_created.store(Global.round()),
        For(i.store(Int(0)), i.load() < Int(64), i.store(i.load() + Int(1))).Do(Seq([
            If(App.globalGet(Concat(GlobalState.Variables.ASSET_PREFIX, Itob(i.load()))) != Int(0),
                Seq([
                    curr_slot.store(App.globalGet(Concat(GlobalState.Variables.ASSET_PREFIX, Itob(i.load())))),
                    curr_i_round.store(Btoi(Extract(curr_slot.load(), ADDRESS_LEN + ASSET_ID_LEN, ROUND_LEN))),
                    If(And(curr_i_round.load() <= Global.round() - FREE_AFTER_ROUNDS, curr_i_round.load() < oldest_created.load()), Seq([
                        asa_slot.store(i.load()),
                        oldest_created.store(curr_i_round.load())
                    ]))
                ]),
                Seq([
                    asa_slot.store(i.load()),
                    Break()
                ])
            ),
        ])),
        If(asa_slot.load() == Int(64), Err()),
        Assert(Txn.application_args.length() == Int(7)),
        # Create an asset
        InnerTxnBuilder.Begin(),
        InnerTxnBuilder.SetFields({
            TxnField.type_enum: TxnType.AssetConfig,
            TxnField.config_asset_total: Btoi(Txn.application_args[1]),
            TxnField.config_asset_decimals: Int(0),
            TxnField.config_asset_unit_name: Txn.application_args[2],
            TxnField.config_asset_name: Txn.application_args[3],
            TxnField.config_asset_url: Txn.application_args[4],
            TxnField.config_asset_metadata_hash: Txn.application_args[5],
            TxnField.config_asset_manager: Txn.application_args[6],
            TxnField.config_asset_reserve: Txn.application_args[6],
            TxnField.config_asset_freeze: Txn.application_args[6],
            TxnField.config_asset_clawback: Txn.application_args[6]
        }),
        InnerTxnBuilder.Submit(),
        App.globalPut(Concat(GlobalState.Variables.ASSET_PREFIX, Itob(asa_slot.load())), Concat(Txn.application_args[0], Itob(InnerTxn.created_asset_id()), Itob(Global.round()))),
        Return(Int(1))
    ])

    get_user_asa_info = Seq([
        asa_slot.store(Int(64)),
        For(i.store(Int(0)), i.load() < Int(64), i.store(i.load() + Int(1))).Do(Seq([
            If(App.globalGet(Concat(GlobalState.Variables.ASSET_PREFIX, Itob(i.load()))) != Int(0),
                Seq([
                    curr_slot.store(App.globalGet(Concat(GlobalState.Variables.ASSET_PREFIX, Itob(i.load())))),
                    If(Extract(curr_slot.load(), Int(0), ADDRESS_LEN) == Txn.sender(), Seq([
                        asa_slot.store(i.load()),
                        Break()
                    ]))
                ])
            ),
        ])),
        If(asa_slot.load() == Int(64), Err())
    ])
    asset_balance = AssetHolding.balance(App.id(), Txn.assets[0])

    claim_asa = Seq([
        get_user_asa_info,
        asset_balance,
        Assert(asset_balance.hasValue()),
        # Send the assest to the specified address
        InnerTxnBuilder.Begin(),
        InnerTxnBuilder.SetFields({
            TxnField.type_enum: TxnType.AssetTransfer,
            TxnField.asset_receiver: Txn.sender(),
            TxnField.asset_amount: asset_balance.value(),
            # Must be in the assets array sent as part of the application call
            TxnField.xfer_asset: Btoi(Extract(App.globalGet(Concat(GlobalState.Variables.ASSET_PREFIX, Itob(asa_slot.load()))), ADDRESS_LEN, ASSET_ID_LEN)),
        }),
        InnerTxnBuilder.Submit(),
        Return(Int(1))
    ])

    handle_noop = If(Global.creator_address() ==
                     Txn.sender(), create_asa, claim_asa)

    # TODO: Can set a required fee for users to pay when they opt-in, could be used for one-time payments for a service
    handle_optin = Seq([
        Return(Int(1))
    ])

    handle_closeout = Seq([
        Return(Int(1))
    ])

    handle_updateapp = If(Global.creator_address() ==
                          Txn.sender(), Return(Int(1)), Err())

    handle_deleteapp = If(Global.creator_address() ==
                          Txn.sender(), Return(Int(1)), Err())

    program = Cond(
        [Txn.on_completion() == OnComplete.NoOp, handle_noop],
        [Txn.on_completion() == OnComplete.OptIn, handle_optin],
        [Txn.on_completion() == OnComplete.CloseOut, handle_closeout],
        [Txn.on_completion() == OnComplete.UpdateApplication, handle_updateapp],
        [Txn.on_completion() == OnComplete.DeleteApplication, handle_deleteapp]
    )
    return program

def clear_program():
    program = Return(Int(1))
    return program

def create(API_TOKEN, API_URL, MNEMONIC, upload=False):
    # Compile to TEAL
    with open(os.path.join(pathlib.Path(__file__).parent, 'blockin_global_approval.teal'), 'w') as f:
        approval_teal = compileTeal(approval_program(), Mode.Application, version=5)
        f.write(approval_teal)

    with open(os.path.join(pathlib.Path(__file__).parent, 'blockin_global_clear.teal'), 'w') as f:
        clear_teal = compileTeal(clear_program(), Mode.Application, version=5)
        f.write(clear_teal)

    if upload:

        # initialize an algodClient
        algod_client = algod.AlgodClient(API_TOKEN, API_URL, {
            "x-api-key": API_TOKEN})

        # Create & Deploy the Application
        return deploy(algod_client, get_private_key_from_mnemonic(MNEMONIC),
               approval_teal, clear_teal, GlobalState.SCHEMA, LocalState.SCHEMA)

max_stars_repo_name: matt-davison/blockin

max_issues_repo_name: matt-davison/blockin

max_forks_repo_name: matt-davison/blockin

reason: from numpy


text:
################################################################################
# infer segments from breakpoint pairs
################################################################################
#------------------------------------------------------------------------------#
import numpy as np
import pandas as pd
import math
import scipy.stats as stats
import statsmodels.api as sm
from itertools import groupby
from itertools import combinations
from operator import itemgetter
from collections import defaultdict
from statsmodels.formula.api import glm
import matplotlib.pyplot as plt

#------------------------------------------------------------------------------#
import cbs
import misc
from grange import GRange
from bpjoint import join_breakpoint

#------------------------------------------------------------------------------#
# get bias in binsize
#------------------------------------------------------------------------------#
def get_bin_bias(bin_count, biasfile, binsize = 10000):
    """calculate bias factor in binsize
    """
    df = pd.read_csv(biasfile, sep = '\t')
    df.columns = ['Chrom', 'Start', 'End', 'GCcontent', 'Mappability']
    df.Start = (np.floor(df.Start/binsize)*binsize).astype('int')
    df_agg = df.groupby(['Chrom', 'Start'])[['GCcontent', 'Mappability']]\
        .mean().reset_index().rename(columns = {"Start":"Coord"})
    dfm = pd.merge(bin_count, df_agg, on = ['Chrom', 'Coord'])
    return dfm

#------------------------------------------------------------------------------#
# glm remove GC and mappability biasa
#------------------------------------------------------------------------------#
def norm_glm(dfm, cutoff = [0, 0], method = 'nb', fold_bound = 5, plot = False):
    """fit glm to remove biases from bin counts

    Parameters
    ----------
    cutoff : select data points used in the model, default use all
    method : 'nb' - negative binomial
             'poisson' - poisson
    fold_bound : bound of bias factor to consider
        to avoid getting super small or large normalized counts
    """
    keep = (dfm.GCcontent >= cutoff[0]) & (dfm.Mappability >= cutoff[1])
    if (method == 'nb'):
        fit = glm('Count ~ GCcontent + Mappability', data = dfm,
            subset = keep, family = sm.families.NegativeBinomial()).fit()
    else:
        fit = glm('Count ~ GCcontent + Mappability', data = dfm,
            subset = keep, family = sm.families.Poisson()).fit()

    coef = fit.params
    fitted_y = np.exp(coef[0] + coef[1] * dfm.GCcontent +
                      coef[2] * dfm.Mappability)
    fac = fitted_y / np.median(fitted_y)

    fac[fac > fold_bound] = fold_bound
    fac[fac < 1/fold_bound] = 1/fold_bound

    norm_y = round(dfm.Count / fac, 4)

    bin_norm = dfm.iloc[:,0:4].copy()
    bin_norm['Count'] = norm_y + 1
    bin_norm['CN'] = 2 * bin_norm['Count'] / get_mode(bin_norm['Count'])

    if (plot):
        plt.plot(dfm.GCcontent, dfm.CN, 'o', alpha = 0.01)
        plt.yscale('log')
        plt.show()

        plt.plot(dfm.GCcontent, norm_y, 'o', alpha = 0.01)
        plt.yscale('log')
        plt.show()

        plt.plot(dfm.Mappability, dfm.CN, 'o', alpha = 0.01)
        plt.yscale('log')
        plt.show()

        plt.plot(dfm.Mappability, norm_y, 'o', alpha = 0.01)
        plt.yscale('log')
        plt.show()

    return bin_norm

#------------------------------------------------------------------------------#
# 2D bias plot
#------------------------------------------------------------------------------#
def plot_bias_2D(x, y, z, npts = 100, cmap = 'gist_rainbow',
                 xlab = 'GCcontent', ylab = 'Mappability'):
    """ 2D bias plot ('viridis')
    """

    # bin data and take average
    zi = stats.binned_statistic_2d(y, x, z, 'mean', bins=npts).statistic
    print(np.nansum(abs(zi-np.nanmedian(zi))))

    # make heatmap
    plt.imshow(zi, cmap=cmap, origin = 'lower', extent=(0,1,0,1))
    plt.colorbar()
    plt.xlabel(xlab)
    plt.ylabel(ylab)
    plt.show()

#------------------------------------------------------------------------------#
# apply triangular smoothing on the normalized bin counts
#------------------------------------------------------------------------------#
def smooth_cn(bin_norm, k = 10):
    """smooth normalized bin counts
    """
    df = bin_norm.copy()
    # calculate log2 ratio
    df['Log2Ratio'] = np.log2(np.array(df.CN / 2))

    df_smooth = []
    for chr in misc.unique(df['Chrom']):
        dfsub = df[df['Chrom'] == chr]
        x = np.array(dfsub['Log2Ratio'])
        x_smooth = trismooth_1d(x, k = k)
        dfsub = dfsub.assign(Log2Ratio = x_smooth)
        df_smooth.append(dfsub)
    return pd.concat(df_smooth)

#------------------------------------------------------------------------------#
def trismooth_1d(x, k = 10, outlier_sd = 5,
                 edge_repeat = False):
    """triangular smoothing on a vector, only smooth outliers
    k : # of points to consider on the left and the right of a point
    outlier_sd : mean +- outlier_sd * sd is the outlier in the region
    edge_repeat : False - boundaries extended by filling nan
        True - boundaries extended by repeating edge points
    """

    x = np.array(x)
    y = np.zeros((len(x), 2*k + 1), dtype = x.dtype)
    y[:, k] = x
    for i in range(k):
        j = k - i
        y[j:,i] = x[:-j]
        y[:-j,-(i+1)] = x[j:]
        if (edge_repeat):
            y[:j,i] = x[0]
            y[-j:,-(i+1)] = x[-1]
        else:
            y[:j,i] = np.nan
            y[-j:,-(i+1)] = np.nan

    # identify outlier
    ym = y[:,k]
    ys = y[:,[i for i in range(2*k+1) if i != k]]
    ys_mean = np.nanmean(ys, axis = 1)
    ys_sd = np.nanstd(ys, ddof=1, axis = 1)
    ys_upper = ys_mean + outlier_sd * ys_sd
    ys_lower = ys_mean - outlier_sd * ys_sd

    # calculate on outlier
    outlier = (ym > ys_upper) | (ym < ys_lower)
    y[outlier, k] = np.nan
    y_outlier = y[outlier,]

    # weighted average on outliers
    tt = list(range(1, k+1))
    w = tt +[k+1] + tt[::-1]
    masked_y = np.ma.masked_array(y_outlier, np.isnan(y_outlier))
    avg = np.ma.average(masked_y, axis = 1, weights = w)
    avg = avg.filled(np.nan)

    # assign back
    smooth_x = x
    smooth_x[outlier] = avg

    return smooth_x

#------------------------------------------------------------------------------#
def trismooth_1d_plus(x, k = 10, outlier_sd = 5,
                      trim_quantile = 0.1, edge_repeat = False):
    """triangular smoothing on a vector, only smooth outliers
    k : # of points to consider on the left and the right of a point
    outlier_sd : mean +- outlier_sd * sd is the outlier in the region
    trim_quantile : remove trim_quantile extreme data points
                    to calcuate mean and sd
                    for both upper and lower quantile
    edge_repeat : False - boundaries extended by filling nan
        True - boundaries extended by repeating edge points
    """

    x = np.array(x)
    y = np.zeros ((len (x), 2*k + 1), dtype = x.dtype)
    y[:, k] = x
    for i in range(k):
        j = k - i
        y[j:,i] = x[:-j]
        y[:-j,-(i+1)] = x[j:]
        if (edge_repeat):
            y[:j,i] = x[0]
            y[-j:,-(i+1)] = x[-1]
        else:
            y[:j,i] = np.nan
            y[-j:,-(i+1)] = np.nan

    # weights
    tt = list(range(1, k+1))
    w = np.array(tt +[k+1] + tt[::-1])

    # smooth
    smooth_x = []
    for i in range(len(x)):
        z = y[i,:]
        zm = z[k]
        z[k] = np.nan

        with np.errstate(invalid='ignore'):
            tf = (z < np.nanquantile(z, trim_quantile)) |\
                 (z > np.nanquantile(z, 1 - trim_quantile))
        z[tf] = np.nan

        z_mean = np.nanmean(z)
        z_sd = np.nanstd(z, ddof=1)
        z_upper = z_mean + outlier_sd * z_sd
        z_lower = z_mean - outlier_sd * z_sd

        if (zm > z_upper or zm < z_lower):
            idx = np.where(~np.isnan(z))[0]
            smooth_x.append(np.average(z[idx], weights = w[idx]))
        else:
            smooth_x.append(zm)

    return smooth_x

#------------------------------------------------------------------------------#
def get_mode(x):
    """get mode from density
    """
    x = np.array(x)
    x = x[np.isfinite(x)]
    density = stats.kde.gaussian_kde(x)
    v = np.linspace(0, np.quantile(x, 0.99), 500)
    density = density(v)
    return v[density.argmax()]

#------------------------------------------------------------------------------#
# output function
#------------------------------------------------------------------------------#
def get_cn_segment(bin_count_clean, binsize = 10000, nperm = 1000, p = 0.01,
        method = 'cutoff', min_cn = None, std_scale = 8):
    """ get copy number segmentations bin count data
    output a dataframe of amplicons

    Parameters
    ----------
    df : counts in genomic bins
    method : method used to perform segmentation
        'cutoff' : use hard cutoff of min_cn to choose amplified regions
        'cbs' : use cbs method to choose amplified regions

    std_scale : mean + std_scale * std as min_cn
    binsize : target genomic bin size to calculate read counts

    output
    ----------
    segmented regions in dataframe
    """
    df = bin_count_clean.copy()

    # set min_cn by mean + n * std
    if (min_cn is None):
        xx = np.array(df['CN'])
        min_cn = np.mean(xx) + std_scale * np.std(xx, ddof = 1)
        print(min_cn)

    # calculate log2 ratio if use cbs
    if (method == 'cbs'):
        df['L2R'] = np.log2(np.array(df.CN / 2))

    # call copy number segmentations in each chromosome
    cnsegs = []
    for chr in set(df['Chrom']):
        # print(chr)
        dfsub = df[df['Chrom'] == chr]
        # seg = []
        if (method == 'cutoff'):
            x = np.array(dfsub['CN'])
            seg = simple_amplicon_segment(x, min_cn)
            seg = [list(dfsub.iloc[i[0],0:2]) +
                   list(dfsub.iloc[i[1]-1,1:2]+binsize) +
                   [i[2]] for i in seg]
        else:
            x = np.array(dfsub['L2R'])
            seg = cbs.cbs_segment(x, nperm=nperm, p=p)
            seg = cbs.cbs_recheck(x, seg, nperm=nperm, p=p)
            # cbs.cbs_plot_segment(x, seg) # plot function
            seg = [list(dfsub.iloc[i[0],0:2]) +
                   list(dfsub.iloc[i[1]-1,1:2]+binsize) +
                   [2**i[2]*2] for i in seg]
        cnsegs.append(pd.DataFrame(seg,
            columns=['Chrom', 'Start', 'End', 'CN']))
        seg_df = pd.concat(cnsegs)

    return seg_df[seg_df.CN >= min_cn]

#------------------------------------------------------------------------------#
# simple segmentation
#------------------------------------------------------------------------------#
def simple_seg(bin_count, binsize = 10000, min_cn = 5, std_scale = 8):
    """ get copy number segmentations bin count data
    output a dataframe of amplicons
    use hard cutoff of min_cn to choose amplified regions

    Parameters
    ----------
    bin_count : counts in genomic bins
    std_scale : mean + std_scale * std as min_cn
    binsize : target genomic bin size to calculate read counts
    min_cn : minimal copy number to consider

    output
    ----------
    segmented regions in dataframe
    """
    df = bin_count.copy()

    # set min_cn by mean + n * std
    if (min_cn is None):
        xx = np.array(df['CN'])
        min_cn = np.mean(xx) + std_scale * np.std(xx, ddof = 1)
        print(min_cn)

    # call copy number segmentations in each chromosome
    cnsegs = []
    for chr in set(df['Chrom']):
        # print(chr)
        dfsub = df[df['Chrom'] == chr]

        # segmentation
        x = np.array(dfsub['CN'])
        seg = simple_amplicon_segment(x, min_cn)
        seg = [list(dfsub.iloc[i[0],0:2]) +
               list(dfsub.iloc[i[1]-1,1:2]+binsize) +
               [i[2]] for i in seg]

        cnsegs.append(pd.DataFrame(seg,
            columns=['Chrom', 'Start', 'End', 'CN']))
        seg_df = pd.concat(cnsegs)

    return seg_df[seg_df.CN >= min_cn]

#------------------------------------------------------------------------------#
# simple calling of amplicons
#------------------------------------------------------------------------------#
def simple_amplicon_segment(x, min_cn = 5, min_dist_to_merge = 2):
    """
    Simple way to identify amplified segments
    score of genomic bins higher than a certain score
    let's say bins with more than 5 copies

    min_dist_to_merge: merge closeby intervals
    """

    xx = np.where(x > min_cn)[0].tolist()
    segments =[]

    # get the continuous numbers in a list with a tolerance
    grp = groupby(zip(xx, xx[1:] + xx[:0]), \
        lambda x: x[1] - x[0] <= min_dist_to_merge)
    for k,g in grp:
        if (k == True):
            gg = [i for i in g]
            segments.append([gg[0][0], gg[-1][1]+1])

    # get the segmentation mean table
    seg = []
    for j in segments:
        seg_start = j[0]
        seg_end = j[1]
        seg_mean = np.mean(x[seg_start:seg_end])
        seg.append([seg_start, seg_end, seg_mean])

    return seg

#------------------------------------------------------------------------------#
def segment_filter_by_size(cn_amp, binsize = 10000, fold = 5):
    """
    filter copy number segments by size (remove small spikes)

    """
    # old version
    # return cn_amp[cn_amp['End'] - cn_amp['Start'] >= fold * binsize]

    cn_amp_merged = misc.merge_bed(cn_amp, gap = 100000)
    cn_amp_drop = pd.DataFrame([
        row for row in cn_amp_merged if (row[2] - row[1] < fold * binsize
        )], columns = cn_amp.columns[0:3])
    df = pd.merge(cn_amp, cn_amp_drop, indicator = True, how = 'left'
        ).query('_merge == "left_only"').drop('_merge', axis = 1)
    return df

#------------------------------------------------------------------------------#
def bin_filter_by_blacklist(bin_count, blacklistfile, \
    binsize = 10000, ext = 0, set_value = 2):
    """
    set the bin CN to 2 (diploid), if overlap with blacklist

    ext - extend both end on ext * binsize
    """
    bin_count_ = bin_count.copy()

    # convert bed file to gr
    gr = GRange(blacklistfile, 'bedfile')

    # init blacklist bins
    bl_bin = defaultdict(lambda: False)

    # save blacklist bins to dict
    for _gr in gr.gr:
        _start = math.floor(_gr[1].start/binsize) - ext
        _end = math.floor(_gr[1].stop/binsize) + 1 + ext
        for i in range(_start, _end):
            bl_bin[(_gr[0], i*binsize)] = True

    # get the bool vector of blacklist bins
    tf = [bl_bin[row[1:3]] for row in bin_count_.itertuples()]

    # set blacklist bins to 0
    # bin_count_.loc[tf, 'Count'] = set_value
    # bin_count_['CN'] = 2*bin_count_['Count']/np.mean(bin_count_['Count'])

    # set blacklist bins CN to 2, left count unchanged
    bin_count_.loc[tf, 'CN'] = set_value

    return bin_count_

#------------------------------------------------------------------------------#
# identify cn segment boundaries by using clip reads
#------------------------------------------------------------------------------#
def finemap_cn_segment_boundary(cn_amp, bp_all, binsize = 10000):
    """
    search break points for single base resolution
    boundaries of amplified segments
    cn_amp: copy number high amplicons
    bp_all: bp_all
    """
    cn_amp = cn_amp[['Chrom', 'Start', 'End', 'CN']]
    seg_bdry = []
    for idx, seg in cn_amp.iterrows():
        bdry = list(seg) + [0, 0]

        # left boundary
        left = bp_all['left']
        left = left[left['Chrom'] == seg[0]]
        left = left.assign(dist = abs(seg[1]-left['Coord']))
        cand = left[left['dist'] <= binsize]
        if (not cand.empty):
            cand = cand.sort_values(['dist'])
            cand = cand.sort_values('Count', ascending=False)
            bdry[1] = list(cand['Coord'])[0]
            bdry[4] = list(cand['Count'])[0]

        # right boundary
        right = bp_all['right']
        right = right[right['Chrom'] == seg[0]]
        right = right.assign(dist = abs(seg[2]-right['Coord']))
        cand = right[right['dist'] <= binsize]
        if (not cand.empty):
            cand = cand.sort_values(['dist'])
            cand = cand.sort_values('Count', ascending=False)
            bdry[2] = list(cand['Coord'])[0]
            bdry[5] = list(cand['Count'])[0]

        # output
        seg_bdry.append(bdry)

    colnames = ['Chrom', 'Start', 'End', 'CN', 'ClipLeft', 'ClipRight']
    return pd.DataFrame(seg_bdry, columns=colnames)

#------------------------------------------------------------------------------#
def get_bppair_from_cnamp(cn_amp, bamfile, min_nt = 4):
    """
    get the bppairs from the refined cnamp
    as addition to the bppairs identifed by bp_cand
    """
    op = []
    for row in cn_amp.itertuples():
        op.append([row[1], row[2], 'L'])
        op.append([row[1], row[3], 'R'])
    op

    tt = [pair[0] + pair[1] + list(
            join_breakpoint(bamfile, *pair, min_nt = min_nt)) \
        for pair in combinations(op, 2)]

    t2 = [row[0:6] + [row[6][1] + row[7][1]] \
          for row in tt \
              if row[6][1] > 0 and row[7][1] > 0 \
              and row[6][0] == row[7][0]]

    colnames = ["Chrom1", "Coord1", "Clip1", "Chrom2", "Coord2", "Clip2", 'Count']
    bp_pair_extra = pd.DataFrame(t2, columns = colnames)

    return bp_pair_extra

#------------------------------------------------------------------------------#
def segment_annotation(cn_amp, refgenef):
    """
    annotate the copy number segments by gene
    """
    return GRange.map_gene(cn_amp, refgenef)

#------------------------------------------------------------------------------#
# filter copy number segments by blacklist intervals (currently not in use)
#------------------------------------------------------------------------------#
def segment_filter_by_blacklist(cn_amp, blacklistfile, f_sort=True):
    """
    filter copy number segments by blacklist intervals

    """
    def filter_blacklist(df, blacklistfile):
        gr1 = GRange(df, 'dataframe_hasend')
        gr2 = GRange(blacklistfile, 'bedfile')
        # extend 100 bp on both end of the blacklist
        gr = gr1.intersect(gr2, a_extend = 0, b_extend = 100, invert = True)

        return pd.DataFrame(\
            [[row[0]] + [row[1].start, row[1].stop] + list(row[2]) for row in gr.gr]\
            , columns=df.columns)

    df = filter_blacklist(cn_amp, blacklistfile)

    if (f_sort):
        df = df.sort_values(['Chrom', 'Start', 'End'])

    return df

max_stars_repo_name: mthjwu/ask

max_issues_repo_name: mthjwu/ask

max_forks_repo_name: mthjwu/ask

reason: import numpy,import scipy,import statsmodels,from statsmodels