text:
"""Create parameters for the problem instance"""

from core import extract_info as ext
from classes.searcher import MySearcher
from classes.belief import MyBelief
from classes.target import MyTarget
from classes.solver_data import MySolverData
from classes.inputs import MyInputs

# external packages
from igraph import *
import numpy as np
import pickle
import random

# init inputs class
def default_specs():
    """Initialize pre-set parameters
    If needed, change parameters here using MyInputs() class functions
    Return: specs"""

    specs = MyInputs()
    # -------------------
    # call here MyInputs() functions to change default specs
    # -------------------

    return specs

# init classes from specs (exp_inputs)
def create_solver_data(specs):
    """Initialize solver data class from specs (exp_inputs)"""
    g = specs.graph

    # planning stuff
    deadline = specs.deadline
    theta = specs.theta
    horizon = specs.horizon
    solver_type = specs.solver_type
    timeout = specs.timeout

    solver_data = MySolverData(horizon, deadline, theta, g, solver_type, timeout)

    return solver_data

def create_belief(specs):
    """Initialize belief
    either from distribution (uniform, split among random vertices)
    or from given b0
    """

    if specs.b0 is None:
        g = specs.graph
        # type of distribution (default: uniform)
        type_distribution = specs.belief_distribution
        if specs.start_target_v_list is None:
            # does not have a list of possible target start vertices
            v_list = placement_list(specs, 't')
        else:
            v_list = specs.start_target_v_list
        # create b_0 (list with probabilities, b_c = b[0] = 0)
        b_0 = set_initial_belief(g, v_list, type_distribution)
    else:
        # user-defined initial belief (with b_c)
        b_0 = specs.b0

    # create belief obj
    belief = MyBelief(b_0)

    # save on specs
    v_list = v_list_from_belief(b_0)
    specs.set_start_target(v_list)

    return belief

def create_searchers(specs):

    # unpack from specs
    g = specs.graph
    capture_range = specs.capture_range
    zeta = specs.zeta
    m = specs.size_team

    if specs.start_searcher_v is None:
        # if initial position was not defined by user
        v_list = placement_list(specs, 's')
        if specs.searcher_together:
            v_list = searchers_start_together(m, v_list)
        # len(v0) = m
        v0 = v_list
        specs.set_start_searchers(v0)
    else:
        # if it was, use that
        v0 = specs.start_searcher_v

    # get graph vertices
    V, n = ext.get_set_vertices(g)
    if any(v0) not in V:
        print("Vertex out of range, V = {1, 2...n}")
        return None

    searchers = create_dict_searchers(g, v0, capture_range, zeta)

    return searchers

def create_target(specs):
    """Needs to be called AFTER belief class"""
    # TODO expand for starting target in specific position
    # for now: if wants to start target in place, set specs.start_target_v_list,
    # with the first item being the true position

    # unpack from target
    motion_rule = specs.target_motion
    g = specs.graph
    my_seed = specs.target_seed

    v_list = specs.start_target_v_list
    true_position = specs.start_target_true

    # generate motion matrix M
    motion_matrix = set_motion_matrix(g, motion_rule)

    target = MyTarget(v_list, motion_matrix, true_position, my_seed)

    return target

# instance parameters
def draw_v_random(g_or_n, q=1, my_seed=None):
    """Choose possible random vertices for the starting point of the target
    positions are given in model indexing (1,2...)
    :param g_or_n : graph
    :param q : number of possible vertices
    :param my_seed : random seed generator (optional)
    return: v_target --> possible initial vertices of the target
     v_left: target from the graph that are 'free' """
    # get set of vertices
    V, n = ext.get_set_vertices(g_or_n)

    v_target = []

    if my_seed is None:
        v_left = V

        # randomly pick vertices
        for pos in range(0, q):
            my_v = int(random.choice(v_left))
            v_target.append(my_v)

            # take out that vertex from list of possible vertices
            v_left.remove(my_v)
    else:
        v_target = pick_pseudo_random(V, my_seed, q)
        v_left = ext.get_v_left(n, v_target)

    return v_target, v_left

def pick_pseudo_random(my_list: list, my_seed: int, qty: int, replace_opt=False):

    # set seed
    np.random.seed(my_seed)

    # idx_list = np.random.randint(low=0, high=last_idx, size=qty)
    random_list = np.random.choice(my_list, qty, replace=replace_opt).tolist()

    return random_list

def placement_list(specs, op='s'):
    """Call only when the initial position is random (no list defined by user)
    Make sure belief and searchers' start vertices are far away (out of reach)
    so that b_c(0) = 0
    :param specs : inputs
    :param op : 's' to place searchers, 't' to get list of vertices for belief """

    g = specs.graph

    # placing searchers or belief?
    if op == 's':
        v_input = specs.start_target_v_list
        my_seed = specs.searcher_seed
        # check if searchers are starting together
        if specs.searcher_together:
            q = 1
        else:
            q = specs.size_team
    else:
        v_input = specs.start_searcher_v
        # quantity of possible nodes (probability > 0)
        q = specs.qty_possible_nodes
        # seed
        my_seed = specs.target_seed

    if v_input is None:
        # draw q random nodes
        v_list, v_left = draw_v_random(g, q, my_seed)
    else:
        # check if it's out of reach at t = 0
        v_taken = v_input
        v_list = []
        out_of_reach = False
        # keep drawing until is out of reach
        while out_of_reach is False:
            v_list = draw_v_random(g, q, my_seed)[0]
            specs.change_seed(my_seed, 't')
            out_of_reach = check_reachability(g, specs.capture_range, v_list, v_taken)
            my_seed += 500

    # check for bug (not a vertex in graph)
    V = ext.get_set_vertices(g)[0]
    if any(v_list) not in V:
        print("Vertex out of range, V = {1, 2...n}")
        exit()

    return v_list

# instance parameters (belief)
def set_initial_belief(g_or_n, v_list: list, type_distribution='uniform'):
    """ Make initial belief vector (assign probabilities to list of vertices) based on
    :param g_or_n : graph or n, to obtain number of total vertices (n)
    :param v_list : list of vertices # [1, 3...] to be assigned non zero probability
    :param type_distribution : type of distribution, default is uniform """

    n = ext.get_set_vertices(g_or_n)[1]

    # create belief vector of zeros
    b_0 = np.zeros(n + 1)

    # get number of possible initial vertices
    q = len(v_list)

    # if the initial_prob is not a list with pre-determined initial probabilities
    # for each possible initial vertex
    # just consider equal probability between possible vertices
    if type_distribution == 'uniform':
        prob_init = 1/q
    else:
        prob_init = type_distribution

    # initial target belief (vertices in v_list will change value, others remain zero)
    # b[0] = 0
    b_0[v_list] = prob_init

    return list(b_0)

# instance parameters (target)
def my_target_motion(g, init_vertex: list,  init_prob='uniform',  motion_rule='random'):
    """Create target motion model matrix and initial belief:
    g is the graph
    init_vertex" list with the nodes numbers (starting at 1)
    init_prob: probability for each none in init_vertex.
    If left blank probability will be equal between all vertices in init_vertex """

    b_0 = set_initial_belief(g, init_vertex, init_prob)

    # target motion model - Markovian, random
    M = set_motion_matrix(g, motion_rule)

    return b_0, M

def set_motion_matrix(g, motion_rule='random'):
    """Define Markovian motion matrix for target
    unless specified, motion is random (uniform) according to neighbor vertices"""

    V_, n = ext.get_idx_vertices(g)

    M = np.zeros((n, n))

    if motion_rule == 'random':
        for v in V_:
            delta_prime = [v] + g.vs[v]["neighbors"]
            n_nei = len(delta_prime)
            prob_v = 1/n_nei
            M[v, delta_prime] = prob_v
    elif motion_rule == 'static':
            M = np.identity(n)
    else:
        exit(print('No other option at this moment'))

    my_M = M.tolist()
    return my_M

def v_list_from_belief(b0: list):

    v_list = []
    b_v = b0[1:]

    v = 0
    for el in b_v:
        v += 1
        if el > 0.0001:
            v_list.append(v)

    return v_list

def target_restricted(g, init_possible=4, my_seed=None):

    # target
    # corners
    v_init_target = [1, 8, 71, 78]
    n_y_target = 3

    v_rtd = random_corners(v_init_target, n_y_target)

    v_target = []

    for i in range(0, init_possible):
        print(v_target)
        my_v = pick_pseudo_random(v_rtd[i], my_seed, 1)
        v_target.append(my_v[0])

    # center
    # searchers
    v_init_searcher = [34]
    n_y_searcher = 4
    v_dict = random_corners(v_init_searcher, n_y_searcher)
    v_possible = v_dict[0]

    # print(v_possible)

    return v_target, v_possible

def random_corners(v_init, n_y, n_columns=10):

    V_rtd = {}

    i = 0
    while i < len(v_init):
        # "groups"
        V_rtd[i] = []
        j = 0
        while j < n_y:
            # rows
            v_it = v_init[i] + j
            v_numbers = list(range(v_it, v_it + (n_y * n_columns), n_columns))
            for el in v_numbers:
                # append
                V_rtd[i].append(el)

            j += 1
            print(V_rtd)
        i += 1

    return V_rtd

def list_random_q(turns=200):

    list_qty = []

    for i in range(0, turns):
        random_n = np.random.randint(low=2, high=15)
        list_qty.append(random_n)

    return list_qty

# instance parameters (searchers)
def searchers_start_together(m: int, v):
    """Place all searchers are one vertex
    :param m : number of searchers
    :param v : integer or list"""

    if isinstance(v, int):
        my_v = v
    else:
        my_v = v[0]

    v_searchers = []
    for i in range(m):
        v_searchers.append(my_v)

    return v_searchers

def create_dict_searchers(g, v0: list, capture_range=0, zeta=None):
    """Create searchers (dictionary with id number as keys).
            Nested: initial position, capture matrices for each vertex"""

    # set of searchers S = {1,..m}
    S = ext.get_set_searchers(v0)[0]
    # create dict
    searchers = {}
    for s_id in S:
        v = ext.get_v0_s(v0, s_id)
        cap_s = ext.get_capture_range_s(capture_range, s_id)
        zeta_s = ext.get_zeta_s(zeta, s_id)

        # create each searcher
        s = MySearcher(s_id, v, g, cap_s, zeta_s)

        # store in dictionary
        searchers[s_id] = s

    return searchers

def check_same_vertex(v_target: list, v_searchers: list):
    """Check if the searchers and target are in different initial vertices"""
    common_el = set(v_target) - (set(v_target) - set(v_searchers))
    if len(common_el) > 0:
        print("searchers can't be in same vertex as target at t = 0")
        return False
    else:
        return True

def check_reachability(g, capture_range, v_target, v_searchers):
    """Check if the target is within reach of any one of the searchers"""

    init_is_ok = True

    for v in v_searchers:
        for vt in v_target:
            distance = ext.get_node_distance(g, vt, v)
            if distance <= capture_range:
                init_is_ok = False
                return init_is_ok

    return init_is_ok

# ----------------------------------------------------------------------------------------------------------------------
# Create graphs
def my_graph(number_vertex: int, ref: str, graph_opt=1, w=None, h=None):
    """Function to create graph according to user inputs"""
    # Graph representing environment
    # create new graph
    g = Graph(directed=False)
    g.add_vertices(number_vertex)

    if graph_opt == 1:
        g.add_edges([(0, 1), (0, 2), (1, 3), (1, 4), (2, 4), (4, 5), (5, 6)])
    elif graph_opt == 2:
        g.add_edges([(0, 1), (1, 2)])
    elif graph_opt == 3:
        g.add_edges([(0, 1), (0, 2), (1, 3), (1, 4), (2, 4), (4, 5), (5, 6), (6, 7)])
    elif graph_opt == 4:
        # create graph from Hollinger 2009 - MUSEUM
        g.add_edges([(0, 55), (1, 55), (2, 55), (3, 55),
                     (55, 56), (55, 54), (55, 34), (55, 33), (55, 35),
                     (54, 32), (54, 31), (54, 30), (54, 53),
                     (53, 29), (53, 28), (53, 27), (53, 26), (53, 39), (53, 52), (53, 38), (53, 37), (53, 36),
                     (52, 57), (52, 25), (52, 50),
                     (50, 51), (50, 49),
                     (25, 24), (24, 43),
                     (57, 56), (57, 4), (57, 5), (57, 6), (57, 40), (57, 41), (57, 43),
                     (49, 48), (48, 43), (48, 59),
                     (43, 42), (43, 44), (43, 22),
                     (44, 7), (44, 8), (44, 9), (44, 21), (44, 23), (44, 45),
                     (21, 45), (21, 22),
                     (46, 45), (45, 10), (46, 11), (46, 12), (46, 13), (46, 47), (46, 20), (47, 14),
                     (58, 47), (58, 15), (58, 59), (58, 18), (58, 19), (59, 17), (59, 16)])

    elif graph_opt == 5:
        g = add_edge_grid(g, w, h)

    elif graph_opt == 7:
        # create graph from Hollinger 2009 - OFFICE
        edges = [(1, 2), (1, 5), (1, 9), (1, 7),
                 (2, 1), (2, 3), (2, 4),
                 (3, 2),
                 (4, 2),
                 (5, 1),
                 (6, 7),
                 (7, 6), (7, 8), (7, 1),
                 (8, 7),
                 (9, 1), (9, 10),
                 (10, 9), (10, 11), (10, 12), (10, 13), (10, 18), (10, 19), (10, 20), (10, 27),
                 (11, 10), (11, 14), (11, 28),
                 (12, 10), (12, 15), (12, 16),
                 (13, 10), (13, 17),
                 (14, 11), (14, 15),
                 (15, 14), (15, 12),
                 (16, 12), (16, 17),
                 (17, 16), (17, 13),
                 (18, 10), (18, 19), (18, 26), (18, 70),
                 (19, 10), (19, 18), (19, 20), (19, 25), (19, 26),
                 (20, 10), (20, 19), (20, 21), (20, 24),
                 (21, 20), (21, 22),
                 (22, 21), (22, 23),
                 (23, 22), (23, 24),
                 (24, 23), (24, 20), (24, 25),
                 (25, 24), (25, 19), (25, 26),
                 (26, 25), (26, 19), (26, 18),
                 (27, 10), (27, 49),
                 (28, 11), (28, 29),
                 (29, 28), (29, 34),
                 (30, 31), (30, 33),
                 (31, 30), (31, 32),
                 (32, 33), (32, 31), (32, 37),
                 (33, 32), (33, 36),
                 (34, 29), (34, 35),
                 (35, 34), (35, 36), (35, 40), (35, 49),
                 (36, 35), (36, 33), (36, 39), (36, 37),
                 (37, 36), (37, 32), (37, 38),
                 (38, 37), (38, 42),
                 (39, 36), (39, 40), (39, 41),
                 (40, 35), (40, 39), (40, 43),
                 (41, 39), (41, 42),
                 (42, 41), (42, 38),
                 (43, 40), (43, 44), (43, 47),
                 (44, 43), (44, 45),
                 (45, 44), (45, 46),
                 (46, 45), (46, 47),
                 (47, 43), (47, 46), (47, 48),
                 (48, 47), (48, 50),
                 (49, 35), (49, 27), (49, 58), (49, 50),
                 (50, 48), (50, 49), (50, 51), (50, 52),
                 (51, 50),
                 (52, 50), (52, 56), (52, 53),
                 (53, 52), (53, 54), (53, 55), (53, 56),
                 (54, 53), (54, 55),
                 (55, 54), (55, 53),
                 (56, 52), (56, 53), (56, 64),
                 (57, 58), (57, 63),
                 (58, 49), (58, 57), (58, 59), (58, 62),
                 (59, 60), (59, 61), (59, 58),
                 (60, 70), (60, 61), (60, 59),
                 (61, 60), (61, 59), (61, 69), (61, 68),
                 (62, 58), (62, 67),
                 (63, 57), (63, 66), (63, 64),
                 (64, 63), (64, 56), (64, 65),
                 (65, 64), (65, 66),
                 (66, 65), (66, 67), (66, 63),
                 (67, 62), (67, 66), (67, 68),
                 (68, 61), (68, 67), (68, 69),
                 (69, 61), (69, 68),
                 (70, 18), (70, 60)]

        edges = map(lambda x: (x[0] - 1, x[1] - 1), edges)
        g.add_edges(edges)

        g.simplify(multiple=True, loops=False)

    V_, n = ext.get_idx_vertices(g)
    V = ext.get_set_vertices(g)[0]

    # label starting at 1
    g.vs["label"] = V

    # find shortest path length between any two vertices
    short_paths = g.shortest_paths_dijkstra()
    g["path_len"] = short_paths

    # find neighbors to each vertex
    for v in V_:
        nei = g.neighbors(v)
        g.vs[v]["neighbors"] = nei

    # name graph
    g["name"] = ref
    return g

def create_graph(n: int, edges: list, base_name: str):
    """Create iGraph based on
    :param n : number of vertices
    :param edges : list of connections, python index: [(0, 1), (1,2)...]
    :param base_name : """

    # create new graph
    g = Graph(directed=False)
    g.add_vertices(n)

    # make sure it's indexed right
    v1 = [e[0] for e in edges]
    v2 = [e[1] for e in edges]
    if min(v1) > 0 or min(v2) > 0:
        edges = [(e[0] - 1, e[1] - 1) for e in edges]

    g.add_edges(edges)

    V = ext.get_set_vertices(g)[0]

    # label starting at 1
    g.vs["label"] = V

    # find shortest path length between any two vertices
    short_paths = g.shortest_paths_dijkstra()
    g["path_len"] = short_paths

    # find neighbors to each vertex
    for v in V:
        vidx = ext.get_python_idx(v)
        nei = g.neighbors(vidx)
        g.vs[vidx]["neighbors"] = nei

    # name graph
    g["name"] = base_name
    return g

def create_grid_graph(w: int, h: int):
    """wxh graph"""
    n_vertex = w*h
    ref = "G" + str(n_vertex) + 'V' + '_grid'
    graph_opt = 5
    g = my_graph(n_vertex, ref, graph_opt, w, h)
    save_graph(g, ref)
    plot_simple_graph(g)
    save_pdf(g, ref)

    return g

def add_edge_grid(g, w, h):

    n = len(g.vs)
    vertices = list(range(0, n))

    for vertex in vertices:
        # not the first column
        if vertex % w != 0:
            g.add_edge(vertex, vertex - 1)
        # not the last row
        if vertex < w*(h - 1):
            g.add_edge(vertex, vertex + w)

    return g

def save_and_plot_graph(g, ref: str, v_searchers: list, v_target: list):
    save_graph(g, ref)
    plot_start_config(g, ref, v_searchers, v_target)
    save_pdf(g, ref)

def save_graph(g, ref):
    """Save created graph as pickle file"""
    graph_name = str(ref) + ".p"
    file_path = ext.get_whole_path(graph_name, 'graphs')
    pickle.dump(g, open(file_path, "wb"))

def plot_start_config(g, ref, v_searchers: list, v_target: list):
    """Plot my graph and show me where the searchers started"""

    g.vs["color"] = "green"
    start_searcher = ext.get_python_idx(v_searchers)
    start_target = ext.get_python_idx(v_target)
    for st in start_searcher:
        g.vs[st]["color"] = "blue"
    for st in start_target:
        g.vs[st]["color"] = "red"
    layout = g.layout("kk")
    name_file = "G" + ref + ".pdf"
    plot(g, name_file, layout=layout, bbox=(300, 300), margin=20)

def save_pdf(g, ref):
    """Plot my graph and show me where the searchers started"""

    g.vs["color"] = "gray"
    if 'grid' in ref:
        my_layout = g.layout("grid")
    else:
        my_layout = g.layout("kk")

    name_file = ref + ".pdf"
    file_path = ext.get_whole_path(name_file, 'graphs')

    plot(g, file_path, layout=my_layout, bbox=(300, 300), margin=20)

def plot_simple_graph(g, my_layout='grid'):
    graph_file = g["name"] + ".pdf"
    name_file = ext.get_whole_path(graph_file, 'graphs')
    g.vs["color"] = "gray"
    plot(g, name_file, layout=my_layout)

def create_office_graph():
    """MESPP paper, Fig 2, first graph (MUSEUM)"""
    n_vertex = 60
    ref = 'G' + str(n_vertex) + 'V'
    graph_opt = 4
    g = my_graph(n_vertex, ref, graph_opt)
    save_graph(g, ref)

    return g

def create_museum_graph():
    """MESPP paper, Fig 2, second graph (OFFICE)"""
    n_vertex = 70
    ref = 'G' + str(n_vertex) + 'V'
    graph_opt = 7
    g = my_graph(n_vertex, ref, graph_opt)
    save_graph(g, ref)
    plot_simple_graph(g, 'kk')

    return g

def create_school_graph(edges):
    """School graph based on list of connections"""

    n = 46
    base_name = 'G' + str(n) + 'V_ss2'
    g = create_graph(n, edges, base_name)
    save_graph(g, base_name)
    plot_simple_graph(g, 'kk')

    return g

def create_test_graph():
    n_vertex = 7
    ref = 'G' + str(n_vertex) + 'V' + '_test'
    graph_opt = 1
    g = my_graph(n_vertex, ref, graph_opt)
    save_graph(g, ref)
    save_pdf(g, ref)

if __name__ == "__main__":

    # create graphs 1-5, 7
    create_test_graph()
    create_office_graph()
    create_grid_graph(10, 10)
    create_grid_graph(16, 16)
    create_grid_graph(8, 8)
    create_grid_graph(3, 3)
    create_museum_graph()


text:
# Authors Royce Burningham & Carter Green

import Class.Adf24Tx2Rx4 as Adf24Tx2Rx4
import Class.RadarProc as RadarProc
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import test_stepper as step
import RPi.GPIO as GPIO
import sar
import datetime as datetime
import sys

STEP_DISTANCE = 3.4636e-3
c0 = 3e8

# Setup Connection
Brd = Adf24Tx2Rx4.Adf24Tx2Rx4()
Brd.BrdRst()

# Configure Receiver
Brd.RfRxEna()
TxPwr = 100

# Configure Transmitter (Antenna 0 - 3, Pwr 0 - 31)
Brd.RfTxEna(1, TxPwr)

# Configure Up-Chirp
dCfg = {
    "fs": 1.0e6,
    "fStrt": 23.9e9,
    "fStop": 24.3e9,
    "TRampUp": 260 / 1.0e6,
    "Tp": 300 / 1.0e6,
    "N": 256,
    "StrtIdx": 0,
    "StopIdx": 1,
    "MimoEna": 0,
    "Frms": 1,
}

Brd.RfMeas("Adi", dCfg)

kf = (dCfg["fStop"] - dCfg["fStrt"]) / dCfg["TRampUp"]
fs = Brd.RfGet("fs")
FuSca = Brd.Get("FuSca")

dRpCfg = {
    "RemoveMean": 1,
    "FFT": 2 ** 9,
    "FuSca": FuSca,
    "fs": fs,
    "kf": kf,
    "RMin": 0,
    "RMax": 4,
    "dB": 0,
    "Abs": 0,
    "Ext": 1,
    "NrFrms": 128,
}

# Configure Range Profile
Proc = RadarProc.RadarProc()
Proc.CfgRangeProfile(dRpCfg)

fStrt = Brd.Adf_Pll.fStrt
fStop = Brd.Adf_Pll.fStop
TRampUp = Brd.Adf_Pll.TRampUp
n = np.arange(int(dCfg["N"]))

Range = Proc.GetRangeProfile("Range")

# Configure SAR Features
if len(sys.argv) < 2:
    points = 166
    sum_size = 9
else:
    distance = float(sys.argv[1])
    points = int(distance / STEP_DISTANCE)
    if int(sys.argv[2]) % 2 == 0:
        sum_size = int(sys.argv[2]) + 1
    else:
        sum_size = int(sys.argv[2])

DataCube = np.zeros((dCfg["N"], points), dtype="int16")

# Collect Data
for i in range(0, points):
    Data = Brd.BrdGetData()
    DataCube[:, i] = np.flip(Data[:, 0])
    step.step(8)
    # print(i)

# Process Data
Out = sar.sar(DataCube, sum_size, window=False) / FuSca
Out2 = sar.focused_sar(DataCube, sum_size, kf, window=False) / FuSca

X, Y = np.meshgrid(np.arange(points - sum_size), Range)
lam = c0 / (dCfg["fStrt"] + dCfg["fStop"]) / 2
CR = X * lam * Y / (2 * STEP_DISTANCE * sum_size)

# Plot SAR Image
fig, ax1 = plt.subplots(1, 2, figsize=(8.5, 8.5))

beep = ax1[0].pcolormesh(
    X * lam / (2 * STEP_DISTANCE * sum_size),
    Y,
    np.abs(Out)[0 : len(Range), 0 : points - sum_size],
    cmap=cm.Greys,
)
fig.colorbar(beep, ax=ax1[1])
ax1[0].set_xlabel("Cross-Range (m)")
ax1[0].set_ylabel("Down-Range (m)")
ax1[0].set_title("Unfocused Sar")

ax1[1].pcolormesh(
    X * lam / (2 * STEP_DISTANCE * sum_size),
    Y,
    np.abs(Out2)[0 : len(Range), 0 : points - sum_size],
    cmap=cm.Greys,
)
ax1[1].set_xlabel("Cross-Range (m)")
ax1[1].set_ylabel("Down-Range (m)")
ax1[1].set_title("Focused Sar")

fig.suptitle("Distance: {} Array Size: {}".format(sys.argv[1], sys.argv[2]))
fig.tight_layout()

# Save Plot
fig_path = "/home/pi/Desktop/DemoRadSoftware/Software/Python/Figures/"
today = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
plt.savefig(fig_path + today + ".png", dpi=300)

# Save Data
my_path = "/home/pi/Desktop/DemoRadSoftware/Software/Python/SavedData/"
np.save(my_path + today + ".npy", DataCube)

# Set GPIO ports to 0
GPIO.cleanup()

# Let user know image and data are saved
print(
    "\nDone! at "
    + today
    + "\n Distance: {} Array Size: {}".format(sys.argv[1], sys.argv[2])
)

ext: py

max_stars_repo_name: carterjgreen/undergrad-thesis

max_issues_repo_name: carterjgreen/undergrad-thesis

max_forks_repo_name: carterjgreen/undergrad-thesis

reason: import numpy


text:
import matplotlib.pyplot as plt
import analytics.forecast.forecast as fc
import analytics.location.path as ap
import analytics.definitions as adef
import pvlib as pv
import pandas as pd
import numpy as np
import pandas.plotting
import pytz
import mplcursors
from loguru import logger
from tqdm import tqdm
import datetime
from typing import Optional, Union
from pint import UnitRegistry

def main():
    # tomorrow = datetime.date.today() + datetime.timedelta(days=1)
    # time = datetime.datetime.combine(tomorrow, time=datetime.time(hour=16),
    #                                  tzinfo=pytz.utc)
    # plot_forecast_area(top_right=pv.location.Location(41.5, -139.5),
    #                    bottom_left=pv.location.Location(40.5, -140.5),
    #                    time=time, metric="low_clouds")

    points = [
        pv.location.Location(36, 141),
        pv.location.Location(42, 162),
        pv.location.Location(44, 168),
        # pv.location.Location(44, 175),
        # pv.location.Location(44, -175),
        # pv.location.Location(44, -165),
        # pv.location.Location(43, -155),
        # pv.location.Location(40, -145),
        # pv.location.Location(37, -135),
        # pv.location.Location(35, -128),
        # pv.location.Location(32, -118)
    ]
    start = datetime.datetime(day=1, month=1, year=2020, tzinfo=pytz.utc)
    end = start + datetime.timedelta(days=1)
    delta = (end - start) / len(points)
    timestamps = [start + (delta * i) for i in range(len(points))]
    path: Optional[ap.SegmentedPath] = None
    for p, t in zip(points, timestamps):
        if path is None:
            path = ap.SegmentedPath.create(p, t)
        else:
            path.append_point(p, t)
    plot_forecast_path(path, "ghi_raw")

    # ureg = UnitRegistry()
    # eff = (50/3.2) * ureg.watt / (ureg.feet ** 2)
    # eff.ito(ureg.watt / (ureg.meters ** 2))
    # print(eff)
    # area = 68000 * ureg.feet ** 2  # square feet
    # area.ito(ureg.meter ** 2)
    # i = calc_integration("forecast-2020-01-19T21-39-46.csv", "ghi_raw")
    # hours = i.astype('timedelta64[h]')
    # hours = hours.astype('int')
    # epm2 = hours * ureg.hour * ureg.watt / (ureg.meter ** 2)  # energy per m^2
    # p = epm2 * area
    # print(p)
    # print(p.to(ureg.kilowatt * ureg.hour))

def plot_forecast_path(path: ap.Path, metric: str):
    pandas.plotting.register_matplotlib_converters()
    df = fc.get_forecast_path(path)

    timestamp = datetime.datetime.now().isoformat()
    timestamp = timestamp[:-7].replace(":", "-")
    with open(adef.get_output_path(f"forecast-{timestamp}.csv"), mode="w") as f:
        df.to_csv(path_or_buf=f)

    plt.plot(df[metric])
    plt.show()

def calc_integration(path: str, column: str) -> Union[float, np.array]:
    df = pd.read_csv(path, index_col=0, parse_dates=True)
    sel_df = df[column].fillna(value=0)
    integral = np.trapz(sel_df, x=df.index, axis=0)
    return integral

def plot_forecast_area(top_right: pv.location.Location, bottom_left: pv.location.Location,
                       time: datetime, metric: str,
                       imshow_kwargs: dict = None, resolution: int = 16, use_tqdm=True):
    """
    :param top_right: Location with the greatest latitude and longitude
    :param bottom_left: Location with the least latitude and longitude
    :param time: Time of forecast
    :param metric: The metric to forecast; any of pvlib's forecast metrics
    :param imshow_kwargs: kwargs to pass to imshow
    :param resolution: number of discrete points to break up each degree into
    :param use_tqdm: boolean, whether to show tqdm loading bar or not
    :return: None
    Valid metrics:
    mid_clouds, convect_clouds, boundary_clouds, ghi_raw, high_clouds, total_clouds
    low_clouds, temp_air, wind_speed_gust, wind_speed_u, wind_speed_v
    """
    pandas.plotting.register_matplotlib_converters()

    # Default key-word arguments (kwargs)
    if imshow_kwargs is None:
        imshow_kwargs = dict()  # use empty dictionary

    min_lat = bottom_left.latitude
    max_lat = top_right.latitude
    min_lon = bottom_left.longitude
    max_lon = top_right.longitude
    n_lat = int((max_lat - min_lat) * resolution)
    n_lon = int((max_lon - min_lon) * resolution)

    lats = np.linspace(min_lat, max_lat, n_lat)
    lons = np.linspace(min_lon, max_lon, n_lon)

    locations = [pv.location.Location(lat, lon) for lat in lats for lon in lons]
    results = []

    it = tqdm if use_tqdm else iter

    for loc in it(locations):
        forecast = fc.get_forecast_single(time, loc)
        res = float(forecast[metric])
        results.append(res)

    coverage = np.array(results)
    grid = coverage.reshape((n_lat, n_lon))

    use_imshow_kwargs = dict(extent=(min_lon, max_lon, min_lat, max_lat),
                             origin="lower",  # plot starting at the lower left corner
                             interpolation='nearest', cmap="Blues_r", vmin=0, vmax=100)
    use_imshow_kwargs.update(imshow_kwargs)
    im = plt.imshow(grid, **use_imshow_kwargs)
    ax = plt.gca()
    fig = plt.gcf()
    cbar = fig.colorbar(im, ax=ax)
    cbar.minorticks_on()
    ax.set_xticks(np.arange(min_lon, max_lon, 0.25))
    ax.set_yticks(np.arange(min_lat, max_lat, 0.25))
    ax.set_title(f"{metric} {time}")
    plt.show()

if __name__ == "__main__":
    main()

max_stars_repo_name: panzer/capstone-mppt

max_issues_repo_name: panzer/capstone-mppt

max_forks_repo_name: panzer/capstone-mppt


text:
import random
import numpy as np

class BatAlgorithm():
    def __init__(self, D, NP, N_Gen, A, r, Qmin, Qmax, Lower, Upper, function):
        self.D = D  #dimension
        self.NP = NP  #population size
        self.N_Gen = N_Gen  #generations
        self.A = A  #loudness
        self.r = r  #pulse rate
        self.Qmin = Qmin  #frequency min
        self.Qmax = Qmax  #frequency max
        self.Lower = Lower  #lower bound
        self.Upper = Upper  #upper bound

        self.f_min = 0.0  #minimum fitness

        self.Lb = [0] * self.D  #lower bound
        self.Ub = [0] * self.D  #upper bound
        self.Q = [0] * self.NP  #frequency

        self.v = [[0 for i in range(self.D)] for j in range(self.NP)]  #velocity
        self.Sol = [[0 for i in range(self.D)] for j in range(self.NP)]  #population of solutions
        self.Fitness = [0] * self.NP  #fitness
        self.best = [0] * self.D  #best solution
        self.Fun = function

    def best_bat(self):
        i = 0
        j = 0
        for i in range(self.NP):
            if self.Fitness[i] < self.Fitness[j]:
                j = i
        for i in range(self.D):
            self.best[i] = self.Sol[j][i]
        self.f_min = self.Fitness[j]

    def init_bat(self):
        for i in range(self.D):
            self.Lb[i] = self.Lower
            self.Ub[i] = self.Upper

        for i in range(self.NP):
            self.Q[i] = 0
            for j in range(self.D):
                rnd = np.random.uniform(0, 1)
                self.v[i][j] = 0.0
                self.Sol[i][j] = self.Lb[j] + (self.Ub[j] - self.Lb[j]) * rnd
            self.Fitness[i] = self.Fun(self.D, self.Sol[i])
        self.best_bat()

    def simplebounds(self, val, lower, upper):
        if val < lower:
            val = lower
        if val > upper:
            val = upper
        return val

    def move_bat(self):
        S = [[0.0 for i in range(self.D)] for j in range(self.NP)]

        self.init_bat()

        for t in range(self.N_Gen):
            for i in range(self.NP):
                rnd = np.random.uniform(0, 1)
                self.Q[i] = self.Qmin + (self.Qmax - self.Qmin) * rnd
                for j in range(self.D):
                    self.v[i][j] = self.v[i][j] + (self.Sol[i][j] -
                                                   self.best[j]) * self.Q[i]
                    S[i][j] = self.Sol[i][j] + self.v[i][j]

                    S[i][j] = self.simplebounds(S[i][j], self.Lb[j],
                                                self.Ub[j])

                rnd = np.random.random_sample()

                if rnd > self.r:
                    for j in range(self.D):
                        S[i][j] = self.best[j] + 0.001 * random.gauss(0, 1)
                        S[i][j] = self.simplebounds(S[i][j], self.Lb[j],
                                                self.Ub[j])

                Fnew = self.Fun(self.D, S[i])

                rnd = np.random.random_sample()

                if (Fnew <= self.Fitness[i]) and (rnd < self.A):
                    for j in range(self.D):
                        self.Sol[i][j] = S[i][j]
                    self.Fitness[i] = Fnew

                if Fnew <= self.f_min:
                    for j in range(self.D):
                        self.best[j] = S[i][j]
                    self.f_min = Fnew

        print(self.f_min)

max_stars_repo_name: buma/BatAlgorithm

max_issues_repo_name: buma/BatAlgorithm

max_forks_repo_name: buma/BatAlgorithm


text:
import netCDF4 as nc
import numpy as np
import pandas as pd

from smrf.data.gridded_input import GriddedInput
from smrf.utils.utils import apply_utm

def metadata_name_from_index(index):
    return 'grid_y{:d}_x{:d}'.format(index[0], index[1])

class InputNetcdf(GriddedInput):

    DATA_TYPE = 'netcdf'

    VARIABLES = [
        'air_temp',
        'vapor_pressure',
        'precip',
        'wind_speed',
        'wind_direction',
        'cloud_factor'
    ]

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

        self.variables = list(self.config.keys())
        self.variables.remove('data_type')
        self.variables.remove('netcdf_file')

    def load(self):
        """
        Load the data from a generic netcdf file
        """

        self._logger.info('Reading data coming from netcdf: {}'.format(
            self.config['netcdf_file'])
        )

        f = nc.Dataset(self.config['netcdf_file'], 'r')

        mlat = f.variables['lat'][:]
        mlon = f.variables['lon'][:]
        mhgt = f.variables['elev'][:]

        if mlat.ndim != 2 & mlon.ndim != 2:
            [mlon, mlat] = np.meshgrid(mlon, mlat)

        # get the values that are in the modeling domain
        ind = (mlat >= self.bbox[1]) & \
            (mlat <= self.bbox[3]) & \
            (mlon >= self.bbox[0]) & \
            (mlon <= self.bbox[2])

        mlat = mlat[ind]
        mlon = mlon[ind]
        mhgt = mhgt[ind]

        # create some fake station names based on the index
        a = np.argwhere(ind)
        primary_id = [metadata_name_from_index(i) for i in a]
        self._logger.debug('{} grid cells within model domain'.format(len(a)))

        # create a metadata dataframe to store all the grid info
        metadata = pd.DataFrame(index=primary_id,
                                columns=('X', 'Y', 'latitude',
                                         'longitude', 'elevation'))

        metadata['latitude'] = mlat.flatten()
        metadata['longitude'] = mlon.flatten()
        metadata['elevation'] = mhgt.flatten()
        metadata = metadata.apply(apply_utm,
                                  args=(self.topo.zone_number,),
                                  axis=1)

        self.metadata = metadata

        t = f.variables['time']
        time = nc.num2date(t[:].astype(int), t.getncattr(
            'units'), t.getncattr('calendar'))
        # Drop milliseconds and prepare to use as pandas DataFrame index
        time = pd.DatetimeIndex(
            [str(tm.replace(microsecond=0)) for tm in time], tz=self.time_zone
        )

        for variable in self.VARIABLES:
            v_file = self.config[variable]
            self._logger.debug(
                'Loading variable {} from netcdf field {}'.format(
                    variable, v_file)
            )

            df = pd.DataFrame(index=time, columns=primary_id)
            for i in a:
                g = metadata_name_from_index(i)
                df[g] = f.variables[v_file][:, i[0], i[1]]

            # deal with any fillValues
            try:
                fv = f.variables[v_file].getncattr('_FillValue')
                df.replace(fv, np.nan, inplace=True)
            except Exception:
                pass

            # Set variable and subset by start and end time
            setattr(self, variable, df)

        f.close()

max_stars_repo_name: UofU-Cryosphere/smrf

max_issues_repo_name: micah-prime/smrf

max_forks_repo_name: UofU-Cryosphere/smrf


text:
import os
import math
import subprocess
import random
import numpy

f = open("standard_model.church")
ch_model = [l.rstrip() for l in f.readlines()]
f.close()

# get priors from slider experiment (23)
#pfile = open("smoothed_15marbles_priors_withnames.txt")
pfile = open("/Users/titlis/cogsci/projects/stanford/projects/thegricean_sinking-marbles/experiments/12_sinking-marbles-prior15/results/data/smoothed_15marbles_priors_epsilon.txt")
priors_15 = [l.rstrip() for l in pfile.readlines()[1:]]
pfile.close()

wpriors = [.3, .5, .6, .7, .8, .9]

rfile = open("modelresults/results_rsa_epsilonsmooth.txt.tmp.results","w")

pline = ch_model.index("(define theprior (list .1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1))")
iline = ch_model.index("(define item 'blabla)")

for k,p in enumerate(priors_15):
	psplit = p.split("\t")
	item = psplit[0]
	print item
	prior = " ".join(psplit[1:])
	print prior
	for ww in wpriors:
		ch_model[pline] = "(define theprior (list %s))"%prior
		ch_model[iline] = "(define item '(%s))"%item
		ofile = open("m.church","w")
		ofile.write("\n".join(ch_model))
		ofile.close()
		subresults = subprocess.Popen(['church',ofile.name], stdout=subprocess.PIPE)
		subresults.communicate()[0]
		rfile.write("".join(open("modelresults/raw_rsa_epsilonsmooth_results.txt.tmp.results").readlines()))

rfile.close()

max_stars_repo_name: thegricean/sinking-marbles

max_issues_repo_name: thegricean/sinking-marbles

max_forks_repo_name: thegricean/sinking-marbles


text:
"""
Selection Detail Example
========================
This example shows a selection that links two views of data: the left panel
contains one point per object, and the right panel contains one line per
object. Clicking on either the points or lines will select the corresponding
objects in both views of the data.

The challenge lies in expressing such hierarchical data in a way that Altair
can handle. We do this by merging the data into a "long form" dataframe, and
aggregating identical metadata for the final plot.
"""
# category: interactive charts
import altair as alt
import pandas as pd
import numpy as np

np.random.seed(0)

n_objects = 20
n_times = 50

# Create one (x, y) pair of metadata per object
locations = pd.DataFrame({
    'id': range(n_objects),
    'x': np.random.randn(n_objects),
    'y': np.random.randn(n_objects)
})

# Create a 50-element time-series for each object
timeseries = pd.DataFrame(np.random.randn(n_times, n_objects).cumsum(0),
                          columns=locations['id'],
                          index=pd.RangeIndex(0, n_times, name='time'))

# Melt the wide-form timeseries into a long-form view
timeseries = timeseries.reset_index().melt('time')

# Merge the (x, y) metadata into the long-form view
data = pd.merge(timeseries, locations, on='id')

# Data is prepared, now make a chart

selector = alt.selection_single(empty='all', fields=['id'])

points = alt.Chart(data).mark_point(filled=True, size=200).encode(
    x='mean(x)',
    y='mean(y)',
    color=alt.condition(selector, 'id:O', alt.value('lightgray'), legend=None),
).properties(
    selection=selector,
    width=250, height=250
).interactive()

timeseries = alt.Chart(data).mark_line().encode(
    x='time',
    y=alt.Y('value', scale=alt.Scale(domain=(-15, 15))),
    color=alt.Color('id:O', legend=None)
).transform_filter(
    selector
).properties(
    selection=selector,
    width=250, height=250
)

points | timeseries

max_stars_repo_name: dhuppenkothen/altair

max_issues_repo_name: dhuppenkothen/altair

max_forks_repo_name: dhuppenkothen/altair


text:
"""A variety of utility functions, especially for visualising the results."""

import os
import numpy as np
import cv2
import matplotlib.pyplot as plt

def import_image(image_path):
    """Import an image from the supplied path.
    Returns the image name and the image.
    """
    img = cv2.imread(image_path)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    #img_name = image_path.split('/')[-1].split('.')[0]
    return img

def plot_image(img, title=''):
    """Plot the supplied image and the corresponding title."""
    if len(img.shape) == 2:
        plt.imshow(img, cmap='gray')
    else:
        plt.imshow(img)
    plt.axis('off')
    plt.title(title, fontsize=20)

def compare_images(img_org, img_undist, titles=None):
    """Display an image comparison in a subplot."""
    if titles is None:
        titles = ('Image Before', 'Image After')

    plt.subplots(1, 2, figsize=(10, 5), dpi=80)
    plt.subplot(1, 2, 1)
    plot_image(img_org, titles[0])
    plt.subplot(1, 2, 2)
    plot_image(img_undist, titles[1])

def overlay_image(img, overlay_img, opacity=1.0):
    """Reliably combine two images based on the opacity. Black pixels are
    treated as transparent. The method also takes in a grayscale base image.
    """
    if len(img.shape) < 3:
        img = np.stack((img, img, img), axis=-1)
    if len(overlay_img.shape) < 3:
        overlay_img = np.stack((overlay_img, overlay_img, overlay_img), axis=-1)

    img_out = np.zeros(img.shape)
    overlay_img = (overlay_img * opacity).astype(int)
    faded_img = (img * (1 - opacity)).astype(int)

    for i in range(0, img.shape[2]):
        channel_out = np.where(overlay_img[:, :, i] == 0,
                               img[:, :, i],
                               faded_img[:, :, i] + overlay_img[:, :, i])

        img_out[:, :, i] = np.where(channel_out > 255, 255, channel_out)

    return img_out.astype('uint8')

def quick_rectangle(img, corners, color='green', opacity=0.9,
                    thickness=4, filled=False):
    """Draws a rectangle on the input image."""
    colors = {'green': (30, 255, 120),
              'blue': (20, 104, 229),
              'red': (224, 52, 0),
              'orange': (252, 163, 9),
              'yellow': (252, 228, 10)}

    if color.lower() not in colors:
        color = 'green'
        print('Warning unknown color, using green instead.  Please choose from:\
              green, blue, red, orange and yellow.')
    else:
        color = color.lower()

    thickness = int(thickness)
    outline = np.zeros(img.shape)
    corners = tuple(map(tuple, corners))
    outline = cv2.rectangle(outline, corners[0], corners[1], colors[color],
                            thickness=thickness)

    if filled:
        fill = np.zeros(img.shape)
        fill = cv2.rectangle(fill, corners[0], corners[1], colors[color],
                             thickness=-1)
        img = overlay_image(img, fill, opacity=opacity*0.3)

    return overlay_image(img, outline, opacity=opacity)

def save_image(image, output_folder='./output_images/',
               name='Current_Image', colorspace='RGB'):
    """Save image to a file."""
    save_path = os.path.join(output_folder, name) + '.jpg'

    if colorspace != 'BGR':
        converter = getattr(cv2, "COLOR_" + colorspace + "2BGR")
        image = cv2.cvtColor(image, converter)

    cv2.imwrite(save_path, image)

max_stars_repo_name: MarkDekker/CarND-Vehicle-Detection

max_issues_repo_name: MarkDekker/CarND-Vehicle-Detection

max_forks_repo_name: MarkDekker/CarND-Vehicle-Detection


text:
"""
Tests simple_setup_problem
"""
import numpy as np

from pymgrit.core.vector import Vector
from pymgrit.core.application import Application
from pymgrit.core.simple_setup_problem import simple_setup_problem

class VectorSimple(Vector):
    def __init__(self):
        super(VectorSimple, self).__init__()

    def __add__(self, other):
        pass

    def __sub__(self, other):
        pass

    def __mul__(self, other):
        pass

    def norm(self):
        pass

    def clone(self):
        pass

    def clone_zero(self):
        pass

    def clone_rand(self):
        pass

    def set_values(self, values):
        pass

    def get_values(self):
        pass

    def pack(self):
        pass

    def unpack(self):
        pass

class ApplicationSimple(Application):
    def __init__(self, *args, **kwargs):
        super(ApplicationSimple, self).__init__(*args, **kwargs)
        self.vector_template = VectorSimple()  # Setting the class which is used for each time point
        self.vector_t_start = VectorSimple()

    def step(self, u_start: VectorSimple, t_start: float, t_stop: float) -> VectorSimple:
        return 1

def test_simple_setup_problem():
    """
    Test simple_setup_problem
    """
    app = ApplicationSimple(t_start=0, t_stop=1, nt=101)

    problem = simple_setup_problem(problem=app, level=3, coarsening=2)

    np.testing.assert_equal(True, isinstance(problem[0], ApplicationSimple))
    np.testing.assert_equal(True, isinstance(problem[1], ApplicationSimple))
    np.testing.assert_equal(True, isinstance(problem[2], ApplicationSimple))

    np.testing.assert_equal(problem[0].nt, 101)
    np.testing.assert_equal(problem[0].t_start, 0)
    np.testing.assert_equal(problem[0].t_end, 1)
    np.testing.assert_equal(problem[0].t, np.linspace(0, 1, 101))

    np.testing.assert_equal(problem[1].nt, 51)
    np.testing.assert_equal(problem[1].t_start, 0)
    np.testing.assert_equal(problem[1].t_end, 1)
    np.testing.assert_equal(problem[1].t, np.linspace(0, 1, 101)[::2])

    np.testing.assert_equal(problem[2].nt, 26)
    np.testing.assert_equal(problem[2].t_start, 0)
    np.testing.assert_equal(problem[2].t_end, 1)
    np.testing.assert_equal(problem[2].t, np.linspace(0, 1, 101)[::4])

def test_simple_setup_problem_warning():
    """
    Test simple_setup_problem
    """
    app = ApplicationSimple(t_start=0, t_stop=1, nt=2)
    problem = simple_setup_problem(problem=app, level=2, coarsening=2)

    np.testing.assert_equal(True, isinstance(problem[0], ApplicationSimple))
    np.testing.assert_equal(True, isinstance(problem[1], ApplicationSimple))

    np.testing.assert_equal(problem[0].nt, 2)
    np.testing.assert_equal(problem[0].t_start, 0)
    np.testing.assert_equal(problem[0].t_end, 1)
    np.testing.assert_equal(problem[0].t, np.linspace(0, 1, 2))

    np.testing.assert_equal(problem[1].nt, 1)
    np.testing.assert_equal(problem[1].t_start, 0)
    np.testing.assert_equal(problem[1].t_end, 0)
    np.testing.assert_equal(problem[1].t, np.linspace(0, 1, 2)[::2])

max_stars_repo_name: masumis/pymgrit

max_issues_repo_name: masumis/pymgrit

max_forks_repo_name: masumis/pymgrit


text:
from mpi4py import MPI
import numpy as np

def psum(a):
    locsum = np.sum(a)
    rcvBuf = np.array(0.0, "d")
    MPI.COMM_WORLD.Allreduce(
        [locsum, MPI.DOUBLE], [rcvBuf, MPI.DOUBLE], op=MPI.SUM
    )
    return rcvBuf

max_stars_repo_name: hofmannedv/jupyter-tutorial

max_issues_repo_name: hofmannedv/jupyter-tutorial

max_forks_repo_name: hofmannedv/jupyter-tutorial


text:
# Lib
import argparse
from argparse import RawTextHelpFormatter
import logging
from pathlib import Path
import sys
import numpy as np
import pandas as pd
# App
from .probes.filters import list_problem_probes, exclude_probes, exclude_sex_control_probes
from .samples.postprocessQC import mean_beta_plot, beta_density_plot, cumulative_sum_beta_distribution, beta_mds_plot
from .reports import controls_report, ReportPDF

LOGGER = logging.getLogger(__name__)

class DefaultParser(argparse.ArgumentParser):
    def error(self, message):
        self._print_message('[Error]:\n')
        self._print_message(f'{message}\n\n')
        self.print_help()
        self.exit(status=2)

def detect_array(df, returns='name', on_lambda=False):
    """Determines array type using number of probes columns in df. Returns array string.
    Note: this is different from methylprep.models.arrays.ArrayType.from_probe_count, which looks at idat files.

    returns (name | filepath)
        default is 'name' -- returns a string
        if 'filepath', this also returns the filepath to the array, using ArrayType and
        methylprep.files.manifests ARRAY_TYPE_MANIFEST_FILENAMES.
    on_lambda (True | False)
        looks for manifest files in /tmp instead of ~/.methylprep_manifest_files

    returns one of: {27k, 450k, epic, epic+, mouse}
    """

    if returns == 'filepath':
        # get manifest data from .methylprep_manifest_files
        try:
            from methylprep.files.manifests import MANIFEST_DIR_PATH, MANIFEST_DIR_PATH_LAMBDA, ARRAY_TYPE_MANIFEST_FILENAMES, ARRAY_FILENAME
            from methylprep.models.arrays import ArrayType
        except ImportError:
            raise ImportError("this function requires `methylprep` be installed (to read manifest array files).")

        def get_filename(array_name):
            if on_lambda:
                man_path = Path(MANIFEST_DIR_PATH_LAMBDA).expanduser()
            else:
                man_path = Path(MANIFEST_DIR_PATH).expanduser()
            man_filename = ARRAY_FILENAME[array_name]
            man_filepath = Path(man_path, man_filename)
            return man_filepath

    # shape: should be wide, with more columns than rows. The larger dimension is the probe count.
    if df.shape[0] > df.shape[1]:
        # WARNING: this will need to be transposed later.
        col_count = (df.shape[0]) #does the index count in shape? assuming it doesn't.
    else:
        col_count = (df.shape[1])
    if 26000 <= col_count <= 28000:
        return '27k' if returns == 'name' else (ArrayType('27k'), get_filename('27k'))
    elif 440000 <= col_count <= 490000: # actual: 485512
        return '450k' if returns == 'name' else (ArrayType('450k'), get_filename('450k'))
    elif 868001 <= col_count <= 869335: # actual: 868578
        return 'epic+' if returns == 'name' else (ArrayType('epic+'), get_filename('epic+'))
    elif 860000 <= col_count <= 868000: # actual: 865860
        return 'epic' if returns == 'name' else (ArrayType('epic'), get_filename('epic'))
    elif 200000 <= col_count <= 300000: # actual: 236685 C20, 262812 v2, 284762 V3-V4, MM285 293220 --- testing, changed max from 270k to 290k for v1.4.6
        return 'mouse' if returns == 'name' else (ArrayType('mouse'), get_filename('mouse'))
    else:
        raise ValueError(f'Unsupported Illumina array type. Your data file contains {col_count} rows for probes.')

def build_parser():
    parser = DefaultParser(
        prog='methylcheck',
        description="""Visualization and quality control functions for methylation, using Illumina IDAT files processed with methylprep or sesame.

Usage
-----
To run the full battery of tests:
    `python -m methylcheck qc -d <data_file> -p all --exclude_all`
or for a standard QC REPORT
    `python -m methylcheck controls -d <data_dir>`

Commands
--------
qc
    - exclude_sex -- remove sex-chromosome-linked probes from a batch of samples.
    - exclude_control -- remove Illumina control probes  from a batch of samples.
    - exclude_probes -- remove less reliable probes, based on lists of problem probes from literature.
    look at the function help for details on how you fine tune the exclusion list.
    - plot -- various kinds of plots to reveal sample quality.
    - array_type -- QC requires the type of array.
    - verbose

controls
    - data_dir -- folder where control_probes.pkl and other methylprep outputs are found. [required]
    - outfilepath=None
        optional: specify a different folder to save XLSX file output
    - bg_offset=3000
        Illumina's defaults add 3000 to control probe intensities. You can adjust that here.
    - colorblind=False
        Set to True if you'd prefer an alternative color-coding to the
        red/yellow/green traffic light for pass/fail in XLSX file.
    - cutoff_adjust=1.0
        A value of 1.0 will mirror Illumina's recommended pass/fail thresholds. Setting it higher
        results in tests that are MORE conservative. For example, setting it to 2 or 3
        will increase the criteria 2X or 3X, respectively. Fractional values are OK.
        """,
        formatter_class=RawTextHelpFormatter,
    )

    #parser.add_argument(
    #    '-v', '--verbose',
    #    help='Provide more detailed processing messages.',
    #    action='store_true',
    #    default=False,
    #)

    parser.add_argument(
        '-d', '--debug',
        help='Provide VERY detailed processing messages.',
        action='store_true',
        default=False,
    )

    subparsers = parser.add_subparsers(dest='command', required=True)
    #subparsers.required = True # this is a python3.4-3.7 bug; cannot specify in the call above.

    qc_parser = subparsers.add_parser('qc', help="Apply filters to exclude probes based on academic refs, or sex-linked, or control probes")
    qc_parser.set_defaults(func=cli_qc)

    controls_parser = subparsers.add_parser('controls', help='Generate an XLSX QC report of control probe performance')
    controls_parser.set_defaults(func=cli_controls_report)

    report_parser = subparsers.add_parser('report', help="Run methylcheck's ReportPDF pipeline and generate a PDF.")
    report_parser.set_defaults(func=cli_ReportPDF)

    #-- when this has subparsers, use -- args = parser.parse_known_args(sys.argv[1:])
    parsed_args, func_args = parser.parse_known_args(sys.argv[1:])
    if parsed_args.debug:
        logging.basicConfig(level=logging.DEBUG)
    elif 'verbose' in func_args:
        logging.basicConfig(level=logging.INFO)
    else:
        logging.basicConfig(level=logging.WARNING)

    parsed_args.func(func_args)
    return parser

def cli_qc(cmd_args):
    parser = DefaultParser(
        prog='methylcheck qc',
        description='Run the QC pipeline on files in a given beta values file.',
    )

    parser.add_argument(
        '-d', '--data_file',
        required=True,
        type=Path,
        help="""path to a file containing sample matrix of beta or m_values. \
You can create this output from `methylprep` using the --betas flag.""",
    )

    parser.add_argument(
        '-v', '--verbose',
        help='Provide more detailed processing messages.',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '-a', '--array_type',
        choices=['27k','450k','EPIC','EPIC+'],
        required=False,
        help='Type of array being processed. If omitted, methylcheck will autodetect based on shape of data.',
    )

    parser.add_argument(
        '--exclude_sex',
        help='filters out probes on sex-chromosomes',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '--exclude_control',
        help='filters out illumina control probes',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '--exclude_probes',
        help='This will exclude all problem probes, based on several lists from academic publications. \
If you want fine-tuning-control, use the exclude_probes() function in a jupyter notebook instead. \
Also see methylcheck.list_problem_probes for more details.',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '--exclude_all',
        help='This will exclude all problem probes, sex-linked probes, and control probes at once.',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '-p', '--plot',
        nargs='*', # -p type1 type2 ... zero or more plots.
        choices=['mean_beta_plot', 'beta_density_plot',
            'cumulative_sum_beta_distribution', 'beta_mds_plot', 'all'],
        help='Select which plots to generate. Note that if you omit this, the default setting will run all plots at once. `-p all`',
        default='all',
    )

    parser.add_argument(
        '-s', '--save',
        help='By default, each plot will only appear on screen, but you can save them as png files with this option.',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '--silent',
        help='Suppress plots and requestss user inputs; for unit testing or automated pipelines.',
        action='store_true',
        default=False,
    )

    #### run_pipeline happens here ###
    #set logging/verbose
    args = parser.parse_args(cmd_args) #sys.argv[1:])

    # load the data file
    if not Path(args.data_file).is_file():
        raise ValueError("Could not find your data_file.")
    if args.data_file.suffix == '.npy':
        npy = np.load(Path(args.data_file))
        df = pd.DataFrame(npy)
    elif args.data_file.suffix == '.pkl':
        df = pd.read_pickle(args.data_file)
    else:
        raise FileNotFoundError("Could not find/read your data file. Must be .pkl or .npy file.")
    # methylprep data will be long format, with samples in columns and probes in rows. MDS transposes this.'

    # determine array type
    if args.array_type is None:
        args.array_type = detect_array(df)

    # apply some filters
    if args.exclude_all:
        df = exclude_sex_control_probes(df, args.array_type, verbose=args.verbose)
        sketchy_probe_list = list_problem_probes(args.array_type)
        #print(len(sketchy_probe_list),'sketchy probes',sketchy_probe_list[:10])
        df = exclude_probes(df, sketchy_probe_list)

    elif args.exclude_sex and args.exclude_control:
        df = exclude_sex_control_probes(df, args.array_type, verbose=args.verbose)
    elif args.exclude_sex:
        df = exclude_sex_control_probes(df, args.array_type, no_sex=True, no_control=False, verbose=args.verbose)
    elif args.exclude_control:
        df = exclude_sex_control_probes(df, args.array_type, no_sex=False, no_control=True, verbose=args.verbose)

    if args.exclude_probes and not args.exclude_all:
        sketchy_probe_list = list_problem_probes(args.array_type)
        #print(type(sketchy_probe_list),len(sketchy_probe_list))
        #print(df.shape, list(df.index)[:20])
        df = exclude_probes(df, sketchy_probe_list)

    # run calculations and plots
    if 'all' in args.plot:
        mean_beta_plot(df, verbose=args.verbose, save=args.save, silent=args.silent)
        beta_density_plot(df, verbose=args.verbose, save=args.save, silent=args.silent)
        wide_df = df.copy().transpose()
        cumulative_sum_beta_distribution(wide_df, verbose=args.verbose, save=args.save, silent=args.silent)
        beta_mds_plot(wide_df, verbose=args.verbose, save=args.save, silent=args.silent)
    else:
        if 'mean_beta_plot' in args.plot:
            mean_beta_plot(df, verbose=args.verbose, save=args.save, silent=args.silent)
        if 'beta_density_plot' in args.plot:
            beta_density_plot(df, verbose=args.verbose, save=args.save, silent=args.silent)
        if 'cumulative_sum_beta_distribution' in args.plot:
            wide_df = df.copy().transpose()
            cumulative_sum_beta_distribution(wide_df, verbose=args.verbose, save=args.save, silent=args.silent)
        if 'beta_mds_plot' in args.plot:
            df_filtered = beta_mds_plot(df, verbose=args.verbose, save=args.save, silent=args.silent)
            # also has silent and filter_stdev params to pass in.

def cli_controls_report(cmd_args):
    parser = DefaultParser(
        prog='methylcheck controls',
        description='Run the controls_report QC pipeline on files in a given folder.',
    )

    parser.add_argument(
        '-d', '--data_dir',
        required=True,
        type=Path,
        help="""Path where input files are located.""",
    )

    parser.add_argument(
        '-o', '--out_dir',
        required=False,
        type=Path,
        help="""Path where output files will be saved.""",
    )
    parser.add_argument(
        '--bg_offset',
        help='Adjust the baseline background offset value. Default is 3000.',
        type=int,
        default=3000,
    )

    parser.add_argument(
        '--colorblind',
        help='Use an alternate color palette in output XLSX file.',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '--cutoff_adjust',
        help='Adjust the pass/fail cutoff. 1.0 uses the standard limits. 2.0 will increase minimum acceptable ratio (more conservative / lower tolerance)',
        type=float,
        default=1.0,
    )

    parser.add_argument(
        '--pval_off',
        help='Flag to disable use of probe failures (pvalue) in report. Note that Illumina Legacy compatability mode does not include pval.',
        action='store_true',
        default=False,
    )

    parser.add_argument(
        '--pval_sig',
        help='p-value failure significance threshold (default <0.05). Ignore if not using pval.',
        type=float,
        default=0.05,
    )

    parser.add_argument(
        '-l', '--legacy',
        help='Change XLSX columns to exclude sex check and overall sample quality score, and match the format of Illumina quality control output files.',
        action='store_true',
        default=False,
    )

    args = parser.parse_args(cmd_args) #sys.argv[1:])
    return controls_report(
        filepath=args.data_dir,
        outfilepath=args.out_dir,
        bg_offset=args.bg_offset,
        colorblind=args.colorblind,
        cutoff_adjust=args.cutoff_adjust,
        legacy=args.legacy,
        pval=(not args.pval_off), # needed to invert kwarg, because default is to run it, and flag is only to disable it
        pval_sig=args.pval_sig)

def cli_ReportPDF(cmd_args):
    parser = DefaultParser(
        prog='methylcheck report',
        description='Create a battery of control probe performance plots, saved as a PDF.',
    )

    parser.add_argument(
        '-d', '--data_dir',
        required=False,
        type=Path,
        default='.',
        help="""Path where input files are located. Default is current folder if omitted.""",
    )

    parser.add_argument(
        '-o', '--out_dir',
        required=False,
        type=Path,
        help="""Path where output files will be saved. Defaults to same folder as 'data_dir'.""",
    )

    parser.add_argument(
        '-f', '--filename',
        required=False,
        default='multipage_pdf.pdf',
        help="""Filename of the PDF generated. Defaults to 'multipage_pdf.pdf' if omitted.""",
    )

    parser.add_argument(
        '--pval_min_percent',
        help='Set the minimum percent of probes that must pass pOOBAH for a sample to be retained. Default is 80%.',
        type=int,
        default=80,
    )

    parser.add_argument(
        '--pval_sig',
        help='p-value failure significance level. Default alpha is 0.05.',
        type=float,
        default=0.05,
    )

    args = parser.parse_args(cmd_args)
    if args.out_dir is None:
        args.out_dir = args.data_dir
    return ReportPDF(
        path=args.data_dir,
        outpath=args.out_dir,
        filename=args.filename,
        pval_min_percent=args.pval_min_percent,
        pval_cutoff=args.pval_sig,
        runme=True) # <--- this autocompletes report and closes it at end of __init__ in one step.

def cli_app():
    build_parser()

max_stars_repo_name: FoxoTech/methylcheck

max_issues_repo_name: FoxoTech/methylcheck

max_forks_repo_name: LifeEGX/methQC


text:
# Standard library imports
import os
from typing import Dict, Iterator, List, Optional

# Third-party imports
import numpy as np
import pandas as pd
from pathlib import Path

# First-party imports
from gluonts.core.component import validated
from gluonts.dataset.common import Dataset
from gluonts.time_feature.seasonality import get_seasonality
from gluonts.model.forecast import SampleForecast
from gluonts.model.predictor import RepresentablePredictor

# https://stackoverflow.com/questions/25329955/check-if-r-is-installed-from-python
from subprocess import Popen, PIPE

proc = Popen(["which", "R"], stdout=PIPE, stderr=PIPE)
R_IS_INSTALLED = proc.wait() == 0

try:
    import rpy2.robjects.packages as rpackages
    from rpy2 import rinterface, robjects
    from rpy2.rinterface import RRuntimeError
except ImportError as e:
    rpy2_error_message = str(e)
    RPY2_IS_INSTALLED = False
else:
    RPY2_IS_INSTALLED = True

USAGE_MESSAGE = """
The RHierarchicalForecastPredictor is a thin wrapper for calling the R HTS package.
In order to use it you need to install R and run

pip install rpy2  # version: rpy2>=2.9.*,<3.*

R -e 'install.packages(c("hts"), repos="https://cloud.r-project.org")'
"""

POINT_HIERARCHICAL_METHODS = [
    "naive_bottom_up", "mint", "erm", "ermParallel"
]

SAMPLE_HIERARCHICAL_METHODS = [
    "depbu_mint"
]

SUPPORTED_METHODS = POINT_HIERARCHICAL_METHODS + SAMPLE_HIERARCHICAL_METHODS

class RHierarchicalForecastPredictor(RepresentablePredictor):
    """
    The RHierarchicalForecastPredictor is a thin wrapper for calling the R HTS package.
    In order to use it you need to install R and run

    pip install rpy2

    R -e 'install.packages(c("hts"), repos="https://cloud.r-project.org")'

    Parameters
    ----------
    freq
    prediction_length
    target_dim
    num_bottom_ts
    nodes
    tags
    is_hts
    method_name
        hierarchical forecast or reconciliation method to be used; mutst be one of:
        "mint", "naive_bottom_up", "erm", "mint_ols", "depbu_mint", "ermParallel"
    fmethod
    algorithm
    covariance
    nonnegative
    period
    correction
    seasonal
    stationary
    params
    lead_time
    """

    @validated()
    def __init__(
        self,
        freq: str,
        prediction_length: int,
        target_dim: int,
        num_bottom_ts: int,
        nodes: List,
        tags: np.ndarray,
        is_hts: bool,
        method_name: str = "mint",
        fmethod: str = "arima",
        algorithm: str = "cg",
        covariance: str = "shr",
        nonnegative: bool = False,
        period: int = None,
        correction: bool = True,
        seasonal: bool = True,
        stationary: bool = False,
        params: Optional[Dict] = None,
        lead_time: int = 0,
    ) -> None:

        super().__init__(freq=freq, prediction_length=prediction_length)

        try:
            from rpy2 import robjects, rinterface
            import rpy2.robjects.packages as rpackages
        except ImportError as e:
            raise ImportError(str(e) + USAGE_MESSAGE) from e

        if not R_IS_INSTALLED:
            raise ImportError("R is not Installed! \n " + USAGE_MESSAGE)

        if not RPY2_IS_INSTALLED:
            raise ImportError(rpy2_error_message + USAGE_MESSAGE)

        self._robjects = robjects
        self._rinterface = rinterface
        self._rinterface.initr()
        self._rpackages = rpackages

        this_dir = os.path.dirname(os.path.realpath(__file__))
        this_dir = this_dir.replace("\\", "/")  # for windows
        r_files = [
            n[:-2] for n in os.listdir(f"{this_dir}/R/") if "hierarchical" in n
        ]

        for n in r_files:
            try:
                path = Path(this_dir, "R", f"{n}.R")
                robjects.r(f'source("{path}")'.replace("\\", "\\\\"))
            except RRuntimeError as er:
                raise RRuntimeError(str(er) + USAGE_MESSAGE) from er

        assert (
            method_name in SUPPORTED_METHODS
        ), f"method {method_name} is not supported please use one of {SUPPORTED_METHODS}"

        self.method_name = method_name

        self._stats_pkg = rpackages.importr('stats')
        self._hts_pkg = rpackages.importr('hts')
        self._r_method = robjects.r[method_name]

        self.prediction_length = prediction_length
        self.freq = freq
        self.period = period if period is not None else get_seasonality(freq)
        self.target_dim = target_dim
        self.num_bottom_ts = num_bottom_ts
        self.nodes = nodes
        self.tags = tags
        self.params = {
            'prediction_length': self.prediction_length,
            'output_types': ['samples'],
            'frequency': self.period,
            "fmethod": fmethod,
            "algorithm": algorithm,
            "covariance": covariance,
            "nonnegative": nonnegative,
            "seasonal": seasonal,
            "stationary": stationary,
            "correction": correction,
        }
        if params is not None:
            self.params.update(params)
        self.lead_time = lead_time
        self.is_hts = is_hts

    def _unlist(self, l):
        if type(l).__name__.endswith("Vector") and type(l).__name__ != "ListVector":
            return [self._unlist(x) for x in l]
        elif type(l).__name__ == "ListVector":
            return [self._unlist(x) for x in l]
        elif type(l).__name__ == "Matrix":
            return np.array(l)
        else:
            return l

    def _run_r_forecast(self, d, params, save_info):
        buf = []

        def save_to_buf(x):
            buf.append(x)

        def dont_save(x):
            pass

        f = save_to_buf if save_info else dont_save

        # save output from the R console in buf
        self._rinterface.set_writeconsole_regular(f)
        self._rinterface.set_writeconsole_warnerror(f)

        params["freq_period"] = self.period
        print("method_name:", self.method_name)

        r_params = self._robjects.vectors.ListVector(params)

        if self.is_hts:
            nodes = self._robjects.ListVector.from_length(len(self.nodes))
            for idx, elem in enumerate(self.nodes):
                if not isinstance(elem, list):
                    elem = [elem]
                nodes[idx] = self._robjects.IntVector(elem)
        else:
            nodes_temp = []
            for idx, elem in enumerate(self.nodes):
                nodes_temp.extend(elem)
            nodes = self._robjects.r.matrix(self._robjects.IntVector(nodes_temp),
                                            ncol=len(self.nodes), nrow=len(self.nodes[0]), byrow=False).transpose()

        num_ts, nobs = d["target"].shape
        tags_r = self._robjects.r.matrix(self._robjects.StrVector(self.tags.flatten()),
                                           ncol=self.tags.shape[1], nrow=self.tags.shape[0], byrow=True)

        y_bottom = self._robjects.r.matrix(self._robjects.FloatVector(d["target"].flatten()),
                                      ncol=nobs, nrow=num_ts, byrow=True).transpose()

        y_bottom_ts = self._stats_pkg.ts(y_bottom, frequency=self.period)

        if self.is_hts:
            hier_ts = self._hts_pkg.hts(y_bottom_ts, nodes=nodes)  # get nodes correctly here
        else:
            hier_ts = self._hts_pkg.gts(y_bottom_ts, groups=nodes)  # get nodes correctly here

        if self.method_name == "depbu_mint":
            forecast = self._r_method(y_bottom_ts, tags_r, r_params)
        else:
            forecast = self._r_method(hier_ts, r_params)

        all_forecasts = list(forecast)
        if self.method_name in POINT_HIERARCHICAL_METHODS:
            assert len(all_forecasts) == self.target_dim * self.prediction_length
            hier_point_forecasts = np.reshape(
                list(forecast), (self.target_dim, self.prediction_length)
            ).transpose()
            hier_forecasts = np.tile(
                hier_point_forecasts, (params["num_samples"], 1, 1)
            )
        else:
            hier_forecasts = np.reshape(
                list(forecast), (params["num_samples"], self.prediction_length, self.target_dim), order='F'
            )
        forecast_dict = dict(
            samples=hier_forecasts
        )

        self._rinterface.set_writeconsole_regular(
            self._rinterface.consolePrint
        )
        self._rinterface.set_writeconsole_warnerror(
            self._rinterface.consolePrint
        )
        return forecast_dict, buf

    def predict(
        self, dataset: Dataset, num_samples=100, save_info=True, **kwargs
    ) -> Iterator[SampleForecast]:
        for data in dataset:
            data["target"] = data["target"][-self.num_bottom_ts:, :]

            if self.method_name == "depbu_mint":
                train_length = data["target"].shape[1]

                # this method expects train + test target; add dummy values for the "test" part
                data["target"] = np.concatenate(
                    (data["target"], data["target"][:, -self.prediction_length:]),
                    axis=1
                )
                self.params["train_length"] = train_length

            params = self.params.copy()
            params['num_samples'] = num_samples

            forecast_dict, console_output = self._run_r_forecast(
                data, params, save_info=save_info
            )

            if self.method_name == "depbu_mint":
                # remove the extra `prediction_length` target that was added
                data["target"] = data["target"][:, :-self.prediction_length]
                assert data['target'].shape[1] == train_length

            samples = np.array(forecast_dict['samples'])
            expected_shape = (params['num_samples'], self.prediction_length, self.target_dim)
            assert (samples.shape == expected_shape),\
                f"Expected shape {expected_shape} but found {samples.shape}"

            info = (
                {'console_output': '\n'.join(console_output)}
                if save_info
                else None
            )

            yield SampleForecast(
                samples,
                # TODO: to use gluonts.support.pandas.forecast_start(data), after updating it so that it understands
                #  multivariate target!
                pd.date_range(data["start"], periods=data["target"].shape[-1] + 1, freq=self.freq)[-1],
                self.freq,
                info=info,
                item_id=data.get("item_id", None),
            )


text:
from tool import Tool
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, r2_score

new_class = Tool()

df = pd.read_csv('resources/postcodes_sampled_15000.csv') # train model with this

df2 = pd.read_csv('resources/test_postcodes_5000.csv') # postcodes we want to find

df3 = pd.read_csv('resources/postcodes_sampled_5000.csv') #true value of those 5000 postcodes

locations = df2['postcode'].tolist()

y_true = df3['medianPrice']

y_pred = new_class.get_median_house_price_estimate(locations, method=2)

<<<<<<< HEAD
prices = new_class.get_annual_flood_risk(locations)
=======
print(np.sqrt(mean_squared_error(y_true.values, y_pred.values)))
print(r2_score(y_true.values, y_pred.values))
>>>>>>> d070fedf959f9067c7faf77db05d0ff9fb44d763