text:
import tkinter as tk
import seaborn as sns
from tkinter.filedialog import askopenfilename
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tkinter import filedialog, messagebox, ttk
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg,FigureCanvas
root=tk.Tk()
sns.set_theme(color_codes=True)

def file2():
    return malai

def file1():
    global malai
    malai=filedialog.askopenfilename(initialdir="/",
                                              title="Select A File",
                                              filetype=(("CSC files", "*.csv"),("xlsx files", "*.xlsx"),("All Files", "*.*")))
    return malai

def get():

    # initalise the tkinter GUI
    root = tk.Tk()

    root.geometry("500x500") # set the root dimensions
    root.pack_propagate(False) # tells the root to not let the widgets inside it determine its size.
    root.resizable(0, 0) # makes the root window fixed in size.

    # Frame for TreeView
    frame1 = tk.LabelFrame(root, text="Excel Data")
    frame1.place(height=250, width=500)

    # Frame for open file dialog
    file_frame = tk.LabelFrame(root, text="Open File")
    file_frame.place(height=100, width=400, rely=0.65, relx=0)

    # Buttons
    button1 = tk.Button(file_frame, text="Browse A File", command=lambda: File_dialog())
    button1.place(rely=0.65, relx=0.50)

    button2 = tk.Button(file_frame, text="Load File", command=lambda: Load_excel_data())
    button2.place(rely=0.65, relx=0.30)

    # The file/file path text
    label_file = ttk.Label(file_frame, text="No File Selected")
    label_file.place(rely=0, relx=0)

    ## Treeview Widget
    tv1 = ttk.Treeview(frame1)
    tv1.place(relheight=1, relwidth=1) # set the height and width of the widget to 100% of its container (frame1).

    treescrolly = tk.Scrollbar(frame1, orient="vertical", command=tv1.yview) # command means update the yaxis view of the widget
    treescrollx = tk.Scrollbar(frame1, orient="horizontal", command=tv1.xview) # command means update the xaxis view of the widget
    tv1.configure(xscrollcommand=treescrollx.set, yscrollcommand=treescrolly.set) # assign the scrollbars to the Treeview Widget
    treescrollx.pack(side="bottom", fill="x") # make the scrollbar fill the x axis of the Treeview widget
    treescrolly.pack(side="right", fill="y") # make the scrollbar fill the y axis of the Treeview widget

    def File_dialog():
        """This Function will open the file explorer and assign the chosen file path to label_file"""
        filename = file1()
        label_file["text"] = filename
        return None

    def Load_excel_data():
        """If the file selected is valid this will load the file into the Treeview"""
        file_path = label_file["text"]
        try:
            excel_filename = r"{}".format(file_path)
            if excel_filename[-4:] == ".csv":
                df = pd.read_csv(excel_filename)
            else:
                df = pd.read_excel(excel_filename)

        except ValueError:
            tk.messagebox.showerror("Information", "The file you have chosen is invalid")
            return None
        except FileNotFoundError:
            tk.messagebox.showerror("Information", f"No such file as {file_path}")
            return None

        clear_data()
        tv1["column"] = list(df.columns)
        tv1["show"] = "headings"
        for column in tv1["columns"]:
            tv1.heading(column, text=column) # let the column heading = column name

        df_rows = df.to_numpy().tolist() # turns the dataframe into a list of lists
        for row in df_rows:
            tv1.insert("", "end", values=row) # inserts each list into the treeview. For parameters see https://docs.python.org/3/library/tkinter.ttk.html#tkinter.ttk.Treeview.insert
        return None

    def clear_data():
        tv1.delete(*tv1.get_children())
        return None

    root.mainloop()

def show():
    filename=file2()
    my_data=pd.read_csv(filename)
    root1=tk.Tk()
    figure1=plt.figure(figsize=(5,4),dpi=100)
    ax1=figure1.add_subplot(211)
    bar1=FigureCanvas(figure1,root1)
    bar1.get_tk_widget().pack(side=tk.LEFT,fill=tk.BOTH)
    df1=my_data[[str(ans2.get()),str(ans3.get())]].groupby(str(ans2.get())).sum()
    df1.plot(kind="line",ax=ax1,legend=True)
    ax1.set_title(str(ans2.get())+" vs "+str(ans3.get()))

    figure2=plt.figure(figsize=(5,5),dpi=100)
    ax2=figure2.add_subplot(111)
    line=FigureCanvas(figure2,root1)
    line.get_tk_widget().pack(side=tk.LEFT,fill=tk.BOTH)
    df1=my_data[[str(ans2.get()),str(ans3.get())]].groupby(str(ans2.get())).sum()
    df1.plot(kind="hist",ax=ax2,legend=True)
    ax1.set_title(str(ans2.get())+" vs "+str(ans3.get()))

    figure3=plt.figure(figsize=(5,4),dpi=100)
    ax3=figure3.add_subplot(111)
    ax3.scatter(my_data[str(ans2.get())],my_data[str(ans3.get())],color="g")
    scatter3=FigureCanvas(figure3,root1)
    scatter3.get_tk_widget().pack(side=tk.LEFT,fill=tk.BOTH)

    ax3.set_title(str(ans2.get())+" vs "+str(ans3.get()))

    ax4=figure1.add_subplot(212)
    bar1=FigureCanvas(figure1,root1)
    df1=my_data[[str(ans2.get()),str(ans3.get())]].groupby(str(ans2.get())).sum()
    df1.plot(kind="box",ax=ax4,legend=True)

    root1.mainloop()

canvas=tk.Canvas(root,height=500,width=500)
canvas.pack()

label=tk.Label(root,text="PlotData",bg="gray90")
canvas.create_window(110,100,window=label)

entry1=tk.Entry(root,width=27)
canvas.create_window(225,100,window=entry1)
button=tk.Button(root,text="show",command=get,bg="gray90")
canvas.create_window(225,130,window=button)

button2=tk.Button(root,text="plot",command=show,bg="gray90")
canvas.create_window(350,230,window=button2)

label2=tk.Label(root,text="Column1",bg="gray90")
canvas.create_window(300,160,window=label2)
ans2=tk.Entry(root,width=13)
canvas.create_window(370,160,window=ans2)

label3=tk.Label(root,text="Column2",bg="gray90")
canvas.create_window(300,190,window=label3)
ans3=tk.Entry(root,width=13)
canvas.create_window(370,190,window=ans3)

root.mainloop()

max_stars_repo_name: Annamalaisaravanan/Plotapp

max_issues_repo_name: Annamalaisaravanan/Plotapp

max_forks_repo_name: Annamalaisaravanan/Plotapp


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Apply bias field correction to DWI. This script doesn't compute the bias
field itself. It ONLY applies an existing bias field. Use the ANTs
N4BiasFieldCorrection executable to compute the bias field
"""

from __future__ import division
from past.utils import old_div
import argparse

import nibabel as nib
import numpy as np

from scilpy.io.utils import (
    add_overwrite_arg, assert_inputs_exist, assert_outputs_exist)

def _build_arg_parser():
    parser = argparse.ArgumentParser(
        description=__doc__, formatter_class=argparse.RawTextHelpFormatter)
    parser.add_argument('dwi', help='DWI Nifti image')
    parser.add_argument('bias_field', help='Bias field Nifti image')
    parser.add_argument('output', help='Corrected DWI Nifti image')
    parser.add_argument('--mask',
                        help='Apply bias field correction only in the region '
                             'defined by the mask')
    add_overwrite_arg(parser)
    return parser

def _rescale_intensity(val, slope, in_max, bc_max):
    return in_max - slope * (bc_max - val)

# https://github.com/stnava/ANTs/blob/master/Examples/N4BiasFieldCorrection.cxx
def _rescale_dwi(in_data, bc_data, mask_data=None):
    nz_in_data = in_data
    nz_bc_data = bc_data
    nz_mask_data = None

    if mask_data is not None:
        nz_mask_data = np.nonzero(mask_data)
        nz_in_data = in_data[nz_mask_data]
        nz_bc_data = bc_data[nz_mask_data]

    in_min = np.amin(nz_in_data)
    in_max = np.amax(nz_in_data)
    bc_min = np.amin(nz_bc_data)
    bc_max = np.amax(nz_bc_data)

    slope = old_div((in_max - in_min), (bc_max - bc_min))

    rescale_func = np.vectorize(_rescale_intensity, otypes=[np.float])
    rescaled_data = rescale_func(nz_bc_data, slope, in_max, bc_max)

    if mask_data is not None:
        bc_data[nz_mask_data] = rescaled_data
    else:
        bc_data = rescaled_data

    return bc_data

def main():
    parser = _build_arg_parser()
    args = parser.parse_args()

    assert_inputs_exist(parser, [args.dwi, args.bias_field], args.mask)
    assert_outputs_exist(parser, args, args.output)

    dwi_img = nib.load(args.dwi)
    dwi_data = dwi_img.get_data()

    bias_field_img = nib.load(args.bias_field)
    bias_field_data = bias_field_img.get_data()

    mask_data = nib.load(args.mask).get_data() if args.mask else None
    nuc_dwi_data = np.divide(dwi_data, bias_field_data[..., np.newaxis])
    rescaled_nuc_data = _rescale_dwi(dwi_data, nuc_dwi_data, mask_data)

    nib.save(nib.Nifti1Image(rescaled_nuc_data, dwi_img.affine,
                             dwi_img.header),
             args.output)

if __name__ == "__main__":
    main()

ext: py

max_stars_repo_name: fullbat/scilpy

max_issues_repo_name: fullbat/scilpy

max_forks_repo_name: fullbat/scilpy

reason: import numpy


text:
"""
   Data cleaning and formating module
"""

from os.path import join
from warnings import simplefilter
import re

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import (
    Polygon,
    MultiPolygon
)
from shapely.wkt import loads

from utilities import cwd
from database import DataBase

# inputs
COUNTY_SHAPES = join(cwd(), 'shapes', 'counties_500k', 'cb_2018_us_county_500k.shx')
STATE_SHAPES = join(cwd(), 'shapes', 'states_500k', 'cb_2018_us_state_500k.shx')

# outputs
US_MAP_TABLE = 'us_map'
STATE_MAP_TABLE = 'state_map'

def remove_islands(map_file, min_area=100000000):
    """Remove small polygons

    Arguments:
        map_file {GeoDataFrame} -- data with a 'geometry' column and
                                   expected crs = 'EPSG:3395

    Keyword Arguments:
        min_area {int} -- polygon area to remove (default: {100000000})

    Returns:
        GeoDataFrame -- data with small polygon removed
    """
    records = map_file.to_dict(orient='records')
    for record in records:
        if record['geometry'].type == 'MultiPolygon':
            multipoly = list(record['geometry'])
            polies = [poly for poly in multipoly if poly.area > min_area]
            if polies:
                record['geometry'] = MultiPolygon(polies)
            else:
                record['geometry'] = max(
                    record['geometry'], key=lambda a: a.area)

    new_map = gpd.GeoDataFrame(records)

    return new_map

def mround(match):
    """Return formated string
    """
    return "{:.2f}".format(float(match.group()))

def transform_alaska_hawaii(us_map, state_map):
    """[summary]

    Arguments:
        us_map {GeoDataFrame} -- us counties geodataframe
        state_map {[type]} -- us states geodataframe

    Returns:
        tuple -- tuple of dataframes with alaska and hawaii transformed
    """
    # transform alaska and hawaii
    texas = state_map['STATEFP'] == '48'
    [(xc_tx, yc_tx)] = list(state_map.loc[texas, 'geometry'].iat[0].centroid.coords)

    alaska = state_map['STATEFP'] == '02'
    xc_ak = state_map.loc[alaska, 'xc'].iat[0]
    yc_ak = state_map.loc[alaska, 'yc'].iat[0]

    hawaii = state_map['STATEFP'] == '15'
    xc_hi = state_map.loc[hawaii, 'xc'].iat[0]
    yc_hi = state_map.loc[hawaii, 'yc'].iat[0]

    state_map.loc[alaska, 'geometry'] = state_map[alaska].scale(
        xfact=0.19, yfact=0.19, zfact=1.0, origin='centroid')
    state_map.loc[alaska, 'geometry'] = state_map[alaska].translate(
        xoff=xc_tx - xc_ak - 2.7e6, yoff=yc_tx - yc_ak - 5.0e5, zoff=0.0)
    state_map.loc[hawaii, 'geometry'] = state_map[hawaii].translate(
        xoff=xc_tx - xc_hi - 1.5e6, yoff=yc_tx - yc_hi - 7.5e5, zoff=0.0)

    alaska_ct = us_map['STATEFP'] == '02'
    hawaii_ct = us_map['STATEFP'] == '15'
    us_map.loc[alaska_ct, 'geometry'] = us_map[alaska_ct].scale(
        xfact=0.19, yfact=0.19, zfact=1.0, origin=(xc_ak, yc_ak, 0))
    us_map.loc[alaska_ct, 'geometry'] = us_map[alaska_ct].translate(
        xoff=xc_tx - xc_ak - 2.7e6, yoff=yc_tx - yc_ak - 5.0e5, zoff=0.0)
    us_map.loc[hawaii_ct, 'geometry'] = us_map[hawaii_ct].translate(
        xoff=xc_tx - xc_hi - 1.5e6, yoff=yc_tx - yc_hi - 7.5e5, zoff=0.0)

    return us_map, state_map

def group_nyc_counties(us_map):
    """Group NYC counties as New York City

    Arguments:
        us_map {GeoDataFrame} -- county metadata

    Returns:
        {GeoDataFrame} -- county metadata
    """
    # new york city counties
    nyc_counties = {'Queens': '36081',
                    'Bronx': '36005',
                    'Richmond': '36085',
                    'New York': '36061',
                    'Kings': '36047'}

    nyc = us_map.loc[us_map['county_id'].isin(nyc_counties.values()), :]
    nyc = nyc.copy(deep=True)

    nyc['county_id'] = '36061'
    nyc = nyc.dissolve(by='county_id', aggfunc='sum').reset_index()
    nyc['name'] = 'New York City'
    nyc['state_id'] = '36'

    counties = us_map.loc[~us_map['county_id'].isin(nyc_counties.values()), :]
    counties = counties.copy(deep=True)
    counties = counties.append(nyc, ignore_index=True)

    return counties

def get_maps(us_map, state_map):
    """Shapes transformation of county and state maps

    Arguments:
        us_map {GeoDataFrame} -- us county map
        state_map {[type]} -- us state map

    Returns:
        tuple -- transformed us_map and state_map
    """
    # remove unwanted states
    drop = ['72', '78', '69', '66', '60']
    us_map = us_map[~us_map['STATEFP'].isin(drop)].copy(deep=True)
    state_map = state_map[~state_map['STATEFP'].isin(drop)].copy(deep=True)

    # coordinate reference system
    simplefilter(action='ignore', category=FutureWarning)
    us_map = us_map.to_crs('EPSG:3395')
    state_map = state_map.to_crs('EPSG:3395')
    simplefilter(action='default')

    # remove islands (first pass)
    us_map = remove_islands(us_map)
    state_map = remove_islands(state_map)

    # add centroids to state map
    points = list(state_map['geometry'].centroid)
    state_map['xc'] = [point.coords[0][0] for point in points]
    state_map['yc'] = [point.coords[0][1] for point in points]

    # scale Alaska and move Alaska and Hawaii under California
    us_map, state_map = transform_alaska_hawaii(us_map, state_map)

    # remove small island (second pass)
    us_map = remove_islands(us_map)
    state_map = remove_islands(state_map)

    # smooth out polygons
    dist = 3500
    us_map['geometry'] = gpd.GeoSeries([x if x.type != 'Polygon' else Polygon(
        x.buffer(dist, join_style=1).buffer(-1 * dist, join_style=1))
                                        for x in us_map['geometry']])

    state_map['geometry'] = gpd.GeoSeries([x if x.type != 'Polygon' else Polygon(
        x.buffer(dist, join_style=1).buffer(-1 * dist, join_style=1))
                                           for x in state_map['geometry']])

    # simplify geometry
    tolerance = 2000
    us_map['geometry'] = us_map['geometry'].simplify(tolerance)
    state_map['geometry'] = state_map['geometry'].simplify(tolerance)

    # dissolve us map by county
    us_map['DISSOLVE'] = us_map['NAME'] + us_map['STATEFP']
    us_map = us_map.dissolve('DISSOLVE')
    us_map.reset_index(drop=True, inplace=True)

    # dissolve state map by state
    state_map = state_map.dissolve('STATEFP')
    state_map.reset_index(inplace=True)

    # remove holes from polygons
    us_map['geometry'] = gpd.GeoSeries(
        [x if x.type != 'Polygon' else Polygon(x.exterior) for x in us_map['geometry']])
    state_map['geometry'] = gpd.GeoSeries(
        [x if x.type != 'Polygon' else Polygon(x.exterior) for x in state_map['geometry']])

    # create unique ids
    us_map['fips'] = us_map['GEOID'].astype(int)

    # add population data from lookup
    pop = pd.read_csv(join(cwd(), 'input', 'fips_county_pop.csv'))
    pop = pop.set_index('fips')['population'].to_dict()
    us_map['pop'] = us_map['fips'].map(pop)
    us_map['pop'].fillna(0)

    pop = pd.read_csv(join(cwd(), 'input', 'state_name_pop.csv'))
    pop = pop.set_index('state')['pop'].to_dict()
    state_map['pop'] = state_map['NAME'].map(pop)

    # select columns, covert to lowercase, and change column names
    keep = ['GEOID', 'STATEFP', 'NAME', 'ALAND', 'AWATER', 'pop', 'geometry']
    us_map = us_map[keep].dropna().copy(deep=True)
    us_map = us_map.rename(columns={'GEOID': 'county_id', 'STATEFP': 'state_id'})
    us_map.columns = [name.lower() for name in us_map.columns]

    keep = ['STATEFP', 'NAME', 'STUSPS', 'ALAND', 'AWATER', 'pop', 'geometry']
    state_map = state_map[keep].dropna().copy(deep=True)
    state_map = state_map.rename(columns={'STATEFP': 'state_id', 'STUSPS': 'abbr'})
    state_map.columns = [name.lower() for name in state_map.columns]

    # convert land area (aland) and water area (awater) to square km
    us_map[['aland', 'awater']] = np.int32(us_map[['aland', 'awater']]/1e6)
    state_map[['aland', 'awater']] = np.int32(state_map[['aland', 'awater']]/1e6)

    # group new york city counties as one (new york city)
    us_map = group_nyc_counties(us_map)

    # finally round geometry coords to two decimal points
    us_map['geometry'] = us_map['geometry'].apply(
        lambda x: loads(re.sub(r'\d*\.\d+', mround, x.wkt)))
    state_map['geometry'] = state_map['geometry'].apply(
        lambda x: loads(re.sub(r'\d*\.\d+', mround, x.wkt)))

    return us_map, state_map

def maps_to_database():
    """Refresh database with counties and states map
    """
    us_map = gpd.read_file(COUNTY_SHAPES, encoding='UTF-8')
    state_map = gpd.read_file(STATE_SHAPES, encoding='UTF-8')

    # format shapes files
    us_map, state_map = get_maps(us_map, state_map)

    # to database
    _db = DataBase()
    _db.add_geotable(US_MAP_TABLE, us_map.set_index('county_id'))
    _db.add_geotable(STATE_MAP_TABLE, state_map.set_index('state_id'))
    _db.close()

if __name__ == "__main__":

    maps_to_database()


text:
'Test of ALU and load instructions'

import numpy as np
from random import getrandbits

from videocore.assembler import qpu
from videocore.driver import Driver

@qpu
def boilerplate(asm, f, nin, nout):
    setup_dma_load(nrows=nin)
    start_dma_load(uniform)
    wait_dma_load()
    setup_vpm_read(nrows=nin)
    setup_vpm_write()

    f(asm)

    setup_dma_store(nrows=nout)
    start_dma_store(uniform)
    wait_dma_store()
    exit()

def run_code(code, X, output_shape, output_type):
    with Driver() as drv:
        X = drv.copy(X)
        Y = drv.alloc(output_shape, dtype=output_type)
        drv.execute(
                n_threads=1,
                program=drv.program(boilerplate, code, X.shape[0],
                                    output_shape[0]),
                uniforms=[X.address, Y.address]
                )
        return np.copy(Y)

#============================= Integer arithmetic =============================

@qpu
def int_ops(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    for op in ['iadd', 'isub', 'imin', 'imax']:
        getattr(asm, op)(r2, r0, r1)
        mov(vpm, r2)

def test_int_ops():
    X = np.array(
            [getrandbits(32) for i in range(2*16)]
            ).reshape(2, 16).astype('int32')
    Y = run_code(int_ops, X, (4, 16), 'int32')
    assert all(Y[0] == X[0] + X[1])
    assert all(Y[1] == X[0] - X[1])
    assert all(Y[2] == np.minimum(X[0], X[1]))
    assert all(Y[3] == np.maximum(X[0], X[1]))

@qpu
def imul24_op(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    imul24(r2, r0, r1)
    mov(vpm, r2)

def test_imul24_op():
    X = np.random.randint(0, 2**24, (2, 16)).astype('uint32')
    Y = run_code(imul24_op, X, (1, 16), 'uint32')
    assert all((X[0]*X[1]).astype('uint32') == Y[0])

#========================== Bitwise binary operation ==========================

@qpu
def bitwise_bin_ops(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    for op in ['band', 'bor', 'bxor']:
        getattr(asm, op)(r2, r0, r1)
        mov(vpm, r2)

def test_bitwise_bin_ops():
    X = np.array(
            [getrandbits(32) for i in range(2*16)]
            ).reshape(2, 16).astype('uint32')
    Y = run_code(bitwise_bin_ops, X, (3, 16), 'uint32')
    assert all(Y[0] == X[0] & X[1])
    assert all(Y[1] == X[0] | X[1])
    assert all(Y[2] == X[0] ^ X[1])

#========================== Bitwise unary operation ==========================

def count_leading_zeros(n):
    bit = 0x80000000
    count = 0
    while not (n & bit):
        count += 1
        bit >>= 1
    return count

@qpu
def bitwise_un_ops(asm):
    mov(r1, vpm)
    for op in ['bnot', 'clz']:
        getattr(asm, op)(r2, r1)
        mov(vpm, r2)

def test_bitwise_un_ops():
    X = np.array(
            [getrandbits(32) for i in range(16)]
            ).reshape(1, 16).astype('uint32')
    Y = run_code(bitwise_un_ops, X, (2, 16), 'uint32')
    assert all(Y[0] == ~X[0])
    assert all(Y[1] == np.vectorize(count_leading_zeros)(X[0]))

#============================== Shift operation ===============================

def rotate_right(n, s):
    return ((n << (32-s)) | (n >> s)) & 0xffffffff

@qpu
def shift_ops(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    for op in ['shr', 'asr', 'ror', 'shl']:
        getattr(asm, op)(r2, r0, r1)
        mov(vpm, r2)

def test_shift_ops():
    X = np.zeros((2, 16), dtype='uint32')
    X[0] = np.array([getrandbits(32) for i in range(16)])
    X[1] = np.random.randint(0, 32, 16)
    Y = run_code(shift_ops, X, (4, 16), 'uint32')
    assert all(X[0] >> X[1] == Y[0])
    assert all(X[0].astype('int32') >> X[1] == Y[1].astype('int32'))
    assert all(np.vectorize(rotate_right)(X[0], X[1]) == Y[2])
    assert all(X[0] << X[1] == Y[3])

#========================= Floating-point arithmetic ==========================

@qpu
def float_ops(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    for op in ['fadd', 'fsub', 'fmul', 'fmin', 'fmax', 'fminabs', 'fmaxabs']:
        getattr(asm, op)(r2, r0, r1)
        mov(vpm, r2)

def test_float_ops():
    X = np.random.randn(2, 16).astype('float32')
    Y = run_code(float_ops, X, (7, 16), 'float32')

    assert np.allclose(X[0] + X[1], Y[0], rtol=1e-6)
    assert np.allclose(X[0] - X[1], Y[1], rtol=1e-6)
    assert np.allclose(X[0] * X[1], Y[2], rtol=1e-6)
    assert all(np.minimum(X[0], X[1]) == Y[3])
    assert all(np.maximum(X[0], X[1]) == Y[4])
    assert all(np.minimum(np.abs(X[0]), np.abs(X[1])) == Y[5])
    assert all(np.maximum(np.abs(X[0]), np.abs(X[1])) == Y[6])

#============================== Type conversion ===============================

@qpu
def type_conv(asm):
    mov(r0, vpm)

    ftoi(r1, r0)
    itof(r2, r1)

    mov(vpm, r1)
    mov(vpm, r2)

def test_type_conv():
    X = np.random.randn(1, 16).astype('float32') * 1000

    Y = run_code(type_conv, X, (2, 16), 'int32')

    assert all(X[0].astype('int32') == Y[0])
    assert all(X[0].astype('int32') == np.ndarray(16, 'float32', Y[1]))

#======================== 8-bit saturation arithmetic =========================

@qpu
def v8sat_ops(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    for op in ['v8adds', 'v8subs', 'v8muld', 'v8min', 'v8max']:
        getattr(asm, op)(r2, r0, r1)
        mov(vpm, r2)

def test_v8sat_ops():
    X = np.array(
            [getrandbits(8) for i in range(2*4*16)]
            ).reshape(2, 4*16).astype('uint8')
    Y = run_code(v8sat_ops, X, (5, 4*16), 'uint8')

    u_X = X
    s_X = X.astype('int32')
    assert all(np.clip(s_X[0] + s_X[1], 0, 255) == Y[0])
    assert all(np.clip(s_X[0] - s_X[1], 0, 255) == Y[1])
    assert all(np.abs((s_X[0]*s_X[1]/255.0).astype('int32') - Y[2]) <= 1)
    assert all(np.minimum(u_X[0], u_X[1]) == Y[3])
    assert all(np.maximum(u_X[0], u_X[1]) == Y[4])

#============================== Small immediate ===============================

@qpu
def small_imm_int(asm):
    mov(r0, vpm)
    for imm in range(-16, 16):
        iadd(vpm, r0, imm)

def test_small_imm_int():
    X = np.array(
            [getrandbits(32) for i in range(16)]
            ).reshape(1, 16).astype('int32')
    Y = run_code(small_imm_int, X, (32, 16), 'int32')
    for i, imm in enumerate(range(-16, 16)):
        assert all(X[0] + imm == Y[i])

@qpu
def small_imm_float(asm):
    mov(r0, vpm)
    for e in range(-8, 8):
        fmul(vpm, r0, 2.0**e)

def test_small_imm_float():
    X = np.random.randn(1, 16).astype('float32')
    Y = run_code(small_imm_float, X, (16, 16), 'float32')
    for i, e in enumerate(range(-8, 8)):
        assert np.allclose(X[0] * 2.0**e, Y[i], rtol=1e-6)

#=============================== Load operation ===============================

@qpu
def load_two_dest(asm):
    mov(r1, vpm)
    ldi(r2, r3, 0x12345678)
    iadd(vpm, r1, r2)
    iadd(vpm, r1, r3)

def test_load_two_dest():
    X = np.ones((1, 16), dtype='uint32')
    Y = run_code(load_two_dest, X, (2, 16), 'uint32')

    assert np.all(X[0] + 0x12345678 == Y[0])
    assert np.all(X[0] + 0x12345678 == Y[1])

#=========================== Per-element immediate ============================

PER_ELMT_UNSIGNED_VALUES = np.random.randint(0, 4, 16)
PER_ELMT_SIGNED_VALUES = np.random.randint(-2, 2, 16)

@qpu
def per_elmt_imm(asm):
    mov(r0, vpm)
    ldi(r1, PER_ELMT_UNSIGNED_VALUES)
    iadd(vpm, r0, r1)
    ldi(r2, PER_ELMT_SIGNED_VALUES)
    iadd(vpm, r0, r2)

def test_per_elmt_imm():
    X = np.array(
            [getrandbits(32) for i in range(16)]
            ).reshape(1, 16).astype('int32')
    Y = run_code(per_elmt_imm, X, (2, 16), 'int32')
    assert all(X[0] + PER_ELMT_UNSIGNED_VALUES == Y[0])
    assert all(X[0] + PER_ELMT_SIGNED_VALUES == Y[1])

#============================== Vector rotatoin ===============================

@qpu
def vector_rotation(asm):
    mov(r0, vpm)
    nop()
    for shift in range(1, 16):
        rotate(vpm, r0, shift)
    for shift in range(-16, 16):
        ldi(broadcast, shift)
        nop()
        rotate(vpm, r0, r5)

def test_vector_rotation():
    X = np.array(
            [getrandbits(32) for i in range(16)]
            ).reshape(1, 16).astype('uint32')
    Y = run_code(vector_rotation, X, (15 + 32, 16), 'uint32')

    for i, shift in enumerate(range(1, 16)):
        assert all(np.roll(X[0], shift) == Y[i])
    for i, shift in enumerate(range(-16, 16)):
        assert all(np.roll(X[0], shift) == Y[i+15])

@qpu
def add_and_rotation(asm):
    mov(r0, vpm)
    mov(r1, vpm)
    rotate(r2, r0, -1)
    fadd(r0, r0, r2).rotate(r3, r1, -1)
    nop()
    fadd(r1, r1, r3).rotate(r2, r0, -2)
    nop()
    fadd(r0, r0, r2).rotate(r3, r1, -2)
    nop()
    fadd(r1, r1, r3).rotate(r2, r0, -4)
    nop()
    fadd(r0, r0, r2).rotate(r3, r1, -4)
    nop()
    fadd(r1, r1, r3).rotate(r2, r0, -8)
    nop()
    fadd(r0, r0, r2).rotate(r3, r1, -8)
    fadd(r1, r1, r3)

    mov(vpm, r0)
    mov(vpm, r1)

def test_add_and_rotation():
    X = np.random.randn(2, 16).astype('float32')
    Y = run_code(add_and_rotation, X, (2, 16), 'float32')
    assert np.allclose(Y[0], X[0].sum(), rtol=1e-4)
    assert np.allclose(Y[1], X[1].sum(), rtol=1e-4)

#=============================== Miscellaneous ================================

@qpu
def shared_small_imm(asm):
    mov(r0, vpm)
    fadd(r1, r0, 2.0).fmul(r2, r0, 2.0)
    mov(vpm, r1)
    mov(vpm, r2)

def test_shared_small_imm():
    X = np.random.randn(1, 16).astype('float32')
    Y = run_code(shared_small_imm, X, (2, 16), 'float32')
    assert np.allclose(X[0] + 2.0, Y[0], rtol=1e-6)
    assert np.allclose(X[0] * 2.0, Y[1], rtol=1e-6)

max_stars_repo_name: long-long-float/py-videocore

max_issues_repo_name: long-long-float/py-videocore

max_forks_repo_name: long-long-float/py-videocore


text:
import os
os.chdir("C:/Users/tprabhan/Documents/My_Books/RSADBE_2.0/R_and_Python_Programs/Python/Chapter_03")
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
MyPath/Python/Chapter_03/
'''
Bar Charts
'''
sc_data = pd.read_csv("Data/SC.csv",delimiter=',')
BT = sc_data.Bug_Type
BFreq = sc_data.Frequency
Ypos = np.arange(len(BT))
plt.barh(Ypos,BFreq,align='center',alpha=0.5)
plt.yticks(Ypos, BT)
plt.xlabel('Severity Counts')
plt.title('Bar graph for the Bug Metrics dataset')
plt.show()

'''
Dot Chart
'''
plt.plot(BFreq,Ypos,'o')
plt.yticks(Ypos, BT)
plt.xlabel('Severity Counts')
plt.title('Dot Chart for the Bug Metrics dataset')
plt.show()

'''
Boxplots
'''
resistivity = pd.read_csv("Data/Resistivity.csv",delimiter=',')
P1 = resistivity['Process.1']
P2 = resistivity['Process.2']
resistivity.boxplot(column=['Process.1','Process.2'])
plt.show()
resistivity.boxplot(column=['Process.1','Process.2'],notch=True)
plt.show()

morley = pd.read_csv("Data/Morley.csv",delimiter=',')
morley.boxplot(column='Speed',by='Expt')
plt.show()

'''
Histograms
'''
Galton = pd.read_csv("Data/Galton.csv",delimiter=',')
Galton.shape
Galton.head(10)
Galton.quantile([0,0.25,0.5,0.75,1])
Galton.mean()
Galton.median()

plt.subplot(1,2,1)
plt.hist(Galton['parent'],bins='fd')
plt.xlim(60,75)
plt.title("Height of Parent")
plt.subplot(1,2,2)
plt.hist(Galton['child'],bins='fd')
plt.xlim(60,75)
plt.title("Height of Child")
plt.show()

plt.subplot(1,2,1)
plt.hist(Galton['parent'],bins='sturges')
plt.xlim(60,75)
plt.title("Height of Parent")
plt.subplot(1,2,2)
plt.hist(Galton['child'],bins='sturges')
plt.xlim(60,75)
plt.title("Height of Child")
plt.show()

Octane = pd.read_csv("Data/Octane.csv",delimiter=',')
plt.subplot(2,2,1)
plt.hist(Octane['Method_1'],color="magenta")
plt.xlabel("Ratings Under Method I")
plt.title("Histogram of Octane Ratings for Method I")
plt.subplot(2,2,2)
plt.hist(Octane['Method_2'],color="cyan")
plt.xlabel("Ratings Under Method II")
plt.title("Histogram of Octane Ratings for Method II")
Samplez = pd.read_csv("Data/Samplez.csv",delimiter=',')
plt.subplot(2,2,3)
plt.hist(Samplez['Sample_1'],color="magenta")
plt.xlabel("Sample 1")
plt.title("Histogram of Sample 1")
plt.subplot(2,2,4)
plt.hist(Samplez['Sample_2'],color="cyan")
plt.xlabel("Sample 2")
plt.title("Histogram of Sample 2")
plt.show()

'''
Matrix of Scatter Plots
'''
Gasoline = pd.read_csv("Data/Gasoline.csv",delimiter=',')
Gasoline2 = Gasoline
del Gasoline2['x11']
from pandas.tools.plotting import scatter_matrix
scatter_matrix(Gasoline2)
plt.show()


text:
#!/usr/bin/env python
import sys
import numpy as np
import telescope_1d
import os
ndishes = int(sys.argv[1])
npix = int(sys.argv[2])
redundant = bool(sys.argv[3]=="1")
redstr = 'red' if redundant else 'nred'

t = None

for seed in range(30):
    for sigt in 'sig point gauss unif'.split():
        sig = None
        for te in [0,0.1,1.,10.,100.]:
            outfname = f"out/{ndishes}_{npix}_{redstr}_{seed}_{sigt}_{te}.npy"
            if os.path.isfile (outfname):
                print (f"{outfname} exists.")
                continue
            if t is None:
                if ndishes<16:
                    Nfreq = 512
                elif ndishes <32:
                    Nfreq = 512
                else:
                    Nfreq = 1024

                t = telescope_1d.Telescope1D(Ndishes=ndishes, Npix_fft=npix, redundant=redundant, Nfreq=Nfreq, seed=22)
            if sig is None:
                if sigt == 'sig':
                    sig = t.get_signal(seed=seed)
                elif sigt == 'point':
                    sig =  t.get_point_source_sky(seed=seed)
                elif sigt =='gauss':
                    sig = t.get_gaussian_sky(seed=seed)
                elif sigt =='unif':
                    sig = t.get_uniform_sky(seed=seed)
                else:
                    print ("Shit!")
                    stop
                uvsig = t.observe_image(sig)

            print (f"Working {outfname}", sig.sum())
            uvplane, uvplane_f, uvplane_1 = t.get_obs_uvplane(uvsig, time_error_sigma=te*1e-12, filter_FG=True)

            np.save(outfname,(uvplane, uvplane_f, uvplane_1))


text:
# %% [code]
import pandas as pd
import numpy as np

def dist(x1,x2):
    return np.sqrt(sum((x1-x2)**2)) # calculating distance

# main algo
def knn(X,Y,queryPoint,k=5):

    vals = [] # creating list to append all distances
    m = X.shape[0]

    for i in range(m):
        d = dist(queryPoint,X[i])
        vals.append((d,Y[i])) #appending all distances

    #sorting the list
    vals = sorted(vals)
    # choose first k distances
    vals = vals[:k]

    vals = np.array(vals)

    new_vals = np.unique(vals[:,1],return_counts=True)

    index = new_vals[1].argmax()
    pred = new_vals[0][index]

    return pred

## For testing Purposes
'''
## Importing libraries

import sklearn.datasets
import matplotlib.pyplot as plt

## creating dataset

x,y = sklearn.datasets.make_classification(n_samples=1000, n_classes=2,
n_clusters_per_class=1, n_features=2,n_informative=2, n_redundant=0, n_repeated=0)

## Visualization

query_p = np.array([0.5,0.5])
plt.scatter(query_p[0],query_p[1],c = 'r') ## plot the query point
plt.scatter(x[:,0],x[:,1],c = y)
plt.show()

## testing the algorithm

result = knn(x,y,query_p)    ### query point ==> x = 0.5,y = 0.5
print(result)
'''

max_stars_repo_path: K Nearest Neighbors/k-nearest neighbors (KNN).py

max_issues_repo_name: imraghavagr/Algo-ScriptML

max_forks_repo_name: imraghavagr/Algo-ScriptML


text:
"""

Code to generate large grids of cloudy models under different assumptions

Written by Michele Fumagalli in Durham, Summer 2015
michele.fumagalli@durham.ac.uk

Copyright (C) 2015 Michele Fumagalli

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/>.

"""

#general import
import numpy as np
from astropy import constants as const
from scipy import interpolate
# import mypython as mp ##comment out cbw
from scipy import integrate
from scipy.interpolate import interp1d
import os
import re

from pkg_resources import resource_filename

class Cldyion():

    """
    Define a class to run ionization modeling to do cloudy grids of elements
    By default, use the HM2005. Other UVB that can be included are: HM2012, FG09

    """

    def __init__(self,uvb='HM05'):

        #This will be constants
        self.column=False
        self.metal=False
        self.redshift=False
        self.density=False
        self.jgal=False
        self.jqso=False
        self.temperature=False
        self.fstar=False
        self.dustscale=False
        self.root=''
        self.path=''
        self.uvbtype=uvb
        self.HM05write=''
        self.subdir=''

        #this is stuff in output of a model (dictionaries)
        self.ins={}  #mean ion fraction
        self.clm={}  #column densities
        self.prof={} #radial profiles
        self.abund={} #abundances
        self.temperature_out=False #abundances

        #A dictionary that stores values of frequency log [Hz] and mean intensity (log 4 pi J_nu [erg/s/Hz/cm^2])
        self.source={'lgnu':0,'lgfpiJ':0}

        #A dictionary that stores correction fraction compared to the
        #Asplund 2009 abundances (1 is default, <1 is depleted)
        #Runs over the elements one wants to deplete
        self.abundance=False

    def minimal(self,column,metal,redshift,density):

        """
        Specify a minimal cloudy model (see Fumagalli et al. 2015)

        column     -> the HI column density [log]
        metal      -> the log metallicity in solar unit [X/H]
        redshift   -> the redshift
        density    -> the total hydrogen volume density [log]

        """

        #Start by defining defaults
        self.column=column
        self.metal=metal
        self.redshift=redshift
        self.density=density
        self.root="n{0}_r{1}_c{2}_Z{3}".format(self.density,self.redshift,self.column,self.metal)

        #leave elements as they are in minimal model
        #Initialize the HM05 UVB at the given redshift
        self.inituvb()

    def cie(self,column,metal,redshift,density,temperature):

        """
        Specify a cloudy model with constant temperature [CIE] (see Fumagalli et al. 2015)

        column     -> the HI column density [log]
        metal      -> the log metallicity in solar unit [X/H]
        redshift   -> the redshift
        density    -> the total hydrogen volume density [log]
        temperature -> the log of the temperature in K

        """

        #first init the minimal model
        self.minimal(column,metal,redshift,density)

        #update as needed
        self.temperature=temperature
        self.root="n{0}_r{1}_c{2}_Z{3}_T{4}".format(self.density,self.redshift,self.column,self.metal,self.temperature)

    def dust(self,column,metal,redshift,density,fstar):

        """
        Specify a cloudy model with dust depletion (see Fumagalli et al. 2015)

        column     -> the HI column density [log]
        metal      -> the log metallicity in solar unit [X/H]
        redshift   -> the redshift
        density    -> the total hydrogen volume density [log]
        fstar      -> the strength of the dust depletion [linear]

        """

        #first init the minimal model
        self.minimal(column,metal,redshift,density)

        #update as needed global variable
        self.fstar=fstar
        self.root="n{0}_r{1}_c{2}_Z{3}_f{4}".format(self.density,self.redshift,self.column,self.metal,self.fstar)

        #now compute the dust scaling
        if(self.fstar > 0):
            #If fstar > 0 ,add dust
            self.dustscale=self.fstar*(10**self.metal)
        else:
            #Else, do not include dust
            self.dustscale=False

        #now compute the abundances depletion following the method by Jenkins et al.
        #Going to f = -1 brings the corrections close to 0
        ejelm=['Carbon','Nitrogen','Oxygen','Magnesium','Silicon','Phosphorus','Chlorine','Titanium',\
                        'Chromium','Manganese','Iron','Nickel','Copper','Zinc']
        ejAx=[-0.101,-0.000,-0.225,-0.997,-1.136,-0.945,-1.242,-2.048,-1.447,-0.857,-1.285,-1.490,-0.710,-0.610]
        ejBx=[-0.193,-0.109,-0.145,-0.800,-0.570,-0.166,-0.314,-1.957,-1.508,-1.354,-1.513,-1.829,-1.102,-0.279]
        ejzx=[0.803,0.550,0.598,0.531,0.305,0.488,0.609,0.430,0.470,0.520,0.437,0.599,0.711,0.555]

        #init depletion
        self.abundance={}

        #compute values
        for ee in range(len(ejelm)):
            depl=ejBx[ee]+ejAx[ee]*(self.fstar-ejzx[ee])
            if(depl < 0):
                self.abundance[ejelm[ee]]=10**depl
            else:
                self.abundance[ejelm[ee]]=1.0

    def modsource(self,column,metal,redshift,density,jgal,jqso):

        """
        Specify a cloudy model with varying SED properties

        column     -> the HI column density [log]
        metal      -> the log metallicity in solar unit [X/H]
        redshift   -> the redshift
        density    -> the total hydrogen volume density [log]
        jgal       -> weight to add in the galaxy contribution
        jqso       -> weight to add in the qso contribution

        """

        #import numpy as np
        #from astropy import constants as const
        #from scipy import interpolate
        #import mypython as mp

        #first init the minimal model
        self.minimal(column,metal,redshift,density)

        #make copy of uvb
        uvb_nu=self.source['lgnu']
        uvb_f=self.source['lgfpiJ']

        #next update coefficient and name
        self.jgal=jgal
        self.jqso=jqso
        self.root="n{0}_r{1}_c{2}_Z{3}_g{4}_q{5}".format(self.density,self.redshift,self.column,self.metal,self.jgal,self.jqso)

        ##################
        ###  quasar    ###
        ##################

        #now load qso sed
        qso_nu=[]
        qso_f=[]

        richardsqsosed = resource_filename("pyigm", "/data/euvb/richards_qsosed.txt")
        fil=open(richardsqsosed)

        for ll in fil:
            fields=ll.split()
            try:
                qso_nu.append(float(fields[0])) #log Hz
                qso_f.append(float(fields[1]))  #log erg/s
            except:
                #skip header
                hhd=0

        fil.close()

        #extend to max frequency for UVB with alpha=-1
        qso_nu.append(np.max(self.source['lgnu']))
        qso_f.append(qso_f[-1])

        #go to np arrays
        qso_f=np.array(qso_f)
        qso_nu=np.array(qso_nu)

        #make qso sed comparable to the UVB in units
        #Cuba is in log of Hz; log of erg/s/cm2/Hz

        #Assume inner radius of 100 kpc for flux conversion
        inrad=1e5 #pc
        qso_f=qso_f-qso_nu-np.log10(4*np.pi)-2*np.log10(inrad*const.pc.to('cm').value)

        #Now apply scaling
        qso_f=qso_f+self.jqso

        #add to model the quasar contrinution
        qsoint=interpolate.interp1d(qso_nu,qso_f,fill_value=-99,bounds_error=False)
        qso_f_int=qsoint(self.source['lgnu'])
        self.source['lgfpiJ']=np.log10(10**self.source['lgfpiJ']+10**qso_f_int)

        ##################
        ###  galaxy    ###
        ##################
        # CBW # gal=mp.seds.sb99.Sb99()
        gal.loadmod('sfr1_ge_z014',path='./Scripts/')

        #grab the last time
        gal_lam=np.array(gal.spectrum[-1]['WAVE']) #A
        gal_f=np.array(gal.spectrum[-1]['TOTAL']) #log erg/s/A

        #convert in unit similar to qso
        gal_nu=np.log10(const.c.to('cm/s').value*1e8/gal_lam) #log Hz
        gal_f=np.log10(10**gal_f*gal_lam/10**gal_nu) #log erg/s/Hz

        #go to flux at given distance
        gal_f=gal_f-np.log10(4*np.pi)-2*np.log10(inrad*const.pc.to('cm').value) #log erg/s/Hz/cm2

        #apply scaling
        gal_f=gal_f+self.jgal

        #add to model
        galint=interpolate.interp1d(gal_nu,gal_f,fill_value=-99,bounds_error=False)
        gal_f_int=galint(self.source['lgnu'])
        self.source['lgfpiJ']=np.log10(10**self.source['lgfpiJ']+10**gal_f_int)

        ##DEBUG sequence
        #import matplotlib as mpl
        #import matplotlib.pyplot as pl
        #pl.plot(uvb_nu,uvb_f,linestyle=':', linewidth=2)
        #pl.plot(self.source['lgnu'],qso_f_int,linestyle=':', linewidth=2)
        #pl.plot(self.source['lgnu'],gal_f_int,linestyle=':', linewidth=2)
        #pl.plot(self.source['lgnu'],self.source['lgfpiJ'], linewidth=4)
        #pl.show()
        #exit()

        #pl.plot(self.source['lgnu'],10**(self.source['lgfpiJ']-uvb_f))
        #pl.show()

        #import pickle
        #src={'lgnu':self.source['lgnu'],'all':self.source['lgfpiJ'],'uvb':uvb_f,'gal':gal_f_int,'qso':qso_f_int}
        #fil=open('source_debug.pkl','w')
        #pickle.dump(src,fil)
        #fil.close()
        ## END of DEBUG sequence

    def writein(self,basedir):

        """
        Write the init file

        """
        #from scipy import integrate
        #from astropy import constants as const
        #from scipy.interpolate import interp1d
        #import numpy as np

        ##CBW separate into sub-directories
        self.subdir="n{0}/r{1}/".format(self.density,self.redshift)
        self.path=basedir+"/"+self.subdir
        if not os.path.isdir(self.path+"/"):
            os.makedirs(self.path+"/")

        ##CBW don't write to stdout every time
        # print('Write {}'.format(self.root))

        fil=open(self.path+"/"+self.root+".in","w")

        #Write header
        fil.write("##Header\n")
        fil.write("# Metallicity {0}\n".format(self.metal))
        fil.write("# Density {0}\n".format(self.density))
        fil.write("# Redshift {0}\n".format(self.redshift))
        fil.write("# Column {0}\n".format(self.column))
        fil.write("# Jgal {0}\n".format(self.jgal))
        fil.write("# Jqso {0}\n".format(self.jqso))
        fil.write("# Temp {0}\n".format(self.temperature))
        fil.write("# Fstar {0}\n".format(self.fstar))

        #set the metallicity
        fil.write("##Abundances\n")
        #turn on default
        fil.write("abundances GASS10 no grains\n")
        #scale all elements
        fil.write("metals {0} log \n".format(self.metal))
        #loop over elemnts for individual corrections
        if(self.abundance):
            for dep in self.abundance:
                fil.write("element scale factor {0} {1}\n".format(dep,self.abundance[dep]))

        #if needed specify dust
        if(self.dustscale):
            fil.write("grains ism {}\n".format(self.dustscale))
            fil.write("no grain qheat\n")

        #density/geometry
        fil.write("##Geometry/density\n")
        fil.write("hden {0}\n".format(self.density))

        #Include constant temperature for CIE
        if(self.temperature):
            fil.write("constant temperature {0} log\n".format(self.temperature))

        #radiation
        fil.write("##Radiation\n")
        fil.write("cmb, z={0}\n".format(self.redshift))

        #if using HM05, write the call to the spectrum here
        ##else write the continuum by hand
        if(self.HM05write != ""):
            fil.write(self.HM05write)
        else:
            #compute integral over range free of difficulties
            #Define min-max in Hz
            maxnu=(const.c.to('cm/s')/(1250*1e-8)).value
            ionnu=(const.c.to('cm/s')/(2000*1e-8)).value

            #make sure interval is in range of model, and convert to Ry
            ionnu = np.max((ionnu,10**np.min(self.source['lgnu'])))
            maxry=(4.1356675e-15*maxnu/13.605698066)
            minry=(4.1356675e-15*ionnu/13.605698066)

            #integrate
            fint = interp1d(10**self.source['lgnu'],10**self.source['lgfpiJ'])
            integral = np.log10(integrate.quad(fint,ionnu,maxnu))

            #now write
            self.writecontinuum(fil)
            fil.write("intensity {0}, range {1} to {2} Ryd\n".format(integral[0],minry,maxry))

        #stop/save
        fil.write("##StopSave\n")
        fil.write("stop neutral column density {0}\n".format(self.column))
        #force no stop temperature
        fil.write("stop temperature off\n")
        #Add stop for CIE
        if(self.temperature):
            fil.write("stop column density 24\n")
        fil.write("iterate to convergence\n")
        fil.write('set save prefix "{0}"\n'.format(self.root))
        #fil.write('save last continuum units Angstrom ".con"\n')
        fil.write('save last column densities ".col"\n')
        fil.write('save last ionization means ".ion"\n')
        #fil.write('save last overview ".ovr"\n')
        fil.write('save last abundances ".abn"\n')
        #fil.write('save last radius ".rad"\n')

        fil.close()

    def inituvb(self):

        """
        This function is a dispatcher of the different UVB models
        currently supported by the code

        """

        #Initialize the HM12 UVB at the given redshift
        if(self.uvbtype == 'HM12'):
            self.initcuba()

        #Initialize the HM05 UVB by printing the
        #  command to the Cloudy *.in file
        elif(self.uvbtype == 'HM05'):
            self.initHM05()

        #Initialize the HM05 UVB that creates the UVB at
        #  a given redshift, used with getionisation.py
        elif(self.uvbtype == 'HM05logU'):
            self.initHM05logU()

        #Initialize the FC09 UVB at the given redshift
        elif(self.uvbtype == 'FG09'):
            self.initcafg()

        else:
            print('Unknown UVB!')
            exit()

    def initcuba(self):

        """
        Initialise the CUBA HM12 UVB at a given redshift

        """

        #import numpy as np
        #from astropy import constants as const

        #load cuba
        input_spectrum = resource_filename("pyigm", "/data/euvb/UVB_HM12.dat")
        uvb=open(input_spectrum)
        flag=0

        #store the input cuba
        cuba_red=[]
        cuba_wave=[]

        #pick the flux at three redshifts near to the z of interest
        cuba_one=[]
        cuba_two=[]
        cuba_three=[]

        #loop cuba file
        for line in uvb:
            if not "#" in line:

                #parse cuba redshift first
                if(flag == 0):
                    line=line.strip()
                    line=line.split()
                    for l in line:
                        cuba_red.append(float(l))
                    cuba_red=np.array(cuba_red)
                    flag = flag+1

                #now read the left and right spectrum for input redshift
                else:

                    #wave
                    line=line.strip()
                    line=line.split()
                    cuba_wave.append(float(line[0]))

                    #check edge of the grid
                    index=np.argmin(abs(self.redshift - cuba_red))

                    #if edge of the grid - stay at that redshift
                    if((index < 1) | (index > cuba_red.size-2)):
                        cuba_one.append(float(line[index+1]))
                        cuba_two.append(float(line[index+1]))
                        cuba_three.append(float(line[index+1]))
                        cuba_int_red=[cuba_red[index],cuba_red[index],cuba_red[index]]

                    else:
                        #prepare for interpolation by loading z +/- 1 setp
                        cuba_one.append(float(line[index]))
                        cuba_two.append(float(line[index+1]))
                        cuba_three.append(float(line[index+2]))
                        cuba_int_red=[cuba_red[index-1],cuba_red[index],cuba_red[index+1]]

        uvb.close()

        #turn to arrays - interpolate in log space
        cuba_wave=np.array(cuba_wave)
        cuba_one=np.log10(np.array(cuba_one))
        cuba_two=np.log10(np.array(cuba_two))
        cuba_three=np.log10(np.array(cuba_three))
        cuba_int_red=np.array(cuba_int_red)
        cuba_int=np.zeros(cuba_one.size)

        #interpolate to give redshift
        for ii in range(cuba_wave.size):
            cuba_int[ii]=np.interp([self.redshift],cuba_int_red,[cuba_one[ii],cuba_two[ii],cuba_three[ii]])
            #print(cuba_wave[ii], cuba_int[ii])

        #cuba has the same wave twice at discontinuities.
        #cloudy does not like that, so perturbe these values
        for ii in range(cuba_wave.size-1):
            if(cuba_wave[ii]==cuba_wave[ii+1]):
                cuba_wave[ii]=cuba_wave[ii]-0.005

        #import matplotlib.pyplot as plt
        #cuba_wave=np.log10(cuba_wave)
        #plt.plot(cuba_wave,cuba_one,color='red')
        #plt.plot(cuba_wave,cuba_two,color='green')
        #plt.plot(cuba_wave,cuba_three,color='blue')
        #plt.plot(cuba_wave,cuba_int,color='black')
        #plt.show()

        #store the interpolated UVB
        lgnu=const.c.to('cm/s')/(cuba_wave*1e-8)
        lgnu=np.log10(lgnu.value) #Log of Hz
        cuba_int=cuba_int+np.log10(4*np.pi) #get rid of /str and go to Log of erg/s/cm2/Hz
        self.source['lgnu']=lgnu
        self.source['lgfpiJ']=cuba_int

    def initHM05(self):

        self.HM05write="table HM05 redshift {0}\n".format(self.redshift)
        self.HM05write="table HM05, z={0}\n".format(self.redshift)

    def initHM05logU(self):

        """
        Initialise the HM05 UVB at a given redshift

        """

        #import numpy as np
        #from astropy import constants as const

        #load hm05cbw
        input_spectrum = resource_filename("pyigm", "/data/euvb/UVB_HM05.dat")
        uvb=open(input_spectrum)
        flag=0

        #store the input hm05cbw
        hm05cbw_red=[]
        hm05cbw_wave=[]

        #pick the flux at three redshifts near to the z of interest
        hm05cbw_one=[]
        hm05cbw_two=[]
        hm05cbw_three=[]

        #loop hm05cbw file
        for line in uvb:
            if not "#" in line:

                #parse hm05cbw redshift first
                if(flag == 0):
                    line=line.strip()
                    line=line.split()
                    for l in line:
                        hm05cbw_red.append(float(l))
                    hm05cbw_red=np.array(hm05cbw_red)
                    flag = flag+1

                #now read the left and right spectrum for input redshift
                else:

                    #wave
                    line=line.strip()
                    line=line.split()
                    hm05cbw_wave.append(float(line[0]))

                    #check edge of the grid
                    index=np.argmin(abs(self.redshift - hm05cbw_red))

                    #if edge of the grid - stay at that redshift
                    if((index < 1) | (index > hm05cbw_red.size-2)):
                        hm05cbw_one.append(float(line[index+1]))
                        hm05cbw_two.append(float(line[index+1]))
                        hm05cbw_three.append(float(line[index+1]))
                        hm05cbw_int_red=[hm05cbw_red[index],hm05cbw_red[index],hm05cbw_red[index]]

                    else:
                        #prepare for interpolation by loading z +/- 1 setp
                        hm05cbw_one.append(float(line[index]))
                        hm05cbw_two.append(float(line[index+1]))
                        hm05cbw_three.append(float(line[index+2]))
                        hm05cbw_int_red=[hm05cbw_red[index-1],hm05cbw_red[index],hm05cbw_red[index+1]]

        uvb.close()

        #turn to arrays - interpolate in log space
        hm05cbw_wave=np.array(hm05cbw_wave)
        hm05cbw_one=np.log10(np.array(hm05cbw_one))
        hm05cbw_two=np.log10(np.array(hm05cbw_two))
        hm05cbw_three=np.log10(np.array(hm05cbw_three))
        hm05cbw_int_red=np.array(hm05cbw_int_red)
        hm05cbw_int=np.zeros(hm05cbw_one.size)

        #interpolate to give redshift
        for ii in range(hm05cbw_wave.size):
            hm05cbw_int[ii]=np.interp([self.redshift],hm05cbw_int_red,[hm05cbw_one[ii],hm05cbw_two[ii],hm05cbw_three[ii]])
            #print(hm05cbw_wave[ii], hm05cbw_int[ii])

        #hm05cbw has the same wave twice at discontinuities.
        #cloudy does not like that, so perturbe these values
        for ii in range(hm05cbw_wave.size-1):
            if(hm05cbw_wave[ii]==hm05cbw_wave[ii+1]):
                hm05cbw_wave[ii]=hm05cbw_wave[ii]-0.005

        #import matplotlib.pyplot as plt
        #hm05cbw_wave=np.log10(hm05cbw_wave)
        #plt.plot(hm05cbw_wave,hm05cbw_one,color='red')
        #plt.plot(hm05cbw_wave,hm05cbw_two,color='green')
        #plt.plot(hm05cbw_wave,hm05cbw_three,color='blue')
        #plt.plot(hm05cbw_wave,hm05cbw_int,color='black')
        #plt.show()

        #store the interpolated UVB
        lgnu=const.c.to('cm/s')/(hm05cbw_wave*1e-8)
        lgnu=np.log10(lgnu.value) #Log of Hz
        hm05cbw_int=hm05cbw_int+np.log10(4*np.pi) #get rid of /str and go to Log of erg/s/cm2/Hz
        self.source['lgnu']=lgnu
        self.source['lgfpiJ']=hm05cbw_int

    def initcafg(self):

        """
        Initialise the FG09 UVB at a given redshift

        """

        #all available redshift for UVB
        avail_z=np.arange(0,10,0.05)

        #pick closest in redshift
        zindex=np.argmin(abs(avail_z-self.redshift))
        nameuvb='fg_uvb_dec11/fg_uvb_dec11_z_{}.dat'.format(avail_z[zindex])
        try:
            fgtab=np.loadtxt(nameuvb,dtype={'names': ('nu','f'),'formats':('f10','f10')})
        except:
            print('Could not read UVB file {}'.format(nameuvb))
            exit()

        #from Ry to Hz
        self.source['lgnu']=np.log10(fgtab['nu'])+np.log10(3.2898419499e15)

        #from (10^-21 erg s^-1 cm^-2 Hz^-1 sr^-1) to Log of erg/s/cm2/Hz
        self.source['lgfpiJ']=np.log10(fgtab['f'])-21+np.log10(4*np.pi)

        #now reverse freq to have high freq first and low freq last (as in HM12)
        #This is nothing special, but it is needed for how the writing procedure
        #fixes the boundary of the SED
        self.source['lgnu'] = self.source['lgnu'][::-1]
        self.source['lgfpiJ'] = self.source['lgfpiJ'][::-1]

    def writecontinuum(self,fil):

        """
        Parse the continuum for an input cloudy file

        """

        #import numpy as np

        #make copies
        spec = 10**np.array(self.source['lgfpiJ'])
        specclean = np.copy(spec)
        logfreq=np.copy(self.source['lgnu'])

        #first point
        specclean[spec == 0.0] = np.amin(spec[spec > 0])*1e-6
        logL_nu = np.log10(specclean)

        #put a floor at -40 when writing to avoid cloudy complaints
        fil.write("interpolate")
        fil.write(" ({0:f} {1:f})".format(7.51,
                                          max(np.amin(logL_nu)-6,-40)))
        fil.write(" ({0:f} {1:f})".format(logfreq[-1]-0.01,
                                          max(np.amin(logL_nu)-6,-40)))
        #stuff in between
        for i in range(len(logL_nu)):
            fil.write("\ncontinue")
            fil.write(" ({0:f} {1:f})".format(logfreq[-i-1],
                                              max(logL_nu[-i-1],-40)))
        #last point
        fil.write("\ncontinue ({0:f} {1:f})".
                  format(logfreq[0]+0.01, max(np.amin(logL_nu)-6,-40)))

        fil.write(" ({0:f} {1:f})\n".format(22.4,max(np.amin(logL_nu)-6,-40)))

    def checkrun(self,stop=False):

        """
        Check if a model has run successfully
        Check stop check if the stop was due to lower column

        """

        #import os

        #check if file exists
        status = False
        warning = True
        filname=self.path+'/'+self.root+'.out'
        filthere=os.path.isfile(filname)

        #check exist status
        linestop=''
        if filthere:
            status=True
            #print('Parse output for warnings')
            file1=open(filname,"r")
            #loop entire file
            for line in file1:
                #check exit condition
                if ("Cloudy exited OK" in line):
                    warning = False
                #grab stopping line
                if "Calculation stopped because" in line:
                    linestop=line
            file1.close()

            #if exit condition not clean, raise warning
            if warning:
                print("Warning: cloudy issues for {}".format(self.root))

            #check if stopped because of unwanted reasons:
            if stop:
                if ("H0 column dens reached" not in linestop) and ("H0-H0+H2 column dens reached" not in linestop):
                    print("Stop: condition not met in {0}".format(self.root))

        return status

    def getion(self):

        """
        Get the ionization fractions from a model

        """
        #import re
        #import numpy as np

        #first check if everything is fine
        status=self.checkrun()

        #if so proceed
        if status:

            #open safe
            filname=self.path+'/'+self.root+'.ion'
            try:
                file1=open(filname,"r")
            except:
                print("Ion not found")
                return

            #now parse
            for line in file1:

                #check if a good data
                if '-' in line:

                    #split it
                    fields=re.split(' |-|\n',line)
                    #remove unwanted and make it array numbers
                    fields=filter(None, fields)
                    values=np.array([-1*float(x) for x in fields[1:] if '.' in x])
                    self.ins[fields[0]]=values

            #close and return
            file1.close()

        else:
            print("Check model... Something is wrong!")
            return

    def getcolumn(self):

        """
        Get the column density from a model

        """

        #import re
        #import numpy as np

        #first check if everything is fine
        status=self.checkrun()

        #if so proceed
        if status:

            #open safe
            filname=self.path+'/'+self.root+'.col'
            try:
                file1=open(filname,"r")
            except:
                print("Column not found")
                return

            #now parse
            header=file1.readline()
            value=file1.readline()

            #extract
            clmd=np.array([float(x) for x in value.split()])
            header=np.array((header.split())[2:])

            #parity check
            if (clmd.size != header.size):
                print("Something wrong in parser column")
                return

            for ii in range(header.size):
                if(clmd[ii] > 0):
                    self.clm[header[ii]]=np.log10(clmd[ii])

            #close and return
            file1.close()

        else:
            print("Check model... Something is wrong!")
            return

    def getprofile(self):

        """
        Get the radial overview from a model

        """

        #import re
        #import numpy as np

        #first check if everything is fine
        status=self.checkrun()

        #if so proceed
        if status:

            #open safe
            filname=self.path+'/'+self.root+'.ovr'
            try:
                file1=open(filname,"r")
            except:
                print("Overview not found")
                return

            #parse header
            header=file1.readline()
            header=header.replace("#","")
            header=header.split()

            #init dictionary
            for ff in header:
                self.prof[ff]=np.empty(0)

            #now parse
            for line in file1:
                line=line.split()

                ii=0
                for ff in header:
                    self.prof[ff]=np.append(self.prof[ff],float(line[ii]))
                    ii=ii+1

            #close and return
            file1.close()

        else:
            print("Check model... Something is wrong!")
            return

    def getabundance(self):

        """
        Get the gas phase abundance in each zone

        """

        #import re
        #import numpy as np

        #first check if everything is fine
        status=self.checkrun()

        #if so proceed
        if status:

            #open safe
            filname=self.path+'/'+self.root+'.abn'
            try:
                file1=open(filname,"r")
            except:
                print("Abundance not found")
                return

            #parse header
            header=file1.readline()
            header=header.replace("#","")
            header=header.split()
            header=header[1:]

            #init dictionary
            for ff in header:
                self.abund[ff]=np.empty(0)

            #now parse
            for line in file1:
                line=line.split()

                ii=0
                for ff in header:
                    self.abund[ff]=np.append(self.abund[ff],float(line[ii]))
                    ii=ii+1

            #close and return
            file1.close()

        else:
            print("Check model... Something is wrong!")
            return

    def gettemperature(self):

        """
        Get the overall/average physical temperature

        """

        ##CBW wrote this

        #import re
        #import numpy as np

        #first check if everything is fine
        status=self.checkrun()

        #if so proceed
        if status:

            #open safe
            filname=self.path+'/'+self.root+'.out'
            try:
                file1=open(filname,"r")
            except:
                print("Temperature not found")
                return

            search_string = "Te      Te(Ne)   Te(NeNp)"

            results = []
            for line in file1:
                if search_string in line:
                    results.append(file1.next())

            res = results[-1]
            self.temperature_out = np.log10(float(res.split()[1]))

            #close and return
            file1.close()

        else:
            print("Check model... Something is wrong!")
            return

max_stars_repo_name: pyigm/pyigm

max_issues_repo_name: pyigm/pyigm

max_forks_repo_name: pyigm/pyigm

reason: import numpy,from scipy,from astropy


text:
import argparse
import os
import shutil
import sys
import time
import warnings
from random import sample

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn import metrics
from torch.autograd import Variable
from torch.optim.lr_scheduler import MultiStepLR

from ogcnn.data import CIFData
from ogcnn.data import collate_pool, get_train_val_test_loader
from ogcnn.model import OrbitalCrystalGraphConvNet

parser = argparse.ArgumentParser(description='Orbital Graph Convolutional Neural Networks')
parser.add_argument('data_options', metavar='OPTIONS', nargs='+',
                    help='dataset options, started with the path to root dir, '
                         'then other options')
parser.add_argument('--task', choices=['regression', 'classification'],
                    default='regression', help='complete a regression or '
                                                   'classification task (default: regression)')
parser.add_argument('--disable-cuda', action='store_true',
                    help='Disable CUDA')
parser.add_argument('-j', '--workers', default=0, type=int, metavar='N',
                    help='number of data loading workers (default: 0)')
parser.add_argument('--epochs', default=100, type=int, metavar='N',
                    help='number of total epochs to run (default: 30)')
parser.add_argument('--start-epoch', default=0, type=int, metavar='N',
                    help='manual epoch number (useful on restarts)')
parser.add_argument('-b', '--batch-size', default=32, type=int,
                    metavar='N', help='mini-batch size (default: 256)')
parser.add_argument('--lr', '--learning-rate', default=0.0005, type=float,
                    metavar='LR', help='initial learning rate (default: '
                                       '0.0005)')
parser.add_argument('--lr-milestones', default=[100], nargs='+', type=int,
                    metavar='N', help='milestones for scheduler (default: '
                                      '[100])')
parser.add_argument('--momentum', default=0.9, type=float, metavar='M',
                    help='momentum')
parser.add_argument('--weight-decay', '--wd', default=0, type=float,
                    metavar='W', help='weight decay (default: 0)')
parser.add_argument('--print-freq', '-p', default=50, type=int,
                    metavar='N', help='print frequency (default: 50)')
parser.add_argument('--resume', default='', type=str, metavar='PATH',
                    help='path to latest checkpoint (default: none)')
train_group = parser.add_mutually_exclusive_group()
train_group.add_argument('--train-ratio', default=None, type=float, metavar='N',
                    help='number of training data to be loaded (default none)')
train_group.add_argument('--train-size', default=None, type=int, metavar='N',
                         help='number of training data to be loaded (default none)')
valid_group = parser.add_mutually_exclusive_group()
valid_group.add_argument('--val-ratio', default=0.1, type=float, metavar='N',
                    help='percentage of validation data to be loaded (default '
                         '0.1)')
valid_group.add_argument('--val-size', default=None, type=int, metavar='N',
                         help='number of validation data to be loaded (default '
                              '1000)')
test_group = parser.add_mutually_exclusive_group()
test_group.add_argument('--test-ratio', default=0.1, type=float, metavar='N',
                    help='percentage of test data to be loaded (default 0.1)')
test_group.add_argument('--test-size', default=None, type=int, metavar='N',
                        help='number of test data to be loaded (default 1000)')

parser.add_argument('--optim', default='SGD', type=str, metavar='SGD',
                    help='choose an optimizer, SGD or Adam, (default: SGD)')
parser.add_argument('--atom-fea-len', default=1536, type=int, metavar='N',
                    help='number of hidden atom features in conv layers')
parser.add_argument('--hot-fea-len', default=768, type=int, metavar='N',
                    help='number of hidden atom features in decoder')
parser.add_argument('--h-fea-len', default=128, type=int, metavar='N',
                    help='number of hidden features after pooling')
parser.add_argument('--n-conv', default=3, type=int, metavar='N',
                    help='number of conv layers')
parser.add_argument('--n-h', default=1, type=int, metavar='N',
                    help='number of hidden layers after pooling')

args = parser.parse_args(sys.argv[1:])

args.cuda = not args.disable_cuda and torch.cuda.is_available()

if args.task == 'regression':
    best_mae_error = 1e10
else:
    best_mae_error = 0.

def main():
    global args, best_mae_error

    # load data
    dataset = CIFData(*args.data_options)
    collate_fn = collate_pool
    train_loader, val_loader, test_loader = get_train_val_test_loader(
        dataset=dataset,
        collate_fn=collate_fn,
        batch_size=args.batch_size,
        train_ratio=args.train_ratio,
        num_workers=args.workers,
        val_ratio=args.val_ratio,
        test_ratio=args.test_ratio,
        pin_memory=args.cuda,
        train_size=args.train_size,
        val_size=args.val_size,
        test_size=args.test_size,
        return_test=True)

    # obtain target value normalizer
    if args.task == 'classification':
        normalizer = Normalizer(torch.zeros(2))
        normalizer.load_state_dict({'mean': 0., 'std': 1.})
    else:
        if len(dataset) < 500:
            warnings.warn('Dataset has less than 500 data points. '
                          'Lower accuracy is expected. ')
            sample_data_list = [dataset[i] for i in range(len(dataset))]
        else:
            sample_data_list = [dataset[i] for i in
                                sample(range(len(dataset)), 500)]
        _, sample_target, _ = collate_pool(sample_data_list)
        normalizer = Normalizer(sample_target)

    # build model
    structures, _, _ = dataset[0]

    orig_atom_fea_len = structures[0][0].shape[-1]     #92 features per atom
    ofm_fea= int(structures[0][1].shape[0]/structures[0][0].shape[0])
    nbr_fea_len = structures[1].shape[-1]      # 41 features per atom neighbors'
    orig_hot_fea_len = structures[0][1].shape[1]*ofm_fea       # 1056 features per crystal
    atom_fea_len = args.atom_fea_len
    hot_fea_len = args.hot_fea_len
    n_conv = args.n_conv
    h_fea_len = args.h_fea_len
    orig_atom_fea_len = orig_atom_fea_len + orig_hot_fea_len

    model = OrbitalCrystalGraphConvNet(orig_atom_fea_len, nbr_fea_len,orig_hot_fea_len,
                                 atom_fea_len=atom_fea_len,
                                 hot_fea_len=hot_fea_len,
                                 n_conv=n_conv,
                                 h_fea_len=h_fea_len,
                                 n_h=args.n_h,
                                 classification=True if args.task ==
                                                        'classification' else False)

    if args.cuda:
        model.cuda()

    # define loss func and optimizer
    if args.task == 'classification':
        criterion = nn.NLLLoss()
    else:
        criterion = nn.MSELoss()
    if args.optim == 'SGD':
        optimizer = optim.SGD(model.parameters(), args.lr,
                              momentum=args.momentum,
                              weight_decay=args.weight_decay)
    elif args.optim == 'Adam':
        optimizer = optim.Adam(model.parameters(), args.lr,
                               weight_decay=args.weight_decay)
    else:
        raise NameError('Only SGD or Adam is allowed as --optim')

    # optionally resume from a checkpoint
    if args.resume:
        if os.path.isfile(args.resume):
            print("=> loading checkpoint '{}'".format(args.resume))
            checkpoint = torch.load(args.resume)
            args.start_epoch = checkpoint['epoch']
            best_mae_error = checkpoint['best_mae_error']
            model.load_state_dict(checkpoint['state_dict'])
            optimizer.load_state_dict(checkpoint['optimizer'])
            normalizer.load_state_dict(checkpoint['normalizer'])
            print("=> loaded checkpoint '{}' (epoch {})"
                  .format(args.resume, checkpoint['epoch']))
        else:
            print("=> no checkpoint found at '{}'".format(args.resume))

    scheduler = MultiStepLR(optimizer, milestones=args.lr_milestones,
                            gamma=0.1)

    for epoch in range(args.start_epoch, args.epochs):
        train(train_loader, model, criterion, optimizer, epoch, normalizer)
        mae_error = validate(val_loader, model, criterion, normalizer)

        if mae_error != mae_error:
            print('Exit due to NaN')
            sys.exit(1)

        scheduler.step()

        # remember the best mae_eror and save checkpoint
        if args.task == 'regression':
            is_best = mae_error < best_mae_error
            best_mae_error = min(mae_error, best_mae_error)
        else:
            is_best = mae_error > best_mae_error
            best_mae_error = max(mae_error, best_mae_error)
        save_checkpoint({
            'epoch': epoch + 1,
            'state_dict': model.state_dict(),
            'best_mae_error': best_mae_error,
            'optimizer': optimizer.state_dict(),
            'normalizer': normalizer.state_dict(),
            'args': vars(args)
        }, is_best)

    # test best model
    print('---------Evaluate Model on Test Set---------------')
    best_checkpoint = torch.load('model_best.pth.tar')
    model.load_state_dict(best_checkpoint['state_dict'])
    validate(test_loader, model, criterion, normalizer, test=True)

def train(train_loader, model, criterion, optimizer, epoch, normalizer):
    batch_time = AverageMeter()
    data_time = AverageMeter()
    losses = AverageMeter()
    if args.task == 'regression':
        mae_errors = AverageMeter()
    else:
        accuracies = AverageMeter()
        precisions = AverageMeter()
        recalls = AverageMeter()
        fscores = AverageMeter()
        auc_scores = AverageMeter()

    # switch to train mode
    model.train()
    epoch_loss = 0

    end = time.time()
    for i, (input, target, batch_cif_ids) in enumerate(train_loader):
        data_time.update(time.time() - end)

        if args.cuda:
            input_var = (Variable(torch.cat((input[0][0],input[0][1].reshape(int(input[0][1].shape[0]/32),1056)),dim=1).cuda(non_blocking=True)),
                         Variable(input[1].cuda(non_blocking=True)),
                         input[2].cuda(non_blocking=True),
                         [crys_idx.cuda(non_blocking=True) for crys_idx in input[3]])
        else:
            input_var = (Variable(torch.cat((input[0][0],input[0][1].reshape(int(input[0][1].shape[0]/32),1056)),dim=1).cuda(non_blocking=True)),
                         Variable(input[1]),
                         input[2],
                         input[3])
        # normalize target
        if args.task == 'regression':
            target_normed = normalizer.norm(target)
        else:
            target_normed = target.view(-1).long()
        if args.cuda:
            target_var = Variable(target_normed.cuda(non_blocking=True))
        else:
            target_var = Variable(target_normed)

        # compute output
        output = model(*input_var)
        loss = criterion(output, target_var)

        if args.task == 'regression':
            mae_error = mae(normalizer.denorm(output.data.cpu()), target)
            losses.update(loss.data.cpu(), target.size(0))
            mae_errors.update(mae_error, target.size(0))
        else:
            accuracy, precision, recall, fscore, auc_score = \
                class_eval(output.data.cpu(), target)
            losses.update(loss.data.cpu().item(), target.size(0))
            accuracies.update(accuracy, target.size(0))
            precisions.update(precision, target.size(0))
            recalls.update(recall, target.size(0))
            fscores.update(fscore, target.size(0))
            auc_scores.update(auc_score, target.size(0))

        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        batch_time.update(time.time() - end)
        end = time.time()

        if i % args.print_freq == 0:
            if args.task == 'regression':
            	print('done')
            	print('Epoch: [{0}][{1}/{2}]\t'
                      'Time {batch_time.val:.3f} ({batch_time.avg:.3f})\t'
                      'Data {data_time.val:.3f} ({data_time.avg:.3f})\t'
                      'Loss {loss.val:.4f} ({loss.avg:.4f})\t'
                      'MAE {mae_errors.val:.3f} ({mae_errors.avg:.3f})'.format(
                    epoch, i, len(train_loader), batch_time=batch_time,
                    data_time=data_time, loss=losses, mae_errors=mae_errors)
                )
            else:
                print('Epoch: [{0}][{1}/{2}]\t'
                      'Time {batch_time.val:.3f} ({batch_time.avg:.3f})\t'
                      'Data {data_time.val:.3f} ({data_time.avg:.3f})\t'
                      'Loss {loss.val:.4f} ({loss.avg:.4f})\t'
                      'Accu {accu.val:.3f} ({accu.avg:.3f})\t'
                      'Precision {prec.val:.3f} ({prec.avg:.3f})\t'
                      'Recall {recall.val:.3f} ({recall.avg:.3f})\t'
                      'F1 {f1.val:.3f} ({f1.avg:.3f})\t'
                      'AUC {auc.val:.3f} ({auc.avg:.3f})'.format(
                    epoch, i, len(train_loader), batch_time=batch_time,
                    data_time=data_time, loss=losses, accu=accuracies,
                    prec=precisions, recall=recalls, f1=fscores,
                    auc=auc_scores)
                )
    torch.cuda.empty_cache()
    del input
    del target
    del batch_cif_ids

def validate(val_loader, model, criterion, normalizer, test=False):
    batch_time = AverageMeter()
    losses = AverageMeter()
    if args.task == 'regression':
        mae_errors = AverageMeter()
    else:
        accuracies = AverageMeter()
        precisions = AverageMeter()
        recalls = AverageMeter()
        fscores = AverageMeter()
        auc_scores = AverageMeter()
    if test:
        test_targets = []
        test_preds = []
        test_cif_ids = []

    # switch to evaluate mode
    model.eval()

    end = time.time()

    epoch_loss = 0
    for i, (input, target, batch_cif_ids) in enumerate(val_loader):
        if args.cuda:
            with torch.no_grad():
                input_var = (Variable(torch.cat((input[0][0],input[0][1].reshape(int(input[0][1].shape[0]/32),1056)),dim=1).cuda(non_blocking=True)),
                             Variable(input[1].cuda(non_blocking=True)),
                             input[2].cuda(non_blocking=True),
                             [crys_idx.cuda(non_blocking=True) for crys_idx in input[3]])
        else:
            with torch.no_grad():
                input_var = (Variable(torch.cat((input[0][0],input[0][1].reshape(int(input[0][1].shape[0]/32),1056)),dim=1)),
                             Variable(input[1]),
                             input[2],
                             input[3])
        if args.task == 'regression':
            target_normed = normalizer.norm(target)
        else:
            target_normed = target.view(-1).long()
        if args.cuda:
            with torch.no_grad():
                target_var = Variable(target_normed.cuda(non_blocking=True))
        else:
            with torch.no_grad():
                target_var = Variable(target_normed)

        # compute output
        output = model(*input_var)
        loss = criterion(output, target_var)

        # measure accuracy and record loss
        if args.task == 'regression':
            mae_error = mae(normalizer.denorm(output.data.cpu()), target)
            losses.update(loss.data.cpu().item(), target.size(0))
            mae_errors.update(mae_error, target.size(0))
            if test:
                test_pred = normalizer.denorm(output.data.cpu())
                test_target = target
                test_preds += test_pred.view(-1).tolist()
                test_targets += test_target.view(-1).tolist()
                test_cif_ids += batch_cif_ids
        else:
            accuracy, precision, recall, fscore, auc_score = \
                class_eval(output.data.cpu(), target)
            losses.update(loss.data.cpu().item(), target.size(0))
            accuracies.update(accuracy, target.size(0))
            precisions.update(precision, target.size(0))
            recalls.update(recall, target.size(0))
            fscores.update(fscore, target.size(0))
            auc_scores.update(auc_score, target.size(0))
            if test:
                test_pred = torch.exp(output.data.cpu())
                test_target = target
                assert test_pred.shape[1] == 2
                test_preds += test_pred[:, 1].tolist()
                test_targets += test_target.view(-1).tolist()
                test_cif_ids += batch_cif_ids

        # measure elapsed time
        batch_time.update(time.time() - end)
        end = time.time()

        if i % args.print_freq == 0:
            if args.task == 'regression':
                print('Test: [{0}/{1}]\t'
                      'Time {batch_time.val:.3f} ({batch_time.avg:.3f})\t'
                      'Loss {loss.val:.4f} ({loss.avg:.4f})\t'
                      'MAE {mae_errors.val:.3f} ({mae_errors.avg:.3f})'.format(
                    i, len(val_loader), batch_time=batch_time, loss=losses,
                    mae_errors=mae_errors))
            else:
                print('Test: [{0}/{1}]\t'
                      'Time {batch_time.val:.3f} ({batch_time.avg:.3f})\t'
                      'Loss {loss.val:.4f} ({loss.avg:.4f})\t'
                      'Accu {accu.val:.3f} ({accu.avg:.3f})\t'
                      'Precision {prec.val:.3f} ({prec.avg:.3f})\t'
                      'Recall {recall.val:.3f} ({recall.avg:.3f})\t'
                      'F1 {f1.val:.3f} ({f1.avg:.3f})\t'
                      'AUC {auc.val:.3f} ({auc.avg:.3f})'.format(
                    i, len(val_loader), batch_time=batch_time, loss=losses,
                    accu=accuracies, prec=precisions, recall=recalls,
                    f1=fscores, auc=auc_scores))
    del input
    del batch_cif_ids
    del target
    torch.cuda.empty_cache()

    if test:
        star_label = '**'
        import csv
        with open('test_results.csv', 'w') as f:
            writer = csv.writer(f)
            for cif_id, target, pred in zip(test_cif_ids, test_targets,
                                            test_preds):
                writer.writerow((cif_id, target, pred))
    else:
        star_label = '*'
    if args.task == 'regression':
    	print(' {star} MAE {mae_errors.avg:.3f}'.format(star=star_label, mae_errors=mae_errors))
    	return mae_errors.avg
    else:
        print(' {star} AUC {auc.avg:.3f}'.format(star=star_label,
                                                 auc=auc_scores))
        return auc_scores.avg

class Normalizer(object):
    """Normalize a Tensor and restore it later. """

    def __init__(self, tensor):
        """tensor is taken as a sample to calculate the mean and std"""
        self.mean = torch.mean(tensor)
        self.std = torch.std(tensor)

    def norm(self, tensor):
        return (tensor - self.mean) / self.std

    def denorm(self, normed_tensor):
        return normed_tensor * self.std + self.mean

    def state_dict(self):
        return {'mean': self.mean,
                'std': self.std}

    def load_state_dict(self, state_dict):
        self.mean = state_dict['mean']
        self.std = state_dict['std']

def mae(prediction, target):
    """
    Computes the mean absolute error between prediction and target

    Parameters
    ----------

    prediction: torch.Tensor (N, 1)
    target: torch.Tensor (N, 1)
    """
    return torch.mean(torch.abs(target - prediction))

def class_eval(prediction, target):
    prediction = np.exp(prediction.numpy())
    target = target.numpy()
    pred_label = np.argmax(prediction, axis=1)
    target_label = np.squeeze(target)
    if not target_label.shape:
        target_label = np.asarray([target_label])
    if prediction.shape[1] == 2:
        precision, recall, fscore, _ = metrics.precision_recall_fscore_support(
            target_label, pred_label, average='binary')
        auc_score = metrics.roc_auc_score(target_label, prediction[:, 1])
        accuracy = metrics.accuracy_score(target_label, pred_label)
    else:
        raise NotImplementedError
    return accuracy, precision, recall, fscore, auc_score

class AverageMeter(object):
    """Computes and stores the average and current value"""

    def __init__(self):
        self.reset()

    def reset(self):
        self.val = 0
        self.avg = 0
        self.sum = 0
        self.count = 0

    def update(self, val, n=1):
        self.val = val
        self.sum += val * n
        self.count += n
        self.avg = self.sum / self.count

def save_checkpoint(state, is_best, filename='checkpoint.pth.tar'):
    torch.save(state, filename)
    if is_best:
        shutil.copyfile(filename, 'model_best.pth.tar')

def adjust_learning_rate(optimizer, epoch, k):
    """Sets the learning rate to the initial LR decayed by 10 every k epochs"""
    assert type(k) is int
    lr = args.lr * (0.1 ** (epoch // k))
    for param_group in optimizer.param_groups:
        param_group['lr'] = lr

if __name__ == '__main__':
    main()

max_stars_repo_name: thazhemadam/OGCNN

max_issues_repo_name: thazhemadam/OGCNN

max_forks_repo_name: thazhemadam/OGCNN


text:
import functools
import re
from datetime import datetime
from tempfile import TemporaryFile
from typing import Any, List, Tuple

import matplotlib.pyplot as plt
import numpy as np
from discord.ext.commands import Cog, Context, group
from tabulate import tabulate

from xythrion.bot import Xythrion
from xythrion.constants import WeatherAPIs
from xythrion.utils import DefaultEmbed, Graph, c2f, check_for_subcommands, http_get, k2c, k2f

EARTH_URL = "https://api.openweathermap.org/data/2.5/forecast?zip={0},{1}&appid={2}"
MARS_URL = f"https://api.nasa.gov/insight_weather/?api_key={WeatherAPIs.MARS}&feedtype=json&ver=1.0"

ZIP_CODE_PATTERN = re.compile(r"^\d{5}$")

class Weather(Cog):
    """Weather for different planets."""

    def __init__(self, bot: Xythrion) -> None:
        self.bot = bot

    @group()
    async def weather(self, ctx: Context) -> None:
        """Getting Weather for different planets."""
        if ctx.invoked_subcommand is None:
            await check_for_subcommands(ctx)

    @weather.command()
    async def earth(self, ctx: Context, zip_code: str, country_code: str = "US") -> None:
        """Getting weather for the planet of Earth."""
        if (zip_code := re.fullmatch(ZIP_CODE_PATTERN, zip_code)) is not None:
            zip_code = zip_code.group(0)

        else:
            await ctx.send(f"`Zip code must be 5 integers long. Received zip code: '{zip_code}'`")
            return

        json = await http_get(ctx, EARTH_URL.format(zip_code, country_code.upper(), WeatherAPIs.EARTH))

        lst, dates = [], []
        for i in json["list"]:
            lst.append(
                [
                    k2f(i["main"]["temp"]),
                    k2c(i["main"]["temp"]),
                    i["main"]["humidity"],
                    i["wind"]["speed"],
                ]
            )
            dates.append(datetime.fromtimestamp(i["dt"]).strftime("%a: %H%p").lower().title())

        titles = ["°F", "°C", "Humidity (%)", "Wind (m/s)"]

        func = functools.partial(self.create_weather_graph_and_table, ctx, lst, titles, dates, "Time")
        embed, table = await self.bot.loop.run_in_executor(None, func)

        embed.title = "**Weather on Earth.**"

        await ctx.send(file=embed.file, embed=embed, content=table)

    @weather.command()
    async def mars(self, ctx: Context) -> None:
        """Getting weather for the planet of Mars."""
        json = await http_get(ctx, MARS_URL)
        sols = json["sol_keys"]
        lst = []
        titles = ["°F", "°C", "Pressure (Pa)", "Wind (m/s)"]

        for sol in sols:
            try:
                i = json[sol]
                c = i["AT"]["av"]
                lst.append([c2f(c), c, i["PRE"]["av"], i["HWS"]["av"]])

            except KeyError:
                break

        func = functools.partial(self.create_weather_graph_and_table, ctx, lst, titles, sols, "Sol")
        embed, table = await self.bot.loop.run_in_executor(None, func)

        embed.title = f"**Weather on Mars sols {sols[0]}-{sols[-1]}.**"

        await ctx.send(file=embed.file, embed=embed, content=table)

    def create_weather_graph_and_table(
        self,
        ctx: Context,
        data_lst: List[List[float]],
        titles: List[str],
        days: List[str],
        day_title: str,
    ) -> Tuple[DefaultEmbed, str]:
        """Manipulating JSON data from weather APIs."""
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)
        axes = [ax1, ax2, ax3, ax4]

        lst = np.array(data_lst)

        lst = [lst[:, i] for i in range(lst.shape[1])]

        for i, (ax, title) in enumerate(zip(axes, titles)):
            ax.plot(lst[i])
            ax.set_title(title)
            ax.set_xticklabels(days, rotation=30)

        with TemporaryFile(suffix=".png") as buffer:
            with Graph(ctx, buffer, fig=fig, ax=axes) as embed:
                return embed, self.create_table(days, day_title, titles, data_lst)

    @staticmethod
    def create_table(days: List[str], day_title: str, titles: List[str], lst: List[Any]) -> str:
        """Creates a table from the tabulate module."""
        table = tabulate(
            [[days[i]] + x for i, x in enumerate(lst)],
            [day_title, *titles],
            tablefmt="simple",
            showindex=False,
            numalign="left",
            stralign="right",
            floatfmt=".2f",
        )

        return f"```py\n{table}```"


text:
# -*- coding: UTF-8 -*-
"""
此脚本用于随机生成线性模型数据、定义模型以及其他工具
"""

import numpy as np
import tensorflow as tf
import os
import math

def generate_linear_data(dimension, num):
    """
    随机产生线性模型数据
    参数
    ----
    dimension ：int，自变量个数
    num ：int，数据个数
    返回
    ----
    x ：np.array，自变量
    y ：np.array，因变量
    """
    np.random.seed(1024)
    beta = np.array(range(dimension)) + 1
    x = np.random.random((num, dimension))
    epsilon = np.random.random((num, 1))
    # 将被预测值写成矩阵形式，会极大加快速度
    y = x.dot(beta).reshape((-1, 1)) + epsilon
    return x, y

def create_linear_model(dimension):
    """
    搭建模型，包括数据中的自变量，应变量和损失函数
    参数
    ----
    dimension : int，自变量的个数
    返回
    ----
    model ：dict，里面包含模型的参数，损失函数，自变量，应变量
    """
    np.random.seed(1024)
    # 定义自变量和应变量
    x = tf.placeholder(tf.float64, shape=[None, dimension], name='x')
    # 将被预测值写成矩阵形式，会极大加快速度
    y = tf.placeholder(tf.float64, shape=[None, 1], name="y")
    # 定义参数估计值和预测值
    beta_pred = tf.Variable(np.random.random([dimension, 1]))
    y_pred = tf.matmul(x, beta_pred, name="y_pred")
    # 定义损失函数
    loss = tf.reduce_mean(tf.square(y_pred - y))
    model = {"loss_function": loss, "independent_variable": x,
             "dependent_variable": y, "prediction": y_pred, "model_params": beta_pred}
    return model

def create_summary_writer(log_path):
    """
    检查所给路径是否已存在，如果存在删除原有日志。并创建日志写入对象
    参数
    ----
    logPath ：string，日志存储路径
    返回
    ----
    summaryWriter ：FileWriter，日志写入器
    """
    if tf.gfile.Exists(log_path):
        tf.gfile.DeleteRecursively(log_path)
    summary_writer = tf.summary.FileWriter(log_path, graph=tf.get_default_graph())
    return summary_writer

def gradient_descent(X, Y, model, learning_rate=0.01, max_iter=10000, tol=1.e-6):
    """
    利用梯度下降法训练模型。
    参数
    ----
    X : np.array, 自变量数据
    Y : np.array, 因变量数据
    model : dict, 里面包含模型的参数，损失函数，自变量，应变量。
    """
    # 确定最优化算法
    method = tf.train.GradientDescentOptimizer(learning_rate=learning_rate)
    optimizer = method.minimize(model["loss_function"])
    # 增加日志
    tf.summary.scalar("loss_function", model["loss_function"])
    tf.summary.histogram("params", model["model_params"])
    tf.summary.scalar("first_param", tf.reduce_mean(model["model_params"][0]))
    tf.summary.scalar("last_param", tf.reduce_mean(model["model_params"][-1]))
    summary = tf.summary.merge_all()
    # 在程序运行结束之后，运行如下命令，查看日志
    # tensorboard --logdir logs/
    # Windows下的存储路径与Linux并不相同
    if os.name == "nt":
        summary_writer = create_summary_writer("logs\\gradient_descent")
    else:
        summary_writer = create_summary_writer("logs/gradient_descent")
    # tensorflow开始运行
    sess = tf.Session()
    # 产生初始参数
    init = tf.global_variables_initializer()
    # 用之前产生的初始参数初始化模型
    sess.run(init)
    # 迭代梯度下降法
    step = 0
    prev_loss = np.inf
    diff = np.inf
    # 当损失函数的变动小于阈值或达到最大循环次数，则停止迭代
    while (step < max_iter) & (diff > tol):
        _, summary_str, loss = sess.run(
            [optimizer, summary, model["loss_function"]],
            feed_dict={model["independent_variable"]: X,
                       model["dependent_variable"]: Y})
        # 将运行细节写入目录
        summary_writer.add_summary(summary_str, step)
        # 计算损失函数的变动
        diff = abs(prev_loss - loss)
        prev_loss = loss
        step += 1
    summary_writer.close()

def stochastic_gradient_descent(X, Y, model, method,
                                mini_batch_fraction=0.01, epoch=10000, tol=1.e-6):
    """
    利用随机梯度下降法训练模型。
    参数
    ----
    X : np.array, 自变量数据
    Y : np.array, 因变量数据
    model : dict, 里面包含模型的参数，损失函数，自变量，应变量
    """
    # 确定最优化算法
    optimizer = method.minimize(model["loss_function"])
    # 增加日志
    tf.summary.scalar("loss_function", model["loss_function"])
    tf.summary.histogram("params", model["model_params"])
    tf.summary.scalar("first_param", tf.reduce_mean(model["model_params"][0]))
    tf.summary.scalar("last_param", tf.reduce_mean(model["model_params"][-1]))
    summary = tf.summary.merge_all()
    # 在程序运行结束之后，运行如下命令，查看日志
    # tensorboard --logdir logs/
    # Windows下的存储路径与Linux并不相同
    if os.name == "nt":
        summary_writer = create_summary_writer("logs\\stochastic_gradient_descent")
    else:
        summary_writer = create_summary_writer("logs/stochastic_gradient_descent")
    # tensorflow开始运行
    sess = tf.Session()
    # 产生初始参数
    init = tf.global_variables_initializer()
    # 用之前产生的初始参数初始化模型
    sess.run(init)
    # 迭代梯度下降法
    step = 0
    batch_size = int(X.shape[0] * mini_batch_fraction)
    batch_num = int(math.ceil(1 / mini_batch_fraction))
    prev_loss = np.inf
    diff = np.inf
    # 当损失函数的变动小于阈值或达到最大训练轮次，则停止迭代
    while (step < epoch) & (diff > tol):
        for i in range(batch_num):
            # 选取小批次训练数据
            batch_x = X[i * batch_size: (i + 1) * batch_size]
            batch_y = Y[i * batch_size: (i + 1) * batch_size]
            # 迭代模型参数
            sess.run([optimizer],
                     feed_dict={model["independent_variable"]: batch_x,
                                model["dependent_variable"]: batch_y})
            # 计算损失函数并写入日志
            summary_str, loss = sess.run(
                [summary, model["loss_function"]],
                feed_dict={model["independent_variable"]: X,
                           model["dependent_variable"]: Y})
            # 将运行细节写入目录
            summary_writer.add_summary(summary_str, step * batch_num + i)
            # 计算损失函数的变动
            diff = abs(prev_loss - loss)
            prev_loss = loss
            if diff <= tol:
                break
        step += 1
    summary_writer.close()


text:
import cv2
import numpy as np
import random
import os

threshold = 0.85

#changed how many pixels are drawn at once
pixelsToShow = 100
#causes pixel fill to bloom outwards rather than scanlines. This is slower
randomizeStack = False

filename = "pic (5)"

img_rgb = cv2.imread("clean images/" + filename + ".png")
img_gray = cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY)

areas = []

def fill(data, original, sx, sy, fill_value):
    global pixelsToShow, randomizeStack
    ysize, xsize, col = data.shape

    stack = []
    stack.append((sx, sy))

    timer = pixelsToShow
    while stack:
        x = 0
        y = 0
        if randomizeStack:
            i = random.randint(0,len(stack)-1)
            x, y = stack[i]
            del stack[i]
        else:
            x, y = stack.pop()

        px = original[y, x]
        if px[0] > 50 and px[1] > 50 and px[2] > 50:

            data[y, x] = (fill_value[0] * (px[0]/255),
                          fill_value[1] * (px[0]/255),
                          fill_value[2] * (px[0]/255))

            original[y,x] = (0,0,0)

            timer = timer - 1
            if (timer == 0):
                cv2.imshow("rgb", img_rgb)
                cv2.waitKey(1);
                timer = pixelsToShow

            if x > 0:
                stack.append((x - 1, y))
            if x < (xsize - 1):
                stack.append((x + 1, y))
            if y > 0:
                stack.append((x, y - 1))
            if y < (ysize - 1):
                stack.append((x, y + 1))

    return data

path, dirs, files = next(os.walk('clean images/numbers/' + filename + '/'))
numbersInFile = len(files) + 1
for i in range(1,numbersInFile):
    number = template = cv2.imread('clean images/numbers/' + filename + '/' + '(' + str(i) + ').png',0)
    w, h = template.shape[::-1]

    num = cv2.matchTemplate(img_gray,template,cv2.TM_CCOEFF_NORMED)
    loc = np.where( num >= threshold)
    for pt in zip(*loc[::-1]):
        cv2.rectangle(img_rgb, pt, (pt[0] + w, pt[1] + h), (0,0,255), 1)
        areas.append((i, pt[0], pt[1]))
        #uncomment the 4 lines below to see it select one number from the image at a time
#    cv2.rectangle(img_rgb,(0,0),(50,40),0,-1)
#    img_rgb = cv2.putText(img_rgb, str(i), (20,25), cv2.FONT_HERSHEY_COMPLEX,
#                   1, (0,0,255), 1, cv2.LINE_AA)
#    cv2.imshow("rgb", img_rgb)
#    cv2.waitKey();
    for pt in zip(*loc[::-1]):
        cv2.rectangle(img_rgb, pt, (pt[0] + w, pt[1] + h), (255,255,255), -1)
    img_gray = cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY)

random.shuffle(areas)
#print(areas)

colorfile = open('clean images/colors/' + filename +'.txt', 'r')
lines = colorfile.readlines()
colorArray = []

# Strips the newline character
for i in range(1,len(lines)):
    line = lines[i].strip().split(",")
    line = list(map(int, line))
    colorArray.append(line)

print(colorArray)

outlines = np.copy(img_rgb)
for i in range(len(areas)):
    print(i)
    color = colorArray[areas[i][0]-1]

    img_rgb = fill(img_rgb, outlines, areas[i][1]+5, areas[i][2]+5, color)

print("DONE")
cv2.imshow("rgb", img_rgb)
cv2.waitKey();
cv2.destroyAllWindows();

cv2.imwrite('output/' + filename +'.png', img_rgb)

max_stars_repo_path: Color Images.py

max_issues_repo_path: Color Images.py

max_forks_repo_path: Color Images.py