text:
import networkx as nx
import math
from logging import getLogger
from libcity.model.abstract_traffic_tradition_model import AbstractTraditionModel
from libcity.utils.GPS_utils import radian2angle, R_EARTH, angle2radian, dist, init_bearing
import numpy as np

class IVMM(AbstractTraditionModel):
    """
    IVMM
    """

    def __init__(self, config, data_feature):

        super().__init__(config, data_feature)

        # logger
        self._logger = getLogger()

        # model param
        self.k = config.get('k', 5)  # maximum number of candidates of any sampling points
        self.r = config.get('r', 100)  # radius of road segment candidate
        self.mu = config.get('mu', 0)
        self.sigma = config.get('sigma', 10)
        self.beta = config.get('beta', 40)
        self.window_size = config.get('window_size', 40)

        # data param
        self.with_time = data_feature.get('with_time', True)
        self.with_rd_speed = data_feature.get('with_rd_speed', True)
        self.delta_time = data_feature.get('delta_time', True)

        # data cache
        self.rd_nwk = None
        self.tra_id = None
        self.trajectory = None
        self.candidates = list()
        self.res_dct = dict()

        # To facilitate the search of candidate points, the road network is indexed using a grid.
        self.lon_r = None
        self.lat_r = None

    def run(self, data):
        # road network & trajectory
        self.rd_nwk = data['rd_nwk']
        trajectory = data['trajectory']

        # set lon_radius and lat_radius based on the first Node of rd_nwk
        self._set_lon_lat_radius(
            self.rd_nwk.nodes[list(self.rd_nwk.nodes)[0]]['lon'],
            self.rd_nwk.nodes[list(self.rd_nwk.nodes)[0]]['lat']
        )

        # deal with every trajectories
        for usr_id, usr_value in trajectory.items():
            self.usr_id = usr_id
            for traj_id, value in usr_value.items():
                self._logger.info('begin map matching, usr_id:%d traj_id:%d' % (usr_id, traj_id))
                self.traj_id = traj_id
                self.trajectory = value
                self._run_one_tra()
                self._logger.info('finish map matching, usr_id:%d traj_id:%d' % (usr_id, traj_id))
        return self.res_dct

    def _run_one_tra(self):
        """
        run ST-Matching for one trajectory
        self.trajectory self.rd_nwk
        Returns:
        """
        self._get_candidates()
        self._logger.info('finish getting candidates')
        self._observation_probability()
        self._logger.info('finish calculating observation probability')
        self._transmission_probability()
        self._logger.info('finish calculating transmission probability')
        if self.with_rd_speed and self.with_time:
            self._temporal_analysis()
            self._logger.info('finish spatial analysis')
        self._score_matrix()
        self._logger.info('finish building score matrix')
        self._interactive_voting()
        self._logger.info('finish interactive voting')
        self._find_matched_sequence()
        self._logger.info('finish finding matched sequence')

    def _set_lon_lat_radius(self, lon, lat):
        """
        get longitude range & latitude range (because radius is actually achieved by a grid search)
        Args:
            lon: longitude local
            lat: latitude local
            self.r
        Returns:
            self.lon_r
            self.lat_r
        """

        # lat_r
        self.lat_r = radian2angle(self.r / R_EARTH)

        # lon_r
        r_prime = R_EARTH * math.cos(angle2radian(lat))
        self.lon_r = radian2angle(self.r / r_prime)

    def _point_edge_dist(self, lon, lat, edge):
        lat_origin = angle2radian(self.rd_nwk.nodes[edge[0]]['lat'])
        lon_origin = angle2radian(self.rd_nwk.nodes[edge[0]]['lon'])
        lat_dest = angle2radian(self.rd_nwk.nodes[edge[1]]['lat'])
        lon_dest = angle2radian(self.rd_nwk.nodes[edge[1]]['lon'])

        a = dist(angle2radian(lat), angle2radian(lon), lat_origin, lon_origin)
        b = dist(angle2radian(lat), angle2radian(lon), lat_dest, lon_dest)
        c = dist(lat_origin, lon_origin, lat_dest, lon_dest)

        # if origin point is the closest
        if b ** 2 > a ** 2 + c ** 2:
            return a, edge[0]  # distance, point

        # if destination point is the closest
        elif a ** 2 > b ** 2 + c ** 2:
            return b, edge[1]

        if c == 0:
            return a, edge[0]
        # otherwise, calculate the Vertical length
        p = (a + b + c) / 2
        s = math.sqrt(p * math.fabs(p - a) * math.fabs(p - b) * math.fabs(p - c))
        return 2 * s / c, None

    def _get_candidates(self):
        """
        get candidates of each GPS sample with given road network
        Returns:
            self.candidates: a list of list.
                In each list are tuples (edge, distance, node)
        """

        # get trajectory without time
        traj_lon_lat = self.trajectory[:, 1:3]
        assert traj_lon_lat.shape[1] == 2

        # for every GPS sample
        for i in range(traj_lon_lat.shape[0]):

            candidate_i = set()
            lon, lat = traj_lon_lat[i, :]

            # for every edge
            for j in self.rd_nwk.edges:
                origin, dest = j[:2]
                lat_origin = self.rd_nwk.nodes[origin]['lat']
                lon_origin = self.rd_nwk.nodes[origin]['lon']
                lat_dest = self.rd_nwk.nodes[dest]['lat']
                lon_dest = self.rd_nwk.nodes[dest]['lon']
                if lat - self.lat_r <= lat_origin <= lat + self.lat_r \
                        and lon - self.lon_r <= lon_origin <= lon + self.lon_r \
                        or lat - self.lat_r <= lat_dest <= lat + self.lat_r \
                        and lon - self.lon_r <= lon_dest <= lon + self.lon_r:
                    candidate_i.add((origin, dest))
                elif lat - self.lat_r <= lat_origin / 2 + lat_dest / 2 <= lat + self.lat_r \
                        and lon - self.lon_r <= lon_origin / 2 + lon_dest / 2 <= lon + self.lon_r:
                    candidate_i.add((origin, dest))
            candidate_i_m = list()  # (edge, distance, point)
            for edge in candidate_i:
                distance, node = self._point_edge_dist(lon, lat, edge)
                candidate_i_m.append((edge, distance, node))
            candidate_i_m.sort(key=lambda a: a[1])  # asc
            candidate_i_m = candidate_i_m[:min(self.k, len(candidate_i_m))]
            candidate_i_k = dict()
            for edge, distance, node in candidate_i_m:
                candidate_i_k[edge] = {'distance': distance, 'node': node, "voted": 0, "pre_set": []}
            self.candidates.append(candidate_i_k)

    def _observation_probability(self):
        """

        Returns:

        """

        # for candidates of every node
        for candidate_i in self.candidates:
            for edge, dct in candidate_i.items():
                candidate_i[edge]['N'] = (1 / math.sqrt(2 * math.pi) / self.sigma * math.exp(
                    - (dct['distance'] - self.mu) ** 2 / (2 * self.sigma ** 2)))

    def _transmission_probability(self):
        """

        Returns:

        """
        i = 1
        while i < len(self.candidates):
            # j and k
            j = i - 1
            if len(self.candidates[i]) == 0:
                k = i + 1
                while len(self.candidates[k]) == 0 and k < len(self.candidates):
                    k += 1
                if k == len(self.candidates):
                    break
            else:
                k = i
            d = dist(
                angle2radian(self.trajectory[j][2]),
                angle2radian(self.trajectory[j][1]),
                angle2radian(self.trajectory[k][2]),
                angle2radian(self.trajectory[k][1])
            )
            for edge_j, dct_j in self.candidates[j].items():
                for edge_k, dct_k in self.candidates[k].items():
                    brng_jk = init_bearing(
                        angle2radian(self.trajectory[j][2]),
                        angle2radian(self.trajectory[j][1]),
                        angle2radian(self.trajectory[k][2]),
                        angle2radian(self.trajectory[k][1])
                    )
                    brng_edge_j = init_bearing(
                        angle2radian(self.rd_nwk.nodes[edge_j[0]]['lat']),
                        angle2radian(self.rd_nwk.nodes[edge_j[0]]['lon']),
                        angle2radian(self.rd_nwk.nodes[edge_j[1]]['lat']),
                        angle2radian(self.rd_nwk.nodes[edge_j[1]]['lon']),
                    )
                    try:
                        if dct_j['node'] is not None and dct_k['node'] is not None:
                            result = d / nx.astar_path_length(self.rd_nwk, dct_j['node'], dct_k['node'],
                                                              weight='distance')
                        elif dct_j['node'] is not None:
                            nd2_origin = edge_k[0]
                            lon, lat = self.rd_nwk.nodes[nd2_origin]['lon'], self.rd_nwk.nodes[nd2_origin]['lat']
                            path_len = nx.astar_path_length(self.rd_nwk, dct_j['node'], nd2_origin, weight='distance')
                            path_len += math.sqrt(
                                math.fabs(
                                    dist(
                                        angle2radian(self.trajectory[k][2]),
                                        angle2radian(self.trajectory[k][1]),
                                        angle2radian(lat),
                                        angle2radian(lon)
                                    ) ** 2 - dct_k['distance'] ** 2
                                )
                            )
                            if edge_j[1] == dct_j['edge']:
                                path_len += self.rd_nwk[edge_j[0]][edge_j[1]]['distance'] * 2
                            result = d / path_len

                        elif dct_k['node'] is not None:
                            nd1_destination = edge_j[1]
                            lon, lat = self.rd_nwk.nodes[nd1_destination]['lon'], self.rd_nwk.nodes[nd1_destination][
                                'lat']
                            path_len = nx.astar_path_length(self.rd_nwk, nd1_destination, dct_k['node'],
                                                            weight='distance')
                            path_len += math.sqrt(
                                math.fabs(
                                    dist(
                                        angle2radian(self.trajectory[j][2]),
                                        angle2radian(self.trajectory[j][1]),
                                        angle2radian(lat),
                                        angle2radian(lon)
                                    ) ** 2 - dct_j['distance'] ** 2
                                )
                            )
                            if edge_k[1] == dct_k['node']:
                                path_len += self.rd_nwk[edge_k[0]][edge_k[1]]['distance'] * 2
                            result = d / path_len
                        else:
                            if edge_j == edge_k and math.fabs(brng_edge_j - brng_jk) < 90:
                                result = 1
                            else:
                                nd1_destination = edge_j[1]
                                lon1, lat1 = self.rd_nwk.nodes[nd1_destination]['lon'], \
                                    self.rd_nwk.nodes[nd1_destination]['lat']
                                nd2_origin = edge_k[0]
                                lon2, lat2 = self.rd_nwk.nodes[nd2_origin]['lon'], self.rd_nwk.nodes[nd2_origin]['lat']
                                result = d / (
                                        nx.astar_path_length(self.rd_nwk, nd1_destination, nd2_origin,
                                                             weight='distance')
                                        + math.sqrt(
                                            math.fabs(
                                                dist(
                                                    angle2radian(self.trajectory[j][2]),
                                                    angle2radian(self.trajectory[j][1]),
                                                    angle2radian(lat1),
                                                    angle2radian(lon1)
                                                ) ** 2 - dct_j['distance'] ** 2
                                            )
                                        )
                                        + math.sqrt(
                                            math.fabs(
                                                dist(
                                                    angle2radian(self.trajectory[k][2]),
                                                    angle2radian(self.trajectory[k][1]),
                                                    angle2radian(lat2),
                                                    angle2radian(lon2)
                                                ) ** 2 - dct_k['distance'] ** 2
                                            )
                                        )
                                )
                        if 'V' in dct_j.keys():
                            dct_j['V'][edge_k] = min(result, 1)
                        else:
                            dct_j['V'] = {edge_k: min(result, 1)}
                    except:
                        if 'V' in dct_j.keys():
                            dct_j['V'][edge_k] = 0
                        else:
                            dct_j['V'] = {edge_k: 0}
            i += 1

    def _temporal_analysis(self):
        """

        Returns:

        """
        # TODO

    def _score_matrix(self):
        sub_static_matrices = list()
        # for every GPS sample
        for i in range(len(self.candidates)):
            # current candidates: self.candidates[i]
            # prev candidates: self.candidates[j]
            # no current candidates
            if i == 0:
                continue
            if len(self.candidates[i]) == 0:
                sub_static_matrix = np.zeros((1, 1))
                sub_static_matrices.append(sub_static_matrix)
                continue
            # if there are current candidates, find prev index j
            j = i - 1
            while j >= 0 and len(self.candidates[j]) == 0:
                j -= 1
            # j >= 0, calculate score of candidates of GPS sample i
            if j < 0:
                sub_static_matrix = np.zeros((1, 1))
                sub_static_matrices.append(sub_static_matrix)
                continue
            sub_static_matrix = np.zeros((len(self.candidates[j].items()), len(self.candidates[i].items())))
            n = 0
            for edge, dct in self.candidates[i].items():
                m = 0
                for edge_pre, dct_pre in self.candidates[j].items():
                    sub_static_matrix[m][n] = dct['N'] * dct_pre['V'][edge] * (
                        1 if 'T' not in dct_pre.keys() else dct_pre['T'][edge])
                    m += 1
                n += 1
            sub_static_matrices.append(sub_static_matrix)
        traj_lon_lat = self.trajectory[:, 1:3]
        self.traj_points = []
        for i in range(len(traj_lon_lat)):
            W = []
            j = 0
            while j < len(traj_lon_lat):
                if i == j:
                    j += 1
                    continue
                lon_i, lat_i = traj_lon_lat[i, :]
                lon_j, lat_j = traj_lon_lat[j, :]
                euclid_distance = dist(angle2radian(lat_i), angle2radian(lon_i),
                                       angle2radian(lat_j), angle2radian(lon_j))
                W.append(pow(math.e, -(euclid_distance ** 2 / self.beta ** 2)))
                j += 1
            fai = []
            self.sub_static_matrices = sub_static_matrices
            self.W = W
            for k in range(len(W)):
                fai.append(sub_static_matrices[k] * W[k])
            traj_point = {"time": self.trajectory[i][0], "lon": self.trajectory[i][1],
                          "lat": self.trajectory[i][2], "fai": fai, "candidates": self.candidates[i]}
            self.traj_points.append(traj_point)  # Φ(fai) is the weighted score matrix

    def _interactive_voting(self):
        for traj_point in self.traj_points:
            fai = traj_point["fai"]
            for candidate in traj_point["candidates"].values():
                self._find_sequence(fai, candidate)

    def _find_sequence(self, fai, candidate):
        for i in range(len(self.candidates)):
            if candidate not in self.candidates[i].values():
                i += 1
            else:
                break
        current_p = i
        j = 0
        for c in self.candidates[i].values():
            if candidate == c:
                current_c = j
                break
            else:
                j += 1
        # get the location of current traj_point
        i += 1
        if i == len(self.candidates):
            # traj_point is in the last column
            last_point = None
            last_score = 0
        while i < len(self.candidates):
            if len(self.candidates[i]) == 0:
                i += 1
                if i == len(self.candidates):
                    last_point = None
                    last_score = 0
                    continue
            # ensure i has candidates
            j = i - 1
            while j > current_p and len(self.candidates[j]) == 0:
                j -= 1
            # ensure j has candidates
            k = 0
            for edge, dct in self.candidates[i].items():
                max_score = -float("inf")
                if j == current_p:
                    tmp = fai[i - 1][current_c][k]
                    if tmp > max_score:
                        dct["score"] = tmp
                else:
                    m = 0
                    for edge_pre, dct_pre in self.candidates[j].items():
                        tmp = dct_pre['score'] + fai[i - 1][m][k]
                        if tmp > max_score:
                            max_score = tmp
                            pre = dct_pre
                        m += 1
                    dct["score"] = max_score
                    for c in pre["pre_set"]:
                        dct["pre_set"].append(c)
                    dct["pre_set"].append(pre)
                k += 1
            i += 1
            if i == len(self.candidates):
                last_point = dct
                last_score = dct["score"]
        i = current_p - 1
        if i == -1:
            first_point = None
            first_score = 0
        while i >= 0:
            if len(self.candidates[i]) == 0:
                i -= 1
                if i == -1:
                    first_point = dct
                    first_score = dct["score"]
                continue
            j = i + 1
            while j < current_p and len(self.candidates[j]) == 0:
                j += 1
            k = 0
            for edge, dct in self.candidates[i].items():
                max_score = -float("inf")
                if j == current_p:
                    tmp = fai[j - 1][k][current_c]
                    if tmp > max_score:
                        dct["score"] = tmp
                else:
                    m = 0
                    for edge_pre, dct_pre in self.candidates[j].items():
                        tmp = dct_pre['score'] + fai[j - 1][k][m]
                        if tmp > max_score:
                            max_score = tmp
                            pre = dct_pre
                        m += 1
                    dct["score"] = max_score
                    for c in pre["pre_set"]:
                        dct["pre_set"].append(c)
                    dct["pre_set"].append(pre)
                k += 1
            i -= 1
            if i == -1:
                first_point = dct
                first_score = dct["score"]
        candidate["fValue"] = first_score + last_score
        if first_point is not None:
            for point in first_point["pre_set"]:
                point["voted"] += 1
        if last_point is not None:
            for point in last_point["pre_set"]:
                point["voted"] += 1
        candidate["voted"] += 1
        for traj_point in self.traj_points:
            for c in traj_point["candidates"].values():
                c["pre_set"] = []

    def _find_matched_sequence(self):
        res = []
        for traj_point in self.traj_points:
            if len(traj_point["candidates"]) != 0:
                max_voted = 0
                fValue = 0
                for key in traj_point["candidates"].keys():
                    candidate = traj_point["candidates"][key]
                    tmp = candidate["voted"]
                    if tmp > max_voted:
                        max_voted = tmp
                        fValue = candidate["fValue"]
                        max_voted_key = key
                    elif tmp == max_voted:
                        if candidate["fValue"] > fValue:
                            max_voted = tmp
                            fValue = candidate["fValue"]
                            max_voted_key = key
            else:
                max_voted_key = None
            res.append(max_voted_key)

        # to geo_id
        res_lst_rel = np.array(list(map(lambda x: self.rd_nwk.edges[x]['geo_id'] if x is not None else None, res)))
        dyna_id_lst = self.trajectory[:, 0].astype(int)
        if self.with_time:
            time_lst = self.trajectory[:, 3]
            res_all = np.stack([dyna_id_lst, res_lst_rel, time_lst], axis=1)
        else:
            res_all = np.stack([dyna_id_lst, res_lst_rel], axis=1)

        # set self.res_dct
        if self.usr_id in self.res_dct.keys():
            self.res_dct[self.usr_id][self.traj_id] = res_all

        else:
            self.res_dct[self.usr_id] = {self.traj_id: res_all}

reason: import numpy,import networkx


text:
import argparse
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt

# Local Moduels
from asiCD.asiCD_utils import img_from_file
from asiCD.asiCD_utils import load_json

def gen_flat_img(img_inst, ann_inst):

    img_flat = np.zeros([img_inst.shape[0]*img_inst.shape[1], 4])

    for i in range(4):
        if i == 3:
            img_flat[:, i] = ann_inst[:, :, 0].flatten()
        else:
            img_flat[:, i] = img_inst[:, :, i].flatten()

    return img_flat

def create_pixel_dist(img_files, ann_files):
    print(f"Image File: {img_files[0]}")
    print(f"Mask File: {ann_files[0]}")
    img_inst = img_from_file(img_files[96])/255
    ann_inst = img_from_file(ann_files[96], load_grey=True)/255

    img_flat = gen_flat_img(img_inst, ann_inst)

    # Pixels within masked area
    img_flat = img_flat[img_flat[:, 3] == 1]

    print(img_flat.shape)

    plt.figure(figsize=(12, 6))

    plt.subplot(1, 2, 1)
    plt.title("Green vs. Red")
    plt.scatter(x=img_flat[:, 0], y=img_flat[:, 1],
                c=img_flat[:, 0:3], alpha=0.1)

    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.xlabel("R")
    plt.ylabel("G")

    plt.subplot(1, 2, 2)
    plt.title("Blue vs. Red")
    plt.scatter(x=img_flat[:, 0], y=img_flat[:, 2],
                c=img_flat[:, 0:3], alpha=0.1)

    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.xlabel("R")
    plt.ylabel("B")

    plt.show()

    return None

def main():

    # construct the argument parse and parse the arguments
    ap = argparse.ArgumentParser(
        description='Detect and encircle the sun for a provided image')
    ap.add_argument("-i", "--image_id", type=int, default=0,
                    help="index of image file")
    ap.add_argument("-r", "--radius", type=int, default=120,
                    help="radius of sun circle for sun_remover_v1")
    ap.add_argument("-t", "--low_thres", type=int, default=253,
                    help="lower thershold for sun_remover_v2")
    args = vars(ap.parse_args())

    # Fetch image files
    img_files = list(Path("dataset/swimseg-2/train").glob("**/*.png"))
    ann_files = list(Path("dataset/swimseg-2/train_labels").glob("**/*.png"))
    config = load_json("config.json")

    print(f"Image Files: {len(img_files)}")
    print(f"Mask Files: {len(ann_files)}")

    create_pixel_dist(img_files, ann_files)

    return None

if __name__ == "__main__":
    main()

ext: py

max_stars_repo_name: vinkmr/asiCD

max_issues_repo_name: vinkmr/asiCD

max_forks_repo_name: vinkmr/asiCD

reason: import numpy


text:
import numpy as np
from scipy.stats import norm

class GMM:
    def __init__(self):
        self.k = None
        self.w = None
        self.mu = None
        self.sigma = None
        self.phi = None
        self.nu = None
        self.gmm_pdfs = None
        self.n_params = None
        self.type = 'GMM'

        self._finite_obs = None
        self._state_props = None

    def set_params(self, config):

        params = {'k', 'w', 'mu', 'sigma', 'phi', 'nu', 'n_params'}
        self.__dict__.update((param, np.array(value)) for param, value in config.items() if param in params)

        self.setup_gmm_pdfs()

    def initialize(self, k, stats):

        self.k = k

        self.w = np.tile(1 / k, (2, k))  # Uniform weights
        # Symmetrically place initial distributions around the median, one higher than the other
        self.mu = np.vstack((np.linspace(stats['P60'], stats['P75'], k), np.linspace(stats['P25'], stats['P40'], k)))
        self.sigma = np.tile(stats['continuous_variance'], (2, k))  # Use overall variance
        self.phi = np.array((0.05, 0.05))  # P(NaN)
        self.nu = np.array((0.1, 0.1))  # P(zero)

        self._finite_obs = stats['finite_obs']
        self._state_props = np.array((0.5, 0.5))

        self.setup_gmm_pdfs()
        self.n_params = 6 * k + 2

    def setup_gmm_pdfs(self):
        self.gmm_pdfs = [[], []]
        for i in range(self.k):
            self.gmm_pdfs[0].append(norm(loc=self.mu[0, i], scale=np.sqrt(self.sigma[0, i])))
            self.gmm_pdfs[1].append(norm(loc=self.mu[1, i], scale=np.sqrt(self.sigma[1, i])))

    def compute_emissions(self, transcript, reference=False):

        # Peform computations for weighted kernels
        self.wnormpdfs(transcript)

        if reference:
            return

        # Correct kernel emissions for phi and nu
        # transcript.gamma_gmm_k *= (1 - self.nu[:, np.newaxis, np.newaxis] - self.phi[:, np.newaxis, np.newaxis])

        b = np.sum(transcript.gamma_gmm_k, axis=1)  # Sum emissions of all kernels
        b *= (1 - self.nu[:, np.newaxis] - self.phi[:, np.newaxis])

        # Assign emission likelihoods for discrete cases
        b[:, transcript.mask_0] = self.nu[:, np.newaxis]
        b[:, transcript.mask_nan] = self.phi[:, np.newaxis]

        transcript.B = b

    @staticmethod
    def post_process(transcript):
        transcript.gamma_gmm_k *= transcript.gamma[:, np.newaxis, :] / \
                                  np.sum(transcript.gamma_gmm_k, axis=1)[:, np.newaxis, :]

    def m_step(self, transcript):

        params = dict()

        gamma_gmm_k_sum = np.nansum(transcript.gamma_gmm_k, axis=2)

        # Overall unpaired proportion
        params["unpaired_prop"] = np.sum(transcript.gamma, axis=1)

        # Re-estimate phi and nu
        params["phi"] = np.sum(transcript.gamma[:, transcript.mask_nan], axis=1)
        params["nu"] = np.sum(transcript.gamma[:, transcript.mask_0], axis=1)
        params["phi_nu_norm"] = np.sum(transcript.gamma, axis=1)

        # Re-estimate GMM means and variances
        params["mu"] = np.nansum(transcript.gamma_gmm_k * transcript.obs[np.newaxis, np.newaxis, :], axis=2)
        sq_residual = (transcript.obs[np.newaxis, np.newaxis, :] - self.mu[:, :, np.newaxis]) ** 2
        params["sigma"] = np.nansum(transcript.gamma_gmm_k * sq_residual, axis=2)
        params["mu_sigma_norm"] = gamma_gmm_k_sum

        # Re-estimate the weight of the Gaussian mixture components - w
        params["w"] = gamma_gmm_k_sum  # Excluding observations not emitted from the GMM
        params["w_norm"] = np.sum(gamma_gmm_k_sum, axis=1)[:, np.newaxis]

        return params

    def update_from_pseudocounts(self, pseudocounts, nan=False):
        """
        Update emission model parameters given the relevant pseudocount sums over all transcripts.
        """

        self.mu = pseudocounts['mu'] / pseudocounts['mu_sigma_norm']
        self.sigma = pseudocounts['sigma'] / pseudocounts['mu_sigma_norm']
        self.w = pseudocounts['w'] / pseudocounts['w_norm']
        self.phi = pseudocounts['phi'] / pseudocounts['phi_nu_norm'] if nan else np.array((0.05, 0.05))
        self.nu = pseudocounts['nu'] / pseudocounts['phi_nu_norm']
        self._state_props = pseudocounts['unpaired_prop'] / pseudocounts['unpaired_prop'].sum()
        self.setup_gmm_pdfs()  # Construct new distribution objects

    def wnormpdfs(self, transcript):
        transcript.gamma_gmm_k = np.tile(np.nan, (2, self.k, len(transcript.obs)))
        for k in range(self.k):
            transcript.gamma_gmm_k[0, k, :] = self.w[0, k, np.newaxis] * self.gmm_pdfs[0][k].pdf(transcript.obs)
            transcript.gamma_gmm_k[1, k, :] = self.w[1, k, np.newaxis] * self.gmm_pdfs[1][k].pdf(transcript.obs)

        transcript.gamma_gmm_k[:, :, transcript.mask_finite] += 1e-20  # Padding to avoid numerical blow-ups
        transcript.gamma_gmm_k[:, :, transcript.mask_0] = np.nan  # Override zeros (would have been mapped to -Infinity)

    def snapshot(self):
        text = ""
        text += "{}:\n{}\n".format('w', np.array2string(self.w))
        text += "{}:\n{}\n".format('mu', np.array2string(self.mu))
        text += "{}:\n{}\n".format('sigma', np.array2string(self.sigma))
        text += "{}:\n{}\n".format('phi', np.array2string(self.phi))
        text += "{}:\n{}\n".format('nu', np.array2string(self.nu))
        return text

    def serialize(self):
        """
        Return a dictionary containing all of the parameters needed to describe the emission model.
        """
        return {'type': self.type,
                'k': self.k,
                'w': self.w.tolist(),
                'mu': self.mu.tolist(),
                'sigma': self.sigma.tolist(),
                'phi': self.phi.tolist(),
                'nu': self.nu.tolist(),
                'n_params': self.n_params}

    def reset(self):
        self.k = None
        self.w = None
        self.mu = None
        self.sigma = None
        self.phi = None
        self.nu = None
        self.gmm_pdfs = None
        self.n_params = None

max_stars_repo_name: AviranLab/patteRNA

max_issues_repo_name: AviranLab/patteRNA

max_forks_repo_name: AviranLab/patteRNA

reason: import numpy,from scipy


text:
import copy
import json
import logging
import math
import os
import pickle
import random

import numpy as np
import nni
import torch
import torch.nn as nn
import torch.optim as optim
from scipy import stats
from nni.nas.pytorch.utils import AverageMeterGroup
from torch.utils.tensorboard import SummaryWriter

from configs import Nb101Parser
from datasets.cifar10 import dataloader_cifar
from space.nb101 import Nb101Mutator, Nb101Network

from trainers.nb101 import train, validate
from trainers.utils import (
    AuxiliaryCrossEntropyLoss, CyclicIterator, Timer, accuracy, load_checkpoint,
    save_checkpoint, set_running_statistics, write_tensorboard
)

logger = logging.getLogger(__name__)

class Nb101RestartParser(Nb101Parser):
    def default_params(self):
        return {
            **super().default_params(),
            "restart_every": -1,
        }

    def validate_args(self, args):
        if args.restart_every < 0:
            args.restart_every = args.epochs
        assert args.epochs % args.restart_every == 0
        return super().validate_args(args)

def main():
    args = Nb101RestartParser.parse_configs()

    if args.cifar_split == "40k":
        train_split, valid_split = "train", "val"
    else:
        train_split, valid_split = "augment", "test"

    train_loader = dataloader_cifar("data/cifar10", train_split, args)
    valid_loader = dataloader_cifar("data/cifar10", valid_split, args)
    sanitize_loader = CyclicIterator(dataloader_cifar("data/cifar10", train_split, args, batch_size=args.bn_sanitize_batch_size))
    writer = SummaryWriter(args.tb_dir)

    model = Nb101Network(args)
    mutator = Nb101Mutator(model, validation_size=args.num_valid_arch, seed=args.seed, pruned=args.pruned)
    if args.resume_checkpoint:
        load_checkpoint(model, args.resume_checkpoint)
    model.cuda()

    if args.aux_weight > 0:
        criterion = AuxiliaryCrossEntropyLoss(args.aux_weight)
    else:
        criterion = nn.CrossEntropyLoss()
    optimizer = optim.SGD(model.parameters(), lr=args.initial_lr, momentum=args.momentum, weight_decay=args.weight_decay)
    scheduler = optim.lr_scheduler.CosineAnnealingWarmRestarts(optimizer, args.restart_every * len(train_loader), eta_min=args.ending_lr)

    validate(model, mutator, sanitize_loader, valid_loader, criterion, writer, args, 0)
    for epoch in range(1, args.epochs + 1):
        train(model, mutator, train_loader, criterion, optimizer, scheduler, writer, args, epoch)
        if (args.eval_every and epoch % args.eval_every == 0) or epoch == args.epochs:
            validate(model, mutator, sanitize_loader, valid_loader, criterion, writer, args, epoch)
    save_checkpoint(args, model, os.path.join(args.output_dir, "checkpoints", "final.pth.tar"))

if __name__ == "__main__":
    main()


text:
import unittest
import numpy as np
from megatron.layers.shaping import Cast, AddDim, OneHotRange, OneHotLabels, Reshape, Flatten
from megatron.layers.shaping import SplitDict, TimeSeries, Concatenate, Slice, Filter

class test_Cast(unittest.TestCase):
    def test_init(self):
        self.assertRaises(TypeError, Cast, 0)
        assert Cast(int).kwargs['new_type'] == int

    def test_transform(self):
        # incompatible cast (whole input)
        self.assertRaises(ValueError, Cast(int).transform, np.array(['a', 'b', 'c']))
        # incompatible cast (one element of input)
        self.assertRaises(ValueError, Cast(int).transform, np.array([1., 2., 'c']))
        # valid cast
        assert all(isinstance(v, np.integer) for v in Cast(int).transform(np.array(['1', '2', '3'])))

class test_AddDim(unittest.TestCase):
    def test_transform(self):
        # axis out of range
        self.assertRaises(ValueError, AddDim(axis=2).transform, np.ones(5))
        # axis at front
        assert AddDim(axis=0).transform(np.ones(5)).shape == (1, 5)
        # axis at back
        assert AddDim().transform(np.ones(5)).shape == (5, 1)
        # axis in middle
        assert AddDim(axis=1).transform(np.ones((5, 5))).shape == (5, 1, 5)

class test_OneHotRange(unittest.TestCase):
    def setUp(self):
        self.transformer = OneHotRange()
        self.transformer_strict = OneHotRange(strict=True)

    def test_partial_fit(self):
        self.transformer.partial_fit(np.array([-2, 4]))
        assert self.transformer.metadata['min_val'] == -2 and self.transformer.metadata['max_val'] == 4
        self.transformer.partial_fit(np.array([-3, 5, 4]))
        assert self.transformer.metadata['min_val'] == -3 and self.transformer.metadata['max_val'] == 5
        self.transformer.partial_fit(np.array([[1, 2], [1, 2]]))
        assert self.transformer.metadata['min_val'] == -3 and self.transformer.metadata['max_val'] == 5

    def test_fit(self):
        self.transformer.fit(np.array([-2, 4]))
        assert self.transformer.metadata['min_val'] == -2 and self.transformer.metadata['max_val'] == 4
        self.transformer.fit(np.array([[1, 2], [1, 2]]))
        assert self.transformer.metadata['min_val'] == 1 and self.transformer.metadata['max_val'] == 2

    def test_transform(self):
        self.assertRaises(TypeError, OneHotRange().fit, np.array(['a', 'b', 'c']))
        self.assertRaises(TypeError, OneHotRange().fit, np.array([1.5, 2, 3]))

        self.transformer_strict = OneHotRange(strict=True)
        self.transformer = OneHotRange()

        # value out of range (strict)
        self.transformer_strict.fit(np.array([-2, 1, 2, 3]))
        self.assertRaises(ValueError, self.transformer_strict.transform, np.array([-2, 0, 2, 4]))

        # value out of range (non-strict)
        self.transformer.fit(np.array([-2, 1, 2, 3]))
        assert np.array_equal(self.transformer.transform(np.array([-2, 0, 2, 4]))[-1], np.zeros(6))

        # 1D full range
        self.transformer.fit(np.array([1., 2, 3]))
        output = self.transformer.transform(np.array([1., 2, 3]))
        correct_output = np.zeros((3, 3))
        for i, j in enumerate([0, 1, 2]):
            correct_output[i, j] = 1
        assert np.array_equal(output, correct_output)

        # 1D sparse
        self.transformer.fit(np.array([-1, 2, 5, 7]))
        output = self.transformer.transform(np.array([5, 2, -1, 7]))
        correct_output = np.zeros((4, 9))
        for i, j in enumerate([6, 3, 0, 8]):
            correct_output[i, j] = 1
        assert np.array_equal(output, correct_output)

        # 2D
        self.transformer.fit(np.random.randint(-2, 3, (50, 50)))
        assert self.transformer.transform(np.random.randint(-2, 3, (40, 40))).shape == (40, 40, 5)

class test_OneHotLabels(unittest.TestCase):
    def setUp(self):
        self.transformer = OneHotLabels()
        self.transformer_strict = OneHotLabels(strict=True)

    def test_partial_fit(self):
        self.transformer.partial_fit(np.array(['a', 2]))
        assert set(self.transformer.metadata['categories']) == {'a', '2'}
        self.transformer.partial_fit(np.array(['a', 'b', 3]))
        assert set(self.transformer.metadata['categories']) == {'a', '2', 'b', '3'}
        self.transformer.partial_fit(np.array([['a', 'b'], [4, 2]]))
        assert set(self.transformer.metadata['categories']) == {'a', '2', 'b', '3', '4'}

    def test_fit(self):
        self.transformer.fit(np.array(['a', 2]))
        assert set(self.transformer.metadata['categories']) == {'a', '2'}
        self.transformer.fit(np.array([['a', 'b'], [3, 'a']]))
        assert set(self.transformer.metadata['categories']) == {'a', 'b', '3'}

    def test_transform(self):
        self.transformer_strict = OneHotLabels(strict=True)
        self.transformer = OneHotLabels()

        # value out of set (strict)
        self.transformer_strict.fit(['a', 1, 0.5, 'b'])
        self.assertRaises(ValueError, self.transformer_strict.transform, np.array(['a', '1', '0.5', 'b', 'c']))

        # value out of set (non-strict)
        self.transformer.fit(np.array(['a', 1, 0.5, 'c']))
        assert np.array_equal(self.transformer.transform(np.array(['a', '1', '0.5', 'b']))[-1], np.zeros(4))

        self.transformer.fit(np.array(['a', 'b', 'c']))
        output = self.transformer.transform(np.array(['a', 'b', 'c']))
        correct_output = np.zeros((3, 3))
        for i, j in enumerate([0, 1, 2]):
            correct_output[i, j] = 1
        assert np.array_equal(output, correct_output)

        self.transformer.fit(np.array([['a', 'b'], ['c', 'd']]))
        output = self.transformer.transform(np.array([['a', 'b'], ['c', 'd']]))
        assert output.shape == (2, 2, 4)

class test_Reshape(unittest.TestCase):
    def test_transform(self):
        self.assertRaises(ValueError, Reshape((3,2)).transform, np.ones(4))
        assert Reshape((2, 2)).transform(np.ones(4)).shape == (2, 2)

class test_Flatten(unittest.TestCase):
    def test_transform(self):
        assert Flatten().transform(np.ones((2, 2))).shape == (4,)

class test_SplitDict(unittest.TestCase):
    def test_transform(self):
        transformer = SplitDict(['a', 'b'])
        output = transformer.transform(np.array([{'a': 1, 'b': 2, 'c': 3}, {'a': 2, 'b': 3, 'c': 4}]))
        assert all(np.array_equal(a, b) for a, b in zip(output, [np.array([1, 2]), np.array([2, 3])]))

class test_TimeSeries(unittest.TestCase):
    def test_init(self):
        self.assertRaises(ValueError, TimeSeries, 1)
        assert TimeSeries(window_size=3).kwargs['window_size'] == 3

    def fit(self):
        transformer = TimeSeries(window_size=3)
        X = np.arange(100).reshape((20, 5))
        transformer.fit(X)
        assert transformer.metadata['shape'] == (5,)
        assert transformer.metadata['previous'] == np.zeros((2, 5))

    def test_transform(self):
        transformer = TimeSeries(window_size=3, reverse=True)
        X = np.arange(15).reshape((3, 5))
        transformer.fit(X)
        output = transformer.transform(X)
        correct_output = np.array([
            [[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 1, 2, 3, 4]],
            [[0, 0, 0, 0, 0], [0, 1, 2, 3, 4], [5, 6, 7, 8, 9]],
            [[0, 1, 2, 3, 4], [5, 6, 7, 8, 9], [10, 11, 12, 13, 14]]
        ])
        assert np.array_equal(output, correct_output)

class test_Concatenate(unittest.TestCase):
    def test_transform(self):
        X = [np.ones(100), np.ones(100), np.ones(100)]
        assert Concatenate().transform(*X).shape == (100, 3)

        X = [np.ones((100, 5)), np.ones((100, 3)), np.ones((100, 2))]
        assert Concatenate(axis=1).transform(*X).shape == (100, 10)
        self.assertRaises(ValueError, Concatenate(axis=0).transform, X[0], X[1], X[2])

class test_Slice(unittest.TestCase):
    def test_transform(self):
        self.assertRaises(ValueError, Slice(2, 2).transform, np.ones(100))

        output = Slice((2, 5), 3).transform(np.arange(100).reshape((10, 10)))
        correct_output = np.array([23, 33, 43])
        assert np.array_equal(output, correct_output)

        output = Slice(3, (0, 2), (0, 3, 2)).transform(np.arange(64).reshape((4, 4, 4)))
        correct_output = np.array([[48, 50], [52, 54]])
        assert np.array_equal(output, correct_output)

class test_Filter(unittest.TestCase):
    def test_transform(self):
        self.assertRaises(TypeError, Filter().transform, np.ones((3, 5)), np.array(['a', 'b', 'c']))

        output = Filter().transform(np.arange(5), np.array([1, 0, 1, 0, 1]))
        correct_output = np.array([0, 2, 4])
        assert np.array_equal(output, correct_output)

        output = Filter().transform(np.arange(16).reshape((4, 4)), np.array([1, 0, 0, 1]))
        correct_output = np.array([[0, 1, 2, 3], [12, 13, 14, 15]])
        assert np.array_equal(output, correct_output)

max_stars_repo_name: ntaylorwss/megatron

max_issues_repo_name: ntaylorwss/megatron

max_forks_repo_name: ntaylorwss/megatron


text:
#!/usr/bin/python3
import time
from typing import Dict, Any
import numpy as np
from wsl.networks.medinet.utils import regression_accuracy
from sklearn.metrics import r2_score
from scipy.stats import spearmanr
from monai.metrics import compute_roc_auc, compute_confusion_metric
import torch

def engine(loader: Any, checkpoint: Dict[str, Any], batchsize: int,
           classes: int, variable_type: str, error_range: int, is_train: bool):

    overall_loss = []
    all_preds = torch.zeros((0, classes))
    all_labels = torch.zeros((0, classes))
    labels_onehot = torch.FloatTensor(batchsize, classes).cuda()
    start = time.time()
    sigmoid = torch.nn.Sigmoid()

    with torch.set_grad_enabled(is_train):
        for iter_num, data in enumerate(loader):
            # name = data[0]
            imgs = data[1].cuda().float()
            labels = data[2].cuda()

            predicted = checkpoint['model'](imgs)
            loss = checkpoint['criterion'](predicted, labels)
            predicted, labels = predicted.detach(), labels.detach()

            if is_train:
                loss.backward()
                checkpoint['optimizer'].step()
                checkpoint['optimizer'].zero_grad()

            overall_loss.append(float(loss.item()))
            all_preds = torch.cat((predicted, all_preds))

            if variable_type == 'categorical':
                if labels_onehot.shape[0] != labels.shape[0]:
                    labels_onehot = torch.FloatTensor(labels.shape[0], classes).cuda()
                labels_onehot.zero_()
                labels_onehot.scatter_(1, labels.unsqueeze(dim=1), 1)
                all_labels = torch.cat((labels_onehot.float(), all_labels))
                predicted = predicted.max(dim=1)[1]  # for correct printing
            else:
                all_labels = torch.cat((labels, all_labels))

            speed = batchsize * iter_num // (time.time() - start)
            print('Epoch:', checkpoint['epoch'], 'Iter:', iter_num,
                  'Pred:', round(predicted.float().mean().item(), 3),
                  'Label:', round(labels.float().mean().item(), 3),
                  'Loss:', round(np.mean(overall_loss), 3),
                  'Speed:', int(speed), 'img/s', end='\r', flush=True)

    loss = np.mean(overall_loss)
    if variable_type == 'continous':
        all_labels, all_preds = all_labels.cpu(), all_preds.cpu()
        rmetric = r2_score(all_labels, all_preds)
        acc = regression_accuracy(all_labels, all_preds, error_range)
        spear, pvalue = spearmanr(all_preds, all_labels)
        summary = (f'Epoch Summary - Loss:{round(loss, 3)} Spearman:{round(spear, 2)} PValue:{round(pvalue, 3)} ' +
                   f'R2:{round(rmetric, 1)} Accuracy(at {error_range}):{round(100 * acc, 1)}')

    else:
        rmetric = compute_roc_auc(all_preds, all_labels, other_act=sigmoid)
        sens = compute_confusion_metric(all_preds, all_labels,
                                        activation=sigmoid, metric_name='sensitivity')
        spec = compute_confusion_metric(all_preds, all_labels,
                                        activation=sigmoid, metric_name='specificity')
        summary = (f'Epoch Summary- Loss:{round(loss, 3)}  ROC:{round(rmetric * 100, 1)} ' +
                   f'Sensitivity:{round(100 * sens, 1)}  Specificity: {round(100 * spec, 1)}')

    print(summary)
    return loss, rmetric, summary

max_stars_repo_name: nishanthta/wsl

max_issues_repo_name: nishanthta/wsl

max_forks_repo_name: nishanthta/wsl


text:
# -------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
# -------------------------------------------------------------
import numpy as np
from systemds.context import SystemDSContext

# create a random array
m1 = np.array(np.random.randint(100, size=5 * 5) + 1.01, dtype=np.double)
m1.shape = (5, 5)
# create another random array
m2 = np.array(np.random.randint(5, size=5 * 5) + 1, dtype=np.double)
m2.shape = (5, 5)

# Create a context
with SystemDSContext() as sds:
    # element-wise matrix multiplication, note that nothing is executed yet!
    m_res = sds.from_numpy(m1) * sds.from_numpy(m2)
    # lets do the actual computation in SystemDS! The result is an numpy array
    m_res_np = m_res.compute()
    print(m_res_np)

max_stars_repo_name: mdbloice/systemds

max_issues_repo_name: mdbloice/systemds

max_forks_repo_name: mdbloice/systemds


text:
from assertion import equals
from functions import LinearFunction
from initializers import ConstShapeInitializer
from mlp import Layer

import numpy as np

__all__ = ['should_be_success_calculate_output']

def should_be_success_calculate_output():
    layer = Layer(
        input_dimension=2,
        output_dimension=3,
        activation_function=LinearFunction(),
        weights_initializer=ConstShapeInitializer(
            np.asarray([
                [1., 2., 3.],
                [1., 2., 3.]
            ])
        ),
        biases_initializer=ConstShapeInitializer(
            np.asarray(
                [1., 2., 3.]
            )
        )
    )
    expected = np.asarray(
        [4., 8, 12.]
    )
    equals(expected, layer([1, 2]))


text:
import unittest

import numpy as np

from print_views import BlockBuilder, build_block, Projector, TOP_VIEW, \
    FRONT_VIEW, RIGHT_VIEW, \
    BOTTOM_VIEW, BACK_VIEW, LEFT_VIEW

__author__ = "Igor Mironov"
__copyright__ = "Copyright 2019, Igor Mironov"
__license__ = "Apache v2.0"

T = 1

F = 0

def assert_array_equal(a, b):
    np.testing.assert_array_equal(a, b)

class BlockBuilderTest(unittest.TestCase):

    def __init__(self, *args, **kwargs):
        super(BlockBuilderTest, self).__init__(*args, **kwargs)

        self.input_3 = '000001010011100101110111000b'
        self.array_3 = [
            [[F, F, F], [F, F, T], [F, T, F]],
            [[F, T, T], [T, F, F], [T, F, T]],
            [[T, T, F], [T, T, T], [F, F, F]]
        ]

        self.input_2 = '00011011b'
        self.array_2 = [
            [[F, F], [F, T]],
            [[T, F], [T, T]]
        ]

    def test_build_block(self):
        bits = self.input_3[:-1]
        dims = [3, 3, 3]
        assert_array_equal(self.array_3, build_block(bits, dims))

    def test_parse_cube(self):
        builder = BlockBuilder()
        blocks = builder.parse_blocks(f"({self.input_3}; {self.input_2})")
        assert_array_equal(self.array_3, next(blocks))
        assert_array_equal(self.array_2, next(blocks))

    def test_parse_rect(self):
        builder = BlockBuilder([2, 4])
        blocks = builder.parse_blocks(f"{self.input_2} {self.input_3}")
        array_a = [[F, F, F, T], [T, F, T, T]]
        assert_array_equal(array_a, next(blocks))
        array_b = [[F, F, F, F], [F, T, F, T]]  # truncated input_3
        assert_array_equal(array_b, next(blocks))

    @staticmethod
    def test_parse_short_input():
        builder = BlockBuilder([2, 3, 4])
        blocks = builder.parse_blocks('0010101b')  # will replicate input
        a = [[[F, F, T, F], [T, F, T, F], [F, T, F, T]],
             [[F, T, F, F], [T, F, T, F], [T, F, F, T]]]
        assert_array_equal(a, next(blocks))

class StandardViewpointTest(unittest.TestCase):

    def __init__(self, *args, **kwargs):
        super(StandardViewpointTest, self).__init__(*args, **kwargs)
        self.projector = Projector()
        self.block = build_block('000000101000110111100000000b', [3, 3, 3])

    def create_viewpoint(self, view_name):
        return self.projector.get_view(view_name)

    def test_get_view(self):
        self.assert_view([[T, F, F], [T, T, T], [T, F, T]], TOP_VIEW)
        self.assert_view([[T, F, F], [T, T, F], [T, T, T]], FRONT_VIEW)
        self.assert_view([[T, F, F], [F, T, F], [F, T, T]], RIGHT_VIEW)
        self.assert_view([[T, F, T], [T, T, T], [T, F, F]], BOTTOM_VIEW)
        self.assert_view([[F, F, T], [F, T, T], [T, T, T]], BACK_VIEW)
        self.assert_view([[F, F, T], [F, T, F], [T, T, F]], LEFT_VIEW)

    def test_projector(self):
        i = self.projector.get_images(self.block, TOP_VIEW)
        image = next(i)
        self.assertEqual("XX  XX\nXXXXXX\nXX    ", image)

    def assert_view(self, expected_view, view_name):
        viewpoint = self.create_viewpoint(view_name)
        view = viewpoint.project(self.block)
        assert_array_equal(view, expected_view)

class ProjectorTest(unittest.TestCase):
    def test_parse_aspect(self):
        projector = Projector()
        self.assertRaises(RuntimeError, projector.parse_aspect, None)
        self.assertRaises(RuntimeError, projector.parse_aspect, '')
        self.assertRaises(RuntimeError, projector.parse_aspect, 'a b+')
        self.assertEqual(projector.parse_aspect('foo'), ('foo', None))
        self.assertEqual(projector.parse_aspect('bar-1'), ('bar', '-1'))
        self.assertEqual(projector.parse_aspect('baz+123'), ('baz', '+123'))

if __name__ == '__main__':
    unittest.main()


text:
import numpy as np

def label2rgb(label_field, image, kind='avg', bg_label=-1, bg_color=(0, 0, 0)):

    #std_list = list()
    out = np.zeros_like(image)
    labels = np.unique(label_field)
    bg = (labels == bg_label)
    if bg.any():
        labels = labels[labels != bg_label]
        mask = (label_field == bg_label).nonzero()
        out[mask] = bg_color
    for label in labels:
        mask = (label_field == label).nonzero()
        #std = np.std(image[mask])
        #std_list.append(std)
        if kind == 'avg':
            color = image[mask].mean(axis=0)
        elif kind == 'median':
            color = np.median(image[mask], axis=0)
        elif kind == 'mix':
            std = np.std(image[mask])
            if std < 20:
                color = image[mask].mean(axis=0)
            elif 20 < std < 40:
                mean = image[mask].mean(axis=0)
                median = np.median(image[mask], axis=0)
                color = 0.5*mean + 0.5*median
            elif 40 < std:
                color = np.median(image[mask], axis=0)
        out[mask] = color
    return out


text:
from functions import functions_roc_curves
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
from sklearn.metrics import average_precision_score
import os
import pymongo
from pymongo import MongoClient

def main():
	#Connection to Mongo DB
	try:
		conn=pymongo.MongoClient()
		print "Connected successfully!!!\n"
	except pymongo.errors.ConnectionFailure, e:
		print "Could not connect to MongoDB: %s" % e

	#Usar base de datos
	db = conn.DUDE
	#Usar una coleccion
	collection = db.usrcat
	coll = db.electroshape

	cursor_usrcat = collection.find()
	cursor_electroshape = coll.find()

	try:
		for i in cursor_electroshape:
			for mol in i["similarities"]:
				if mol["molid"].startswith("CHEMBL"):
					molid_ref = "CHEMBL"
					break
				elif mol["molid"].startswith("ZINC"):
					molid_ref = "ZINC"
					break
				else:
					continue

		cursor_electroshape.rewind()

		# Calcul lists with usrcat and electroshape
		if molid_ref == "CHEMBL":
			usrcat_list = np.array([(1, mol["similarity"]) if mol["molid"].startswith("CHEMBL") else (0, mol["similarity"]) for i in cursor_usrcat for mol in i["similarities"] ])
			electroshape_list = np.array([(1, mol["similarity"]) if mol["molid"].startswith("CHEMBL") else (0, mol["similarity"]) for i in cursor_electroshape for mol in i["similarities"] ])
		elif molid_ref == "ZINC":
			usrcat_list = np.array([(0, mol["similarity"]) if mol["molid"].startswith("ZINC") else (1, mol["similarity"]) for i in cursor_usrcat for mol in i["similarities"] ])
			electroshape_list = np.array([(0, mol["similarity"]) if mol["molid"].startswith("ZINC") else (1, mol["similarity"]) for i in cursor_electroshape for mol in i["similarities"] ])

		y_test, y_score = np.hsplit(usrcat_list, 2)
		print average_precision_score(y_test, y_score)

		y_test_e, y_score_e = np.hsplit(electroshape_list, 2)
		print average_precision_score(y_test_e, y_score_e)

		# Compute ROC curve and ROC area for the target for usrcat
		fpr, tpr, roc_auc = functions_roc_curves.curve(usrcat_list)

		# Compute ROC curve and ROC area for the target for electroshape
		fpr_e, tpr_e, roc_auc_e = functions_roc_curves.curve(electroshape_list)

		# List roc curves for representation
		list_roc_curves = [];
		_label = ""

		_label="ROC average curve for usrcat"
		list_roc_curves.append((fpr,tpr,roc_auc,_label))

		_label="ROC average curve for electroshape"
		list_roc_curves.append((fpr_e,tpr_e,roc_auc_e,_label))

		# Plot of a ROC curve for a specific target in usrcat and electroshape
		path = "./graphics/"
		target = "Average"

		#functions_roc_curves.graphic(target, path, list_roc_curves)

	except Exception as e:
		print "error"

if __name__ == '__main__':
	main()

max_stars_repo_name: albags/VirtualScreeningOnDecoys

max_issues_repo_name: albags/VirtualScreeningOnDecoys

max_forks_repo_name: albags/VirtualScreeningOnDecoys


text:
from lensit.clusterlens.profile import  profile
import numpy as np
from lensit import ell_mat, pbs
from lensit.sims import ffs_phas, ffs_cmbs, ffs_maps
from lensit.misc.misc_utils import enumerate_progress, gauss_beam
from lensit.ffs_deflect import ffs_deflect
import lensit as li
from os.path import join as opj
import os
from camb import CAMBdata

def get_cluster_libdir(cambinifile, profilename, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp):
    return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename,
    'npix%s_lpix_%samin_lmaxsky%s' % (npix, lpix_amin, ellmax_sky), 'M200_%.6E_z%s' % (M200, z), '%s_sims' % nsims)

class cluster_maps(object):
    def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, cmb_exp='5muKamin_1amin', cache_maps=False):
        """Library for flat-sky CMB simulations lensed by a galaxy cluster.

        Args:
            lib_dir: various things will be cached there
            npix: number of pixels on one side of the square box
            lpix_amin: physical size (in arcmin) of the pixels
            nsims: number of CMB maps to simulate
            cosmo: Instantiated camb.results.CAMBdata object
            profparams: dict containing the parameters defining the profile
            profilename: string defining the density profile of the cluster (e.g. nfw)
            ellmax_sky: maximum multipole of CMB spectra used to generate the CMB maps
    """
        self.libdir = libdir
        self.cosmo = cosmo
        if profilename == 'nfw':
            self.M200 = profparams['M200c']
            self.z = profparams['z']
        self.npix = npix
        self.lpix_amin = lpix_amin
        shape = (self.npix, self.npix)
        lpix_rad = self.lpix_amin * np.pi / 180 / 60
        lsides = (lpix_rad*self.npix, lpix_rad*self.npix)
        ellmatdir = opj(li._get_lensitdir()[0], 'temp', 'ellmats', 'ellmat_npix%s_lpix_%samin' % (self.npix, self.lpix_amin))
        self.ellmat = ell_mat.ell_mat(ellmatdir, shape, lsides)

        assert ellmax_sky < cosmo.Params.max_l, "Ask for a higher ellmax in CAMB, see l_max_scalar in param inifile or camb.CAMBparams.set_for_lmax()"
        self.ellmax_sky = ellmax_sky

        num_threads=int(os.environ.get('OMP_NUM_THREADS', 1))
        self.lib_skyalm =  ell_mat.ffs_alm_pyFFTW(self.ellmat, num_threads=num_threads, filt_func=lambda ell: ell <= self.ellmax_sky)

        nfields = 3

        camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T
        cls_cmb = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]}

        # Generate the CMB random phases
        skypha_libdir = opj(self.libdir,  'len_alms', 'skypha')
        skypha = ffs_phas.ffs_lib_phas(skypha_libdir, nfields, self.lib_skyalm, nsims_max=nsims)
        if not skypha.is_full() and pbs.rank == 0:
            for i, idx in enumerate_progress(np.arange(nsims, dtype=int), label='Generating CMB phases'):
                skypha.get_sim(int(idx))
        pbs.barrier()

        self.haloprofile = profile(self.cosmo, profilename)

        lencmbs_libdir = opj(self.libdir, 'len_alms')

        # Get the noise and beam of the experiment
        sN_uKamin, sN_uKaminP, Beam_FWHM_amin, ellmin, ellmax = li.get_config(cmb_exp)
        self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, cls_cmb, self.M200, self.z, self.haloprofile, lib_pha=skypha)

        # The lmax of the data maps is set by the pixel size
        lmax_data = self.ellmat._get_ellmax()
        cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=lmax_data)

        # The resolution of the data map could be lower than the resolution of the sky map (=the true map)
        self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= lmax_data, num_threads=num_threads)

        vcell_amin2 = np.prod(self.lib_datalm.ell_mat.lsides) / np.prod(self.lib_datalm.ell_mat.shape) * (180 * 60. / np.pi) ** 2
        nTpix = sN_uKamin / np.sqrt(vcell_amin2)
        nPpix = sN_uKaminP / np.sqrt(vcell_amin2)

        # Generate the noise random phases
        pixpha_libdir = opj(self.libdir, 'pixpha')
        pixpha = ffs_phas.pix_lib_phas(pixpha_libdir, 3, self.lib_datalm.ell_mat.shape, nsims_max=nsims)
        if not pixpha.is_full() and pbs.rank == 0:
            for _i, idx in enumerate_progress(np.arange(nsims), label='Generating Noise phases'):
                pixpha.get_sim(idx)
        pbs.barrier()
        maps_libdir = opj(self.libdir, 'maps')

        self.maps_lib = ffs_maps.lib_noisemap(maps_libdir, self.lib_datalm, self.len_cmbs, cl_transf, nTpix, nPpix, nPpix,
                                        pix_pha=pixpha, cache_sims=cache_maps)

    def get_unl_map(self, idx, field='t'):
        """Get the unlensed CMB map
        Args:
            idx: index of the simulation to return
            field: 't', 'e' or 'b' for temperature, E and B modes of polarization respectively
        Returns:
            unlensed map: numpy array of shape self.lib_skyalm.shape
        """
        if field in ['t', 'e', 'b']:
            return self.lib_skyalm.alm2map(self.len_cmbs.unlcmbs.get_sim_alm(idx, field))
        elif field in ['q', 'u']:
            i = 0 if field =='q' else 1
            return self.lib_skyalm.alm2map(self.len_cmbs.unlcmbs.get_sim_qulm(idx)[i])

    def get_len_map(self, idx, field='t'):
        """Get the lensed CMB map
        Args:
            idx: index of the simulation to return
            field: 't', 'e' or 'b' for temperature, E and B modes of polarization respectively
        Returns:
            lensed map: numpy array of shape self.lib_skyalm.shape
        """
        return self.lib_skyalm.alm2map(self.len_cmbs.get_sim_alm(idx, field))

    def get_obs_map(self, idx, field='t'):
        """Get the observed CMB map, inclusing the noise and the bean of the instrument
        Args:
            idx: index of the simulation to return
            field: 't', 'e', 'b', 'q' or 'u' for temperature, E, B, Q and U modes of polarization respectively
        Returns:
            lensed map: numpy array of shape self.lib_datalm.shape
        """
        if field=='t':
            return self.maps_lib.get_sim_tmap(idx)
        elif field=='q':
            return self.maps_lib.get_sim_qumap(idx)[0]
        elif field =='u':
            return self.maps_lib.get_sim_qumap(idx)[1]
        elif field in ['e', 'b']:
            qmap, umap = self.maps_lib.get_sim_qumap(idx)
            qlm = self.lib_datalm.map2alm(qmap)
            ulm = self.lib_datalm.map2alm(umap)
            elm, blm = self.lib_datalm.QUlms2EBalms(np.array([qlm, ulm]))
            if field == 'e':
                return self.lib_datalm.alm2map(elm)
            elif field == 'b':
                return self.lib_datalm.alm2map(blm)

    def get_noise_map(self, idx, field='t'):
        """Get the noise map
        Args:
            idx: index of the simulation to return
            field: 't', 'q' or 'u' for temperature, Q and U modes of polarization respectively
        Returns:
            noise map: numpy array of shape self.lib_datalm.shape
        """
        assert field in ['t', 'q', 'u'], "The noise maps are generated for t, q and u fields, not {}".format(field)
        if field =='t':
            return self.maps_lib.get_noise_sim_tmap(idx)
        elif field =='q':
            return self.maps_lib.get_noise_sim_qmap(idx)
        elif field =='u':
            return self.maps_lib.get_noise_sim_umap(idx)

    def get_kappa_map(self, M200, z):
        """Get the convergence map of the cluster
        Args:
            M200: mass (in Msol) of the cluster in a sphere of density
                200 times the critical density of the universe
            z: redshift of the cluster
        Returns:
            convergence map: numpy array of shape self.lib_skyalm.shape
        """
        return self.haloprofile.kappa_map(M200, z, self.lib_skyalm.shape, self.lib_skyalm.lsides)

class sim_cmb_len_cluster(ffs_cmbs.sims_cmb_len):
    def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_cmb:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False):
        super(sim_cmb_len_cluster, self).__init__(lib_dir, lib_skyalm, cls_cmb, lib_pha=lib_pha, use_Pool=0, cache_lens=False)

        kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides)
        self.defl_map = haloprofile.kmap2deflmap(kappa_map, lib_skyalm.shape, lib_skyalm.lsides)

        assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_cmb dictionnary to avoid errors."

    def _get_f(self, idx=None):
        """We refine the displacement field using the convergence map of the cluster"""
        return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides)

max_stars_repo_name: samuelsimko/LensIt

max_issues_repo_name: samuelsimko/LensIt

max_forks_repo_name: samuelsimko/LensIt


text:
#Ridge regression individually on each of the 5 target variables
def my_score(model,X,y_true):
    import numpy as np
    y_true = np.array(y_true)
    y_pred=model.predict(X)
    return np.mean(np.sum(np.abs(y_true - y_pred), axis=0)/np.sum(y_true, axis=0))
from sklearn.linear_model import Ridge
alpha=0.17
m=Ridge(alpha,normalize=True)
m.fit(train_df[:3634],y_train_df[:3634]['age'])
np.round(my_score(m,train_df[3634:],y_train_df[3634:]['age']),5)
#FOR AGE gives 0.1483 score

max_stars_repo_name: nizamphoenix/kaggle

max_issues_repo_name: nizamphoenix/kaggle

max_forks_repo_name: nizamphoenix/kaggle


text:
#!/usr/bin/env python2

from __future__ import print_function
import numpy as np
import matplotlib
matplotlib.use('PDF')
import matplotlib.pyplot as plt
import sys
import csv

from optparse import OptionParser

class CmdLineParser(object):
	def __init__(self):
		self.parser = OptionParser(usage='usage: python2 parserScores.py scores.csv outputfile.csv')

	def addOption(self, *args, **kwargs):
		self.parser.add_option(*args, **kwargs)

	def parseArgs(self):
		(options, args) = self.parser.parse_args()
		options.input = args[0]
		options.output = args[1]
		return options

def parseScores(fp, outputfile):
	"""
	read a CSV of scores
	the format is as follows:
	instance, seed, model, score, time, run_number
	"""
	firstRow = [0] # n_models, model1, model2, ...
	instances = [] # name of the instances(with the seed)
	runs = []	# rows

	for row in fp:
		row = row.strip().split(', ')
		instance = row[0] + "_" + row[5]
		model = row[2]
		time = row[4]

		try: #search for the model's index
			model_index = firstRow.index(model)-1
		except ValueError:
			#print("new model found", model)
			firstRow.append(model)
			firstRow[0] += 1
			model_index = firstRow[0]-1

		try: #search for the instance's index
			instance_index = instances.index(instance)
			l = len(runs[instance_index])
			if (model_index >= l):
				# add extra cells
				runs[instance_index].extend([None]* (model_index-l+1))
			runs[instance_index][model_index] = time
		except ValueError:
			#print("new instance:", instance)
			instances.append(instance)
			instance_index = instances.index(instance)
			runs.append([time])

	# create table
	with open(outputfile, mode='a') as file:
		writer = csv.writer(file, delimiter=',')
		writer.writerow(firstRow)
		# write runs
		for i in range(len(instances)):
			writer.writerow([instances[i]] + runs[i])

def main():
	parser = CmdLineParser()
	opt = parser.parseArgs()
	parseScores(open(opt.input, 'r'), opt.output)

if __name__ == '__main__':
	main()

max_stars_repo_name: mikelemilia/travelling-salesman-problem

max_issues_repo_name: mikelemilia/travelling-salesman-problem

max_forks_repo_name: mikelemilia/travelling-salesman-problem


text:
# -*- coding: utf-8 -*-
#
# test_stdp_triplet_synapse.py
#
# This file is part of NEST.
#
# Copyright (C) 2004 The NEST Initiative
#
# NEST is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# NEST is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NEST.  If not, see <http://www.gnu.org/licenses/>.

# This script tests the stdp_triplet_synapse in NEST.

import nest
import unittest
from math import exp
import numpy as np

@nest.check_stack
class STDPTripletConnectionTestCase(unittest.TestCase):
    """Check stdp_triplet_connection model properties."""

    def setUp(self):
        nest.set_verbosity('M_WARNING')
        nest.ResetKernel()

        # settings
        self.dendritic_delay = 1.0
        self.decay_duration = 5.0
        self.synapse_model = "stdp_triplet_synapse"
        self.syn_spec = {
            "model": self.synapse_model,
            "delay": self.dendritic_delay,
            # set receptor 1 post-synaptically, to not generate extra spikes
            "receptor_type": 1,
            "weight": 5.0,
            "tau_plus": 16.8,
            "tau_plus_triplet": 101.0,
            "Aplus": 0.1,
            "Aminus": 0.1,
            "Aplus_triplet": 0.1,
            "Aminus_triplet": 0.1,
            "Kplus": 0.0,
            "Kplus_triplet": 0.0,
            "Wmax": 100.0,
        }
        self.post_neuron_params = {
            "tau_minus": 33.7,
            "tau_minus_triplet": 125.0,
        }

        # setup basic circuit
        self.pre_neuron = nest.Create("parrot_neuron")
        self.post_neuron = nest.Create(
            "parrot_neuron", 1, params=self.post_neuron_params)
        nest.Connect(self.pre_neuron, self.post_neuron, syn_spec=self.syn_spec)

    def generateSpikes(self, neuron, times):
        """Trigger spike to given neuron at specified times."""
        delay = 1.
        gen = nest.Create("spike_generator", 1, {
                          "spike_times": [t - delay for t in times]})
        nest.Connect(gen, neuron, syn_spec={"delay": delay})

    def status(self, which):
        """Get synapse parameter status."""
        stats = nest.GetConnections(
            self.pre_neuron, synapse_model=self.synapse_model)
        return nest.GetStatus(stats, [which])[0][0]

    def decay(self, time, Kplus, Kplus_triplet, Kminus, Kminus_triplet):
        """Decay variables."""
        Kplus *= exp(-time / self.syn_spec["tau_plus"])
        Kplus_triplet *= exp(-time / self.syn_spec["tau_plus_triplet"])
        Kminus *= exp(-time / self.post_neuron_params["tau_minus"])
        Kminus_triplet *= exp(-time /
                              self.post_neuron_params["tau_minus_triplet"])
        return (Kplus, Kplus_triplet, Kminus, Kminus_triplet)

    def facilitate(self, w, Kplus, Kminus_triplet):
        """Facilitate weight."""
        Wmax = self.status("Wmax")
        return np.sign(Wmax) * (abs(w) + Kplus * (
            self.syn_spec["Aplus"] +
            self.syn_spec["Aplus_triplet"] * Kminus_triplet)
        )

    def depress(self, w, Kminus, Kplus_triplet):
        """Depress weight."""
        Wmax = self.status("Wmax")
        return np.sign(Wmax) * (abs(w) - Kminus * (
            self.syn_spec["Aminus"] +
            self.syn_spec["Aminus_triplet"] * Kplus_triplet)
        )

    def assertAlmostEqualDetailed(self, expected, given, message):
        """Improve assetAlmostEqual with detailed message."""
        messageWithValues = "%s (expected: `%s` was: `%s`" % (
            message, str(expected), str(given))
        self.assertAlmostEqual(given, expected, msg=messageWithValues)

    def test_badPropertiesSetupsThrowExceptions(self):
        """Check that exceptions are thrown when setting bad parameters."""

        def setupProperty(property):
            bad_syn_spec = self.syn_spec.copy()
            bad_syn_spec.update(property)
            nest.Connect(self.pre_neuron, self.post_neuron,
                         syn_spec=bad_syn_spec)

        def badPropertyWith(content, parameters):
            self.assertRaisesRegexp(
                nest.NESTError, "BadProperty(.+)" + content,
                setupProperty, parameters
            )

        badPropertyWith("Kplus", {"Kplus": -1.0})
        badPropertyWith("Kplus_triplet", {"Kplus_triplet": -1.0})

    def test_varsZeroAtStart(self):
        """Check that pre and post-synaptic variables are zero at start."""
        self.assertAlmostEqualDetailed(
            0.0, self.status("Kplus"), "Kplus should be zero")
        self.assertAlmostEqualDetailed(0.0, self.status(
            "Kplus_triplet"), "Kplus_triplet should be zero")

    def test_preVarsIncreaseWithPreSpike(self):
        """Check that pre-synaptic variables (Kplus, Kplus_triplet) increase
        after each pre-synaptic spike."""

        self.generateSpikes(self.pre_neuron, [2.0])

        Kplus = self.status("Kplus")
        Kplus_triplet = self.status("Kplus_triplet")

        nest.Simulate(20.0)
        self.assertAlmostEqualDetailed(
            Kplus + 1.0,
            self.status("Kplus"),
            "Kplus should have increased by 1")
        self.assertAlmostEqualDetailed(
            Kplus_triplet + 1.0,
            self.status("Kplus_triplet"),
            "Kplus_triplet should have increased by 1")

    def test_preVarsDecayAfterPreSpike(self):
        """Check that pre-synaptic variables (Kplus, Kplus_triplet) decay
        after each pre-synaptic spike."""

        self.generateSpikes(self.pre_neuron, [2.0])
        # trigger computation
        self.generateSpikes(self.pre_neuron, [2.0 + self.decay_duration])

        (Kplus, Kplus_triplet, _, _) = self.decay(
            self.decay_duration, 1.0, 1.0, 0.0, 0.0)
        Kplus += 1.0
        Kplus_triplet += 1.0

        nest.Simulate(20.0)
        self.assertAlmostEqualDetailed(
            Kplus, self.status("Kplus"), "Kplus should have decay")
        self.assertAlmostEqualDetailed(Kplus_triplet, self.status(
            "Kplus_triplet"), "Kplus_triplet should have decay")

    def test_preVarsDecayAfterPostSpike(self):
        """Check that pre-synaptic variables (Kplus, Kplus_triplet) decay
        after each post-synaptic spike."""

        self.generateSpikes(self.pre_neuron, [2.0])
        self.generateSpikes(self.post_neuron, [3.0, 4.0])
        # trigger computation
        self.generateSpikes(self.pre_neuron, [2.0 + self.decay_duration])

        (Kplus, Kplus_triplet, _, _) = self.decay(
            self.decay_duration, 1.0, 1.0, 0.0, 0.0)
        Kplus += 1.0
        Kplus_triplet += 1.0

        nest.Simulate(20.0)
        self.assertAlmostEqualDetailed(
            Kplus, self.status("Kplus"), "Kplus should have decay")
        self.assertAlmostEqualDetailed(Kplus_triplet, self.status(
            "Kplus_triplet"), "Kplus_triplet should have decay")

    def test_weightChangeWhenPrePostSpikes(self):
        """Check that weight changes whenever a pre-post spike pair happen."""

        self.generateSpikes(self.pre_neuron, [2.0])
        self.generateSpikes(self.post_neuron, [4.0])
        self.generateSpikes(self.pre_neuron, [6.0])  # trigger computation

        Kplus = self.status("Kplus")
        Kplus_triplet = self.status("Kplus_triplet")
        Kminus = 0.0
        Kminus_triplet = 0.0
        weight = self.status("weight")
        Wmax = self.status("Wmax")

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0, Kplus, Kplus_triplet, Kminus, Kminus_triplet)
        weight = self.depress(weight, Kminus, Kplus_triplet)
        Kplus += 1.0
        Kplus_triplet += 1.0

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0 + self.dendritic_delay, Kplus, Kplus_triplet,
            Kminus, Kminus_triplet
        )
        weight = self.facilitate(weight, Kplus, Kminus_triplet)
        Kminus += 1.0
        Kminus_triplet += 1.0

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0 - self.dendritic_delay, Kplus, Kplus_triplet,
            Kminus, Kminus_triplet
        )
        weight = self.depress(weight, Kminus, Kplus_triplet)

        nest.Simulate(20.0)
        self.assertAlmostEqualDetailed(weight, self.status(
            "weight"), "weight should have decreased")

    def test_weightChangeWhenPrePostPreSpikes(self):
        """Check that weight changes whenever a pre-post-pre spike triplet
        happen."""

        self.generateSpikes(self.pre_neuron, [2.0, 6.0])
        self.generateSpikes(self.post_neuron, [4.0])
        self.generateSpikes(self.pre_neuron, [8.0])  # trigger computation

        Kplus = self.status("Kplus")
        Kplus_triplet = self.status("Kplus_triplet")
        Kminus = 0.0
        Kminus_triplet = 0.0
        weight = self.status("weight")
        Wmax = self.status("Wmax")

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0, Kplus, Kplus_triplet, Kminus, Kminus_triplet)
        weight = self.depress(weight, Kminus, Kplus_triplet)
        Kplus += 1.0
        Kplus_triplet += 1.0

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0 + self.dendritic_delay, Kplus, Kplus_triplet,
            Kminus, Kminus_triplet
        )

        weight = self.facilitate(weight, Kplus, Kminus_triplet)
        Kminus += 1.0
        Kminus_triplet += 1.0

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0 - self.dendritic_delay, Kplus, Kplus_triplet,
            Kminus, Kminus_triplet
        )
        weight = self.depress(weight, Kminus, Kplus_triplet)
        Kplus += 1.0
        Kplus_triplet += 1.0

        (Kplus, Kplus_triplet, Kminus, Kminus_triplet) = self.decay(
            2.0, Kplus, Kplus_triplet, Kminus, Kminus_triplet)
        weight = self.depress(weight, Kminus, Kplus_triplet)

        nest.Simulate(20.0)
        self.assertAlmostEqualDetailed(weight, self.status(
            "weight"), "weight should have decreased")

    def test_maxWeightStaturatesWeight(self):
        """Check that setting maximum weight property keep weight limited."""

        limited_weight = self.status("weight") + 1e-10
        limited_syn_spec = self.syn_spec.copy()
        limited_syn_spec.update({"Wmax": limited_weight})
        nest.Connect(self.pre_neuron, self.post_neuron,
                     syn_spec=limited_syn_spec)

        self.generateSpikes(self.pre_neuron, [2.0])
        self.generateSpikes(self.pre_neuron, [3.0])  # trigger computation

        nest.Simulate(20.0)
        self.assertAlmostEqualDetailed(limited_weight, self.status(
            "weight"), "weight should have been limited")

@nest.check_stack
class STDPTripletInhTestCase(STDPTripletConnectionTestCase):

    def setUp(self):
        nest.set_verbosity('M_WARNING')
        nest.ResetKernel()

        # settings
        self.dendritic_delay = 1.0
        self.decay_duration = 5.0
        self.synapse_model = "stdp_triplet_synapse"
        self.syn_spec = {
            "model": self.synapse_model,
            "delay": self.dendritic_delay,
            # set receptor 1 post-synaptically, to not generate extra spikes
            "receptor_type": 1,
            "weight": -5.0,
            "tau_plus": 16.8,
            "tau_plus_triplet": 101.0,
            "Aplus": 0.1,
            "Aminus": 0.1,
            "Aplus_triplet": 0.1,
            "Aminus_triplet": 0.1,
            "Kplus": 0.0,
            "Kplus_triplet": 0.0,
            "Wmax": -100.0,
        }
        self.post_neuron_params = {
            "tau_minus": 33.7,
            "tau_minus_triplet": 125.0,
        }

        # setup basic circuit
        self.pre_neuron = nest.Create("parrot_neuron")
        self.post_neuron = nest.Create("parrot_neuron", 1,
                                       params=self.post_neuron_params)
        nest.Connect(self.pre_neuron, self.post_neuron, syn_spec=self.syn_spec)

def suite_inh():
    return unittest.makeSuite(STDPTripletInhTestCase, "test")

def suite():
    return unittest.makeSuite(STDPTripletConnectionTestCase, "test")

def run():
    runner = unittest.TextTestRunner(verbosity=2)
    runner.run(suite())
    runner.run(suite_inh())

if __name__ == "__main__":
    run()

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
"""
 *
 * Copyright (c) 2021 Manuel Yves Galliker
 *               2021 Autonomous Systems Lab ETH Zurich
 * All rights reserved.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 * 3. Neither the name Data Driven Dynamics nor the names of its contributors may be
 *    used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *

The model class contains properties shared between all models and shgall simplyfy automated checks and the later
export to a sitl gazebo model by providing a unified interface for all models. """

__author__ = "Manuel Yves Galliker"
__maintainer__ = "Manuel Yves Galliker"
__license__ = "BSD 3"

from progress.bar import Bar
import pandas as pd
import math
import time
import yaml
import numpy as np
import matplotlib.pyplot as plt
import os
from src.models.model_config import ModelConfig
from src.tools.ulog_tools import load_ulog, pandas_from_topic
from src.tools.dataframe_tools import compute_flight_time, resample_dataframe_list
from src.tools.quat_utils import quaternion_to_rotation_matrix

class DataHandler(object):

    visual_dataframe_selector_config_dict = {
        "x_axis_col": "timestamp",
        "sub_plt1_data": ["q0", "q1", "q2", "q3"],
        "sub_plt2_data": ["u0", "u1", "u2", "u3"]}

    def __init__(self, config_file):
        print("===============================================================================")
        print("                              Data Processing                                  ")
        print("===============================================================================")
        self.config = ModelConfig(config_file)
        config_dict = self.config.dynamics_model_config

        assert type(
            config_dict) is dict, 'req_topics_dict input must be a dict'
        assert bool(config_dict), 'req_topics_dict can not be empty'
        self.config_dict = config_dict
        self.resample_freq = config_dict["resample_freq"]
        self.estimate_angular_acceleration = config_dict["estimate_angular_acceleration"]
        print("Resample frequency: ", self.resample_freq, "Hz")
        self.req_topics_dict = config_dict["data"]["required_ulog_topics"]

        self.req_dataframe_topic_list = config_dict["data"]["req_dataframe_topic_list"]

        self.estimate_forces = config_dict["estimate_forces"]
        self.estimate_moments = config_dict["estimate_moments"]

        # used to generate a dict with the resulting coefficients later on.
        self.coef_name_list = []
        self.result_dict = {}

    def loadLogs(self, rel_data_path):
        self.rel_data_path = rel_data_path
        if (os.path.isdir(rel_data_path)):
            self.data_df = pd.DataFrame()
            for filename in os.listdir(rel_data_path):
                self.loadLogFile(os.path.join(rel_data_path, filename))

        else:
            if not self.loadLogFile(rel_data_path):
                raise TypeError("File extension needs to be either csv or ulg")

    def loadLogFile(self, rel_data_path):

        if (rel_data_path.endswith(".csv")):
            print("Loading CSV file: ", rel_data_path)
            self.data_df = pd.read_csv(rel_data_path, index_col=0)
            print("Loading topics: ", self.req_dataframe_topic_list)
            for req_topic in self.req_dataframe_topic_list:
                assert(
                    req_topic in self.data_df), ("missing topic in loaded csv: " + str(req_topic))
            return True

        elif (rel_data_path.endswith(".ulg")):
            print("Loading uLog file: ", rel_data_path)
            ulog = load_ulog(rel_data_path)
            print("Loading topics:")
            for req_topic in self.req_topics_dict:
                print(req_topic)
            self.check_ulog_for_req_topics(ulog)
            self.data_df = self.compute_resampled_dataframe(ulog)
            return True

        else:
            return False

    def check_ulog_for_req_topics(self, ulog):
        for topic_type in self.req_topics_dict.keys():
            try:
                topic_type_data = ulog.get_dataset(topic_type)
            except:
                print("Missing topic type: ", topic_type)
                exit(1)
            topic_type_data = topic_type_data.data
            ulog_topic_list = self.req_topics_dict[topic_type]["ulog_name"]
            for topic_index in range(len(ulog_topic_list)):
                try:
                    topic = ulog_topic_list[topic_index]
                    topic_data = (topic_type_data[topic])
                except:
                    print("Missing topic: ", topic_type,
                          ulog_topic_list[topic_index])
                    exit(1)
        return

    def compute_resampled_dataframe(self, ulog):
        print("Starting data resampling of topic types: ",
              self.req_topics_dict.keys())
        # setup object to crop dataframes for flight data
        df_list = []
        topic_type_bar = Bar('Resampling', max=len(
            self.req_topics_dict.keys()))

        # getting data
        fts = []
        for topic_type in self.req_topics_dict.keys():
            topic_dict = self.req_topics_dict[topic_type]

            if "id" in topic_dict.keys():
                id = topic_dict["id"]
                curr_df = pandas_from_topic(ulog, [topic_type], id)
            else:
                curr_df = pandas_from_topic(ulog, [topic_type])

            if topic_type == "actuator_outputs":
                fts = compute_flight_time(curr_df)
            elif topic_type == "actuator_controls_0":
                fts = compute_flight_time(curr_df)

            curr_df = curr_df[topic_dict["ulog_name"]]
            if "dataframe_name" in topic_dict.keys():
                assert (len(topic_dict["dataframe_name"]) == len(topic_dict["ulog_name"])), (
                    'could not rename topics of type', topic_type, "due to rename list not having an entry for every topic.")
                curr_df.columns = topic_dict["dataframe_name"]
            topic_type_bar.next()
            if (topic_type == "vehicle_angular_velocity" and self.estimate_angular_acceleration):
                ang_vel_mat = curr_df[["ang_vel_x",
                                      "ang_vel_y",  "ang_vel_z"]].to_numpy()
                time_in_secods_np = (curr_df[["timestamp"]].to_numpy()/1000000)
                time_in_secods_np = time_in_secods_np.flatten()
                ang_acc_np = np.gradient(
                    ang_vel_mat, time_in_secods_np, axis=0)
                topic_type_bar.next()
                curr_df[["ang_acc_b_x", "ang_acc_b_y",
                         "ang_acc_b_z"]] = ang_acc_np

            df_list.append(curr_df)

        topic_type_bar.finish()

        # Check if actuator topics are empty
        if not fts :
            print("could not select flight time due to missing actuator topic")
            exit(1)
        # set start and end time of flight duration

        resampled_df = resample_dataframe_list(
            df_list, fts, self.resample_freq)
        if (self.estimate_angular_acceleration):
            ang_vel_mat = resampled_df[["ang_vel_x",
                                       "ang_vel_y",  "ang_vel_z"]].to_numpy()
            for i in range(3):
                ang_vel_mat[:, i] = np.convolve(
                    ang_vel_mat[:, i], np.ones(33), mode='same') / 33

            # Alternate forward differentiation version
            # ang_vel_mat_1 = np.roll(ang_vel_mat, -1, axis=0)
            # diff_angular_acc_mat = (
            #     ang_vel_mat_1 - ang_vel_mat) * self.resample_freq
            # resampled_df[["ang_acc_b_x", "ang_acc_b_y",
            #               "ang_acc_b_z"]] = diff_angular_acc_mat

            time_in_secods_np = (
                resampled_df[["timestamp"]].to_numpy()/1000000)
            time_in_secods_np = time_in_secods_np.flatten()
            ang_acc_np = np.gradient(ang_vel_mat, time_in_secods_np, axis=0)
            topic_type_bar.next()
            resampled_df[["ang_acc_b_x", "ang_acc_b_y",
                          "ang_acc_b_z"]] = ang_acc_np
        return resampled_df.dropna()

    def visually_select_data(self, plot_config_dict=None):
        print("===============================================================================")
        print("                           Data Selection Enabled                              ")
        print("===============================================================================")
        from visual_dataframe_selector.data_selector import select_visual_data
        print("Number of data samples before cropping: ",
              self.data_df.shape[0])
        self.data_df = select_visual_data(
            self.data_df, self.visual_dataframe_selector_config_dict)

    def get_dataframes(self):
        return self.data_df

    def visualize_data(self):
        def plot_scatter(ax, title, dataframe_x, dataframe_y, dataframe_z, color='blue'):
            ax.scatter(self.data_df[dataframe_x], self.data_df[dataframe_y],
                       self.data_df[dataframe_z], s=10, facecolor=color, lw=0, alpha=0.1)
            ax.set_title(title)
            ax.set_xlabel(dataframe_x)
            ax.set_ylabel(dataframe_y)
            ax.set_zlabel(dataframe_z)

        num_plots = 2
        fig = plt.figure("Data Visualization")
        ax1 = fig.add_subplot(num_plots, 2, 1, projection='3d')
        plot_scatter(ax1, 'Local Velocity', 'vx', 'vy', 'vz')

        ax2 = fig.add_subplot(num_plots, 2, 2, projection='3d')
        plot_scatter(ax2, 'Body Acceleration', 'acc_b_x',
                     'acc_b_y', 'acc_b_z', 'red')

        ax3 = fig.add_subplot(num_plots, 2, 3, projection='3d')
        plot_scatter(ax3, 'Body Angular Velocity', 'ang_vel_x',
                     'ang_vel_y', 'ang_vel_z', 'red')

        ax4 = fig.add_subplot(num_plots, 2, 4, projection='3d')
        plot_scatter(ax4, 'Body Angular Acceleration',
                     'ang_acc_b_x', 'ang_acc_b_y', 'ang_acc_b_z', 'red')
        plt.show(block=False)

max_stars_repo_name: ethz-asl/data-driven-dynamics

max_issues_repo_name: ethz-asl/data-driven-dynamics

max_forks_repo_name: ethz-asl/data-driven-dynamics