text:
import numpy as np

def IOU(x, centroids):
    similarities = []
    k = len(centroids)
    for centroid in centroids:
        c_w, c_h = centroid
        w, h = x
        if c_w >= w and c_h >= h:
            similarity = w * h / (c_w * c_h)
        elif c_w >= w and c_h <= h:
            similarity = w * c_h / (w * h + (c_w - w) * c_h)
        elif c_w <= w and c_h >= h:
            similarity = c_w * h / (w * h + c_w * (c_h - h))
        else:  # means both w,h are bigger than c_w and c_h respectively
            similarity = (c_w * c_h) / (w * h)
        similarities.append(similarity)  # will become (k,) shape
    return np.array(similarities)

def avg_IOU(X, centroids):
    n, d = X.shape
    sum = 0.
    for i in range(X.shape[0]):
        # note IOU() will return array which contains IoU for each centroid and X[i] // slightly ineffective, but I am too lazy
        sum += max(IOU(X[i], centroids))
    return sum / n

def write_anchors_to_file(centroids, X, anchor_file):
    f = open(anchor_file, 'w')

    anchors = centroids.copy()
    print(anchors.shape)

    for i in range(anchors.shape[0]):
        anchors[i][0] *= images_width / 32.
        anchors[i][1] *= images_height / 32.

    widths = anchors[:, 0]
    sorted_indices = np.argsort(widths)

    print('Anchors = ', anchors[sorted_indices])

    for i in sorted_indices[:-1]:
        f.write('%0.2f,%0.2f, ' % (anchors[i, 0], anchors[i, 1]))

    # there should not be comma after last anchor, that's why
    f.write('%0.2f,%0.2f\n' % (anchors[sorted_indices[-1:], 0], anchors[sorted_indices[-1:], 1]))

    f.write('%f\n' % (avg_IOU(X, centroids)))
    print()

def kmeans(X, centroids, eps, anchor_file):
    N = X.shape[0]
    iterations = 0
    k, dim = centroids.shape
    prev_assignments = np.ones(N) * (-1)
    iter = 0
    old_D = np.zeros((N, k))

    while True:
        D = []
        iter += 1
        for i in range(N):
            d = 1 - IOU(X[i], centroids)
            D.append(d)
        D = np.array(D)  # D.shape = (N,k)

        print("iter {}: dists = {}".format(iter, np.sum(np.abs(old_D - D))))

        # assign samples to centroids
        assignments = np.argmin(D, axis=1)

        if (assignments == prev_assignments).all():
            print("Centroids = ", centroids)
            write_anchors_to_file(centroids, X, anchor_file)
            return

        # calculate new centroids
        centroid_sums = np.zeros((k, dim), np.float)
        for i in range(N):
            centroid_sums[assignments[i]] += X[i]
        for j in range(k):
            centroids[j] = centroid_sums[j] / (np.sum(assignments == j))

        prev_assignments = assignments.copy()
        old_D = D.copy()


text:
def get_scorecard_objects_and_metrics(scorecard, path_scorecard):
    import os, sys
    import pandas as pd
    import xlrd
    import numpy as np
    from openpyxl import load_workbook
    import xlsxwriter
    import itertools
    import math
    import json

    '''Here we get all lists from all metrics and fields'''
    METRIC_NAMES_BY_SHEET = []
    TYPES_BY_SHEET = []
    UNITS_BY_SHEET = []
    GREATORLESSS_BY_SHEET = []
    POST_RESULTS_BY_SHEET = []
    FLAG_BY_SHEET = []
    OUTOFSPEC_BY_SHEET = []
    DESCRIPT_BY_SHEET = []
    CALC_BY_SHEET = []

    '''Remove all of the fields we dont want'''
    sheets = scorecard.sheet_names
    for sheet in sheets:
        df = pd.read_excel(path_scorecard, sheet_name= sheet)  # get data from sheet
        column_names = list(df.columns)
        metric_col = column_names[1] #metric sheet will be the reference for deleting empty rows or label rows
        length_col = len(df[metric_col]) #original number of rows in the excel sheet

        if sheet == sheets[0]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Quick Assessment' \
                        or df[metric_col][i] == 'System Margins' or (df[metric_col][i]) == 'Precision' \
                        or (df[metric_col][i]) == 'EVRs - Out of Spec' or (df[metric_col][i]) == 'Entry - Out of Spec'\
                        or (df[metric_col][i]) == 'Parachute Descent - Out of Spec' \
                        or (df[metric_col][i]) == 'Powered Flight - Out of Spec' \
                        or (df[metric_col][i]) == 'Heatshield Separation Constraints'\
                        or (df[metric_col][i]) == 'Additional Critical Metrics':
                    df.drop([i], inplace = True)
                    DF_SummaryOutofspec = df

            METRIC_NAMES_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[7]]))
            DESCRIPT_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[8]]))
            CALC_BY_SHEET.append(list(DF_SummaryOutofspec[column_names[9]]))

        if sheet == sheets[1]:
            for i in range(length_col):
                metric_col = column_names[1]
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Quick Assessment' \
                        or df[metric_col][i] == 'Terrain Interaction Failure Rate' or (df[metric_col][i]) == 'Precision' \
                        or (df[metric_col][i]) == 'EVRs - Out of Spec' or (df[metric_col][i]) == 'System Margins'\
                        or (df[metric_col][i]) == 'Precision':
                    df.drop([i], inplace = True)
                    DF_DecisionMeeting = df

            METRIC_NAMES_BY_SHEET.append(list(DF_DecisionMeeting[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_DecisionMeeting[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_DecisionMeeting[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_DecisionMeeting[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_DecisionMeeting[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_DecisionMeeting[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_DecisionMeeting[column_names[7]]))
            DESCRIPT_BY_SHEET.append(list(DF_DecisionMeeting[column_names[8]]))
            CALC_BY_SHEET.append(list(DF_DecisionMeeting[column_names[9]]))

        if sheet == sheets[2]:
            for i in range(length_col):
                algo = df[metric_col][i]
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric':
                    df.drop([i], inplace = True)
                    DF_EDLMETRICS = df

            METRIC_NAMES_BY_SHEET.append(list(DF_EDLMETRICS[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_EDLMETRICS[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_EDLMETRICS[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_EDLMETRICS[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_EDLMETRICS[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_EDLMETRICS[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_EDLMETRICS[column_names[7]]))
            DESCRIPT_BY_SHEET.append(list(DF_EDLMETRICS[column_names[9]]))
            CALC_BY_SHEET.append(list(DF_EDLMETRICS[column_names[10]]))

        if sheet == sheets[3]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Terrain Relative Navigation' \
                        or df[metric_col][i] == 'TRN Error Budget' \
                        or (df[metric_col][i]) == 'Landing Performance':
                    df.drop([i], inplace = True)
                DF_TRN = df

            METRIC_NAMES_BY_SHEET.append(list(DF_TRN[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_TRN[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_TRN[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_TRN[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_TRN[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_TRN[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_TRN[column_names[7]]))
            CALC_BY_SHEET.append(list(DF_TRN[column_names[9]]))
            DESCRIPT_BY_SHEET.append(list(DF_TRN[column_names[8]]))

        if sheet == sheets[4]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Rover Mass Properties' \
                        or (df[metric_col][i]) == 'Precision' or (df[metric_col][i]) == 'PDV Dynamic State (1%-tile)' \
                        or (df[metric_col][i]) == 'PDV Dynamic State (99%-tile)' \
                        or (df[metric_col][i]) == 'Close Clearances (mm)' \
                        or (df[metric_col][i]) == 'Close Clearances (in)' \
                        or (df[metric_col][i]) == 'Margined Loads (N), LUF = 1.2':
                    df.drop([i], inplace=True)
                    DF_BUDMETRICS = df

            calculation7 = list(DF_BUDMETRICS[column_names[8]])
            n = len(calculation7)
            description7 = [0] * n

            METRIC_NAMES_BY_SHEET.append(list(DF_BUDMETRICS[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_BUDMETRICS[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_BUDMETRICS[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_BUDMETRICS[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_BUDMETRICS[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_BUDMETRICS[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_BUDMETRICS[column_names[7]]))
            CALC_BY_SHEET.append(list(DF_BUDMETRICS[column_names[8]]))
            DESCRIPT_BY_SHEET.append(description7)

        if sheet == sheets[5]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Entry Environments' \
                        or (df[metric_col][i]) == 'Entry Guidance' or (df[metric_col][i]) == 'Entry Control':
                    df.drop([i], inplace = True)
                    DF_Entry = df

            metrics3 = list(DF_Entry[column_names[1]])
            types3 = list(DF_Entry[column_names[2]])
            units3 = list(DF_Entry[column_names[3]])
            post_results3 = list(DF_Entry[column_names[4]])
            great_or_less3 = list(DF_Entry[column_names[5]])
            flags3 = list(DF_Entry[column_names[6]])
            out_of_spec3 = list(DF_Entry[column_names[7]])
            description3 = list(DF_Entry[column_names[8]])
            calculation3 = list(DF_Entry[column_names[9]])

            METRIC_NAMES_BY_SHEET.append(metrics3)
            TYPES_BY_SHEET.append(types3)
            UNITS_BY_SHEET.append(units3)
            POST_RESULTS_BY_SHEET.append(post_results3)
            GREATORLESSS_BY_SHEET.append(great_or_less3)
            FLAG_BY_SHEET.append(flags3)
            OUTOFSPEC_BY_SHEET.append(out_of_spec3)
            CALC_BY_SHEET.append(calculation3)
            DESCRIPT_BY_SHEET.append(description3)

        if sheet == sheets[6]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Parachute Deployment Conditions' \
                        or (df[metric_col][i]) == 'Parachute Opening Load Indicators (New)' \
                        or (df[metric_col][i]) == 'Pre-HSS' or df[metric_col][i] == 'Heatshield Separation Conditions'\
                        or df[metric_col][i] == 'Post-HSS' or df[metric_col][i] == 'Backshell Separation Conditions':
                    df.drop([i], inplace = True)
                    DF_ParachuteDescent = df

            metrics4 = list(DF_ParachuteDescent[column_names[1]])
            types4 = list(DF_ParachuteDescent[column_names[2]])
            units4 = list(DF_ParachuteDescent[column_names[3]])
            post_results4 = list(DF_ParachuteDescent[column_names[4]])
            great_or_less4 = list(DF_ParachuteDescent[column_names[5]])
            flags4 = list(DF_ParachuteDescent[column_names[6]])
            out_of_spec4 = list(DF_ParachuteDescent[column_names[7]])
            description4 = list(DF_ParachuteDescent[column_names[8]])
            calculation4 = list(DF_ParachuteDescent[column_names[9]])

            METRIC_NAMES_BY_SHEET.append(metrics4)
            TYPES_BY_SHEET.append(types4)
            UNITS_BY_SHEET.append(units4)
            POST_RESULTS_BY_SHEET.append(post_results4)
            GREATORLESSS_BY_SHEET.append(great_or_less4)
            FLAG_BY_SHEET.append(flags4)
            OUTOFSPEC_BY_SHEET.append(out_of_spec4)
            CALC_BY_SHEET.append(calculation4)
            DESCRIPT_BY_SHEET.append(description4)

        if sheet == sheets[7]:
            for i in range(length_col):
                metric_col = column_names[2]
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace = True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Backshell Sep' \
                        or (df[metric_col][i]) == 'Powered Approach' or (df[metric_col][i]) == 'Constant Velocity' \
                        or df[metric_col][i] == 'Constant Decel' or df[metric_col][i] == 'Throttle Down' \
                        or df[metric_col][i] == 'Rover Separation' or df[metric_col][i] == 'Touchdown' \
                        or df[metric_col][i] == 'Flyaway':
                    df.drop([i], inplace = True)
                    DF_Powereddesent = df

            METRIC_NAMES_BY_SHEET.append(list(DF_Powereddesent[column_names[2]]))
            TYPES_BY_SHEET.append(list(DF_Powereddesent[column_names[3]]))
            UNITS_BY_SHEET.append(list(DF_Powereddesent[column_names[4]]))
            POST_RESULTS_BY_SHEET.append(list(DF_Powereddesent[column_names[5]]))
            GREATORLESSS_BY_SHEET.append(list(DF_Powereddesent[column_names[6]]))
            FLAG_BY_SHEET.append(list(DF_Powereddesent[column_names[7]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_Powereddesent[column_names[8]]))
            CALC_BY_SHEET.append(list(DF_Powereddesent[column_names[10]]))
            DESCRIPT_BY_SHEET.append(list(DF_Powereddesent[column_names[9]]))

        if sheet == sheets[8]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'EVRs - Out of Spec' \
                        or (df[metric_col][i]) == 'Precision' or (df[metric_col][i]) == 'Timeline Margins':
                    df.drop([i], inplace=True)
                    DF_TIMELINE = df

            METRIC_NAMES_BY_SHEET.append(list(DF_TIMELINE[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_TIMELINE[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_TIMELINE[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_TIMELINE[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_TIMELINE[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_TIMELINE[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_TIMELINE[column_names[7]]))
            CALC_BY_SHEET.append(list(DF_TIMELINE[column_names[9]]))
            DESCRIPT_BY_SHEET.append(list(DF_TIMELINE[column_names[8]]))

        if sheet == sheets[9]:
            for i in range(length_col):
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'EDLCOMM MFSK Tones (R9.4)':
                    df.drop([i], inplace=True)
                    DF_COMM = df

            METRIC_NAMES_BY_SHEET.append(list(DF_COMM[column_names[1]]))
            TYPES_BY_SHEET.append(list(DF_COMM[column_names[2]]))
            UNITS_BY_SHEET.append(list(DF_COMM[column_names[3]]))
            POST_RESULTS_BY_SHEET.append(list(DF_COMM[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list(DF_COMM[column_names[5]]))
            FLAG_BY_SHEET.append(list(DF_COMM[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_COMM[column_names[7]]))
            CALC_BY_SHEET.append(list(DF_COMM[column_names[9]]))
            DESCRIPT_BY_SHEET.append(list(DF_COMM[column_names[8]]))

        if sheet == sheets[10]:
            for i in range(length_col):
                metric_col = column_names[2]
                if type(df[metric_col][i]) == float and df[column_names[3]][i] != str('50%-tile') \
                        and df[column_names[3]][i] != str('99%-tile'):
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Received Power' \
                        or df[metric_col][i] == 'Received Power Watermarks' or df[metric_col][i] == 'Orbiter Tracking to MSL' \
                        or df[metric_col][i] == 'Event Times Relative to Nominal Entry Interface Time (time = 540 s)':
                    df.drop([i], inplace=True)
                DF_TELECOM = df

            METRIC_NAMES_BY_SHEET.append(list(DF_TELECOM[column_names[2]]))
            TYPES_BY_SHEET.append(list(DF_TELECOM[column_names[3]]))
            UNITS_BY_SHEET.append(list(DF_TELECOM[column_names[4]]))
            POST_RESULTS_BY_SHEET.append(list(DF_TELECOM[column_names[5]]))
            GREATORLESSS_BY_SHEET.append(list(DF_TELECOM[column_names[6]]))
            FLAG_BY_SHEET.append(list(DF_TELECOM[column_names[7]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_TELECOM[column_names[8]]))
            CALC_BY_SHEET.append(list(DF_TELECOM[column_names[10]]))
            DESCRIPT_BY_SHEET.append(list(DF_TELECOM[column_names[9]]))

        if sheet == sheets[11]:
            for i in range(length_col):
                metric_col = column_names[2]
                if type(df[metric_col][i]) == float and df[column_names[5]][i] != str('Scalar'):
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'ID':
                    df.drop([i], inplace=True)
                DF_EFC = df

            METRIC_NAMES_BY_SHEET.append(list(DF_EFC[column_names[4]]))
            TYPES_BY_SHEET.append(list(DF_EFC[column_names[5]]))
            UNITS_BY_SHEET.append(list(DF_EFC[column_names[6]]))
            POST_RESULTS_BY_SHEET.append(list(DF_EFC[column_names[7]]))
            GREATORLESSS_BY_SHEET.append(list(DF_EFC[column_names[8]]))
            FLAG_BY_SHEET.append(list(DF_EFC[column_names[9]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_EFC[column_names[10]]))
            CALC_BY_SHEET.append(list(DF_EFC[column_names[12]]))
            DESCRIPT_BY_SHEET.append(list(DF_EFC[column_names[11]]))

        if sheet == sheets[12]:
            for i in range(length_col):
                metric_col = column_names[3]
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Monte Carlo' \
                        or df[metric_col][i] == 'Landing Accuracy' or df[metric_col][i] == 'Touchdown Conditions' \
                        or df[metric_col][i] == 'Propellant Consumption' or df[metric_col][i] == 'Timeline Margin' \
                        or df[metric_col][i] == 'Aeroheating' or df[metric_col][i] == 'Loads' \
                        or df[metric_col][i] == 'Heatshield Separation Constraints':
                    df.drop([i], inplace=True)
                DF_Summary = df

            METRIC_NAMES_BY_SHEET.append(list(DF_Summary[column_names[3]]))
            TYPES_BY_SHEET.append(list(DF_Summary[column_names[8]]))
            UNITS_BY_SHEET.append(list(DF_Summary[column_names[6]]))
            POST_RESULTS_BY_SHEET.append(list(DF_Summary[column_names[7]]))
            GREATORLESSS_BY_SHEET.append(list(DF_Summary[column_names[4]]))
            n = DF_Summary.shape[0]
            FLAG_BY_SHEET.append(list([str('nan')]*n))
            OUTOFSPEC_BY_SHEET.append(list(DF_Summary[column_names[5]]))
            CALC_BY_SHEET.append(list(DF_Summary[column_names[9]]))
            DESCRIPT_BY_SHEET.append(list([str('nan')]*n))

        if sheet == sheets[13]:
            for i in range(length_col):
                metric_col = column_names[3]
                if type(df[metric_col][i]) == float and df[column_names[5]][i] != str('Scalar'):
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Failure Mechanism':
                    df.drop([i], inplace=True)
                DF_Failure = df

            METRIC_NAMES_BY_SHEET.append(list(DF_Failure[column_names[3]]))
            TYPES_BY_SHEET.append(list(DF_Failure[column_names[5]]))
            n = DF_Failure.shape[0]
            UNITS_BY_SHEET.append(list([str('nan')]*n))
            POST_RESULTS_BY_SHEET.append(list(DF_Failure[column_names[4]]))
            GREATORLESSS_BY_SHEET.append(list([str('nan')]*n))
            FLAG_BY_SHEET.append(list([str('nan')]*n))
            OUTOFSPEC_BY_SHEET.append(list([str('nan')]*n))
            CALC_BY_SHEET.append(list(DF_Failure[column_names[6]]))
            DESCRIPT_BY_SHEET.append(list(DF_Failure[column_names[7]]))

        if sheet == sheets[14]:
            for i in range(length_col):
                metric_col = column_names[3]
                if type(df[metric_col][i]) == float:
                    df.drop([i], inplace=True)
                elif (df[metric_col][i]) == 'Metric' or df[metric_col][i] == 'Landing Accuracy' \
                        or df[metric_col][i] == 'Touchdown Conditions' or df[metric_col][i] == 'Propellant Consumption' \
                        or df[metric_col][i] == 'Timeline' or df[metric_col][i] == 'EDL Comm ' \
                        or df[metric_col][i] == 'Aeroheating' or df[metric_col][i] == 'Entry Loads (Does not include PD)' \
                        or df[metric_col][i] == 'Parachute Deploy Constraints' \
                        or df[metric_col][i] == 'Heatshield Separation Constraints' \
                        or df[metric_col][i] == 'Radar Constraints' or df[metric_col][i] == 'Rover Wind Loading' \
                        or df[metric_col][i] == 'Long Term Recontact' or df[metric_col][i] == 'Downrange Distance for MARDI' \
                        or df[metric_col][i] == 'Backshell Separation Conditions' or df[metric_col][i] == 'Accordion Values' \
                        or df[metric_col][i] == 'Rover Separation Conditions' or df[metric_col][i] == 'Sky Crane Dynamics' \
                        or df[metric_col][i] == 'Flyaway Start State ':
                    df.drop([i], inplace=True)
                DF_Detailed = df

            METRIC_NAMES_BY_SHEET.append(list(DF_Detailed[column_names[3]]))
            TYPES_BY_SHEET.append(list(DF_Detailed[column_names[9]]))
            UNITS_BY_SHEET.append(list(DF_Detailed[column_names[7]]))
            POST_RESULTS_BY_SHEET.append(list(DF_Detailed[column_names[8]]))
            GREATORLESSS_BY_SHEET.append(list(DF_Detailed[column_names[4]]))
            FLAG_BY_SHEET.append(list(DF_Detailed[column_names[6]]))
            OUTOFSPEC_BY_SHEET.append(list(DF_Detailed[column_names[5]]))
            CALC_BY_SHEET.append(list(DF_Detailed[column_names[10]]))
            n = DF_Detailed.shape[0]
            DESCRIPT_BY_SHEET.append(list([str('nan')]*n))

    '''Create objects and classes'''
    class ScoreCardCategory(object):
        def __init__(self, name = 0,  entries = 0):
            self.name = name# these are the tab names (categories)
            self.entries= entries# These are all the metrics in each
        def __call__(self, name, entries):
            print(name)
            print(entries)

    # A = edl_metric.name #How we call class objects
    # B = edl_metric.entries

    class ScoreCardMetrics(object):
        def __init__(self, name = 0, type = 0 , units=0, POST_results = 0, GreatOrLess = 0, Flag = 0, OutOfSpec = 0, Enum = 0,
                     Description = 0, Calculation = 0, CalculationCheckedBy = 0, MetricCheckedBy = 0, FlagOutOfSpecOwner = 0):
            self.name = name
            self.type = type
            self.units = units
            self.POST_results = POST_results
            self.GreatOrLess = GreatOrLess
            self.Flag = Flag
            self.OutOfSpec = OutOfSpec
            self.Enum = Enum
            self.Description = Description
            self.Calculation = Calculation
            self.CalculationCheckedBy = CalculationCheckedBy
            self.MetricCheckedBy = MetricCheckedBy
            self.FlagOutOfSpecOwner = FlagOutOfSpecOwner
        def __call__(self, name, type, units, POST_results, GreatOrLess, Flag, OutOfSpec, Enum,
                     Description, Calculation, CalculationCheckedBy, MetricCheckedBy, FlagOutOfSpecOwner):
            print(name)
            print(type)
            print(units)
            print(POST_results)
            print(GreatOrLess)
            print(Flag)
            print(OutOfSpec)
            print(Enum)
            print(Description)
            print(Calculation)
            print(CalculationCheckedBy)
            print(MetricCheckedBy)
            print(FlagOutOfSpecOwner)

    METRIC_OBJECTS = []
    METRIC_OBJECTS_GROUPED = []

    for j in range(len(METRIC_NAMES_BY_SHEET)): # number of columns (number of excel sheets)
        METRIC_OBJECTS = []
        for i in range(len(METRIC_NAMES_BY_SHEET[j])): # number of metrics in each sheet
            METRIC_OBJECTS.append(ScoreCardMetrics(name = METRIC_NAMES_BY_SHEET[j][i], type=TYPES_BY_SHEET[j][i],
                                                     units=UNITS_BY_SHEET[j][i], POST_results=POST_RESULTS_BY_SHEET[j][i],
                                                     GreatOrLess=GREATORLESSS_BY_SHEET[j][i], Flag=FLAG_BY_SHEET[j][i],
                                                     OutOfSpec=OUTOFSPEC_BY_SHEET[j][i], Enum=0,
                                                     Description=DESCRIPT_BY_SHEET[j][i],
                                                     Calculation=CALC_BY_SHEET[j][i], CalculationCheckedBy=0,
                                                     MetricCheckedBy=0, FlagOutOfSpecOwner=0))
            METRIC_OBJECTS = list(METRIC_OBJECTS)
        METRIC_OBJECTS_GROUPED.append(METRIC_OBJECTS)

    '''Make single list of all of the metric objects'''
    METRIC_OBJECTS = [item for items in METRIC_OBJECTS_GROUPED for item in items]
    '''Make a list of all of the POST results'''
    print(METRIC_OBJECTS[:])

    summaryOutofspec_objects = METRIC_OBJECTS_GROUPED[0]
    decisionmeeting_objects = METRIC_OBJECTS_GROUPED[1]
    edl_metric_objects = METRIC_OBJECTS_GROUPED[2]
    trn_objects = METRIC_OBJECTS_GROUPED[3]
    bud_metrics_objects = METRIC_OBJECTS_GROUPED[4]
    entry_objects = METRIC_OBJECTS_GROUPED[5]
    parachute_descent_objects = METRIC_OBJECTS_GROUPED[6]
    powered_flight_objects = METRIC_OBJECTS_GROUPED[7]
    timeline_objects = METRIC_OBJECTS_GROUPED[8]
    comm_objects = METRIC_OBJECTS_GROUPED[9]
    telecom_objects = METRIC_OBJECTS_GROUPED[10]
    efc_objects = METRIC_OBJECTS_GROUPED[11]
    summary_objects = METRIC_OBJECTS_GROUPED[12]
    failures_objects = METRIC_OBJECTS_GROUPED[13]
    detailed_objects = METRIC_OBJECTS_GROUPED[14]

    '''Here we have the list of metrics for the entry, parachute descent, powered descent, TRN and BUD metrics from the scorecard'''
    list_of_metrics = [item for items in METRIC_NAMES_BY_SHEET for item in items]
    '''Save the list of metrics to excel'''

    '''
    import csv
    df = pd.DataFrame(list_of_metrics)
    listofmetrics = pd.ExcelWriter("list_of_metrics_complete.xlsx", engine='xlsxwriter')
    df.to_excel(listofmetrics,'Sheet 1')
    listofmetrics.save()
    '''

    '''After modifications, get the new list of metrics '''

    df = pd.read_excel('/Users/ssantini/Code/ExtractDataMatlab/ExtractScorecardData/list_of_metrics_complete_ver2.xlsx')  # get data from sheet
    col_name = list(df.columns)
    list_of_metrics = list(df[0])

    summary_outofspec = ScoreCardCategory(scorecard.sheet_names[0],  summaryOutofspec_objects)
    decision_meeting = ScoreCardCategory(scorecard.sheet_names[1],  decisionmeeting_objects)
    edl_metric = ScoreCardCategory(scorecard.sheet_names[2], edl_metric_objects)
    trn= ScoreCardCategory(scorecard.sheet_names[3], trn_objects)
    bud_metrics = ScoreCardCategory(scorecard.sheet_names[4], bud_metrics_objects)
    entry = ScoreCardCategory(scorecard.sheet_names[5], entry_objects)
    parachute_descent= ScoreCardCategory(scorecard.sheet_names[6], parachute_descent_objects)
    powered_flight = ScoreCardCategory(scorecard.sheet_names[7], powered_flight_objects)
    edltimeline = ScoreCardCategory(scorecard.sheet_names[8], timeline_objects)
    comms = ScoreCardCategory(scorecard.sheet_names[9], comm_objects)
    telecom = ScoreCardCategory(scorecard.sheet_names[10], telecom_objects)
    efc = ScoreCardCategory(scorecard.sheet_names[11], efc_objects)
    summary = ScoreCardCategory(scorecard.sheet_names[12], summary_objects)
    failure = ScoreCardCategory(scorecard.sheet_names[13], failures_objects)
    detailed = ScoreCardCategory(scorecard.sheet_names[14], detailed_objects)

    return (summary_outofspec, decision_meeting, edl_metric, trn, bud_metrics, entry, parachute_descent, powered_flight,
        edltimeline, comms, telecom, efc, summary, failure, detailed, list_of_metrics)
#
#
# [summary_outofspec, decision_meeting, edl_metric, trn, bud_metrics, entry, parachute_descent, powered_flight,
#         edltimeline, comms, telecom, efc, summary, failure, detailed, list_of_metrics] = get_scorecard_objects_and_metrics(scorecard, path_scorecard)
#
# print(get_scorecard_objects_and_metrics.get_scorecard_objects_and_metrics())
# i = 1

ext: py

reason: import numpy


text:
#!/usr/bin/env python3
#
# hausdorff_subsampling.py: investigates the behaviour of the Hausdorff
# distance under random subsampling.

import argparse
import sys

import numpy as np
import matplotlib.pyplot as plt

def diameter(X):
    '''
    Calculates the diameter of a given set of points and returns it.
    '''

    # Calculate pairwise distance matrix of all points; this is required
    # as a prerequisite of the diameter calculation.
    D = np.sum((X[None, :] - X[:, None])**2, -1)**0.5

    return np.amax(D)

def pairwise_distances(X):
    '''
    Calculates pairwise distance matrix of a given finite metric space and
    returns it in the form of a vector.
    '''

    D = np.sum((X[None, :] - X[:, None])**2, -1)**0.5
    return D.ravel()

def hausdorff_distance(X, Y):
    '''
    Calculates the Hausdorff distance between two finite metric spaces,
    i.e. two finite point clouds.
    '''

    d_x_Y = 0

    for x in X:
        d_x_y = sys.float_info.max
        for y in Y:
            d = np.linalg.norm(x - y, ord=None)
            d_x_y = min(d_x_y, d)

        d_x_Y = max(d_x_Y, d_x_y)

    d_y_X = 0

    for y in Y:
        d_y_x = sys.float_info.max
        for x in X:
            d = np.linalg.norm(x - y, ord=None)
            d_y_x = min(d_y_x, d)

        d_y_X = max(d_y_X, d_y_x)

    return max(d_x_Y, d_y_X)

if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument('-d', '--dimension', type=int, nargs='+', default=2)
    parser.add_argument('-n', '--n_points', type=int, default=100)
    parser.add_argument('-s', '--n_subsamples', type=int, default=20)

    args = parser.parse_args()

    n_points = args.n_points
    d = args.dimension
    n_subsamples = args.n_subsamples

    for d in args.dimension:

        print(f'# dim = {d}')

        np.random.seed(42)

        X = np.random.normal(size=(n_points, d))
        X_diam = diameter(X)

        X_pairwise_distances = pairwise_distances(X)

        X_mean_distance = np.mean(X_pairwise_distances)
        X_mean_distance_upper_bound = np.sqrt(2 * d)

        for m in range(1, n_points + 1):

            hausdorff_distances = []
            Y_diameters = []

            for _ in range(n_subsamples):
                Y = X[np.random.choice(X.shape[0], m, replace=False)]

                Y_diameters.append(diameter(Y))
                hausdorff_distances.append(hausdorff_distance(X, Y))

            print(
                m,
                np.mean(hausdorff_distances),
                np.std(hausdorff_distances)
            )

            sys.stdout.flush()

        print()

max_stars_repo_name: duxiaodan/topological-autoencoders

max_issues_repo_name: duxiaodan/topological-autoencoders

max_forks_repo_name: duxiaodan/topological-autoencoders


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 19 12:32:11 2019

@author: gpang
"""

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
import time
#from SALib.sample import sobol_sequence
import scipy as sci
import scipy.io as sio
#

class one_GP:

    def __init__(self):
        pass

    def model(self, dt, a, b, c, Nu, Nf, kernel_type = 'SE', is_noise = True):
        self.dt = dt
        self.a = a
        self.b = b
        self.c = c
        self.Nu = Nu
        self.Nf = Nf
        self.kernel_type = kernel_type

        self.is_noise = is_noise

    def u_exact(self, x,t, u_exa, t_exa, x_exa, dim):
        if dim == 1:
            tt = np.ndarray.flatten(t_exa)
            uu1 =  np.ndarray.flatten(u_exa[0,:])
            uu2 = np.ndarray.flatten(u_exa[-1,:])
            f1 = sci.interpolate.interp1d(tt,uu1,kind='cubic')
            f2 = sci.interpolate.interp1d(tt,uu2,kind='cubic')
            u1 = f1(t)
            u2 = f2(t)
            return np.array([[u1],[u2]],dtype=np.float64)
        elif dim == 2:
            t = t*np.ones((x.shape[0],1),dtype=np.float64)
            [tt, xx] = np.meshgrid(t_exa,x_exa)
            ttt = tt.reshape((-1,1))
            xxx = xx.reshape((-1,1))
            uuu = u_exa.reshape((-1,1))
            return sci.interpolate.griddata(np.concatenate((ttt,xxx),axis=1),uuu, np.concatenate((t,x),axis=1), fill_value = 0.0, method='cubic')

    def f_exact(self, x,t):

        return np.zeros((x.shape[0],1),dtype=np.float64)

    def kernel(self, X, Y, NX, NY, t1, equal=False, diag=False):
        if self.kernel_type == 'SE':
            Y = tf.reshape(Y,(1,-1))
            if diag == False:
                return tf.exp(-0.5* (X-Y)**2/t1**2)
            else:
                return tf.ones((NX,1),dtype=tf.float64)
        elif self.kernel_type == 'Matern1':
            dist = tf.sqrt(self.square_dist(X,Y,t1,equal))
            return (1.0+3.0**0.5*dist)*tf.exp(-3.0**0.5*dist)
        elif self.kernel_type == 'Matern2':
            dist = tf.sqrt(self.square_dist(X,Y,t1,equal))
            return (1.0+5.0**0.5*dist+5.0/3.0*dist**2)*tf.exp(-5.0**0.5*dist)

    def kx(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag == False:
            return (Y-X)/t1**2*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return tf.zeros((NX,1),dtype=tf.float64)

    def ky(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag == False:
            return (X-Y)/t1**2*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return tf.zeros((NX,1),dtype=tf.float64)

    def kxx(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag==False:
            return (-1.0/t1**2+(X-Y)**2/t1**4)*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return -1.0/t1**2 * tf.ones((NX,1),dtype=tf.float64)

    def kyy(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag==False:
            return (-1.0/t1**2+(X-Y)**2/t1**4)*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return -1.0/t1**2 * tf.ones((NX,1),dtype=tf.float64)

    def kxy(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag==False:
            return (1.0/t1**2-(X-Y)**2/t1**4)*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return 1.0/t1**2*tf.ones((NX,1),dtype=tf.float64)

    def kyyx(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag==False:
            return (3*(X-Y)/t1**4-(X-Y)**3/t1**6)*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return tf.zeros((NX,1),dtype=tf.float64)

    def kyxx(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag==False:
            return (3*(Y-X)/t1**4+(X-Y)**3/t1**6)*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return tf.zeros((NX,1),dtype=tf.float64)

    def kxxyy(self, X, Y, NX, NY, t1, diag=False):
        Y = tf.reshape(Y,(1,-1))
        if diag==False:
            return (3.0/t1**4-6*(X-Y)**2/t1**6+(X-Y)**4/t1**8)*tf.exp(-0.5*(X-Y)**2/t1**2)
        else:
            return 3.0/t1**4*tf.ones((NX,1),dtype=tf.float64)

    def Lap2_kernel(self, X, Y, NX, NY, t1, lambda1, lambda2, un_x, un_y, equal=False, diag=False):
        unx = tf.reshape(un_x,[-1])
        uny = tf.reshape(un_y,[-1])
        unx = tf.diag(unx)
        uny = tf.diag(uny)
        if self.kernel_type == 'SE':
            if diag == False:
                k = lambda1**2*tf.matmul(tf.matmul(unx,self.kxy(X,Y,NX, NY, t1,diag)),uny)-lambda1*lambda2*tf.matmul(unx,self.kyyx(X,Y,NX, NY, t1,diag))\
                  -lambda1*lambda2*tf.matmul(self.kyxx(X,Y,NX, NY, t1,diag),uny)+lambda2**2*self.kxxyy(X,Y,NX, NY, t1,diag)
            else:
                k  = lambda1**2* un_x**2*self.kxy(X,Y,NX, NY, t1,diag)-lambda1*lambda2*un_x*self.kyyx(X,Y,NX, NY, t1,diag)\
                  -lambda1*lambda2*un_y*self.kyxx(X,Y,NX, NY, t1,diag)+lambda2**2*self.kxxyy(X,Y,NX, NY, t1,diag)
            return k

    def Lap1_kernel(self, X, Y, NX, NY, t1, lambda1, lambda2, un_x, un_y, equal=False, diag=False): ## -\Delta rather than \Delta

      if self.kernel_type == 'SE':
            unx = tf.reshape(un_x,[-1])
            uny = tf.reshape(un_y,[-1])
            unx = tf.diag(unx)
            uny = tf.diag(uny)
            if diag == False:
                 k = lambda1*tf.matmul(unx,self.kx(X,Y, NX, NY, t1,diag))-lambda2*self.kxx(X,Y,NX, NY, t1,diag)

            else:
                 k = lambda1*un_x*self.kx(X,Y,NX, NY, t1,diag)-lambda2*self.kxx(X,Y,NX, NY, t1,diag)
            return k

    def Lap1_kernel_prime(self, X, Y, NX, NY, t1, lambda1, lambda2, un_x, un_y, equal=False, diag=False): ## -\Delta rather than \Delta

        if self.kernel_type == 'SE':
            unx = tf.reshape(un_x,[-1])
            uny = tf.reshape(un_y,[-1])
            unx = tf.diag(unx)
            uny = tf.diag(uny)
            if diag == False:
                  k = lambda1*tf.matmul(self.ky(X,Y,NX, NY, t1,diag),uny)-lambda2*self.kyy(X,Y,NX, NY, t1,diag)

            else:
                  k = lambda1*un_y*self.ky(X,Y,NX, NY, t1,diag)-lambda2*self.kyy(X,Y,NX, NY, t1,diag)

            return k

    def kernel_uf_train(self, Xu, Xf, Nu, Nf, t1, t3, t5, a, b, c, lambda1, lambda2, un_u, un_f, dt, diag=False):

        if self.kernel_type == 'SE':
            if diag == False:
                ku3u3 = self.kernel(Xu, Xu, Nu, Nu, t1, equal=True)
                ku2u2 = self.kernel(Xu, Xu, Nu, Nu, t3, equal=True)
                ku1u1 = self.kernel(Xu, Xu, Nu, Nu, t5, equal=True)
                kf3f3 = self.kernel(Xf, Xf, Nf, Nf, t1, equal=True)  \
                        + dt**2*b[0]**2*self.Lap2_kernel(Xf, Xf, Nf, Nf, t5, lambda1, lambda2, un_f, un_f, equal=True) \
                        + dt**2*b[1]**2*self.Lap2_kernel(Xf, Xf, Nf, Nf, t3, lambda1, lambda2, un_f, un_f, equal=True)
                kf2f2 = self.kernel(Xf, Xf, Nf, Nf, t3, equal=True) \
                        + dt*a[1,1]*self.Lap1_kernel(Xf, Xf, Nf, Nf, t3, lambda1, lambda2, un_f, un_f, equal=True) \
                        + dt*a[1,1]*self.Lap1_kernel_prime(Xf, Xf, Nf, Nf, t3, lambda1, lambda2, un_f, un_f, equal=True)\
                        +dt**2*a[1,0]**2*self.Lap2_kernel(Xf, Xf, Nf, Nf, t5, lambda1, lambda2, un_f, un_f, equal=True) \
                        +dt**2*a[1,1]**2*self.Lap2_kernel(Xf, Xf, Nf, Nf, t3, lambda1, lambda2, un_f, un_f, equal=True)
                kf1f1 = self.kernel(Xf, Xf, Nf, Nf, t5, equal=True) \
                        +dt*a[0,0]*self.Lap1_kernel(Xf, Xf, Nf, Nf, t5, lambda1, lambda2, un_f, un_f, equal=True)\
                        +dt*a[0,0]*self.Lap1_kernel_prime(Xf, Xf, Nf, Nf, t5, lambda1, lambda2, un_f, un_f, equal=True)\
                        +dt**2*a[0,1]**2*self.Lap2_kernel(Xf, Xf, Nf, Nf, t3, lambda1, lambda2, un_f, un_f, equal=True) \
                        +dt**2*a[0,0]**2*self.Lap2_kernel(Xf, Xf, Nf, Nf, t5, lambda1, lambda2, un_f, un_f, equal=True)

                kf3u3 = self.kernel(Xf, Xu, Nf, Nu, t1)
                kf3u2 = dt*b[1]*self.Lap1_kernel(Xf, Xu, Nf, Nu, t3, lambda1, lambda2, un_f, un_u)
                kf2u2 = self.kernel(Xf, Xu, Nf, Nu, t3) + dt*a[1,1]*self.Lap1_kernel(Xf,Xu,Nf, Nu, t3,lambda1, lambda2, un_f, un_u)
                kf1u2 = dt*a[0,1]*self.Lap1_kernel(Xf, Xu, Nf, Nu, t3, lambda1, lambda2, un_f, un_u)
                kf3u1 = dt*b[0]*self.Lap1_kernel(Xf, Xu,Nf,Nu, t5, lambda1, lambda2, un_f, un_u)
                kf2u1 = dt*a[1,0]*self.Lap1_kernel(Xf, Xu, Nf,Nu,t5, lambda1, lambda2, un_f, un_u)
                kf1u1 = self.kernel(Xf, Xu, Nf,Nu,t5) + dt*a[0,0]*self.Lap1_kernel(Xf, Xu, Nf,Nu,t5, lambda1, lambda2, un_f, un_u)
                kf2f3 = dt*b[1]*self.Lap1_kernel_prime(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f) \
                       +dt**2*b[0]*a[1,0]*self.Lap2_kernel(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f) \
                       +dt**2*b[1]*a[1,1]*self.Lap2_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f)
                kf1f3 = dt*b[0]*self.Lap1_kernel_prime(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f) \
                        + dt**2*b[0]*a[0,0]*self.Lap2_kernel(Xf, Xf,Nf,Nf, t5, lambda1, lambda2, un_f, un_f) \
                        + dt**2*b[1]*a[0,1]*self.Lap2_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f)
                kf1f2 = dt*a[0,1]*self.Lap1_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f) \
                        +dt*a[1,0]*self.Lap1_kernel_prime(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f) \
                        + dt**2*a[1,0]*a[0,0]*self.Lap2_kernel(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f)\
                        + dt**2*a[1,1]*a[0,1]*self.Lap2_kernel(Xf, Xf, Nf, Nf, t3, lambda1, lambda2, un_f, un_f)

                zuu = tf.zeros((Nu,Nu),dtype=tf.float64)
                zuf = tf.zeros((Nu,Nf),dtype=tf.float64)
                zfu = tf.zeros((Nf,Nu),dtype=tf.float64)

                k1 = tf.concat( (ku3u3, zuu, zuu, tf.transpose(kf3u3), zuf, zuf),axis=1)
                k2 = tf.concat( (zuu, ku2u2, zuu, tf.transpose(kf3u2), tf.transpose(kf2u2), tf.transpose(kf1u2)),axis=1)
                k3 = tf.concat( (zuu, zuu, ku1u1, tf.transpose(kf3u1), tf.transpose(kf2u1), tf.transpose(kf1u1)),axis=1)
                k4 = tf.concat( (kf3u3, kf3u2, kf3u1, kf3f3, tf.transpose(kf2f3), tf.transpose(kf1f3)),axis=1)
                k5 = tf.concat( (zfu, kf2u2, kf2u1, kf2f3, kf2f2, tf.transpose(kf1f2)),axis=1)
                k6 = tf.concat( (zfu, kf1u2, kf1u1, kf1f3, kf1f2, kf1f1),axis=1)

                k = tf.concat((k1,k2,k3,k4,k5,k6),axis=0)

                return k
            else:
                ku3u3 = self.kernel(Xu, Xu, Nu,Nu,t1, diag=True)
                ku2u2 = self.kernel(Xu, Xu, Nu,Nu,t3, diag=True)
                ku1u1 = self.kernel(Xu, Xu, Nu,Nu,t5, diag=True)
                kf3f3 = self.kernel(Xf, Xf, Nf,Nf,t1, diag=True)  \
                        + dt**2*b[0]**2*self.Lap2_kernel(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f, diag=True) \
                        + dt**2*b[1]**2*self.Lap2_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f, diag=True)
                kf2f2 = self.kernel(Xf, Xf, Nf,Nf,t3,  diag=True) \
                        + 2.0*dt*a[1,1]*self.Lap1_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f, diag=True) \
                        +dt**2*a[1,0]**2*self.Lap2_kernel(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f, diag=True) \
                        +dt**2*a[1,1]**2*self.Lap2_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f, diag=True)
                kf1f1 = self.kernel(Xf, Xf, Nf,Nf,t5, diag=True) \
                        +2.0*dt*a[0,0]*self.Lap1_kernel(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f, diag=True)\
                        +dt**2*a[0,1]**2*self.Lap2_kernel(Xf, Xf, Nf,Nf,t3, lambda1, lambda2, un_f, un_f, diag=True) \
                        +dt**2*a[0,0]**2*self.Lap2_kernel(Xf, Xf, Nf,Nf,t5, lambda1, lambda2, un_f, un_f, diag=True)

                return tf.concat((ku3u3,ku2u2,ku1u1,kf3f3, kf2f2, kf1f1),axis=0)

    def kernel_u_test(self,  Xt, Xu, Xf, Nt, Nu,Nf,t1, t3, t5, a, b, c, lambda1, lambda2, un_u, un_f, un_t, dt):

        if self.kernel_type == 'SE':
            ku3u3 = self.kernel(Xt, Xu, Nt,Nu,t1)
            ku2u2 = self.kernel(Xt, Xu, Nt,Nu,t3)
            ku1u1 = self.kernel(Xt, Xu, Nt,Nu,t5)

            ku3f3 = self.kernel(Xt, Xf, Nt,Nf,t1)

            ku2f3 = dt*b[1]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t3,lambda1, lambda2, un_t, un_f )
            ku2f2 = self.kernel(Xt, Xf, Nt,Nf,t3) + dt*a[1,1]*self.Lap1_kernel_prime(Xt,Xf,Nt,Nf,t3,lambda1, lambda2, un_t, un_f)
            ku2f1 = dt*a[0,1]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)
            ku1f3 = dt*b[0]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)
            ku1f2 = dt*a[1,0]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)
            ku1f1 = self.kernel(Xt, Xf, Nt,Nf, t5) + dt*a[0,0]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)

            zuu = tf.zeros((Nt,Nu),dtype=tf.float64)
            zuf = tf.zeros((Nt,Nf),dtype=tf.float64)
#            zfu = tf.zeros((Xf.shape[0],Xt.shape[0]),dtype=tf.float64)

            k1 = tf.concat( (ku3u3, zuu, zuu, ku3f3, zuf, zuf),axis=1)
            k2 = tf.concat( (zuu, ku2u2, zuu, ku2f3, ku2f2, ku2f1),axis=1)
            k3 = tf.concat( (zuu, zuu, ku1u1, ku1f3, ku1f2, ku1f1),axis=1)

            k = tf.concat((k1,k2,k3),axis=0)
            return k

    def kernel_f_test(self, Xt, Xu, Xf, Nt,Nu,Nf,t1, t3, t5, a, b, c, lambda1, lambda2, un_u, un_f, un_t, dt):

#        sess1 = tf.Session()
#        sess1.run(tf.global_variables_initializer())
#

        if self.kernel_type == 'SE':
            kf3f3 = self.kernel(Xt, Xf, Nt,Nf,t1)  \
                    + dt**2*b[0]**2*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                    + dt**2*b[1]**2*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)

            kf2f2 = self.kernel(Xt, Xf, Nt,Nf,t3) \
                     + dt*a[1,1]*self.Lap1_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f) \
                    + dt*a[1,1]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf, t3, lambda1, lambda2, un_t, un_f)\
                    +dt**2*a[1,0]**2*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                    +dt**2*a[1,1]**2*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)
            kf1f1 = self.kernel(Xt, Xf, Nt, Nf, t5) \
                    +dt*a[0,0]*self.Lap1_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)\
                    +dt*a[0,0]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)\
                    +dt**2*a[0,1]**2*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f) \
                    +dt**2*a[0,0]**2*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)

            kf3u3 = self.kernel(Xt, Xu, Nt,Nu,t1)
            kf3u2 = dt*b[1]*self.Lap1_kernel(Xt, Xu, Nt,Nu,t3, lambda1, lambda2, un_t, un_u)
            kf2u2 = self.kernel(Xt, Xu, Nt,Nu, t3) + dt*a[1,1]*self.Lap1_kernel(Xt,Xu,Nt,Nu,t3,lambda1, lambda2, un_t, un_u)
            kf1u2 = dt*a[0,1]*self.Lap1_kernel(Xt, Xu, Nt,Nu, t3, lambda1, lambda2, un_t, un_u)
            kf3u1 = dt*b[0]*self.Lap1_kernel(Xt, Xu, Nt,Nu,t5, lambda1, lambda2, un_t, un_u)
            kf2u1 = dt*a[1,0]*self.Lap1_kernel(Xt, Xu, Nt,Nu,t5, lambda1, lambda2, un_t, un_u)
            kf1u1 = self.kernel(Xt, Xu, Nt, Nu, t5) + dt*a[0,0]*self.Lap1_kernel(Xt, Xu, Nt, Nu, t5, lambda1, lambda2, un_t, un_u)

            kf2f3 = dt*b[1]*self.Lap1_kernel_prime(Xt, Xf, Nt, Nf, t3,lambda1, lambda2, un_t, un_f) \
                   +dt**2*b[0]*a[1,0]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                   +dt**2*b[1]*a[1,1]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)
            kf1f3 = dt*b[0]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                    + dt**2*b[0]*a[0,0]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                    + dt**2*b[1]*a[0,1]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)
            kf1f2 = dt*a[0,1]*self.Lap1_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f) \
                        +dt*a[1,0]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                        + dt**2*a[1,0]*a[0,0]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)\
                        + dt**2*a[1,1]*a[0,1]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)

            kf3f2 = dt*b[1]*self.Lap1_kernel(Xt, Xf, Nt,Nf,t3,lambda1, lambda2, un_t, un_f) \
                   +dt**2*b[0]*a[1,0]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                   +dt**2*b[1]*a[1,1]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)
            kf3f1 = dt*b[0]*self.Lap1_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                    + dt**2*b[0]*a[0,0]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                    + dt**2*b[1]*a[0,1]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)
            kf2f1 = dt*a[0,1]*self.Lap1_kernel_prime(Xt, Xf, Nt,Nf, t3, lambda1, lambda2, un_t, un_f) \
                        +dt*a[1,0]*self.Lap1_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f) \
                        + dt**2*a[1,0]*a[0,0]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t5, lambda1, lambda2, un_t, un_f)\
                        + dt**2*a[1,1]*a[0,1]*self.Lap2_kernel(Xt, Xf, Nt,Nf,t3, lambda1, lambda2, un_t, un_f)

            zfu = tf.zeros((Nt,Nu),dtype=tf.float64)

            k4 = tf.concat( (kf3u3, kf3u2, kf3u1, kf3f3, kf3f2,kf3f1),axis=1)
            k5 = tf.concat( (zfu, kf2u2, kf2u1, kf2f3, kf2f2, kf2f1),axis=1)
            k6 = tf.concat( (zfu, kf1u2, kf1u1, kf1f3, kf1f2, kf1f1),axis=1)

            k = tf.concat((k4,k5,k6),axis=0)

            return k

#    def nlml(self,Xu,Xf,Yu,Yf,dt, hyp1, hyp2, sig_n, kernel_type, jitter=1.0e-10): # negative logarithm marginal-likelihood
#
##        sess1 = tf.Session()
##        sess1.run(tf.global_variables_initializer())
#
#        N = Xu.shape[0] + Xf.shape[0]
#        Nu = Xu.shape[0]
#        K = self.kernel_uf_train(Xu,Xf,hyp1,hyp2,dt) + (sig_n**2+jitter)*tf.eye(N,dtype=tf.float64)
#        self.L = tf.cholesky(K)
#        r = np.concatenate((Yu,Yf),axis=0) - np.concatenate((np.zeros((Nu,1),dtype=np.float64), self.prior_mean_train),axis=0)
#        self.alpha = tf.cholesky_solve(self.L, r)
#        self.sig2_tf = tf.matmul(r, self.alpha, transpose_a=True)/N
#        return 0.5 * N * tf.log(2.0*np.pi*self.sig2_tf)\
#                +tf.reduce_sum(tf.log(tf.diag_part(self.L))) \
#                + N/2.0

    def nlml(self,Xu,Xf,Yu1, Yu2, Yu3, Yf,dt, hyp1, hyp3, hyp5, sig_n, lambda1, lambda2, un_u, un_f, kernel_type, jitter=1.0e-10): # negative logarithm marginal-likelihood

#        sess1 = tf.Session()
#        sess1.run(tf.global_variables_initializer())

#        xf_train = np.linspace(-8.0,8.0,self.Nf+2)[1:-1].reshape((-1,1))
#
#
#        yf_train = self.u_exact(xf_train,0.0, u_exa, t_exa, x_exa, 2)
#        #yf_train = yf_train+np.linalg.cholesky(previous_cov_mat[:Nf,:Nf])@ np.random.randn(Nf,1)
#
#        xu_train = np.array([[-8.0], [8.0]],dtype=np.float64)
#
#        Nu = xu_train.shape[0]
#        Nf = xf_train.shape[0]
#
#
#
#        un_u = self.u_exact(xu_train,init_time,u_exa, t_exa, x_exa, 1)
#        un_f = yf_train

        N = 3*(self.Nu + self.Nf)
        self.K0 = self.kernel_uf_train(Xu,Xf,self.Nu, self.Nf, hyp1,hyp3,hyp5,self.a, self.b, self.c, lambda1, lambda2, un_u, un_f, dt)
#        self.K0 = self.kernel_uf_train(xu_train,xf_train,self.Nu, self.Nf, hyp1,hyp3,hyp5,self.a, self.b, self.c, lambda1, lambda2, un_u, un_f, dt)

        K = self.K0 + (sig_n**2+jitter)*tf.eye(N,dtype=tf.float64)

        self.L = tf.cholesky(K)
        r = tf.concat((Yu1,Yu2,Yu3,Yf,Yf,Yf),axis=0)
        self.alpha = tf.cholesky_solve(self.L, r)
        self.sig2_tf = tf.matmul(r, self.alpha, transpose_a=True)/N
        return 0.5 * N * tf.log(2.0*np.pi*self.sig2_tf)\
                +tf.reduce_sum(tf.log(tf.diag_part(self.L))) \
                + N/2.0

    def training(self, optimizer = 'Adam', num_iter=10001, learning_rate = 5.0e-4, jitter = 1.0e-10):

        tf.reset_default_graph()

        self.hyp1 = tf.exp(tf.Variable(0.0,dtype=np.float64))
        self.hyp3 = tf.exp(tf.Variable(0.0,dtype=np.float64))
        self.hyp5 = tf.exp(tf.Variable(0.0,dtype=np.float64))

        self.lambda1 = tf.Variable(1.0,dtype=np.float64, trainable=False)
        self.lambda2 = tf.Variable(0.1,dtype=np.float64, trainable=False)

        self.xu_train = tf.placeholder(tf.float64, shape=(None,1))
        self.xu_test = tf.placeholder(tf.float64, shape=(None,1))
        self.xf_train = tf.placeholder(tf.float64, shape=(None,1))
        self.xf_test = tf.placeholder(tf.float64, shape=(None,1))
        self.yu_train1 = tf.placeholder(tf.float64, shape=(None,1))
        self.yu_train2 = tf.placeholder(tf.float64, shape=(None,1))
        self.yu_train3 = tf.placeholder(tf.float64, shape=(None,1))
        self.yf_train = tf.placeholder(tf.float64, shape=(None,1))
        self.un_u = tf.placeholder(tf.float64, shape=(None,1))
        self.un_f = tf.placeholder(tf.float64, shape=(None,1))
        self.un_t = tf.placeholder(tf.float64, shape=(None,1))

#        self.hyp = tf.Variable(self.hyp)

        if self.is_noise:
            self.sig_n = tf.exp(tf.Variable(np.log(1.0e-4),dtype=tf.float64))
        else:
            self.sig_n = tf.Variable(0.0,dtype=tf.float64, trainable=False)

#        sess1 = tf.Session()
#        sess1.run(tf.global_variables_initializer())
#        k1,k4 = self.kernel_uf_train(self.xu_train, self.xf_train, self.hyp1, self.hyp3, self.hyp5, self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.dt)
#        kk1, kk4 = self.kernel_f_test(self.xf_test, self.xu_train, self.xf_train, self.hyp1,  self.hyp3,  self.hyp5,  self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.un_t, self.dt)
#        aaa = sess1.run(k4-kk4)
        self.num_iter = num_iter
        self.jitter = jitter
        self.optimizer = optimizer
        self.learning_rate = learning_rate

        init_time = 0.0

        u_simulation = sio.loadmat('burgers.mat')
        u_exa = np.real(u_simulation['usol'])
        t_exa = u_simulation['t'].reshape((-1,1))
        x_exa = u_simulation['x'].reshape((-1,1))

        tt0 = time.time()

#    Nu = 2
#        Nf = 70
#        self.dt = 1.0e-1

        xf_train = np.linspace(-8.0,8.0,self.Nf+2)[1:-1].reshape((-1,1))
        xf_test = x_exa # np.linspace(-8.0,8.0,1000).reshape((-1,1))

        xu_test = np.concatenate((xf_train,xf_test),axis=0)
        yf_train = self.u_exact(xf_train,init_time, u_exa, t_exa, x_exa, 2)
        #yf_train = yf_train+np.linalg.cholesky(previous_cov_mat[:Nf,:Nf])@ np.random.randn(Nf,1)

        xu_train = np.array([[-8.0], [8.0]],dtype=np.float64)

        Nu = xu_train.shape[0]
        Nf = xf_train.shape[0]
        Nt = xf_test.shape[0]

        self.previous_cov_mat = 1e-4*np.eye(3*Nf,dtype=np.float64)

        un_u = self.u_exact(xu_train,init_time,u_exa, t_exa, x_exa, 1)
        un_f = yf_train
        un_t = self.u_exact(xf_test,init_time,u_exa,t_exa,x_exa,2)

        yu_train = []
        yu_train.append(self.u_exact(xu_train,dt*1,u_exa, t_exa, x_exa, 1))
        yu_train.append(self.u_exact(xu_train,dt*(1-1)+self.c[1]*dt,u_exa, t_exa, x_exa, 1))
        yu_train.append(self.u_exact(xu_train,dt*(1-1)+self.c[0]*dt, u_exa, t_exa, x_exa, 1))

        nlml_tf \
             = self.nlml(self.xu_train,self.xf_train, self.yu_train1,self.yu_train2,self.yu_train3, self.yf_train, self.dt, self.hyp1, self.hyp3, self.hyp5, self.sig_n, self.lambda1, self.lambda2, self.un_u, self.un_f, self.kernel_type, self.jitter)
        self.sign_var = self.sig2_tf * self.sig_n**2

        self.K_train = self.kernel_uf_train(self.xu_train, self.xf_train, Nu, Nf, self.hyp1, self.hyp3, self.hyp5, self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.dt)
        self.m_train = tf.matmul(self.K_train,self.alpha)
        L1 =  tf.concat((tf.zeros((3*Nf,3*Nu),dtype=tf.float64),self.previous_cov_mat),axis=1)
        L1 = tf.concat((tf.zeros((3*Nu,3*Nu+3*Nf),dtype=tf.float64),L1),axis=0)
        V1 = tf.linalg.triangular_solve(self.L,tf.transpose(self.K_train))
        V2 = tf.cholesky_solve(self.L, tf.transpose(self.K_train))
        self.var_train = self.sig2_tf*(self.kernel_uf_train(self.xu_train, self.xf_train, Nu, Nf, self.hyp1, self.hyp3, self.hyp5, self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.dt, diag=True)\
                                  - tf.reshape(tf.reduce_sum(V1*V1,axis=0),(-1,1)))
        self.var_train = self.var_train + tf.reshape( tf.diag_part(tf.matmul(tf.matmul(tf.transpose(V2),L1),V2)),(-1,1))

        self.var_train = tf.maximum(self.var_train,tf.zeros((3*Nu+3*Nf,1),dtype=tf.float64) )

#

        k_test_u = self.kernel_u_test(self.xu_test[Nf:,:], self.xu_train, self.xf_train, Nt, Nu, Nf, self.hyp1, self.hyp3, self.hyp5, self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.un_t, self.dt)
        self.m_test_u = tf.matmul(k_test_u,self.alpha)
        V1_test_u = tf.linalg.triangular_solve(self.L,tf.transpose(k_test_u))
        V2_test_u = tf.cholesky_solve(self.L, tf.transpose(k_test_u))
        self.var_test_u = self.sig2_tf * (1.0 - tf.reshape(tf.reduce_sum(V1_test_u*V1_test_u,axis=0),(-1,1))) +self.sign_var
        self.var_test_u = self.var_test_u + tf.reshape( tf.diag_part(tf.matmul(tf.matmul(tf.transpose(V2_test_u),L1),V2_test_u)),(-1,1))

        self.var_test_u = tf.maximum(self.var_test_u,tf.zeros((3*Nt,1),dtype=tf.float64) )

        k_test_u0 = self.kernel_u_test(self.xu_test[:Nf,:], self.xu_train, self.xf_train, Nf, Nu, Nf, self.hyp1, self.hyp3, self.hyp5,self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.un_f,  self.dt)
        self.m_test_u0 = tf.matmul(k_test_u0[:Nf,:],self.alpha)
        V1_test_u0 = tf.linalg.triangular_solve(self.L,tf.transpose(k_test_u0[:Nf,:]))
        V2_test_u0 = tf.cholesky_solve(self.L, tf.transpose(k_test_u0[:Nf,:]))
        self.var_test_u0 = self.sig2_tf * (self.kernel(self.xu_test[:Nf,:],self.xu_test[:Nf,:],Nf, Nf, self.hyp1,equal=True)\
                    - tf.matmul(tf.transpose(V1_test_u0),V1_test_u0)) + self.sign_var* tf.eye(Nf,dtype=tf.float64)
        self.var_test_u0 = self.var_test_u0 + tf.reshape( tf.diag_part(tf.matmul(tf.matmul(tf.transpose(V2_test_u0),L1),V2_test_u0)),(-1,1))
        self.var_test_u0 = tf.maximum(self.var_test_u0,tf.zeros((Nf,Nf),dtype=tf.float64) )

        k_test_f = self.kernel_f_test(self.xf_test, self.xu_train, self.xf_train, Nt, Nu, Nf, self.hyp1,  self.hyp3,  self.hyp5,  self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_u, self.un_f, self.un_t, self.dt)
        self.m_test_f = tf.matmul(k_test_f,self.alpha)
        V1_test_f = tf.linalg.triangular_solve(self.L,tf.transpose(k_test_f))
        V2_test_f = tf.cholesky_solve(self.L, tf.transpose(k_test_f))
        self.var_test_f = self.sig2_tf * (self.kernel_uf_train(self.xf_test, self.xf_test, Nt, Nt, self.hyp1,  self.hyp3, self.hyp5,self.a, self.b, self.c, self.lambda1, self.lambda2, self.un_t, self.un_t,  self.dt,diag=True)[3*Nt:,0:1] \
                   - tf.reshape(tf.reduce_sum(V1_test_f*V1_test_f,axis=0),(-1,1)))  + self.sign_var
        self.var_test_f = self.var_test_f + tf.reshape( tf.diag_part(tf.matmul(tf.matmul(tf.transpose(V2_test_f),L1),V2_test_f)),(-1,1))

        self.var_test_f = tf.maximum(self.var_test_f,tf.zeros((3*Nt,1),dtype=tf.float64) )

        num_t = int(5.0/self.dt)

        u_pred = np.zeros((x_exa.shape[0],num_t+1),dtype=np.float64)
        u_pred[:,0:1] = self.u_exact(x_exa,0.0,u_exa,t_exa,x_exa,2)

        u_interp = np.zeros((x_exa.shape[0],num_t+1),dtype=np.float64)
        u_interp[:,0:1] = self.u_exact(x_exa,0.0,u_exa,t_exa,x_exa,2)

        sig_pred = np.zeros((x_exa.shape[0],num_t+1),dtype=np.float64)

        sig_pred[:,0] = 1.0e-2

        if optimizer == 'Adam':
            optimizer_Adam = tf.train.AdamOptimizer(learning_rate)
            train_op_Adam = optimizer_Adam.minimize(nlml_tf)

            grad1 = tf.gradients(nlml_tf,self.hyp1)[0]
            grad2 = tf.gradients(nlml_tf,self.hyp3)[0]
            grad3 = tf.gradients(nlml_tf,self.hyp5)[0]
            gradn = tf.gradients(nlml_tf,self.sig_n)[0]
            std_train = tf.sqrt(self.var_train)
            std_test_u = tf.sqrt(self.var_test_u)
            std_test_f = tf.sqrt(self.var_test_f)
            std_signal = tf.sqrt(self.sig2_tf)
            std_noise = tf.sqrt(self.sign_var)

        with tf.Session() as sess:
            sess.run(tf.global_variables_initializer())

            for k in np.arange(1,num_t+1):
                print ('\n      t = '+ str(dt*k)+ '  *********************')
                nlml_min = 1.0e16

                yu_train = []
                yu_test = []

                u_interp[:,k:(k+1)] = self.u_exact(xf_test,dt*k, u_exa, t_exa, x_exa, 2)
                yf_test = self.u_exact(xf_test,dt*(k-1),u_exa,t_exa,x_exa,2)

                yu_train.append(self.u_exact(xu_train,dt*k,u_exa, t_exa, x_exa, 1))
                yu_train.append(self.u_exact(xu_train,dt*(k-1)+self.c[1]*dt,u_exa, t_exa, x_exa, 1))
                yu_train.append(self.u_exact(xu_train,dt*(k-1)+self.c[0]*dt, u_exa, t_exa, x_exa, 1))

                yu_test.append(self.u_exact(xf_train,dt*k, u_exa, t_exa, x_exa, 2))
                yu_test.append(u_interp[:,k:(k+1)])
                yu_test.append(self.u_exact(xf_test,dt*(k-1)+self.c[1]*dt, u_exa, t_exa, x_exa, 2))
                yu_test.append(self.u_exact(xf_test,dt*(k-1)+self.c[0]*dt, u_exa, t_exa, x_exa, 2))

                feed_dict = {self.xu_train: xu_train, self.yu_train1: yu_train[0], \
                             self.yu_train2: yu_train[1], self.yu_train3: yu_train[2],\
                            self.xu_test:  xu_test,   \
                            self.xf_train: xf_train, self.yf_train: yf_train,  \
                           self.xf_test: xf_test, self.un_u: un_u, self.un_f: un_f, self.un_t: un_t}

                if k == 1:
                    self.num_iter = 5001
                    print_skip = 5000
                else:
                    self.num_iter = np.maximum(int(5000*self.dt)+1,50)
                    print_skip = self.num_iter-1

                for i in range(self.num_iter):
                    sess.run(train_op_Adam, feed_dict)
                    if i % print_skip == 0:
                        nlml_temp = sess.run(nlml_tf, feed_dict)
                        if nlml_temp < nlml_min:
                            nlml_min = nlml_temp
                            self.mm_train = sess.run(self.m_train, feed_dict)
                            self.ss_train = sess.run(std_train,feed_dict)
                            self.mm_test_u = sess.run(self.m_test_u, feed_dict)
                            self.ss_test_u = sess.run(std_test_u, feed_dict)
                            self.mm_test_f = sess.run(self.m_test_f, feed_dict)
                            self.ss_test_f = sess.run(std_test_f, feed_dict)
                            self.mm_test_u0 = sess.run(self.m_test_u0, feed_dict)
                            self.previous_cov_mat = np.tile(sess.run(self.var_test_u0, feed_dict),(3,3))
                            lambda1_val, lambda2_val, nlml_val, hyp1_val, hyp3_val, hyp5_val, sig_f, sig_n, grad_f1, grad_f3, grad_f5, grad_n= \
                            sess.run([self.lambda1, self.lambda2, nlml_tf, self.hyp1, self.hyp3,  self.hyp5,  std_signal, \
                                        std_noise,grad1,grad2,grad3,gradn],feed_dict)

                            print ('*************************\n')
                            print ('Iter: ', i, '  nlml =', nlml_min, '\n')
                            print ('nlml:   '   , nlml_val)
                            print ('hyp:   ' , [hyp1_val,hyp3_val, hyp5_val])
                            print ('signal std:   ', sig_f)
                            print ('noise std:   ',sig_n)
                            print('grads of nlml over hyp ', [grad_f1, grad_f3, grad_f5])
                            print ('lambda: ', [lambda1_val, lambda2_val])
                            print ('grad of nlml over sig_n', grad_n)

                            print ('Training_err_u3:', np.linalg.norm(self.mm_train[:Nu,0:1]-yu_train[0],2)/np.linalg.norm(yu_train[0],2))
                            print ('Training_err_f3:', np.linalg.norm(self.mm_train[(3*Nu):(3*Nu+Nf),0:1]-yf_train,2)/np.linalg.norm(yf_train,2))

                            print ('Training_err_u2:', np.linalg.norm(self.mm_train[Nu:(2*Nu),0:1]-yu_train[1],2)/np.linalg.norm(yu_train[1],2))
                            print ('Training_err_f2:', np.linalg.norm(self.mm_train[(3*Nu+Nf):(3*Nu+2*Nf),0:1]-yf_train,2)/np.linalg.norm(yf_train,2))

                            print ('Training_err_u1:', np.linalg.norm(self.mm_train[(2*Nu):(3*Nu),0:1]-yu_train[2],2)/np.linalg.norm(yu_train[2],2))
                            print ('Training_err_f1:', np.linalg.norm(self.mm_train[(3*Nu+2*Nf):(3*Nu+3*Nf),0:1]-yf_train,2)/np.linalg.norm(yf_train,2))

                            print ('Test_err_u0:', np.linalg.norm(self.mm_test_u0-yu_test[0],2)/np.linalg.norm(yu_test[0],2))

                            print ('Test_err_u3:', np.linalg.norm(self.mm_test_u[:Nt,0:1]-yu_test[1],2)/np.linalg.norm(yu_test[1],2))
                            print ('Test_err_f3:', np.linalg.norm(self.mm_test_f[:Nt,0:1]-yf_test,2)/np.linalg.norm(yf_test,2))

                            print ('Test_err_u2:', np.linalg.norm(self.mm_test_u[Nt:(2*Nt),0:1]-yu_test[2],2)/np.linalg.norm(yu_test[2],2))
                            print ('Test_err_f2:', np.linalg.norm(self.mm_test_f[Nt:(2*Nt),0:1]-yf_test,2)/np.linalg.norm(yf_test,2))

                            print ('Test_err_u1:', np.linalg.norm(self.mm_test_u[(2*Nt):(3*Nt),0:1]-yu_test[3],2)/np.linalg.norm(yu_test[3],2))
                            print ('Test_err_f1:', np.linalg.norm(self.mm_test_f[(2*Nt):(3*Nt),0:1]-yf_test,2)/np.linalg.norm(yf_test,2))

                yf_train = self.mm_test_u0 #+ np.linalg.cholesky(previous_cov_mat)@np.random.randn(Nf,1)
                u_pred[:,k:(k+1)] = self.mm_test_u[:xf_test.shape[0]]
                sig_pred[:,k:(k+1)] = self.ss_test_u[:xf_test.shape[0]]

                un_u = yu_train[0]
                un_f = self.mm_test_u0
                un_t = self.mm_test_u[:Nt,0:1]

            fig = plt.figure()
            plt.contourf(np.linspace(0.0,5.0,num_t+1), np.ndarray.flatten(x_exa), u_interp, 100, cmap='jet')
            plt.colorbar()
            plt.xlabel('t')
            plt.ylabel('x')
            plt.title('1D Burgers\' equation: Exact solution')
            plt.tight_layout()
            plt.savefig('D-GP-FW-FIG/Exact-Burgers-dt-'+str(100*self.dt)+'.png',dpi=1000)
        #plt.show()
            plt.close(fig)

            fig = plt.figure()
            plt.contourf(np.linspace(0.0,5.0,num_t+1), np.ndarray.flatten(x_exa), u_pred, 100, cmap='jet')
            plt.colorbar()
            plt.xlabel('t')
            plt.ylabel('x')
            plt.title('1D Burgers\' equation: Discrete time GP (solution)')
            plt.tight_layout()
            plt.savefig('D-GP-FW-FIG/D-GP-Burgers-solution-dt-'+str(100*self.dt)+'.png',dpi=1000)
            #plt.show()
            plt.close(fig)

            fig = plt.figure()
            plt.contourf(np.linspace(0.0,5.0,num_t+1), np.ndarray.flatten(x_exa), sig_pred, 100, cmap='jet')
            plt.colorbar()
            plt.xlabel('t')
            plt.ylabel('x')
            plt.title('1D Burgers\' equation: Discrete time GP (std)')
            plt.tight_layout()
            plt.savefig('D-GP-FW-FIG/D-GP-Burgers-std-dt-'+str(100*self.dt)+'.png',dpi=1000)
            #plt.show()
            plt.close(fig)

            fig = plt.figure()
            plt.contourf(np.linspace(0.0,5.0,num_t+1), np.ndarray.flatten(x_exa), np.abs(u_interp-u_pred), 100, cmap='jet')
            plt.colorbar()
            plt.xlabel('t')
            plt.ylabel('x')
            plt.title('1D Burgers\' equation: Discrete time GP (absolute error)')
            plt.tight_layout()
            plt.savefig('D-GP-FW-FIG/D-GP-Burgers-Ab-err-dt-'+str(100*self.dt)+'.png',dpi=1000)
            #plt.show()
            plt.close(fig)

            np.savetxt('D-GP-FW-FIG/exact_u.txt', u_exa, fmt='%10.5e')
            np.savetxt('D-GP-FW-FIG/interp_u.txt', u_interp, fmt='%10.5e')

            np.savetxt('D-GP-FW-FIG/predicted_u_dt_'+str(100*self.dt)+'.txt', u_pred, fmt='%10.5e')
            np.savetxt('D-GP-FW-FIG/predicted_std_dt_'+str(100*self.dt)+'.txt', sig_pred, fmt='%10.5e')

            u_error = np.linalg.norm(u_interp.reshape((-1,1))-u_pred.reshape((-1,1)))/np.linalg.norm(u_interp.reshape((-1,1)))
            print('u_error= ', u_error)
            np.savetxt('D-GP-FW-FIG/u_error_dt_'+str(100*self.dt)+'.txt', [u_error], fmt='%10.5e' )
            tt1 = time.time()

            print ('CPU time ', tt1-tt0)

a = np.array([[0.25, 0.25-np.sqrt(3.0)/6.0], [0.25+np.sqrt(3.0)/6.0, 0.25]],dtype=np.float64)
b = np.array([0.5,0.5],dtype=np.float64)
c = np.array([0.5-np.sqrt(3.0)/6.0, 0.5+np.sqrt(3.0)/6.0],dtype=np.float64)
Nu = 2
Nf = 70
dt = 5.0e-3

GP_instance = one_GP()
GP_instance.model(dt,a,b,c,Nu,Nf,is_noise=True)
GP_instance.training(num_iter=5001,jitter=0.0)

del GP_instance

reason: import numpy,import scipy


text:
from sympy import primefactors
from fpack import primelist_till_x as primo
from time import time

def phi_n(x):
    l = primefactors(x)
    if len(l) == 2:  # because we want to maximize phi(n)
        phi_n = x
        for ele in l:
            phi_n *= 1 - 1 / ele
        return int(phi_n)
    else:
        return 0

def p(x, y):
    return sorted(list(str(x))) == sorted(list(str(y)))

t = time()
# i = 2
# low = 10**7
# num = 0
# while i < 10**7:
#     print(i)
#     if p(i,phi_n(i)) and i/phi_n(i) < low:
#         low = i/phi_n(i)
#         num = i
#     i += 1
i = 10**7 - 1
low = 10**7
num = 0
while True:
    print(i)
    if len(primefactors(i)) == 2:
        if (
            p(i, (primefactors(i)[0] - 1) * (primefactors(i)[1] - 1))
            and primefactors(i)[0]
            * primefactors(i)[1]
            / (primefactors(i)[0] - 1)
            * (primefactors(i)[1] - 1)
            < low
        ):
            low = (
                primefactors(i)[0]
                * primefactors(i)[1]
                / (primefactors(i)[0] - 1)
                * (primefactors(i)[1] - 1)
            )
            num = i
            print("===============>", i)
    i -= 2

# print(primies, time() - t)

reason: from sympy


text:
from dlib import *
import numpy as np

import sys
sys.path = ['./superfast/build'] + sys.path
import superfast
import tools

# img = load_grayscale_image(sys.argv[1])

img = load_grayscale_image('./images/find_page/paper01.jpg')

# What about this image?  Need to do something to fix it
# img = load_grayscale_image('./images/find_page/tissue_01.jpg')

ht = hough_transform(300)

img = sub_image(gaussian_blur(img,3))
img = resize_image(img, ht.size, ht.size)

win1 = image_window(img)

img = equalize_histogram(img)
x,y = sobel_edge_detector(img)
edges = suppress_non_maximum_edges(x,y)
win3 = image_window(edges)
normalize_image_gradients(x,y)
superfast.discard_wacky_edge_groups(edges, x, y)
edges = hysteresis_threshold(edges, partition_pixels(edges), 50)
win4 = image_window(edges)
himg = ht(edges)

hits = ht.find_strong_hough_points(himg, hough_count_thresh=ht.size/5, angle_nms_thresh=15, radius_nms_thresh=10)

print("hough hits found: ", len(hits))
# keep just the 50 best lines
if (len(hits) > 50):
    hits = hits[0:50]

# NEW!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Use tools.find_hough_boxes_less_simple(), which now also checks if all the corners are connected.
line_pixels = ht.find_pixels_voting_for_lines(edges, hits, 4,4)
boxes = tools.find_hough_boxes_less_simple(ht, hits, line_pixels)

if len(boxes) > 0:
    c1,c2,c3,c4,area,idx1,idx2,idx3,idx4 = boxes[0]
    win1.add_overlay(line(c1,c2))
    win1.add_overlay(line(c2,c3))
    win1.add_overlay(line(c3,c4))
    win1.add_overlay(line(c4,c1))

    page = extract_image_4points(img, [c1,c2,c3,c4], 200,200)
    win_page = image_window(page)

    # NEW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    # plot the pixels that make up this box.
    cimg = convert_image(edges, dtype='rgb_pixel')
    for l in [line_pixels[i] for i in [idx1, idx2, idx3, idx4]]:
        for p in l:
            cimg[p.y,p.x,:] = (255,0,0)
    win4.set_image(cimg)

input("hit enter to exit")


text:
'''
Run SExtractor on all the original images to identify saturating sources.

Find the XY position of the possible crosstalk and convert to WCS coordinates.

create a file including the world coords of possible regions affected by crosstalk.
'''

# Author:  Igor Andreoni <andreoni@caltech.edu>,
# License: MIT.

import numpy as np
import subprocess
from subprocess import call
import os
from astropy.io import ascii
from astropy.io import fits
from astropy import wcs
from astropy.table import Table
import pdb

def flag_crosstalk(path_images, image_name, callsex, path_config,
                   set_saturation=True, saturation_adu=5000, scale = 0.2637):

    # Initialise the lists:
    listRA = []
    listDEC = []
    listRADII = []

    # Initialise a list of saturating sources (different catalog)
    listSATra = []
    listSATdec = []

    # Run SExtractor to find the sources flagged as saturating.
    catname = f"{path_images}/{image_name.replace('.fits','_xtalk_temp.cat')}"

    sexcommand = f"{callsex} {path_images}/{image_name} -c {path_config}/config_crosstalk.sex -CATALOG_NAME {catname} -PARAMETERS_NAME {path_config}/sex_crosstalk.param"

    # Fix the saturation level in case the user does not want to use the one given in the header:
    if set_saturation:
        sexcommand = sexcommand + f" -SATUR_LEVEL {saturation_adu} -SATUR_KEY  NONE "
    else:
        sexcommand = sexcommand + " -SATUR_KEY SATURATION "

    # Call the SExtractor command:
    os.system(sexcommand)

    # Select the saturating sources (flag = 4,5,6,7)
    saturation_tbl = ascii.read(catname, format='sextractor')
    saturation_tbl = saturation_tbl[saturation_tbl['FLAGS'] >= 4]
    saturation_tbl = saturation_tbl[saturation_tbl['FLAGS'] <= 7]

    #Return if there are no saturating sources
    if len(saturation_tbl) == 0:
        print('No saturating sources. The crosstalk flagging file will be empty.')
        return

    # Read in the image (fits) and get exact size information
    image_hdul = fits.open(f"{path_images}/{image_name}")
    image_data = image_hdul[0].data
    NAXIS2, NAXIS1 = len(image_data[:,0]), len(image_data[0,:])

    # Get the WCS solution
    w = wcs.WCS(image_hdul[0].header)

    # Initialize new table for the crosstalk
    crosstalk_tbl = Table([[],[],[],[],[]], names = ('RA_sat','Dec_sat', 'RA_xtalk', 'Dec_xtalk', 'radius'))

    # For each saturating star, compute the location of the crosstalk
    for Xsat, Ysat, RAsat, DECsat, kr, a_image in zip(saturation_tbl['X_IMAGE'], saturation_tbl['Y_IMAGE'],
                                                      saturation_tbl['X_WORLD'], saturation_tbl['Y_WORLD'],
                                                      saturation_tbl['KRON_RADIUS'],
                                                      saturation_tbl['A_IMAGE']):
        xcross = (NAXIS1 / 2.) + ((NAXIS1 / 2.) - Xsat)
        ycross = Ysat
        world_coords_xtalk = w.wcs_pix2world(np.array([[xcross, ycross]]), 1)
        RAcross, DECcross = world_coords_xtalk[0][0], world_coords_xtalk[0][1]
        newradius = kr * a_image
        # Add a row to the output table, converting the radius from px to arcsec
        crosstalk_tbl.add_row([RAsat, DECsat, RAcross, DECcross, newradius*scale])
        print(Xsat, Ysat, xcross,ycross)

    # File that will include RA, Dec, RADIUS of the regions to flag as possible crosstalk, along with bright sources.
    crosstalk_filename=f"{path_images}/{image_name.replace('.fits','_crosstalk.csv')}"
    crosstalk_tbl.write(crosstalk_filename, format='csv')

    return crosstalk_tbl

if __name__ == '__main__':
    path_images = 'test'
    image_name = '4hr.g.ut151219.503748_23.fits'
    callsex = '/Users/igor/Software/ureka/Ureka/bin/sex'
    path_config = '.'
    crosstalk = flag_crosstalk(path_images, image_name, callsex, path_config,
                               set_saturation = True, saturation_adu=5000,
                               scale=0.2637)

max_stars_repo_name: igorandreoni/DECam-codes

max_issues_repo_name: igorandreoni/DECam-codes

max_forks_repo_name: igorandreoni/DECam-codes

reason: import numpy,from astropy


text:
import numpy as np
from scipy import sparse as sp
# note that we use the scikit-learn bundled version of joblib
from sklearn.externals import joblib

# --------------------------------------------------------------

def build_centroids( X, partition, k ):
	"""
	Build a set of K centroids based on the specified partition memberships.
	"""
	# NB: need to convert to a numpy array before we do this...
	memberships = np.array(partition)
	n_features = X.shape[1]
	centroids = np.empty((k, n_features), dtype=np.float64)
	for i in range(k):
		center_mask = memberships == i
		if sp.issparse(X):
			center_mask = np.where(center_mask)[0]
		centroids[i] = X[center_mask].mean(axis=0)
	return centroids

def clustermap_to_partition( cluster_map, doc_ids ):
	"""
	Convert a dictionary representing a clustering into a partition.
	"""
	cluster_names = list(cluster_map.keys())
	cluster_names.sort()
	# build document map
	partition = []
	doc_map = {}
	for i in range( len(doc_ids) ):
		partition.append( -1 )
		doc_map[ doc_ids[i] ] = i
	# now create partition
	for cluster_index in range( len(cluster_map) ):
		for doc_id in cluster_map[cluster_names[cluster_index]]:
			partition[doc_map[doc_id]] = cluster_index
	return partition

# --------------------------------------------------------------

def save_term_rankings( out_path, term_rankings, labels = None ):
	"""
	Save a list of multiple term rankings using Joblib.
	"""
	# no labels? generate some standard ones
	if labels is None:
		labels = []
		for i in range( len(term_rankings) ):
			labels.append( "C%02d" % (i+1) )
	joblib.dump((term_rankings,labels), out_path )

def load_term_rankings( in_path ):
	"""
	Load a list of multiple term rankings using Joblib.
	"""
	(term_rankings,labels) = joblib.load( in_path )
	return (term_rankings,labels)

def save_nmf_factors( out_path, W, H, doc_ids ):
	"""
	Save a NMF factorization result using Joblib.
	"""
	joblib.dump((W,H,doc_ids), out_path )

def load_nmf_factors( in_path ):
	"""
	Load a NMF factorization result using Joblib.
	"""
	(W,H,doc_ids) = joblib.load( in_path )
	return (W,H,doc_ids)

def save_partition( out_path, partition, doc_ids ):
	"""
	Save a disjoint partition (clustering) result using Joblib.
	"""
	joblib.dump((partition,doc_ids), out_path )

def load_partition( in_path):
	"""
	Load a disjoint partition (clustering) result using Joblib.
	"""
	(partition,doc_ids) = joblib.load( in_path )
	return (partition,doc_ids)

reason: import numpy,from scipy


text:
import numpy as np
import math
from struct import unpack
#import pylab as pl

from functools import wraps
import functools
from collections import namedtuple
from typing import NamedTuple
from pprint import pprint
import logging

from .logging import AddLocationFilter
from .constants import eodmarker,eotmarker,magicmarker

logger = logging.getLogger(__name__)
logflt = AddLocationFilter()
logger.addFilter(logflt)

class decode:
	"""Decorator decoder class for 32-bit data words from TRAPconfig

	The constructor takes a description of a TRAP data word in the format
	used in the TRAP User Manual. It can then be used to decorate functions
	to help with the parsing of data words according to	this format. If the
	parsing succeeds, the function is called with an additional argument that
	contains the extracted fields as a namedtuple. An assertion error is
	raised dif the parsing fails."""

	def __init__(self, pattern):

		fieldinfo = dict()

		# calculate bitmask and required shift for every character
		for i,x in enumerate(x for x in pattern if x not in ": "):
			p = 31-i # bit position
			if x not in fieldinfo:
				fieldinfo[x] = dict(mask=1<<p, shift=p)
			else:
				fieldinfo[x]['mask'] |= 1<<p
				fieldinfo[x]['shift'] = p

		# remember bits marked as '0' or '1' for validation of the dword
		zero_mask = fieldinfo.pop('0')['mask'] if '0' in fieldinfo else 0
		one_mask  = fieldinfo.pop('1')['mask'] if '1' in fieldinfo else 0
		self.validate_value = one_mask
		self.validate_mask = zero_mask | one_mask

		# create a named tuple based on the field names
		self.dtype = namedtuple("decoded_fields", fieldinfo.keys() )

		# flatten the fields info to a tuple for faster retrieval
		# the field name is not stritcly necessary, but might help to debug
		self.fields = tuple(
		  (k,fieldinfo[k]['mask'],fieldinfo[k]['shift']) for k in fieldinfo )

	def __call__(self,func):

		@wraps(func)
		def wrapper(ctx, dword):
			assert( (dword & self.validate_mask) == self.validate_value)
			return func(ctx,dword,self.decode(dword))

		return wrapper

	def decode(self,dword):
		return self.dtype(*[ (dword & x[1]) >> x[2] for x in self.fields ])

class describe:
	"""Decorator to generate messages about dwords

	Probably this function is overkill and should be replaced with a single
	log statement in the decorated functions."""

	def __init__(self, fmt):
		self.format = fmt

	def __call__(self,func):

		@wraps(func)
		def wrapper(ctx,dword,fields=None):

			if fields is None:
				retval = func(ctx,dword)
				fielddata = {}

			else:
				retval = func(ctx,dword,fields)
				fielddata = fields._asdict()

			if retval is True or retval is None:
				retval = dict()

			if 'description' not in retval:
				logger.info(self.format.format(dword=dword, **fielddata, ctx=ctx))
				# retval['description'] = msg

			return retval

		return wrapper

ParsingContext = namedtuple('ParsingContext', [
  'major', 'minor', 'nhw', 'sm', 'stack', 'layer', 'side', #from HC0
  'ntb', 'bc_counter', 'pre_counter', 'pre_phase', # from HC1
  'SIDE', 'HC', 'VER', 'det', ## derived from HCx
  'rob', 'mcm', ## from MCM header
  'store_digits', ## links to helper functions/functors
  'event', ## event number
])

# ------------------------------------------------------------------------
# Generic dwords

@describe("SKP ... skip parsing ...")
def skip_until_eod(ctx, dword):
	assert(dword != eodmarker)
	return dict(readlist=[[parse_eod, skip_until_eod]])

@describe("TRK tracklet")
def parse_tracklet(state, dword):
	assert(dword != eotmarker)
	return dict(readlist=[[parse_tracklet, parse_eot]])

@describe("EOT")
def parse_eot(ctx, dword):
	assert(dword == eotmarker)
	return dict(readlist=[[parse_eot, parse_hc0]])

@describe("EOD")
def parse_eod(ctx, dword):
	assert(dword == eodmarker)
	return dict(readlist=[[parse_eod]])

# ------------------------------------------------------------------------
# Half-chamber headers

@decode("xmmm : mmmm : nnnn : nnnq : qqss : sssp : ppcc : ci01")
@describe("HC0 {ctx.HC} ver=0x{m:X}.{n:X} nw={q}")
def parse_hc0(ctx, dword, fields):

	ctx.major = fields.m  # (dword >> 24) & 0x7f
	ctx.minor = fields.n  # (dword >> 17) & 0x7f
	ctx.nhw   = fields.q  # (dword >> 14) & 0x7
	ctx.sm    = fields.s  # (dword >>  9) & 0x7
	ctx.layer = fields.p  # (dword >>  6) & 0x1f
	ctx.stack = fields.c  # (dword >>  3) & 0x3
	ctx.side  = fields.i  # (dword >>  2) & 0x1

	ctx.det = 18*ctx.sm + 6*ctx.stack + ctx.layer

	# An alternative to update the context - which one is easier to read?
	# (ctx.major,ctx.minor,ctx.nhw,ctx.sm,ctx.layer,ctx.stack,ctx.side) = fields

	# Data corruption seen with configs around svn r5930 -> no major/minor info
	# This is a crude fix, and the underlying problem should be solved ASAP
	if ctx.major==0 and ctx.minor==0 and ctx.nhw==0:
		ctx.major = 0x20 # ZS
		ctx.minor = 0
		ctx.nhw   = 2

	# set an abbreviation for further log messages
	side = 'A' if fields.i==0 else 'B'
	ctx.HC   = f"{fields.s:02}_{fields.c}_{fields.p}{side}"

	readlist = list()
	for i in range(ctx.nhw):
		# check additional HC header in with HC1 last, because HC2 and HC3
		# appear like HC1 with the (invalid) phase >= 12. This order avoids
		# this ambiguity.
		readlist.append([parse_hc3, parse_hc2, parse_hc1])

	readlist.append([parse_mcmhdr])
	return dict(readlist=readlist)

@decode("tttt : ttbb : bbbb : bbbb : bbbb : bbpp : pphh : hh01")
@describe("HC1 tb={t} bc={b} ptrg={p} phase={h}")
def parse_hc1(ctx, dword, fields):

	ctx.ntb         = fields.t  # (dword >> 26) & 0x3f
	ctx.bc_counter  = fields.b  # (dword >> 10) & 0xffff
	ctx.pre_counter = fields.p  # (dword >>  6) & 0xF
	ctx.pre_phase   = fields.h  # (dword >>  2) & 0xF

@decode("pgtc : nbaa : aaaa : xxxx : xxxx : xxxx : xx11 : 0001")
@describe("HC2 - filter settings")
def parse_hc2(ctx, dword, fields):
	pass

@decode("ssss : ssss : ssss : saaa : aaaa : aaaa : aa11 : 0101")
@describe("HC3 - svn version {s} {a}")
def parse_hc3(ctx, dword, fields):
	pass

# ------------------------------------------------------------------------
# MCM headers

@decode("1rrr : mmmm : eeee : eeee : eeee : eeee : eeee : 1100")
@describe("MCM {r}:{m:02} event {e}")
def parse_mcmhdr(ctx, dword, fields):

	ctx.rob = fields.r
	ctx.mcm = fields.m

	if ctx.major & 0x20:   # Zero suppression
		return dict(readlist=[[parse_adcmask]])

	else:  # No ZS -> read 21 channels, then expect next MCM header or EOD
		readlist = list()
		for i in range ( 21 * (ctx.ntb // 3) ):
			readlist.append([parse_adcdata])

		readlist.append([parse_mcmhdr, parse_eod])
		return dict(readlist=readlist)

@decode("nncc : cccm : mmmm : mmmm : mmmm : mmmm : mmmm : 1100")
def parse_adcmask(ctx, dword, fields):
	desc = "MSK "
	count = 0
	readlist = list()

	adcdata = np.zeros(ctx.ntb, dtype=np.uint16)

	for ch in range(21):
		if ch in [9,19]:
			desc += " "

		if fields.m & (1<<ch):
			count += 1
			desc += str(ch%10)
			for tb in range ( 0, ctx.ntb , 3 ):
				readlist.append([parse_adcdata(channel=ch, timebin=tb, adcdata=adcdata)])
		else:
			desc += "."

	desc += f"  ({~fields.c & 0x1F} channels)"
	readlist.append([parse_mcmhdr, parse_eod])
	assert( count == (~fields.c & 0x1F) )

	logger.info(desc)
	return dict(readlist=readlist)
	# return dict(description=desc, readlist=readlist)

# ------------------------------------------------------------------------
# Raw data
class parse_adcdata:
	"""ADC data parser

	To parse ADC data, we need to know the channel number and the timebins
	in this dword. I don't think this data should be kept in the context.
	The parser for the MCM header / adcmask therefore stores it in the parser
	for the ADC data word. This parser therefore has to be a callable object.
	"""

	def __init__(self, channel, timebin, adcdata=None):
		self.channel = channel
		self.timebin = timebin
		self.adcdata = adcdata
		self.__name__ = "parse_adcdata"

	# @decode("xxxx:xxxx:xxyy:yyyy:yyyy:zzzz:zzzz:zzff")
	def __call__(self, ctx, dword):
		x = (dword & 0xFFC00000) >> 22
		y = (dword & 0x003FF000) >> 12
		z = (dword & 0x00000FFC) >>  2
		f = (dword & 0x00000003) >>  0

		msg = "ADC "
		msg += f"ch {self.channel:2} " if self.timebin==0 else " "*6
		msg += f"tb {self.timebin:2} (f={f})   {x:4}  {y:4}  {z:4}"

		logger.info(msg)

		if self.adcdata is not None and ctx.store_digits is not None:
			# store the ADC values in the reserved array
			for i,adc in enumerate((x,y,z)):
				if self.timebin+i < len(self.adcdata):
					self.adcdata[self.timebin+i] = adc

			# if this is the last dword for this channel -> store the digit
			if self.timebin+3 >= len(self.adcdata):
				ctx.store_digits(ctx.event, ctx.det, ctx.rob, ctx.mcm,
				                 self.channel, self.adcdata)

		return dict()

# ------------------------------------------------------------------------
class LinkParser:

	# end_of_tracklet = 0x10001000
	# end_of_data     = 0x00000000

	#Defining the initial variables for class
	def __init__(self, store_digits = None):
		self.ctx = ParsingContext
		self.ctx.event = 0
		self.ctx.store_digits = store_digits

	def next_event(self):
		self.ctx.event += 1

	def process(self,linkdata, linkpos=-1):
		'''Initialize parsing context and process data from one optical link.

		Parameter: linkdata = iterable list of uint32 data words
        '''

		self.ctx.current_linkpos = linkpos

		self.readlist = [ list([parse_tracklet, parse_eot]) ]
		self.process_linkdata(linkdata)

	def process_linkdata(self, linkdata):

		for dword in linkdata:

			self.ctx.current_linkpos += 1
			self.ctx.current_dword = dword

			logflt.where = f"{self.ctx.current_linkpos:06x} {dword:08x}  "

			# Debugging:
			# self.dump_readlist()

			try:
				# for fct in self.readlist[i]:
				for fct in self.readlist.pop(0):

					# The function can raise an AssertionError to signal that
					# it does not understand the dword
					try:
						 result = fct(self.ctx,dword)

					except AssertionError as ex:
						continue

					if not isinstance(result, dict):
						break

					if 'readlist' in result:
						self.readlist.extend(result['readlist'])

					# the function handled the dword -> we are done
					# if 'description' in result:
					# 	print(f"{ctx.current_dword:06x} {dword:08x}  ", end="")
					# 	print(result['description'])

					break

				else:
					logger.error(logflt.where + "NO MATCH - testing possibilities")
					# check_dword(dword)

					# skip everything until EOD
					self.readlist.extend([[parse_eod, skip_until_eod]])
					continue

			except IndexError:
				logger.error(logflt.where + "extra data after end of readlist")
				break

	def dump_readlist(self):
		for j,l in enumerate(self.readlist):
			print( [ f.__name__ for f in self.readlist[j] ] )

def check_dword(dword):

	ctx = dict()

	parsers = [ parse_tracklet, parse_eot, parse_eod,
	  parse_hc0, parse_hc1, parse_hc2, parse_hc3,
	  parse_mcmhdr, parse_adcmask, parse_adcdata(-1,-1) ]

	for p in parsers:
		try:
			p(ctx,dword)
		except AssertionError:
			continue

		# log.debug("Match:", p.__name__)


text:
import cv2
import pickle
import czifile
import numpy as np
import os
from time import time

def read_preprocessed(file_name, cell_id, image_number):
    image = cv2.imread(f"images/{file_name}_{cell_id}/Images{image_number}.jpg")
    return image

def get_landmarks(image):
    saturation_epsilon = 20
    no_grayscale = image.copy()
    for i in range(0, no_grayscale.shape[0]):
        for j in range(0, no_grayscale.shape[1]):

            distance = abs(int(no_grayscale[i][j][0]) - int(no_grayscale[i][j][1]))
            distance2 = abs(int(no_grayscale[i][j][1]) - int(no_grayscale[i][j][2]))

            if no_grayscale[i][j][0] == no_grayscale[i][j][1] == no_grayscale[i][j][2]:
                no_grayscale[i][j][0], no_grayscale[i][j][1], no_grayscale[i][j][2] = (0, 0, 0)
            elif distance <= saturation_epsilon and distance2 <= saturation_epsilon and no_grayscale[i][j][2] <= 160:
                no_grayscale[i][j][0], no_grayscale[i][j][1], no_grayscale[i][j][2] = (0, 0, 0)
            else:
                #print(no_grayscale[i][j][0], no_grayscale[i][j][1],no_grayscale[i][j][2])
                no_grayscale[i][j][0], no_grayscale[i][j][1],no_grayscale[i][j][2] = (0, 0, 255)

    #cv2.imshow('noname', no_grayscale)
    landmarks = cv2.Canny(no_grayscale, 100, 200)
    return landmarks

def get_centroid(image_shape):
    avg_x = 0.0
    avg_y = 0.0
    points_count = 0
    for i in range(0, image_shape.shape[0]):
        for j in range(0, image_shape.shape[1]):
            if image_shape[i][j] > 0:
                avg_x += i
                avg_y += j
                points_count += 1
    print(int(avg_x/points_count), int(avg_y/points_count))
    return (int(avg_x/points_count), int(avg_y/points_count))

def translate(to, centroid, image_shape):

    for i in range(0, image_shape.shape[0]):
        for j in range(0, image_shape.shape[1]):
            if image_shape[i][j] > 0:
                image_shape[i][j] = 0
                new_i = int(i - centroid[0]) + to[0]
                new_j = int(j - centroid[1]) + to[1]
                image_shape[new_i][new_j] = 255

    return image_shape

# distance between pixels
def in_rangeof(landmarks, point, distance=3):
    for landmark in landmarks:
        if abs(landmark[0]-point[0])+ abs(landmark[1]-point[1]) < distance:
            return True
    return False

if __name__ == "__main__":

    "(start_slice, cell_id, data_file)"
    images_data = [(29, 0, "p3.dat"),
                (29, 1, "p3.dat"),
                (25, 2, "p3.dat"),
                (35, 0, "p6.dat")]
    z_resolution = 4.0
    first = True
    recommended = 100
    shapes = []
    for element in images_data:
        print("---NEW image data---")
        shape_data = []
        start_slice = element[0]
        cell_id = element[1]
        file_name = element[2]

        preprocessed_image = read_preprocessed(file_name, cell_id, start_slice)

        land_marks = get_landmarks(preprocessed_image)
        centroid = get_centroid(land_marks)
        cv2.imwrite(f"images/landmarks/image{start_slice}.jpg", land_marks)

        translated = translate((300, 300), centroid, land_marks)
        new_imagesize = (600, 600)
        cropped = translated[0:new_imagesize[0], 0:new_imagesize[1]]

        reduced_landmarks = np.zeros((600, 600))
        current_landmarks = []

        default_distance = 12
        while len(current_landmarks) != recommended:
            for i in range(0, cropped.shape[0]):
                for j in range(0, cropped.shape[1]):
                    if cropped[i][j] > 0:
                        point = (i, j)
                        if not in_rangeof(current_landmarks, point, default_distance):
                            current_landmarks.append((i,j, 0*z_resolution))
                            reduced_landmarks[i][j] = 255
                            if len(current_landmarks) == recommended:
                                break
                if len(current_landmarks) == recommended:
                    break
            if len(current_landmarks) == recommended:
                break
            default_distance -= 1
            if default_distance < 0:
                print("fail")
                break

        cv2.imwrite(f"images/reduced/image{start_slice}{time()}.jpg", reduced_landmarks)
        print("landmarks: " + str(len(current_landmarks)))

        current_landmarks.sort()
        print(current_landmarks)
        shape_data += current_landmarks

        slices_count = 15
        for k in range(start_slice+1, start_slice+slices_count):

            cropped = None
            if not os.path.exists(f"images/cropped_cache/{file_name}_{cell_id}_r{k}_cache.dat"):
                preprocessed_image = read_preprocessed(file_name, cell_id, k)

                land_marks = get_landmarks(preprocessed_image)
                cv2.imwrite(f"images/landmarks/image{k}.jpg", land_marks)

                translated = translate((300, 300), centroid, land_marks)

                new_imagesize = (600, 600)
                cropped = translated[0:new_imagesize[0], 0:new_imagesize[1]]
                with open(f"images/cropped_cache/{file_name}_{cell_id}_r{k}_cache.dat", 'wb+') as f:
                    pickle.dump(cropped, f)
            else:
                with open(f"images/cropped_cache/{file_name}_{cell_id}_r{k}_cache.dat", mode="rb+") as f:
                    cropped = pickle.load(f)

            reduced_landmarks = np.zeros((600, 600))
            current_landmarks = []

            default_distance = 12
            while len(current_landmarks) != recommended:
                for i in range(0, cropped.shape[0]):
                    for j in range(0, cropped.shape[1]):
                        if cropped[i][j] > 0:
                            point = (i, j)
                            if not in_rangeof(current_landmarks, point, default_distance):
                                current_landmarks.append((i,j, (k-start_slice)*z_resolution))
                                reduced_landmarks[i][j] = 255
                                if len(current_landmarks) == recommended:
                                    break
                    if len(current_landmarks) == recommended:
                        break
                if len(current_landmarks) == recommended:
                    break
                default_distance -= 1
                if default_distance < 0:
                    print("fail")
                    break

            print("landmarks: " + str(len(current_landmarks)))
            current_landmarks.sort()
            print(current_landmarks)
            shape_data += current_landmarks
            cv2.imwrite(f"images/reduced/image{k}{time()}.jpg", reduced_landmarks)

        shapes.append(shape_data)

    if True:
        with open("shapes.dat", mode="wb+") as f_shapes:
            pickle.dump(shapes, f_shapes)

max_stars_repo_name: PascalSavage/procedural-cell-generation

max_issues_repo_name: PascalSavage/procedural-cell-generation

max_forks_repo_name: PascalSavage/procedural-cell-generation