text:
'''Generate images, plots and tables for bachelor's thesis paper.'''

import csv
import os

import matplotlib.pyplot as plt
import numpy as np
from skimage.color import rgb2gray
from skimage.filters import threshold_otsu
from skimage.io import imread, imsave
from skimage.transform import rescale
from skimage.util import invert, img_as_ubyte
import tikzplotlib
from sklearn.model_selection import train_test_split
from tensorflow import keras

from blob_analysis import (count_blobs_with_all_methods, patch_plot_legend,
                           plot_blob_stat)
from blob_detection_compare_demo import compare_detection
from blob_series_tracker import (find_blob_series,
                                 ratio_of_remaining_blobs_in_stages)
from img_processing import (crop_ui, default_img_set, full_prepare,
                            load_img_series)
from neural_network import (default_grain_classifier_model,
                            network_cross_validation, mean_confusion_matrix)

def temp_bounds_imgs_gen():
    img = imread('img/103_E5R_1.jpg')
    img = rgb2gray(img)
    img = invert(img)

    bounds = (((6, 24), (283, 318)),
              ((219, 236), (283, 318)))

    for bound in bounds:
        bound_img = img[slice(*bound[0]), slice(*bound[1])]
        bound_img = rescale(bound_img, 4, anti_aliasing=True)
        thr = threshold_otsu(bound_img)
        img_txt = bound_img > thr

        plt.figure()
        plt.imshow(bound_img, cmap="gray")
        plt.axis('off')
        imsave('exports/temp_bounds_scale.png', img_as_ubyte(bound_img))
        plt.figure()
        plt.axis('off')
        plt.imshow(img_txt, cmap="gray")
        imsave('exports/temp_bounds_bin.png', img_as_ubyte(img_txt))

def grain_samples_imgs_gen():
    samples_names = ('104_E5R', '117_E6R')
    for name in samples_names:
        imgs = load_img_series('img/' + name)
        for i, img in enumerate(imgs):
            img = rgb2gray(img)
            imsave('exports/' + name + '_' + str(i) + '.png',
                   img_as_ubyte(crop_ui(img)))

def blob_detection_compare_plots_gen():
    img = imread('img/104_E5R_0.jpg')
    img_crop = crop_ui(rgb2gray(img))
    img_prep = full_prepare(img)
    blobs_list = compare_detection(img_prep)

    suffixes = ('LoG', 'DoG', 'DoH')

    for blobs, suffix in zip(blobs_list, suffixes):
        _, ax = plt.subplots()
        plt.title('Liczba wykrytych detali: {}'.format(len(blobs)))
        plt.imshow(img_crop, cmap=plt.get_cmap('gray'))
        for blob in blobs:
            y, x, r = blob
            c = plt.Circle((x, y), r, color='r', fill=False)
            ax.add_patch(c)

        ax.set_axis_off()
        tikzplotlib.save('exports/blob_detection_compare_' + suffix)

def blob_count_plots_gen():
    imgs = load_img_series('img/104_E5R')
    imgs_prep = [full_prepare(img) for img in imgs]
    imgs_crop = [crop_ui(rgb2gray(img)) for img in imgs]

    stages_all = find_blob_series(imgs_prep, only_remaining=False)
    stages_rem = find_blob_series(imgs_prep)

    # Map stages on first image
    colors = ('blue', 'blueviolet', 'magenta', 'crimson', 'red')
    fig = plt.figure(frameon=False)
    ax = fig.add_axes([0, 0, 1, 1])
    plt.imshow(imgs_crop[0], cmap=plt.get_cmap('gray'))
    for stage, color in zip(stages_rem, colors):
        for blob in stage:
            y, x, r = blob
            c = plt.Circle((x, y), r, color=color, fill=False)
            ax.add_patch(c)
    ax.set_axis_off()
    plt.savefig('exports/blob_tracker', dpi=300)

    # Show two methods combined to compare
    loop_set = enumerate(zip(stages_rem, stages_all, imgs_crop))
    for i, (stage_rem, stage_all, img) in loop_set:
        fig = plt.figure(frameon=False)
        ax = fig.add_axes([0, 0, 1, 1])
        plt.imshow(img, cmap=plt.get_cmap('gray'))
        for blob_all in stage_all:
            y, x, r = blob_all
            c = plt.Circle((x, y), r, color='b', fill=False)
            ax.add_patch(c)
        for blob_rem in stage_rem:
            y, x, r = blob_rem
            c = plt.Circle((x, y), r, color='r', fill=False)
            ax.add_patch(c)

        ax.set_axis_off()
        plt.savefig('exports/blob_tracker_min_' + str(i))

def blob_ratio_table_gen():
    sample_names = ('104_E5R', '106_E11R', '107_E6R', '111_E16R')

    with open('exports/neural_network_comparison.csv', 'w') as csvfile:
        filewriter = csv.writer(csvfile, delimiter=';')
        # Header
        filewriter.writerow(('Próbka', 'Minuta 0', 'Minuta 1', 'Minuta 2',
                             'Minuta 3', 'Minuta 4'))
        for name in sample_names:
            imgs = load_img_series('img/' + name)
            imgs_prep = [full_prepare(img) for img in imgs]
            stages_rem = find_blob_series(imgs_prep)
            ratios = ratio_of_remaining_blobs_in_stages(stages_rem)
            ratios = (round(ratio, 2) for ratio in ratios)
            filewriter.writerow((name, *ratios))

def blob_analysis_plots_gen():
    X, y = default_img_set()
    X = [[full_prepare(img) for img in same_sample] for same_sample in X]

    Xa, Xr, Xp = count_blobs_with_all_methods(X)

    colors = ('r', 'g', 'b', 'y')
    labels = ('E5R', 'E11R', 'E6R', 'E16R')

    plot_blob_stat(Xa, y, colors)
    plt.title('Liczba wszystkich detali')
    plt.xlabel('Minuty')
    plt.ylabel('Liczba detali')
    patch_plot_legend(colors, labels)
    tikzplotlib.save('exports/blob_analysis_all')

    plot_blob_stat(Xr, y, colors)
    plt.title('Liczba śledzonych detali')
    plt.xlabel('Minuty')
    plt.ylabel('Liczba detali')
    patch_plot_legend(colors, labels)
    tikzplotlib.save('exports/blob_analysis_remaining')

    plot_blob_stat(Xp, y, colors)
    plt.title('Pozostały procent śledzonych detali')
    plt.xlabel('Minuty')
    plt.ylabel('Liczba detali')
    patch_plot_legend(colors, labels)
    tikzplotlib.save('exports/blob_analysis_ratio')

def neural_network_trainig_plots_gen():
    X, y = default_img_set()
    X = [[full_prepare(img) for img in same_sample] for same_sample in X]
    Xs = count_blobs_with_all_methods(X)

    files_suffixes = ('all', 'remaining', 'ratio')

    for X, suffix in zip(Xs, files_suffixes):
        X = np.array(X)
        y = np.array(y)

        X_train, _, y_train, _ = train_test_split(
            X, y, stratify=y, test_size=0.33, random_state=1)

        model = default_grain_classifier_model()
        model.compile(
            optimizer='adam',
            loss='sparse_categorical_crossentropy',
            metrics=['accuracy'])
        history = model.fit(X_train, y_train, epochs=300, verbose=0)

        plt.figure()
        plt.title('Historia treningu modelu')
        plt.xlabel('Epoka')

        plt.plot(history.history['accuracy'], c='b')
        plt.plot(history.history['loss'], c='r')

        plt.legend(('Dokładność', 'Błąd'))

        tikzplotlib.save('exports/neural_network_trainig_' + suffix)

def neural_network_test_table_gen():
    X, y = default_img_set()
    X = [[full_prepare(img) for img in same_sample] for same_sample in X]
    Xs = count_blobs_with_all_methods(X)
    Xs = [np.array(X_count) for X_count in Xs]
    y = np.array(y)

    row_names = ('Wszystkie detale', 'Śledzone detale',
                 'Stosunek śledzonych detali')

    with open('exports/neural_network_test.csv', 'w') as csvfile:
        filewriter = csv.writer(csvfile, delimiter=';')
        # Header
        filewriter.writerow(('Metoda zliczania detali', 'Wskaźnik',
                             'wskaźnik'))
        filewriter.writerow(('Metoda zliczania detali', 'Błąd', 'Dokładność'))

        for X, name in zip(Xs, row_names):
            X = np.array(X)
            X_train, X_test, y_train, y_test = train_test_split(
                X, y, stratify=y, test_size=0.33, random_state=1)

            model = default_grain_classifier_model()
            model.compile(
                optimizer='adam',
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])
            model.fit(X_train, y_train, epochs=300, verbose=0)

            score = model.evaluate(X_test, y_test, verbose=0)
            filewriter.writerow((name, *score))

def neural_network_validation_table_gen():
    X, y = default_img_set()
    X = [[full_prepare(img) for img in same_sample] for same_sample in X]
    Xs = count_blobs_with_all_methods(X)
    Xs = [np.array(X_count) for X_count in Xs]
    y = np.array(y)

    row_names = ('Wszystkie detale', 'Śledzone detale',
                 'Stosunek śledzonych detali')

    with open('exports/neural_network_validation.csv', 'w') as csvfile:
        filewriter = csv.writer(csvfile, delimiter=';')
        # Header
        filewriter.writerow(('Metoda zliczania detali', 'Wskaźnik',
                     'wskaźnik'))
        filewriter.writerow(('Metoda zliczania detali', 'Błąd', 'Dokładność'))

        for X, name in zip(Xs, row_names):
            model = default_grain_classifier_model()
            model.compile(
                optimizer='adam',
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])

            scores = np.array(network_cross_validation(model, X, y, 3))
            score = np.round(scores.mean(axis=0), 2)
            filewriter.writerow((name, *score))

def network_comparison_table_gen():
    X, y = default_img_set()
    X = [[full_prepare(img) for img in same_sample] for same_sample in X]
    X = [
        ratio_of_remaining_blobs_in_stages(find_blob_series(img_series))
        for img_series in X
    ]
    X = np.array(X)
    y = np.array(y)

    with open('exports/neural_network_comparison.csv', 'w') as csvfile:
        filewriter = csv.writer(csvfile, delimiter=';')
        # Header
        filewriter.writerow(('Parametr', 'Wartość', 'Błąd', 'Dokładność'))

        # Activation functions
        activation_funcs = ('sigmoid', 'relu', 'elu', 'tanh')
        for func in activation_funcs:
            model = keras.Sequential([
                keras.layers.Dense(5, activation=func),
                keras.layers.Dense(256, activation=func),
                keras.layers.Dense(128, activation=func),
                keras.layers.Dense(4, activation='softmax')
            ])

            model.compile(
                optimizer='adam',
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])
            scores = np.array(network_cross_validation(model, X, y, 3))
            score = np.round(scores.mean(axis=0), 2)

            filewriter.writerow(('Funkcja aktywacji', func, *score))

        # Number of hidden layers
        models = []
        models.append(
            keras.Sequential([
                keras.layers.Dense(5, activation='tanh'),
                keras.layers.Dense(512, activation='tanh'),
                keras.layers.Dense(4, activation='softmax')
            ]))
        models.append(
            keras.Sequential([
                keras.layers.Dense(5, activation='tanh'),
                keras.layers.Dense(256, activation='tanh'),
                keras.layers.Dense(128, activation='tanh'),
                keras.layers.Dense(4, activation='softmax')
            ]))

        for model, i in zip(models, (1, 2)):
            model.compile(
                optimizer='adam',
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])
            scores = np.array(network_cross_validation(model, X, y, 3))
            score = np.round(scores.mean(axis=0), 2)

            filewriter.writerow(('Liczba warstw ukrytych', i, *score))

        # Number of neurons in hidden layers
        neurons_num = ((128, 64), (256, 128), (512, 126))

        for num in neurons_num:
            model = keras.Sequential([
                keras.layers.Dense(5, activation='tanh'),
                keras.layers.Dense(num[0], activation='tanh'),
                keras.layers.Dense(num[1], activation='tanh'),
                keras.layers.Dense(4, activation='softmax')
            ])

            model.compile(
                optimizer='adam',
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])
            scores = np.array(network_cross_validation(model, X, y, 3))
            score = np.round(scores.mean(axis=0), 2)

            filewriter.writerow(('Liczba neuronów w warstwach ukrytych',
                                 '{} i {}'.format(num[0], num[1]), *score))

        # Optimizer
        model = default_grain_classifier_model()
        optimizers = ('sgd', 'adam')
        for opt in optimizers:
            model.compile(
                optimizer=opt,
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])
            scores = np.array(network_cross_validation(model, X, y, 3))
            score = np.round(scores.mean(axis=0), 2)

            filewriter.writerow(('Algorytm uczenia', opt, *score))

def confusion_matrix_table_gen():
    X, y = default_img_set()
    X = [[full_prepare(img) for img in same_sample] for same_sample in X]
    X = count_blobs_with_all_methods(X)[2]
    X = np.array(X)
    y = np.array(y)

    model = default_grain_classifier_model()
    model.compile(
        optimizer='adam',
        loss='sparse_categorical_crossentropy',
        metrics=['accuracy'])

    mcm = mean_confusion_matrix(model, X, y, 3)
    np.savetxt(
        "exports/mean_confusion_matrix_ratio.csv",
        mcm,
        fmt='%.2f',
        delimiter=";")

def clear_dir(dir_to_clear):
    for root, _, files in os.walk(dir_to_clear):
        for file_to_rm in files:
            os.remove(os.path.join(root, file_to_rm))

if __name__ == '__main__':
    if not os.path.exists('exports'):
        os.makedirs('exports')
    else:
        clear_dir('exports')

    temp_bounds_imgs_gen()
    grain_samples_imgs_gen()
    blob_detection_compare_plots_gen()
    blob_ratio_table_gen()
    blob_count_plots_gen()
    blob_analysis_plots_gen()
    neural_network_test_table_gen()
    neural_network_validation_table_gen()
    neural_network_trainig_plots_gen()
    network_comparison_table_gen()
    confusion_matrix_table_gen()

    plt.show()

max_stars_repo_name: MaciejZj/Copper-grains-classification

max_issues_repo_name: MaciejZj/Copper-grains-classification

max_forks_repo_name: MaciejZj/Copper-grains-classification


text:
from pathlib import Path
import os

import numpy as np
import pygeos as pg
import rasterio
from rasterio.features import rasterize
from pyogrio import read_dataframe

from analysis.constants import DATA_CRS, MASK_FACTOR
from analysis.lib.io import write_raster
from analysis.lib.input_areas import get_input_area_mask, get_input_area_boundary
from analysis.lib.raster import add_overviews, create_lowres_mask, extract_window
from analysis.lib.speedups import remap
from analysis.lib.pygeos_util import to_dict_all, dissolve, explode

src_dir = Path("source_data/natures_network")
data_dir = Path("data/inputs")
out_dir = data_dir / "indicators/natures_network"
outfilename = out_dir / "nn_priority.tif"

if not out_dir.exists():
    os.makedirs(out_dir)

### Get input area mask
print("Extracting Nature's Network input area mask...")
mask, transform, window = get_input_area_mask("nn")
nodata = 255

print("Reading and warping Nature's Network...")
with rasterio.open(src_dir / "NaturesNetwork_conservdesign_180625.tif") as src:
    # note: raster does not have nodata set; 0 indicates NODATA (outside extent)
    # and 0 values
    data = extract_window(src, window, transform, nodata=nodata)

# apply input area mask
data = np.where(mask == 1, data, nodata).astype("uint8")

print("Reclassifying data...")
# Remap the raw values to priorities and categories
table = read_dataframe(
    src_dir / "NaturesNetwork_conservdesign_180625.tif.vat.dbf",
    columns=["Value", "Priority", "Descrpt"],
    read_geometry=False,
)
table = table.loc[table.Value > 0].copy()
table.Priority = table.Priority.astype("uint8")
table["category"] = table.Descrpt.str[0].astype("uint8")

remap_table = table[["Value", "Priority"]].values.astype("uint8")
priority_data = remap(data, remap_table, nodata=nodata)

write_raster(outfilename, priority_data, transform, DATA_CRS, nodata=nodata)

print("Adding overviews and masks...")
add_overviews(outfilename)

create_lowres_mask(
    outfilename,
    str(outfilename).replace(".tif", "_mask.tif"),
    factor=MASK_FACTOR,
    ignore_zero=False,
)

### Process indicators

# Extract cores from the categories
table["aquatic_cores"] = table.Descrpt.str.contains("Aquatic")
table["terrestrial_cores"] = table.Descrpt.str.contains("Terrestrial")
table["imperiled_cores"] = table.Descrpt.str.contains("Imperiled")

for core_type in ["aquatic_cores", "terrestrial_cores", "imperiled_cores"]:
    print(f"Processing {core_type}")
    cores_data = remap(
        data, table[["Value", core_type]].values.astype("uint8"), nodata=nodata
    )

    outfilename = out_dir / f"{core_type}.tif"
    write_raster(outfilename, cores_data, transform, DATA_CRS, nodata=nodata)
    add_overviews(outfilename)
    create_lowres_mask(
        outfilename,
        str(outfilename).replace(".tif", "_mask.tif"),
        factor=MASK_FACTOR,
        ignore_zero=True,
    )

### Process habitat condition
print("Extracting habitat condition...")
with rasterio.open(src_dir / "indicators/HABITAT_CONDITION_web.tif") as src:
    data = extract_window(src, window, transform, nodata=nodata)

# apply input area mask
data = np.where(mask == 1, data, nodata).astype("uint8")

# Classes used for symbology are 0-47 (lowest), 48-88, 89-200 (highest)
# see the *.sd file for this dataset
data[data <= 47] = 0
data[(data > 47) & (data <= 88)] = 1
data[(data > 88) & (data <= 200)] = 2

outfilename = out_dir / "habitat_condition.tif"
write_raster(outfilename, data, transform, DATA_CRS, nodata=nodata)
add_overviews(outfilename)
create_lowres_mask(
    outfilename,
    str(outfilename).replace(".tif", "_mask.tif"),
    factor=MASK_FACTOR,
    ignore_zero=False,
)

### Process Habitat Importance
print("Extracting habitat importance...")
with rasterio.open(src_dir / "indicators/habitat_importance.tif") as src:
    dtype = src.dtypes[0]
    data = extract_window(src, window, transform, nodata=src.nodata).astype(dtype)

table = read_dataframe(
    src_dir / "indicators/habitat_importance.tif.vat.dbf", read_geometry=False
)

# Classes on Impt_sum used for symbology are 0-22 (highest), 23-77, 78-200 (lowest)
# see the *.sd file for this dataset
table["new_value"] = 0
table.loc[table.Impt_sum <= 77, "new_value"] = 1
table.loc[table.Impt_sum <= 22, "new_value"] = 2

data = remap(data, table[["Value", "new_value"]].values.astype(dtype), nodata=nodata)

# apply input area mask
data = np.where(mask == 1, data, nodata).astype("uint8")

outfilename = out_dir / "habitat_importance.tif"
write_raster(outfilename, data, transform, DATA_CRS, nodata=nodata)
add_overviews(outfilename)
create_lowres_mask(
    outfilename,
    str(outfilename).replace(".tif", "_mask.tif"),
    factor=MASK_FACTOR,
    ignore_zero=False,
)

### Process Aquatic Ecological Integrity
print("Processing aquatic integrity")
with rasterio.open(
    src_dir / "indicators/AquaticIEI_v3_2/DSL_aquaiei-r_2010_v3.2.tif"
) as src:
    dtype = src.dtypes[0]
    data = extract_window(src, window, transform, nodata=src.nodata)

    # apply input area mask
    data = np.where(mask == 1, data, src.nodata).astype(dtype)

# this is on a continuous scale from 0 to 1; bin into 10 classes
data[data == src.nodata] = 255
data[(data >= 0) & (data <= 0.1)] = 0
data[(data > 0.1) & (data <= 0.2)] = 1
data[(data > 0.2) & (data <= 0.3)] = 2
data[(data > 0.3) & (data <= 0.4)] = 3
data[(data > 0.4) & (data <= 0.5)] = 4
data[(data > 0.5) & (data <= 0.6)] = 5
data[(data > 0.6) & (data <= 0.7)] = 6
data[(data > 0.7) & (data <= 0.8)] = 6
data[(data > 0.8) & (data <= 0.9)] = 8
data[(data > 0.9) & (data <= 1)] = 9

data = data.astype("uint8")

outfilename = out_dir / "aquatic_integrity.tif"
write_raster(outfilename, data, transform, DATA_CRS, nodata=nodata)
add_overviews(outfilename)
create_lowres_mask(
    outfilename,
    str(outfilename).replace(".tif", "_mask.tif"),
    factor=MASK_FACTOR,
    ignore_zero=True,
)

### Process terrestrial integrity
print("Processing terrestrial integrity")
with rasterio.open(src_dir / "indicators/IEI_2010_r_v3_2/iei-r_2010_v3.2.tif") as src:
    dtype = src.dtypes[0]
    data = extract_window(src, window, transform, nodata=src.nodata)

    # apply input area mask
    data = np.where(mask == 1, data, src.nodata).astype(dtype)

# this is on a continuous scale from 0 to 1; bin into 10 classes
data[data == src.nodata] = 255
data[(data >= 0) & (data <= 0.1)] = 0
data[(data > 0.1) & (data <= 0.2)] = 1
data[(data > 0.2) & (data <= 0.3)] = 2
data[(data > 0.3) & (data <= 0.4)] = 3
data[(data > 0.4) & (data <= 0.5)] = 4
data[(data > 0.5) & (data <= 0.6)] = 5
data[(data > 0.6) & (data <= 0.7)] = 6
data[(data > 0.7) & (data <= 0.8)] = 6
data[(data > 0.8) & (data <= 0.9)] = 8
data[(data > 0.9) & (data <= 1)] = 9

data = data.astype("uint8")

outfilename = out_dir / "integrity.tif"
write_raster(outfilename, data, transform, DATA_CRS, nodata=nodata)
add_overviews(outfilename)
create_lowres_mask(
    outfilename,
    str(outfilename).replace(".tif", "_mask.tif"),
    factor=MASK_FACTOR,
    ignore_zero=False,
)

### Process freshwater resilience by watershed

print("Processing freshwater resilience")

bnd = get_input_area_boundary("nn")
df = read_dataframe(
    src_dir
    / "indicators/Freshwater_Resilience/FW_resilience_highesthigh_watersheds.shp",
    columns=["RES_CLASS"],
).to_crs(DATA_CRS)

tree = pg.STRtree(df.geometry.values.data)
ix = tree.query(bnd, predicate="intersects")
df = df.iloc[ix].copy()

# remap values so that they are 1-4; 0 is fill value, 255 nodata
df["value"] = df.RES_CLASS.map(
    {
        "Complex: Highest Relative Resilience": 4,
        "Complex: High Relative Resilience": 3,
        "Non-Complex: Highest Relative Score": 2,
        "Non-Complex: High Relative Score": 2,
    }
)

# dissolve by value
df = dissolve(explode(df), by=["value"])

# rasterize
df["g"] = to_dict_all(df.geometry.values.data)
shapes = df[["g", "value"]].values.tolist()

data = rasterize(shapes, mask.shape, transform=transform, dtype="uint8", fill=0)

# apply input area mask
nodata = 255
data = np.where(mask == 1, data, nodata).astype("uint8")

outfilename = out_dir / "freshwater_resilience.tif"
write_raster(outfilename, data, transform, DATA_CRS, nodata=nodata)
add_overviews(outfilename)
create_lowres_mask(
    outfilename,
    str(outfilename).replace(".tif", "_mask.tif"),
    factor=MASK_FACTOR,
    ignore_zero=False,
)

ext: py

max_stars_repo_name: astutespruce/secas-blueprint

max_issues_repo_name: astutespruce/secas-blueprint

max_forks_repo_name: astutespruce/secas-blueprint

reason: import numpy


text:
import pickle
import numpy as np
import pandas as pd
import os
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score,confusion_matrix

df1 = pickle.load(open('../models/en-all-users-mbfc-fasttext-embed_df_1.pkl','rb'))
df2 = pickle.load(open('../models/en-all-users-mbfc-fasttext-embed_df_2.pkl','rb'))
df1['feature_vector']=df1.values.tolist()
df2['feature_vector']=df2.values.tolist()
df1 = df1.reset_index()
df2 = df2.reset_index()

df1 = df1.rename(columns={'index': 'user'})
df2 = df2.rename(columns={'index': 'user'})

df1=df1[['user','feature_vector']]
df2=df2[['user','feature_vector']]

#Overall
new_df=df1.append(df2,ignore_index=True)
mbfc_scores=pd.read_csv('../data/domain-score.csv')
subdf=new_df[new_df.user.isin(mbfc_scores['user'].tolist())]

combined_df=subdf.merge(mbfc_scores,on='user')
#Consider the top and bottom 30% of the data
s1 = combined_df['science'].quantile(0.3)
s2 = combined_df['science'].quantile(0.7)
p1 = combined_df['political'].quantile(0.3)
p2 = combined_df['political'].quantile(0.7)
m1 = combined_df['moderacy'].quantile(0.3)
m2 = combined_df['moderacy'].quantile(0.7)

hard_sci=[]
hard_politics=[]
hard_moderacy=[]

for i in range(len(combined_df)):
    if combined_df['science'].iloc[i]<=s1:
        hard_sci.append(0)
    elif combined_df['science'].iloc[i]>=s2:
        hard_sci.append(1.0)
    else:
        hard_sci.append(-100)
    if combined_df['political'].iloc[i]<=p1:
        hard_politics.append(0)
    elif combined_df['political'].iloc[i]>=p2:
        hard_politics.append(1.0)
    else:
        hard_politics.append(-100)
    if combined_df['moderacy'].iloc[i]<=m1:
        hard_moderacy.append(0)
    elif combined_df['moderacy'].iloc[i]>=m2:
        hard_moderacy.append(1.0)
    else:
        hard_moderacy.append(-100)
combined_df['hard_sci']=hard_sci
combined_df['hard_politics']=hard_politics
combined_df['hard_moderacy']=hard_moderacy

tmp = combined_df[combined_df['hard_sci']!=-100]
X_new=tmp['feature_vector']
Y_new=tmp['hard_sci']

X_train, X_test, y_train, y_test = train_test_split(X_new, Y_new,test_size=0.33,shuffle=True,random_state=23)

clf = LogisticRegression(random_state=23,max_iter=700).fit(np.asarray(X_train.tolist()), y_train.tolist())
y_pred=clf.predict(np.asarray(X_test.tolist()))
print("Science-Predictions")
print("Accuracy",accuracy_score(np.asarray(y_test.tolist()),y_pred))
print("Precision",precision_score(np.asarray(y_test.tolist()),y_pred))
print("Recall",recall_score(np.asarray(y_test.tolist()),y_pred))
print("F1-Score",f1_score(np.asarray(y_test.tolist()),y_pred))
tn, fp, fn, tp = confusion_matrix(np.asarray(y_test.tolist()), y_pred).ravel()

all_X = new_df['feature_vector'].tolist()
all_Y_pred = clf.predict(all_X)
new_df['science']=list(all_Y_pred)

tmp = combined_df[combined_df['hard_politics']!=-100]
X_new=tmp['feature_vector']
Y_new=tmp['hard_politics']

X_train, X_test, y_train, y_test = train_test_split(X_new, Y_new,test_size=0.33,shuffle=True,random_state=23)
print()
print("Political-Predictions")
clf = LogisticRegression(random_state=23,max_iter=700).fit(np.asarray(X_train.tolist()), y_train.tolist())
y_pred=clf.predict(np.asarray(X_test.tolist()))
print("Accuracy",accuracy_score(np.asarray(y_test.tolist()),y_pred))
print("Precision",precision_score(np.asarray(y_test.tolist()),y_pred))
print("Recall",recall_score(np.asarray(y_test.tolist()),y_pred))
print("F1-Score",f1_score(np.asarray(y_test.tolist()),y_pred))
tn, fp, fn, tp = confusion_matrix(np.asarray(y_test.tolist()), y_pred).ravel()

all_X = new_df['feature_vector'].tolist()
all_Y_pred = clf.predict(all_X)
new_df['political']=list(all_Y_pred)

tmp = combined_df[combined_df['hard_moderacy']!=-100]
X_new=tmp['feature_vector']
Y_new=tmp['hard_moderacy']

X_train, X_test, y_train, y_test = train_test_split(X_new, Y_new,test_size=0.33,shuffle=True,random_state=23)

print()
print("Moderacy-Predictions")
clf = LogisticRegression(random_state=23,max_iter=700).fit(np.asarray(X_train.tolist()), y_train.tolist())
y_pred=clf.predict(np.asarray(X_test.tolist()))
print("Accuracy",accuracy_score(np.asarray(y_test.tolist()),y_pred))
print("Precision",precision_score(np.asarray(y_test.tolist()),y_pred))
print("Recall",recall_score(np.asarray(y_test.tolist()),y_pred))
print("F1-Score",f1_score(np.asarray(y_test.tolist()),y_pred))
tn, fp, fn, tp = confusion_matrix(np.asarray(y_test.tolist()), y_pred).ravel()

all_X = new_df['feature_vector'].tolist()
all_Y_pred = clf.predict(all_X)
new_df['moderacy']=list(all_Y_pred)

del new_df['feature_vector']
new_df.to_csv('../results/Final-Predictions.csv',index=False)


text:
#!/usr/bin/env python

## Alright people, prepare yourself.
## I am going to:

## Populate my mseed files with event and station information
## Write them to SaC files
## Time shift the SAC files (I KNOW RIGHT)
## Remove the instrument response, detrend, remove mean, filter etc...
## ROTATE THE TRACES
## Then leave them all in the "processed" directory
## Good stuff

#### -----This needs to be run from the "raw" directory as given in the ObspyDMT directory structure ----- ####

## import useful packages:

import obspy
import os
from glob import *

## To take in arguments from the terminal
import argparse

# To get relative path
from pathlib import Path

## to only read certain lines in file:
from itertools import islice

##distance baz calculator
from utilities_Array import deg_km_az_baz

# To dictate headers
from obspy.core.util import AttribDict

## Web service packages
from obspy.clients.iris import Client

# TauP
from obspy.taup import TauPyModel

# my model of choice is prem :D
model = TauPyModel(model="prem")

# Make client - just for giggles
client = Client()

import shutil
import numpy

import matplotlib.pyplot as plt

import argparse

parser = argparse.ArgumentParser(description='Preprocessing script for data retrieved from obspy DMT')

parser.add_argument("-fmin","--minimuum_frequency", help="Enter the lowest frequency to be bandpass filtered.", type=float, required=False, action="store", default = 0.05)

parser.add_argument("-fmax","--maximuum_frequency", help="Enter the highest frequency to be bandpass filtered.", type=float, required=False, action="store", default = 1.0)

parser.add_argument("-p","--phase_list", help="Enter the phases you want the travel time predictions for.", type=str, nargs='+', required=True, action="store", default=[])

parser.add_argument("-tap","--taper_percentage", help="Enter the percentage of the trace at each end you want to taper as a decimal.", type=float, required=False, action="store", default=0.05)

parser.add_argument("-rot","--rotate", help="Enter whether you want the traces to be rotated or not.", required=False, action="store_true", default=False)

parser.add_argument("-split","--splitting", help="Enter whether you want the traces stored in a file and be compativle with SHEBA.", type=bool, required=False, action="store", default=False)

parser.add_argument("-detrend","--detrend_T_F", help="Enter whether you want to detrend the trace or not.", type=bool, required=False, action="store", default=True)

parser.add_argument("-v", "--verbose", help="Increase verbosity of output, just --verbose is enough to turn on.",action="store_true")

args = parser.parse_args()

min_freq = args.minimuum_frequency
max_freq = args.maximuum_frequency
phases = args.phase_list[:]
taper_percent = args.taper_percentage
Q_rotate = args.rotate
Q_detrend = args.detrend_T_F
Q_splitting = args.splitting
major=True
print(min_freq, max_freq,phases,taper_percent,Q_rotate,Q_detrend)

### Find the files with the station and event information in them

## Event dir name:
ev_relpath = "EVENTS-INFO/catalog.txt"

## full path of current directory
pwd = os.getcwd()

## path to directory with the Event directory
base_path = str(Path(pwd).parents[1]) + "/"

ev_path = os.path.join(base_path, ev_relpath)

## path to directory with station XML files :D

## gets the dierctory outside the current one:
base_dir = os.path.split(pwd)[0]
station_dir = os.path.join(base_dir,"resp/")
#print(station_dir)

## get event id
ev_name = os.path.basename(base_dir)
#print(ev_name)

## make a list for the stations
station_temp = []

### BEGIN THE LOOP

for M_file in glob('*BH*'):
	print(M_file)
	bits = M_file.split('.')
	print('bits', len(bits))
	station_name=bits[1]
	network_name=bits[0]

### stream the file
	stream = obspy.read(M_file)
## extract the trace
	trace = stream[0]
## copy the trace for no real reason.
	s_time = trace.stats.starttime
	e_time = trace.stats.endtime
	trace2= trace.trim(starttime=s_time)
## open the xml file as an inventory
	inv = obspy.read_inventory("%sSTXML.%s" %(station_dir, M_file))

	network=inv[0]
	station=network[0]
	channel=station[0]

## For each of the channels the station has, if the code (BHZ BHE BHN) is the same
# take the information
# this was the bane of my existance for an hour
	for chnl in station:

		if chnl.code==trace.stats.channel:
			channel_name = chnl.code
			response=channel.response
			dip= chnl.dip
			inclination = chnl.dip + 90

	## take the coordinates of the station
	stla=station.latitude
	stlo=station.longitude
	stel=station.elevation
	station_temp.append(station.code)

	## MORE SANITY CHECKS
	with open(ev_path) as f:
	## only read lines from line 6
	## This will remove all of the nonesense at the start of the file

		for line in islice(f,5,None):
			if not line.startswith('-') and ev_name in line:
				number,event_id,datetime,latitude,longitude,depth,magnitude,magnitude_type,author,flynn_region,mrr,mtt,mpp,mrt,mrp,mtp,stf_func,stf_duration,t1,t2 = line.split(',')
				evla = latitude
				evlo = longitude
				evdp = depth
				mag = magnitude
				event_time = obspy.UTCDateTime(datetime)

	#distances = client.distaz(stalat=stla, stalon=stlo, evtlat=evla, evtlon=evlo)
	dist_deg,dist_km,az,baz = deg_km_az_baz(lat1=float(evla),lon1=float(evlo),lat2=float(stla),lon2=float(stlo))

	#### Want to find the origin time relative to the trace period
	### i.e. we want to find the difference between the event and
	## the start of the trace
	O = obspy.UTCDateTime(event_time) - obspy.UTCDateTime(trace.stats.starttime)
	### I think O is equal to zero...

	### Going to need to loop over all of the phases ###
	## I don't think this is the most efficient way of doing things but I'm ##
	## going to make a list of lists, where i have [phase name, travel time, number (to be put on t) ##
	# good luck my friend #

	ph_time_list = []

#### Populate information into SAC part ####

	trace2.stats.sac = AttribDict({'baz': float(baz), 'az': float(az), 'kstnm': str(station),'cmpinc': inclination, 'gcarc': str(dist_deg),
							'dist': str(dist_km), 'evla': evla, 'evlo': evlo, 'evdp': evdp, 'stla': stla, 'stlo': stlo, 'stel': stel,'o': O,
							'stdp': dip, 'iztype':9, 'mag': mag})

	for t,phase in enumerate(phases):
		time_prediction =  model.get_travel_times(source_depth_in_km=float(evdp), distance_in_degree=float(dist_deg), receiver_depth_in_km=0.0, phase_list=[phase])

		## if the phase is present add the origin time to it
		if time_prediction:
			time_pred_total = O + time_prediction[0].time
		else:
			time_pred_total = None

		print(phase, time_pred_total)
		### add these values to a list ###
		ph_time_list.append([phase, time_pred_total])

		tn = str('t'+str(t+1))
		ktn = str('k'+tn)

		print(tn,ktn)

		# i.e. for SKKS major arc.
		if major == True:
			if len(time_prediction) > 1:
				# only use the first major arc arrival
				time_pred_total_2 = O + time_prediction[1].time
				phase_label2 = "%s_Major" %phase
				ph_time_list.append([phase_label2, time_pred_total_2])

	### for each loop se tthe attributes that I want ###
	# set predicted travel time
	for t,pt_list in enumerate(ph_time_list):
		time_pred = pt_list[1]
		phase_label = pt_list[0]
		tn = str('t'+str(t+1))
		ktn = str('k'+tn)
		setattr(trace2.stats.sac, tn, time_pred)
		setattr(trace2.stats.sac,ktn, phase_label)

	if trace2.stats.channel == "BHZ":
		trace2.stats.sac.cmpaz = 0.0

	print(trace2.stats)

#### process the data, filter etc...

	trace_rr = trace2.copy()

	#question=input("Velocity(v) or Displacement(d)?")

	#trace_rr.remove_response(inventory=inv, output='VEL')

	trace_rr.remove_response(inventory=inv, output='DISP')

	trace_detrend = trace_rr.copy()
	if Q_detrend:
	#### detrend/demean stuff
		trace_detrend.detrend(type='demean')

	#### SO, SO, SO, SO. I'll now filter the trace with the response removed! - Cheers Josh
	#### taper it
	trace_taper = trace_detrend.copy()

	# This gets the same results as the SAC function taper.
	trace_taper.taper(max_percentage=taper_percent) ## Defaults to both sides and 'hann' type.

	trace_filt = trace_taper.copy()

	#trace_filt.filter('bandpass', freqmin=0.05, freqmax=0.3, corners=4, zerophase=True)

	trace_filt.filter('bandpass', freqmin=min_freq, freqmax=max_freq, corners=4, zerophase=True)

	trace_filt.write("%s.%s.%s.SAC" %(network_name,station_name,channel_name), format="SAC")

	print("------------Next file -------------")

## Time shift the SAC file
import obspy.io.sac.sactrace

for sac_file in glob('*SAC'):

	sac=obspy.io.sac.sactrace.SACTrace.read(sac_file)
	FILE_NAME = os.path.splitext(sac_file)
	print(FILE_NAME[0])
	sac.reftime= event_time

	print(sac.t1, sac.b, sac.o)

	tr2 = sac.to_obspy_trace()
	print(tr2.stats)
	tr2.write("%s.SAC" %FILE_NAME[0], format="SAC")

print("YAY - It worked, I now have SAC files... WITH POPULATED INFORMATION")

print("I hope this is living up to your expectations")
print("If not, you can leave and write your own script")
print("Then send me it.")
print("It'll be better and more efficient than mine...")

## make a list with only no repeated station name

stations = list(set(station_temp))

print(stations)

# Q_rotate = True
#### ROTATE! ####
print("rotate:", Q_rotate)
if Q_rotate == True:

	## Now we rotate ## if the user wants to ##

	## use the lists made at the start of the code

	for station in stations:
	## print station, because, why not?
		print(station)
		## make temporary list for the sac files
		temp_list = []

		## loop over all the sac files, because efficiency is for losers
		for sac_file in glob('*SAC'):
		## print the name, because, see above comment about efficiency and losers
			## take the BHN and BHE files and shove them in the temporary list
			if "%s." %station in sac_file and 'BHN' in sac_file:
				temp_list.append(sac_file)
			if "%s." %station in sac_file and 'BHE' in sac_file:
				temp_list.append(sac_file)
			#print(temp_list)

	## make a stream of those two lists only
		stream = obspy.read(temp_list[0])
		stream += obspy.read(temp_list[1])

		## find back azimuth and inclination
		baz = stream[0].stats.sac.baz
		inc = stream[0].stats.sac.cmpinc

	### check if they are the same length, if not, trim the longer one by the difference in data points

		npts1 = stream[0].stats.npts
		npts2 = stream[1].stats.npts
		samp_rate = stream[0].stats.delta

		s_time1 = stream[0].stats.starttime
		s_time2 = stream[1].stats.starttime

		e_time1 = stream[0].stats.endtime
		e_time2 = stream[1].stats.endtime

		# print(stream)

		## correct for different trace start times

		if s_time1 != s_time2:
			if s_time1 > s_time2:
				stream[0].trim(starttime=s_time1)
				stream[1].trim(starttime=s_time1)
			if s_time1 < s_time2:
				stream[0].trim(starttime=s_time2)
				stream[1].trim(starttime=s_time2)
		if e_time1 != e_time2:
			if e_time1 > e_time2:
				stream[0].trim(endtime=e_time2)
				stream[1].trim(endtime=e_time2)
			if e_time1 < e_time2:
				stream[0].trim(endtime=e_time1)
				stream[1].trim(endtime=e_time1)

		stream.rotate(method = 'NE->RT', back_azimuth = baz, inclination = inc)

	### Get the traces from the stream ###
		trace1 = stream[0]
		trace2 = stream[1]

		s_time3 = trace1.stats.starttime
		s_time4 = trace2.stats.starttime

		e_time3 = trace1.stats.endtime
		e_time4 = trace2.stats.endtime

		## correct for different trace start times

		if s_time3 != s_time4:
			if s_time3 > s_time4:
				trace1.trim(starttime=s_time3)
				trace2.trim(starttime=s_time3)
			if s_time3 < s_time4:
				trace1.trim(starttime=s_time4)
				trace2.trim(starttime=s_time4)
		if e_time3 != e_time4:
			if e_time3 > e_time4:
				trace1.trim(endtime=e_time4)
				trace2.trim(endtime=e_time4)
			if e_time1 < e_time2:
				trace1.trim(endtime=e_time3)
				trace2.trim(endtime=e_time3)

		print(stream)

		ch1 = trace1.stats.channel
		ch2 = trace2.stats.channel

		print(ch1,ch2)

		if ch1[-1] == "R":
			cmpaz_1 = float(trace1.stats.sac.baz) + 180
		elif ch1[-1] == "T":
			cmpaz_1 = float(trace1.stats.sac.baz) + 270
		else:
			pass

		if ch2[-1] == "R":
			cmpaz_2 = float(trace1.stats.sac.baz) + 180
		elif ch2[-1] == "T":
			cmpaz_2 = float(trace1.stats.sac.baz) + 270
		else:
			pass

		### write component information into the file
		trace1.stats.sac['cmpaz'] = float(cmpaz_1)
		trace2.stats.sac['cmpaz'] = float(cmpaz_2)

		print("CMPAZ: ", trace1.stats.sac.cmpaz)
	###### get info from the name of the first trace #####
		network_name, station_name,channel_name,SAC = temp_list[0].split(".")

	## write the dude ##

		trace1.write("%s.%s.%s.SAC" %(str(trace1.stats.network), station_name,str(trace1.stats.channel)), format = "SAC")
		trace2.write("%s.%s.%s.SAC" %(str(trace2.stats.network), station_name,str(trace2.stats.channel)), format = "SAC")

	## Store the BHN and BHE in case I want them for whatever reason.
	# make directory
	H_E_Dir="../processed/BHN_BHE_processed"
	if not os.path.exists(H_E_Dir):
	    os.makedirs(H_E_Dir)
	# move files
	for old_file in glob('*BHN*.SAC'):
		if os.path.exists(str(H_E_Dir) + "/" + str(old_file)):
			os.remove(str(H_E_Dir) + "/" + str(old_file))
		shutil.move(old_file,H_E_Dir)

	for old_file in glob('*BHE*.SAC'):
		if os.path.exists(str(H_E_Dir) + "/" + str(old_file)):
			os.remove(str(H_E_Dir) + "/" + str(old_file))
			shutil.move(old_file,H_E_Dir)

## should be left with only BHR/BHT/BHZ files

### We are, in fact, not done...
## need to move the files to the "processed" directory

if Q_splitting:
### Going to do my splitting work within a seperate directory
	Splitting_Dir = "../Splitting"
	if not os.path.exists(Splitting_Dir):
	    os.makedirs(Splitting_Dir)
else:
	pass

for pro_file in glob('*SAC'):
	pro_file_First = pro_file.split(".")[0]
	pro_file_Second = pro_file.split(".")[1]
	pro_file_no_SAC = pro_file_First + "." + pro_file_Second
	if os.path.exists("../processed/%s" %pro_file):
		os.remove("../processed/%s" %pro_file)
	shutil.copy(pro_file, "../processed")

	if Q_splitting:
		if os.path.exists("%s/%s" %(Splitting_Dir,pro_file)):
			os.remove("%s/%s" %(Splitting_Dir,pro_file))
		shutil.move(pro_file, Splitting_Dir)
		os.rename("%s/%s" %(Splitting_Dir, pro_file), "%s/%s" %(Splitting_Dir,pro_file_no_SAC))
	else:
		pass

max_stars_repo_name: eejwa/Data-Preprocessing

max_issues_repo_name: eejwa/Data-Preprocessing

max_forks_repo_name: eejwa/Data-Preprocessing


text:
""" Module to correct pulsar and FRB DMs for the MW ISM """

from ne2001 import ne_io, density #ne2001 ism model
import pygedm #ymw ism model
import numpy as np
import pandas as pd
from astropy import units as u
from astropy.coordinates import SkyCoord, Galactic
import logging

logging.basicConfig(format='%(asctime)s - %(message)s', datefmt='%d-%b-%y %H:%M:%S', level=logging.INFO)
ne = density.ElectronDensity()

def find_delta_dm(transient_type,transient_data,ism_model,b_val,mc_deg=5,save_df=True):
    """
    Find pulsar/FRB DMs corrected for by the MW ISM DM and remove observations in complex DM regions.
    Returns array of DMs
    FRB data is available as a csv in the FRBs/FRB/frb/data/FRBs repo (FRB catalogue [Petroff et al. 2017])
    Pulsar data is avaiable as a csv in the FRBs/pulsars/pulsars/data/atnf_cat repo (v1.61 ATNF pulsar catalogue [Manchester et al. 2005])

    Arguments:
        transient_type (str):
            Accepts 'frb' or 'pulsar'.
        transient_data (str):
            Path to data (in .csv format).
        ism_model (str):
            Model used to calculated the MW halo DM.
            Accepts 'ymw16'  [Yao et al. 2017] or 'ne2001' [Cordes & Lazio 2003].
        b_val (int):
            Galactic latitude considered (b>b_val, b<-b_val).
        mc_deg (int):
            Number of degrees from Magellanic clouds within which transients are removed.
        save_df (str, optional):
            Save transient DMs and coords to csv.

    Outputs:
    """
    # Sort data and get coords
    if transient_type=='frb':
        transcat_df = pd.read_csv(transient_data, skiprows=1, usecols= [0,5,6,7], names=['Name','l','b','dm'])
        transcat_df['dm'] = transcat_df['dm'].str.split('&').str[0].astype(float).values
        coords = SkyCoord(l=transcat_df['l'], b=transcat_df['b'], unit=(u.degree),frame=Galactic)
    elif transient_type=='pulsar':
        transcat_df = pd.read_csv(transient_data, skiprows=2, usecols = [1,2,3,9,10], names=['Name','Pref','dm','RAJD','DECJD'])
        transcat_df = transcat_df[~transcat_df['dm'].str.contains('*', regex=False)].reset_index(drop=True)
        transcat_df['dm'] = transcat_df['dm'].astype(float)
        c_icrs = SkyCoord(ra=transcat_df['RAJD'], dec=transcat_df['DECJD'], unit=(u.degree), frame='icrs')
        transcat_df['l'] = pd.DataFrame(c_icrs.galactic.l.value)
        transcat_df['b'] = pd.DataFrame(c_icrs.galactic.b.value)
        coords = SkyCoord(l=transcat_df['l'], b=transcat_df['b'], unit=(u.degree),frame=Galactic)

    # Find transients in line of sight of MCs
    logging.info('Removing transients near Magellanic clouds...')
    # LMC
    lmc_distance = 50*u.kpc
    lmc_coord = SkyCoord('J052334.6-694522',unit=(u.hourangle, u.deg),distance=lmc_distance)
    close_to_lmc = lmc_coord.separation(coords) < mc_deg*u.deg
    lmc_trans = list(transcat_df[close_to_lmc]['Name'])
    # SMC
    smc_distance = 61*u.kpc
    smc_coord = SkyCoord('J005238.0-724801',unit=(u.hourangle, u.deg),distance=smc_distance)
    close_to_smc = smc_coord.separation(coords) < mc_deg*u.deg
    smc_trans = list(transcat_df[close_to_smc]['Name'])

    transcat_df = transcat_df[~transcat_df['Name'].isin(lmc_trans)].reset_index(drop=True)
    transcat_df = transcat_df[~transcat_df['Name'].isin(smc_trans)].reset_index(drop=True)
    if transient_type=='pulsar':
        transcat_df = transcat_df[~transcat_df['Pref'].str.contains('mfl+06', regex=False)].reset_index(drop=True)
    elif transient_type=='frb':
        pass

    # Remove transients with low Galactic lattitudes
    logging.info('Removing transients with low Galactic lattitudes...')
    transcat_df = pd.concat([transcat_df[transcat_df.b > b_val], transcat_df[transcat_df.b < -b_val]], ignore_index=True)

    # ISM model
    logging.info('Correcting transient DMs for ISM...')
    trans_ism = []
    if ism_model=='ymw16':
        for i in range(len(transcat_df['dm'])):
            trans_ism_ = pygedm.dist_to_dm(transcat_df['l'].iloc[i], transcat_df['b'].iloc[i], 100000)[0].value
            trans_ism = np.append(trans_ism,trans_ism_)
    elif ism_model=='ne2001':
        for i in range(len(transcat_df['dm'])):
            trans_ism_ = ne.DM(transcat_df['l'].iloc[i], transcat_df['b'].iloc[i], 100.).value
            trans_ism = np.append(trans_ism,trans_ism_)

    transcat_df['trans_ism'] = pd.DataFrame(trans_ism)
    transcat_df['deltaDM'] = pd.DataFrame(transcat_df['dm']-transcat_df['trans_ism'])

    if save_df==True:
        transcat_df.to_csv('transient_data/'+transient_type+'cat_df_'+ism_model+'_'+str(int(b_val))+'.csv')
        logging.info('Transient data saved to csv.')
    else:
        pass

    return np.array(transcat_df['deltaDM'])

max_stars_repo_name: Lachimax/FRB

max_issues_repo_name: Lachimax/FRB

max_forks_repo_name: Lachimax/FRB

reason: import numpy,from astropy


text:
from sklearn import datasets
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

X, y = datasets.load_iris(return_X_y=True, as_frame=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42, test_size=0.3)


text:
#######################################################################
# Plotting script for Wetropolis Au dynamics - test cases
#######################################################################

'''
Plots simulation data from main_wetro_tests.py

Specify in spec below:
> config#0 for steady state test
> config#1 for floodwave test

Output: interactive video of simulations.
'''

##################################################################
# GENERIC MODULES REQUIRED
##################################################################
import numpy as np
import scipy as sp
import os
import errno
import sys
import importlib.util
import matplotlib.pyplot as plt
from matplotlib import animation
import animatplot as amp # see ex_animatplot.py for example script.

##################################################################
# CUSTOM MODULES REQUIRED
##################################################################
from cross_sections import xsec_hAs, xsec_Ahs, plot_xsec_hAs

##################################################################
# IMPORT PARAMETERS FROM CONFIGURATION FILE
##################################################################
#spec = importlib.util.spec_from_file_location("config", sys.argv[1])
spec = importlib.util.spec_from_file_location("config","configs/config#0.py")
config = importlib.util.module_from_spec(spec)
spec.loader.exec_module(config)

## config pars
hr = config.hr
wr = config.wr
hf = config.hf
hc = hr+hf
wf = config.wf
wc = config.wc
tana = config.tana
LR1 = config.LR1
LR2 = config.LR2
LR3 = config.LR3
LR11 = config.LR11
LR22 = config.LR22
tr = config.tr
Nk = config.Nk
s_r = config.s_r
s_m = config.s_m
dbds = config.dbds
g = config.g
Cm = config.Cm
Neq = config.Neq
ic = config.ic
cfl = config.cfl
BC = config.BC

##################################################################
# Set up dirs
##################################################################
outdir = config.outdir
cwd = os.getcwd()
dirn = str(cwd+outdir)
try:
    os.makedirs(dirn)
except OSError as exception:
    if exception.errno != errno.EEXIST:
        raise

##################################################################
# Set up grid
##################################################################

## set up grid
L=LR3 #length of domain
# Nk=25*L #number of gridcells (excluding ghost)
# Nk = int(Nk)
Nf=Nk+1 #number of nodes
Kk=L/Nk #length of cell

s = np.linspace(0, L, Nk+1)
sBC = np.linspace(-Kk, L+Kk, Nk+3)  #node loc with ghosts

# locating floodplain/city
index_fp = np.where((s < LR1) | (s > LR2))
index_fp = np.array(index_fp)[0]
index_city = np.where((s >= LR1) & (s <= LR2))
index_city = np.array(index_city)[0]

##################################################################
# Define time parameters
##################################################################

tn = config.tn
wd = config.wd #wetropolis day
tmax = config.tmax
Nmeas = config.Nmeas
dtmeasure = tmax/Nmeas
tmeasure = dtmeasure
timevec = np.linspace(tn,tmax,Nmeas+1)

print(' Loading simulation data from:', dirn)

U_array = np.load(str(dirn+'/U_array.npy'))
h_array = np.load(str(dirn+'/h_array.npy'))
Z_array = np.load(str(dirn+'/Z_array.npy'))

fp = index_fp[50]
ct = index_city[5]

##################################################################
# Animation using animatplot
##################################################################

h =  h_array[:,:,:][0]
A =  U_array[0,:,:]
Au =  U_array[1,:,:]
fp = index_fp[50]
ct = index_city[5]

fig, axes = plt.subplots(2, 1, figsize=(10,6))

axes[0].plot([s[fp], s[fp]],[0,0.04],'k:')
axes[0].plot([s[ct], s[ct]],[0,0.04],'k:')
axes[0].fill([config.LR1, config.LR2,config.LR2,config.LR1],[0,0,config.hc,config.hc],'r',alpha=0.1,linestyle='None')
axes[0].set_ylim([0,0.03])
axes[0].set_xlim([0,L])
axes[0].set_ylabel('$h(s,t)$',fontsize=14)
# axes[0].set_xlabel('$s$',fontsize=14)

axes[1].set_ylim([0,0.0006])
axes[1].set_xlim([0,L])
axes[1].set_ylabel('$Au(s,t)$',fontsize=14)
axes[1].set_xlabel('$s$',fontsize=14)

hline_block = amp.blocks.Line(s,h[1:,:],ax=axes[0], t_axis=1, linewidth = 1.0, drawstyle='steps-post')

Auline_block = amp.blocks.Line(s,Au[1:,:],ax=axes[1], t_axis=1, linewidth = 1.0, drawstyle='steps-post')

timeline = amp.Timeline(timevec, fps=10)

# now to construct the animation
anim = amp.Animation([hline_block, Auline_block], timeline)
anim.controls()

fname = dirn+'/simulation'
anim.save_gif(fname)

plt.show()

reason: import numpy,import scipy


text:
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

from bias_detector.BiasConfidence import BiasConfidence
from bias_detector.BiasMetric import get_bias_metric_description, get_bias_metric, get_bias_metric_short_description
from bias_detector.common import get_p_others, get_other_groups_names, races

class BiasReport:

    def __init__(self,
                 bias_metrics_results: pd.DataFrame,
                 estimated_groups_sizes: pd.Series,
                 y_true: pd.Series,
                 y_pred: pd.Series,
                 y_scores: pd.Series,
                 privileged_race: str,
                 p_groups: pd.DataFrame,
                 detect_gender_bias: bool,
                 detect_race_bias: bool) -> None:
        self.bias_metrics_results = bias_metrics_results
        self.estimated_groups_sizes = estimated_groups_sizes
        self.y_true = y_true
        self.y_pred = y_pred
        self.y_scores = y_scores
        self.privileged_race = privileged_race
        self.p_groups = p_groups
        self.detect_gender_bias = detect_gender_bias
        self.detect_race_bias = detect_race_bias

    def plot_groups(self) -> None:
        if self.detect_race_bias:
            race_group_sizes = self.estimated_groups_sizes[races + ['race_unknown']]
            race_group_percentages = pd.Series(race_group_sizes.values / race_group_sizes.sum()).apply(lambda fraction: '{fraction:.2%}'.format(fraction=fraction))
            if len(race_group_sizes) > 0:
                plt.title('Race estimated groups sizes')
                plt.pie(race_group_sizes.values)
                plt.legend(labels=race_group_sizes.index.values + ' - ' + race_group_percentages)
                plt.show()
        if self.detect_gender_bias:
            gender_group_sizes = self.estimated_groups_sizes[['male', 'female', 'gender_unknown']]
            gender_group_percentages = pd.Series(gender_group_sizes.values / gender_group_sizes.sum()).apply(lambda fraction: '{fraction:.2%}'.format(fraction=fraction))
            if len(gender_group_sizes) > 0:
                plt.title('Gender estimated groups sizes')
                plt.pie(gender_group_sizes.values)
                plt.legend(labels=gender_group_sizes.index.values + ' - ' + gender_group_percentages)
                plt.show()

    def plot_y_pred(self) -> None:
        for group_name in self.bias_metrics_results.columns:
            if group_name == 'female':
                continue
            group = self.p_groups[group_name]
            others = get_p_others(self.p_groups, group_name, self.privileged_race)
            group_title, others_title = self.get_titles(group_name)
            if self.y_pred is not None:
                plt.title('y_pred distribution by {group_title}/{others_title} classes'.format(
                    group_title=group_title, others_title=others_title))
                pos_pred_group_size = (self.y_pred == 1).multiply(group).sum()
                pos_pred_others_size = (self.y_pred == 1).multiply(others).sum()
                neg_pred_group_size = (self.y_pred == 0).multiply(group).sum()
                neg_pred_others_size = (self.y_pred == 0).multiply(others).sum()
                plt.bar(x=[0, 1], height=[neg_pred_group_size/group.sum(), pos_pred_group_size/group.sum()], label=group_title, alpha=0.5, width=0.1)
                plt.bar(x=[0, 1], height=[neg_pred_others_size/others.sum(), pos_pred_others_size/others.sum()], label=others_title, alpha=0.5, width=0.1)
                plt.legend()
                plt.xlabel('y_pred')
                plt.ylabel('density')
                plt.show()

    def plot_y_scores(self) -> None:
        for group_name in self.bias_metrics_results.columns:
            if group_name == 'female':
                continue
            group = self.p_groups[group_name]
            others = get_p_others(self.p_groups, group_name, self.privileged_race)
            group_title, others_title = self.get_titles(group_name)
            if self.y_scores is not None:
                plt.title('y_scores distribution by {group_label}/{others_label} classes'.format(
                    group_label=group_title, others_label=others_title))
                bucket_count = 20
                bucket_size = 1 / bucket_count
                buckets = [i * bucket_size for i in range(0, bucket_count)]
                normalized_group_sizes = []
                normalized_others_sizes = []
                for bucket in buckets:
                    bucket_scores = ((self.y_scores >= bucket) & (self.y_scores < bucket + bucket_size))
                    normalized_group_sizes.append(bucket_scores.multiply(group).sum()/group.sum())
                    normalized_others_sizes.append(bucket_scores.multiply(others).sum()/others.sum())
                plt.bar(x=np.array(buckets) + bucket_size / 2, height=normalized_group_sizes, label=group_title, alpha=0.5, width=bucket_size)
                plt.bar(x=np.array(buckets) + bucket_size / 2, height=normalized_others_sizes, label=others_title, alpha=0.5, width=bucket_size)
                plt.legend()
                plt.xlabel('y_score')
                plt.ylabel('density')
                plt.show()

    def get_titles(self, group_name):
        group_title = group_name.replace('_', ' ').title()
        other_groups_names = get_other_groups_names(group_name, self.privileged_race)
        others_title = other_groups_names[0].replace('_', ' ').title() if len(other_groups_names) == 1 else 'Others'
        group_title = group_title if group_title != 'Api' else 'API'
        others_title = others_title if others_title != 'Api' else 'API'
        return group_title, others_title

    def plot_summary(self) -> dict:
        for bias_metric in self.bias_metrics_results.index:
            bias_metric_title = bias_metric.replace('_', ' ').title()
            results = self.bias_metrics_results.loc[bias_metric]
            x_labels = []
            bar_results = []
            for group_name in results.index:
                if group_name == 'female':
                    continue
                group_title, others_title = self.get_titles(group_name)
                x_labels.append(group_title + '\n' + others_title)
                bar_results.append(results[group_name])

            fig, ax = plt.subplots()
            bar_width = 0.2
            index = np.arange(len(bar_results))
            group_colors = []
            others_colors = []
            for i, r in enumerate(bar_results):
                confidence = r.bias_confidence
                if confidence is None or confidence.p_value > BiasConfidence.alpha:
                    group_colors.append('lightsteelblue')
                    others_colors.append('bisque')
                else:
                    group_colors.append('tab:blue')
                    others_colors.append('tab:orange')
            group_bar = ax.bar(x=index, height=[r.group_p for r in bar_results], width=bar_width, color=group_colors)
            others_bar = ax.bar(x=index+bar_width, height=[r.others_p for r in bar_results], width=bar_width, color=others_colors)
            ax.set_xticks(index + bar_width / 2)
            ax.set_xticklabels(x_labels)
            ax.set_ylabel(get_bias_metric_short_description(get_bias_metric(bias_metric)))
            biased = np.any([r.bias_confidence is not None and r.bias_confidence.p_value <= BiasConfidence.alpha for r in bar_results])
            for i, r in enumerate(bar_results):
                confidence = r.bias_confidence
                if confidence is None or confidence.p_value > BiasConfidence.alpha:
                    continue
                interval = confidence.get_interval()
                plt.plot((i+(bar_width/2), i+(bar_width/2)), (min(r.group_p, r.others_p) + interval/2, max(r.group_p, r.others_p) - interval/2), 'r_-', color='black')
            if biased:
                percents = int((1 - BiasConfidence.alpha) * 100)
                plt.title(bias_metric_title + '\n\n' + 'Intervals represents the {percents}% confidence min difference'.format(percents=percents))
            else:
                plt.title(bias_metric_title + '\n\n' + 'No bias detected')
            plt.xticks(rotation=90)
            plt.show()

    def get_summary(self) -> dict:
        summary = {}
        for bias_metric in self.bias_metrics_results.index:
            groups_alerts = []
            bias_metric_formatted = bias_metric.replace('_', ' ').title()
            summary[bias_metric_formatted] = groups_alerts
            results = self.bias_metrics_results.loc[bias_metric]
            for group_name in results.index:
                if group_name == 'female':
                    continue
                result = results[group_name]
                confidence = result.bias_confidence
                if result is not None and confidence is not None and type(confidence) == BiasConfidence and confidence.p_value <= BiasConfidence.alpha:
                    group_title, others_title = self.get_titles(group_name)
                    group_metric_description = get_bias_metric_description(get_bias_metric(bias_metric), group_title)
                    others_metric_description = get_bias_metric_description(get_bias_metric(bias_metric), others_title)
                    group_alert = '{group_metric_description}-{others_metric_description}={group_p:.2}-{others_p:.2}={diff:.2}±{ci:.2} (p-value={p_value:.2})'.format(
                        group_metric_description=group_metric_description,
                        others_metric_description=others_metric_description,
                        diff=result.get_diff(),
                        ci=confidence.get_interval(),
                        group_p=result.group_p,
                        others_p=result.others_p,
                        p_value=confidence.p_value)
                    groups_alerts.append(group_alert)
        return summary

    def print_summary(self):
        from IPython.core.display import display, HTML
        display(HTML(self.get_summary_html()))

    def get_summary_html(self):
        summary = self.get_summary()
        formatted_summary = '''Report Summary:'''
        formatted_summary += '''
'''
        for bias_metric_title, groups_alerts in summary.items():
            formatted_summary += '''
{bias_metric_title}:
'''.format(bias_metric_title=bias_metric_title)
            if len(groups_alerts) == 0:
                formatted_summary += '''No bias detected.'''
            else:
                formatted_summary += 'We observed the following statistically significant differences (𝛼={alpha:.2}):'\
                    .format(alpha=BiasConfidence.alpha)
                formatted_summary += '''
'''
                for group_alert in groups_alerts:
                    formatted_summary += '''
{group_alert}
'''.format(group_alert=group_alert)
                formatted_summary += '''
'''
        formatted_summary += '''
'''
        return formatted_summary

    def get_bias_metrics_results(self):
        return self.bias_metrics_results

    def get_estimated_groups_sizes(self):
        return self.estimated_groups_sizes

    def get_p_groups(self):
        return self.p_groups

max_stars_repo_name: yairhoresh/bias-detector

max_issues_repo_name: yairhoresh/bias-detector


text:
#CUDA Code
from pycuda.compiler import SourceModule
gpuSource = SourceModule("""
#include <pycuda-complex.hpp>
#include <stdio.h>

typedef pycuda::complex<double> dcmplx;

const int vectorSize = 10;
const int spinComps = 5;
const double pi = 3.14159265358979323846264338328;

/* QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize] */

__global__ void test(dcmplx *QField, dcmplx *vortField, dcmplx *rhoField, int *lattice)
{

}

__device__ __forceinline__ double phase (dcmplx z)
{
    double R = atan2(z.imag(), z.real());
    return R + pi;

}

__device__ __forceinline__ double phaseSubtract (double phase1, double phase2)
{
    bool edgeCase1 = (phase1 < pi/2. && phase2 > 3.*pi/2.);
    bool edgeCase2 = (phase2 < pi/2. && phase1 > 3.*pi/2.);
    double phaseDiff = phase1 - phase2;
    if (edgeCase1){
      return 2.*pi + phaseDiff;
    } else if (edgeCase2){
      return phaseDiff - 2.*pi;
    } else {
      return phaseDiff;
    }

}

__global__ void aveVorticity(dcmplx *vortField, dcmplx *vortFieldAve, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int xL = (x-1+xSize)%(xSize);
  int xR = (x+1)%(xSize);
  int yL = (y-1+ySize)%(ySize);
  int yR = (y+1)%(ySize);
  vortFieldAve[y + x*ySize] =  1./9.*(vortField[yL + xL*ySize] +
                                      vortField[yL + xR*ySize] +
                                      vortField[yL + x*ySize] +
                                      vortField[yR + xL*ySize] +
                                      vortField[yR + xR*ySize] +
                                      vortField[yR + x*ySize] +
                                      vortField[y + xL*ySize] +
                                      vortField[y + xR*ySize] +
                                      vortField[y + x*ySize]);
}

__global__ void getPlotDetailsPhaseAndDensity(dcmplx *QField, dcmplx *vortField, dcmplx *rhoField, dcmplx *bosonField, double *phaseField, dcmplx *VxField, dcmplx *VyField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  /* For Vorticity */
  double phiRightAbove[ spinComps ];
  double phiRightBelow[ spinComps ];
  double phiLeftAbove[ spinComps ];
  double phiLeftBelow[ spinComps ];

  /* For Velocity */
  double phiRight[ spinComps ];
  double phiLeft[ spinComps ];
  double phiAbove[ spinComps ];
  double phiBelow[ spinComps ];

  double Vx = 0.;
  double Vy = 0.;
  double VxAbove = 0.;
  double VxBelow = 0.;
  double VyLeft = 0.;
  double VyRight = 0.;
  dcmplx i(0.,1.);
  int xLeft = (x-1+xSize)%(xSize);
  int xRight = (x+1)%(xSize);
  int yBelow = (y-1+ySize)%(ySize);
  int yAbove = (y+1)%(ySize);

  phaseField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  for(n = 0; n < spinComps; n++){
    phi[n] = QField[2*n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];

    /* For Vorticity */
    phiRightAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiRightBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeftAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiLeftBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);

    /* For Velocity */
    phiRight[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeft[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]);
    phiBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]);

    /* For Phase */
    phaseField[y+x*ySize] += phase(phi[n]);
  }

  while(phaseField[y+x*ySize] > 2.*pi){
    phaseField[y+x*ySize] -= 2.*pi;
  }

  rhoField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }
  bosonField[y+x*ySize] = 0;
  for(n = 0; n < spinComps; n++){
    bosonField[y+x*ySize] += phi[n]*conj(phi[n]);
  }
  for(n = 0; n < spinComps; n++){
      Vx += (1./2.)*phaseSubtract(phiRight[n], phiLeft[n]);
      Vy += (1./2.)*phaseSubtract(phiAbove[n], phiBelow[n]);
      VxAbove += (1./2.)*phaseSubtract(phiRightAbove[n], phiLeftAbove[n]);
      VxBelow += (1./2.)*phaseSubtract(phiRightBelow[n], phiLeftBelow[n]);
      VyLeft += (1./2.)*phaseSubtract(phiLeftAbove[n], phiLeftBelow[n]);
      VyRight += (1./2.)*phaseSubtract(phiRightAbove[n], phiRightBelow[n]);
  }
  VxField[z+y*zSize+x*zSize*ySize] = Vx;
  VyField[z+y*zSize+x*zSize*ySize] = Vy;
  vortField[z+y*zSize+x*zSize*ySize] = (1./2.)*((VyRight - VyLeft)-(VxAbove-VxBelow));

}

__global__ void getPlotDetailsVorticity(dcmplx *QField, dcmplx *vortField, dcmplx *rhoField, dcmplx *bosonField, double *phaseField, dcmplx *VxField, dcmplx *VyField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  /* For Vorticity */
  double phiRightAbove[ spinComps ];
  double phiRightBelow[ spinComps ];
  double phiLeftAbove[ spinComps ];
  double phiLeftBelow[ spinComps ];

  /* For Velocity */
  double phiRight[ spinComps ];
  double phiLeft[ spinComps ];
  double phiAbove[ spinComps ];
  double phiBelow[ spinComps ];

  double Vx [ spinComps ];
  double Vy [ spinComps ];
  double VxAbove [ spinComps ];
  double VxBelow [ spinComps ];
  double VyLeft [ spinComps ];
  double VyRight [ spinComps ];
  dcmplx i(0.,1.);
  int xLeft = (x-1+xSize)%(xSize);
  int xRight = (x+1)%(xSize);
  int yBelow = (y-1+ySize)%(ySize);
  int yAbove = (y+1)%(ySize);

  phaseField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  for(n = 0; n < spinComps; n++){
    phi[n] = QField[2*n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];

    /* For Vorticity */
    phiRightAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiRightBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeftAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiLeftBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);

    /* For Velocity */
    phiRight[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeft[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]);
    phiBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]);

    /* For Phase */
    phaseField[y+x*ySize] += phase(phi[n]);
  }

  while(phaseField[y+x*ySize] > 2.*pi){
    phaseField[y+x*ySize] -= 2.*pi;
  }

  rhoField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }
  for(n = 0; n < spinComps; n++){
  bosonField[n+y*spinComps+x*ySize*spinComps] = 0;
  }
  for(n = 0; n < spinComps; n++){
    bosonField[n+y*spinComps+x*ySize*spinComps] = phi[n]*conj(phi[n]);
  }
  for(n = 0; n < spinComps; n++){
      Vx[n]= (1./2.)*phaseSubtract(phiRight[n], phiLeft[n]);
      Vy[n]= (1./2.)*phaseSubtract(phiAbove[n], phiBelow[n]);
      VxAbove[n] += (1./2.)*phaseSubtract(phiRightAbove[n], phiLeftAbove[n]);
      VxBelow[n]= (1./2.)*phaseSubtract(phiRightBelow[n], phiLeftBelow[n]);
      VyLeft[n]= (1./2.)*phaseSubtract(phiLeftAbove[n], phiLeftBelow[n]);
      VyRight[n]= (1./2.)*phaseSubtract(phiRightAbove[n], phiRightBelow[n]);
  }
  VxField[z+y*zSize+x*zSize*ySize] = Vx[0];
  VyField[z+y*zSize+x*zSize*ySize] = Vy[0];
  for (n = 0; n < spinComps; n++){
    vortField[n + y*spinComps+x*spinComps*ySize] = (1./2.)*((VyRight[n] - VyLeft[n])-(VxAbove[n]-VxBelow[n]));
  }
}

__global__ void getPlotDetailsMayavi(dcmplx *QField, dcmplx *rhoField, double *phaseField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int comp = lattice[6];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  dcmplx i(0.,1.);

  phaseField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  for(n = 0; n < spinComps; n++){
    phi[n] = QField[2*n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }
  /* For Phase */
  phaseField[y+x*ySize] = phase(phi[comp]);

  while(phaseField[y+x*ySize] > 2.*pi){
    phaseField[y+x*ySize] -= 2.*pi;
  }

  rhoField[y+x*ySize] = psi[2*comp]*conj(psi[2*comp]) + psi[2*comp+1]*conj(psi[2*comp+1]);

}

__global__ void getPlotDetails(dcmplx *QField, dcmplx *vortField, dcmplx *rhoField, double *phaseField, dcmplx *VxField, dcmplx *VyField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  /* For Vorticity */
  double phiRightAbove[ spinComps ];
  double phiRightBelow[ spinComps ];
  double phiLeftAbove[ spinComps ];
  double phiLeftBelow[ spinComps ];

  /* For Velocity */
  double phiRight[ spinComps ];
  double phiLeft[ spinComps ];
  double phiAbove[ spinComps ];
  double phiBelow[ spinComps ];

  double Vx = 0.;
  double Vy = 0.;
  double VxAbove = 0.;
  double VxBelow = 0.;
  double VyLeft = 0.;
  double VyRight = 0.;
  dcmplx i(0.,1.);
  int xLeft = (x-1+xSize)%(xSize);
  int xRight = (x+1)%(xSize);
  int yBelow = (y-1+ySize)%(ySize);
  int yAbove = (y+1)%(ySize);

  phaseField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  for(n = 0; n < spinComps; n++){
    phi[n] = QField[2*n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];

    /* For Vorticity */
    phiRightAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiRightBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeftAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiLeftBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);

    /* For Velocity */
    phiRight[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeft[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]);
    phiBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]);

    /* For Phase */
    phaseField[y+x*ySize] += phase(phi[n]);
  }

  while(phaseField[y+x*ySize] > 2.*pi){
    phaseField[y+x*ySize] -= 2.*pi;
  }

  rhoField[y+x*ySize] = 0;
  for(n = 0; n < vectorSize; n++){
    rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }

  for(n = 0; n < spinComps; n++){
      Vx += (1./2.)*phaseSubtract(phiRight[n], phiLeft[n]);
      Vy += (1./2.)*phaseSubtract(phiAbove[n], phiBelow[n]);
      VxAbove += (1./2.)*phaseSubtract(phiRightAbove[n], phiLeftAbove[n]);
      VxBelow += (1./2.)*phaseSubtract(phiRightBelow[n], phiLeftBelow[n]);
      VyLeft += (1./2.)*phaseSubtract(phiLeftAbove[n], phiLeftBelow[n]);
      VyRight += (1./2.)*phaseSubtract(phiRightAbove[n], phiRightBelow[n]);
  }
  VxField[z+y*zSize+x*zSize*ySize] = Vx;
  VyField[z+y*zSize+x*zSize*ySize] = Vy;
  vortField[z+y*zSize+x*zSize*ySize] = (1./2.)*((VyRight - VyLeft)-(VxAbove-VxBelow));

}

__global__ void getPlotDetailsDS(dcmplx *QField, dcmplx *vortField, dcmplx *rhoField, double *phaseField, dcmplx *VxField, dcmplx *VyField, dcmplx *rhoFieldComp, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  /* For Vorticity */
  double phiRightAbove[ spinComps ];
  double phiRightBelow[ spinComps ];
  double phiLeftAbove[ spinComps ];
  double phiLeftBelow[ spinComps ];

  /* For Velocity */
  double phiRight[ spinComps ];
  double phiLeft[ spinComps ];
  double phiAbove[ spinComps ];
  double phiBelow[ spinComps ];

  double Vx = 0.;
  double Vy = 0.;
  double VxAbove = 0.;
  double VxBelow = 0.;
  double VyLeft = 0.;
  double VyRight = 0.;
  dcmplx i(0.,1.);
  int xLeft = (x-1+xSize)%(xSize);
  int xRight = (x+1)%(xSize);
  int yBelow = (y-1+ySize)%(ySize);
  int yAbove = (y+1)%(ySize);

  phaseField[y+x*ySize] = 0;

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  for(n = 0; n < spinComps; n++){
    phi[n] = QField[2*n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];

    /* For Vorticity */
    phiRightAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiRightBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeftAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiLeftBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);

    /* For Velocity */
    phiRight[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
    phiLeft[n] = phase(QField[2*n+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
    phiAbove[n] = phase(QField[2*n+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]);
    phiBelow[n] = phase(QField[2*n+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]
              +QField[2*n+1+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]);

    /* For Phase */
    phaseField[y+x*ySize] += phase(phi[n]);
  }

  while(phaseField[y+x*ySize] > 2.*pi){
    phaseField[y+x*ySize] -= 2.*pi;
  }

  rhoField[y+x*ySize] = 0;
  for(n = 0; n < vectorSize; n++){
    rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }

  rhoFieldComp[n+y*vectorSize+x*ySize*vectorSize] = 0;
  for(n = 0; n < vectorSize; n++){
    rhoFieldComp[n+y*vectorSize+x*ySize*vectorSize] += psi[n]*conj(psi[n]);
  }

  for(n = 0; n < spinComps; n++){
      Vx += (1./2.)*phaseSubtract(phiRight[n], phiLeft[n]);
      Vy += (1./2.)*phaseSubtract(phiAbove[n], phiBelow[n]);
      VxAbove += (1./2.)*phaseSubtract(phiRightAbove[n], phiLeftAbove[n]);
      VxBelow += (1./2.)*phaseSubtract(phiRightBelow[n], phiLeftBelow[n]);
      VyLeft += (1./2.)*phaseSubtract(phiLeftAbove[n], phiLeftBelow[n]);
      VyRight += (1./2.)*phaseSubtract(phiRightAbove[n], phiRightBelow[n]);
  }
  VxField[z+y*zSize+x*zSize*ySize] = Vx;
  VyField[z+y*zSize+x*zSize*ySize] = Vy;
  vortField[z+y*zSize+x*zSize*ySize] = (1./2.)*((VyRight - VyLeft)-(VxAbove-VxBelow));

}

__global__ void getPlotDetailsForComponent(dcmplx *QField, dcmplx *vortField, dcmplx *rhoField, dcmplx *phaseField, dcmplx *VxField, dcmplx *VyField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int component = lattice[6];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  /* For Vorticity */
  double phiRightAbove[ spinComps ];
  double phiRightBelow[ spinComps ];
  double phiLeftAbove[ spinComps ];
  double phiLeftBelow[ spinComps ];

  /* For Velocity */
  double phiRight[ spinComps ];
  double phiLeft[ spinComps ];
  double phiAbove[ spinComps ];
  double phiBelow[ spinComps ];

  double Vx = 0.;
  double Vy = 0.;
  double VxAbove = 0.;
  double VxBelow = 0.;
  double VyLeft = 0.;
  double VyRight = 0.;
  dcmplx i(0.,1.);
  int xLeft = (x-1+xSize)%(xSize);
  int xRight = (x+1)%(xSize);
  int yBelow = (y-1+ySize)%(ySize);
  int yAbove = (y+1)%(ySize);

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  phi[component] = QField[2*component+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];

  /* For Vorticity */
  phiRightAbove[component] = phase(QField[2*component+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+yAbove*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
  phiRightBelow[component] = phase(QField[2*component+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+yBelow*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
  phiLeftAbove[component] = phase(QField[2*component+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+yAbove*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
  phiLeftBelow[component] = phase(QField[2*component+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+yBelow*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);

  /* For Velocity */
  phiRight[component] = phase(QField[2*component+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+y*vectorSize*zSize+xRight*zSize*ySize*vectorSize]);
  phiLeft[component] = phase(QField[2*component+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+y*vectorSize*zSize+xLeft*zSize*ySize*vectorSize]);
  phiAbove[component] = phase(QField[2*component+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+yAbove*vectorSize*zSize+x*zSize*ySize*vectorSize]);
  phiBelow[component] = phase(QField[2*component+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+yBelow*vectorSize*zSize+x*zSize*ySize*vectorSize]);

  /* For Phase */
  phaseField[y+x*ySize] = phase(phi[component]);

  rhoField[y+x*ySize] = 0;
  for(n = 2*component; n <= 2*component + 1; n++){
  rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }

  Vx = (1./2.)*phaseSubtract(phiRight[component], phiLeft[component]);
  Vy = (1./2.)*phaseSubtract(phiAbove[component], phiBelow[component]);
  VxAbove = (1./2.)*phaseSubtract(phiRightAbove[component], phiLeftAbove[component]);
  VxBelow = (1./2.)*phaseSubtract(phiRightBelow[component], phiLeftBelow[component]);
  VyLeft = (1./2.)*phaseSubtract(phiLeftAbove[component], phiLeftBelow[component]);
  VyRight = (1./2.)*phaseSubtract(phiRightAbove[component], phiRightBelow[component]);

  VxField[z+y*zSize+x*zSize*ySize] = Vx;
  VyField[z+y*zSize+x*zSize*ySize] = Vy;
  vortField[z+y*zSize+x*zSize*ySize] = (1./2.)*((VyRight - VyLeft)-(VxAbove-VxBelow));

}

__global__ void getPlotDetailsForMF(dcmplx *QField, dcmplx *rhoField, double *phaseField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int component = lattice[6];
  //printf("mf = %d", component);
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];

  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }

  phi[component] = QField[2*component+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize]
            +QField[2*component+1+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];

  /* For Phase */
  phaseField[y+x*ySize] = phase(phi[component]);

  rhoField[y+x*ySize] = 0;
  for(n = 2*component; n <= 2*component + 1; n++){
  rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }

}

__global__ void getTotalRho(dcmplx *QField, dcmplx *rhoField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int component = lattice[6];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];
  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }
  rhoField[y+x*ySize] = 0;
  for(n = 0; n <= vectorSize ; n++){
  rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }
  }

__global__ void getTotalRhoComp(dcmplx *QField, dcmplx *rhoField, int *lattice)
{
  int xSize = lattice[0];
  int ySize = lattice[2];
  int zSize = lattice[4];
  int component = lattice[6];
  int x = blockIdx.x * blockDim.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;
  int z = blockIdx.z * blockDim.z + threadIdx.z;
  int n;
  dcmplx phi[ spinComps ];
  dcmplx psi[ vectorSize ];
  for(n = 0; n < vectorSize; n++){
    psi[n] = QField[n+z*vectorSize+y*vectorSize*zSize+x*zSize*ySize*vectorSize];
  }
  rhoField[y+x*ySize] = 0;
  for(n = 2*component; n <= 2*component + 1 ; n++){
    rhoField[y+x*ySize] += psi[n]*conj(psi[n]);
  }

}

""")

reason: from pycuda


text:
import pandas as pd
import unittest
import numpy as np

from typing import List, Dict, Union, Any

import sys, os

cwd = os.path.abspath(os.curdir)
sys.path.insert(0, cwd + "/src/")

from GeneralProcess import Base
import tempfile, logging
# ---------------------------------------------------------------------------- #

class TestLogging(unittest.TestCase):

    def test_empty_path(self):
        out = Base.createLogger(log_path=None)
        self.assertEqual(None, out)

    def test_instantiation(self):
        with tempfile.TemporaryDirectory() as tmpdirname:
            log_path = os.path.abspath(tmpdirname)

            log = Base.createLogger(log_path=log_path)
            self.assertIsNotNone(log)

            log.info("This is an info line")
            log.debug("This is a debug line")
            log.warning("This is a warning")
            log.error("This is an error")

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

max_stars_repo_name: haganenoneko/EphysAnalysisTools

max_issues_repo_name: haganenoneko/EphysAnalysisTools

max_forks_repo_name: haganenoneko/EphysAnalysisTools


text:
# Copyright 2021 Huawei Technologies Co., Ltd
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================
"""
Test case for fully connected layer module
"""

import pytest
import numpy as np
import mindspore.nn as nn
from mindspore import Tensor
from mindspore_rl.network.fully_connected_net import FullyConnectedLayers

class Net(nn.Cell):
    """
    Ground truth network
    """

    def __init__(self):
        super().__init__()
        self.fc1 = nn.Dense(10, 20)
        self.fc2 = nn.Dense(20, 3)
        self.relu = nn.ReLU()

    def construct(self, x):
        x = self.fc1(x)
        x = self.relu(x)
        x = self.fc2(x)
        return self.relu(x)

@pytest.mark.level0
@pytest.mark.platform_x86_gpu_training
@pytest.mark.platform_x86_ascend_training
@pytest.mark.platform_arm_ascend_training
@pytest.mark.platform_x86_cpu
@pytest.mark.env_onecard
def test_fully_connected_layer():
    '''
    Feature: Test test_fully_connected_layer class
    Description: Test test_fully_connected_layer class
    Expectation: success.
    '''
    net = Net()
    fc_layers = FullyConnectedLayers(fc_layer_params=[10, 20, 3])
    net_fc1_weight = net.trainable_params()[0]
    net_fc2_weight = net.trainable_params()[2]
    for k, v in fc_layers.parameters_dict().items():
        if 'weight' in k and '0' in k:
            v.set_data(net_fc1_weight)
        if 'weight' in k and '1' in k:
            v.set_data(net_fc2_weight)

    input_x = Tensor(np.random.random([5, 10]).astype(np.float32))
    expected = net(input_x)
    actual = fc_layers(input_x)

    assert np.allclose(expected.asnumpy(), actual.asnumpy(), atol=1e-5)

if __name__ == '__main__':
    test_fully_connected_layer()

max_stars_repo_name: mindspore-ai/reinforcement

max_issues_repo_name: mindspore-ai/reinforcement

max_forks_repo_name: mindspore-ai/reinforcement