text:
import numpy as np
import re
from sys import argv
import sys
import csv
import glob
import os
import cv2
import operator
import matplotlib.pyplot as plt

"""
Convert a rgb image to grayscale.
"""
def rgb_to_gray(rgb):
    return np.dot(rgb[... , :3] , [0.299 , 0.587, 0.114])

"""
Order strings in natural order.
"""
def stringSplitByNumbers(x):
    r = re.compile('(\d+)')
    l = r.split(x)
    return [int(y) if y.isdigit() else y for y in l]

"""
Remove duplicated line from csv file saved by Biolapse interface.
"""
def remove_duplicate_csv(source):
	# # Create the output file as input with _deduped before .csv extension
	# source = argv[1]
	destination = source.replace('.csv', '_deduped.csv')
	data = open(source, 'r')
	target = open(destination, 'w')
	# Let the user know you are starting, in case you are de-dupping a huge file
	print("\nRemoving duplicates from %r" % source)

	# Initialize variables and counters
	unique_lines = set()
	source_lines = 0
	duplicate_lines = 0

	# Loop through data, write uniques to output file, skip duplicates.
	for line in data:
		source_lines += 1
		# Strip out the junk for an easy set check, also saves memory
		line_to_check = line.strip('\r\n')
		if line_to_check in unique_lines: # Skip if line is already in set
			duplicate_lines += 1
			continue
		else: # Write if new and append stripped line to list of seen lines
			target.write(line)
			unique_lines.add(line_to_check)

	# Be nice and close out the files
	target.close()
	data.close()

"""
Extract information from csv file saved by Biolapse interface.
"""
def get_labels(path_to_crop):
    # Load labels
    all_img = glob.glob(os.path.join(path_to_crop, 'zoom','*', '*.png'))
    all_img = sorted(all_img,key=stringSplitByNumbers)
    labels = []
    idx = []
    csv_file_ = os.path.join(path_to_crop, 'cells.csv')
    remove_duplicate_csv(csv_file_)
    csv_file_ = os.path.join(path_to_crop, 'cells_deduped.csv')
    with open(csv_file_) as csv_file:
        csv_reader = csv.reader(csv_file, delimiter=',')
        name_cell_ = ''
        first_run = True
        for row in csv_reader:
            idx_matching = [k for k in range(
                len(all_img)) if name_cell_ in all_img[k]]
            name_cell = row[0]
            if name_cell == name_cell_:
                # add previous cell labels
                if not(row[1] == ' '):
                    crt_time = int(row[1])
                    for k in range(prev_time, crt_time):
                        lab_cell.append(0)
                    prev_time = crt_time
                elif not(row[3] == ' '):
                    crt_time = int(row[3])
                    for k in range(prev_time, crt_time):
                        lab_cell.append(0)
                    prev_time = crt_time
                elif not(row[5] == ' '):
                    crt_time = int(row[5])
                    for k in range(prev_time, crt_time):
                        lab_cell.append(1)
                    prev_time = crt_time
                elif not(row[7] == ' '):
                    # import ipdb; ipdb.set_trace()
                    crt_time = int(row[7])
                    for k in range(prev_time, crt_time):
                        lab_cell.append(2)
                    prev_time = crt_time
                elif not(row[9] == ' '):
                    crt_time = int(row[9])
                    for k in range(prev_time, crt_time):
                        lab_cell.append(3)
                    prev_time = crt_time
                elif not(row[11] == ' '):
                    crt_time = int(row[11])
                    for k in range(prev_time, crt_time+1):
                        lab_cell.append(4)
                    prev_time = crt_time
                    # lab_cell.append(-1)
                    if not(len(lab_cell) == int(row[11])-start_time+1):
                        print('Warning: inconsistance in time acquisition.')
                    if len(idx_matching) != len(lab_cell):
                        # import ipdb; ipdb.set_trace()
                        print(
                            'Error: not same number of cell than number of labels found.')
                else:
                    print('Warning: no phase mentionned.')
                    print(row)
                    name_cell_ = ''
            else:
                # add corresponding index in image sequence
                if not(first_run):
                    if len(lab_cell) != 0:
                        labels.append(lab_cell)
                        idx.append(idx_matching)
                # create new label list
                lab_cell = []
                if not(row[1] == ' '):
                    start_time = int(row[1])
                    prev_time = start_time
                    first_run = False
                elif not(row[3] == ' '):
                    start_time = int(row[3])
                    prev_time = start_time
                    first_run = False
                elif not(row[5] == ' '):
                    start_time = int(row[5])
                    prev_time = start_time
                    first_run = False
                elif not(row[7] == ' '):
                    start_time = int(row[7])
                    prev_time = start_time
                    first_run = False
                elif not(row[9] == ' '):
                    start_time = int(row[9])
                    prev_time = start_time
                    first_run = False
                elif not(row[11] == ' '):
                    print('Warning: end before it starts.')
                    name_cell_ = ''
                    # start_time = int(row[1])
                    # lab_cell.append(-1)
                else:
                    print('Warning: no phase mentionned.')
                    print(row)
                    name_cell_ = ''
            name_cell_ = name_cell
        idx_matching = [k for k in range(
            len(all_img)) if name_cell_ in all_img[k]]
        if len(lab_cell) != 0:
            labels.append(lab_cell)
            idx.append(idx_matching)
    return labels, idx

## NN
import argparse
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.utils import data
from torchvision import datasets, transforms
import imageio

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.conv1 = nn.Conv2d(1, 10, kernel_size = 5, stride=1, padding=2, dilation=1)
        self.conv2 = nn.Conv2d(10, 20, kernel_size = 5, stride=1, padding=2, dilation=1)
        self.conv3 = nn.Conv2d(20, 30, kernel_size = 5, stride=1, padding=2, dilation=1)
        self.fc1 = nn.Linear(16*16*30, 500)
        self.fc2 = nn.Linear(500, 4)

    def forward(self, x):
        x = torch.squeeze(x,0)
        x = F.relu(self.conv1(x))
        x = F.max_pool2d(x, 2, 2)
        # x = nn.Dropout(p=0.5)(x)
        x = F.relu(self.conv2(x))
        x = F.max_pool2d(x, 2, 2)
        # x = nn.Dropout(p=0.5)(x)
        x = F.relu(self.conv3(x))
        x = F.max_pool2d(x, 2, 2)
        # x = nn.Dropout(p=0.5)(x)
        x = x.view(-1, self.num_flat_features(x))
        x = F.relu(self.fc1(x))
        # x = nn.Dropout(p=0.5)(x)
        x = x.view(-1, self.num_flat_features(x))
        x = self.fc2(x)
        return F.softmax(x, dim=1)

    def num_flat_features(self, x):
        size = x.size()[1:]  # all dimensions except the batch dimension
        num_features = 1
        for s in size:
            num_features *= s
        return num_features

"""
Apply blur operator onto points at position on the grid.
Operator is assumed to act separately on two parts of the image: im[:n/2,:] and im[n/2:,:].
Image makes operator always identifiable: minimum 2 sources id different part of the image.
"""
class CropGenertor(torch.utils.data.Dataset):
    def __init__(self, im, feature):
        self.L=feature.shape[0]
        self.feature=feature
        self.im=im

    def __len__(self):
        return self.L

    def __getitem__(self, i):
        return np.expand_dims(self.im[i],0), self.feature[i]

from torch.autograd import Variable
def train(model, device, train_loader, test_loader, optimizer, epoch, nbatch=20):
    model.train()
    use_cuda=torch.cuda.is_available()

    # data, target = train_loader
    # data_test, target_test = test_loader
    # nbImages = data.shape[0]

    cpt=0
    cum_loss_train=0
    for im, feat in train_loader:
        if use_cuda:
            im = Variable(im).type(torch.float32).cuda()
        else:
            im = Variable(im).type(torch.float32)

    # for batch_idx in range(nit):
        # randi = np.random.randint(0,nbImages,size=nbatch)
        # dataloc, targetloc = data[randi,:,:,:] ,target[randi]
        # dataloc.to(device)
        # targetloc.to(device)

        optimizer.zero_grad()
        # output = model(dataloc)
        output = model(im)
        if use_cuda:
            loss = nn.BCELoss()(output,Variable(feat.float()).cuda())
        else:
            loss = nn.BCELoss()(output,Variable(feat.float()))
        # loss = nn.BCELoss()(output,targetloc)
        # loss = nn.CrossEntropyLoss().cuda()(output,targetloc)
        loss.backward()
        optimizer.step()
        cum_loss_train+=loss.item()
        if cpt % 50 == 0:
            print('Train Epoch: {} [{}]\tLoss: {:.6f}'.format(
                epoch, cpt,loss.item()))
        cpt+=1

    loss_train = cum_loss_train/cpt
    # # Test
    # randi = np.random.randint(0,data_test.shape[0],size=nbatch)
    # dataloc_test, targetloc_test = data_test[randi,:], target_test[randi,:]
    # dataloc_test.to(device)
    # targetloc_test.to(device)
    # output = model(dataloc_test)
    # loss_test = nn.BCELoss()(output,targetloc_test)

    cpt=0
    cum_loss_test=0
    for image_test, label_test in test_loader:
        if use_cuda:
            image_test = Variable(image_test).type(torch.float32).cuda()
        else:
            image_test = Variable(image_test).type(torch.float32)
        output = model(image_test)
        if use_cuda:
            loss = nn.BCELoss()(output, Variable(label_test.float()).cuda())
        else:
            loss = nn.BCELoss()(output, Variable(label_test.float()))
        cum_loss_test += loss.item()
        cpt+=1
    loss_test=cum_loss_test/cpt

    # loss_test = nn.CrossEntropyLoss().cuda()(output,targetloc_test)
    print('Test loss: {:.6f}'.format(loss_test))
    return loss_train, loss_test

class Dataset():

    """Dataset class
        Attributes:
        dataPath (TYPE): Position of the folder containing all the images
        lenDataset (TYPE): Number of images present in this folder
    """

    def __init__(self, dataPath,lenDataset):
        'Initialization'
        self.dataPath = dataPath
        self.lenDataset = lenDataset

    def __len__(self):
        'Denotes the total number of samples'
        return self.lenDataset

    def __getitem__(self, index):
        'Generates one sample of data'
        images = np.zeros((self.lenDataset,1,128,128))
        labs = np.load(self.dataPath+os.sep+'feature.npy')
        labs = labs[:self.lenDataset]
        for i in range(self.lenDataset):
            images[i,0,:,:] = np.array(imageio.imread(self.dataPath+os.sep+str(i).zfill(5)+'.png'),dtype=float)/255.

        return torch.Tensor(images), torch.Tensor(labs)
    def getALLDATA(self):
        images = np.zeros((self.lenDataset,1,128,128))
        labs = np.load(self.dataPath+os.sep+'feature.npy')
        labs = labs[:self.lenDataset]
        for i in range(self.lenDataset):
            images[i,0,:,:] = np.array(imageio.imread(self.dataPath+os.sep+str(i).zfill(5)+'.png'),dtype=float)/255.

        return (torch.Tensor(images).cuda(), torch.Tensor(labs).cuda())

"""
Viterbi algorithm to predict the most likely Markov chain based on the probabilities of each time.

INPUT:
    - proba_pred: array of size (nSamples,nStates). nSamples is the number of sample and nStates is the number of possible states. Contains
    probability of each sample to be in each states.
    - nclasses: number of phase to predict.
    - transition: average probability of transition between classes. This possibly encode the restriction on the state path.
"""
def viterbi(proba_pred, nclasses, transition):
    # Check shape
    nSamples = proba_pred.shape[0]
    nStates = proba_pred.shape[1]
    bestPath = np.zeros(nSamples)
    psi = np.zeros((nStates,nSamples))
    c = np.zeros(nSamples)
    viterbi = np.zeros((nStates, nSamples))

    viterbi[:, 0] = proba_pred[0].T
    c[0] = 1.0/np.sum(viterbi[:, 0])
    viterbi[:, 0] = c[0] * viterbi[:, 0]  # apply the scaling factor
    psi[0] = 0

    for t in range(1, nSamples):
        for s in range(0, nStates):
            trans_p = viterbi[:, t-1] * transition[:, s]
            psi[s, t], viterbi[s, t] = max(
                enumerate(trans_p), key=operator.itemgetter(1))
            # viterbi[s,t] = viterbi[s,t]*probaPredictor(model,featMark[t])[s]
            viterbi[s,t] = viterbi[s,t]*proba_pred[t,s]

        if np.sum(viterbi[:,t])==0:
            import ipdb; ipdb.set_trace()
        c[t] = 1.0/np.sum(viterbi[:,t])
        viterbi[:,t] = c[t] * viterbi[:,t]

    bestPath[nSamples-1] =  viterbi[:,nSamples-1].argmax()
    for t in range(nSamples-1,0,-1):
        bestPath[t-1] = psi[int(bestPath[t]),t]

    return bestPath

def buildTransitionMatrix(vectorEpsilon):
    n = vectorEpsilon.shape[0]
    transitionMat = np.zeros([n,n])
    for i in range(n-1):
        transitionMat[i,i] = 1 - vectorEpsilon[i]
        transitionMat[i,i+1] = vectorEpsilon[i]

    # transitionMat[n-1,n-1] = 1
    transitionMat[n-1,n-1] = 1-vectorEpsilon[-1]
    transitionMat[n-1,0] = vectorEpsilon[-1]
    return transitionMat

"""
Detect connected components with minimal area defined by 'minimalSize'.
Return list containing the number of the component, its barycenter and its area. It also returns
'compo', an image where each connected componenent is defined by the number given in localCells.
"""
def find_connected_components(masks,minimalSize):
	binary = cv2.threshold(masks, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]
	num,compo=cv2.connectedComponents(binary)
	localCells = []
	for chunckNo in range(1,num):
		if len(np.where(compo==chunckNo)[0]) < minimalSize:
			compo[np.where(compo==chunckNo)] = 0
		else:
			localCells += [[chunckNo,(np.where(compo==chunckNo)[0].mean(),np.where(compo==chunckNo)[1].mean()),len(np.where(compo==chunckNo)[0])]]
	return localCells, compo

"""
Display train and testing dataset for classification
"""
def display_classification(im,feature_pred=[],feature_true=[],leg=[],save=True,name='visu',L=200):
    if len(feature_true)!=0:
        l=feature_true.shape[1]
    elif len(feature_pred)!=0:
        l=feature_pred.shape[1]
    else:
        l=-1
    if len(leg)!=l:
        leg=[]
    K=np.minimum(im.shape[0],9)
    plt.figure(1)
    for k in range(K):
        plt.subplot(331+k)
        plt.imshow(im[k])
        if len(feature_true)==0 and len(feature_pred)!=0:
            if len(leg)!=0:
                plt.title('Image {0} -- Prediction: {1}'.format(k,(leg[np.where(feature_pred[k])[0][0]])))
            else:
                plt.title('Image {0} -- Prediction: {1}'.format(k,(np.where(feature_pred[k])[0][0])))
        if len(feature_true)!=0 and len(feature_pred)==0:
            if len(leg)!=0:
                plt.title('Image {0} -- True: {1}'.format(k,leg[np.where(feature_true[k])[0][0]]))
            else:
                plt.title('Image {0} -- True: {1}'.format(k,(np.where(feature_true[k])[0][0])))
        if len(feature_true)!=0 and len(feature_pred)!=0:
            if len(leg)!=0:
                plt.title('Image {0} -- Prediction: {1} -- True: {2}'.format(k,(leg[np.where(feature_pred[k])[0][0]]),(leg[np.where(feature_true[k])[0][0]])))
            else:
                plt.title('Image {0} -- Prediction: {1} -- True: {2}'.format(k,(np.where(feature_pred[k])[0][0]),(np.where(feature_true[k])[0][0])))
        if len(feature_true)==0 and len(feature_pred)==0:
            raise ValueError('No feature found.')

    if save:
        if not os.path.exists(os.path.join('..','Data','Classification',name)):
            os.makedirs(os.path.join('..','Data','Classification',name))

        for k in range(np.minimum(im.shape[0],L)):
            plt.figure(2)
            plt.imshow(im[k])
            if len(feature_true)==0 and len(feature_pred)!=0:
                if len(leg)!=0:
                    plt.title('Image {0} -- Prediction: {1}'.format(k,(leg[np.where(feature_pred[k])[0][0]])))
                else:
                    plt.title('Image {0} -- Prediction: {1}'.format(k,(np.where(feature_pred[k])[0][0])))
            if len(feature_true)!=0 and len(feature_pred)==0:
                if len(leg)!=0:
                    plt.title('Image {0} -- True: {1}'.format(k,(leg[np.where(feature_true[k])[0][0]])))
                else:
                    plt.title('Image {0} -- True: {1}'.format(k,(np.where(feature_true[k])[0][0])))
            if len(feature_true)!=0 and len(feature_pred)!=0:
                if len(leg)!=0:
                    plt.title('Image {0} -- Prediction: {1} -- True: {2}'.format(k,(leg[np.where(feature_pred[k])[0][0]]),(leg[np.where(feature_true[k])[0][0]])))
                else:
                    plt.title('Image {0} -- Prediction: {1} -- True: {2}'.format(k,(np.where(feature_pred[k])[0][0]),(np.where(feature_true[k])[0][0])))
            if len(feature_true)==0 and len(feature_pred)==0:
                raise ValueError('No feature found.')
            plt.savefig(os.path.join('..','Data','Classification',name,'im_'+str(k).zfill(5)+'.png'))

max_stars_repo_name: lebrat/Biolapse

max_issues_repo_name: lebrat/Biolapse

max_forks_repo_name: lebrat/Biolapse


text:
#Programmer: Chris Tralie
#Purpose: To extract cover song alignments for use in the GUI
import numpy as np
import sys
sys.path.append("../")
sys.path.append("../SequenceAlignment")
import os
import glob
import scipy.io as sio
import skimage
import skimage.io
import time
import matplotlib.pyplot as plt
from CSMSSMTools import *
from BlockWindowFeatures import *
from Onsets import *
from AudioIO import *
import json
import argparse
import base64

def getBase64File(filename):
    fin = open(filename, "rb")
    b = fin.read()
    b = base64.b64encode(b)
    fin.close()
    return b

def getBase64PNGImage(D, cmapstr):
    c = plt.get_cmap(cmapstr)
    D = np.round(255.0*D/np.max(D))
    C = c(np.array(D, dtype=np.int32))
    C = np.array(np.round(C*255), dtype=np.uint8)
    skimage.io.imsave("temp.png", C)
    b = getBase64File("temp.png")
    os.remove("temp.png")
    return b

#http://stackoverflow.com/questions/1447287/format-floats-with-standard-json-module
class PrettyFloat(float):
    def __repr__(self):
        return '%.4g' % self
def pretty_floats(obj):
    if isinstance(obj, float):
        return PrettyFloat(obj)
    elif isinstance(obj, dict):
        return dict((k, pretty_floats(v)) for k, v in obj.items())
    elif isinstance(obj, (list, tuple)):
        return map(pretty_floats, obj)
    return obj

def compareTwoSongsJSON(filename1, TempoBias1, filename2, TempoBias2, hopSize, FeatureParams, CSMTypes, Kappa, outfilename, song1name = 'Song 1', song2name = 'Song 2'):
    print("Getting features for %s..."%filename1)
    (XAudio, Fs) = getAudioLibrosa(filename1)
    (tempo, beats1) = getBeats(XAudio, Fs, TempoBias1, hopSize, filename1)
    (Features1, O1) = getBlockWindowFeatures((XAudio, Fs, tempo, beats1, hopSize, FeatureParams))

    print("Getting features for %s..."%filename2)
    (XAudio, Fs) = getAudioLibrosa(filename2)
    (tempo, beats2) = getBeats(XAudio, Fs, TempoBias2, hopSize, filename2)
    (Features2, O2) = getBlockWindowFeatures((XAudio, Fs, tempo, beats2, hopSize, FeatureParams))

    print("Feature Types: ", Features1.keys())

    beats1 = beats1*hopSize/float(Fs)
    beats2 = beats2*hopSize/float(Fs)
    Results = {'song1name':song1name, 'song2name':song2name, 'hopSize':hopSize, 'FeatureParams':FeatureParams, 'Kappa':Kappa, 'CSMTypes':CSMTypes, 'beats1':pretty_floats(beats1.tolist()), 'beats2':pretty_floats(beats2.tolist())}
    #Do each feature individually
    FeatureCSMs = {}
    for FeatureName in Features1:
        print("Doing %s..."%FeatureName)
        res =  getCSMSmithWatermanScores(Features1[FeatureName], O1, Features2[FeatureName], O2, Kappa, CSMTypes[FeatureName], True)
        CSMs = {}
        CSMs['D'] = getBase64PNGImage(res['D'], 'afmhot')
        CSMs['CSM'] = getBase64PNGImage(res['CSM'], 'afmhot')
        CSMs['DBinary'] = getBase64PNGImage(1-res['DBinary'], 'gray')
        CSMs['score'] = res['score']
        FeatureCSMs[FeatureName] = CSMs

    #Do OR Merging
    print("Doing OR Merging...")
    res = getCSMSmithWatermanScoresORMerge(Features1, O1, Features2, O2, Kappa, CSMTypes, True)
    CSMs = {}
    CSMs['D'] = getBase64PNGImage(res['D'], 'afmhot')
    CSMs['CSM'] = getBase64PNGImage(1-res['DBinary'], 'gray')
    CSMs['DBinary'] = CSMs['CSM']
    CSMs['score'] = res['score']
    CSMs['FeatureName'] = 'ORFusion'
    FeatureCSMs['ORFusion'] = CSMs

    #Do cross-similarity fusion
    print("Doing similarity network fusion...")
    K = 20
    NIters = 3
    res = getCSMSmithWatermanScoresEarlyFusionFull(Features1, O1, Features2, O2, Kappa, K, NIters, CSMTypes, True)
    CSMs = {}
    CSMs['D'] = getBase64PNGImage(res['D'], 'afmhot')
    CSMs['CSM'] = getBase64PNGImage(res['CSM'], 'afmhot')
    CSMs['DBinary'] = getBase64PNGImage(1-res['DBinary'], 'gray')
    CSMs['score'] = res['score']
    FeatureCSMs['SNF'] = CSMs

    Results['FeatureCSMs'] = FeatureCSMs
    print("Saving results...")
    Results['beats1'] = list(Results['beats1'])
    Results['beats2'] = list(Results['beats2'])
    for feature in ['MFCCs', 'SSMs', 'Chromas', 'ORFusion', 'SNF']:
        for imgtype in ['D', 'CSM', 'DBinary']:
            Results['FeatureCSMs'][feature][imgtype] = "data:image/png;base64, " + Results['FeatureCSMs'][feature][imgtype].decode('ascii')

    #Add music as base64 files
    path, ext = os.path.splitext(filename1)
    Results['file1'] = "data:audio/%s;base64, "%ext[1::] + getBase64File(filename1).decode('ascii')
    path, ext = os.path.splitext(filename2)
    Results['file2'] = "data:audio/%s;base64, "%ext[1::] + getBase64File(filename2).decode('ascii')
    fout = open(outfilename, "w")
    fout.write(json.dumps(Results))
    fout.close()

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--filename1', type=str, required=True, help="Path to audio file for first song")
    parser.add_argument('--filename2', type=str, required=True, help="Path to audio file for first song")
    parser.add_argument('--jsonfilename', type=str, required=True, help="Output name of json file")
    parser.add_argument('--artist1', type=str, default="Artist 1", help="Name of the artist of the first song")
    parser.add_argument('--artist2', type=str, default="Artist 2", help="Name of the artist of the second song")
    parser.add_argument('--tempobias1', type=int, default=120, help="Tempo bias of the beat tracker for the first song")
    parser.add_argument('--tempobias2', type=int, default=120, help="Tempo bias of the beat tracker for the second song")
    parser.add_argument('--hopsize', type=int, default=512, help="Hop size to use for the features")
    parser.add_argument('--kappa', type=float, default=0.1, help="Nearest neighbor threshold")
    opt = parser.parse_args()

    FeatureParams = {'MFCCBeatsPerBlock':20, 'MFCCSamplesPerBlock':200, 'DPixels':50, 'ChromaBeatsPerBlock':20, 'ChromasPerBlock':40}
    CSMTypes = {'MFCCs':'Euclidean', 'SSMs':'Euclidean', 'SSMsDiffusion':'Euclidean', 'Geodesics':'Euclidean', 'Jumps':'Euclidean', 'Curvs':'Euclidean', 'Tors':'Euclidean', 'CurvsSS':'Euclidean', 'TorsSS':'Euclidean', 'D2s':'EMD1D', 'Chromas':'CosineOTI'}

    compareTwoSongsJSON(opt.filename1, opt.tempobias1, opt.filename2, opt.tempobias2, opt.hopsize, FeatureParams, CSMTypes, opt.kappa, opt.jsonfilename, opt.artist1, opt.artist2)

ext: py

reason: import numpy,import scipy


text:
import math
from abc import ABCMeta, abstractmethod
from enum import Enum
from pathlib import Path

import audioread
import librosa
import numpy
import os
from lazy import lazy
from numpy import ndarray, mean, std, vectorize, dot
from typing import List, Optional, Tuple, Callable

from speechless.tools import name_without_extension, mkdir, write_text, log

# This is the place where pyaudio is heavily used

# how we devide the frequency bank, linear or mf
class SpectrogramFrequencyScale(Enum):
    linear = "linear"
    mel = "mel"

# how we calculate the power of a frequency band
class SpectrogramType(Enum):
    power = "power"
    amplitude = "amplitude"
    power_level = "power level"

def z_normalize(array: ndarray) -> ndarray:
    return (array - mean(array)) / std(array)

class PositionalLabel:
    def __init__(self, labeled_sections: List[Tuple[str, Tuple[float, float]]]):
        if not labeled_sections:
            raise ValueError("Sections must be specified.")

        if any(range is None for label, range in labeled_sections):
            raise ValueError("Range must be specified.")

        self.labeled_sections = labeled_sections
        self.labels = [word for word, range in labeled_sections]
        self.label = " ".join(word for word in self.labels)

    def convert_range_to_seconds(self, original_sample_rate: int) -> 'PositionalLabel':
        # start: the start frame (in time domain)
        return PositionalLabel(list(
            (word, (start_end[0] / original_sample_rate, start_end[1] / original_sample_rate))
            for word, start_end in self.labeled_sections))

    def with_corrected_labels(self, correction: Callable[[str], str]) -> 'PositionalLabel':
        # section is the label of that section
        return PositionalLabel([(correction(section), range) for section, range in self.labeled_sections])

    def serialize(self) -> str:
        return "\n".join("{}|{}|{}".format(label, start, end) for label, (start, end) in
                         self.labeled_sections)

    @staticmethod
    def deserialize(serialized: str) -> 'PositionalLabel':
        return PositionalLabel(list((label, (float(start), float(end)))
                                    for label, start, end in
                                    map(lambda item: item.split("|"), serialized.splitlines())))

# a transformed example in frequency domain
class LabeledSpectrogram:
    __metaclass__ = ABCMeta

    def __init__(self, id: str, label: str):
        self.label = label
        self.id = id

    @abstractmethod
    def z_normalized_transposed_spectrogram(self) -> ndarray: raise NotImplementedError

class LabeledExample(LabeledSpectrogram):
    def __init__(self,
                 get_raw_audio: Callable[[], ndarray],
                 sample_rate: int = 16000,
                 id: Optional[str] = None,
                 label: Optional[str] = "nolabel",
                 # 32 frames
                 fourier_window_length: int = 512,
                 # or 128/16 = 8 ms
                 hop_length: int = 128,
                 mel_frequency_count: int = 128,
                 label_with_tags: str = None,
                 positional_label: Optional[PositionalLabel] = None):
        super().__init__(id=id, label=label)

        # The default values for hop_length and fourier_window_length are powers of 2 near the values specified in the wave2letter paper.
        self.get_raw_audio = get_raw_audio
        self.sample_rate = sample_rate
        self.fourier_window_length = fourier_window_length
        self.hop_length = hop_length
        self.mel_frequency_count = mel_frequency_count
        self.label_with_tags = label_with_tags
        self.positional_label = positional_label

    def tag_count(self, tag: str) -> int:
        return self.label_with_tags.count(tag)

    def _power_spectrogram(self) -> ndarray:
        return self._amplitude_spectrogram() ** 2

    def _amplitude_spectrogram(self) -> ndarray:
        return abs(self._complex_spectrogram())

    def _complex_spectrogram(self) -> ndarray:
        return librosa.stft(y=self.get_raw_audio(), n_fft=self.fourier_window_length, hop_length=self.hop_length)

    def mel_frequencies(self) -> List[float]:
        # according to librosa.filters.mel code
        return librosa.mel_frequencies(self.mel_frequency_count + 2, fmax=self.sample_rate / 2)

    def _convert_spectrogram_to_mel_scale(self, linear_frequency_spectrogram: ndarray) -> ndarray:
        return dot(
            librosa.filters.mel(sr=self.sample_rate, n_fft=self.fourier_window_length, n_mels=self.mel_frequency_count),
            linear_frequency_spectrogram)

    def highest_detectable_frequency(self) -> float:
        return self.sample_rate / 2

    def spectrogram(self, type: SpectrogramType = SpectrogramType.power_level,
                    frequency_scale: SpectrogramFrequencyScale = SpectrogramFrequencyScale.linear) -> ndarray:
        def spectrogram_by_type():
            if type == SpectrogramType.power:
                return self._power_spectrogram()
            if type == SpectrogramType.amplitude:
                return self._amplitude_spectrogram()
            if type == SpectrogramType.power_level:
                return self._power_level_from_power_spectrogram(self._power_spectrogram())

            raise ValueError(type)

        s = spectrogram_by_type()
        # TODO: simplify it:
        # time_series, sample_rate = librosa.core.load("7a.wav",sr=20000)
        # librosa.feature.melspectrogram(time_series, sr=20000, n_mels=23, n_fft=500, hop_length=200)
        return self._convert_spectrogram_to_mel_scale(s) if frequency_scale == SpectrogramFrequencyScale.mel else s

    # each sample is normalized independently, across frames
    def z_normalized_transposed_spectrogram(self):
        """
        :return: Array with shape (time, frequencies)
        """
        # transposed: time steps * frame (frequency bins), normalization inside a frame
        return z_normalize(self.spectrogram(frequency_scale=SpectrogramFrequencyScale.mel).T)
    # shape of the spectrogram : number of bank * time steps
    def frequency_count_from_spectrogram(self, spectrogram: ndarray) -> int:
        return spectrogram.shape[0]

    def time_step_count(self) -> int:
        return self.spectrogram().shape[1]

    def time_step_rate(self) -> float:
        return self.time_step_count() / self.duration_in_s

    @staticmethod
    def _power_level_from_power_spectrogram(spectrogram: ndarray) -> ndarray:
        # default value for min_decibel found by experiment (all values except for 0s were above this bound)
        def power_to_decibel(x, min_decibel: float = -150) -> float:
            if x == 0:
                return min_decibel
            l = 10 * math.log10(x)
            return min_decibel if l < min_decibel else l

        return vectorize(power_to_decibel)(spectrogram)

    def reconstructed_audio_from_spectrogram(self) -> ndarray:
        return librosa.istft(self._complex_spectrogram(), win_length=self.fourier_window_length,
                             hop_length=self.hop_length)

    @lazy
    def duration_in_s(self) -> float:
        return len(self.get_raw_audio()) / self.sample_rate

    def __str__(self) -> str:
        return self.id + (": {}".format(self.label) if self.label else "")

class LabeledExampleFromFile(LabeledExample):
    def __init__(self,
                 audio_file: Path,
                 id: Optional[str] = None,
                 sample_rate_to_convert_to: int = 16000,
                 label: Optional[str] = "nolabel",
                 fourier_window_length: int = 512,
                 hop_length: int = 128,
                 mel_frequency_count: int = 128,
                 label_with_tags: str = None,
                 positional_label: Optional[PositionalLabel] = None):
        # The default values for hop_length and fourier_window_length are powers of 2 near the values specified in the wave2letter paper.

        if id is None:
            id = name_without_extension(audio_file)

        self.audio_file = audio_file

        super().__init__(
            id=id, get_raw_audio=lambda: librosa.load(str(self.audio_file), sr=self.sample_rate)[0],
            label=label, sample_rate=sample_rate_to_convert_to,
            fourier_window_length=fourier_window_length, hop_length=hop_length, mel_frequency_count=mel_frequency_count,
            label_with_tags=label_with_tags, positional_label=positional_label)

    @property
    def audio_directory(self):
        return Path(self.audio_file.parent)

    @lazy
    def original_sample_rate(self) -> int:
        return LabeledExampleFromFile.file_sample_rate(self.audio_file)

    @staticmethod
    def file_sample_rate(audio_file: Path) -> int:
        with audioread.audio_open(os.path.realpath(str(audio_file))) as input_file:
            return input_file.samplerate

    @lazy
    def duration_in_s(self) -> float:
        try:
            return librosa.get_duration(filename=str(self.audio_file))
        except Exception as e:
            log("Failed to get duration of {}: {}".format(self.audio_file, e))
            return 0

    def sections(self) -> Optional[List[LabeledExample]]:
        if self.positional_label is None:
            return None

        audio = self.get_raw_audio()

        def section(label, start, end):
            return LabeledExample(
                get_raw_audio=lambda: audio[int(start * self.sample_rate):int(end * self.sample_rate)],
                label=label,
                sample_rate=self.sample_rate,
                fourier_window_length=self.fourier_window_length, hop_length=self.hop_length,
                mel_frequency_count=self.mel_frequency_count)
        # fine grained samples, each section corresponds to one sample
        return [section(label, start, end) for label, (start, end) in self.positional_label.labeled_sections]

# use pre-compuated spectrogram to speed up the training
class CachedLabeledSpectrogram(LabeledSpectrogram):
    def __init__(self, original: LabeledSpectrogram, spectrogram_cache_directory: Path):
        super().__init__(id=original.id, label=original.label)
        self.original = original
        self.spectrogram_cache_file = spectrogram_cache_directory / "{}.npy".format(original.id)

    def z_normalized_transposed_spectrogram(self) -> ndarray:
        if not self.is_cached():
            return self._calculate_and_save_spectrogram()

        return self._load_from_cache()

    def _load_from_cache(self):
        try:
            return numpy.load(str(self.spectrogram_cache_file))
        except ValueError:
            log("Recalculating cached file {} because loading failed.".format(self.spectrogram_cache_file))
            return self._calculate_and_save_spectrogram()

    def _calculate_and_save_spectrogram(self):
        spectrogram = self.original.z_normalized_transposed_spectrogram()
        self._save_to_cache(spectrogram)
        return spectrogram

    def _save_to_cache(self, spectrogram: ndarray):
        numpy.save(str(self.spectrogram_cache_file), spectrogram)

    def is_cached(self):
        return self.spectrogram_cache_file.exists()

    def repair_cached_file_if_incorrect(self):
        if not self.is_cached():
            self._calculate_and_save_spectrogram()
            return

        from_cache = self._load_from_cache()
        calculated = self.original.z_normalized_transposed_spectrogram()
        try:
            numpy.testing.assert_almost_equal(calculated, from_cache, decimal=1)
        except AssertionError as e:
            self.move_incorrect_cached_file_to_backup_location_and_save_error(str(e))
            self._save_to_cache(calculated)

    def move_incorrect_cached_file_to_backup_location_and_save_error(self, error_text: str):
        parent_directory = Path(self.spectrogram_cache_file.parent)
        incorrect_cached_backup_directory = Path(parent_directory.parent / (parent_directory.name + "-incorrect"))
        mkdir(incorrect_cached_backup_directory)
        incorrect_backup_file = incorrect_cached_backup_directory / self.spectrogram_cache_file.name
        incorrect_backup_message_file = incorrect_cached_backup_directory / (
            name_without_extension(self.spectrogram_cache_file) + "-error.txt")
        write_text(incorrect_backup_message_file, error_text)
        self.spectrogram_cache_file.rename(incorrect_backup_file)

reason: import numpy,from numpy


text:
#from sipyco.pc_rpc import Client

from artiq.experiment import *
from artiq.master.databases import DeviceDB
from artiq.master.worker_db import DeviceManager, DatasetManager
#from artiq.master.worker_impl import CCB, Scheduler

import numpy as np

class ARTIQ_api(object):
    """
    An API for the ARTIQ box.
    Directly accesses the hardware on the box without having to use artiq_master.
    """

    def __init__(self):
        self.devices = DeviceDB('C:\\Users\\EGGS1\\Documents\\ARTIQ\\artiq-master\\device_db.py')
        self.device_manager = DeviceManager(self.devices)
        self._getDevices()
        self._initializeDevices()

    # SETUP
    def _getDevices(self):
        # get core
        self.core = self.device_manager.get("core")
        self.core_dma = self.device_manager.get("core_dma")

        # store devices in dictionary where device name is key
        # and device itself is value
        self.device_db = self.devices.get_device_db()
        self.ttlout_list = {}
        self.ttlin_list = {}
        self.dds_list = {}
        self.urukul_list = {}
        self.dac = None
        self.sampler = None

        # assign names and devices
        for name, params in self.device_db.items():
            # only get devices with named class
            if 'class' not in params:
                continue
            # set device as attribute
            devicetype = params['class']
            device = self.device_manager.get(name)
            if devicetype == 'TTLInOut':
                self.ttlin_list[name] = device
            elif devicetype == 'TTLOut':
                self.ttlout_list[name] = device
            elif devicetype == 'AD9910':
                self.dds_list[name] = device
            elif devicetype == 'CPLD':
                self.urukul_list[name] = device
            elif devicetype == 'Zotino':
                self.dac = device
            elif devicetype == 'Sampler':
                self.sampler = device

    def _initializeDevices(self):
        """
        Initialize devices that need to be initialized.
        """
        # set ttlinout devices to be input
        # self.core.reset()
        # for device in self.ttlin_list.values():
        #     try:
        #         device.input()
        #     except RTIOUnderflow:
        #         self.core.break_realtime()
        #         device.input()
        # initialize DDSs
        # self.initializeDDSAll()
        # one-off device init
        self.initializeDAC()

    # PULSE SEQUENCING
    @kernel
    def record(self, sequencename):
        self.core.reset()
        with self.core_dma.record(sequencename):
            for i in range(50):
                with parallel:
                    self.ttlout_list[0].pulse(1*ms)
                    self.ttlout_list[1].pulse(1*ms)
                delay(1.0*ms)

    def record2(self, ttl_seq, dds_seq, sequencename):
        """
        Processes the TTL and DDS sequence into a format
        more easily readable and processable by ARTIQ.
        """
        # TTL
        ttl_times = list(ttl_seq.keys())
        ttl_commands = list(ttl_seq.values())
        # DDS
        dds_times = list(dds_seq.keys())
        dds_commands = list(dds_seq.values())
        # send to kernel
        self._record(ttl_times, ttl_commands, dds_times, dds_commands, sequencename)

    @kernel
    def eraseSequence(self, sequencename):
        """
        Erase the given pulse sequence from DMA.
        """
        self.core.reset()
        self.core_dma.erase(sequencename)

    # TTL
    def setTTL(self, ttlname, state):
        """
        Manually set the state of a TTL.
        """
        try:
            dev = self.ttlout_list[ttlname]
        except KeyError:
            raise Exception('Invalid device name.')
        self._setTTL(dev, state)

    @kernel
    def _setTTL(self, dev, state):
        self.core.reset()
        if state:
            dev.on()
        else:
            dev.off()

    # DDS
    @kernel
    def initializeDDSAll(self):
        '''
        Initialize all DDS channels and their CPLDs.
        '''
        #todo: fix, don't use dicts
        self.core.reset()
        # initialize urukul cpld
        for device in list(self.urukul_list.values()):
            try:
                device.init()
            except RTIOUnderflow:
                self.core.break_realtime()
                device.init()
        # initialize each DDS channel
        self.core.reset()
        for device in self.dds_list.values():
            try:
                device.init()
            except RTIOUnderflow:
                self.core.break_realtime()
                device.init()

    def initializeDDS(self, dds_name):
        dev = self.dds_list[dds_name]
        self._initializeDDS(dev)

    @kernel
    def _initializeDDS(self, dev):
        self.core.reset()
        dev.init()

    def toggleDDS(self, dds_name, state):
        """
        Toggle a DDS using the RF switch.
        """
        dev = self.dds_list[dds_name]
        self._toggleDDS(dev, state)

    @kernel
    def _toggleDDS(self, dev, state):
        self.core.reset()
        dev.cfg_sw(state)

    def setDDS(self, dds_name, param, val):
        """
        Manually set the state of a DDS.
        """
        dev = self.dds_list[dds_name]
        if param == 0:
            self._setDDSFreq(dev, val)
        elif param == 1:
            self._setDDSAmpl(dev, val)
        elif param == 2:
            self._setDDSPhase(dev, val)

    @kernel
    def _setDDSFreq(self, dev, ftw):
        self.core.reset()
        dev.set_ftw(ftw)

    @kernel
    def _setDDSAmpl(self, dev, asf):
        self.core.reset()
        dev.set_asf(asf)

    @kernel
    def _setDDSPhase(self, dev, pow):
        self.core.reset()
        dev.set_pow(pow)

    def setDDSAtt(self, dds_name, att_mu):
        """
        Set the DDS attenuation.
        """
        dev = self.dds_list[dds_name]
        self._setDDSAtt(dev, att_mu)

    @kernel
    def _setDDSAtt(self, dev, att_mu):
        self.core.reset()
        dev.set_att_mu(att_mu)

    def readDDS(self, dds_name, reg, length):
        """
        Read the value of a DDS register.
        """
        dev = self.dds_list[dds_name]
        if length == 16:
            return self._readDDS16(dev, reg)
        elif length == 32:
            return self._readDDS32(dev, reg)
        elif length == 64:
            return self._readDDS64(dev, reg)

    @kernel
    def _readDDS16(self, dev, reg):
        self.core.reset()
        return dev.read16(reg)

    @kernel
    def _readDDS32(self, dev, reg):
        self.core.reset()
        return dev.read32(reg)

    @kernel
    def _readDDS64(self, dev, reg):
        self.core.reset()
        return dev.read64(reg)

    # DAC
    @kernel
    def initializeDAC(self):
        """
        Initialize the DAC.
        """
        self.core.reset()
        self.dac.init()

    @kernel
    def setDAC(self, channel_num, volt_mu):
        """
        Set the voltage of a DAC register.
        """
        self.core.reset()
        self.dac.write_dac_mu(channel_num, volt_mu)
        self.dac.load()

    @kernel
    def setDACGain(self, channel_num, gain_mu):
        """
        Set the gain of a DAC channel.
        """
        self.core.reset()
        self.dac.write_gain_mu(channel_num, gain_mu)
        self.dac.load()

    @kernel
    def setDACOffset(self, channel_num, volt_mu):
        """
        Set the voltage of a DAC offset register.
        """
        self.core.reset()
        self.dac.write_offset_mu(channel_num, volt_mu)
        self.dac.load()

    @kernel
    def setDACGlobal(self, word):
        """
        Set the OFSx registers on the AD5372.
        """
        self.core.reset()
        self.dac.write_offset_dacs_mu(word)

    @kernel
    def readDAC(self, channel_num, address):
        """
        Read the value of one of the DAC registers.
        :param channel_num: Channel to read from
        :param address: Register to read from
        :return: the value of the register
        """
        self.core.reset()
        reg_val = self.dac.read_reg(channel_num, address)
        return reg_val

    @kernel
    def initializeSampler(self):
        """
        Initialize the Sampler.
        """
        self.core.reset()
        self.sampler.init()

    @kernel
    def setSamplerGain(self, channel_num, gain_mu):
        """
        Set the gain for a sampler channel.
        :param channel_num: Channel to set
        :param gain_mu: Register to read from
        :return: the value of the register
        """
        self.core.reset()
        self.sampler.set_gain_mu(channel_num, gain_mu)

    @kernel
    def getSamplerGains(self):
        """
        Get the gain of all sampler channels.
        :return: the sample channel gains.
        """
        self.core.reset()
        return self.sampler.get_gains_mu()

    @kernel
    def readSampler(self, sampleArr):
        """
        Set the gain of
        :param channel_num: Channel to set
        :param gain_mu: Register to read from
        :return: the value of the register
        """
        self.core.reset()
        return self.sampler.sample_mu(sampleArr)

reason: import numpy


text:
#######################################################################
#
#  import_bobbi.py
#
#  Code: Robin van Emden - robin at pavlov.tech - 2017
#
#  Totem Bobbi: http://www.totemopenhealth.com/
#
#######################################################################

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import os
from pathlib import Path

###################### Set csv filename ###############################

csv_file = "data/8-sdataHR.csv"

#######################################################################

# Retrieve some file info
f = Path(csv_file)                                      # path object
f_no_ext = os.path.join(str(f.parent), f.stem)          # remove extension

# Save Pickled DF - faster to save and load later
is_cached = Path(f_no_ext +".df")                       # check for cache file
if is_cached.is_file():                                 # cache exists?
    dataset = pd.read_pickle(f_no_ext +".df")           # yes - load cache

else:
    dataset = pd.read_csv( f_no_ext +".csv",            # no  - load csv
                index_col=False,
                sep=';',
                dtype={'ms': np.int32,
                       'heartrate': np.int16,
                       'BPM': np.int16 } )

    dataset['hart'] = dataset['heartrate']              # set hart
    # dataset['hart'] = -dataset['hart']+1000           # if needed, flip
    dataset.to_pickle(f_no_ext +".df")                  # save to cache

# Save small subset of data as a sample
dataset = dataset[(dataset['ms'] >= 104000)
                    & (dataset['ms'] <= 115000)]

dataset = dataset.reset_index(drop=True)                # Start index at 0
dataset.to_csv(f_no_ext + "_sample.csv"
               ,sep=';',index=False)                    # Save CSV sample
dataset.to_pickle(f_no_ext +"_sample.df")               # Save DF sample

# Lets plot our sample, see if its ok
plt.title("Heart Rate Signal")                          #The title of our plot
plt.plot(dataset.ms,dataset.hart)                       #Draw the plot object
plt.show()                                              #Display the plot

max_stars_repo_name: robinvanemden/TotemBobbiScripts

max_issues_repo_name: robinvanemden/TotemBobbiScripts

max_forks_repo_name: robinvanemden/TotemBobbiScripts


text:
"""
LANGUAGE: PYTHON
AUTHOR: cdefga
GITHUB: https://github.com/cdefga
"""

import numpy as np

def n_dimensional_distance(a: np.ndarray, b: np.ndarray) -> float:
    assert a.shape == b.shape, 'dimension of input must be equal'
    dist = np.sqrt(np.sum((a - b)**2))
    return dist

if __name__ == '__main__':
    a = np.array([1,2,3,4])
    b = np.array([1,2,3,5])
    dist = n_dimensional_distance(a, b)
    print(dist)


text:
# --------------------------------------------------------
# Faster R-CNN
# Licensed under The MIT License [see LICENSE for details]
# Written by Xinlei Chen
# --------------------------------------------------------
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
from model.config import cfg
from model.bbox_transform import bbox_transform_inv, clip_boxes
import numpy.random as npr

import torch
from torch.autograd import Variable

def proposal_top_layer(rpn_cls_prob, rpn_bbox_pred, im_info, _feat_stride, anchors, num_anchors):
  """A layer that just selects the top region proposals
     without using non-maximal suppression,
     For details please see the technical report
  """
  rpn_top_n = cfg.TEST.RPN_TOP_N

  scores = rpn_cls_prob[:, :, :, num_anchors:]

  rpn_bbox_pred = rpn_bbox_pred.view(-1, 4)
  scores = scores.contiguous().view(-1, 1)

  length = scores.size(0)
  if length < rpn_top_n:
    # Random selection, maybe unnecessary and loses good proposals
    # But such case rarely happens
    top_inds = torch.from_numpy(npr.choice(length, size=rpn_top_n, replace=True)).long().cuda()
  else:
    top_inds = scores.sort(0, descending=True)[1]
    top_inds = top_inds[:rpn_top_n]
    top_inds = top_inds.view(rpn_top_n)

  # Do the selection here
  anchors = anchors[top_inds.data, :].contiguous()
  rpn_bbox_pred = rpn_bbox_pred[top_inds.data, :].contiguous()
  scores = scores[top_inds.data].contiguous()

  # Convert anchors into proposals via bbox transformations
  proposals = bbox_transform_inv(anchors, rpn_bbox_pred)

  # Clip predicted boxes to image
  proposals = clip_boxes(proposals, im_info[:2])

  # Output rois blob
  # Our RPN implementation only supports a single input image, so all
  # batch inds are 0
  batch_inds = Variable(proposals.data.new(proposals.size(0), 1).zero_())
  blob = torch.cat((batch_inds, proposals), 1)
  return blob, scores

max_stars_repo_name: zhbli/pytorch-faster-rcnn

max_issues_repo_name: zhbli/pytorch-faster-rcnn

max_forks_repo_name: zhbli/pytorch-faster-rcnn


text:
import os,sys,inspect
current_dir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
parent_dir = os.path.dirname(current_dir)
sys.path.insert(0, parent_dir)

import numpy as np
import matplotlib.pyplot as plt
from pdfsolver import PdfSolver, PdfGrid
from Learning import PDElearn
from datamanage import DataIO
from visualization import Visualize
from scipy.signal import savgol_filter
from sklearn.metrics import mean_squared_error
import time
import pdb
from __init__ import *

def advection():
    case = 'advection_marginal'

    dt = 0.05
    t0 = 0
    tend = 5

    dx = 0.05
    x0 = -3.0
    xend = 3.0

    dk = 0.1
    k0 = -1.3
    kend = 2.3

    du = 0.05
    u0 = -5
    uend = 3

    sigu0=1.1
    mink=-0.5
    maxk=0.5
    a = 1.0
    b = 0.0

    muk = 0.5
    sigk= 1.0
    mux = 0.5
    sigx = 0.7
    muU = 0
    sigU = 1.2
    rho = 0

    gridvars = {'u': [u0, uend, du], 'k': [k0, kend, dk], 't': [t0, tend, dt], 'x':[x0, xend, dx]}
    ICparams = {'fu0':'compact_gaussian',
              'fu0param': [mux, sigx, muU, sigU, rho],
              'fk':'uniform',
              'fkparam': [mink, maxk]}

    grid = PdfGrid(gridvars)
    S = PdfSolver(grid, ICparams, save=True, case=case)
    S.solve() # no need to return anything

def reaction():
    case = 'reaction_linear'

    dt = 0.01
    t0 = 0
    tend = 5

    du = 0.005
    u0 = -7
    uend = 7

    umean = 0.2
    sigu0 = 1.1
    k = 0.5

    gridvars = {'u': [u0, uend, du], 't': [t0, tend, dt]}
    ICparams = {'u0': umean,
              'fu0':'gaussian',
              'fu0param': sigu0,
              'k': k}

    grid = PdfGrid(gridvars)
    S = PdfSolver(grid, ICparams, save=True, case=case)
    S.solve() # no need to return anything

if __name__ == "__main__":

    if len(sys.argv)>1:
        if sys.argv[1] == 'reaction':
            reaction()
        elif sys.argv[1] == 'advection':
            advection()
        else:
            raise exception("wrong option")
    else:
        advection()

max_stars_repo_name: josephbakarji/learning-pdf

max_issues_repo_name: josephbakarji/learning-pdf

max_forks_repo_name: josephbakarji/learning-pdf

reason: import numpy,from scipy


text:
"""
Generic test utilities.

"""

import os
import re
import sys
import numpy as np
import inspect

__all__ = ['PytestTester', 'check_free_memory', '_TestPythranFunc']

class FPUModeChangeWarning(RuntimeWarning):
    """Warning about FPU mode change"""
    pass

class PytestTester:
    """
    Pytest test runner entry point.
    """

    def __init__(self, module_name):
        self.module_name = module_name

    def __call__(self, label="fast", verbose=1, extra_argv=None, doctests=False,
                 coverage=False, tests=None, parallel=None):
        import pytest

        module = sys.modules[self.module_name]
        module_path = os.path.abspath(module.__path__[0])

        pytest_args = ['--showlocals', '--tb=short']

        if doctests:
            raise ValueError("Doctests not supported")

        if extra_argv:
            pytest_args += list(extra_argv)

        if verbose and int(verbose) > 1:
            pytest_args += ["-" + "v"*(int(verbose)-1)]

        if coverage:
            pytest_args += ["--cov=" + module_path]

        if label == "fast":
            pytest_args += ["-m", "not slow"]
        elif label != "full":
            pytest_args += ["-m", label]

        if tests is None:
            tests = [self.module_name]

        if parallel is not None and parallel > 1:
            if _pytest_has_xdist():
                pytest_args += ['-n', str(parallel)]
            else:
                import warnings
                warnings.warn('Could not run tests in parallel because '
                              'pytest-xdist plugin is not available.')

        pytest_args += ['--pyargs'] + list(tests)

        try:
            code = pytest.main(pytest_args)
        except SystemExit as exc:
            code = exc.code

        return (code == 0)

class _TestPythranFunc:
    '''
    These are situations that can be tested in our pythran tests:
    - A function with multiple array arguments and then
      other positional and keyword arguments.
    - A function with array-like keywords (e.g. `def somefunc(x0, x1=None)`.
    Note: list/tuple input is not yet tested!

    `self.arguments`: A dictionary which key is the index of the argument,
                      value is tuple(array value, all supported dtypes)
    `self.partialfunc`: A function used to freeze some non-array argument
                        that of no interests in the original function
    '''
    ALL_INTEGER = [np.int8, np.int16, np.int32, np.int64, np.intc, np.intp]
    ALL_FLOAT = [np.float32, np.float64]
    ALL_COMPLEX = [np.complex64, np.complex128]

    def setup_method(self):
        self.arguments = {}
        self.partialfunc = None
        self.expected = None

    def get_optional_args(self, func):
        # get optional arguments with its default value,
        # used for testing keywords
        signature = inspect.signature(func)
        optional_args = {}
        for k, v in signature.parameters.items():
            if v.default is not inspect.Parameter.empty:
                optional_args[k] = v.default
        return optional_args

    def get_max_dtype_list_length(self):
        # get the the max supported dtypes list length in all arguments
        max_len = 0
        for arg_idx in self.arguments:
            cur_len = len(self.arguments[arg_idx][1])
            if cur_len > max_len:
                max_len = cur_len
        return max_len

    def get_dtype(self, dtype_list, dtype_idx):
        # get the dtype from dtype_list via index
        # if the index is out of range, then return the last dtype
        if dtype_idx > len(dtype_list)-1:
            return dtype_list[-1]
        else:
            return dtype_list[dtype_idx]

    def test_all_dtypes(self):
        for type_idx in range(self.get_max_dtype_list_length()):
            args_array = []
            for arg_idx in self.arguments:
                new_dtype = self.get_dtype(self.arguments[arg_idx][1],
                                           type_idx)
                args_array.append(self.arguments[arg_idx][0].astype(new_dtype))
            self.pythranfunc(*args_array)

    def test_views(self):
        args_array = []
        for arg_idx in self.arguments:
            args_array.append(self.arguments[arg_idx][0][::-1][::-1])
        self.pythranfunc(*args_array)

    def test_strided(self):
        args_array = []
        for arg_idx in self.arguments:
            args_array.append(np.repeat(self.arguments[arg_idx][0],
                                        2, axis=0)[::2])
        self.pythranfunc(*args_array)

def _pytest_has_xdist():
    """
    Check if the pytest-xdist plugin is installed, providing parallel tests
    """
    # Check xdist exists without importing, otherwise pytests emits warnings
    from importlib.util import find_spec
    return find_spec('xdist') is not None

def check_free_memory(free_mb):
    """
    Check *free_mb* of memory is available, otherwise do pytest.skip
    """
    import pytest

    try:
        mem_free = _parse_size(os.environ['SCIPY_AVAILABLE_MEM'])
        msg = '{0} MB memory required, but environment SCIPY_AVAILABLE_MEM={1}'.format(
            free_mb, os.environ['SCIPY_AVAILABLE_MEM'])
    except KeyError:
        mem_free = _get_mem_available()
        if mem_free is None:
            pytest.skip("Could not determine available memory; set SCIPY_AVAILABLE_MEM "
                        "variable to free memory in MB to run the test.")
        msg = '{0} MB memory required, but {1} MB available'.format(
            free_mb, mem_free/1e6)

    if mem_free < free_mb * 1e6:
        pytest.skip(msg)

def _parse_size(size_str):
    suffixes = {'': 1e6,
                'b': 1.0,
                'k': 1e3, 'M': 1e6, 'G': 1e9, 'T': 1e12,
                'kb': 1e3, 'Mb': 1e6, 'Gb': 1e9, 'Tb': 1e12,
                'kib': 1024.0, 'Mib': 1024.0**2, 'Gib': 1024.0**3, 'Tib': 1024.0**4}
    m = re.match(r'^\s*(\d+)\s*({0})\s*$'.format('|'.join(suffixes.keys())),
                 size_str,
                 re.I)
    if not m or m.group(2) not in suffixes:
        raise ValueError("Invalid size string")

    return float(m.group(1)) * suffixes[m.group(2)]

def _get_mem_available():
    """
    Get information about memory available, not counting swap.
    """
    try:
        import psutil
        return psutil.virtual_memory().available
    except (ImportError, AttributeError):
        pass

    if sys.platform.startswith('linux'):
        info = {}
        with open('/proc/meminfo', 'r') as f:
            for line in f:
                p = line.split()
                info[p[0].strip(':').lower()] = float(p[1]) * 1e3

        if 'memavailable' in info:
            # Linux >= 3.14
            return info['memavailable']
        else:
            return info['memfree'] + info['cached']

    return None

max_stars_repo_name: lorentzenchr/scipy

max_issues_repo_name: lorentzenchr/scipy

max_forks_repo_name: lorentzenchr/scipy


text:
"""
Functions for forward modeling of thermochronometric ages
"""
import warnings

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from scipy.integrate import romb,romberg
from scipy.linalg import solve_banded
from scipy.interpolate import griddata

def tridiag(a, b, c, diag_length):
    """
    Set up a tridiagonal matrix from the values for the 3 diagonals (a,b,c)

    Parameters
    ----------
    a : float
        First diagonal value
    b : float
        Second diagonal value
    c : float
        Third diagonal value
    diag_length : flaot
        Length of principal diagonal

    Returns
    -------
    tridiag_matrix : NumPy array
        Tridiagonal matri

    """
    a_array = np.ones(diag_length-1)*a
    b_array = np.ones(diag_length)*b
    c_array = np.ones(diag_length-1)*c

    tridiag_matrix = (
        np.diag(a_array, -1) + np.diag(b_array, 0) + np.diag(c_array, 1)
        )

    return(tridiag_matrix)

def tridiag_banded(a, b, c, diag_length):
    """
    Set up a tridiagonal matrix in banded form
    from the values for the 3 diagonals (a,b,c).

    Parameters
    ----------
    a : float
        First diagonal value
    b : float
        Second diagonal value
    c : float
        Third diagonal value
    diag_length : flaot
        Length of principal diagonal

    Returns
    -------
    tridiag_matrix : NumPy array
        Tridiagonal matri

    """
    a_array = np.ones(diag_length)*a
    b_array = np.ones(diag_length)*b
    c_array = np.ones(diag_length)*c

    banded_matrix = np.vstack((a_array,b_array,c_array))
    banded_matrix[0,0] = 0
    banded_matrix[-1,-1] = 0

    return(banded_matrix)

def get_parameters(system='AHe'):
    """
    Get appropriate diffusion parameters for particular isotopic system.

    Frequency factors and activation energies after Reiners and Brandon, 2006
    and references therein. Stopping distances after Ketcham et al., 2011.

    Parameters
    ----------
    system : string, optional
        Shorthand name of the system. Options include AHe, ZHe, BtAr, MsAr,
        HbAr, and KsAr. The default is 'AHe'.

    Returns
    -------
    freq_factor : float
        Frequency factor (um^2*yr^-1)
    activ_energy : float
        Activation energy (J*mol^-1)
    stopping_distances : tuple
        Set of stopping distances (um) for U238, U235, and Th235. Only for
        AHe and ZHe.

    """
    if system == 'AHe':
        # Frequency Factor(um^2*yr^-1)
        freq_factor=50e8*3.154e7;
        # Activation Energy (J*mol^-1)
        activ_energy=138000
        # Stopping Distances (um)
        R_238U = 18.81
        R_235U = 21.80
        R_232Th = 22.25

    elif system == 'ZHe':
        # Frequency Factor(um^2*yr^-1)
        freq_factor=0.46e8*3.154e7;
        # Activation Energy (J*mol^-1)
        activ_energy=169000
        # Stopping Distances (um)
        R_238U = 15.55
        R_235U = 18.05
        R_232Th = 18.43

    elif system == 'BtAr':
        # Frequency Factor(um^2*yr^-1)
        freq_factor=7.5e-2*1e8*3.154e7;
        # Activation Energy (J*mol^-1)
        activ_energy=197000

    elif system == 'MsAr':
        # Frequency Factor(um^2*yr^-1)
        freq_factor=9.8e-3*1e8*3.154e7;
        # Activation Energy (J*mol^-1)
        activ_energy=18000

    elif system == 'HbAr':
        # Frequency Factor(um^2*yr^-1)
        freq_factor=6e-2*1e8*3.154e7;
        # Activation Energy (J*mol^-1)
        activ_energy=268000

    elif system == 'KsAr':
        # Frequency Factor(um^2*yr^-1)
        freq_factor=9.8e-3*1e8*3.154e7;
        # Activation Energy (J*mol^-1)
        activ_energy=183000

    else:
        raise Exception('Isotopic System Not Found')

    if (system == 'AHe') or (system == 'ZHe'):
        stopping_distances = np.array([R_238U,R_235U,R_232Th])
        return(freq_factor,activ_energy,stopping_distances)
    else:
        return(freq_factor,activ_energy)

def calculate_diffusivity(T,freq_factor,activ_energy,R=8.3144598):
    """
    Calculate diffusivity (kappa) based on temperature and system diffusion
    parameters

    After Reiners and Brandon, 2006, with PVa term assumed to be 0.

    Parameters
    ----------
    T : float
        Temperature (K)
    freq_factor : float
        Frequency factor (um^2*yr^-1)
    activ_energy : float
        Activation energy (J*mol^-1)
    R : float, optional
        Gas constant (J*K^-1*mol^-1). The default is 8.3144598.

    Returns
    -------
    kappa : float
        Diffusivity

    """
    kappa = freq_factor*np.exp(-activ_energy/(R*T))
    return(kappa)

def calculate_beta(diffusivity,node_spacing,time_interval):
    """
    Calculate beta, a substitution term based on diffusivity, the spacing
    of nodes within the modeled grain, and the timestep. After Ketcham, 2005.

    Parameters
    ----------
    diffusivity : float
        Diffusivity.
    node_spacing : float
        Spacing of nodes in the modeled crystal (um)
    time_interval : float
        Timestep in the thermal model (yr)

    Returns
    -------
    beta : float
        Beta, after Ketcham, 2005.

    """
    beta = (2*(node_spacing**2))/(diffusivity*time_interval)
    return(beta)

def UTh_ppm2molg(U,Th):
    """
    Convert concentrations of U and Th from ppm to mol/g

    Parameters
    ----------
    U : float
        U concentration (ppm)
    Th : float
        Th concentration (ppm)

    Returns
    -------
    U238_molg : float
        U238 (mol/g)
    U235_molg : float
        U235 (mol/g)
    Th_molg : float
        Th232 (mol/g)
    """

    U238_ppm = (137.88/(1+137.88))*U
    U235_ppm = (1/(1+137.88))*U
    Th_ppm = Th

    U238_molg = U238_ppm*1e-6/238
    U235_molg = U235_ppm*1e-6/235
    Th_molg = Th_ppm*1e-6/232

    return(U238_molg,U235_molg,Th_molg)

def calculate_He_production_rate(U238_molg,U235_molg,Th_molg):
    """
    Calculate He production rate as a function of U and Th.

    Parameters
    ----------
    U238_molg : float
        U238 (mol/g)
    U235_molg : float
        U235 (mol/g)
    Th_molg : float
        Th232 (mol/g)

    Returns
    -------
    He_production : float
        He production (mol/g)

    """

    lambda238 = -np.log(1/2)/4.468e9
    lambda235 = -np.log(1/2)/7.04e8
    lambda232 = -np.log(1/2)/1.40e10

    term238 = 8*lambda238*U238_molg
    term235 = 7*lambda235*U235_molg
    term232 = 6*lambda232*Th_molg

    He_production = term238+term235+term232

    return(He_production)

def calculate_node_positions(node_spacing,radius):
    """
    Calculate node positions given spacing and radius.

    Follows Ketcham, 2005.

    Parameters
    ----------
    node_spacing : float
        Distance between nodes in the crystal (um)
    radius : float
                Radius of the grain (um)

    Returns
    -------
    node_positions : NumPy array
        Radial positions of each modeled node (um)

    """
    node_positions = np.arange(node_spacing/2,radius,node_spacing)

    return(node_positions)

def sum_He_shells(x,node_positions,radius,nodes,print_He=True):
    """
    Sum He produced within all nodes of the modeled crystal.

    Uses substition for He volume after Ketcham, 2005. Converts radial profile
    to system of shells, so He is weighted by volume of shell. Also prints out
    He volume for comparison with real data.

    Parameters
    ----------
    x : NumPy array
        Matrix x solved for using matrices A and B after Ketcham, 2005.
        Equivalent to the volume times the node position.
    node_positions : NumPy array
        Radial positions of each modeled node (um)
    radius : float
            Radius of the grain (um)
    nodes : float
            Number of modeled nodes in the crystal

    Returns
    -------
    He_molg : float
        Total amount (mol/g) of He within the modeled crystal.
    v : NumPy array
        Radial profile of He (mol/g)

    """

    # Back-substitute u=vr to get radial profile
    v = (x/node_positions)

    # Get volumes of spheres at each node
    sphere_volumes = (node_positions)**3 * (4*np.pi/3)

    # Get total volume of the sphere
    total_volume = radius**3 * (4*np.pi/3)

    # Calculate volumes for the shell corresponding to each node
    shell_volumes = np.empty(sphere_volumes.size)
    for x in range(sphere_volumes.size):
        if x==0:
            shell_volumes[x]=sphere_volumes[x]
        else:
            shell_volumes[x] = sphere_volumes[x] - sphere_volumes[x-1]

    # Get shell as fraction of total volume
    shell_fraction = shell_volumes/total_volume

    # Scale He within radial profile by shell fraction
    v_shells = v*shell_fraction

    # Integrate weighted radial profile
    He_molg = romb(v_shells)
    He_nccg = He_molg * 22.4e12

    if print_He==True:
        print('He (ncc/g): ',He_nccg)
        print('He (nmol/g): ',He_molg*1e9)

    return(He_molg,v)

def calculate_age(He_molg,U238_molg,U235_molg,Th_molg):
    """
    Calculate (U-Th)/He age.

    Note that no alpha correction is applied here. Instead, the alpha
    correction is applied to the amounts of each parent isotope fed into
    this function, following Ketcham et al., 2011.

    Parameters
    ----------
    He_molg : float
        Amount of He (mol/g)
    U238_molg : float
        Amount of U238 (mol/g)
    U235_molg : float
        Amount of U235 (mol/g)
    Th_molg : float
        Amount of Th232 (mol/g)

    Returns
    -------
    age_Ma : float
        Calculated (U-Th)/He age

    """

    lambda238 = -np.log(1/2)/4.468e9
    lambda235 = -np.log(1/2)/7.04e8
    lambda232 = -np.log(1/2)/1.40e10

    ageterm_238 = 8*U238_molg
    ageterm_235 = 7*U235_molg
    ageterm_232 = 6*Th_molg

    def age_equation(t):
        root = (
            ageterm_238*(np.exp(lambda238*t)-1)
            + ageterm_235*(np.exp(lambda235*t)-1)
            + ageterm_232*(np.exp(lambda232*t)-1) - He_molg
            )

        return(root)

    warnings.filterwarnings('ignore','The iteration is not making good progress')

    age = float(fsolve(age_equation,1e6))
    age_Ma = age/1e6

    return(age_Ma)

def alpha_correction(stopping_distance,radius):
    """
    Calculate alpha ejection correction factor, after Ketcham et al., 2011.

    Parameters
    ----------
    stopping_distance : float
        Stopping distance for particular isotopic system (um).
    radius : float
        Radius of the grain (um)

    Returns
    -------
    tau : float
        Alpha correction factor

    """
    volume = (4/3) * np.pi * radius**3
    surface_area = 4 * np.pi * radius**2

    tau = 1 - 0.25*((surface_area*stopping_distance)/volume)

    return(tau)

def model_alpha_ejection(node_positions,stopping_distance,radius):
    """
    Model retained fraction of He after alpha ejection, given node position,
    stopping distances, and radius.

    Parameters
    ----------
    node_positions : NumPy array
        Radial positions of each modeled node (um)
    stopping_distance : float
        Stopping distance for particular isotopic system (um).
    radius : float
        Radius of the grain (um)

    Returns
    -------
    retained_fraction_edge : NumPy array
        Fraction of He retained after alpha ejection for each node position

    """

    # Find edge nodes based on stopping distance and radius
    edge_nodes = node_positions >= radius-stopping_distance

    # Calculate intersection planes for all nodes
    intersection_planes = (
        (node_positions**2 + radius**2 - stopping_distance**2)/(2*node_positions)
        )

    # Calculate retained fractions for all nodes hypothetically
    retained_fractions_all = (
        0.5 + (intersection_planes-node_positions)/(2*stopping_distance)
        )

    # Only apply retained fraction to edge nodes
    retained_fractions_edge = np.where(edge_nodes,retained_fractions_all,1)

    return(retained_fractions_edge)

def forward_model(U,Th,radius,temps,time_interval,system,nodes=513,
                  initial_He=np.nan,calc_age=True,print_age=True):
    """
    Forward model a (U-Th)/He age for a particular time-temperature path.

    Uses finite difference method for diffusion within a sphere as described
    in Ketcham, 2005. Applies alpha ejection correction after Ketcham et al.,
    2011. Returns corrected age but uncorrected age also printed.

    Parameters
    ----------
    U : float
        U concentration (ppm)
    Th : float
        Th concentration (ppm)
    radius : float
        Radius of the grain (um)
    temps : NumPy array
        List of temperatures (K) for the time-temperature path
    time_interval : float
        Time (yr) between each temperature in the time-temperature path.
    system : string
        Isotopic system. Current options are 'AHe' and 'ZHe'.
    nodes : float, optional
        Number of nodes to model within the crystal. The default is 513.

    Returns
    -------
    age_corrected : float
        (U-Th)/He age (Ma) corrected for alpha ejection

    """

    # Find node spacing and time interval based on radius and T-t path
    node_spacing = radius/nodes

    node_positions = calculate_node_positions(node_spacing,radius)

    # Get parameters for the appropriate mineral
    freq_factor,activ_energy,stop_distances = get_parameters(system)

    # Get mol/g of U,Th
    U238_molg,U235_molg,Th_molg = UTh_ppm2molg(U,Th)

    # Package parameters to pass to He profile

    node_information = (nodes,node_spacing,node_positions)
    UTh_molg = (U238_molg,U235_molg,Th_molg)
    system_parameters = (freq_factor,activ_energy,stop_distances)

    # Calculate He profile
    x = He_profile(temps,time_interval,node_information,radius,UTh_molg,
                   system_parameters,initial_He)

    if calc_age==True:

        # Calculate corrected age age
        age,vol,pos = profile2age(x,node_positions,radius,nodes,UTh_molg,
                                  stop_distances,print_age)

        return(age,vol,pos,x)

    else:
        return(x)

def He_profile(temps,time_interval,node_information,radius,UTh_molg,
               system_parameters,initial_He=np.nan):

    # Unpack parameters
    nodes,node_spacing,node_positions = node_information
    U238_molg,U235_molg,Th_molg = UTh_molg
    freq_factor,activ_energy,stop_distances = system_parameters

    # Modify mol/g of U,Th for alpha ejection
    U238_alpha = (
        U238_molg*model_alpha_ejection(node_positions,stop_distances[0],
                                                radius)
        )

    U235_alpha = (
        U235_molg*model_alpha_ejection(node_positions,stop_distances[1],
                                                radius)
        )

    Th_alpha = (
        Th_molg*model_alpha_ejection(node_positions,stop_distances[2],
                                                radius)
        )

    # Calculate He production based on U and Th, adjusted for alpha ejection.
    He_production = calculate_He_production_rate(U238_alpha,U235_alpha,Th_alpha)

    if np.all(np.isnan(initial_He)):
        # Set initial x (He) equal to 0
        x = np.zeros(nodes)
    else:
        x = initial_He

    # Loop through each step of the T-t path

    for temp in temps:

        # Use temperature to calculate diffusivity and beta

        diffusivity = calculate_diffusivity(temp,freq_factor,activ_energy)
        beta = calculate_beta(diffusivity,node_spacing, time_interval)

       # Calculate A (banded) and use x to calculate new B

        AB = tridiag_banded(1,-2-beta,1,nodes)
        AB[1,0]=-3-beta # Following Guenthner Matlab script - Neumann
        B = np.zeros(nodes)
        # Loop through each node
        for j in range(len(node_positions)):
            # For first node, use boundary condition (Neumann)
            if j==0:
                B[j] = (
                    x[j] + (2-beta)*x[j] - x[j+1]
                    - He_production[j]*node_positions[j]*beta*time_interval
                    )

            # For last node, use boundary condition (Drichlet)
            elif j==(len(B)-1):
                B[j] = (
                    -x[j-1] + (2-beta)*x[j] - 0
                    - He_production[j]*node_positions[j]*beta*time_interval
                    )
            # Use previous and subsequent nodes for remaining nodes
            else:
                B[j] = (
                    -x[j-1] + (2-beta)*x[j] - x[j+1]
                    - He_production[j]*node_positions[j]*beta*time_interval
                    )

        # Solve for x using banded A and B
        x = solve_banded((1,1),AB,B)

    return(x)

def profile2age(x,node_positions,radius,nodes,UTh_molg,stop_distances,
                print_age=True,):

    He_molg,volumes = sum_He_shells(x,node_positions,radius,nodes,
                                    print_He=print_age)

    U238_molg,U235_molg,Th_molg = UTh_molg

    # Because alpha ejection modeled, model age is an "uncorrected" age.
    age_uncorrected = calculate_age(He_molg,U238_molg,U235_molg,Th_molg)

    # "Corrected" age uses alpha-adjusted U-Th values

    # Get array of correction values (238,235,232)
    tau = alpha_correction(stop_distances,radius)

    # Correct U and Th accordingly
    U238_corr = U238_molg*tau[0]
    U235_corr = U235_molg*tau[1]
    Th_corr = Th_molg*tau[2]

    age_corrected = calculate_age(He_molg,U238_corr,U235_corr,Th_corr)

    if print_age==True:
        print('Age (Ma) Uncorrected: ',age_uncorrected)
        print('Age (Ma) Corrected: ',age_corrected)

    # Make diffusional profile
    volumes_normalized = volumes/np.max(volumes)
    position_normalized = node_positions/radius

    return(age_corrected,volumes_normalized,position_normalized)

def interpolate(x,y,age,**kwargs):
    all_parts = (x,y)

    known_x = x[~np.isnan(age)]
    known_y = y[~np.isnan(age)]
    knowns = (known_x,known_y)

    known_ages = age[~np.isnan(age)]

    # Make negative ages 0
    positive_ages = np.where(known_ages>0,known_ages,0)

    filled_grid = griddata(knowns,positive_ages,all_parts,**kwargs)

    return(filled_grid)

max_stars_repo_name: dyvasey/riftinversion

max_issues_repo_name: dyvasey/riftinversion

max_forks_repo_name: dyvasey/riftinversion


text:
import six
import chainer
import numpy as np
import chainer.links as L
import chainer.functions as F
import nutszebra_chainer
import functools
from collections import defaultdict

class BN_ReLU_Conv(nutszebra_chainer.Model):

    def __init__(self, in_channel, out_channel, filter_size=(3, 3), stride=(1, 1), pad=(1, 1)):
        super(BN_ReLU_Conv, self).__init__(
            conv=L.Convolution2D(in_channel, out_channel, filter_size, stride, pad),
            bn=L.BatchNormalization(in_channel),
        )

    def weight_initialization(self):
        self.conv.W.data = self.weight_relu_initialization(self.conv)
        self.conv.b.data = self.bias_initialization(self.conv, constant=0)

    def __call__(self, x, train=False):
        return self.conv(F.relu(self.bn(x, test=not train)))

    def count_parameters(self):
        return functools.reduce(lambda a, b: a * b, self.conv.W.data.shape)

class ResBlock(nutszebra_chainer.Model):

    def __init__(self, in_channel, out_channel, n=18, stride_at_first_layer=2, multiplier=4):
        super(ResBlock, self).__init__()
        modules = []
        modules += [('skip_bn', L.BatchNormalization(in_channel))]
        modules += [('bn_relu_conv1_1', BN_ReLU_Conv(in_channel, out_channel, 1, stride_at_first_layer, 0))]
        modules += [('bn_relu_conv2_1', BN_ReLU_Conv(out_channel, out_channel))]
        modules += [('bn_relu_conv3_1', BN_ReLU_Conv(out_channel, int(multiplier * out_channel), 1, 1, 0))]
        for i in six.moves.range(2, n + 1):
            modules.append(('bn_relu_conv1_{}'.format(i), BN_ReLU_Conv(int(multiplier * out_channel), out_channel, 1, 1, 0)))
            modules.append(('bn_relu_conv2_{}'.format(i), BN_ReLU_Conv(out_channel, out_channel)))
            modules.append(('bn_relu_conv3_{}'.format(i), BN_ReLU_Conv(out_channel, int(multiplier * out_channel), 1, 1, 0)))
        # register layers
        [self.add_link(*link) for link in modules]
        self.modules = modules
        self.in_channel = in_channel
        self.out_channel = out_channel
        self.n = n
        self.stride_at_first_layer = stride_at_first_layer
        self.multiplier = multiplier

    def weight_initialization(self):
        for i in six.moves.range(1, self.n + 1):
            self['bn_relu_conv1_{}'.format(i)].weight_initialization()
            self['bn_relu_conv2_{}'.format(i)].weight_initialization()
            self['bn_relu_conv3_{}'.format(i)].weight_initialization()

    @staticmethod
    def concatenate_zero_pad(x, h_shape, volatile, h_type):
        _, x_channel, _, _ = x.data.shape
        batch, h_channel, h_y, h_x = h_shape
        if x_channel == h_channel:
            return x
        pad = chainer.Variable(np.zeros((batch, h_channel - x_channel, h_y, h_x), dtype=np.float32), volatile=volatile)
        if h_type is not np.ndarray:
            pad.to_gpu()
        return F.concat((x, pad))

    def maybe_pooling(self, x):
        if self.stride_at_first_layer == 2:
            return F.average_pooling_2d(x, 1, 2, 0)
        return x

    def __call__(self, x, train=False):
        h = self['bn_relu_conv1_1'](x, train=train)
        h = self['bn_relu_conv2_1'](h, train=train)
        h = self['bn_relu_conv3_1'](h, train=train)
        x = h + ResBlock.concatenate_zero_pad(self.maybe_pooling(F.relu(self.skip_bn(x, test=not train))), h.data.shape, h.volatile, type(h.data))
        for i in six.moves.range(2, self.n + 1):
            h = self['bn_relu_conv1_{}'.format(i)](x, train=train)
            h = self['bn_relu_conv2_{}'.format(i)](h, train=train)
            x = self['bn_relu_conv3_{}'.format(i)](h, train=train) + x
        return x

    def count_parameters(self):
        count = 0
        for i in six.moves.range(1, self.n + 1):
            count += self['bn_relu_conv1_{}'.format(i)].count_parameters()
            count += self['bn_relu_conv2_{}'.format(i)].count_parameters()
            count += self['bn_relu_conv3_{}'.format(i)].count_parameters()
        return count

class ResidualNetwork(nutszebra_chainer.Model):

    def __init__(self, category_num, block_num=3, out_channels=(16, 32, 64), N=(18, 18, 18), multiplier=4):
        super(ResidualNetwork, self).__init__()
        # conv
        modules = [('conv1', L.Convolution2D(3, out_channels[0], 3, 1, 1))]
        in_channel = out_channels[0]
        strides = [1] + [2] * (block_num - 1)
        for i, out_channel, n, stride in six.moves.zip(six.moves.range(1, block_num + 1), out_channels, N, strides):
            modules.append(('res_block{}'.format(i), ResBlock(in_channel, out_channel, n=n, stride_at_first_layer=stride)))
            in_channel = int(out_channel * multiplier)
        modules.append(('bn_relu_conv', BN_ReLU_Conv(in_channel, category_num, filter_size=(1, 1), stride=(1, 1), pad=(0, 0))))
        # register layers
        [self.add_link(*link) for link in modules]
        self.modules = modules
        self.category_num = category_num
        self.block_num = block_num
        self.out_channels = out_channels
        self.N = N
        self.name = 'residual_network_{}_{}_{}_{}_{}'.format(category_num, block_num, out_channels, N, multiplier)

    def weight_initialization(self):
        self.conv1.W.data = self.weight_relu_initialization(self.conv1)
        self.conv1.b.data = self.bias_initialization(self.conv1, constant=0)
        for i in six.moves.range(1, self.block_num + 1):
            self['res_block{}'.format(i)].weight_initialization()
        self.bn_relu_conv.weight_initialization()

    def __call__(self, x, train=False):
        h = F.relu(self.conv1(x))
        for i in six.moves.range(1, self.block_num + 1):
            h = self['res_block{}'.format(i)](h, train=train)
        h = self.bn_relu_conv(F.relu(h), train=train)
        num, categories, y, x = h.data.shape
        h = F.reshape(F.average_pooling_2d(h, (y, x)), (num, categories))
        return h

    def count_parameters(self):
        count = 0
        count += functools.reduce(lambda a, b: a * b, self.conv1.W.data.shape)
        for i in six.moves.range(1, self.block_num + 1):
            count += self['res_block{}'.format(i)].count_parameters()
        count += self.bn_relu_conv.count_parameters()
        return count

    def calc_loss(self, y, t):
        loss = F.softmax_cross_entropy(y, t)
        return loss

    def accuracy(self, y, t, xp=np):
        y.to_cpu()
        t.to_cpu()
        indices = np.where((t.data == np.argmax(y.data, axis=1)) == True)[0]
        accuracy = defaultdict(int)
        for i in indices:
            accuracy[t.data[i]] += 1
        indices = np.where((t.data == np.argmax(y.data, axis=1)) == False)[0]
        false_accuracy = defaultdict(int)
        false_y = np.argmax(y.data, axis=1)
        for i in indices:
            false_accuracy[(t.data[i], false_y[i])] += 1
        return accuracy, false_accuracy


text:
import dadi
import dadi.DFE
#import dadi.NLopt_mod
import nlopt
import os, time
import numpy as np
from src.Pdfs import get_dadi_pdf

def infer_dfe(fs, cache1d, cache2d, sele_dist, sele_dist2, theta,
              p0, upper_bounds, lower_bounds, fixed_params, misid, cuda, maxeval, maxtime, seed):

    # Randomize starting parameter values
    if seed != None:
        np.random.seed(seed)
    else:
        ts = time.time()
        seed = int(time.time()) + int(os.getpid())
        np.random.seed(seed)

    # fs = dadi.Spectrum.from_file(fs)

    # theta = ns_s * theta
    # #popt = np.array(open(popt, 'r').readline().rstrip().split(), dtype=float)
    # #theta = ns_s * popt[-1]

    if cache1d != None:
        # spectra1d = pickle.load(open(cache1d, 'rb'))
        func = cache1d.integrate
    if cache2d != None:
        # spectra2d = pickle.load(open(cache2d, 'rb'))
        func = cache2d.integrate

    if sele_dist != None:
        sele_dist = get_dadi_pdf(sele_dist)
    if sele_dist2 != None:
        sele_dist2 = get_dadi_pdf(sele_dist2)
    if (sele_dist == None) and (sele_dist2 != None):
        sele_dist = sele_dist2

    if (cache1d != None) and (cache2d != None):
        func = dadi.DFE.Cache2D_mod.mixture
        func_args = [cache1d, cache2d, sele_dist, sele_dist2, theta]
    else:
        func_args = [sele_dist, theta]

    if misid:
        func = dadi.Numerics.make_anc_state_misid_func(func)

    # print('Inputs\n',p0,'\n',lower_bounds,'\n',upper_bounds,'\n',fixed_params,'\n\n')
    p0_len = len(p0)
    p0 = _convert_to_None(p0, p0_len)
    lower_bounds = _convert_to_None(lower_bounds, p0_len)
    upper_bounds = _convert_to_None(upper_bounds, p0_len)
    fixed_params = _convert_to_None(fixed_params, p0_len)

    # Fit a DFE to the data
    # Initial guess and bounds
    #print(p0)
    # print('Inputs\n',p0,'\n',lower_bounds,'\n',upper_bounds,'\n',fixed_params,'\n\n')
    print(func)
    p0 = dadi.Misc.perturb_params(p0, lower_bound=lower_bounds, upper_bound=upper_bounds)
    popt, _ = dadi.Inference.opt(p0, fs, func, pts=None,
                                func_args=func_args, fixed_params=fixed_params,
                                lower_bound=lower_bounds, upper_bound=upper_bounds,
                                maxeval=maxeval, maxtime=maxtime, multinom=False, verbose=0)
    # print(popt)

    #print('Optimized parameters: {0}'.format(popt))

    # Get expected SFS for MLE
    if (cache1d != None) and (cache2d != None):
        model = func(popt, None, cache1d, cache2d, sele_dist, sele_dist2, theta, None)
    else:
        model = func(popt, None, sele_dist, theta, None)
    #print(model)
    # Likelihood of the data given the model AFS.
    ll_model = dadi.Inference.ll(model, fs)
    #print('Maximum log composite likelihood: {0}'.format(ll_model))

    #with open(output, 'w') as f:
    #    f.write(str(ll_model))
    #    for p in popt:
    #        f.write("\t")
    #        f.write(str(p))
    #    f.write("\n")
    return ll_model, popt, theta

def _convert_to_None(inference_input, p0_len):
    if inference_input == -1: inference_input = [inference_input]*p0_len
    inference_input = list(np.where(np.array(inference_input) == -1, None, np.array(inference_input)))
    return inference_input

max_stars_repo_name: xin-huang/dadi-cli

max_issues_repo_name: xin-huang/dadi-cli

max_forks_repo_name: xin-huang/dadi-cli


text:
import cv2
import numpy as np
import glob
import os
def calibrate(dirpath):
    """ Apply camera calibration operation for images in the given directory path. """
    # prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(8,6,0)
    CHESS_SQUARE_LEN = 0.027

    # termination criteria
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
    # prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
    objp = np.zeros((6*9,3), np.float32)
    objp[:,:2] = np.mgrid[0:9,0:6].T.reshape(-1,2)*CHESS_SQUARE_LEN
    # Arrays to store object points and image points from all the images.
    objpoints = [] # 3d point in real world space
    imgpoints = [] # 2d points in image plane.
    images = glob.glob(os.path.join(dirpath,"*.jpg"))
    gray = None
    for fname in images:
        print(fname)
        img = cv2.imread(fname)
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        # Find the chess board corners
        ret, corners = cv2.findChessboardCorners(gray, (9,6), None)
        # If found, add object points, image points (after refining them)
        if ret == True:
            objpoints.append(objp)
            corners2 = cv2.cornerSubPix(gray,corners, (11,11), (-1,-1), criteria)
            imgpoints.append(corners)
            # Draw and display the corners
            #cv2.drawChessboardCorners(img, (9,6), corners2, ret)
            #cv2.imshow('img', img)
            #cv2.waitKey(500)

    ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1], None, None)

    return [ret, mtx, dist, rvecs, tvecs]

def saveCoefficients(mtx, dist, path):
    """ Save the camera matrix and the distortion coefficients to given path/file. """
    cv_file = cv2.FileStorage(path, cv2.FILE_STORAGE_WRITE)
    cv_file.write("camera_matrix", mtx)
    cv_file.write("dist_coeff", dist)
    # note you *release* you don't close() a FileStorage object
    cv_file.release()

if __name__ == '__main__':
    ret, matrix_coefficients, distortion_coefficients, rvecs, tvecs = calibrate('calib_images')
    if not os.path.exists( "./params"):
        os.makedirs( "./params")
    print('rvecs')
    for item in rvecs:
        print(item)
        print()

    print('tvecs')
    for item in tvecs:
        print(item)
        print()
    saveCoefficients(matrix_coefficients, distortion_coefficients, "./params/calibration_coefficients.yaml")

    #h,  w = gray.shape[:2]
    #newcameramtx, roi = cv2.getOptimalNewCameraMatrix(matrix_coefficients, distortion_coefficients, (w,h), 1, (w,h))
    # undistort
    #dst = cv2.undistort(gray, matrix_coefficients, distortion_coefficients, None, newcameramtx)


text:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import gym
from gym.envs.registration import register
register(id = "FrozenLakeEasy-v0", entry_point = "gym.envs.toy_text:FrozenLakeEnv",
         kwargs = {"is_slippery": False})

################################################################################
def show_q_value(Q):
    """
    Show Q-values for FrozenLake-v0.
    To show each action's evaluation,
    a state is shown as 3 x 3 matrix like following.

    +---+---+---+
    |   | u |   |  u: up value
    | l | m | r |  l: left value, r: right value, m: mean value
    |   | d |   |  d: down value
    +---+---+---+
    """
    env = gym.make("FrozenLake-v0")
    nrow = env.unwrapped.nrow
    ncol = env.unwrapped.ncol
    state_size = 3
    q_nrow = nrow * state_size
    q_ncol = ncol * state_size
    reward_map = np.zeros((q_nrow, q_ncol))
    #
    for r in range(nrow):
        for c in range(ncol):
            s = r * nrow + c
            state_exist = False
            if isinstance(Q, dict) and s in Q:
                state_exist = True
            elif isinstance(Q, (np.ndarray, np.generic)) and s < Q.shape[0]:
                state_exist = True

            if state_exist:
                # 在展示图中，纵向的序号是反转的
                _r = 1 + (nrow - 1 - r) * state_size
                _c = 1 + c * state_size
                reward_map[_r][_c - 1] = Q[s][0]  # LEFT = 0
                reward_map[_r - 1][_c] = Q[s][1]  # DOWN = 1
                reward_map[_r][_c + 1] = Q[s][2]  # RIGHT = 2
                reward_map[_r + 1][_c] = Q[s][3]  # UP = 3
                reward_map[_r][_c] = np.mean(Q[s])  # Center
    #
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    plt.imshow(reward_map, cmap = cm.RdYlGn, interpolation = "bilinear",
               vmax=abs(reward_map).max(), vmin = -abs(reward_map).max())
    ax.set_xlim(-0.5, q_ncol - 0.5)
    ax.set_ylim(-0.5, q_nrow - 0.5)
    ax.set_xticks(np.arange(-0.5, q_ncol, state_size))
    ax.set_yticks(np.arange(-0.5, q_nrow, state_size))
    ax.set_xticklabels(range(ncol + 1))
    ax.set_yticklabels(range(nrow + 1))
    ax.grid(which = "both")
    plt.show()

max_stars_repo_name: KhalilWong/Learn-RL

max_issues_repo_name: KhalilWong/Learn-RL

max_forks_repo_name: KhalilWong/Learn-RL


text:
import cv2
from numpy import array, matrix

# Camera width and height (for a single eye)
WIDTH = 960
HEIGHT = 960

# Camera device index or video file path
# Should be a dual fisheye camera. Left and Right eye views
# are assumed to be side-by-side and of the same pixel size.
VIDEO_INPUT = '../test-video.mp4'

# Path to output video file
VIDEO_OUT = 'out.mp4'

# Camera intrinsics, distortion parameters and stereo parameters computed with calibration.py
leftCameraMatrix = matrix([[418.12009366,   0.        , 474.70168291],
        [  0.        , 418.26107355, 475.82296102],
        [  0.        ,   0.        ,   1.        ]])
leftDistCoeffs = array([[-1.41472770e-01, -5.00821052e-03, -1.36898095e-04,
        -5.31109897e-04,  5.92505609e-03]])
rightCameraMatrix = matrix([[420.70164293,   0.        , 494.0718344 ],
        [  0.        , 421.00224821, 485.6593556 ],
        [  0.        ,   0.        ,   1.        ]])
rightDistCoeffs = array([[-1.47821076e-01,  2.58986389e-03, -5.24606308e-04,
         1.42863352e-04,  3.41432318e-03]])
R = array([[ 0.98734558,  0.02063431, -0.15723527],
       [-0.02137971,  0.99976677, -0.00305063],
       [ 0.15713565,  0.00637367,  0.98755646]])
T = array([[-0.0902457 ],
       [ 0.00017903],
       [-0.00636208]])

# Compute reprojection and correction matrices from above data
imgSize = (WIDTH, HEIGHT)
R1, R2, P1, P2, Q, roi1, roi2 = cv2.stereoRectify(leftCameraMatrix, leftDistCoeffs, rightCameraMatrix, rightDistCoeffs, imgSize, R, T, flags=cv2.CALIB_ZERO_DISPARITY)
leftMapX, leftMapY = cv2.initUndistortRectifyMap(leftCameraMatrix, leftDistCoeffs, R1, P1, imgSize, cv2.CV_32FC1)
rightMapX, rightMapY = cv2.initUndistortRectifyMap(rightCameraMatrix, rightDistCoeffs, R2, P2, imgSize, cv2.CV_32FC1)

cap = cv2.VideoCapture(VIDEO_INPUT)
frameSize = (int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)), int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)))
fourcc = cv2.VideoWriter_fourcc(*'mp4v')
out = cv2.VideoWriter(VIDEO_OUT, fourcc, float(cap.get(cv2.CAP_PROP_FPS)), frameSize)

while cap.isOpened():
  ret, frame = cap.read()
  if ret == True:
    left  = frame[0:HEIGHT, 0:WIDTH]
    right = frame[0:HEIGHT, WIDTH:WIDTH*2]

    left = cv2.remap(left, leftMapX, leftMapY, cv2.INTER_LINEAR)
    right = cv2.remap(right, rightMapX, rightMapY, cv2.INTER_LINEAR)

    frame = cv2.hconcat([left, right])
    cv2.imshow('Frame', frame)
    key = cv2.waitKey(1) & 0xFF
    out.write(frame)
  else:
    break

cap.release()
out.release()

cv2.destroyAllWindows()

max_stars_repo_name: ClearlyClaire/opencv-stereo-tests

max_issues_repo_name: ClearlyClaire/opencv-stereo-tests

max_forks_repo_name: ClearlyClaire/opencv-stereo-tests

reason: from numpy