text:
from . import Layers as L
# from . import network # as Net

import matplotlib.pyplot as plt
import numpy as np

def visualizeNeuronPotentials(network, neuronHistory):

    neuronList = []
    for ix, layer in enumerate(network.layers):
        maxNeurons = (ix, len(layer.postneurons))
        neuronList.append(maxNeurons)
    numNeurons = np.array(neuronList)
    numberLayers = len(neuronHistory[0])
    maxNeuronsPerLayer = np.max(numNeurons[:,1])

    t = np.linspace(0, len(neuronHistory), len(neuronHistory))
    fig, axs = plt.subplots(numberLayers, maxNeuronsPerLayer,
                            sharex='all', sharey='all', figsize=(maxNeuronsPerLayer * 2, numberLayers * 2))

    for layerIx, layer in enumerate(network.layers):
        for neuronId in range(numNeurons[layerIx,1]):
            layerHistory = np.array(neuronHistory[:, layerIx])
            excitatoryData = []
            inhibitoryData = []
            for time in range(len(layerHistory)):
                excitatoryData.append(layerHistory[time][neuronId].neuron[0])
                inhibitoryData.append(layerHistory[time][neuronId].neuron[1])

            excitatoryData = np.array(excitatoryData)
            inhibitoryData = np.array(inhibitoryData)

            title = 'Hidden neuron ' + str(neuronId)
            if layerIx == 0:
                title = 'Input neuron ' + str(neuronId)
            if layerIx == network.numLayers - 1: #len(layer) - 1:
                title = 'Output neuron ' + str(neuronId)

            if numNeurons[layerIx, 1] > 1:
                axs[layerIx, neuronId].set_title(title)
                line1, = axs[layerIx, neuronId].plot(t, excitatoryData, dashes=[6, 2], label='Neuron Potential')
                line2, = axs[layerIx, neuronId].plot(t, inhibitoryData, dashes=[4, 2, 8, 3], label='Neuron Inhibition')
                axs[layerIx, neuronId].set_ylabel('potential')
              #  axs[layerIx, neuronId].legend()

                if neuronId == numNeurons[layerIx,1]:
                    axs[layerIx, neuronId].set_xlabel('time')

                if layerIx == 0 and neuronId == 0:
                    axs[layerIx, neuronId].legend()

            else: # there is only one neuron in this layer
                axs[layerIx].set_title(title)
                line1, = axs[layerIx].plot(t, excitatoryData, dashes=[6, 2], label='Neuron Potential')
                line2, = axs[layerIx].plot(t, inhibitoryData, dashes=[4, 2, 8, 3], label='Neuron Inhibition')
                axs[layerIx].set_ylabel('potential')
              #  axs[layerIx].legend()

                if neuronId == numNeurons[layerIx,1]:
                    axs[layerIx].set_xlabel('time')

                if layerIx == 0 and neuronId == 0:
                    axs[layerIx].legend()

   # plt.show()

def visualizeSynapseMatrix(synapseMatrixHistory, title=""):

    exampleMatrix = synapseMatrixHistory[0]
    noPre, noPost = exampleMatrix.getNumberOfConnectedNeurons()

    ventricleCounts = np.zeros((noPre, noPost, len(synapseMatrixHistory)))
    criftTransmitterCounts = np.zeros((noPre, noPost, len(synapseMatrixHistory)))
    freeRezeptorCounts = np.zeros((noPre, noPost, len(synapseMatrixHistory)))
    acceptedTransmitterCounts = np.zeros((noPre, noPost, len(synapseMatrixHistory)))

    t = np.linspace(0, len(synapseMatrixHistory), len(synapseMatrixHistory))
    fig, axs = plt.subplots(noPre, noPost, sharex='all', sharey='all', figsize=(noPost * 2, noPre * 2))

    fig.suptitle(title)

    # Synapse Data
    for time in range(len(synapseMatrixHistory)):
        for preIx in range(noPre):
            for postIx in range(noPost):
                synapseMatrix = synapseMatrixHistory[time].synapseMatrix
                synapse = synapseMatrix[preIx, postIx].synapse
                ventricleCounts[preIx, postIx,time] = synapse[0]
                criftTransmitterCounts[preIx, postIx,time] = synapse[1]
                freeRezeptorCounts[preIx, postIx,time] = synapse[2]
         #       acceptedTransmitterCounts[preIx, postIx] = synapse[3]

    for preIx in range(noPre):
        for postIx in range(noPost):
            axs[preIx, postIx].set_title('Synapse ' + str(preIx) + ' to ' + str(postIx))  # 0 to 1')
            line1, = axs[preIx, postIx].plot(t, ventricleCounts[preIx, postIx], dashes=[2, 2, 10, 2],
                                             label='Ventricle Count')
            # line1.set_dashes([2, 2, 10, 2])  # 2pt line, 2pt break, 10pt line, 2pt break
            line2, = axs[preIx, postIx].plot(t, criftTransmitterCounts[preIx, postIx], dashes=[6, 2],
                                             label='Clift Transmitter Count')
            line3, = axs[preIx, postIx].plot(t, freeRezeptorCounts[preIx, postIx], dashes=[6, 2, 4, 5],
                                             label='Free Rezeptor Count')
   #         line4, = axs[preIx, postIx].plot(t, acceptedTransmitterCounts[preIx, postIx], dashes=[8, 2, 2, 5],
    #                                         label='Accepted Transmitters Count')
            if preIx == 0 and postIx == 0:
                axs[preIx, postIx].legend()

            if postIx == noPost:
                axs[preIx, postIx].set_xlabel('time')
   # plt.show()

def visualizeNetwork(networkHistory, visualizeSynapses = True, visualizeNeurons = True):

    plt.style.use('bmh')
    #plt.style.use('Solarize_Light2')
    #plt.style.use('fivethirtyeight')
    #plt.style.use('dark_background')

    if visualizeSynapses:
        excitatoryMatrixHistory = np.full((networkHistory[0].numLayers, len(networkHistory)), L.SynapseMatrix)
        inhibitoryMatrixHistory = np.full((networkHistory[0].numLayers, len(networkHistory)), L.SynapseMatrix)
        lateralInhibitoryMatrixHistory = np.full((networkHistory[0].numLayers, len(networkHistory)), L.SynapseMatrix)

    if visualizeNeurons:
        allNeuronsHistory = []

    for time in range(len(networkHistory)):
        network = networkHistory[time]
        if visualizeNeurons:
            allNeuronsHistory.append(network.logNeurons())
        if visualizeSynapses:
            iterLayers = iter(network.layers)
            next(iterLayers)  # skipping the first layer
            for layerId, layer in enumerate(iterLayers):
                # TODO check which kind of layer is to be visualized
                excitatoryMatrixHistory[layerId,time] = layer.excitatoryMatrix
                lateralInhibitoryMatrixHistory[layerId,time] = layer.lateralInhibitoryMatrix

                if (isinstance(layer, L.DuoLateralInhibitoryLayer)):
                    inhibitoryMatrixHistory[layerId, time] = layer.inhibitoryMatrix

    exampleNetwork = networkHistory[0]
    if visualizeNeurons:
        allNeuronsHistory = np.array(allNeuronsHistory)
        visualizeNeuronPotentials(exampleNetwork, allNeuronsHistory)

    if visualizeSynapses:
        iterLayers = iter(exampleNetwork.layers)
        next(iterLayers)  # skipping the first layer
        for layerId, layer in enumerate(iterLayers):
            visualizeSynapseMatrix(excitatoryMatrixHistory[layerId], 'Layer ' + str(layerId) + ': Excitatory Synapses')
            visualizeSynapseMatrix(lateralInhibitoryMatrixHistory[layerId], 'Layer ' + str(layerId) + ': lateral Inhibitory Synapses')

            if (isinstance(layer, L.DuoLateralInhibitoryLayer)):
                visualizeSynapseMatrix(inhibitoryMatrixHistory[layerId],
                                       'Layer ' + str(layerId) + ': Inhibitory Synapses')

    plt.show()

# ----------------------------- STDP - visualisation ---------------------------

def visualiseLearningSignalHistory(networkHistory, calculations, timePerCalculation):
    exampleNetwork = networkHistory[0]
    exitoryLearningSignalsHistory = np.empty((exampleNetwork.numLayers, calculations, 3))
    inhibitoryLearningSignalsHistory = np.empty((exampleNetwork.numLayers, calculations, 3))

    excitatorySubsA = np.empty((exampleNetwork.numLayers, calculations, 1))
    excitatorySubsB = np.empty((exampleNetwork.numLayers, calculations, 1))
    excitatorySubsC = np.empty((exampleNetwork.numLayers, calculations, 1))
    effectiveExcitatoryAdjustment = np.empty((exampleNetwork.numLayers, calculations, 1))

    inhibSubsA = np.empty((exampleNetwork.numLayers, calculations, 1))
    inhibSubsB = np.empty((exampleNetwork.numLayers, calculations, 1))
    inhibSubsC = np.empty((exampleNetwork.numLayers, calculations, 1))
    effectiveInhibitoryAdjustment = np.empty((exampleNetwork.numLayers, calculations, 1))

    for c in range(0, calculations):
        # prepare visualisation of the learning signal
        cal = c + 1
        network = networkHistory[cal * (timePerCalculation - 1)]

        iterLayers = iter(network.layers)
        next(iterLayers)  # skipping the first layer
        for layerId, layer in enumerate(iterLayers):
            excitatorySynapseMatrix = layer.excitatoryMatrix.synapseMatrix
            excitatorySynapse = excitatorySynapseMatrix[0, 0].synapse
            exitoryLearningSignalsHistory[layerId, c] = excitatorySynapse[5]

            inhibitorySynapseMatrix = layer.inhibitoryMatrix.synapseMatrix
            inhibitorySynapse = inhibitorySynapseMatrix[0, 0].synapse
            inhibitoryLearningSignalsHistory[layerId, c] = inhibitorySynapse[5]

            if c == 0:
                excitatorySubsA[layerId, c] = exitoryLearningSignalsHistory[layerId, c, 0]
                excitatorySubsB[layerId, c] = exitoryLearningSignalsHistory[layerId, c, 1]
                excitatorySubsC[layerId, c] = exitoryLearningSignalsHistory[layerId, c, 2]

                inhibSubsA[layerId, c] = inhibitoryLearningSignalsHistory[layerId, c, 0]
                inhibSubsB[layerId, c] = inhibitoryLearningSignalsHistory[layerId, c, 1]
                inhibSubsC[layerId, c] = inhibitoryLearningSignalsHistory[layerId, c, 2]

            if c > 0:
                excitatorySubsA[layerId, c] = \
                    exitoryLearningSignalsHistory[layerId, c, 0] - exitoryLearningSignalsHistory[layerId, c-1, 0]
                excitatorySubsB[layerId, c] = \
                    exitoryLearningSignalsHistory[layerId, c, 1] - exitoryLearningSignalsHistory[layerId, c-1, 1]
                excitatorySubsC[layerId, c] = \
                    exitoryLearningSignalsHistory[layerId, c, 2] - exitoryLearningSignalsHistory[layerId, c-1, 2]

                inhibSubsA[layerId, c] = \
                    inhibitoryLearningSignalsHistory[layerId, c, 0] - inhibitoryLearningSignalsHistory[layerId, c-1, 0]
                inhibSubsB[layerId, c] = \
                    inhibitoryLearningSignalsHistory[layerId, c, 1] - inhibitoryLearningSignalsHistory[layerId, c-1, 1]
                inhibSubsC[layerId, c] = \
                    inhibitoryLearningSignalsHistory[layerId, c, 2] - inhibitoryLearningSignalsHistory[layerId, c-1, 2]

    effectiveExcitatoryAdjustment[layerId] = excitatorySubsC[layerId] - excitatorySubsB[layerId]
    effectiveInhibitoryAdjustment[layerId] = inhibSubsB[layerId] - inhibSubsC[layerId]

    iterLayers = iter(network.layers)
    next(iterLayers)  # skipping the first layer
    for layerId, layer in enumerate(iterLayers):
        cs = np.linspace((int(-calculations / 2)), int(calculations / 2), calculations)
        fig, axs = plt.subplots(3, 1, sharex='col', figsize=(3, 3))

        # ----
        axs[0].set_title('Excitatory Learning Signals')
        line1, = axs[0].plot(cs, excitatorySubsA[layerId], dashes=[2, 2, 10, 2], label='Substance A')
        line2, = axs[0].plot(cs, excitatorySubsB[layerId], dashes=[6, 2], label='Substance B - Weaken')
        line3, = axs[0].plot(cs, excitatorySubsC[layerId], dashes=[6, 2, 4, 2], label='Substance C - Enforce')
        line4, = axs[0].plot(cs, effectiveExcitatoryAdjustment[layerId], dashes=[3, 1, 5, 2], label='Effective Change')

        axs[0].set_ylabel('Intensity of the change') # in the information transmission effect')

        # set(leg1,'Interpreter','latex');
        # set(leg1,'FontSize',17);
        axs[0].legend()

        axs[1].set_title('Inhibitory Learning Signals')
        line5, = axs[1].plot(cs, inhibSubsA[layerId], label='Inhibitory Substance A')
        line5.set_dashes([8, 4, 10, 4])  # 2pt line, 2pt break, 10pt line, 2pt break
        line6, = axs[1].plot(cs, inhibSubsB[layerId], dashes=[6, 2, 2, 2], label='Inhibitory Substance B - Enforce')
        line7, = axs[1].plot(cs, inhibSubsC[layerId], dashes=[6, 4, 6, 2], label='Inhibitory Substance C - Weaken')
        line8, = axs[1].plot(cs, effectiveInhibitoryAdjustment[layerId], dashes=[4, 1, 6, 1, 4, 1],
                             label='Inhibitory Effective Change')
        axs[1].set_ylabel('Intensity of the change')
        axs[1].legend()

        axs[2].set_title('Learning Signals')
        line9, = axs[2].plot(cs, effectiveInhibitoryAdjustment[layerId], dashes=[4, 1, 6, 1, 4, 1],
                             label='Inhibitory Effective Change')
        line10, = axs[2].plot(cs, effectiveExcitatoryAdjustment[layerId], dashes=[3, 1, 5, 2], label='Effective Change')

        axs[2].set_ylabel('Intensity of the change')
        axs[2].set_xlabel('Time difference $\Delta(Spike_{Preneuron}, Spike_{Postneuron})$')
        axs[2].legend()

    plt.show()

# ----------------------------- Iris - visualisation ---------------------------


text:
import numpy as np
from train_classifier import *
from ge_criterion_baseline import *
from utilities import saveToFile, runBaselineTests, getModelAccuracy, getWeakSignalAccuracy
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from experiments import breast_cancer_reader, cardio_reader, obs_network_reader

def train_weak_signals(data, weak_signal_data, num_weak_signal):
    """
    Trains different views of weak signals

    :param data: dictionary of training and test data
    :type data: dict
    :param weak_signal_data: data representing the different views for the weak signals
    :type: array
    :param num_weak_signal: number of weak_signals
    type: in
    :return: dictionary containing of models, probabilities and error bounds of weak signals
    :rtype: dict
    """

    train_data, train_labels = data['training_data']
    val_data, val_labels = data['validation_data']
    test_data, test_labels = data['test_data']

    n, d = train_data.shape

    weak_signal_train_data = weak_signal_data[0]
    weak_signal_val_data = weak_signal_data[1]
    weak_signal_test_data = weak_signal_data[2]

    weak_signals = []
    stats = np.zeros(num_weak_signal)
    w_sig_probabilities = []
    w_sig_test_accuracies = []
    weak_val_accuracy = []

    for i in range(num_weak_signal):
        # fit model
        model = LogisticRegression(solver = "lbfgs", max_iter= 1000)
        model.fit(weak_signal_train_data[i], train_labels)
        weak_signals.append(model)

        # evaluate probability of P(X=1)
        probability = model.predict_proba(weak_signal_val_data[i])[:, 1]
        score = val_labels * (1 - probability) + (1 - val_labels) * probability
        stats[i] = np.sum(score) / score.size
        w_sig_probabilities.append(probability)

        # evaluate accuracy for validation data
        weak_val_accuracy.append(accuracy_score(val_labels, np.round(probability)))

        # evaluate accuracy for test data
        test_predictions = model.predict(weak_signal_test_data[i])
        w_sig_test_accuracies.append(accuracy_score(test_labels, test_predictions))

    model = {}
    model['models'] = weak_signals
    model['probabilities'] = np.array(w_sig_probabilities)
    model['error_bounds'] = stats
    model['validation_accuracy'] = weak_val_accuracy
    model['test_accuracy'] = w_sig_test_accuracies

    return model

def run_experiment(data, weak_signal_data, num_weak_signal, constant_bound=False):
    """
    Runs experiment with the given dataset

    :param data: dictionary of validation and test data
    :type data: dict
    :param weak_signal_data: data representing the different views for the weak signals
    :type: array
    :param num_weak_signal: number of weak signals
    :type num_weak_signal: int
    """

    w_model = train_weak_signals(data, weak_signal_data, num_weak_signal)

    training_data = data['training_data'][0].T
    training_labels = data['training_data'][1]
    val_data, val_labels = data['validation_data']
    val_data = val_data.T
    test_data = data['test_data'][0].T
    test_labels = data['test_data'][1]

    num_features, num_data_points = training_data.shape

    weak_signal_ub = w_model['error_bounds']
    # weak_signal_ub = np.ones(w_model['error_bounds'].shape) * 0.3
    models = w_model['models']
    weak_signal_probabilities = w_model['probabilities']

    weights = np.zeros(num_features)

    print("Running tests...")
    if constant_bound:
        optimized_weights, y = train_all(val_data, weights, weak_signal_probabilities, np.zeros(weak_signal_ub.size) + 0.3, max_iter=10000)
    else:
       optimized_weights, y = train_all(val_data, weights, weak_signal_probabilities, weak_signal_ub, max_iter=10000)

    # calculate validation results
    learned_probabilities = probability(val_data, optimized_weights)
    validation_accuracy = getModelAccuracy(learned_probabilities, val_labels)

    # calculate test results
    learned_probabilities = probability(test_data, optimized_weights)
    test_accuracy = getModelAccuracy(learned_probabilities, test_labels)

    # calculate weak signal results
    weak_val_accuracy = w_model['validation_accuracy']
    weak_test_accuracy = w_model['test_accuracy']

    adversarial_model = {}
    adversarial_model['validation_accuracy'] = validation_accuracy
    adversarial_model['test_accuracy'] = test_accuracy

    weak_model = {}
    weak_model['num_weak_signal'] = num_weak_signal
    weak_model['validation_accuracy'] = weak_val_accuracy
    weak_model['test_accuracy'] = weak_test_accuracy

    print("")
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
    print("We trained %d learnable classifiers with %d weak signals" %(1, num_weak_signal))
    print("The accuracy of learned model on the validatiion data is", validation_accuracy)
    print("The accuracy of weak signal(s) on the validation data is", weak_val_accuracy)
    print("The accuracy of the model on the test data is", test_accuracy)
    print("The accuracy of weak signal(s) on the test data is", weak_test_accuracy)
    print("")
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")

    # calculate baseline
    print("Running tests on the baselines...")
    baselines = runBaselineTests(val_data, weak_signal_probabilities) #remove the transpose to enable it run
    b_validation_accuracy = getWeakSignalAccuracy(val_data, val_labels, baselines)
    b_test_accuracy = getWeakSignalAccuracy(test_data, test_labels, baselines)
    print("The accuracy of the baseline models on test data is", b_test_accuracy)
    print("")
    weak_model['baseline_val_accuracy'] = b_validation_accuracy
    weak_model['baseline_test_accuracy'] = b_test_accuracy
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")

    # calculate ge criteria
    print("Running tests on ge criteria...")
    model = ge_criterion_train(val_data.T, val_labels, weak_signal_probabilities, num_weak_signal)
    ge_validation_accuracy = accuracy_score(val_labels, np.round(probability(val_data, model)))
    ge_test_accuracy = accuracy_score(test_labels, np.round(probability(test_data, model)))
    print("The accuracy of ge criteria on validation data is", ge_validation_accuracy)
    print("The accuracy of ge criteria on test data is", ge_test_accuracy)
    weak_model['gecriteria_val_accuracy'] = ge_validation_accuracy
    weak_model['gecriteria_test_accuracy'] = ge_test_accuracy
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")

    return adversarial_model, weak_model

def bound_experiment(data, weak_signal_data, num_weak_signal, bound):
    """
    Runs experiment with the given dataset

    :param data: dictionary of validation and test data
    :type data: dict
    :param weak_signal_data: data representing the different views for the weak signals
    :type: array
    :param num_weak_signal: number of weak signals
    :type num_weak_signal: int
    :param bound: error bound of the weak signal
    :type bound: int
    :return: outputs from the bound experiments
    :rtype: dict
    """

    w_model = train_weak_signals(data, weak_signal_data, num_weak_signal)

    training_data = data['training_data'][0].T
    training_labels = data['training_data'][1]
    val_data, val_labels = data['validation_data']
    val_data = val_data.T
    test_data = data['test_data'][0].T
    test_labels = data['test_data'][1]

    num_features, num_data_points = training_data.shape

    weak_signal_ub = w_model['error_bounds']
    weak_signal_probabilities = w_model['probabilities']

    weights = np.zeros(num_features)

    print("Running tests...")

    optimized_weights, ineq_constraint = train_all(val_data, weights, weak_signal_probabilities, bound, max_iter=10000)

    # calculate test probabilities
    test_probabilities = probability(test_data, optimized_weights)
    # calculate error bound on test data
    test_probabilities = np.round(test_probabilities)
    error_bound = (test_probabilities.dot(1 - test_labels) + (1 - test_probabilities).dot(test_labels)) / test_labels.size
    test_accuracy = getModelAccuracy(test_probabilities, test_labels)

    print("")
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
    print("The error_bound of learned model on the weak signal is", weak_signal_ub)
    print("The test accuracy of the weak signal is", w_model['test_accuracy'])
    print("The error_bound of learned model on the test data is", error_bound[0])
    print("The accuracy of the model on the test data is", test_accuracy)

    output = {}
    output['weak_signal_ub'] = weak_signal_ub[0]
    output['weak_test_accuracy'] = w_model['test_accuracy'][0]
    output['error_bound'] = error_bound[0]
    output['test_accuracy'] = test_accuracy
    output['ineq_constraint'] = ineq_constraint[0]

    return output

def dependent_error_exp(data, weak_signal_data, num_weak_signal):
    """
    Runs experiment with the given dataset

    :param data: dictionary of validation and test data
    :type data: dict
    :param weak_signal_data: data representing the different views for the weak signals
    :type: array
    :param num_weak_signal: number of weak signals
    :type num_weak_signal: int
    :return: test accuracies of ALL and respective baselines
    :rtype: dict
    """

    w_model = train_weak_signals(data, weak_signal_data, num_weak_signal)

    training_data = data['training_data'][0].T
    training_labels = data['training_data'][1]
    val_data, val_labels = data['validation_data']
    val_data = val_data.T
    test_data = data['test_data'][0].T
    test_labels = data['test_data'][1]

    num_features, num_data_points = training_data.shape

    weak_signal_ub = w_model['error_bounds']
    weak_signal_probabilities = w_model['probabilities']
    weak_test_accuracy = w_model['test_accuracy']

    weights = np.zeros(num_features)

    print("Running tests...")

    optimized_weights, ineq_constraint = train_all(val_data, weights, weak_signal_probabilities, weak_signal_ub, max_iter=5000)

    # calculate test probabilities
    test_probabilities = probability(test_data, optimized_weights)
    # calculate test accuracy
    test_accuracy = getModelAccuracy(test_probabilities, test_labels)

    print("")
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
    print("Experiment %d"%num_weak_signal)
    print("We trained %d learnable classifiers with %d weak signals" %(1, num_weak_signal))
    print("The accuracy of the model on the test data is", test_accuracy)
    print("The accuracy of weak signal(s) on the test data is", weak_test_accuracy)
    print("")

    # calculate ge criteria
    print("Running tests on ge criteria...")
    model = ge_criterion_train(val_data.T, val_labels, weak_signal_probabilities, num_weak_signal)
    ge_test_accuracy = accuracy_score(test_labels, np.round(probability(test_data, model)))
    print("The accuracy of ge criteria on test data is", ge_test_accuracy)
    print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")

    # calculate baseline
    print("Running tests on the baselines...")
    baselines = runBaselineTests(val_data, weak_signal_probabilities)
    b_test_accuracy = getWeakSignalAccuracy(test_data, test_labels, baselines)
    print("The accuracy of the baseline models on test data is", b_test_accuracy)
    print("")

    output = {}
    output['ALL'] = test_accuracy
    output['WS'] = w_model['test_accuracy'][-1]
    output['GE'] = ge_test_accuracy
    output['AVG'] = b_test_accuracy[-1]

    return output

def run_tests():
    """
    Runs experiment.
    :return: None
    """

    print("Running breast cancer experiment...")
    breast_cancer_reader.run_experiment(run_experiment, saveToFile)
    print("Running obs network experiment...")
    obs_network_reader.run_experiment(run_experiment, saveToFile)
    print("Running cardio experiment...")
    cardio_reader.run_experiment(run_experiment, saveToFile)

    # # un-comment to run bounds experimrnt in the paper
    #breast_cancer_reader.run_bounds_experiment(bound_experiment)
    #obs_network_reader.run_bounds_experiment(bound_experiment)

    # # un-comment to run dependency error experiment in the paper
    #print("Running dependent error on cardio experiment...")
    # cardio_reader.run_dep_error_exp(dependent_error_exp)

if __name__ == '__main__':
    run_tests()

ext: py

max_stars_repo_name: VTCSML/Adversarial-Label-Learning

max_issues_repo_name: VTCSML/Adversarial-Label-Learning

max_forks_repo_name: VTCSML/Adversarial-Label-Learning

reason: import numpy


text:
#!/usr/bin/env python
import torch
import math
import torch.nn as nn
import torch.nn.functional as F
import numpy as np

import aux_funcs as af
import model_funcs as mf

from architectures.CNNs.ResNet import ResNet

import pdb

class PolicyNet(nn.Module):
    """docstring for PolicyNet"""
    def __init__(self, data_type, num_classes=200):
        super(PolicyNet, self).__init__()
        params = {}
        params['num_blocks'] = [5,5,5]
        params['num_classes'] = num_classes
        params['augment_training'] = False
        if data_type == 'tiny_imagenet':
            params['input_size'] = 64
        else:
            params['input_size'] = 32
        params['block_type'] = 'basic'

        self.resnet = ResNet(params)

        self.end_layer = nn.Sequential(
            nn.Linear(num_classes, num_classes),
            nn.ReLU(),
            nn.Linear(num_classes, 1),
            nn.Sigmoid()
            )

    def forward(self, x0, yi):
        # x0_embed = self.resnet(x0)
        # merged = torch.cat((x0_embed, yi), 1)
        return self.end_layer(yi)

class SeqNet(nn.Module):
    def __init__(self, A, args, train_post, device):
        super(SeqNet, self).__init__()
        self.A = A
        self.m, self.n = A.shape
        self.args = args
        hidden_dims = args.policy_hidden_dims + '-' + str(1) # output size is one (score)
        self.mlp = MLP(input_dim=2*self.m + self.n + args.post_dim,
                       hidden_dims=hidden_dims,
                       nonlinearity=args.nonlinearity,
                       act_last=None)
        self.post_dim = args.post_dim
        self.train_post = train_post
        self.device = device

    def forward(self, y, xhs):
        pi_t_score = []
        position_enc = position_encoding(len(self.train_post), self.post_dim).to(self.device)
        batch_size = y.shape[0]
        for i, t in self.train_post.items():
            x_hat = xhs[t]
            ax = torch.matmul(x_hat, self.A.t())
            pe = torch.stack([position_enc[i] for _ in range(batch_size)], dim=0)
            f = torch.cat([y, y-ax, x_hat, pe], dim=-1)
            pi_t_score.append(self.mlp(f))

        scores = torch.cat(pi_t_score, dim=-1)
        # p = torch.sigmoid(self.mlp(f).view(-1))
        return scores

class MulticlassNet(nn.Module):
    def __init__(self, args, A, nz_post):
        super(MulticlassNet, self).__init__()
        self.A = A
        self.m, self.n = A.shape
        self.args = args

        self.nz_post = nz_post
        self.num_class = len(nz_post)
        multiclass_dims = args.policy_multiclass_dims + '-' + str(self.num_class)
        hidden_dims = tuple(map(int, args.policy_hidden_dims.split("-")))
        input_dim = hidden_dims[-1] * self.num_class

        self.mlp = MLP(input_dim=2*self.m + self.n,
                       hidden_dims=args.policy_hidden_dims,
                       nonlinearity=args.nonlinearity,
                       act_last=args.nonlinearity)

        self.mlp2 = MLP(input_dim=input_dim,
                       hidden_dims=multiclass_dims,
                       nonlinearity=args.nonlinearity,
                       act_last=None)

    def forward(self, y, xhs):
        class_feature = []
        for t in self.nz_post.values():
            x_hat = xhs[t]
            ax = torch.matmul(x_hat, self.A.t())
            inputs = torch.cat([y, y-ax, x_hat], dim=-1)
            f = self.mlp(inputs)
            class_feature.append(f)
        features = torch.cat(class_feature, dim=-1)
        scores = self.mlp2(features)
        # p = F.softmax(score, dim=-1)
        return scores

class MulticlassNetImage(nn.Module):
    """docstring for MulticlassNetImage"""
    def __init__(self, args, x, internal_fm, nz_post, category=200, pd=20, net_size=1,
        device='cuda', target_channel=128, share=False):
        super(MulticlassNetImage, self).__init__()
        self.args = args
        self.internal_fm = internal_fm
        self.nz_post = nz_post
        self.num_class = len(nz_post)
        self.target_channel = target_channel
        self.device = device
        self.pe = af.position_encoding(self.num_class, pd).to(self.device)
        self.input_shape_list = self._get_input_shapes()
        self.used_shape_list = self._get_used_shapes()
        self.smallest_fm_shape = self._get_smallest_shape()
        self.fm_module_list = nn.ModuleList()
        self.flatten_size = target_channel*self.smallest_fm_shape*self.smallest_fm_shape

        self.x_module = nn.Sequential(
                nn.Conv2d(in_channels=3,
                    out_channels=target_channel,
                    kernel_size=3, padding=1),
                nn.BatchNorm2d(target_channel),
                nn.ReLU(),
                nn.AdaptiveAvgPool2d(self.smallest_fm_shape),
                af.Flatten(),
                nn.Linear(self.flatten_size, self.flatten_size)
                )

        for shape in self.used_shape_list:
            if len(shape)==4:
                fm_module = nn.Sequential(
                    nn.Conv2d(in_channels=shape[1],
                        out_channels=target_channel,
                        kernel_size=3, padding=1),
                    nn.BatchNorm2d(target_channel),
                    nn.ReLU(),
                    nn.AdaptiveAvgPool2d(self.smallest_fm_shape),
                    af.Flatten(),
                    nn.Linear(self.flatten_size, self.flatten_size)
                    )
            if len(shape)==2:
                fm_module = nn.Sequential(
                    nn.Linear(shape[1], self.flatten_size),
                    nn.ReLU()
                    )
            self.fm_module_list.append(fm_module)

        self.total_channels = target_channel*(len(self.used_shape_list)+1)

        self.final_module = nn.Sequential(
            nn.Conv2d(in_channels=self.total_channels,
                out_channels=target_channel,
                kernel_size=3, padding=1),
            nn.BatchNorm2d(target_channel),
            nn.ReLU(),
            nn.Dropout(0.8),
            af.Flatten(),
            nn.Linear(self.flatten_size, self.num_class)
            )

        self.xhs_module_fc = nn.Sequential(
            nn.Linear(category*self.num_class, category),
            nn.ReLU(),
            nn.Dropout(0.7),
            nn.Linear(category, self.num_class)
            )

        self.encoder_layer = nn.TransformerEncoderLayer(category, 4)
        self.transformer_encoder = nn.TransformerEncoder(self.encoder_layer, 2)

        self.xhs_module_conv = nn.Sequential(
            nn.Conv1d(in_channels=category,
                out_channels=50,
                kernel_size=1),
            nn.ReLU(),
            nn.Dropout(0.7),
            af.Flatten(),
            nn.Linear(self.num_class*50, self.num_class)
            )

        self.shared_module = nn.Sequential(
            nn.Linear(category+pd+3+self.flatten_size, 400),
            nn.ReLU(),
            nn.Dropout(0.8),
            nn.Linear(400, 200),
            nn.ReLU(),
            nn.Dropout(0.8)
            )

        self.shared_module_fm = nn.Sequential(
            nn.Linear(self.flatten_size*2+3, 400),
            nn.ReLU(),
            nn.Dropout(0.8),
            nn.Linear(400, 200),
            nn.ReLU(),
            nn.Dropout(0.8)
            )

        self.last_layer = nn.Sequential(
            nn.Linear(200*self.num_class, self.num_class)
            )

    def _get_input_shapes(self):
        input_shape_list = []
        for fm in self.internal_fm:
            input_shape_list.append(fm.shape)
        return input_shape_list

    def _get_used_shapes(self):
        used_shape_list = []
        for i in range(self.num_class):
        # for t in self.nz_post.values():
            t = self.nz_post[i]
            used_shape_list.append(self.internal_fm[t].shape)
        return used_shape_list

    def _get_smallest_shape(self):
        len_list = list(map(lambda x: x[-1], self.used_shape_list))
        return min(len_list)

    def forward(self, x, fm, xhs):
        '''
        xhs: a list of 200dims, t*200 channels
        '''
        # x_map = self.x_module(x)
        # fm_map = list()
        # for i in range(self.num_class):
        #     fm_map.append(self.fm_module_list[i](fm[self.nz_post[i]]))
        # maps = torch.cat([x_map]+fm_map, dim=1)
        # return self.final_module(maps)

        batch_size = x.size(0)
        # pdb.set_trace()

        processed_x = self.x_module(x)
        # processed_x = processed_x.view(batch_size, -1) # batch*flatten_size

        xhs_used = list()
        for i in range(self.num_class):
            # pe = torch.stack([self.pe[i] for _ in range(batch_size)], dim=0)

            fm_processed = self.fm_module_list[i](fm[self.nz_post[i]])
            mme = af.mean_max_entropy(xhs[self.nz_post[i]])
            # go through a shared parameter module
            xhs_used.append(
                self.shared_module_fm(torch.cat([fm_processed, mme, processed_x], dim=-1)))

        # go through the final classification layer

        return self.last_layer(torch.cat(xhs_used, dim=-1))

        # xhs_used = xhs_used.permute(1, 0, -1) # t*batch*200
        # xhs_trans = self.transformer_encoder(xhs_used)
        # xhs_trans = xhs_trans.permute(1, -1, 0) # batch*200*t
        # return self.xhs_module_conv(xhs_trans)

class MNIconfidence(nn.Module):
    """docstring for MNIconfidence"""
    def __init__(self, args, x, internal_fm, nz_post, net_size=1,
        category=200, pd=20, device='cuda', target_channel=128, share=False):
        super(MNIconfidence, self).__init__()
        self.share = share
        # print('Share parameter: ', self.share)
        self.net_size = net_size
        self.args = args
        self.internal_fm = internal_fm
        self.nz_post = nz_post
        self.num_class = len(nz_post)
        self.target_channel = target_channel
        self.device = device
        self.pe = af.position_encoding(self.num_class, pd).to(self.device)
        self.input_shape_list = self._get_input_shapes()
        self.used_shape_list = self._get_used_shapes()
        self.smallest_fm_shape = self._get_smallest_shape()
        self.flatten_size = target_channel*self.smallest_fm_shape*self.smallest_fm_shape

        self.x_module = nn.Sequential(
                nn.Conv2d(in_channels=3,
                    out_channels=target_channel,
                    kernel_size=3, padding=1),
                nn.BatchNorm2d(target_channel),
                nn.ReLU(),
                nn.AdaptiveAvgPool2d(self.smallest_fm_shape),
                af.Flatten(),
                nn.Linear(self.flatten_size, self.flatten_size)
                )

        self.unshared_module_list = nn.ModuleList()

        for i in range(self.num_class):
            confidence_module = nn.Sequential(
                nn.Linear(category+pd+3, 200*net_size),
                nn.ReLU(),
                nn.Dropout(0.7),
                nn.Linear(200*net_size, 100*net_size),
                nn.ReLU(),
                nn.Dropout(0.7)
                )
            self.unshared_module_list.append(confidence_module)

        self.shared_module = nn.Sequential(
            nn.Linear(category+pd+3, 200*net_size),
            nn.ReLU(),
            nn.Dropout(0.8),
            nn.Linear(200*net_size, 100*net_size),
            nn.ReLU(),
            nn.Dropout(0.8)
            )
        self.last_layer = nn.Sequential(
            nn.Linear(100*self.num_class*net_size, self.num_class)
            )

    def _get_input_shapes(self):
        input_shape_list = []
        for fm in self.internal_fm:
            input_shape_list.append(fm.shape)
        return input_shape_list

    def _get_used_shapes(self):
        used_shape_list = []
        for i in range(self.num_class):
        # for t in self.nz_post.values():
            t = self.nz_post[i]
            used_shape_list.append(self.internal_fm[t].shape)
        return used_shape_list

    def _get_smallest_shape(self):
        len_list = list(map(lambda x: x[-1], self.used_shape_list))
        return min(len_list)

    def forward(self, x, fm, xhs):
        batch_size = x.size(0)

        processed_x = self.x_module(x)
        # processed_x = processed_x.view(batch_size, -1) # batch*flatten_size

        xhs_used = list()
        for i in range(self.num_class):
            pe = torch.stack([self.pe[i] for _ in range(batch_size)], dim=0)
            mme = af.mean_max_entropy(xhs[self.nz_post[i]])
            # go through a shared parameter module
            if self.share:
                xhs_used.append(
                    self.shared_module(torch.cat([xhs[self.nz_post[i]], mme, pe], dim=-1)))
            else:
                xhs_used.append(
                    self.unshared_module_list[i](torch.cat([xhs[self.nz_post[i]], mme, pe], dim=-1)))
        # go through the final classification layer

        return self.last_layer(torch.cat(xhs_used, dim=-1))

class Imiconfidence(nn.Module):
    """docstring for MNIconfidence"""
    def __init__(self, args, x, internal_fm, nz_post, net_size=1,
        category=200, pd=20, device='cuda', target_channel=128, share=False):
        super(Imiconfidence, self).__init__()
        self.share = share
        # print('Share parameter: ', self.share)
        self.net_size = net_size
        self.args = args
        self.internal_fm = internal_fm
        self.nz_post = nz_post
        self.num_class = len(nz_post)
        self.target_channel = target_channel
        self.device = device
        self.pe = af.position_encoding(self.num_class, pd).to(self.device)
        self.input_shape_list = self._get_input_shapes()
        self.used_shape_list = self._get_used_shapes()
        self.smallest_fm_shape = self._get_smallest_shape()
        self.flatten_size = target_channel*self.smallest_fm_shape*self.smallest_fm_shape

        self.unshared_module_list = nn.ModuleList()

        for i in range(self.num_class):
            confidence_module = nn.Sequential(
                nn.Linear(category+3, 1),
                # nn.Dropout(0.7)
                )
            self.unshared_module_list.append(confidence_module)

        self.shared_module = nn.Sequential(
            nn.Linear(category+3, 1),
            # nn.Dropout(0.7)
            )

    def _get_input_shapes(self):
        input_shape_list = []
        for fm in self.internal_fm:
            input_shape_list.append(fm.shape)
        return input_shape_list

    def _get_used_shapes(self):
        used_shape_list = []
        for i in range(self.num_class):
        # for t in self.nz_post.values():
            t = self.nz_post[i]
            used_shape_list.append(self.internal_fm[t].shape)
        return used_shape_list

    def _get_smallest_shape(self):
        len_list = list(map(lambda x: x[-1], self.used_shape_list))
        return min(len_list)

    def forward(self, x, fm, xhs):
        batch_size = x.size(0)

        xhs_used = list()
        for i in range(self.num_class):
            pe = torch.stack([self.pe[i] for _ in range(batch_size)], dim=0)

            mme = af.mean_max_entropy(xhs[self.nz_post[i]])
            # go through a shared parameter module
            if self.share:
                xhs_used.append(
                    self.shared_module(torch.cat([xhs[self.nz_post[i]], mme], dim=-1)))
            else:
                xhs_used.append(
                    self.unshared_module_list[i](torch.cat([xhs[self.nz_post[i]], mme], dim=-1)))
        # go through the final classification layer
        return torch.cat(xhs_used, dim=-1)

# model = MulticlassNetImage(cmd_args,x,  internal_fm, nz_post)
# model.to(device)
# tmp = model(x, internal_fm)


text:
"""
Simulation for problem 01 of SP1.
"""
from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import factorial

def R(N: int) -> Tuple[np.ndarray, np.ndarray]:
    # Generate the x array
    x = np.arange(0, N + 1, dtype=np.float64)
    # Probability
    p = 0.5
    # Calculate R
    y = factorial(N - 1) * np.power(p, N - 1) / (factorial(x) * factorial(N - x))

    return x, y

if __name__ == '__main__':
    for N in [6, 12, 24, 170]:
        x, y = R(N)
        plt.plot(x, y, label=f'N={N}')

    plt.legend()
    plt.xlabel(r'$0 \leq x \leq N$')
    plt.ylabel(r'$R(x, N)$')
    plt.savefig('export/problem01.svg')
    plt.show()

reason: import numpy,from scipy


text:
"""
@brief Calculate the gain and noise of a CCD camera system by
examining two flat field images. The calculation is the standard
mean/variance thing.
"""
from __future__ import print_function
from __future__ import absolute_import
import numpy as np
import numpy.random as random
import lsst.eotest.image_utils as imutils
import lsst.geom as lsstGeom
import lsst.afw.image as afwImage
import lsst.pex.exceptions as pexExcept

def flat_gain(image1, image2, count=1000, dx=100, dy=100, binsize=1,
              seed=None):
    """
    Calculate the gain and noise of a CCD camera system by
    examining two flat field images. The calculation is the standard
    mean/variance thing.
    """
    # If seed is None, the seed is generated from /dev/urandom.
    random.seed(seed)

    # Unbias using the mean bias of the two images.
    bmean = (imutils.bias(image1) + imutils.bias(image2))/2.
    image1 -= bmean
    image2 -= bmean

    # Trim prescan and overscan.
    image1 = imutils.trim(image1)
    image2 = imutils.trim(image2)

    # Rebin into binsize x binsize pixels.
    im1 = imutils.rebin(image1, binsize)
    im2 = imutils.rebin(image2, binsize)

    if dx > im1.getWidth():
        dx = im1.getWidth()
    if dy > im1.getHeight():
        dy = im1.getHeight()

    # Sample detector at size=count locations.
    try:
        xarr = random.randint(im1.getWidth() - dx - 1, size=count)
    except ValueError:
        # Rebinned image width smaller than requested dx, so just
        # generate subregions using the full x extent.
        xarr = np.zeros(count, dtype=np.int)
    try:
        yarr = random.randint(im1.getHeight() - dy - 1, size=count)
    except ValueError:
        # Rebinned image height smaller than requested dy, so just
        # generate subregions using the full y extent.
        yarr = np.zeros(count, dtype=np.int)

    gains = []
    ntrial = 0
    exception_count = 0
    for x, y in zip(xarr, yarr):
        bbox = lsstGeom.Box2I(lsstGeom.Point2I(int(x), int(y)),
                              lsstGeom.Extent2I(dx, dy))
        imarr1 = im1.Factory(im1, bbox).getArray()
        imarr2 = im2.Factory(im2, bbox).getArray()
        # Calculate flat ratio and correct subarray of image 2.
        fratio = np.mean(imarr1/imarr2)
        imarr2 *= fratio
        # Calculate the mean value of the flat field images.
        fmean = (np.mean(imarr1) + np.mean(imarr2))/2.
        # Calculate the variance of the flat difference image.
        fvar = np.var(imarr1 - imarr2)/2.
        gains.append(fvar/fmean)
    gain = 1./np.median(gains)  # gain in Ne/DN
    return gain, im1, im2

if __name__ == '__main__':
    import os
    from .sim_tools import simulateFlat

    file1 = 'test_flat1.fits'
    file2 = 'test_flat2.fits'

    hdus = 2
    simulateFlat(file1, 4000, 5, hdus=hdus)
    simulateFlat(file2, 4000, 5, hdus=hdus)

    for amp in range(1, hdus+1):
        image1 = afwImage.ImageF(file1, imutils.dm_hdu(amp))
        image2 = afwImage.ImageF(file2, imutils.dm_hdu(amp))
        gain, im1, im2 = flat_gain(image1, image2, count=1000)
        print(amp, gain)

    os.remove(file1)
    os.remove(file2)

max_stars_repo_name: tguillemLSST/eotest

max_issues_repo_name: tguillemLSST/eotest

max_forks_repo_name: tguillemLSST/eotest


text:
# Exercise 9.3
# Author: Noah Waterfield Price

import numpy as np

class Parabola:

    def __init__(self, c0, c1, c2):
        self.c0 = c0
        self.c1 = c1
        self.c2 = c2

    def __call__(self, x):
        return self.c2 * x ** 2 + self.c1 * x + self.c0

    def table(self, L, R, n):
        """Return a table with n points for L <= x <= R."""
        s = ''
        for x in np.linspace(L, R, n):
            y = self(x)
            s += '%12g %12g\n' % (x, y)
        return s

class SinPlusParabola(Parabola):

    def __init__(self, A, w, a, b, c):
        Parabola.__init__(self, a, b, c)
        self.A = A
        self.w = w

    def __call__(self, x):
        return Parabola.__call__(self, x) + self.A * np.sin(self.w * x)

sinPP = SinPlusParabola(10, 2, 3, 5, 2.5)
print sinPP.table(0, 10, 10)

"""
Sample run:
python sin_plus_quadratic.py
           0            3
     1.11111      19.5942
     2.22222      16.8136
     3.33333       51.186
     4.44444       79.711
     5.55556      98.0049
     6.66667      154.384
     7.77778      194.642
     8.88889      236.195
          10      312.129

"""


text:
import torch
import numpy as np
import os
from torch.utils.data import DataLoader, Dataset
from torchvision import transforms, utils
import pandas as pd
import random
import time
from sklearn import feature_selection, random_projection, svm
from sklearn.metrics import accuracy_score

class Data(Dataset):
    def __init__(self, filepath_n, filepath_t):
        n_x, n_y = Data.loader(filepath_n, 1)
        t_x, t_y = Data.loader(filepath_t, 0)

        self.x = np.concatenate((n_x, t_x))
        self.y = np.concatenate((n_y, t_y))
        # self.loader = Data.preprocess_img()

    def __getitem__(self, index):
        x = self.x[index]
        # x = self.loader(x)
        x = Data.preprocess(x)
        y = self.y[index]
        return x, y

    def __len__(self):
        return len(self.x)

    @staticmethod
    def preprocess(x):
        return torch.from_numpy(x)

    @staticmethod
    def preprocess_img(mean=0.0, std=1.0, normalize=False):
        if normalize:
            return transforms.Compose([
                transforms.ToTensor(),
                transforms.Normalize(mean=mean, std=std)
            ])
        else:
            return transforms.Compose([transforms.ToTensor()])

    @staticmethod
    def loader(filepath, label):
        data = pd.read_csv(filepath, sep="\t", index_col=0)
        data = pd.DataFrame(data.values.T, index=data.columns, columns=data.index)
        for col in list(data.columns[data.isnull().sum() > 0]):
            data[col].fillna(data[col].mean(), inplace=True)
            data[col] = (data[col] - data[col].mean()) / (data[col].std())
        data['label'] = label
        return data.drop('label', axis=1).values, data['label'].values

def floatTensor(cuda):
    return torch.cuda.FloatTensor if cuda else torch.FloatTensor

def longTensor(cuda):
    return torch.cuda.LongTensor if cuda else torch.LongTensor

def save(path):
    os.makedirs(path, exist_ok=True)

if __name__ == '__main__':
    data = Data('../RucPriv/BC-TCGA/BC-TCGA-Normal.txt', '../RucPriv/BC-TCGA/BC-TCGA-Tumor.txt')

max_stars_repo_name: HyperAndy/GANPriv

max_issues_repo_name: HyperAndy/GANPriv

max_forks_repo_name: HyperAndy/GANPriv


text:
__author__ = 'sibirrer'

import numpy as np
import pytest
import numpy.testing as npt
import lenstronomy.Util.param_util as param_util

try:
    import fastell4py
    fastell4py_bool = True
except:
    print("Warning: fastell4py not available, tests will be trivially fulfilled without giving the right answer!")
    fastell4py_bool = False

class TestSPEMD(object):
    """
    tests the Gaussian methods
    """
    def setup(self):
        from lenstronomy.LensModel.Profiles.spemd_smooth import SPEMD_SMOOTH
        self.SPEMD_SMOOT = SPEMD_SMOOTH()
        from lenstronomy.LensModel.Profiles.nie import NIE
        self.NIE = NIE()

    def test_function(self):
        phi_E = 1.
        gamma = 2.
        q = 0.999
        phi_G = 1.
        s_scale = 0.1
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        x = np.array([1., 2])
        y = np.array([2, 0])
        values = self.SPEMD_SMOOT.function(x, y, phi_E, gamma, e1, e2, s_scale)
        if fastell4py_bool:
            values_nie = self.NIE.function(x, y, phi_E, e1, e2, s_scale)
            delta_f = values[0] - values[1]
            delta_f_nie = values_nie[0] - values_nie[1]
            npt.assert_almost_equal(delta_f, delta_f_nie, decimal=5)
        else:
            npt.assert_almost_equal(values, 0, decimal=5)

    def test_derivatives(self):
        x = np.array([1])
        y = np.array([2])
        phi_E = 1.
        gamma = 2.
        q = 0.7
        phi_G = 1.
        s_scale = 0.1
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        f_x, f_y = self.SPEMD_SMOOT.derivatives(x, y, phi_E, gamma, e1, e2, s_scale)
        if fastell4py_bool:
            f_x_nie, f_y_nie = self.NIE.derivatives(x, y, phi_E, e1, e2, s_scale)
            #TODO test with higher precision, convention of theta_E might be slightly different from NIE with SPEMD
            npt.assert_almost_equal(f_x, f_x_nie, decimal=2)
            npt.assert_almost_equal(f_y, f_y_nie, decimal=2)
        else:
            npt.assert_almost_equal(f_x, 0, decimal=7)
            npt.assert_almost_equal(f_y, 0, decimal=7)

    def test_hessian(self):
        x = np.array([1.])
        y = np.array([2.])
        phi_E = 1.
        gamma = 2.
        q = 0.9
        phi_G = 1.
        s_scale = 0.1
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        f_xx, f_yy,f_xy = self.SPEMD_SMOOT.hessian(x, y, phi_E, gamma, e1, e2, s_scale)
        if fastell4py_bool:
            f_xx_nie, f_yy_nie, f_xy_nie = self.NIE.hessian(x, y, phi_E, e1, e2, s_scale)
            npt.assert_almost_equal(f_xx, f_xx_nie, decimal=3)
            npt.assert_almost_equal(f_yy, f_yy_nie, decimal=3)
            npt.assert_almost_equal(f_xy, f_xy_nie, decimal=3)
        else:
            npt.assert_almost_equal(f_xx, 0, decimal=7)
            npt.assert_almost_equal(f_yy, 0, decimal=7)
            npt.assert_almost_equal(f_xy, 0, decimal=7)

    def test_bounds(self):
        theta_E, gamma, q, phi_G, s_scale = self.SPEMD_SMOOT._parameter_constraints(theta_E=-1, s_scale=0, gamma=3, q=2, phi_G=0)
        assert theta_E == 0

    def test_is_not_empty(self):
        func = self.SPEMD_SMOOT.is_not_empty

        assert func(0.1, 0.2)
        assert func([0.1], [0.2])
        assert func((0.1, 0.3), (0.2, 0.4))
        assert func(np.array([0.1]), np.array([0.2]))
        assert not func([], [])
        assert not func(np.array([]), np.array([]))

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

max_stars_repo_name: tvelagapudi/lenstronomy

max_issues_repo_name: tvelagapudi/lenstronomy

max_forks_repo_name: tvelagapudi/lenstronomy


text:
"""
Loader functions for the datasets. You can add yours to here.

The function requires 2 inputs:
- name: Name of dataset as the name of the file without suffix
- path: Path to where the file is located

The functions return 2 outputs:
- dict_snapshots: Dictionary of snapshots, where keys are time stamps and values are list of edges occurring at this
                    time.
- dict_weights: Dictionary of weights according to dict_snapshots - keys are the time stamps and values are list of
                weight edges corresponding to the order of edges in dict_snapshots.

Note that if you change the name of the function "data_loader", you have to change it as well in the file "main.py" in
in line 51.
"""

from datetime import datetime
import os
import numpy as np

def data_loader(name, path):
    """
    Data loader for all datasets used.
    :param name: Name of the file (str, without .txt)
    :param path: Path to where it is
    :return:
    """
    if name == "dblp":
        dict_snapshots, dict_weights = load_dblp(name, path)
    else:
        dict_snapshots, dict_weights, times = {}, {}, []
        with open(os.path.join(path, name + ".txt"), 'r') as filehandle:
            for line in filehandle:
                a = line.split(" ")
                if len(a) == 1:
                    a = line.split(",")
                if "%" in a[0]:
                    continue
                else:
                    node1 = a[0]
                    node2 = a[1]
                    if node1 == node2:
                        continue
                    try:
                        w = float(a[2])
                        t = a[3].split("\n")[0]
                    except:
                        w = 1.
                        t = a[2].split("\n")[0]
                    if t == '0':
                        date = int(t)
                    else:
                        x = datetime.fromtimestamp(int(t))
                        times.append(x)
                        month = x.month
                        year = x.year
                        date = "{}.{}".format(month, year)
                    if dict_snapshots.get(date) is None:
                        dict_snapshots.update({date: [(node1, node2)]})
                    else:
                        dict_snapshots[date].append((node1, node2))
                    if dict_weights.get(date) is None:
                        dict_weights.update({date: [w]})
                    else:
                        dict_weights[date].append(w)
        times.sort()
        sorted_times = sort_keys(times)
        dict_snapshots = sort_dict_snapshots(dict_snapshots, sorted_times)
        if name == "facebook-wosn-wall":
            del dict_snapshots["10.2004"]
            del dict_snapshots["11.2004"]
            del dict_snapshots["12.2004"]
    return dict_snapshots, dict_weights

def load_dblp(name, path, is_weighted=False):
    """
    DBLP dataset
    """
    dict_snapshots = {}
    with open(os.path.join(path, name + ".txt"), 'r') as filehandle:
        for line in filehandle:
            a = line.split(" ")
            node1 = a[0]
            node2 = a[1]
            if is_weighted:
                w = float(a[2])
                t = a[3].split("\n")[0]
            else:
                w = 1.
                t = a[2].split("\n")[0]
            if dict_snapshots.get(t) is None:
                dict_snapshots.update({t: [(node1, node2, w)]})
            else:
                dict_snapshots[t].append((node1, node2, w))
    sorted_dict = dict(sorted(dict_snapshots.items(), key=lambda kv: kv[0]))
    keys = sorted_dict.keys()
    dict_weights = {key: [] for key in keys}
    new_dict = {key: [] for key in keys}
    for key in keys:
        for edge in sorted_dict[key]:
            dict_weights[key].append(edge[2])
            new_dict[key].append((edge[0], edge[1]))
    return new_dict, dict_weights

def sort_keys(times):
    """
    Function to sort the times
    """
    sorted_times = []
    for t in times:
        year = t.year
        month = t.month
        date = "{}.{}".format(month, year)
        if date in sorted_times:
            continue
        else:
            sorted_times.append(date)
    return sorted_times

def sort_dict_snapshots(dict_snapshots, times):
    """
    Sort dictionary of snapshots by time
    :param dict_snapshots:
    :param times:
    :return:
    """
    new_dict = {t: dict_snapshots[t] for t in times}
    return new_dict

def load_embedding(path, file_name):
    """
    Given a .npy file - embedding of a given graph. return the embedding dictionary
    :param path: Where this file is saved.
    :param file_name: The name of the file
    :return: Embedding dictionary
    """
    data = np.load(os.path.join(path, '{}.npy'.format(file_name)), allow_pickle=True)
    dict_embedding = data.item()
    return dict_embedding


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

from __future__ import division
import sys
from optparse import OptionParser
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import blackmanharris
from scipy.fftpack import fft
from common import load, analyze, analyze_function, rms_flat, dB, rms_fft, rms_rfft, setup_graph, takeClosest
#from thd import THDN, THD

import traceback

import os
local_dir = os.path.dirname(os.path.realpath(__file__))
'''
# print current info
filename, codeline, funcName, text = traceback.extract_stack()[-1]
print (filename, codeline, funcName, text)
print traceback.extract_stack(None, 1)[0][2] # print current function name

# print caller info
filename, codeline, funcName, text = traceback.extract_stack()[-2]
print (filename, codeline, funcName, text)
print traceback.extract_stack(None, 2)[0][2] # print current function name
'''

'''
def find_range(f, x):
    """
    Find range between nearest local minima from peak at index x
    """
    for i in np.arange(x+1, len(f)):
        if f[i+1] >= f[i]:
            uppermin = i
            break
    for i in np.arange(x-1, 0, -1):
        if f[i] <= f[i-1]:
            lowermin = i + 1
            break
    return (lowermin, uppermin)

def THDN(signal, sample_rate):
    """
    Measure the THD+N for a signal and print the results

    Prints the estimated fundamental frequency and the measured THD+N.  This is
    calculated from the ratio of the entire signal before and after
    notch-filtering.

    Currently this tries to find the "skirt" around the fundamental and notch
    out the entire thing.  A fixed-width filter would probably be just as good,
    if not better.
    """
    # Get rid of DC and window the signal
    signal -= np.mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    windowed = signal * blackmanharris(len(signal))  # TODO Kaiser?

    # Measure the total signal before filtering but after windowing
    total_rms = rms_flat(windowed)
#    print '\ntotal_rms: %f' % (total_rms)
#    print "\n RMS = %f" % ( np.sqrt(np.mean(windowed**2)) )

    # Find the peak of the frequency spectrum (fundamental frequency), and
    # filter the signal by throwing away values between the nearest local
    # minima
    f = np.fft.rfft(windowed)
    i = np.argmax(abs(f))
#    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed)))  # Not exact
    lowermin, uppermin = find_range(abs(f), i)
    f[lowermin: uppermin] = 0

    # Transform noise back into the signal domain and measure it
    # TODO: Could probably calculate the RMS directly in the frequency domain instead
    noise = np.fft.irfft(f)
    THDN = rms_flat(noise) / total_rms
    print "THD+N:     %.4f%% or %.1f dB" % (THDN * 100, 20 * np.log10(THDN))

def THD(signal, sample_rate):
    """Measure the THD for a signal

    This function is not yet trustworthy.

    Returns the estimated fundamental frequency and the measured THD,
    calculated by finding peaks in the spectrum.

    There are two definitions for THD, a power ratio or an amplitude ratio
    When finished, this will list both

    """
    # Get rid of DC and window the signal
    signal -= np.mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    windowed = signal * blackmanharris(len(signal))  # TODO Kaiser?

    # Find the peak of the frequency spectrum (fundamental frequency), and
    # filter the signal by throwing away values between the nearest local
    # minima
    f = np.fft.rfft(windowed)
    i = np.argmax(abs(f))
#    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed)))  # Not exact
#    print 'fundamental amplitude: %.3f' % abs(f[i])

    # Find the values for the first 15 harmonics.  Includes harmonic peaks
    # only, by definition
    # TODO: Should peak-find near each one, not just assume that fundamental
    # was perfectly estimated.
    # Instead of limited to 15, figure out how many fit based on f0 and
    # sampling rate and report this "4 harmonics" and list the strength of each
#    for x in range(2, 15):
#        print('%.3f' % abs(f[i * x]))

    THD = sum([abs(f[i*x]) for x in range(2, 15)]) / abs(f[i])
    print '\nTHD: %f%%' % (THD * 100)
'''

def SaveGraph(signal, sample_rate):
    # https://github.com/calebmadrigal/FourierTalkOSCON/blob/master/06_FFTInPython.ipynb
#    signal = signal / (2.**15) #convert sound array to float pt. values
    print 'sample_rate = ' + str(sample_rate) # samples per sec

    sample = len(signal) # euq signal.shape[0]
    print 'samples = ' + str(sample)

    total_sampling_time = sample/sample_rate
    print 'total_sampling_time = ' + str(total_sampling_time) + ' seconds'

    num_samples = sample_rate * total_sampling_time # euq signal.shape[0] and euq len(signal)
#    print 'num_samples (samples) = ' + str(num_samples)

    t = np.linspace(0, total_sampling_time, num_samples)
#    '''
    setup_graph(title='Signal output', x_label='time (in seconds)', y_label='amplitude', fig_size=(12,6))
    plt.plot(t, signal)
    plt.savefig('signal.png')
#    '''

    # Take the fft
    fft_output = np.fft.fft(signal)
#    '''
    setup_graph(title='FFT output', x_label='FFT bins', y_label='amplitude', fig_size=(12,6))
    plt.plot(fft_output)
    plt.savefig('fft.png')
#    '''

    # For real-valued input, the fft output is always symmetric.
    rfft_output = np.fft.rfft(signal)
#    '''
    setup_graph(title='rFFT output', x_label='frequency (in Hz)', y_label='amplitude', fig_size=(12,6))
    plt.plot(rfft_output)
    plt.savefig('rfft.png')
#   '''

    # Getting frequencies on x-axis labels
    # the x-axis to represent frequency
    # Frequencies range from 0 to the Nyquist Frequency (sample rate / 2)
    rfreqs = [(i*1.0/num_samples)*sample_rate for i in range(int(num_samples//2+1))]
#    '''
    setup_graph(title='Corrected rFFT Output', x_label='frequency (in Hz)', y_label='amplitude', fig_size=(12,6))
    plt.plot(rfreqs, rfft_output)
    plt.savefig('rfft_corrected.png')
#    '''

    # Getting negative values on y-axis labels
    # the y-axis to represent magnitude
#    '''
    rfft_mag = [np.sqrt(i.real**2 + i.imag**2)/len(rfft_output) for i in rfft_output]
    setup_graph(title='Corrected rFFT Output', x_label='frequency (in Hz)', y_label='magnitude', fig_size=(12,6))
    plt.plot(rfreqs, rfft_mag)
    plt.savefig('rfft_magnitude.png')
#    '''

#    '''
    # take the output of the FFT and perform an Inverse FFT to get back to our original wave (using the Inverse Real FFT - irfft).
    irfft_output = np.fft.irfft(rfft_output)
    setup_graph(title='Inverse rFFT Output', x_label='time (in seconds)', y_label='amplitude', fig_size=(12,6))
    plt.plot(t, irfft_output)
    plt.savefig('rfft_inverse.png')
#    '''
#    plt.show()

def DrawGraph(signal, sample_rate):
    # https://github.com/calebmadrigal/FourierTalkOSCON/blob/master/06_FFTInPython.ipynb
#    signal = signal / (2.**15) #convert sound array to float pt. values
    print 'sample_rate = ' + str(sample_rate) # samples per sec

    sample = len(signal) # euq signal.shape[0]
    print 'samples = ' + str(sample)

    total_sampling_time = sample/sample_rate
    print 'total_sampling_time = ' + str(total_sampling_time) + ' seconds'

    num_samples = sample_rate * total_sampling_time # euq signal.shape[0] and len(signal)
#    print 'num_samples (samples) = ' + str(num_samples)

#    '''
    time_array = np.arange(0, num_samples*1.0/sample_rate, 1.0/sample_rate) #start, stop, step
    setup_graph(title='Time domain Output', x_label='Time (in Sec)', y_label='amplitude', fig_size=(12,4))
    plt.plot(time_array, signal)
#    '''

    rfft_output = np.fft.rfft(signal)

    # Plotting the Frequency Content
    amplitude = [np.sqrt(i.real**2 + i.imag**2)/len(rfft_output) for i in rfft_output]
    frequencies = [(i*1.0/num_samples)*sample_rate for i in range(int(num_samples//2+1))]

    setup_graph(title='Frequency domain Output', x_label='frequency (in Hz)', y_label='amplitude', fig_size=(12,4))
    plt.plot(frequencies, amplitude, color='r')

    power = dB(amplitude)
    setup_graph(title='Frequency domain Output', x_label='frequency (in Hz)', y_label='power (dB)', fig_size=(12,4))
    plt.plot(frequencies, power, color='k')

    setup_graph(title='Spectrogram', x_label='time (in seconds)', y_label='frequency', fig_size=(14,7))
    plt.specgram(signal, Fs=sample_rate)

    plt.show()

def wave_analyzer(files):
    if files:
        for filename in files:
            try:
                load(filename)
            except IOError:
                print 'Couldn\'t analyze "' + filename + '"\n'
            print ''
    else:
        # TODO: realtime analyzer goes here
        sys.exit("You must provide at least one file to analyze:\n"
                 "python wave_analyzer.py filename.wav")

def CalculateBandRMS(signal, sample_rate, band=1, percent=5):
#    print sys._getframe().f_code.co_name # print function name

    # https://gist.github.com/endolith/1257010
    # https://scipy.github.io/devdocs/tutorial/fftpack.html
    num_samples = len(signal)

    rfft_array = np.fft.rfft(signal)
    frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)] #0,(1/n)*rate,(2/n)*rate,....((n/2)/n)*rate

    '''
    rms = np.mean(map(lambda x: x**2, signal)) ** .5
    print rms
    rms_val = np.sqrt(np.mean(signal**2))
    db_val = 20*np.log10(rms_flat(signal))
    print (rms_val, db_val)

    rms_val_raw = rms_flat(signal)
    print('Original RAW RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_raw, dB(rms_val_raw) ))

    fft_signal = np.fft.fft(signal)
    rms_val_fft = rms_fft( fft_signal )
    print('Original FFT RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_fft, dB(rms_val_fft) ))

    rms_val_fft2ifft = rms_flat( np.fft.ifft(fft_signal) )
    print('Original FFT to iFFT RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_fft2ifft, dB(rms_val_fft2ifft) ))

    rfft_signal = np.fft.rfft(signal)
    # Only approximate for odd n:
    rms_val_rfft2irfft = rms_flat( np.fft.irfft(rfft_signal) )
    print('Original rFFT to irFFT RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_rfft2irfft, dB(rms_val_rfft2irfft) ))

    rms_val_rfft = rms_rfft(rfft_signal)
    print('Original rFFT RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_rfft, dB(rms_val_rfft) ))

    # Accurate for odd n:
    rms_val_rfft2irfft_with_size = rms_flat( np.fft.irfft(rfft_signal, num_samples) )
    print('Original rFFT to irFFT RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_rfft2irfft_with_size, dB(rms_val_rfft2irfft_with_size) ))

    rms_val_rfft_with_size = rms_rfft( rfft_signal, num_samples )
    print('Original rFFT RMS: {:+.2f}, {:+.2f} dB'.format( rms_val_rfft_with_size, dB(rms_val_rfft_with_size) ))
    '''

#    irfft_array = np.fft.irfft(rfft_array) # irfft_array must euq signal
    bandwidth = band*1000
    lowcut = 0 # Hz # Remove lower frequencies.
    highcut = 0 # Hz # Remove higher frequencies.

    search_range = ( bandwidth * percent ) / 100

    lowcut = takeClosest(frequencies, (bandwidth-search_range) )
    highcut  = takeClosest(frequencies, (bandwidth+search_range) )
#    print(lowcut, highcut)

    cut_rfft_array = rfft_array.copy()
    cut_rfft_array[:lowcut] = 0 # Hz # Low cut
    cut_rfft_array[highcut:] = 0 # Hz # High cut

    cut_signal = np.fft.irfft(cut_rfft_array)
    rms_val_cut = rms_flat(cut_signal)
#    print('band=%d rms=%.4f power=%.4f' %( band, rms_val_cut, dB(rms_val_cut) ))
    print('SET BAND_%d=%.4f' %( band, dB(rms_val_cut) ))

    '''
    setup_graph(title='Spectrogram (Before)', x_label='time (in seconds)', y_label='frequency', fig_size=(14,7))
    plt.specgram(signal, Fs=sample_rate)

    setup_graph(title='Spectrogram (After)', x_label='time (in seconds)', y_label='frequency', fig_size=(14,7))
    plt.specgram(cut_signal, Fs=sample_rate)

    time_array = np.arange(0, num_samples*1.0/sample_rate, 1.0/sample_rate) #start, stop, step

    rms = rms_flat(signal)
    plt.figure(figsize=(10,5))
    plt.plot(time_array, signal, color="red")
    plt.plot([0, num_samples*1.0/sample_rate], [rms, rms], color="green", linewidth=2, alpha=.8)
    plt.plot([0, num_samples*1.0/sample_rate], [-rms, -rms], color="green", linewidth=2, alpha=.8)
#    plt.plot(signal, color="blue")
#    plt.plot([0, len(signal)], [rms, rms], color="green", linewidth=2, alpha=.8)
#    plt.plot([0, len(signal)], [-rms, -rms], color="green", linewidth=2, alpha=.8)
    plt.title( 'Time domain (Before)\nRMS=' + str(rms) )
    plt.xlabel('time (in seconds)')
    plt.ylabel('amplitude')

    rms = rms_flat(cut_signal)
    plt.figure(figsize=(10,5))
    plt.plot(time_array, cut_signal, color="red")
    plt.plot([0, num_samples*1.0/sample_rate], [rms, rms], color="green", linewidth=2, alpha=.8)
    plt.plot([0, num_samples*1.0/sample_rate], [-rms, -rms], color="green", linewidth=2, alpha=.8)
    plt.title( 'Time domain (After)\nRMS=' + str(rms) )
    plt.xlabel('time (in seconds)')
    plt.ylabel('amplitude')

    plt.show()
    '''

def CalculateBandMAX(signal, sample_rate, band_start=1, band_stop=1, percent=5, left_right="none", filename="none"):
#    print traceback.extract_stack(None, 1)[0][2] # print current function name

    '''
    import traceback
    stack = traceback.extract_stack()
    filename, codeline, funcName, text = stack[-2] # print caller function info
    print (filename, codeline, funcName, text)
    filename, codeline, funcName, text = stack[-1] # print current function info
    print (filename, codeline, funcName, text)
    filename, codeline, funcName, text = stack[0] # print root function info
    print (filename, codeline, funcName, text)
    '''
#    print traceback.extract_stack(None, 2)[0][2] # print caller function name
#    print traceback.extract_stack(None, 1)[0][2] # print current function name

    num_samples = len(signal)

    rfft_array = np.fft.rfft(signal)
    amplitude = [np.sqrt(i.real**2 + i.imag**2)/len(rfft_array) for i in rfft_array]
    frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)] #0,(1/n)*rate,(2/n)*rate,....((n/2)/n)*rate
    power = dB(amplitude)

    fft_signal = np.fft.fft(signal)
    freqs = np.fft.fftfreq(len(fft_signal))
    #print(freqs.min(), freqs.max())
    # (-0.5, 0.499975)

    # Find the peak in the coefficients
    idx = np.argmax(np.abs(fft_signal))
    freq = freqs[idx]
    freq_in_hertz = abs(freq * sample_rate)
    print(freq_in_hertz)

#    DrawGraph(signal, sample_rate)

    '''
    np.savetxt("amplitude.txt", amplitude, fmt='%6.3f')
    np.savetxt("frequencies.txt", frequencies, fmt='%6.3f')
    np.savetxt("power.txt", power, fmt='%6.3f')
    '''

#    '''
#    min_val,idx = min( (l[i],i) for i in xrange(len(l)) )
#    max_val,idx = max( (l[i],i) for i in xrange(len(l)) )

    peaks = []
#    print "{:*^80s}".format('Peak List')
    for x in range(band_start, band_stop+1): #check peak value from band_start to band_end kHz
        kx = 1000 * x
        search_range = (kx*percent)/100
        kx_start = kx - search_range
        kx_stop  = kx + search_range

        kx_idx_start = takeClosest(frequencies, kx_start)
        kx_idx_stop  = takeClosest(frequencies, kx_stop) + 1

        #find a good max amplitude.
#        maxAmplitude = np.max( np.absolute( amplitude[kx_idx_start:kx_idx_stop] ) )
#        peak_idx = closest(amplitude, maxAmplitude)
#        peak_idx = int( np.where( amplitude==np.max(amplitude[kx_idx_start:kx_idx_stop]) )[0] )
        peak_idx = int( np.where( power==np.max(power[kx_idx_start:kx_idx_stop]) )[0] )

        '''
        print 'search_range = ' + str(search_range)
        print 'kx = ' + str(kx) + ' [ ' + str(kx_start) + ' ~ ' + str(kx_stop) + ' ]'
        print 'kx_idx_start = ' + str(kx_idx_start)
        print 'kx_idx_stop  = ' + str(kx_idx_stop)

        print '--> peak_idx = ' + str(peak_idx)

        offset = (abs(frequencies[kx_idx_start] - frequencies[peak_idx])/frequencies[kx_idx_start]) * 100
        print 'Frequency: %.4f Hz amplitude: %.4f, offset of %.4f Hz: %.4%%' %(frequencies[peak_idx], amplitude[peak_idx], frequencies[kx_idx_start], offset)
        '''
#        print 'band=%d frequency=%.4f amplitude=%.4f power=%.4f' %(x, frequencies[peak_idx], amplitude[peak_idx], power[peak_idx] )
        # 1 kilohertz (kHz) to 1000 hertz (Hz) = 1 kHz × 1000 = 1000 Hz

        print 'SET %s=%.4f' %(left_right, power[peak_idx] )
        result(power[peak_idx], left_right, filename)
#    print "{:*^80s}".format('')

def result(audio_dB, left_right, filename):
    set_result = "Fail"

    result_file = filename.split(".")[0] + "_out.bat"
    if audio_dB > -10:
        set_result = "PASS"

    result_path = os.path.join(local_dir, result_file)
    result_string = "Set " + left_right + "dB=" + str(audio_dB) + "\n"
    result_string2 = "Set " + left_right + "=" + set_result + "\n"
    #print(result_string, result_path)

    f = open(result_path, 'a+')
    f.write(result_string)
    f.write(result_string2)
    f.closed

def usage():
    print("Usage:")
    print(sys.argv[0] + ' -w wavefile -f function -b begin -e end -p percent')
    sys.exit(2)

def parse_args():
    # http://python.usyiyi.cn/translate/python_278/library/argparse.html
    parser = OptionParser(usage='%prog -w wavefile -f function -b begin -e end -p percent')
    parser.add_option('-w', '--wavefile',
                        action="store", dest="wavefile",
                        help="Wave Filename", default="None")
    parser.add_option('-f', '--function',
                        action="store", dest="function",
                        help="Run Function", default="None")
    parser.add_option('-b', '--begin',
                        action="store", dest="begin",
                        help="Start Band", default=1, type=int)
    parser.add_option('-e', '--end',
                        action="store", dest="end",
                        help="Stop Band", default=1, type=int)
    parser.add_option('-p', '--percent',
                        action="store", dest="percent",
                        help="Search Percent", default=5, type=int)

    options, args = parser.parse_args()
#    print "Length : %d" % len(options.__dict__)
#    print options

    if not (options.wavefile and options.function):
        usage()

    return options

if __name__ == '__main__':
    '''
    options = parse_args()

    signal, sample_rate, channels = load(options.wavefile)

    if options.function.lower()=='calculatebandmax':
        CalculateBandMAX(signal, sample_rate, options.begin, options.end, options.percent)
    elif options.function.lower()=='calculatebandrms':
        for band in range(options.begin, options.end+1):
            CalculateBandRMS(signal, sample_rate, band, options.percent)
    else:
        print "Support Function:\n\tCalculateBandMAX\n\tCalculateBandRMS"

    '''
    try:
        import sys
        files = sys.argv[1:]
        if files:
            for filename in files:
                try:
                    if(os.path.isfile(filename)):
                        print("file exist")
                    else:
                        sys.exit("File not found")
#                    analyze(filename)
#                    PrintWavHeader(filename)
#                   analyze_function(filename, SaveGraph)
#                    analyze_function(filename, DrawGraph)
#                    analyze_function(filename, THDN)
                    analyze_function(filename, CalculateBandMAX)
#                    analyze_function(filename, CalculateBandRMS)

                except IOError:
                    print 'Couldn\'t analyze "' + filename + '"\n'
                print ''
        else:
            sys.exit("You must provide at least one file to analyze")
    except BaseException as e:
        print('Error:')
        print(e)
        raise
#    finally:
        # Otherwise Windows closes the window too quickly to read
#        raw_input('(Press <Enter> to close)')
#    '''

max_stars_repo_name: auspbro/sound-analysis

max_issues_repo_name: cs-joy/sound-analysis

max_forks_repo_name: cs-joy/sound-analysis