text:
import numpy as np
from normalization import *
from .linear_regression import LinearRegression

class MultivariateLinearRegression(LinearRegression):
    """
        Multivariate Linear Regression with Gradient Descent
            Pros:
                | 1. Scale with number of features n

            Cons:
                | 1. Need to choose alpha
                | 2. Need to iterate
                | 3. Need Feature Scaling
    """

    def __init__(self, feature_scaling="auto"):
        super(MultivariateLinearRegression, self).__init__()
        if feature_scaling == "mean":
            self.fs = MeanNormalization()
        elif feature_scaling == "zscore":
            self.fs = ZScoreNormalization()
        else:
            self.fs = MeanNormalization()

    def fit(self, x, y, alpha=0.1, l=0, tol=0, scale_feature=True):
        """
              Fits a regression model (Theta) on training data

              :param x: numpy array or sparse matrix of shape [n_samples, n_features]
                    Training data
              :param y: numpy array of shape (n_samples,)
                    Target values
              :param alpha: float, default: 0.1
                    Learning rate
              :param l: float, optional - default is zero
                    Lambda value for Regularization term
              :param tol: float, default: 0
                    Tollerance of convergence value
              :param scale_feature: boolean, optional
                    Flag for feature scaling, default is True
              :return self:  returns an instance of self.

        """

        self.fs.flush()
        if scale_feature:
            x = self.fs.scale(x)
        x = self.reshape_training_set(x)

        th = np.zeros(x.shape[1])
        m = len(x)
        beta = alpha / m
        converged = False
        while not converged:
            error = np.dot(x, th) - y
            temp = th - beta * np.dot(error.T, x)
            temp[1:] += l * beta * th[1:]
            diff = abs(temp - th) <= tol
            converged = np.all(diff)

            th = temp

        self._model = th

        return self

    def fit_model(self, x, y, alpha, iterations, l=0.0, scale_feature=True):
        """
            Fits a regression model (Theta) on training data

              :param x: numpy array or sparse matrix of shape [n_samples, n_features]
                  Training data
              :param y: numpy array of shape (n_samples,)
                  Target values
              :param alpha: float
                  Learning rate
              :param iterations: int
                  number of iterations
              :param l: float, optional - default is zero
                  Lambda value for Regularization term
              :param scale_feature: flag for feature scaling
              :return self:  returns an instance of self.

              """

        self.fs.flush()
        if scale_feature:
            x = self.fs.scale(x)

        x = self.reshape_training_set(x)

        th = np.zeros(x.shape[1])
        m = len(x)
        beta = alpha / m
        i = 0

        j_history = np.zeros(iterations)

        while i < iterations:

            error = np.dot(x, th) - y
            th[1:] = (1 - l * beta) * th[1:]
            th -= beta * np.dot(error.T, x)
            j_history[i] = self.compute_cost(x, y, th, l)
            i += 1

        self.j_history = j_history
        self._model = th

        return self

    def predict(self, x, scale_feature=True):
        """
        Predict using the regression model

            :param x: {array-like, sparse matrix}, shape = (n_samples, n_features)
                  Samples.
            :param scale_feature: flag for feature scaling
            :return: target: array, shape = (n_samples,)
                  Returns predicted values.

        """
        if scale_feature:
            x = self.fs.scale(x)
        th = self._model
        # checks whether the example to predict is a vector or a matrix
        if len(x.shape) > 1:
            xn = np.ones((x.shape[0], x.shape[1] + 1))
            xn[:, 1:] = x
        else:
            xn = np.ones(x.shape[0]+1)
            xn[1:] = x
        return np.dot(xn, th)


text:
import streamlit as st
import pandas as pd
import geopandas as gpd
import plotly.express as px
import matplotlib.pyplot as plt
import numpy as np

st.title("Анализ данных авиаиндустрии")

st.markdown("Привет!")
st.markdown("В этом проекте я буду анализировать данные авиаиндустрии: по перелётам, самолётам и аэропортам.")
st.markdown("Это первая часть анализа, выполненная в питоне. Чтобы увидеть вторую часть (R), пройди по ссылке avia-project-knn.heroku.com.")

st.markdown("Список использованных технологий:")
st.markdown("1. Pandas")
st.markdown("2. SQL")
st.markdown("3. GeoPandas")
st.markdown("4. Продвинутый web scraping с помощью Selenium")
st.markdown("5. Streamlit")
st.markdown("6. API")
st.markdown("7. Networkx")
st.markdown("8. Visualisation (plotly)")
st.markdown("Во второй части проекта использованы:")
st.markdown("9. R")
st.markdown("10. Machine Learning (kNN)")
st.markdown("11. Ggplot2, продвинутое использование с многослойными картинками")
st.markdown("12. Дополнение ggiraph к ggplot2")
st.markdown("Всего строк кода порядка 130 в Python и 80 в R")
st.markdown("Весь проект - единое целое")
st.markdown("И, конечно, он крутой!")

st.markdown("## 1. Описание доступа к данным по авиаперелётам")
st.markdown("Мой основной источник данных - OpenSky Network (https://opensky-network.org). Портал щедро предоставил мне доступ к данным, за что им большое спасибо!")
st.markdown("Данные лежат в Impala Shell (это такая оболочка над базой данных Impala). К оболочке можно подключаться и делать к ней SQL-запросы.")
st.markdown("1. Сначала нужно завести аккаунт на OpenSky Network и запросить доступ к данным (https://opensky-network.org/index.php?option=com_users).")
st.markdown("2. Затем нужно подключиться к БД. Для этого нужно запустить в терминале команду `ssh -p 2230 -l USERNAME data.opensky-network.org`. Вместо `USERNAME` нужно написать свой логин (кажется, он совпадает с фамилией, написанной с большой буквы).")
st.markdown("3. Чтобы увидеть таблицы, которые лежат в БД, нужно запустить команду `show tables;`. Получится такой список:")
st.markdown("""`+--------------------------+
| name                     |
+--------------------------+
| acas_data4               |
| allcall_replies_data4    |
| flarm_raw                |
| flights                  |
| flights_data4            |
| identification_data4     |
| operational_status_data4 |
| position_data4           |
| rollcall_replies_data4   |
| sensor_visibility        |
| sensor_visibility_data3  |
| state_vectors            |
| state_vectors_data3      |
| state_vectors_data4      |
| velocity_data4           |
+--------------------------+`""")
st.markdown("4. Я буду пользоваться только одной из них - `flights`. Это данные по всем завершившимся полётам с 1 января 2016 года. Узнать, что там лежит, можно с помощью команды `describe flights;`")
st.markdown("""`+-------------+--------+---------+
| name        | type   | comment |
+-------------+--------+---------+
| icao24      | string |         |
| firstseen   | string |         |
| lastseen    | string |         |
| duration    | bigint |         |
| callsign    | string |         |
| departure   | string |         |
| destination | string |         |
+-------------+--------+---------+`""")
st.markdown("`icao24` - это уникальный 24-битный код трансмиттера, который стоит на каждом уважающем себя воздушном судне. Эти трансмиттеры передают данные о важных характеристиках полёта - местоположении, направлении, скорости, и так далее. Именно благодаря данным с трансмиттеров существует эта БД. Благодаря этим же данным работает какой-нибудь flightradar24.com. Поскольку код уникальный, а трансмиттеры зачастую не заменяются, этот код можно отождествить с конкретным самолётом.")
st.markdown("`duration` - это длительность полёта в минутах.")
st.markdown("`departure` и `destination` - это соответственно аэропорт отправления и аэропорт назначения. Аэропорты для краткости закодированы, и в этих колонках лежат их ICAO-коды (ICAO - это International Civil Aviation Organization, такой международный регулятор/стандартизатор. Они занимаются в том числе унификацией данных - например, общепринятых алиасов аэропортов).")
st.markdown("`firstseen` и `lastseen` - это время начала и окончания полёта.")
st.markdown("`callsign` - это позывной, который авиалиния даёт рейсу. Он не стандартизирован ICAO, поэтому мне он не нужен.")
st.markdown("5. Чтобы достать данные о каком-нибудь полёте, нужно просто исполнить SQL-запрос: `select icao24, firstseen, lastseen, duration, departure, destination from flights limit 1;`")
st.markdown("""`+--------+---------------------+---------------------+----------+-----------+-------------+
| icao24 | firstseen           | lastseen            | duration | departure | destination |
+--------+---------------------+---------------------+----------+-----------+-------------+
| a029b7 | 2017-10-09 20:30:45 | 2017-10-09 21:13:12 | 42       | KFKL      | KJHW        |
+--------+---------------------+---------------------+----------+-----------+-------------+`""")
st.markdown("""Самолёт с трансмиттером `a029b7` вылетел 9 октября 2017 года из `KFKL` в `KJHW` на 42 минуты. При желании можно узнать, что это легкомоторный самолёт (https://opensky-network.org/aircraft-profile?icao24=a029b7), а летел он из
Venango Regional Airport (штат Пенсильвания) в Chautauqua County-Jamestown Airport (штат Нью-Йорк).""")
st.markdown("""6. К сожалению, у меня не получилось интегрировать конкретно эту Impala Shell в питон-скрипт, поэтому придётся делать это руками.
            Буду копировать output терминала, вставлять в файл, а потом парсить средствами питона.""")
st.markdown("7. Теперь я могу собрать нужные мне данные в CSV-файл. Об этом ниже.")

st.markdown("## 2. Выгрузка данных по авиаперелётам")
st.markdown("Возьмём данные за 2019 год. Это достаточно длинный период времени, к тому же не искажённый пандемией.")
st.markdown("Всего в таблице `flights` 118 714 936 записей, а за 2019 год - 32 088 373. Это довольно много, heroku такое не вывезет. Поэтому сделаем несколько файлов с нужными агрегациями.")
st.markdown("Считаем количество перелетов по каждой (упорядоченной) паре аэропортов")
st.markdown("Исполняем SQL-запрос:")
st.markdown("""`select
departure, destination, count(icao24)
from flights
where
firstseen >= "2019"
and lastseen < "2020"
group by departure, destination;`
""")
st.markdown("Сохраняем в файл flights_by_airport_pairs.txt")
st.markdown("Теперь напишем функцию, которая будет парсить этот файл")
with st.echo():
    def parse_hadoop_output(filename):
        f = open(filename, "r")
        for i, line in enumerate(f):
            if i == 1:
                header = parse_line(line)
                ncolumns = len(header)
                columns_data = dict.fromkeys(header, 0)
            elif len(line) == 1:
                continue
            elif line[1] != "-":
                content = parse_line(line)
                if content != header:
                    if columns_data[header[0]] == 0:
                        for ncol in range(ncolumns):
                            columns_data[header[ncol]] = [content[ncol],]
                        continue
                    for ncol in range(ncolumns):
                        columns_data[header[ncol]].append(content[ncol])
        f.close()
        df = pd.DataFrame(data=columns_data).rename(columns={'count(icao24)': 'count'})
        return df

    def parse_line(line):
        content = [x.strip(" ") for x in line.strip("|").split("|")][:-1]
        return content

    df_fl = parse_hadoop_output("flights_by_airports.txt")
    df_fl[:5]
    st.write(df_fl["departure"].nunique())

st.markdown("## 4. Получение данных по аэропортам")
st.markdown("Чтобы сделать красивые визуализации, нужно понять, где находится каждый аэропорт.")
st.markdown("Это отличная возможность использовать **Selenium**!")
st.markdown("Мой комп умрёт парсить 14 000 аэропортов, поэтому тберём только крупные. Скажем, с как минимум 500 рейсами в год.")
st.markdown("Таких всего:")

with st.echo(code_location="below"):
    df = df_fl.copy()
    df = df[pd.notna(df.departure)]
    df = df[pd.notna(df.destination)]
    df_dep = df.departure.value_counts().copy().reset_index()
    df_des = df.destination.value_counts().copy()
    df = df_dep.join(df_des, how='outer', on="index")
    df.fillna(0, inplace=True)
    df['total'] = df.departure + df.destination
    df = df[df.total >= 500]
    df.rename(columns={'index': 'airport'}, inplace=True)
    df.to_csv("major_airports.csv", sep=';')
    st.write(df.shape[0])

st.markdown("Чтобы тебе не пришлось ждать вечность, пока все аэропорты пропарсятся на Heroku, я делаю парсинг у себя, а результаты сохраняю в `airport_data.csv`")
st.markdown("Парсить будем с сайта https://www.avcodes.co.uk/aptcodesearch.asp")
st.markdown("Вот код парсера:")
st.code("""
from selenium.webdriver import Chrome
from selenium.webdriver.common.by import By
import pandas as pd

driver=Chrome("/Users/afonin/PycharmProjects/temp-avia-data/chromedriver")

def get_info(icao):
    driver.get("https://www.avcodes.co.uk/aptcodesearch.asp")

    driver.find_element_by_name('icaoapt').send_keys(icao)
    driver.find_element_by_name('B1').click()

    lon = driver.find_element(By.XPATH, "//*[contains(text(), 'Longitude:')]").get_attribute("innerHTML").splitlines()
    lat = driver.find_element(By.XPATH, "//*[contains(text(), 'Latitude:')]").get_attribute("innerHTML").splitlines()
    name = driver.find_element_by_class_name('tablebg').get_attribute("innerHTML").splitlines()
    country = driver.find_element(By.XPATH, "//*[contains(text(), 'Country:')]").get_attribute("innerHTML").splitlines()

    lon = lon[0][14:]
    if lon[-1] == "W":
        lon = -(int(lon[:3]) + int(lon[4:6])/60 + int(lon[7:9])/3600)
    else:
        lon = int(lon[:3]) + int(lon[4:6]) / 60 + int(lon[7:9]) / 3600

    lat = lat[0][13:]
    if lat[-1] == "S":
        lat = -((int(lat[:3]) + int(lat[4:6])/60 + int(lat[7:9])/3600))
    else:
        lat = int(lat[:3]) + int(lat[4:6]) / 60 + int(lat[7:9]) / 3600

    city = name[0].strip(' ').split(' / ')[0]
    name = name[0][name[0].find('/')+2:].replace(' ', '')

    country = country[0].split('<br>')[1]

    return [lat, lon, city, name, country]

df = pd.read_csv("major_airports.csv", sep=';')
airports = df.airport.to_list()
df = pd.DataFrame(columns=['airport', 'lat', 'lon'])
cnt = 0

for airport in airports[1:]:
    try:
        lat, lon, city, name, country = get_info(airport)
        if (country == "United States of America") and (lon > 0):  # фиксю баг сайта
            lon = -lon
        df = df.append({'airport': airport, 'lat': lat, 'lon': lon, 'name': name, 'city': city, 'country': country},
                       ignore_index=True)
    except:
        cnt += 1
        print(f"{cnt}. {airport}")

df.to_csv("major_airports_data.csv", sep=';')
""")

st.markdown("Парсер нашёл данные только по 401 аэропорту. Остальные аэропорты, видимо, слишком маленькие, и их нет в базе на том сайте.")

with st.echo():
    df_ap = pd.read_csv("major_airports_data.csv", sep=';')
    st.write(df_ap[:5])

st.markdown("## 5. Визуализации")

st.markdown("### Граф связности аэропортов")
st.markdown("""Интересно посмотреть, какие есть ключевые транспортные хабы - аэропорты, которые связаны с многими другими.
Скорее всего, это будут какие-нибудь Heathrow (Лондон), O'Hare (Чикаго), Аэропорт Шанхая, JFK (Нью-Йорк).
В этих аэропортах много трансконтинентальных рейсов, пассажиры которых затем пересаживаются и летят в более мелкие аэропорты.
Проверим интуицию :)""")
st.markdown("Выкинем все перелёты, где одна из двух точек неизвестна, а также перелёты, кончившиеся там же, где и начались")
st.markdown("""Выкинем также все связки между аэропортами, где меньше 1460 перелётов в год в каждую сторону (по 4 в неделю).
Можно было бы поставить отсечку поменьше, но тогда граф будет слишком нечитаемый - точки будут неразличимы друг от друга.
С другой библиотекой (не `networkx`, а `Gephi`) было бы лучше.""")
st.markdown("Подружить библиотеку для отрисовки `networkx` (`pygraphviz`) с Heroku довольно непросто. Поэтому я прилагаю код, который запустил локально, и скриншот графа.")

st.image("networkx-graph.png")

st.code("""
    df_fl_graph = df_fl.copy()
    df_fl_graph = df_fl_graph[df_fl_graph['destination'] != "NULL"]
    df_fl_graph = df_fl_graph[df_fl_graph['departure'] != "NULL"]
    df_fl_graph = df_fl_graph[df_fl_graph.departure != df_fl_graph.destination]
    df_fl_graph['a1'] = df_fl_graph[['destination', 'departure']].min(axis=1)
    df_fl_graph['a2'] = df_fl_graph[['destination', 'departure']].max(axis=1)
    df_fl_graph.drop(['departure', 'destination'], inplace=True, axis=1)
    df_fl_graph['count'] = df_fl_graph['count'].astype(int)
    df_fl_graph = df_fl_graph[df_fl_graph['count'] >= 365*4]
    df_fl_graph = df_fl_graph.groupby(['a1', 'a2']).sum().reset_index()
    G = nx.from_pandas_edgelist(df_fl_graph, source='a1', target='a2', edge_attr="count")
    degrees = sorted(G.degree, key=lambda x: x[1], reverse=True)
    top_nodes = [x[0] for x in degrees[:5]]
    other_nodes = [x[0] for x in degrees[5:]]
    plt.plot(1)
    nx.draw(G, node_size=[5, ]*len(other_nodes) + [300, ]*5,
            pos=nx.nx_agraph.graphviz_layout(G, prog="fdp"),
            node_color=['blue', ]*len(other_nodes) + ['orange', ]*5,
            nodelist=other_nodes + top_nodes,
            edge_color='grey',
            labels=dict(zip(top_nodes, top_nodes)),
            alpha=0.85,
            font_size=4,
            font_family='helvetica',
            width=0.1)
    fig = plt.plot()
    st.set_option('deprecation.showPyplotGlobalUse', False)
    st.pyplot(fig)
""")

st.markdown("Видно, что 5 узловых аэропортов - это EGLL (Heathrow, London), EDDF (Franfkfurt), KORD (O'Hare, Chicago), KDFW (Dallas-Fort Worth), KLAX (Los Angeles). Интуиция плюс-минус не подвела :)")

st.markdown("### Самые крупные аэропорты на карте-1")
st.markdown("Давай посмотрим, где находятся 10 крупнейших аэропортов мира (по количеству рейсов)")

with st.echo(code_location="below"):
    st.set_option('deprecation.showPyplotGlobalUse', False)
    df1 = pd.read_csv("major_airports.csv", sep=';').drop("Unnamed: 0", axis=1)
    df2 = pd.read_csv("major_airports_data.csv", sep=';').drop("Unnamed: 0", axis=1)
    df = df2.join(df1.set_index("airport"), on='airport').drop(columns=['departure', 'destination'])
    df['total'] = df.total.astype('int32')
    df_slice = df.sort_values(by=['total'], ascending=False)[:10]
    gdf = gpd.GeoDataFrame(df_slice, geometry=gpd.points_from_xy(df_slice.lon, df_slice.lat))
    world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
    ax = world.plot(facecolor='w', edgecolor='k')
    plt.plot(2)
    gdf.plot(markersize=np.square(gdf['total'])/1000000, ax=ax)
    fig2 = plt.plot()
    st.pyplot(fig2)

st.markdown("Все в Америке, кто бы мог подумать...")

st.markdown("### Самые крупные аэропорты на карте-2")
st.markdown("У geopandas не самый богатый встроенный функционал. По сути он работает на `matplotlib` без каких-то модификаций, поэтому красоты и интерактивности от него ждать не стоит.")
st.markdown("Хорошо, что есть библиотеки получше! Например, `plotly`")
st.markdown("Оставим только 100 крупнейших аэропортов")
st.markdown("Картинку лучше открыть целиком (для этого есть кнопка в правом верхнем углу).")

with st.echo(code_location="below"):
    df = df.sort_values(by=['total'], ascending=False)[:100]
    fig3 = px.scatter_geo(df, lat='lat', lon='lon', size=df['total']/1000, color='country',
                          hover_data=['country', 'city', 'name', 'total'],
                          hover_name='name', opacity=0.7, scope='world', title='Крупнейшие аэропорты',
                          basemap_visible=True)
    st.write(fig3)

st.markdown("Картинка получилась посимпатичнее, круто!")

st.markdown("## 6. Получение данных по самолётам")
st.markdown("Прежде чем уйти в R, чтобы поработать там с ML, нужно подготовить данные.")
st.markdown("Будем пользоваться API сайта aviationstack.com.")
st.markdown("Там всего 500 реквестов, поэтому я работаю с API локально, а код просто прикрепляю")
st.code("""
import json as js
import pandas as pd
import requests as rq

def execute_request(API_key, entrypoint, offset):
    url = entrypoint + f"&offset={offset}"
    r = rq.get(url)
    j = r.json()
    data = j['data']
    df = pd.DataFrame(data)
    return df

API_key = "YOUR_API_KEY"  # если хочешь, можешь получить свой ключ на сайте
entrypoint = f"http://api.aviationstack.com/v1/airplanes?access_key={API_key}"

for i in range(190):
    try:
        df = execute_request(API_key=API_key, entrypoint=entrypoint, offset=i*100)
    except:
        print(i)
    if i == 0:
        df_fetched = df.copy()
    else:
        df_fetched = df_fetched.append(df, ignore_index=True)

df_fetched.to_csv("aircraft_data.csv", sep=';')
""")

st.markdown("""В датасете, который мы собрали через API, есть icao24-коды самолётов, дата изготовления и модель.
Нужно собрать ещё один численный параметр - например, среднюю длительность полёта. Заберём его из OpenSky Network c
помощью SQL-запроса:""")
st.markdown("`SELECT icao24, avg(duration) from flights where firstseen >= '2019' and lastseen < '2020' group by icao24`")
st.markdown("Осталось только связать два последних датасета между собой.")
st.markdown("Поскольку это нужно для R, а  R лежит в другом приложении, я делаю это локально, а код просто копирую сюда, чтобы был.")

st.code("""
df_fld = parse_hadoop_output("avg_flight_duration.txt")
df_fld['icao24'] = df_fld['icao24'].apply(lambda x: x.upper())
df_ac = pd.read_csv("aircraft_data.csv", sep=';')
df_ac = df_ac[['icao_code_hex', 'iata_type', 'first_flight_date']]
df = df_ac.join(df_fld.set_index('icao24'), on='icao_code_hex', how='inner')
df.rename(columns={'avg(duration)': 'avg_duration'}, inplace=True)
df.to_csv("knn_data.csv", sep=';')
""")

st.markdown("Всего часть с питоном занимает около 130 строчек. А ещё есть часть в R, так что 120 строк я точно набираю :)")

ext: py

reason: import numpy


text:
#   Copyright (c) 2021, Xilinx, Inc.
#   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 of the copyright holder 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 HOLDER 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.
#

################################################################################
#
# Authors: Timoteo Garcia Bertoa <timoteog@xilinx.com>
#
# \file tmrmapping.py
#
# This library contains function definitios for TMR implementation.
#
################################################################################

import numpy as np

def isvalid(pe, ofm_ch, redf, num_red):

	"""
    =======================================================================================

    This function checks if, given a number of PEs and number of redundancies for an OFM,
	it is possible to map the complete channel stack so that at least each repeated channel
	is computed by separate PEs.

	Errors are flagged if the number of PEs suggested is not valid.

	Parameters:

	 - pe      : number of PEs
	 - ofm_ch  : number of OFM channels, not counting the triplications
	 - redf    : redundancy factor, (3 for triplication)
	 - num_red : number of triplicated channels

    =======================================================================================
    """

	if pe < 3:
		print("Error: Number of PEs must be at least 3, to be able to map 1 PE per triplicated channel")
		exit(1)
	else:

		# Number of output channels including triplications
		ofm_ch_red = ofm_ch + (num_red * (redf - 1))

		if pe > ofm_ch_red:
			print("Error: Number of PEs can be reduced to {}".format(ofm_ch_red))
			exit(1)
		elif pe == ofm_ch_red:
			print("Number of PE = {} is valid!".format(pe))
		elif pe < ofm_ch_red:
			if ofm_ch_red % pe == 0:
				fold = ofm_ch_red / pe
				print("Number of PE = {} is valid, neuron fold = {}".format(pe, fold))
			else:
				print("Error: Number of PEs = {} is not valid. Number of total channels should be divisible by number of PE.".format(pe))
				exit(1)

def map(param_in, ofm_ch, pe, redf, num_red, ch_list = []):

	"""
    =======================================================================================

    This function maps the triplications and returns the new set of weights, providing also
	a channel mask value and a vector of indexes.

	The user can provide weights, number of redundancies and redundancy factor, and the
	function will return a new set of weights with the specified channels repeated as many times
	as selected with the redundancy factor.

	The user should specify which channels are desired to be repeated.

	Parameters:

	 - param_in     : matrix of weights, NxM, where N = OFM channels (not inc. trip.), M = tile
	 - ofm_ch       : number of OFM channels, not counting the triplications
	 - pe	        : number of PEs used to compute convolution
	 - redf         : redundancy factor, (3 for triplication)
	 - num_red      : number of triplicated channels
	 - ch_list      : list of channel indexes desired to be repeated

	Returns:

	 - param_out    : matrix of weights, NxM, where N = OFM channels (inc. trip.), M = tile
	 - channel_mask : value to mask repeated (1) and not repeated channels (0)
	 - red_ch_index : vector with the indexes of the repeated channels

    =======================================================================================
    """

    # Flag if the number of triplications don't match
	if len(ch_list) != num_red:
		print("Error: You have specified {} triplications, but {} channel/s to triplicate".format(num_red, len(ch_list)))
		exit(1)

    # Flag if user specifies a wrong channel number
	for x in ch_list:
		if x >= ofm_ch:
			print("Error: You have specified an invalid channel, it should be betwen 0 and {}!".format(ofm_ch - 1))
			exit(1)

	# Number of output channels including triplications
	ofm_ch_red = ofm_ch + (num_red * (redf - 1))

	# Generate channel_mask
	channel_mask = ''
	for k in range(ofm_ch):
		added = 0
		for ch in ch_list:
			if k == ch:
				channel_mask = "111" + channel_mask
				added = 1
		if added == 0:
			channel_mask = "0" + channel_mask
	channel_mask = "0b" + channel_mask
	channel_mask = int(channel_mask, 2)

	red_ch_index = np.zeros((num_red, 1))
	# Generate red_ch_index, and order rows to show LSB index first
	aux = 0
	ch_count = 0
	for ch in range(ofm_ch_red):
		if ch in ch_list:
			red_ch_index[aux] = ch_count
			aux += 1
			ch_count += 2
		ch_count += 1
	red_ch_index = red_ch_index[np.lexsort(np.fliplr(red_ch_index).T)]

	ch_repeats = []
	# Create list with number of repeats per channel, and triplicate weights
	for ch in range(ofm_ch):
		added = 0
		for x in ch_list:
			if ch == x:
				ch_repeats.append(redf)
				added = 1
		if added == 0:
			ch_repeats.append(1)

	# Triplicate channels
	param_out = np.repeat(param_in, ch_repeats, axis = 0)

	# Consider folding case, split in chunks and concatenate the rows a neuron folding number of times
	param_out = np.concatenate(np.array_split(param_out, ofm_ch_red/pe), axis=1)

	return param_out, channel_mask, red_ch_index

max_stars_repo_name: timoteogb/finn-hlslib

max_issues_repo_name: timoteogb/finn-hlslib

max_forks_repo_name: timoteogb/finn-hlslib


text:
import os
import glob
import inspect
import warnings
import datetime
import functools
import contextlib
import collections

import numpy as np
import pandas as pd
import tqdm

from . import getdata
from . import getstatic
from .. import config

def _filter_df(df: pd.DataFrame, query: str) -> pd.DataFrame:
    if '@' in query:
        # In this case, we appear to be using
        # external variables, so we need to
        # pull in the local variables from the caller

        frame = inspect.currentframe()
        caller_locals = frame.f_back.f_locals
        with contextlib.suppress(TypeError):
            # We hop back two frames, so the caller and the caller's caller.
            caller_locals.update(frame.f_back.f_back.f_locals)

    else:
        # Otherwise, leave it empty.
        caller_locals = None

    try:
        return df.query(query, local_dict=caller_locals)
    except TypeError:
        # This may be triggered by an operation unsupported by Numexpr.
        # Try to fall back to Python backend.
        return df.query(query, local_dict=caller_locals, engine='python')
    except pd.core.computation.ops.UndefinedVariableError:
        msg = ('One of the keys does not seem to appear in some of the data. '
               'This may be due to wrong spelling, or simply missing from one of the hours.')
        warnings.warn(RuntimeWarning(msg))
        return pd.DataFrame()

def _clear_duplicates(df: pd.DataFrame, feed: str) -> None:
    if feed == 'TripUpdates' and not df.empty:
        df.sort_values(by='timestamp', ascending=True)

        # Not all feeds provide exactly the same fields, so this filters for it:
        keys = list({'trip_id', 'direction_id', 'stop_sequence', 'stop_id'}.intersection(df.keys()))
        df.drop_duplicates(subset=keys, keep='last', inplace=True)

def get_data_range(feed: str, company: str, start_date: str, start_hour: int = 0, end_date: str = None,
                   end_hour: int = 23, query: (str, None) = None, merge_static: bool = False,
                   clear_duplicates: bool = True) -> pd.DataFrame:
    """Get all data from start_date:start_hour to end_date:end_hour, inclusive, downloading if necessary.
    Date should be given in the format YYYY-MM-DD or YYYY_MM_DD.

    If end_date is not specified, it is the same as start_date.

    Use the parameter `query` to filter the data as it is being read if you want to minimise the memory footprint.

    merge_static includes some information from the static feed.
    """
    if feed not in ('VehiclePositions', 'TripUpdates', 'ServiceAlerts'):
        raise ValueError(f'Feed {feed} not recognised.')
    if company not in ('dt', 'klt', 'otraf', 'skane', 'sl', 'ul', 'varm', 'xt'):
        warnings.warn(RuntimeWarning(f'Company {company} is not recognised. Maybe the API will fail?'))

    frames = []
    if end_date is None:
        end_date = start_date

    start = start_date.replace('_', '-') + f' {start_hour}:00'
    end = end_date.replace('_', '-') + f' {end_hour}:00'

    for date_hour in tqdm.tqdm(pd.date_range(start, end, freq='h'), desc='Loading data'):
        date, time_code = str(date_hour).split()
        hour = int(time_code.split(':')[0])
        full_path = getdata._get_data_path(company, feed=feed, hour=hour, date=date)
        if not os.path.exists(full_path):
            # Download if it doesn't exist'
            try:
                getdata.get_data(date, hour, feed, company)
            except ValueError:
                warnings.warn(RuntimeWarning(f'The API did not return data for {date_hour}'))
                continue

        df = pd.read_feather(full_path)
        if clear_duplicates:
            _clear_duplicates(df, feed)

        if feed == 'VehiclePositions':
            if 'trip_id' in df:
                # For vehicle positions, skip data that does not correspond to routes.
                df = df.query('trip_id.notna()', engine='python')
            else:
                df = pd.DataFrame()

        if merge_static and not df.empty:
            this_static = load_static_data(company, date, remove_unused_stations=True)
            if this_static is None:
                continue

            if feed == 'TripUpdates':
                # Overwrite some parameters already present in the static data
                df.drop(columns=['route_id', 'direction_id'], errors='ignore', inplace=True)

                # Rename some common columns
                df.rename(
                    columns={'arrival_time': 'observed_arrival_time', 'departure_time': 'observed_departure_time'},
                    inplace=True)

                this_static.stop_times.rename(columns={'arrival_time': 'scheduled_arrival_time',
                                                       'departure_time': 'scheduled_departure_time'}, inplace=True)

                df = df.merge(this_static.stop_times, how='left',
                              on=('trip_id', 'stop_id', 'stop_sequence'),
                              validate='m:1')
                df = df.merge(this_static.trips.reset_index(level=['route_id', 'direction_id']), how='left',
                              on='trip_id', validate='m:1')

            elif feed == 'VehiclePositions':
                df.drop(columns=['route_id', 'direction_id'], errors='ignore', inplace=True)
                df = df.merge(this_static.trips.reset_index(level=['route_id', 'direction_id']), how='left',
                              on='trip_id', validate='m:1')

            elif feed == 'ServiceAlerts':
                warnings.warn(RuntimeWarning('ServiceAlerts cannot be merged with static data.'))

        if query is not None:
            df = _filter_df(df, query)

        # Drop the index, since it will be regenerated when concatenated
        df.drop(columns='index', errors='ignore', inplace=True)
        frames.append(df)

    df_merged = pd.concat(frames)
    if 'timestamp' in df_merged.keys():
        df_merged['datetime'] = pd.to_datetime(df_merged.timestamp, unit='s')
    return df_merged

static_data = collections.namedtuple('StaticData', ['stop_times', 'stops', 'trips', 'shapes', 'routes'])

def _remove_unused_stations(stops_data: pd.DataFrame, stops_times: pd.DataFrame) -> (pd.DataFrame, pd.DataFrame):
    parent_stations = stops_data.parent_station.dropna().values
    used_stations = set(stops_times.stop_id).union(parent_stations)
    stops_times = stops_times.query('stop_id in @used_stations')
    stops_data = stops_data.query('stop_id in @used_stations')
    return stops_data, stops_times

@functools.lru_cache(1)
def load_static_data(company: str, date: str, remove_unused_stations: bool = False) -> (static_data, None):
    """Load static data from the cache, downloading it if necessary.

    Date should be given in the format YYYY-MM-DD or YYYY_MM_DD.
    """
    output_folder = getstatic._get_static_data_path(company, date)
    if not os.path.isdir(output_folder):
        try:
            getstatic.get_static_data(company=company, date=date)
        except ValueError:
            warnings.warn(RuntimeWarning(f'Missing static data for {company} at {date}'))
            return

    if not os.path.exists(os.path.join(output_folder, 'stop_times.txt')):
        warnings.warn(RuntimeWarning(f'Missing static data for {company} at {date}'))
        return

    def parse_date(time_code: str):
        time_code = list(map(int, time_code.split(':')))
        if time_code[0] < 24:
            return datetime.datetime(*map(int, date.replace('_', '-').split('-')), *time_code)
        else:
            # If it is over 24 h, roll over the next day.
            time_code[0] -= 24
            dt = datetime.datetime(*map(int, date.replace('_', '-').split('-')), *time_code)
            return dt + datetime.timedelta(days=1)

    stops_times = pd.read_csv(os.path.join(output_folder, 'stop_times.txt'), dtype={'trip_id': 'str', 'stop_id': 'str'},
                              parse_dates=['arrival_time', 'departure_time'], date_parser=parse_date)
    stops_data = pd.read_csv(os.path.join(output_folder, 'stops.txt'),
                             dtype={'stop_id': 'str', 'parent_station': 'str'})
    trips_data = pd.read_csv(os.path.join(output_folder, 'trips.txt'), dtype={'trip_id': 'str', 'route_id': 'str'})
    shapes_data = pd.read_csv(os.path.join(output_folder, 'shapes.txt'))
    routes = pd.read_csv(os.path.join(output_folder, 'routes.txt'), dtype={'route_id': 'str', 'agency_id': 'str',
                                                                           'route_short_name': 'str',
                                                                           'route_long_name': 'str'})

    if remove_unused_stations:
        stops_data, stops_times = _remove_unused_stations(stops_data, stops_times)

    stop_times = pd.merge(stops_times, stops_data, on='stop_id', how='left', validate='m:1')
    trips_data = pd.merge(trips_data, routes, on='route_id', how='left', validate='m:1')

    # Set indexes for faster querying and merging
    trips_data.set_index(['trip_id', 'route_id', 'direction_id'], inplace=True, drop=True, verify_integrity=True)
    stop_times.set_index(['trip_id', 'stop_id', 'stop_sequence'], inplace=True, drop=True, verify_integrity=True)
    shapes_data.set_index(['shape_id', 'shape_pt_sequence'], inplace=True, drop=True, verify_integrity=True)
    stops_data.set_index('stop_id', inplace=True, drop=True, verify_integrity=True)

    data = static_data(stop_times=stop_times, stops=stops_data, trips=trips_data, shapes=shapes_data, routes=routes)
    return data

def clean_cache() -> None:
    """Apply cleaning filters to all files in the cache, to fix issues with arrays downloaded with older versions."""

    for f in tqdm.tqdm(glob.glob(os.path.join(config.CACHE_DIR, '*feather'))):
        df = pd.read_feather(f)
        if 'ServiceAlerts' in f:
            df = getdata.unpack_jsons(df)
        getdata.sanitise_array(df)

        if df.empty:  # Feather does not support a DF without columns, so add a dummy one
            df['_'] = np.zeros(len(df), dtype=np.bool_)
        df.reset_index(drop=True).to_feather(f, compression='zstd', compression_level=9)

max_stars_repo_name: KoDa-project/pykoda

max_issues_repo_name: KoDa-project/pykoda

max_forks_repo_name: KoDa-project/pykoda


text:
import numpy as np
import pandas as pd

from asset_model import geometric_brownian_motion

#TODO
# class CcpiStrategy(InvestmentStrategy):
#
#     def __init__(self, drawdown=None, multiplier=3):
#         self.drawdown = drawdown
#         self.multiplier = multiplier
#
#     def update_portfolio_weighs(self, current_weights, account_value, returns):
#         risky_weight = current_weights[0]
#         safe_weight = current_weights[1]
#         #TODO
#         # if self.drawdown is not None:
#         #     peak = np.maximum(peak, account_value)
#         #     floor_value = peak * (1 - drawdown)
#         cushion = (account_value - floor_value) / account_value
#         risky_weight = self.multiplier * cushion
#         risky_weight = np.minimum(risky_weight, 1)
#         risky_weight = np.maximum(risky_weight, 0)
#         safe_weight = 1 - risky_weight
#         return np.array(risky_weight, safe_weight)

def backtest_cppi(risky_returns,
                  safe_returns=None,
                  risk_free_rate=0.03,
                  multiplier=3,
                  cushion_ratio=0.8,
                  drawdown=None,
                  start_value=1000):
    """
    Run a backtest of the CPPI strategy, given a set of returns for the risky asset
    Returns a dictionary containing: Asset Value History, Risk Budget History, Risky Weight History
    :param risky_returns: history of risky returns
    :param safe_returns: history of safe returns. If None, default to risk free rate
    :param risk_free_rate: Rate of return of the risk free asset
    :param multiplier: multiplier to allocate to risky asset = multiplier*(1-cushion)*wealth
    :param cushion_ratio: ratio of the wealth to protect
    :param drawdown: max drawdown allowed (as ratio)
    :param start_value: Initial monetary value of the account
    """

    # Ensure returns are a data frame
    if isinstance(risky_returns, pd.Series):
        risky_returns = pd.DataFrame(risky_returns, columns=["R"])

    # If no safe asset is specified, default to the risk free rate
    if safe_returns is None:
        safe_returns = pd.DataFrame().reindex_like(risky_returns)
        safe_returns.values[:] = risk_free_rate / 12  # fast way to set all values to a number

    # set up the CPPI parameters
    dates = risky_returns.index
    steps = len(dates)
    account_value = start_value
    floor_value = start_value * cushion_ratio
    peak = account_value

    # set up some DataFrames for saving intermediate values
    account_history = pd.DataFrame().reindex_like(risky_returns)
    risky_w_history = pd.DataFrame().reindex_like(risky_returns)
    cushion_history = pd.DataFrame().reindex_like(risky_returns)
    floorval_history = pd.DataFrame().reindex_like(risky_returns)
    peak_history = pd.DataFrame().reindex_like(risky_returns)

    for step in range(steps):
        # If a drawdown is specified, re-calibrate the floor value
        if drawdown is not None:
            peak = np.maximum(peak, account_value)
            floor_value = peak * (1 - drawdown)
        cushion = (account_value - floor_value) / account_value
        risky_weight = multiplier * cushion
        risky_weight = np.minimum(risky_weight, 1)
        risky_weight = np.maximum(risky_weight, 0)
        safe_weight = 1 - risky_weight
        risky_alloc = account_value * risky_weight
        safe_alloc = account_value * safe_weight
        # recompute the new account value at the end of this step
        account_value = risky_alloc * (1 + risky_returns.iloc[step]) + safe_alloc * (1 + safe_returns.iloc[step])

        # save the histories for analysis and plotting
        cushion_history.iloc[step] = cushion
        risky_w_history.iloc[step] = risky_weight
        account_history.iloc[step] = account_value
        floorval_history.iloc[step] = floor_value
        peak_history.iloc[step] = peak
    risky_wealth = start_value * (1 + risky_returns).cumprod()
    return {
        "wealth": account_history,
        "risky_wealth": risky_wealth,
        "risk_budget": cushion_history,
        "risky_allocation": risky_w_history,
        "peak": peak_history,
        "floor": floorval_history
    }

def cppi_monte_carlo(scenarios=50, mu=0.07, sigma=0.15,
                     multiplier=3, cushion_ratio=0.0, risk_free_rate=0.03, start_value=1000):
    risky_returns = geometric_brownian_motion(scenarios=scenarios, mu=mu, sigma=sigma, prices=False)
    return backtest_cppi(risky_returns=risky_returns, risk_free_rate=risk_free_rate,
                         multiplier=multiplier, start_value=start_value, cushion_ratio=cushion_ratio)

max_stars_repo_name: algorythmist/finance-toolbox

max_issues_repo_name: algorythmist/finance-toolbox

max_forks_repo_name: algorythmist/finance-toolbox


text:
from keras import backend as K
from keras.layers import Dense, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras.layers import Activation, Dropout, BatchNormalization
from keras.models import Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import Adam
from keras.preprocessing.image import ImageDataGenerator
from keras.preprocessing.image import img_to_array
from sklearn.preprocessing import LabelBinarizer
from sklearn.model_selection import train_test_split
import os
import random
import numpy as np
import cv2
import glob
import pickle
import matplotlib.pyplot as plt

# Load Data from folder
current_path = os.getcwd()
All_files = glob.glob(current_path + '/Reduced_Data_further/**/*')
print("Number of total files: " + str(len(All_files)))

######--HYPERPARAMETERS--#######
Epochs = 1000
Learning_Rate = 1e-5
Batch_Size = 5
Image_Dimensions = (120, 84, 3)
model_name = '2-32-2-64-2-128--5'
act_type = "relu"
filter_size = (3, 3)
pad_type = "same"
dropout_rate = 0.5
################################

# data and label intitialization
os.makedirs(model_name)
data = []
labels = []

random.shuffle(All_files)

for image_name in All_files:
    image_data = cv2.imread(image_name)
    image_data = cv2.resize(image_data, (Image_Dimensions[0], Image_Dimensions[1]))
    image_data = img_to_array(image_data)
    data.append(image_data)
    label = image_name.split(os.path.sep)
    labels.append(label[-2])

data = np.array(data, dtype="float") /255.0
labels = np.array(labels)

lb = LabelBinarizer()
labels = lb.fit_transform(labels)

# Data split into training and testing set
(x_train, x_test, y_train, y_test) = train_test_split(data, labels, test_size=0.2)

# Data Augmentation
data_augmentation = ImageDataGenerator(width_shift_range=0.15, height_shift_range=0.15)

# Model inputs
width = Image_Dimensions[0]
height = Image_Dimensions[1]
depth = Image_Dimensions[2]
input_parameters = (height, width, depth)
channel_dimensions = -1
classes = len(lb.classes_)

# Channels depend on the OS used
if K.image_data_format() == "channels_first":
    input_parameters = (depth, height, width)
    channel_dimensions = 1

# Defining model architecture
model = Sequential()
model.add(Conv2D(32, filter_size, padding=pad_type, input_shape=input_parameters))
model.add(Activation(act_type))
model.add(BatchNormalization(axis=channel_dimensions))
model.add(Dropout(dropout_rate))

model.add(Conv2D(32, filter_size, padding=pad_type))
model.add(Activation(act_type))
model.add(BatchNormalization(axis=channel_dimensions))
model.add(MaxPooling2D(pool_size=(3, 3)))
model.add(Dropout(dropout_rate))

model.add(Conv2D(64, filter_size, padding=pad_type))
model.add(Activation(act_type))
model.add(BatchNormalization(axis=channel_dimensions))
model.add(Dropout(dropout_rate))

model.add(Conv2D(64, filter_size, padding=pad_type))
model.add(Activation(act_type))
model.add(BatchNormalization(axis=channel_dimensions))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(dropout_rate))

model.add(Conv2D(128, filter_size, padding=pad_type))
model.add(Activation(act_type))
model.add(BatchNormalization(axis=channel_dimensions))
model.add(Dropout(dropout_rate))

model.add(Conv2D(128, filter_size, padding=pad_type))
model.add(Activation(act_type))
model.add(BatchNormalization(axis=channel_dimensions))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(dropout_rate))

model.add(Flatten())
model.add(Dense(1024, use_bias=False))
model.add(BatchNormalization())
model.add(Activation(act_type))
model.add(Dropout(dropout_rate+0.3))

model.add(Dense(classes))
model.add(Activation("softmax"))

# Model compilation
opt = Adam(lr=Learning_Rate, decay=Learning_Rate/Epochs)
model.compile(loss="categorical_crossentropy", optimizer=opt, metrics=["accuracy"])
filepath = model_name + "/weights-improvement-{epoch:02d}-{val_acc:.2f}.hdf5"
checkpoint = ModelCheckpoint(filepath, monitor='val_acc', verbose=1, save_best_only=True, mode='max')
callbacks_list = [checkpoint]

# Train Network
model_history = model.fit_generator(data_augmentation.flow(x_train, y_train, batch_size=Batch_Size), validation_data=(x_test, y_test),
   steps_per_epoch=(len(x_train)//Batch_Size), epochs=Epochs, verbose=1, callbacks=callbacks_list)
data_augmentation.flow(x_train, y_train, batch_size=Batch_Size)

# Save Model
model.save(model_name+ '/' + model_name+'.model')

# Save Labels
pickle_file = open('labels.pickle', "wb")
pickle_file.write(pickle.dumps(lb))
pickle_file.close()

# Plot loss and accuracy of model
plt.figure()
plt.plot(np.arange(0, Epochs), model_history.history["acc"], label="Training Accuracy")
plt.plot(np.arange(0, Epochs), model_history.history["val_acc"], label="Validation Accuracy")
plt.title("Accuracy vs Validation Accuracy")
plt.xlabel("Epoch Number")
plt.ylabel("Accuracy")
plt.legend(loc="upper left")
plt.savefig(model_name+ '/' + model_name + '-acc.png')

plt.figure()
plt.plot(np.arange(0, Epochs), model_history.history["loss"], label="Training Loss")
plt.plot(np.arange(0, Epochs), model_history.history["val_loss"], label="Validation Loss")
plt.title("Loss vs Validation Loss")
plt.xlabel("Epoch Number")
plt.ylabel("Loss")
plt.legend(loc="upper left")
plt.savefig(model_name+ '/' + model_name + '-loss.png')

# Creating arrays of the data to plot manually
accuracy_history = np.array(model_history.history["acc"])
val_accuracy_history = np.array(model_history.history["val_acc"])
loss_history = np.array(model_history.history["loss"])
val_loss_history = np.array(model_history.history["val_loss"])

# Saving arrays to text files
np.savetxt(model_name+ '/' + model_name + "-accuracy.txt", accuracy_history, delimiter=",")
np.savetxt( model_name+ '/' + model_name + "-val_accuracy.txt", val_accuracy_history, delimiter=",")
np.savetxt(model_name+ '/' + model_name + "-loss.txt", loss_history, delimiter=",")
np.savetxt(model_name+ '/' + model_name + "-val_loss.txt", val_loss_history, delimiter=",")

max_stars_repo_name: BAD-Classifier/Machine-Learning

max_issues_repo_name: BAD-Classifier/Machine-Learning

max_forks_repo_name: BAD-Classifier/Machine-Learning


text:
import shutil
from pathlib import Path
from typing import List, Tuple
from urllib.request import urlopen

import numpy as np

import cv2

def download_remote_data_file(local_dir_name: str,
    remote_dir_url: str, file_name: str) -> Tuple[str, bool]:
    """Downloads data from url if it's not saved locally yet.

    Args:
        local_dir_name (str): The local folder where we'll save the
        downloaded data.
        remote_dir_url (str): The url prefix of the data file we want to download.
        file_name (str): The name of the data file we want to download.
    Returns:
        The path to the local file, and a bool indicating whether the file
        was downloaded or not.
    """
    # Create a data directory if it doesn't exist.
    local_dir_path = Path('.', local_dir_name)
    local_dir_path.mkdir(exist_ok=True)

    # Download the data file if it doesn't exist.
    remote_file_url = remote_dir_url + file_name
    local_file_path = Path(local_dir_path, file_name)
    downloaded = False
    if not local_file_path.exists():
        print(f'Downloading data file {remote_file_url}...')
        with urlopen(remote_file_url) as response:
            with open(local_file_path, 'wb') as data_file:
                shutil.copyfileobj(response, data_file)
        downloaded = True
        print('Done downloading data file.')

    return (local_file_path, downloaded)

def show_video(title: str, video: np.ndarray) -> None:
    """Shows the given video.

    Args:
        video (np.ndarray): The video, of shape (frames, height, width).
    """
    frame_count = video.shape[0]
    for frame_index in range(frame_count):
        frame = video[frame_index]
        cv2.imshow(f'{title} - press "q" to exit', frame)
        key = cv2.waitKey(30) & 0xFF
        if key == ord('q'):
            break

    cv2.destroyAllWindows()

def save_video(file_path: str, video: np.ndarray) -> None:
    """Saves the given video as an AVI file.

    Args:
        file_path (str): The path to the file.
        video (np.ndarray): The video to save, of shape (frames, width, height).
    """
    frame_count = video.shape[0]
    size_y = video.shape[1]
    size_x = video.shape[2]

    fourcc = cv2.VideoWriter_fourcc(*'XVID')
    out = cv2.VideoWriter(
        file_path, fourcc, fps=30.0, frameSize=(size_x, size_y)
    )

    for frame_index in range(frame_count):
        frame = video[frame_index]
        frame = cv2.cvtColor(frame, cv2.COLOR_GRAY2BGR)
        out.write(frame)

    out.release()

max_stars_repo_name: bstollnitz/grad-school-portfolio

max_issues_repo_name: bstollnitz/portfolio

max_forks_repo_name: bstollnitz/portfolio


text:
import os

from urllib.request import urlretrieve
from tqdm import tqdm
import zipfile
import hashlib
import shutil

import pandas as pd
import re
import pickle
import tensorflow as tf

# % matplotlib inline
# % config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import time
import datetime
import numpy as np
from sklearn.model_selection import train_test_split
import random

os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"

# if __name__ == "__main__":
#     pass

# 数据预处理
# 数据预处理
# 数据预处理
def load_data():
    """
    Load Dataset from File
    """
    # 读取User数据
    users_title = ['UserID', 'Gender', 'Age', 'JobID', 'Zip-code']
    users = pd.read_table('./ml-1m/users.dat', sep='::', header=None, names=users_title, engine='python')
    users = users.filter(regex='UserID|Gender|Age|JobID')
    users_orig = users.values
    # 改变User数据中性别和年龄
    gender_map = {'F': 0, 'M': 1}
    users['Gender'] = users['Gender'].map(gender_map)

    age_map = {val: ii for ii, val in enumerate(set(users['Age']))}
    users['Age'] = users['Age'].map(age_map)

    # 读取Movie数据集
    movies_title = ['MovieID', 'Title', 'Genres']
    movies = pd.read_table('./ml-1m/movies.dat', sep='::', header=None, names=movies_title, engine='python')
    movies_orig = movies.values
    # 将Title中的年份去掉
    pattern = re.compile(r'^(.*)\((\d+)\)$')

    title_map = {val: pattern.match(val).group(1) for ii, val in enumerate(set(movies['Title']))}
    movies['Title'] = movies['Title'].map(title_map)

    # 电影类型转数字字典
    genres_set = set()
    for val in movies['Genres'].str.split('|'):
        genres_set.update(val)

    genres_set.add('<PAD>')
    genres2int = {val: ii for ii, val in enumerate(genres_set)}

    # 将电影类型转成等长数字列表，长度是18
    genres_map = {val: [genres2int[row] for row in val.split('|')] for ii, val in enumerate(set(movies['Genres']))}

    for key in genres_map:
        for cnt in range(max(genres2int.values()) - len(genres_map[key])):
            genres_map[key].insert(len(genres_map[key]) + cnt, genres2int['<PAD>'])

    movies['Genres'] = movies['Genres'].map(genres_map)

    # 电影Title转数字字典
    title_set = set()
    for val in movies['Title'].str.split():
        title_set.update(val)

    title_set.add('<PAD>')
    title2int = {val: ii for ii, val in enumerate(title_set)}

    # 将电影Title转成等长数字列表，长度是15
    title_count = 15
    title_map = {val: [title2int[row] for row in val.split()] for ii, val in enumerate(set(movies['Title']))}

    for key in title_map:
        for cnt in range(title_count - len(title_map[key])):
            title_map[key].insert(len(title_map[key]) + cnt, title2int['<PAD>'])

    movies['Title'] = movies['Title'].map(title_map)

    # 读取评分数据集
    ratings_title = ['UserID', 'MovieID', 'ratings', 'timestamps']
    ratings = pd.read_table('./ml-1m/ratings.dat', sep='::', header=None, names=ratings_title, engine='python')
    ratings = ratings.filter(regex='UserID|MovieID|ratings')

    # 合并三个表
    data = pd.merge(pd.merge(ratings, users), movies)

    # 将数据分成X和y两张表
    target_fields = ['ratings']
    features_pd, targets_pd = data.drop(target_fields, axis=1), data[target_fields]

    features = features_pd.values
    targets_values = targets_pd.values

    return title_count, title_set, genres2int, features, targets_values, ratings, users, movies, data, movies_orig, users_orig

title_count, title_set, genres2int, features, targets_values, ratings, users, movies, data, movies_orig, users_orig = load_data()

pickle.dump((title_count, title_set, genres2int, features, targets_values, ratings, users, movies, data, movies_orig, users_orig), open('preprocess.p', 'wb'))

# 显示预处理后的数据
print(users.head())

print(movies.head())

print(movies.values[0])

# 从本地读取数据

title_count, title_set, genres2int, features, targets_values, ratings, users, movies, data, movies_orig, users_orig = pickle.load(open('preprocess.p', mode='rb'))


text:
# <Copyright 2022, Argo AI, LLC. Released under the MIT license.>

"""Unit tests for motion forecasting metrics."""

from typing import Final

import numpy as np
import pytest

import av2.datasets.motion_forecasting.eval.metrics as metrics
from av2.utils.typing import NDArrayFloat, NDArrayNumber

# Build stationary GT trajectory at (0, 0)
test_N: Final[int] = 10
_STATIONARY_GT_TRAJ = np.zeros((test_N, 2))

# Case 1: K=1 forecast stationary at (1, 1)
forecasted_trajectories_stationary_k1 = np.ones((1, test_N, 2))
expected_ade_stationary_k1 = np.full((1,), np.sqrt(2))
expected_fde_stationary_k1 = np.full((1,), np.sqrt(2))

# Case 2: K=6 forecasts stationary at (1, 1)
forecasted_trajectories_stationary_k6 = np.ones((6, test_N, 2))
expected_ade_stationary_k6 = np.full((6,), np.sqrt(2))
expected_fde_stationary_k6 = np.full((6,), np.sqrt(2))

# Case 3: K=1 forecast in straight line on X axis
forecasted_trajectories_straight_k1 = np.stack([np.arange(test_N), np.zeros(test_N)], axis=1)[np.newaxis, ...]  # 1xNx2
expected_ade_straight_k1 = np.full((1,), np.arange(test_N).mean())
expected_fde_straight_k1 = np.full((1,), test_N - 1)

# Case 4: K=6 forecasts in straight line on X axis
forecasted_trajectories_straight_k6: NDArrayFloat = np.concatenate(  # type: ignore
    [forecasted_trajectories_straight_k1] * 6, axis=0
)  # 6xNx2
expected_ade_straight_k6 = np.full((6,), np.arange(test_N).mean())
expected_fde_straight_k6 = np.full((6,), test_N - 1)

# Case 5: K=1 forecast in diagonal line
forecasted_trajectories_diagonal_k1 = np.stack([np.arange(test_N), np.arange(test_N)], axis=1)[np.newaxis, ...]  # 1xNx2
expected_ade_diagonal_k1 = np.full((1,), 6.36396103)
expected_fde_diagonal_k1 = np.full((1,), np.hypot(test_N - 1, test_N - 1))

@pytest.mark.parametrize(
    "forecasted_trajectories, expected_ade",
    [
        (forecasted_trajectories_stationary_k1, expected_ade_stationary_k1),
        (forecasted_trajectories_stationary_k6, expected_ade_stationary_k6),
        (forecasted_trajectories_straight_k1, expected_ade_straight_k1),
        (forecasted_trajectories_straight_k6, expected_ade_straight_k6),
        (forecasted_trajectories_diagonal_k1, expected_ade_diagonal_k1),
    ],
    ids=["stationary_k1", "stationary_k6", "straight_k1", "straight_k6", "diagonal_k1"],
)
def test_compute_ade(forecasted_trajectories: NDArrayNumber, expected_ade: NDArrayFloat) -> None:
    """Test that compute_ade returns the correct output with valid inputs.

    Args:
        forecasted_trajectories: Forecasted trajectories for test case.
        expected_ade: Expected average displacement error.
    """
    ade = metrics.compute_ade(forecasted_trajectories, _STATIONARY_GT_TRAJ)
    np.testing.assert_allclose(ade, expected_ade)  # type: ignore

@pytest.mark.parametrize(
    "forecasted_trajectories, expected_fde",
    [
        (forecasted_trajectories_stationary_k1, expected_fde_stationary_k1),
        (forecasted_trajectories_stationary_k6, expected_fde_stationary_k6),
        (forecasted_trajectories_straight_k1, expected_fde_straight_k1),
        (forecasted_trajectories_straight_k6, expected_fde_straight_k6),
        (forecasted_trajectories_diagonal_k1, expected_fde_diagonal_k1),
    ],
    ids=["stationary_k1", "stationary_k6", "straight_k1", "straight_k6", "diagonal_k1"],
)
def test_compute_fde(forecasted_trajectories: NDArrayNumber, expected_fde: NDArrayFloat) -> None:
    """Test that compute_fde returns the correct output with valid inputs.

    Args:
        forecasted_trajectories: Forecasted trajectories for test case.
        expected_fde: Expected final displacement error.
    """
    fde = metrics.compute_fde(forecasted_trajectories, _STATIONARY_GT_TRAJ)
    assert np.array_equal(fde, expected_fde)

@pytest.mark.parametrize(
    "forecasted_trajectories, miss_threshold_m, expected_is_missed_label",
    [
        (forecasted_trajectories_stationary_k1, 2.0, False),
        (forecasted_trajectories_stationary_k6, 2.0, False),
        (forecasted_trajectories_straight_k1, expected_fde_straight_k1[0] + 0.01, False),
        (forecasted_trajectories_straight_k1, expected_fde_straight_k1[0] - 0.01, True),
        (forecasted_trajectories_diagonal_k1, 2.0, True),
    ],
    ids=["stationary_k1", "stationary_k6", "straight_below_threshold", "straight_above_threshold", "diagonal"],
)
def test_compute_is_missed_prediction(
    forecasted_trajectories: NDArrayNumber, miss_threshold_m: float, expected_is_missed_label: bool
) -> None:
    """Test that compute_is_missed_prediction returns the correct output with valid inputs.

    Args:
        forecasted_trajectories: Forecasted trajectories for test case.
        miss_threshold_m: Minimum distance threshold for final displacement to be considered a miss.
        expected_is_missed_label: Expected is_missed label for test case.
    """
    is_missed_prediction = metrics.compute_is_missed_prediction(
        forecasted_trajectories, _STATIONARY_GT_TRAJ, miss_threshold_m
    )

    # Check that is_missed labels are of the correct shape and have the correct value
    assert is_missed_prediction.shape == forecasted_trajectories.shape[:1]
    assert np.all(is_missed_prediction == expected_is_missed_label)


text:
import numpy as np
import torch
import torch.nn as nn

from encoder import ResBlock, rnn_encoder
from decoder import rnn_decoder

from processor import preprocess, postprocess

class seq2seq(nn.Module):

    def __init__(self, device, teacher_forcing_ratio=0., dropout=0., random_init=False):
        super(seq2seq, self).__init__()
        self.random_init = random_init
        self.teacher_forcing_ratio = teacher_forcing_ratio

        self.encoder = rnn_encoder(device=device, feature_size=250, encoder_hidden=256, dropout=dropout)
        self.decoder = rnn_decoder(device=device, decoder_hidden=256,teacher_forcing_ratio=self.teacher_forcing_ratio,
                                   dropout=dropout, random_init=random_init)

    def forward(self, inputs, targets):
        contexts, encoder_state = self.encoder(inputs)
        outputs, final_state = self.decoder(targets, contexts.transpose(1, 2), encoder_state)
        return outputs

    def sample(self, inputs, return_attns=False):
        contexts, encoder_state = self.encoder(inputs)
        outputs, attns = self.decoder.sample(contexts.transpose(1, 2), encoder_state)
        return outputs if not return_attns else (outputs, attns)

class CNN_en1D(nn.Module):
    def __init__(self, L, ks=71):
        super(CNN_en1D, self).__init__()
        self.L = L
        self.C = 64
        self.cnn_g1 = nn.Sequential(nn.Conv1d(1, self.C, ks, 1, padding=(ks-1)//2),
                                    nn.BatchNorm1d(self.C),
                                    nn.ReLU(self.C),
                                    ResBlock(self.C, ks=ks, norm_func=nn.BatchNorm1d),
                                    ResBlock(self.C, ks=ks, norm_func=nn.BatchNorm1d))

        self.cnn_g2 = nn.Sequential(ResBlock(self.C, ks=ks, norm_func=nn.BatchNorm1d),
                                    ResBlock(self.C, ks=ks, norm_func=nn.BatchNorm1d),
                                    nn.Conv1d(self.C, 1, ks, 1, padding=(ks-1)//2))

        for m in self.modules():
            if isinstance(m, nn.Conv1d):
                nn.init.kaiming_normal_(m.weight)
            elif isinstance(m, nn.BatchNorm1d):
                nn.init.constant_(m.weight, 1)
                nn.init.constant_(m.bias, 0)
            elif isinstance(m, nn.Linear):
                nn.init.xavier_uniform_(m.weight)

    def forward(self, input):
        cnn_out = self.cnn_g1(input)
        cnn_out = self.cnn_g2(cnn_out)
        cnn_out = cnn_out.view(-1, self.L)

        return torch.sigmoid(cnn_out)

class qrs_detector(object):
    def __init__(self, model_path, model_type='fcn', fs=500):
        assert model_type in ['fcn', 'seq2seq']
        self.model_type=model_type
        self.fs = fs
        self.device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')
        if self.model_type=='fcn':
            self.model = CNN_en1D(L=5000).to(self.device)
        else:
            self.model = seq2seq(device=self.device).to(self.device)

        #print('Loading model from {}.'.format(model_path))

        checkpoint = torch.load(model_path)
        self.model.load_state_dict(checkpoint['state_dict'])
        self.model.eval()

        # print('Epoch: {}，Se: {}， pre: {}, F1: {}, acc:{}'.
        #       format(checkpoint['epoch'], checkpoint['se'], checkpoint['pre'], checkpoint['f1'],
        #              checkpoint['cpsc_acc']))

    def detection(self, signal):
        signal = preprocess(signal, fs=self.fs)
        signal_tensor = torch.from_numpy(signal).view(1, 1, -1).float().to(self.device)

        if self.model_type=='fcn':
            outputs = self.model(signal_tensor)  # outputs: 1 x 1 x beats
            beat_locs2 = list(outputs.squeeze().detach().cpu().numpy())
            beat_locs = postprocess(beat_locs2, margin=12)
        else:
            outputs = self.model.sample(signal_tensor)  # outputs: 1 x 1 x beats
            beat_locs = np.rint(outputs.view(-1).detach().cpu().numpy()).astype(np.int) if outputs \
                                                                                is not None else np.array([])
        return beat_locs


text:
import numpy as np
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()
import os
import image_utils
from random import sample
from training_params import batch_size
import training_params

count_by_folders = {}
image_paths = {}
for root, dir, files in os.walk(image_utils.image_dir):
    count = 0
    for file in files:
        if file.split('.')[-1] == 'jpg':
            count += 1
            if root in image_paths:
                    image_paths[root].append(file)
            else:
                    image_paths[root] = [file]
            # display_image(unpreprocess_image(load_image(os.path.join(root, file))))
    count_by_folders[root] = count

singleton_folders = []
multiple_folders = []
for folder in count_by_folders:
        if count_by_folders[folder] == 1:
                singleton_folders.append(folder)
        elif count_by_folders[folder] > 1:
                multiple_folders.append(folder)
        else:
                print("EMPTY FOLDER: ", folder)

def get_batch():
        number_of_pairs = ((batch_size + 1) // 2) // 2
        singles = sample(singleton_folders, k=batch_size - number_of_pairs * 2)
        multiples = sample(multiple_folders, k=number_of_pairs)
        batch = []
        for m in multiples:
                ms = sample(image_paths[m], k=2)
                # print(ms)
                for filename in ms:
                        batch.append(image_utils.load_image(os.path.join(m, filename), training_params.shape))
        for s in singles:
                batch.append(image_utils.load_image(os.path.join(s, image_paths[s][0]), training_params.shape))
        return np.array(batch)

if __name__ == '__main__':
        batch = get_batch()
        print(batch.shape)
        image_utils.display_image(image_utils.unpreprocess_image(batch[0]))
        image_utils.display_image(image_utils.unpreprocess_image(batch[1]))


text:
from ctypes import c_int, c_float, pointer, POINTER
import numpy as np
from .load_libmatting import load_libmatting

library = load_libmatting()

# Need those types for function signatures:
c_int_p = POINTER(c_int)
c_float_p = POINTER(c_float)

# Declare function signatures
_knn = library.knn
_knn.argtypes = [
    c_float_p,
    c_float_p,
    c_int_p,
    c_float_p,
    c_int,
    c_int,
    c_int,
    c_int]

def knn(data_points, query_points, k, overwrite_data_points=False):
    if data_points.dtype == np.float64:
        data_points = data_points.astype(np.float32)

    if query_points.dtype == np.float64:
        query_points = query_points.astype(np.float32)

    assert(len(data_points.shape) == 2)
    assert(len(query_points.shape) == 2)
    assert(data_points.dtype == np.float32)
    assert(query_points.dtype == np.float32)
    assert(len(data_points) >= k)

    n_data_points, data_points_dim = data_points.shape
    n_query_points, query_points_dim = query_points.shape

    assert(data_points_dim == query_points_dim)

    indices = np.empty(n_query_points * k, dtype=np.int32)
    squared_distances = np.empty(n_query_points * k, dtype=np.float32)

    data_points = data_points.ravel()
    query_points = query_points.ravel()

    if not overwrite_data_points:
        data_points = data_points.copy()

    assert(data_points.flags['C_CONTIGUOUS'])
    assert(query_points.flags['C_CONTIGUOUS'])

    _knn(
        np.ctypeslib.as_ctypes(data_points),
        np.ctypeslib.as_ctypes(query_points),
        np.ctypeslib.as_ctypes(indices),
        np.ctypeslib.as_ctypes(squared_distances),
        n_data_points,
        n_query_points,
        data_points_dim,
        k)

    indices = indices.reshape(n_query_points, k)
    squared_distances = squared_distances.reshape(n_query_points, k)

    return indices

def test():
    for d in range(1, 10):
        for k in range(1, 5):
            na = np.random.randint(1, 100)
            nb = np.random.randint(1, 100)

            if k > na:
                continue

            a = np.random.randn(na, d).astype(np.float32)
            b = np.random.randn(nb, d).astype(np.float32)

            indices1 = knn(a, b, k)
            indices2 = np.empty((nb, k), dtype=np.int32)

            for i in range(nb):
                distances = np.linalg.norm(a - b[i].reshape(1, d), axis=1)
                indices2[i] = np.argsort(distances)[:k]

            assert(indices1.shape == indices2.shape)
            assert(np.allclose(indices1, indices2))
    print("tests passed")

if __name__ == "__main__":
    test()

max_stars_repo_name: styanddty/matting

max_issues_repo_name: styanddty/matting


text:
#! /usr/bin/env python

import numpy as np
import time

def czt(x, m = None, w = None, a = None):

    # Copyright (C) 2004 Daniel Gunyan
    #
    # This program 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 3 of the License, or (at your option) any later
    # version.
    #
    # This program 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
    # this program; if not, see <http://www.gnu.org/licenses/>.

    # -*- texinfo -*-
    # @deftypefn  {Function File} {} czt (@var{x})
    # @deftypefnx {Function File} {} czt (@var{x}, @var{m})
    # @deftypefnx {Function File} {} czt (@var{x}, @var{m}, @var{w})
    # @deftypefnx {Function File} {} czt (@var{x}, @var{m}, @var{w}, @var{a})
    # Chirp z-transform.  Compute the frequency response starting at a and
    # stepping by w for m steps.  a is a point in the complex plane, and
    # w is the ratio between points in each step (i.e., radius increases
    # exponentially, and angle increases linearly).
    #
    # To evaluate the frequency response for the range f1 to f2 in a signal
    # with sampling frequency Fs, use the following:
    #
    # @example
    # @group
    # m = 32;                               ## number of points desired
    # w = exp(-j*2*pi*(f2-f1)/((m-1)*Fs));  ## freq. step of f2-f1/m
    # a = exp(j*2*pi*f1/Fs);                ## starting at frequency f1
    # y = czt(x, m, w, a);
    # @end group
    # @end example
    #
    # If you don't specify them, then the parameters default to a fourier
    # transform:
    #     m=length(x), w=exp(-j*2*pi/m), a=1
    #
    # If x is a matrix, the transform will be performed column-by-column.
    # @end deftypefn

    # Algorithm (based on Oppenheim and Schafer, "Discrete-Time Signal
    # Processing", pp. 623-628):
    #   make chirp of length -N+1 to max(N-1,M-1)
    #     chirp => w^([-N+1:max(N-1,M-1)]^2/2)
    #   multiply x by chirped a and by N-elements of chirp, and call it g
    #   convolve g with inverse chirp, and call it gg
    #     pad ffts so that multiplication works
    #     ifft(fft(g)*fft(1/chirp))
    #   multiply gg by M-elements of chirp and call it done

    n = len(x)

    if m is None: m = n
    if w is None: w = np.exp(-2j * np.pi / m)
    if a is None: a = 1

    wExponents1 =   np.arange(    0   , n) ** 2 / 2.0    # n elements         [ 0 .. (n - 1) ]
    wExponents2 = - np.arange(-(n - 1), m) ** 2 / 2.0    # m + n - 1 elements [ -(n - 1) .. +(m - 1) ]
    wExponents3 =   np.arange(    0   , m) ** 2 / 2.0    # m elements         [ 0        ..  (m - 1) ]

    xx = x * a ** -np.arange(n) * w ** wExponents1

    # Determine next-biggest FFT of power-of-two

    nfft = 1
    while nfft < (m + n - 1):
        nfft += nfft

    fxx = np.fft.fft(xx, nfft)

    ww = w ** wExponents2

    fww = np.fft.fft(ww, nfft)

    fyy = fxx * fww

    yy = np.fft.ifft(fyy, nfft)

    # select output

    yy = yy[n - 1 : m + n - 1]

    y = yy * w ** wExponents3

    return y

NUM_POINTS = 3

x = (10.0 + np.arange(NUM_POINTS)) + (20.0 + np.arange(NUM_POINTS) ** 2) * 1j

t1 = time.time()
numpy_fft_result = np.fft.fft(x)
t2 = time.time()
numpy_fft_time = (t2 - t1)
print "numpy_fft_time", numpy_fft_time

t1 = time.time()
czt_result = czt(x)
t2 = time.time()
czt_time = (t2 - t1)

err_czt = czt_result - np.fft.fft(x)
err_czt = sum(err_czt.conj() * err_czt)

print "czt_time", czt_time
print "czt_error", abs(err_czt)

max_stars_repo_name: sidneycadot/FourierBenchmarks

max_issues_repo_name: sidneycadot/FourierBenchmarks

max_forks_repo_name: sidneycadot/FourierBenchmarks


text:
#!/user/bin/env python
"""Python binary to train COVID-19 epidemic model with single CPU."""
from typing import Text
from absl import app
from absl import flags
import io_utils
import json
import logging
import numpy as np
import os
import resample
# pylint: disable=import-not-at-top
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
import combined_model

NUM_STATES = 58
flags.DEFINE_float("learning_rate", 0.01, "Initial learning rate.")
flags.DEFINE_integer("max_epochs", 2000, "Number of steps to run trainer.")
flags.DEFINE_integer("test_duration", 100, "Number of days to predict")
flags.DEFINE_string("metadata", None, "")
flags.DEFINE_string("residuals", None, "Only used in resample jobs.")
flags.DEFINE_string("output_path", None, "Path to store the output data.")
flags.DEFINE_integer("min_t0", 1, "value for which t0 iteration starts.")
flags.DEFINE_integer("max_t0", 21, "value for which t0 iteration ends.")
flags.DEFINE_string("file_index_plus_one", None, "Index for permutation.")
flags.DEFINE_bool(
    "output_to_json", True, "If the export data are saved in json file.")
flags.DEFINE_bool(
    "resample_jobs", False, "If the model estimation is on resampled data.")
flags.DEFINE_bool("flatten_future", True, "Indicates whether or not the "
                  "future prediction relies on the flat death and infection "
                  "rates.")
FLAGS = flags.FLAGS

def read_metadata(fpath: Text, n_states: int):
  with open(fpath, "r") as f:
    jdata = json.load(f)
  if len(jdata.keys()) != n_states:
    raise ValueError(f"The metadata should contains exactly {n_states} records.")
  for k, v in jdata.items():
    if set(v.keys()) != {"state", "knot", "data"}:
      raise ValueError(f"At the record {k}, the json data don't follow the required format, "
                       "expected to have three fields: state, knot, and data.")
  return jdata

def map_keys(file_index: int, n_states: int, resample_jobs: bool = False):
  int_file_index = int(file_index)
  if resample_jobs:
    return str(int_file_index % n_states), str(int_file_index // n_states)
  else:
    if int_file_index >= n_states:
      raise ValueError(f"The file index should be between 0 and {n_states}.")
    return file_index, None

def main(unused_argv):
  logger = logging.getLogger("Covid-19-estimation")
  job_metadata = read_metadata(FLAGS.metadata, NUM_STATES)
  job_file_index = str(int(FLAGS.file_index_plus_one) - 1)
  state_key, resample_key = map_keys(
      job_file_index, NUM_STATES, FLAGS.resample_jobs)
  job_state = job_metadata[state_key]["state"]
  logger.info(f"Create estimation job for the state {job_state}.")

  # Parse the required data & args for running the estimation model.
  job_knots = job_metadata[state_key]["knot"]
  if FLAGS.resample_jobs:
    job_residuals = resample.read_residuals(
        FLAGS.residuals, NUM_STATES)
    job_input_data = resample.get_resampled_input(
        job_residuals[state_key]["fitted"],
        job_residuals[state_key]["residual"],
        int(resample_key))
  else:
    job_input_data = job_metadata[state_key]["data"]
  job_knots_connect = [1] * len(job_knots)
  job_initial_a = [0.2] * (len(job_knots) + 1)

  # Step 1 estimation (infection cases only).
  estimator = combined_model.Covid19CombinedEstimator(
      knots=job_knots,
      knots_connect=job_knots_connect,
      loss_prop=0,
      estimator_args={
          "learning_rate": FLAGS.learning_rate,
          "epochs": FLAGS.max_epochs},
      initial_guess_a=np.array(job_initial_a),
      variable_death_rate_trainable=False
  )
  estimator.fit(data=job_input_data, min_t0=FLAGS.min_t0, max_t0=FLAGS.max_t0)
  # Save the estimated weights from step 1 (will be fixed in step 2).
  stage1_estimated_a, stage1_estimated_t0 = estimator.final_model.a.numpy(), estimator.final_model.t0
  logger.info("First stage estimation done.")

  # Step 2 estimation (death only).
  estimator = combined_model.Covid19CombinedEstimator(
      knots=job_knots,
      knots_connect=job_knots_connect,
      loss_prop=1,
      estimator_args={
          "learning_rate": FLAGS.learning_rate,
          "epochs": FLAGS.max_epochs},
      initial_guess_a=stage1_estimated_a,
      variable_a_trainable=False
  )
  estimator.fit(data=job_input_data, min_t0=stage1_estimated_t0, max_t0=stage1_estimated_t0)
  io_utils.parse_estimated_model(estimator)
  if FLAGS.resample_jobs:
    job_suffix = "state_" + state_key + "resample_" + resample_key
  else:
    job_suffix = "state_" + state_key
  io_utils.export_estimation_and_prediction(
      estimator=estimator,
      test_duration=FLAGS.test_duration,
      output_path=FLAGS.output_path,
      suffix=job_suffix,
      flatten_future=FLAGS.flatten_future,
      to_json=FLAGS.output_to_json
  )
  logger.info("Second stage estimation done.")

if __name__ == "__main__":
  flags.mark_flag_as_required("metadata")
  flags.mark_flag_as_required("output_path")
  flags.mark_flag_as_required("file_index_plus_one")
  app.run(main)


text:
import json
import logging
import math
import queue
import socket
import sys
import threading
import time
from multiprocessing.sharedctypes import Value
from sqlite3 import Time

import cv2
import numpy as np
from cscore import CameraServer
from networktables import NetworkTables, NetworkTablesInstance

#Below link has a barebones version, useful for getting the camera server stuff
#https://docs.wpilib.org/en/stable/docs/software/vision-processing/wpilibpi/basic-vision-example.html

# Queue of Packets
# Thread Safe.. Packets being sent to robot are placed here!
PacketQueue = queue.Queue()

# Creates a socket
#No clue if this will work, just pulled it from last year
class SocketWorker(threading.Thread):
    def __init__(self, q, *args, **kwargs):
        self.queue = q
        super().__init__(*args, **kwargs)

        # Initialize Socket Connect
        SocketHost = "10.1.38.2"
        SocketPort = 5800
        self.sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM)
        self.sock.connect((SocketHost, SocketPort))

    def run(self):
        while True:
            try:
                packet = self.queue.get()
                # Convert Packet to JSON
                # Send Message Over Socket
                try:
                    data_string = json.dumps(packet)
                    self.sock.sendall(data_string.encode())
                except Exception as e:
                    print("Socket Exception " + str(e))
            except Exception as e1:
                pass

def calculateDistanceFeet(targetPixelWidth):
    # d = Tft*FOVpixel/(2*Tpixel*tanΘ)

    #Target width in feet, ball for 2022 is ~9 1/2 inches (so tft is just an estimate to 9 inches)
    Tft = .75
    #170 is our cameras FOV for 2022. PSEye was 75
    tanFOV = math.tan(170/2)
    #For some reason, 640 works for a 320 resolution width, if set to 320 it may remove the need for the arbitrary value
    distEst = Tft * 640 / (2 * targetPixelWidth * tanFOV)
    # Unsure as to what measurement distEst is producing in the above line, but multiplying it by .32 will return your distance in feet
    distEstFeet = distEst * .32
    return (abs(distEstFeet))

if __name__ == "__main__":
    #Avoid touching camera server settings
    print('2022 Ball Vision Blue Starting')

    cameraConfig = ''
    with open('/boot/frc.json', 'r') as f:
        mycamera = json.load(f)
        cameraConfig = mycamera
        camera = cameraConfig['cameras'][0]
        '''
        Create this list of values by going to wpilibpi.local, vision settings, open stream, change settings to whats needed, open "source config JSON",
        then paste it into "Custom Properties Json" in Vision Settings. Save it, then copy the list created, and replace whatevers set in the cameraConfig variable
        Dont try to create this manually, its extremely picky on formatting.
        '''
        cameraConfig = {"fps":120,"height":240,"pixel format":"mjpeg","properties":[{"name":"connect_verbose","value":1},{"name":"raw_brightness","value":-8},{"name":"brightness","value":43},{"name":"raw_contrast","value":0},{"name":"contrast","value":0},{"name":"raw_saturation","value":128},{"name":"saturation","value":100},{"name":"raw_hue","value":0},{"name":"hue","value":50},{"name":"white_balance_temperature_auto","value":True},{"name":"gamma","value":100},{"name":"raw_gain","value":0},{"name":"gain","value":0},{"name":"power_line_frequency","value":1},{"name":"white_balance_temperature","value":4600},{"name":"raw_sharpness","value":2},{"name":"sharpness","value":33},{"name":"backlight_compensation","value":1},{"name":"exposure_auto","value":3},{"name":"raw_exposure_absolute","value":157},{"name":"exposure_absolute","value":3},{"name":"exposure_auto_priority","value":True}],"width":320}

    cs = CameraServer.getInstance()
    cameraSettings = cs.startAutomaticCapture()

    cameraSettings.setConfigJson(json.dumps(cameraConfig))
    input_stream = cs.getVideo()
    output_stream = cs.putVideo('Processed', 320, 240)

    SocketThread = SocketWorker(PacketQueue).start()

    #Numpy creates an array of zeros in the size of the image width/height. Its mentioned in documentation this can be performance intensive, and to do it outside the loop
    #Note that the order is (height, width) in shape
    imgForm = np.zeros(shape=(240, 320, 3), dtype=np.uint8)

    #Blue ball
    blueHue = [81, 122]
    blueSat = [0, 255]
    blueVal = [76, 255]

    #Creating settings for blur filter. Radius should be figured out by testing in GRIP
    radius = 8
    ksize = int(2 * round(radius) + 1)

    #Parameters for targeting, I set these all up here because its easier to go through and change them when tuning with grip
    #500
    cnt_area_low = 300
    #cnt_area_high = 7500
    minimum_perimeter = 9
    width_minimum = 9

    #300
    width_maximum = 500
    height_minimum = 9
    #300
    height_maximum = 500
    solid_Low = 94
    solid_High = 100
    #18
    min_vertices = 5
    max_vertices = 100
    rat_low = 0
    rat_high = 3
    cy = ''
    cx = ''
    solid = 0
    last_cnt_area = 0
    conCount = 0
    cnt_to_process = 0

    #Used to pick the lowest detected contour in the image which is most likely the closest ball
    lowest_y = 1000

    last_contour_x = 1000
    last_contour_y = 1000
    current_frame = 0

    #Create empty dictionary for values we'll send in our packet
    PacketValue = {}

    '''
    #Values used for testing framerate
    curTime = time.time()
    oldTime = ''
    '''
    printCount = 1
    myDistFeet = 0

    print('Blue ball vision setup complete')

    while True:
        #Try covers the following code to make sure we never fail during a match.
        try:
            current_frame += 1
            '''
            if current_frame % 60 == 0:
                oldTime = curTime
                curTime = time.time()
                print(curTime - oldTime)
            '''

            PacketValue = {}
            PacketValue['cameraid'] = 0
            PacketValue['ballColor'] = 'Blue'

            #start_time = time.time() #Use this to get FPS below
            frame_time, input_img = input_stream.grabFrame(imgForm)

            # Notify output of error and skip iteration
            if frame_time == 0:
                output_stream.notifyError(input_stream.getError())
                continue

            if current_frame == 1000:
                cv2.imwrite('blueoriginal.jpeg', input_img)

            #Blur image, change image colorspace to HSV, blur it, and for the 2022 robot we flip the image due to the camera being upside down.
            #cv2.imwrite('original.jpeg', input_img)
            #You have to blur before converting color space!! If you don't, the output will be messed up.
            input_img = cv2.blur(input_img, (ksize, ksize))
            #cv2.imwrite('blur.jpeg', input_img)
            input_img = cv2.cvtColor(input_img, cv2.COLOR_BGR2HSV)
            #cv2.imwrite('2hsv.jpeg', input_img)
            #Color bumper area to black
            #If camera is flipped
            '''
            input_img[0:30, 0:320] = (0,0,0)
            input_img = cv2.flip(input_img, 1)
            '''

            #If camera is not flipped
            input_img[210:240, 0:320] = (0,0,0)

            #Mask out colors that dont fall in the range we'd find the blue ball in
            mask = cv2.inRange(input_img, (blueHue[0], blueSat[0], blueVal[0]),
                                (blueHue[1], blueSat[1], blueVal[1]))
            #cv2.imwrite('masked.jpeg', mask)

            if current_frame == 1000:
                cv2.imwrite('bluethresh.jpeg', input_img)
            _, contours, _ = cv2.findContours(mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_TC89_KCOS)

            # Sort contours by area size (biggest to smallest)
            cntsSorted = sorted(contours, key=lambda x: cv2.contourArea(x), reverse=True)
            #cv2.drawContours(mask, cntsSorted, -1, (0,255,0), 3)
            #cv2.imwrite('contours.jpeg', mask)

            con = []
            for cnt in cntsSorted:
                # Calculate Contour area early as a way to eliminate unnecesarily processing contours that arent viable
                cntArea = cv2.contourArea(cnt)
                if (cntArea > cnt_area_low):# and (cntArea < cnt_area_high)
                    # Get moments of contour; mainly for centroid
                    M = cv2.moments(cnt)
                    # Get convex hull (bounding polygon on contour)
                    hull = cv2.convexHull(cnt)

                    # calculate area of convex hull
                    #hullArea = cv2.contourArea(hull)

                    # Approximate shape - Approximates a polygonal curve(s) with the specified precision.
                    approximateShape = cv2.approxPolyDP(cnt, 0.01 * cv2.arcLength(cnt, True), True)

                    x, y, w, h = cv2.boundingRect(cnt)
                    #Check what amount of the bounding rectangle is filled by the detected object.
                    ratio = float(w) / h

                    perimeter = cv2.arcLength(cnt, True)

                    # finding center point of shape
                    M = cv2.moments(cnt)
                    if M['m00'] != 0.0:
                        x = int(M['m10']/M['m00'])
                        y = int(M['m01']/M['m00'])

                    #Frequently attempts to divide by zero, so a check that they aren't is necessary
                    '''
                    if cntArea != 0 and cntArea != 0:
                        solid = 100 * cntArea / hullArea
                    '''

                    #Filtering out the contours based on tuned values we are looking for
                    validCnt = True
                    validCnt &= (perimeter > minimum_perimeter)
                    validCnt &= (w >= width_minimum) and (w <= width_maximum)
                    validCnt &= (h >= height_minimum) and (h <= height_maximum)
                    '''
                    if solid != 0:
                        validCnt &= (solid > solid_Low) and (solid <= solid_High)
                    '''
                    validCnt &= (len(approximateShape) >= 8)
                    validCnt &= (ratio >= rat_low) and (ratio < rat_high)

                    #List of prints for debugging
                    '''
                    print('Perimeter:', perimeter)
                    print('Width:', w)
                    print('Height:', h)
                    print('Solid:', solid)
                    print('Approximate Shape:', approximateShape)
                    print('Ratio:', ratio)
                    print('Vetices', )
                    '''

                    validCnt &= (cv2.arcLength(cnt, True) < 10000)

                    circularity = 0
                    if(perimeter == 0):
                        validCnt = False
                    else:
                        circularity = 4*math.pi*(cntArea/(perimeter*perimeter))
                        validCnt &= (.5 < circularity < 1.5)

                    if(validCnt) and y < lowest_y :
                        myDistFeet = calculateDistanceFeet(w)
                        #print(cntArea, circularity, ratio)
                        cnt_to_process = cnt

            x, y, w, h = cv2.boundingRect(cnt_to_process)

            M = cv2.moments(cnt_to_process)
            cy = int(M["m01"] / M["m00"])
            cx = int(M["m10"] / M["m00"])

            #Every 100 cnt detections, print X + Y. Mostly just used to guestimate framerate
            if printCount % 100 == 0:
                print('X center:', cx, 'Y center:',cy)
                printCount = 1
            else:
                printCount += 1

            #The /2 and -.6 are entirely arbitrary values. If re-using this code in the future, you will need to re-sample to find those valvues.
            dist = (myDistFeet/2)-.6
            #print(dist)

            PacketValue['BallX'] = cy
            PacketValue['BallY'] = cx
            PacketValue['Dist'] = dist

            last_contour_x = cx
            last_contour_y = cy

            #Send packet, then clear values
            PacketQueue.put_nowait(PacketValue)
            cx = ''
            cy = ''
            last_cnt_area = 0
            cnt_to_process = ''
            lowest_y = 1000

            #print(PacketValue)
        except Exception as e:
            print('Error, likely that a ball wasnt found, error: ')


text:
import math
import RegExService
import random
import numpy
from functools import reduce
import sys
import getopt

alfa = 2
beta = 5
sigm = 3
ro = 0.8
th = 80
fileName = "E-n22-k4.txt"
iterations = 100
ants = 25

def generateGraph():
    capacityLimit, graph, demand, optimalValue = RegExService.getData(fileName)
    vertices = list(graph.keys())
    print(vertices)
    vertices.remove(1)

    edges = {(min(a, b), max(a, b)): numpy.sqrt((graph[a][0] - graph[b][0]) ** 2 + (graph[a][1] - graph[b][1]) ** 2) for
             a in graph.keys() for b in graph.keys()}
    feromones = {(min(a, b), max(a, b)): 1 for a in graph.keys() for b in graph.keys() if a != b}

    return vertices, edges, capacityLimit, demand, feromones, optimalValue

def solutionOfOneAnt(vertices, edges, capacityLimit, demand, feromones):
    solution = list()

    while (len(vertices) != 0):
        path = list()
        city = numpy.random.choice(vertices)
        capacity = capacityLimit - demand[city]
        path.append(city)
        vertices.remove(city)
        while (len(vertices) != 0):
            probabilities = list(map(lambda x: ((feromones[(min(x, city), max(x, city))]) ** alfa) * (
                    (1 / edges[(min(x, city), max(x, city))]) ** beta), vertices))
            probabilities = probabilities / numpy.sum(probabilities)
            city = numpy.random.choice(vertices, p=probabilities)
            capacity = capacity - demand[city]

            if (capacity >= 0):
                path.append(city)
                vertices.remove(city)
            else:
                break
        solution.append(path)
    return solution

def rateSolution(solution, edges):
    s = 0
    for i in solution:
        a = 1
        for j in i:
            b = j
            s = s + edges[(min(a, b), max(a, b))]
            a = b
        b = 1
        s = s + edges[(min(a, b), max(a, b))]
    return s

def updateFeromone(feromones, solutions, bestSolution):
    Lavg = reduce(lambda x, y: x + y, (i[1] for i in solutions)) / len(solutions)
    feromones = {k: (ro + th / Lavg) * v for (k, v) in feromones.items()}
    solutions.sort(key=lambda x: x[1])
    if (bestSolution != None):
        if (solutions[0][1] < bestSolution[1]):
            bestSolution = solutions[0]
        for path in bestSolution[0]:
            for i in range(len(path) - 1):
                feromones[(min(path[i], path[i + 1]), max(path[i], path[i + 1]))] = sigm / bestSolution[1] + feromones[
                    (min(path[i], path[i + 1]), max(path[i], path[i + 1]))]
    else:
        bestSolution = solutions[0]
    for l in range(sigm):
        paths = solutions[l][0]
        L = solutions[l][1]
        for path in paths:
            for i in range(len(path) - 1):
                feromones[(min(path[i], path[i + 1]), max(path[i], path[i + 1]))] = (sigm - (l + 1) / L ** (l + 1)) + \
                                                                                feromones[(
                                                                                min(path[i], path[i + 1]),
                                                                                max(path[i], path[i + 1]))]
    return bestSolution

def main():
    bestSolution = None
    vertices, edges, capacityLimit, demand, feromones, optimalValue = generateGraph()

    for i in range(iterations):
        solutions = list()
        for _ in range(ants):
            solution = solutionOfOneAnt(vertices.copy(), edges, capacityLimit, demand, feromones)
            solutions.append((solution, rateSolution(solution, edges)))
        bestSolution = updateFeromone(feromones, solutions, bestSolution)
        print(str(i) + ":\t" + str(int(bestSolution[1])) + "\t" + str(optimalValue))
    return bestSolution

if __name__ == "__main__":
    argv = sys.argv[1:]

    try:
        opts, args = getopt.getopt(argv, "f:a:b:s:r:t:i:n:", ["fileName=",
                                                              "alpha=", "beta=", "sigma=", "rho=", "theta=",
                                                              "iterations=", "numberOfAnts="])
    except getopt.GetoptError:
        sys.exit(2)
    for opt, arg in opts:
        if (opt in ("-a", "--alpha")):
            alfa = float(arg)
        elif (opt in ("-b", "--beta")):
            beta = float(arg)
        elif (opt in ("-s", "--sigma")):
            sigm = float(arg)
        elif (opt in ("-r", "--rho")):
            ro = float(arg)
        elif (opt in ("-t", "--theta")):
            th = float(arg)
        elif (opt in ("-f", "--fileName", "--file")):
            fileName = str(arg)
        elif (opt in ("-i", "--iterations")):
            iterations = int(arg)
        elif (opt in ("-n", "--numberOfAnts")):
            ants = int(arg)

    print("file name:\t" + str(fileName) +
          "\nalpha:\t" + str(alfa) +
          "\nbeta:\t" + str(beta) +
          "\nsigma:\t" + str(sigm) +
          "\nrho:\t" + str(ro) +
          "\ntheta:\t" + str(th) +
          "\niterations:\t" + str(iterations) +
          "\nnumber of ants:\t" + str(ants))

    solution = main()

    print("Solution: " + str(solution))
    if fileName == "E-n22-k4.txt":
        optimalSolution = ([[18, 21, 19, 16, 13], [17, 20, 22, 15], [14, 12, 5, 4, 9, 11], [10, 8, 6, 3, 2, 7]], 375)
        print("Optimal solution: " + str(optimalSolution))